Mercurial > forge

comparison main/sparse/SuperLU/SRC/dpivotL.c @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial revision
author pkienzle
date Wed, 10 Oct 2001 19:54:49 +0000
parents
children 785c84b30f53
comparison
equal deleted inserted replaced
-1:000000000000 0:6b33357c7561
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