|
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 * File name: zgscon.c |
|
|
12 * History: Modified from lapack routines ZGECON. |
|
|
13 */ |
|
|
14 #include <math.h> |
|
|
15 #include "util.h" |
|
|
16 #include "zsp_defs.h" |
|
|
17 |
|
|
18 void |
|
|
19 zgscon(char *norm, SuperMatrix *L, SuperMatrix *U, |
|
|
20 double anorm, double *rcond, int *info) |
|
|
21 { |
|
|
22 /* |
|
|
23 Purpose |
|
|
24 ======= |
|
|
25 |
|
|
26 ZGSCON estimates the reciprocal of the condition number of a general |
|
|
27 real matrix A, in either the 1-norm or the infinity-norm, using |
|
|
28 the LU factorization computed by ZGETRF. |
|
|
29 |
|
|
30 An estimate is obtained for norm(inv(A)), and the reciprocal of the |
|
|
31 condition number is computed as |
|
|
32 RCOND = 1 / ( norm(A) * norm(inv(A)) ). |
|
|
33 |
|
|
34 See supermatrix.h for the definition of 'SuperMatrix' structure. |
|
|
35 |
|
|
36 Arguments |
|
|
37 ========= |
|
|
38 |
|
|
39 NORM (input) char* |
|
|
40 Specifies whether the 1-norm condition number or the |
|
|
41 infinity-norm condition number is required: |
|
|
42 = '1' or 'O': 1-norm; |
|
|
43 = 'I': Infinity-norm. |
|
|
44 |
|
|
45 L (input) SuperMatrix* |
|
|
46 The factor L from the factorization Pr*A*Pc=L*U as computed by |
|
|
47 zgstrf(). Use compressed row subscripts storage for supernodes, |
|
|
48 i.e., L has types: Stype = SC, Dtype = _Z, Mtype = TRLU. |
|
|
49 |
|
|
50 U (input) SuperMatrix* |
|
|
51 The factor U from the factorization Pr*A*Pc=L*U as computed by |
|
|
52 zgstrf(). Use column-wise storage scheme, i.e., U has types: |
|
|
53 Stype = NC, Dtype = _Z, Mtype = TRU. |
|
|
54 |
|
|
55 ANORM (input) double |
|
|
56 If NORM = '1' or 'O', the 1-norm of the original matrix A. |
|
|
57 If NORM = 'I', the infinity-norm of the original matrix A. |
|
|
58 |
|
|
59 RCOND (output) double* |
|
|
60 The reciprocal of the condition number of the matrix A, |
|
|
61 computed as RCOND = 1/(norm(A) * norm(inv(A))). |
|
|
62 |
|
|
63 INFO (output) int* |
|
|
64 = 0: successful exit |
|
|
65 < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
66 |
|
|
67 ===================================================================== |
|
|
68 */ |
|
|
69 |
|
|
70 /* Local variables */ |
|
|
71 int kase, kase1, onenrm, i; |
|
|
72 double ainvnm; |
|
|
73 doublecomplex *work; |
|
|
74 extern int zrscl_(int *, doublecomplex *, doublecomplex *, int *); |
|
|
75 |
|
|
76 extern int zlacon_(int *, doublecomplex *, doublecomplex *, double *, int *); |
|
|
77 |
|
|
78 |
|
|
79 /* Test the input parameters. */ |
|
|
80 *info = 0; |
|
|
81 onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O"); |
|
|
82 if (! onenrm && ! lsame_(norm, "I")) *info = -1; |
|
|
83 else if (L->nrow < 0 || L->nrow != L->ncol || |
|
|
84 L->Stype != SC || L->Dtype != _Z || L->Mtype != TRLU) |
|
|
85 *info = -2; |
|
|
86 else if (U->nrow < 0 || U->nrow != U->ncol || |
|
|
87 U->Stype != NC || U->Dtype != _Z || U->Mtype != TRU) |
|
|
88 *info = -3; |
|
|
89 if (*info != 0) { |
|
|
90 i = -(*info); |
|
|
91 xerbla_("zgscon", &i); |
|
|
92 return; |
|
|
93 } |
|
|
94 |
|
|
95 /* Quick return if possible */ |
|
|
96 *rcond = 0.; |
|
|
97 if ( L->nrow == 0 || U->nrow == 0) { |
|
|
98 *rcond = 1.; |
|
|
99 return; |
|
|
100 } |
|
|
101 |
|
|
102 work = doublecomplexCalloc( 3*L->nrow ); |
|
|
103 |
|
|
104 |
|
|
105 if ( !work ) |
|
|
106 ABORT("Malloc fails for work arrays in zgscon."); |
|
|
107 |
|
|
108 /* Estimate the norm of inv(A). */ |
|
|
109 ainvnm = 0.; |
|
|
110 if ( onenrm ) kase1 = 1; |
|
|
111 else kase1 = 2; |
|
|
112 kase = 0; |
|
|
113 |
|
|
114 do { |
|
|
115 zlacon_(&L->nrow, &work[L->nrow], &work[0], &ainvnm, &kase); |
|
|
116 |
|
|
117 if (kase == 0) break; |
|
|
118 |
|
|
119 if (kase == kase1) { |
|
|
120 /* Multiply by inv(L). */ |
|
|
121 sp_ztrsv("Lower", "No transpose", "Unit", L, U, &work[0], info); |
|
|
122 |
|
|
123 /* Multiply by inv(U). */ |
|
|
124 sp_ztrsv("Upper", "No transpose", "Non-unit", L, U, &work[0],info); |
|
|
125 |
|
|
126 } else { |
|
|
127 |
|
|
128 /* Multiply by inv(U'). */ |
|
|
129 sp_ztrsv("Upper", "Transpose", "Non-unit", L, U, &work[0], info); |
|
|
130 |
|
|
131 /* Multiply by inv(L'). */ |
|
|
132 sp_ztrsv("Lower", "Transpose", "Unit", L, U, &work[0], info); |
|
|
133 |
|
|
134 } |
|
|
135 |
|
|
136 } while ( kase != 0 ); |
|
|
137 |
|
|
138 /* Compute the estimate of the reciprocal condition number. */ |
|
|
139 if (ainvnm != 0.) *rcond = (1. / ainvnm) / anorm; |
|
|
140 |
|
|
141 SUPERLU_FREE (work); |
|
|
142 return; |
|
|
143 |
|
|
144 } /* zgscon */ |
|
|
145 |
