1
2 /*
3 * -- SuperLU routine (version 2.0) --
4 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
5 * and Lawrence Berkeley National Lab.
6 * November 15, 1997
7 *
8 */
9 #include <math.h>
10 #include "slu_zdefs.h"
11
zgst01(int m,int n,SuperMatrix * A,SuperMatrix * L,SuperMatrix * U,int * perm_c,int * perm_r,double * resid)12 int zgst01(int m, int n, SuperMatrix *A, SuperMatrix *L,
13 SuperMatrix *U, int *perm_c, int *perm_r, double *resid)
14 {
15 /*
16 Purpose
17 =======
18
19 ZGST01 reconstructs a matrix A from its L*U factorization and
20 computes the residual
21 norm(L*U - A) / ( N * norm(A) * EPS ),
22 where EPS is the machine epsilon.
23
24 Arguments
25 ==========
26
27 M (input) INT
28 The number of rows of the matrix A. M >= 0.
29
30 N (input) INT
31 The number of columns of the matrix A. N >= 0.
32
33 A (input) SuperMatrix *, dimension (A->nrow, A->ncol)
34 The original M x N matrix A.
35
36 L (input) SuperMatrix *, dimension (L->nrow, L->ncol)
37 The factor matrix L.
38
39 U (input) SuperMatrix *, dimension (U->nrow, U->ncol)
40 The factor matrix U.
41
42 perm_c (input) INT array, dimension (N)
43 The column permutation from ZGSTRF.
44
45 perm_r (input) INT array, dimension (M)
46 The pivot indices from ZGSTRF.
47
48 RESID (output) DOUBLE*
49 norm(L*U - A) / ( N * norm(A) * EPS )
50
51 =====================================================================
52 */
53
54 /* Local variables */
55 doublecomplex zero = {0.0, 0.0};
56 int i, j, k, arow, lptr,isub, urow, superno, fsupc, u_part;
57 doublecomplex utemp, comp_temp;
58 double anorm, tnorm, cnorm;
59 double eps;
60 doublecomplex *work;
61 SCformat *Lstore;
62 NCformat *Astore, *Ustore;
63 doublecomplex *Aval, *Lval, *Uval;
64 int *colbeg, *colend;
65
66 /* Function prototypes */
67 extern double zlangs(char *, SuperMatrix *);
68
69 /* Quick exit if M = 0 or N = 0. */
70
71 if (m <= 0 || n <= 0) {
72 *resid = 0.f;
73 return 0;
74 }
75
76 work = (doublecomplex *)doublecomplexCalloc(m);
77
78 Astore = A->Store;
79 Aval = Astore->nzval;
80 Lstore = L->Store;
81 Lval = Lstore->nzval;
82 Ustore = U->Store;
83 Uval = Ustore->nzval;
84
85 /* Determine EPS and the norm of A. */
86 eps = dlamch_("Epsilon");
87 anorm = zlangs("1", A);
88 cnorm = 0.;
89
90 /* Compute the product L*U, one column at a time */
91 for (k = 0; k < n; ++k) {
92
93 /* The U part outside the rectangular supernode */
94 for (i = U_NZ_START(k); i < U_NZ_START(k+1); ++i) {
95 urow = U_SUB(i);
96 utemp = Uval[i];
97 superno = Lstore->col_to_sup[urow];
98 fsupc = L_FST_SUPC(superno);
99 u_part = urow - fsupc + 1;
100 lptr = L_SUB_START(fsupc) + u_part;
101 work[L_SUB(lptr-1)].r -= utemp.r;
102 work[L_SUB(lptr-1)].i -= utemp.i;
103 for (j = L_NZ_START(urow) + u_part; j < L_NZ_START(urow+1); ++j) {
104 isub = L_SUB(lptr);
105 zz_mult(&comp_temp, &utemp, &Lval[j]);
106 z_sub(&work[isub], &work[isub], &comp_temp);
107 ++lptr;
108 }
109 }
110
111 /* The U part inside the rectangular supernode */
112 superno = Lstore->col_to_sup[k];
113 fsupc = L_FST_SUPC(superno);
114 urow = L_NZ_START(k);
115 for (i = fsupc; i <= k; ++i) {
116 utemp = Lval[urow++];
117 u_part = i - fsupc + 1;
118 lptr = L_SUB_START(fsupc) + u_part;
119 work[L_SUB(lptr-1)].r -= utemp.r;
120 work[L_SUB(lptr-1)].i -= utemp.i;
121 for (j = L_NZ_START(i)+u_part; j < L_NZ_START(i+1); ++j) {
122 isub = L_SUB(lptr);
123 zz_mult(&comp_temp, &utemp, &Lval[j]);
124 z_sub(&work[isub], &work[isub], &comp_temp);
125 ++lptr;
126 }
127 }
128
129 /* Now compute A[k] - (L*U)[k] (Both matrices may be permuted.) */
130
131 colbeg = intMalloc(n);
132 colend = intMalloc(n);
133 for (i = 0; i < n; i++) {
134 colbeg[perm_c[i]] = Astore->colptr[i];
135 colend[perm_c[i]] = Astore->colptr[i+1];
136 }
137
138 for (i = colbeg[k]; i < colend[k]; ++i) {
139 arow = Astore->rowind[i];
140 work[perm_r[arow]].r += Aval[i].r;
141 work[perm_r[arow]].i += Aval[i].i;
142 }
143
144 /* Now compute the 1-norm of the column vector work */
145 tnorm = 0.;
146 for (i = 0; i < m; ++i) {
147 tnorm += fabs(work[i].r) + fabs(work[i].i);
148 work[i] = zero;
149 }
150 cnorm = SUPERLU_MAX(tnorm, cnorm);
151 }
152
153 *resid = cnorm;
154
155 if (anorm <= 0.f) {
156 if (*resid != 0.f) {
157 *resid = 1.f / eps;
158 }
159 } else {
160 *resid = *resid / (float) n / anorm / eps;
161 }
162
163 SUPERLU_FREE(work);
164 SUPERLU_FREE(colbeg);
165 SUPERLU_FREE(colend);
166 return 0;
167
168 /* End of ZGST01 */
169
170 } /* zgst01_ */
171
172