← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/cpivotL.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/cpivotL.c
 
1
 
 
2
/*
 
3
 * -- SuperLU routine (version 3.0) --
 
4
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 
5
 * and Lawrence Berkeley National Lab.
 
6
 * October 15, 2003
 
7
 *
 
8
 */
 
9
/*
 
10
  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 
11
 
 
12
  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 
13
  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 
14
 
 
15
  Permission is hereby granted to use or copy this program for any
 
16
  purpose, provided the above notices are retained on all copies.
 
17
  Permission to modify the code and to distribute modified code is
 
18
  granted, provided the above notices are retained, and a notice that
 
19
  the code was modified is included with the above copyright notice.
 
20
*/
 
21
 
 
22
#include <math.h>
 
23
#include <stdlib.h>
 
24
#include "csp_defs.h"
 
25
 
 
26
#undef DEBUG
 
27
 
 
28
int
 
29
cpivotL(
 
30
        const int  jcol,     /* in */
 
31
        const float u,      /* in - diagonal pivoting threshold */
 
32
        int        *usepr,   /* re-use the pivot sequence given by perm_r/iperm_r */
 
33
        int        *perm_r,  /* may be modified */
 
34
        int        *iperm_r, /* in - inverse of perm_r */
 
35
        int        *iperm_c, /* in - used to find diagonal of Pc*A*Pc' */
 
36
        int        *pivrow,  /* out */
 
37
        GlobalLU_t *Glu,     /* modified - global LU data structures */
 
38
        SuperLUStat_t *stat  /* output */
 
39
       )
 
40
{
 
41
/*
 
42
 * Purpose
 
43
 * =======
 
44
 *   Performs the numerical pivoting on the current column of L,
 
45
 *   and the CDIV operation.
 
46
 *
 
47
 *   Pivot policy:
 
48
 *   (1) Compute thresh = u * max_(i>=j) abs(A_ij);
 
49
 *   (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
 
50
 *           pivot row = k;
 
51
 *       ELSE IF abs(A_jj) >= thresh THEN
 
52
 *           pivot row = j;
 
53
 *       ELSE
 
54
 *           pivot row = m;
 
55
 * 
 
56
 *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
 
57
 *
 
58
 *   Return value: 0      success;
 
59
 *                 i > 0  U(i,i) is exactly zero.
 
60
 *
 
61
 */
 
62
    complex one = {1.0, 0.0};
 
63
    int          fsupc;     /* first column in the supernode */
 
64
    int          nsupc;     /* no of columns in the supernode */
 
65
    int          nsupr;     /* no of rows in the supernode */
 
66
    int          lptr;      /* points to the starting subscript of the supernode */
 
67
    int          pivptr, old_pivptr, diag, diagind;
 
68
    float       pivmax, rtemp, thresh;
 
69
    complex       temp;
 
70
    complex       *lu_sup_ptr; 
 
71
    complex       *lu_col_ptr;
 
72
    int          *lsub_ptr;
 
73
    int          isub, icol, k, itemp;
 
74
    int          *lsub, *xlsub;
 
75
    complex       *lusup;
 
76
    int          *xlusup;
 
77
    flops_t      *ops = stat->ops;
 
78
 
 
79
    /* Initialize pointers */
 
80
    lsub       = Glu->lsub;
 
81
    xlsub      = Glu->xlsub;
 
82
    lusup      = Glu->lusup;
 
83
    xlusup     = Glu->xlusup;
 
84
    fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
 
85
    nsupc      = jcol - fsupc;          /* excluding jcol; nsupc >= 0 */
 
86
    lptr       = xlsub[fsupc];
 
87
    nsupr      = xlsub[fsupc+1] - lptr;
 
88
    lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
 
89
    lu_col_ptr = &lusup[xlusup[jcol]];  /* start of jcol in the supernode */
 
90
    lsub_ptr   = &lsub[lptr];   /* start of row indices of the supernode */
 
91
 
 
92
#ifdef DEBUG
 
93
if ( jcol == MIN_COL ) {
 
94
    printf("Before cdiv: col %d\n", jcol);
 
95
    for (k = nsupc; k < nsupr; k++) 
 
96
        printf("  lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
 
97
}
 
98
#endif
 
99
    
 
100
    /* Determine the largest abs numerical value for partial pivoting;
 
101
       Also search for user-specified pivot, and diagonal element. */
 
102
    if ( *usepr ) *pivrow = iperm_r[jcol];
 
103
    diagind = iperm_c[jcol];
 
104
    pivmax = 0.0;
 
105
    pivptr = nsupc;
 
106
    diag = EMPTY;
 
107
    old_pivptr = nsupc;
 
108
    for (isub = nsupc; isub < nsupr; ++isub) {
 
109
        rtemp = c_abs1 (&lu_col_ptr[isub]);
 
110
        if ( rtemp > pivmax ) {
 
111
            pivmax = rtemp;
 
112
            pivptr = isub;
 
113
        }
 
114
        if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
 
115
        if ( lsub_ptr[isub] == diagind ) diag = isub;
 
116
    }
 
117
 
 
118
    /* Test for singularity */
 
119
    if ( pivmax == 0.0 ) {
 
120
        *pivrow = lsub_ptr[pivptr];
 
121
        perm_r[*pivrow] = jcol;
 
122
        *usepr = 0;
 
123
        return (jcol+1);
 
124
    }
 
125
 
 
126
    thresh = u * pivmax;
 
127
    
 
128
    /* Choose appropriate pivotal element by our policy. */
 
129
    if ( *usepr ) {
 
130
        rtemp = c_abs1 (&lu_col_ptr[old_pivptr]);
 
131
        if ( rtemp != 0.0 && rtemp >= thresh )
 
132
            pivptr = old_pivptr;
 
133
        else
 
134
            *usepr = 0;
 
135
    }
 
136
    if ( *usepr == 0 ) {
 
137
        /* Use diagonal pivot? */
 
138
        if ( diag >= 0 ) { /* diagonal exists */
 
139
            rtemp = c_abs1 (&lu_col_ptr[diag]);
 
140
            if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
 
141
        }
 
142
        *pivrow = lsub_ptr[pivptr];
 
143
    }
 
144
    
 
145
    /* Record pivot row */
 
146
    perm_r[*pivrow] = jcol;
 
147
    
 
148
    /* Interchange row subscripts */
 
149
    if ( pivptr != nsupc ) {
 
150
        itemp = lsub_ptr[pivptr];
 
151
        lsub_ptr[pivptr] = lsub_ptr[nsupc];
 
152
        lsub_ptr[nsupc] = itemp;
 
153
 
 
154
        /* Interchange numerical values as well, for the whole snode, such 
 
155
         * that L is indexed the same way as A.
 
156
         */
 
157
        for (icol = 0; icol <= nsupc; icol++) {
 
158
            itemp = pivptr + icol * nsupr;
 
159
            temp = lu_sup_ptr[itemp];
 
160
            lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
 
161
            lu_sup_ptr[nsupc + icol*nsupr] = temp;
 
162
        }
 
163
    } /* if */
 
164
 
 
165
    /* cdiv operation */
 
166
    ops[FACT] += 10 * (nsupr - nsupc);
 
167
 
 
168
    c_div(&temp, &one, &lu_col_ptr[nsupc]);
 
169
    for (k = nsupc+1; k < nsupr; k++) 
 
170
        cc_mult(&lu_col_ptr[k], &lu_col_ptr[k], &temp);
 
171
 
 
172
    return 0;
 
173
}
 
174