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 /*! @file dgsrfs.c
13 * \brief Improves computed solution to a system of inear equations
14 *
15 * <pre>
16 * -- SuperLU routine (version 5.1) --
17 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
18 * and Lawrence Berkeley National Lab.
19 * October 15, 2003
20 *
21 * Modified from lapack routine DGERFS
22 * Last modified: December 3, 2015
23 * </pre>
24 */
25 /*
26 * File name: dgsrfs.c
27 * History: Modified from lapack routine DGERFS
28 */
29 #include <math.h>
30 #include "slu_ddefs.h"
31
32 /*! \brief
33 *
34 * <pre>
35 * Purpose
36 * =======
37 *
38 * DGSRFS improves the computed solution to a system of linear
39 * equations and provides error bounds and backward error estimates for
40 * the solution.
41 *
42 * If equilibration was performed, the system becomes:
43 * (diag(R)*A_original*diag(C)) * X = diag(R)*B_original.
44 *
45 * See supermatrix.h for the definition of 'SuperMatrix' structure.
46 *
47 * Arguments
48 * =========
49 *
50 * trans (input) trans_t
51 * Specifies the form of the system of equations:
52 * = NOTRANS: A * X = B (No transpose)
53 * = TRANS: A'* X = B (Transpose)
54 * = CONJ: A**H * X = B (Conjugate transpose)
55 *
56 * A (input) SuperMatrix*
57 * The original matrix A in the system, or the scaled A if
58 * equilibration was done. The type of A can be:
59 * Stype = SLU_NC, Dtype = SLU_D, Mtype = SLU_GE.
60 *
61 * L (input) SuperMatrix*
62 * The factor L from the factorization Pr*A*Pc=L*U. Use
63 * compressed row subscripts storage for supernodes,
64 * i.e., L has types: Stype = SLU_SC, Dtype = SLU_D, Mtype = SLU_TRLU.
65 *
66 * U (input) SuperMatrix*
67 * The factor U from the factorization Pr*A*Pc=L*U as computed by
68 * dgstrf(). Use column-wise storage scheme,
69 * i.e., U has types: Stype = SLU_NC, Dtype = SLU_D, Mtype = SLU_TRU.
70 *
71 * perm_c (input) int*, dimension (A->ncol)
72 * Column permutation vector, which defines the
73 * permutation matrix Pc; perm_c[i] = j means column i of A is
74 * in position j in A*Pc.
75 *
76 * perm_r (input) int*, dimension (A->nrow)
77 * Row permutation vector, which defines the permutation matrix Pr;
78 * perm_r[i] = j means row i of A is in position j in Pr*A.
79 *
80 * equed (input) Specifies the form of equilibration that was done.
81 * = 'N': No equilibration.
82 * = 'R': Row equilibration, i.e., A was premultiplied by diag(R).
83 * = 'C': Column equilibration, i.e., A was postmultiplied by
84 * diag(C).
85 * = 'B': Both row and column equilibration, i.e., A was replaced
86 * by diag(R)*A*diag(C).
87 *
88 * R (input) double*, dimension (A->nrow)
89 * The row scale factors for A.
90 * If equed = 'R' or 'B', A is premultiplied by diag(R).
91 * If equed = 'N' or 'C', R is not accessed.
92 *
93 * C (input) double*, dimension (A->ncol)
94 * The column scale factors for A.
95 * If equed = 'C' or 'B', A is postmultiplied by diag(C).
96 * If equed = 'N' or 'R', C is not accessed.
97 *
98 * B (input) SuperMatrix*
99 * B has types: Stype = SLU_DN, Dtype = SLU_D, Mtype = SLU_GE.
100 * The right hand side matrix B.
101 * if equed = 'R' or 'B', B is premultiplied by diag(R).
102 *
103 * X (input/output) SuperMatrix*
104 * X has types: Stype = SLU_DN, Dtype = SLU_D, Mtype = SLU_GE.
105 * On entry, the solution matrix X, as computed by dgstrs().
106 * On exit, the improved solution matrix X.
107 * if *equed = 'C' or 'B', X should be premultiplied by diag(C)
108 * in order to obtain the solution to the original system.
109 *
110 * FERR (output) double*, dimension (B->ncol)
111 * The estimated forward error bound for each solution vector
112 * X(j) (the j-th column of the solution matrix X).
113 * If XTRUE is the true solution corresponding to X(j), FERR(j)
114 * is an estimated upper bound for the magnitude of the largest
115 * element in (X(j) - XTRUE) divided by the magnitude of the
116 * largest element in X(j). The estimate is as reliable as
117 * the estimate for RCOND, and is almost always a slight
118 * overestimate of the true error.
119 *
120 * BERR (output) double*, dimension (B->ncol)
121 * The componentwise relative backward error of each solution
122 * vector X(j) (i.e., the smallest relative change in
123 * any element of A or B that makes X(j) an exact solution).
124 *
125 * stat (output) SuperLUStat_t*
126 * Record the statistics on runtime and floating-point operation count.
127 * See util.h for the definition of 'SuperLUStat_t'.
128 *
129 * info (output) int*
130 * = 0: successful exit
131 * < 0: if INFO = -i, the i-th argument had an illegal value
132 *
133 * Internal Parameters
134 * ===================
135 *
136 * ITMAX is the maximum number of steps of iterative refinement.
137 *
138 * </pre>
139 */
140 void
dgsrfs(trans_t trans,SuperMatrix * A,SuperMatrix * L,SuperMatrix * U,int * perm_c,int * perm_r,char * equed,double * R,double * C,SuperMatrix * B,SuperMatrix * X,double * ferr,double * berr,SuperLUStat_t * stat,int * info)141 dgsrfs(trans_t trans, SuperMatrix *A, SuperMatrix *L, SuperMatrix *U,
142 int *perm_c, int *perm_r, char *equed, double *R, double *C,
143 SuperMatrix *B, SuperMatrix *X, double *ferr, double *berr,
144 SuperLUStat_t *stat, int *info)
145 {
146
147
148 #define ITMAX 5
149
150 /* Table of constant values */
151 int ione = 1;
152 double ndone = -1.;
153 double done = 1.;
154
155 /* Local variables */
156 NCformat *Astore;
157 double *Aval;
158 SuperMatrix Bjcol;
159 DNformat *Bstore, *Xstore, *Bjcol_store;
160 double *Bmat, *Xmat, *Bptr, *Xptr;
161 int kase;
162 double safe1, safe2;
163 int i, j, k, irow, nz, count, notran, rowequ, colequ;
164 int ldb, ldx, nrhs;
165 double s, xk, lstres, eps, safmin;
166 char transc[1];
167 trans_t transt;
168 double *work;
169 double *rwork;
170 int *iwork;
171 int isave[3];
172
173 extern int dlacon2_(int *, double *, double *, int *, double *, int *, int []);
174 #ifdef _CRAY
175 extern int SCOPY(int *, double *, int *, double *, int *);
176 extern int SSAXPY(int *, double *, double *, int *, double *, int *);
177 #else
178 extern int dcopy_(int *, double *, int *, double *, int *);
179 extern int daxpy_(int *, double *, double *, int *, double *, int *);
180 #endif
181
182 Astore = A->Store;
183 Aval = Astore->nzval;
184 Bstore = B->Store;
185 Xstore = X->Store;
186 Bmat = Bstore->nzval;
187 Xmat = Xstore->nzval;
188 ldb = Bstore->lda;
189 ldx = Xstore->lda;
190 nrhs = B->ncol;
191
192 /* Test the input parameters */
193 *info = 0;
194 notran = (trans == NOTRANS);
195 if ( !notran && trans != TRANS && trans != CONJ ) *info = -1;
196 else if ( A->nrow != A->ncol || A->nrow < 0 ||
197 A->Stype != SLU_NC || A->Dtype != SLU_D || A->Mtype != SLU_GE )
198 *info = -2;
199 else if ( L->nrow != L->ncol || L->nrow < 0 ||
200 L->Stype != SLU_SC || L->Dtype != SLU_D || L->Mtype != SLU_TRLU )
201 *info = -3;
202 else if ( U->nrow != U->ncol || U->nrow < 0 ||
203 U->Stype != SLU_NC || U->Dtype != SLU_D || U->Mtype != SLU_TRU )
204 *info = -4;
205 else if ( ldb < SUPERLU_MAX(0, A->nrow) ||
206 B->Stype != SLU_DN || B->Dtype != SLU_D || B->Mtype != SLU_GE )
207 *info = -10;
208 else if ( ldx < SUPERLU_MAX(0, A->nrow) ||
209 X->Stype != SLU_DN || X->Dtype != SLU_D || X->Mtype != SLU_GE )
210 *info = -11;
211 if (*info != 0) {
212 i = -(*info);
213 input_error("dgsrfs", &i);
214 return;
215 }
216
217 /* Quick return if possible */
218 if ( A->nrow == 0 || nrhs == 0) {
219 for (j = 0; j < nrhs; ++j) {
220 ferr[j] = 0.;
221 berr[j] = 0.;
222 }
223 return;
224 }
225
226 rowequ = strncmp(equed, "R", 1)==0 || strncmp(equed, "B", 1)==0;
227 colequ = strncmp(equed, "C", 1)==0 || strncmp(equed, "B", 1)==0;
228
229 /* Allocate working space */
230 work = doubleMalloc(2*A->nrow);
231 rwork = (double *) SUPERLU_MALLOC( A->nrow * sizeof(double) );
232 iwork = intMalloc(2*A->nrow);
233 if ( !work || !rwork || !iwork )
234 ABORT("Malloc fails for work/rwork/iwork.");
235
236 if ( notran ) {
237 *(unsigned char *)transc = 'N';
238 transt = TRANS;
239 } else if ( trans == TRANS ) {
240 *(unsigned char *)transc = 'T';
241 transt = NOTRANS;
242 } else if ( trans == CONJ ) {
243 *(unsigned char *)transc = 'C';
244 transt = NOTRANS;
245 }
246
247 /* NZ = maximum number of nonzero elements in each row of A, plus 1 */
248 nz = A->ncol + 1;
249 eps = dmach("Epsilon");
250 safmin = dmach("Safe minimum");
251
252 /* Set SAFE1 essentially to be the underflow threshold times the
253 number of additions in each row. */
254 safe1 = nz * safmin;
255 safe2 = safe1 / eps;
256
257 /* Compute the number of nonzeros in each row (or column) of A */
258 for (i = 0; i < A->nrow; ++i) iwork[i] = 0;
259 if ( notran ) {
260 for (k = 0; k < A->ncol; ++k)
261 for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
262 ++iwork[Astore->rowind[i]];
263 } else {
264 for (k = 0; k < A->ncol; ++k)
265 iwork[k] = Astore->colptr[k+1] - Astore->colptr[k];
266 }
267
268 /* Copy one column of RHS B into Bjcol. */
269 Bjcol.Stype = B->Stype;
270 Bjcol.Dtype = B->Dtype;
271 Bjcol.Mtype = B->Mtype;
272 Bjcol.nrow = B->nrow;
273 Bjcol.ncol = 1;
274 Bjcol.Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
275 if ( !Bjcol.Store ) ABORT("SUPERLU_MALLOC fails for Bjcol.Store");
276 Bjcol_store = Bjcol.Store;
277 Bjcol_store->lda = ldb;
278 Bjcol_store->nzval = work; /* address aliasing */
279
280 /* Do for each right hand side ... */
281 for (j = 0; j < nrhs; ++j) {
282 count = 0;
283 lstres = 3.;
284 Bptr = &Bmat[j*ldb];
285 Xptr = &Xmat[j*ldx];
286
287 while (1) { /* Loop until stopping criterion is satisfied. */
288
289 /* Compute residual R = B - op(A) * X,
290 where op(A) = A, A**T, or A**H, depending on TRANS. */
291
292 #ifdef _CRAY
293 SCOPY(&A->nrow, Bptr, &ione, work, &ione);
294 #else
295 dcopy_(&A->nrow, Bptr, &ione, work, &ione);
296 #endif
297 sp_dgemv(transc, ndone, A, Xptr, ione, done, work, ione);
298
299 /* Compute componentwise relative backward error from formula
300 max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
301 where abs(Z) is the componentwise absolute value of the matrix
302 or vector Z. If the i-th component of the denominator is less
303 than SAFE2, then SAFE1 is added to the i-th component of the
304 numerator before dividing. */
305
306 for (i = 0; i < A->nrow; ++i) rwork[i] = fabs( Bptr[i] );
307
308 /* Compute abs(op(A))*abs(X) + abs(B). */
309 if ( notran ) {
310 for (k = 0; k < A->ncol; ++k) {
311 xk = fabs( Xptr[k] );
312 for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
313 rwork[Astore->rowind[i]] += fabs(Aval[i]) * xk;
314 }
315 } else { /* trans = TRANS or CONJ */
316 for (k = 0; k < A->ncol; ++k) {
317 s = 0.;
318 for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) {
319 irow = Astore->rowind[i];
320 s += fabs(Aval[i]) * fabs(Xptr[irow]);
321 }
322 rwork[k] += s;
323 }
324 }
325 s = 0.;
326 for (i = 0; i < A->nrow; ++i) {
327 if (rwork[i] > safe2) {
328 s = SUPERLU_MAX( s, fabs(work[i]) / rwork[i] );
329 } else if ( rwork[i] != 0.0 ) {
330 /* Adding SAFE1 to the numerator guards against
331 spuriously zero residuals (underflow). */
332 s = SUPERLU_MAX( s, (safe1 + fabs(work[i])) / rwork[i] );
333 }
334 /* If rwork[i] is exactly 0.0, then we know the true
335 residual also must be exactly 0.0. */
336 }
337 berr[j] = s;
338
339 /* Test stopping criterion. Continue iterating if
340 1) The residual BERR(J) is larger than machine epsilon, and
341 2) BERR(J) decreased by at least a factor of 2 during the
342 last iteration, and
343 3) At most ITMAX iterations tried. */
344
345 if (berr[j] > eps && berr[j] * 2. <= lstres && count < ITMAX) {
346 /* Update solution and try again. */
347 dgstrs (trans, L, U, perm_c, perm_r, &Bjcol, stat, info);
348
349 #ifdef _CRAY
350 SAXPY(&A->nrow, &done, work, &ione,
351 &Xmat[j*ldx], &ione);
352 #else
353 daxpy_(&A->nrow, &done, work, &ione,
354 &Xmat[j*ldx], &ione);
355 #endif
356 lstres = berr[j];
357 ++count;
358 } else {
359 break;
360 }
361
362 } /* end while */
363
364 stat->RefineSteps = count;
365
366 /* Bound error from formula:
367 norm(X - XTRUE) / norm(X) .le. FERR = norm( abs(inv(op(A)))*
368 ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
369 where
370 norm(Z) is the magnitude of the largest component of Z
371 inv(op(A)) is the inverse of op(A)
372 abs(Z) is the componentwise absolute value of the matrix or
373 vector Z
374 NZ is the maximum number of nonzeros in any row of A, plus 1
375 EPS is machine epsilon
376
377 The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
378 is incremented by SAFE1 if the i-th component of
379 abs(op(A))*abs(X) + abs(B) is less than SAFE2.
380
381 Use DLACON2 to estimate the infinity-norm of the matrix
382 inv(op(A)) * diag(W),
383 where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) */
384
385 for (i = 0; i < A->nrow; ++i) rwork[i] = fabs( Bptr[i] );
386
387 /* Compute abs(op(A))*abs(X) + abs(B). */
388 if ( notran ) {
389 for (k = 0; k < A->ncol; ++k) {
390 xk = fabs( Xptr[k] );
391 for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
392 rwork[Astore->rowind[i]] += fabs(Aval[i]) * xk;
393 }
394 } else { /* trans == TRANS or CONJ */
395 for (k = 0; k < A->ncol; ++k) {
396 s = 0.;
397 for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) {
398 irow = Astore->rowind[i];
399 xk = fabs( Xptr[irow] );
400 s += fabs(Aval[i]) * xk;
401 }
402 rwork[k] += s;
403 }
404 }
405
406 for (i = 0; i < A->nrow; ++i)
407 if (rwork[i] > safe2)
408 rwork[i] = fabs(work[i]) + (iwork[i]+1)*eps*rwork[i];
409 else
410 rwork[i] = fabs(work[i])+(iwork[i]+1)*eps*rwork[i]+safe1;
411
412 kase = 0;
413
414 do {
415 dlacon2_(&A->nrow, &work[A->nrow], work,
416 &iwork[A->nrow], &ferr[j], &kase, isave);
417 if (kase == 0) break;
418
419 if (kase == 1) {
420 /* Multiply by diag(W)*inv(op(A)**T)*(diag(C) or diag(R)). */
421 if ( notran && colequ )
422 for (i = 0; i < A->ncol; ++i) work[i] *= C[i];
423 else if ( !notran && rowequ )
424 for (i = 0; i < A->nrow; ++i) work[i] *= R[i];
425
426 dgstrs (transt, L, U, perm_c, perm_r, &Bjcol, stat, info);
427
428 for (i = 0; i < A->nrow; ++i) work[i] *= rwork[i];
429 } else {
430 /* Multiply by (diag(C) or diag(R))*inv(op(A))*diag(W). */
431 for (i = 0; i < A->nrow; ++i) work[i] *= rwork[i];
432
433 dgstrs (trans, L, U, perm_c, perm_r, &Bjcol, stat, info);
434
435 if ( notran && colequ )
436 for (i = 0; i < A->ncol; ++i) work[i] *= C[i];
437 else if ( !notran && rowequ )
438 for (i = 0; i < A->ncol; ++i) work[i] *= R[i];
439 }
440
441 } while ( kase != 0 );
442
443
444 /* Normalize error. */
445 lstres = 0.;
446 if ( notran && colequ ) {
447 for (i = 0; i < A->nrow; ++i)
448 lstres = SUPERLU_MAX( lstres, C[i] * fabs( Xptr[i]) );
449 } else if ( !notran && rowequ ) {
450 for (i = 0; i < A->nrow; ++i)
451 lstres = SUPERLU_MAX( lstres, R[i] * fabs( Xptr[i]) );
452 } else {
453 for (i = 0; i < A->nrow; ++i)
454 lstres = SUPERLU_MAX( lstres, fabs( Xptr[i]) );
455 }
456 if ( lstres != 0. )
457 ferr[j] /= lstres;
458
459 } /* for each RHS j ... */
460
461 SUPERLU_FREE(work);
462 SUPERLU_FREE(rwork);
463 SUPERLU_FREE(iwork);
464 SUPERLU_FREE(Bjcol.Store);
465
466 return;
467
468 } /* dgsrfs */
469