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