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