1
2 /*! @file dpivotL.c
3 * \brief Performs numerical pivoting
4 *
5 * <pre>
6 * -- SuperLU routine (version 3.0) --
7 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
8 * and Lawrence Berkeley National Lab.
9 * October 15, 2003
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 * </pre>
22 */
23
24
25 #include <math.h>
26 #include <stdlib.h>
27 #include "slu_ddefs.h"
28
29 #undef DEBUG
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
dpivotL(const int jcol,const double u,int * usepr,int * perm_r,int * iperm_r,int * iperm_c,int * pivrow,GlobalLU_t * Glu,SuperLUStat_t * stat)56 dpivotL(
57 const int jcol, /* in */
58 const double u, /* in - diagonal pivoting threshold */
59 int *usepr, /* re-use the pivot sequence given by perm_r/iperm_r */
60 int *perm_r, /* may be modified */
61 int *iperm_r, /* in - inverse of perm_r */
62 int *iperm_c, /* in - used to find diagonal of Pc*A*Pc' */
63 int *pivrow, /* out */
64 GlobalLU_t *Glu, /* modified - global LU data structures */
65 SuperLUStat_t *stat /* output */
66 )
67 {
68
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 int pivptr, old_pivptr, diag, diagind;
74 double pivmax, rtemp, thresh;
75 double temp;
76 double *lu_sup_ptr;
77 double *lu_col_ptr;
78 int *lsub_ptr;
79 int isub, icol, k, itemp;
80 int *lsub, *xlsub;
81 double *lusup;
82 int *xlusup;
83 flops_t *ops = stat->ops;
84
85 /* Initialize pointers */
86 lsub = Glu->lsub;
87 xlsub = Glu->xlsub;
88 lusup = Glu->lusup;
89 xlusup = Glu->xlusup;
90 fsupc = (Glu->xsup)[(Glu->supno)[jcol]];
91 nsupc = jcol - fsupc; /* excluding jcol; nsupc >= 0 */
92 lptr = xlsub[fsupc];
93 nsupr = xlsub[fsupc+1] - lptr;
94 lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
95 lu_col_ptr = &lusup[xlusup[jcol]]; /* start of jcol in the supernode */
96 lsub_ptr = &lsub[lptr]; /* start of row indices of the supernode */
97
98 #ifdef DEBUG
99 if ( jcol == MIN_COL ) {
100 printf("Before cdiv: col %d\n", jcol);
101 for (k = nsupc; k < nsupr; k++)
102 printf(" lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
103 }
104 #endif
105
106 /* Determine the largest abs numerical value for partial pivoting;
107 Also search for user-specified pivot, and diagonal element. */
108 if ( *usepr ) *pivrow = iperm_r[jcol];
109 diagind = iperm_c[jcol];
110 pivmax = 0.0;
111 pivptr = nsupc;
112 diag = EMPTY;
113 old_pivptr = nsupc;
114 for (isub = nsupc; isub < nsupr; ++isub) {
115 rtemp = fabs (lu_col_ptr[isub]);
116 if ( rtemp > pivmax ) {
117 pivmax = rtemp;
118 pivptr = isub;
119 }
120 if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
121 if ( lsub_ptr[isub] == diagind ) diag = isub;
122 }
123
124 /* Test for singularity */
125 if ( pivmax == 0.0 ) {
126 *pivrow = lsub_ptr[pivptr];
127 perm_r[*pivrow] = jcol;
128 *usepr = 0;
129 return (jcol+1);
130 }
131
132 thresh = u * pivmax;
133
134 /* Choose appropriate pivotal element by our policy. */
135 if ( *usepr ) {
136 rtemp = fabs (lu_col_ptr[old_pivptr]);
137 if ( rtemp != 0.0 && rtemp >= thresh )
138 pivptr = old_pivptr;
139 else
140 *usepr = 0;
141 }
142 if ( *usepr == 0 ) {
143 /* Use diagonal pivot? */
144 if ( diag >= 0 ) { /* diagonal exists */
145 rtemp = fabs (lu_col_ptr[diag]);
146 if ( rtemp != 0.0 && rtemp >= thresh ) pivptr = diag;
147 }
148 *pivrow = lsub_ptr[pivptr];
149 }
150
151 /* Record pivot row */
152 perm_r[*pivrow] = jcol;
153
154 /* Interchange row subscripts */
155 if ( pivptr != nsupc ) {
156 itemp = lsub_ptr[pivptr];
157 lsub_ptr[pivptr] = lsub_ptr[nsupc];
158 lsub_ptr[nsupc] = itemp;
159
160 /* Interchange numerical values as well, for the whole snode, such
161 * that L is indexed the same way as A.
162 */
163 for (icol = 0; icol <= nsupc; icol++) {
164 itemp = pivptr + icol * nsupr;
165 temp = lu_sup_ptr[itemp];
166 lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
167 lu_sup_ptr[nsupc + icol*nsupr] = temp;
168 }
169 } /* if */
170
171 /* cdiv operation */
172 ops[FACT] += nsupr - nsupc;
173
174 temp = 1.0 / lu_col_ptr[nsupc];
175 for (k = nsupc+1; k < nsupr; k++)
176 lu_col_ptr[k] *= temp;
177
178 return 0;
179 }
180
181