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