1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
/*! \file
Copyright (c) 2003, The Regents of the University of California, through
Lawrence Berkeley National Laboratory (subject to receipt of any required
approvals from U.S. Dept. of Energy)
All rights reserved.
The source code is distributed under BSD license, see the file License.txt
at the top-level directory.
*/
/*! @file dpivotgrowth.c
* \brief Computes the reciprocal pivot growth factor
*
* <pre>
* -- SuperLU routine (version 2.0) --
* Univ. of California Berkeley, Xerox Palo Alto Research Center,
* and Lawrence Berkeley National Lab.
* November 15, 1997
* </pre>
*/
#include <math.h>
#include "slu_ddefs.h"
/*! \brief
*
* <pre>
* Purpose
* =======
*
* Compute the reciprocal pivot growth factor of the leading ncols columns
* of the matrix, using the formula:
* min_j ( max_i(abs(A_ij)) / max_i(abs(U_ij)) )
*
* Arguments
* =========
*
* ncols (input) int
* The number of columns of matrices A, L and U.
*
* A (input) SuperMatrix*
* Original matrix A, permuted by columns, of dimension
* (A->nrow, A->ncol). The type of A can be:
* Stype = NC; Dtype = SLU_D; Mtype = GE.
*
* L (output) SuperMatrix*
* The factor L from the factorization Pr*A=L*U; use compressed row
* subscripts storage for supernodes, i.e., L has type:
* Stype = SC; Dtype = SLU_D; Mtype = TRLU.
*
* U (output) SuperMatrix*
* The factor U from the factorization Pr*A*Pc=L*U. Use column-wise
* storage scheme, i.e., U has types: Stype = NC;
* Dtype = SLU_D; Mtype = TRU.
* </pre>
*/
double
dPivotGrowth(int ncols, SuperMatrix *A, int *perm_c,
SuperMatrix *L, SuperMatrix *U)
{
NCformat *Astore;
SCformat *Lstore;
NCformat *Ustore;
double *Aval, *Lval, *Uval;
int fsupc, nsupr, luptr, nz_in_U;
int i, j, k, oldcol;
int *inv_perm_c;
double rpg, maxaj, maxuj;
double smlnum;
double *luval;
/* Get machine constants. */
smlnum = dmach("S");
rpg = 1. / smlnum;
Astore = A->Store;
Lstore = L->Store;
Ustore = U->Store;
Aval = Astore->nzval;
Lval = Lstore->nzval;
Uval = Ustore->nzval;
inv_perm_c = (int *) SUPERLU_MALLOC(A->ncol*sizeof(int));
for (j = 0; j < A->ncol; ++j) inv_perm_c[perm_c[j]] = j;
for (k = 0; k <= Lstore->nsuper; ++k) {
fsupc = L_FST_SUPC(k);
nsupr = L_SUB_START(fsupc+1) - L_SUB_START(fsupc);
luptr = L_NZ_START(fsupc);
luval = &Lval[luptr];
nz_in_U = 1;
for (j = fsupc; j < L_FST_SUPC(k+1) && j < ncols; ++j) {
maxaj = 0.;
oldcol = inv_perm_c[j];
for (i = Astore->colptr[oldcol]; i < Astore->colptr[oldcol+1]; ++i)
maxaj = SUPERLU_MAX( maxaj, fabs(Aval[i]) );
maxuj = 0.;
for (i = Ustore->colptr[j]; i < Ustore->colptr[j+1]; i++)
maxuj = SUPERLU_MAX( maxuj, fabs(Uval[i]) );
/* Supernode */
for (i = 0; i < nz_in_U; ++i)
maxuj = SUPERLU_MAX( maxuj, fabs(luval[i]) );
++nz_in_U;
luval += nsupr;
if ( maxuj == 0. )
rpg = SUPERLU_MIN( rpg, 1.);
else
rpg = SUPERLU_MIN( rpg, maxaj / maxuj );
}
if ( j >= ncols ) break;
}
SUPERLU_FREE(inv_perm_c);
return (rpg);
}
|
