1
2 /*! @file ilu_dpivotL.c
3 * \brief Performs numerical pivoting
4 *
5 * <pre>
6 * -- SuperLU routine (version 4.0) --
7 * Lawrence Berkeley National Laboratory
8 * June 30, 2009
9 * </pre>
10 */
11
12
13 #include <math.h>
14 #include <stdlib.h>
15 #include "slu_ddefs.h"
16
17 #ifndef SGN
18 #define SGN(x) ((x)>=0?1:-1)
19 #endif
20
21 /*! \brief
22 *
23 * <pre>
24 * Purpose
25 * =======
26 * Performs the numerical pivoting on the current column of L,
27 * and the CDIV operation.
28 *
29 * Pivot policy:
30 * (1) Compute thresh = u * max_(i>=j) abs(A_ij);
31 * (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
32 * pivot row = k;
33 * ELSE IF abs(A_jj) >= thresh THEN
34 * pivot row = j;
35 * ELSE
36 * pivot row = m;
37 *
38 * Note: If you absolutely want to use a given pivot order, then set u=0.0.
39 *
40 * Return value: 0 success;
41 * i > 0 U(i,i) is exactly zero.
42 * </pre>
43 */
44
45 int
ilu_dpivotL(const int jcol,const double u,int * usepr,int * perm_r,int diagind,int * swap,int * iswap,int * marker,int * pivrow,double fill_tol,milu_t milu,double drop_sum,GlobalLU_t * Glu,SuperLUStat_t * stat)46 ilu_dpivotL(
47 const int jcol, /* in */
48 const double u, /* in - diagonal pivoting threshold */
49 int *usepr, /* re-use the pivot sequence given by
50 * perm_r/iperm_r */
51 int *perm_r, /* may be modified */
52 int diagind, /* diagonal of Pc*A*Pc' */
53 int *swap, /* in/out record the row permutation */
54 int *iswap, /* in/out inverse of swap, it is the same as
55 perm_r after the factorization */
56 int *marker, /* in */
57 int *pivrow, /* in/out, as an input if *usepr!=0 */
58 double fill_tol, /* in - fill tolerance of current column
59 * used for a singular column */
60 milu_t milu, /* in */
61 double drop_sum, /* in - computed in ilu_dcopy_to_ucol()
62 (MILU only) */
63 GlobalLU_t *Glu, /* modified - global LU data structures */
64 SuperLUStat_t *stat /* output */
65 )
66 {
67
68 int n; /* number of columns */
69 int fsupc; /* first column in the supernode */
70 int nsupc; /* no of columns in the supernode */
71 int nsupr; /* no of rows in the supernode */
72 int lptr; /* points to the starting subscript of the supernode */
73 register int pivptr;
74 int old_pivptr, diag, ptr0;
75 register double pivmax, rtemp;
76 double thresh;
77 double temp;
78 double *lu_sup_ptr;
79 double *lu_col_ptr;
80 int *lsub_ptr;
81 register int isub, icol, k, itemp;
82 int *lsub, *xlsub;
83 double *lusup;
84 int *xlusup;
85 flops_t *ops = stat->ops;
86 int info;
87
88 /* Initialize pointers */
89 n = Glu->n;
90 lsub = Glu->lsub;
91 xlsub = Glu->xlsub;
92 lusup = Glu->lusup;
93 xlusup = Glu->xlusup;
94 fsupc = (Glu->xsup)[(Glu->supno)[jcol]];
95 nsupc = jcol - fsupc; /* excluding jcol; nsupc >= 0 */
96 lptr = xlsub[fsupc];
97 nsupr = xlsub[fsupc+1] - lptr;
98 lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
99 lu_col_ptr = &lusup[xlusup[jcol]]; /* start of jcol in the supernode */
100 lsub_ptr = &lsub[lptr]; /* start of row indices of the supernode */
101
102 /* Determine the largest abs numerical value for partial pivoting;
103 Also search for user-specified pivot, and diagonal element. */
104 pivmax = -1.0;
105 pivptr = nsupc;
106 diag = EMPTY;
107 old_pivptr = nsupc;
108 ptr0 = EMPTY;
109 for (isub = nsupc; isub < nsupr; ++isub) {
110 if (marker[lsub_ptr[isub]] > jcol)
111 continue; /* do not overlap with a later relaxed supernode */
112
113 switch (milu) {
114 case SMILU_1:
115 rtemp = fabs(lu_col_ptr[isub] + drop_sum);
116 break;
117 case SMILU_2:
118 case SMILU_3:
119 /* In this case, drop_sum contains the sum of the abs. value */
120 rtemp = fabs(lu_col_ptr[isub]);
121 break;
122 case SILU:
123 default:
124 rtemp = fabs(lu_col_ptr[isub]);
125 break;
126 }
127 if (rtemp > pivmax) { pivmax = rtemp; pivptr = isub; }
128 if (*usepr && lsub_ptr[isub] == *pivrow) old_pivptr = isub;
129 if (lsub_ptr[isub] == diagind) diag = isub;
130 if (ptr0 == EMPTY) ptr0 = isub;
131 }
132
133 if (milu == SMILU_2 || milu == SMILU_3) pivmax += drop_sum;
134
135 /* Test for singularity */
136 if (pivmax < 0.0) {
137 fprintf(stderr, "[0]: jcol=%d, SINGULAR!!!\n", jcol);
138 fflush(stderr);
139 exit(1);
140 }
141 if ( pivmax == 0.0 ) {
142 if (diag != EMPTY)
143 *pivrow = lsub_ptr[pivptr = diag];
144 else if (ptr0 != EMPTY)
145 *pivrow = lsub_ptr[pivptr = ptr0];
146 else {
147 /* look for the first row which does not
148 belong to any later supernodes */
149 for (icol = jcol; icol < n; icol++)
150 if (marker[swap[icol]] <= jcol) break;
151 if (icol >= n) {
152 fprintf(stderr, "[1]: jcol=%d, SINGULAR!!!\n", jcol);
153 fflush(stderr);
154 exit(1);
155 }
156
157 *pivrow = swap[icol];
158
159 /* pick up the pivot row */
160 for (isub = nsupc; isub < nsupr; ++isub)
161 if ( lsub_ptr[isub] == *pivrow ) { pivptr = isub; break; }
162 }
163 pivmax = fill_tol;
164 lu_col_ptr[pivptr] = pivmax;
165 *usepr = 0;
166 #ifdef DEBUG
167 printf("[0] ZERO PIVOT: FILL (%d, %d).\n", *pivrow, jcol);
168 fflush(stdout);
169 #endif
170 info =jcol + 1;
171 } /* if (*pivrow == 0.0) */
172 else {
173 thresh = u * pivmax;
174
175 /* Choose appropriate pivotal element by our policy. */
176 if ( *usepr ) {
177 switch (milu) {
178 case SMILU_1:
179 rtemp = fabs(lu_col_ptr[old_pivptr] + drop_sum);
180 break;
181 case SMILU_2:
182 case SMILU_3:
183 rtemp = fabs(lu_col_ptr[old_pivptr]) + drop_sum;
184 break;
185 case SILU:
186 default:
187 rtemp = fabs(lu_col_ptr[old_pivptr]);
188 break;
189 }
190 if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = old_pivptr;
191 else *usepr = 0;
192 }
193 if ( *usepr == 0 ) {
194 /* Use diagonal pivot? */
195 if ( diag >= 0 ) { /* diagonal exists */
196 switch (milu) {
197 case SMILU_1:
198 rtemp = fabs(lu_col_ptr[diag] + drop_sum);
199 break;
200 case SMILU_2:
201 case SMILU_3:
202 rtemp = fabs(lu_col_ptr[diag]) + drop_sum;
203 break;
204 case SILU:
205 default:
206 rtemp = fabs(lu_col_ptr[diag]);
207 break;
208 }
209 if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
210 }
211 *pivrow = lsub_ptr[pivptr];
212 }
213 info = 0;
214
215 /* Reset the diagonal */
216 switch (milu) {
217 case SMILU_1:
218 lu_col_ptr[pivptr] += drop_sum;
219 break;
220 case SMILU_2:
221 case SMILU_3:
222 lu_col_ptr[pivptr] += SGN(lu_col_ptr[pivptr]) * drop_sum;
223 break;
224 case SILU:
225 default:
226 break;
227 }
228
229 } /* else */
230
231 /* Record pivot row */
232 perm_r[*pivrow] = jcol;
233 if (jcol < n - 1) {
234 register int t1, t2, t;
235 t1 = iswap[*pivrow]; t2 = jcol;
236 if (t1 != t2) {
237 t = swap[t1]; swap[t1] = swap[t2]; swap[t2] = t;
238 t1 = swap[t1]; t2 = t;
239 t = iswap[t1]; iswap[t1] = iswap[t2]; iswap[t2] = t;
240 }
241 } /* if (jcol < n - 1) */
242
243 /* Interchange row subscripts */
244 if ( pivptr != nsupc ) {
245 itemp = lsub_ptr[pivptr];
246 lsub_ptr[pivptr] = lsub_ptr[nsupc];
247 lsub_ptr[nsupc] = itemp;
248
249 /* Interchange numerical values as well, for the whole snode, such
250 * that L is indexed the same way as A.
251 */
252 for (icol = 0; icol <= nsupc; icol++) {
253 itemp = pivptr + icol * nsupr;
254 temp = lu_sup_ptr[itemp];
255 lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
256 lu_sup_ptr[nsupc + icol*nsupr] = temp;
257 }
258 } /* if */
259
260 /* cdiv operation */
261 ops[FACT] += nsupr - nsupc;
262 temp = 1.0 / lu_col_ptr[nsupc];
263 for (k = nsupc+1; k < nsupr; k++) lu_col_ptr[k] *= temp;
264
265 return info;
266 }
267