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