kusano 7d535a
kusano 7d535a
/*! @file sgsrfs.c
kusano 7d535a
 * \brief Improves computed solution to a system of inear equations
kusano 7d535a
 * 
kusano 7d535a
 * 
kusano 7d535a
 * -- SuperLU routine (version 3.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * October 15, 2003
kusano 7d535a
 *
kusano 7d535a
 * Modified from lapack routine SGERFS
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
/*
kusano 7d535a
 * File name:	sgsrfs.c
kusano 7d535a
 * History:     Modified from lapack routine SGERFS
kusano 7d535a
 */
kusano 7d535a
#include <math.h></math.h>
kusano 7d535a
#include "slu_sdefs.h"
kusano 7d535a
kusano 7d535a
/*! \brief
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 *   Purpose   
kusano 7d535a
 *   =======   
kusano 7d535a
 *
kusano 7d535a
 *   SGSRFS improves the computed solution to a system of linear   
kusano 7d535a
 *   equations and provides error bounds and backward error estimates for 
kusano 7d535a
 *   the solution.   
kusano 7d535a
 *
kusano 7d535a
 *   If equilibration was performed, the system becomes:
kusano 7d535a
 *           (diag(R)*A_original*diag(C)) * X = diag(R)*B_original.
kusano 7d535a
 *
kusano 7d535a
 *   See supermatrix.h for the definition of 'SuperMatrix' structure.
kusano 7d535a
 *
kusano 7d535a
 *   Arguments   
kusano 7d535a
 *   =========   
kusano 7d535a
 *
kusano 7d535a
 * trans   (input) trans_t
kusano 7d535a
 *          Specifies the form of the system of equations:
kusano 7d535a
 *          = NOTRANS: A * X = B  (No transpose)
kusano 7d535a
 *          = TRANS:   A'* X = B  (Transpose)
kusano 7d535a
 *          = CONJ:    A**H * X = B  (Conjugate transpose)
kusano 7d535a
 *   
kusano 7d535a
 *   A       (input) SuperMatrix*
kusano 7d535a
 *           The original matrix A in the system, or the scaled A if
kusano 7d535a
 *           equilibration was done. The type of A can be:
kusano 7d535a
 *           Stype = SLU_NC, Dtype = SLU_S, Mtype = SLU_GE.
kusano 7d535a
 *    
kusano 7d535a
 *   L       (input) SuperMatrix*
kusano 7d535a
 *	     The factor L from the factorization Pr*A*Pc=L*U. Use
kusano 7d535a
 *           compressed row subscripts storage for supernodes, 
kusano 7d535a
 *           i.e., L has types: Stype = SLU_SC, Dtype = SLU_S, Mtype = SLU_TRLU.
kusano 7d535a
 * 
kusano 7d535a
 *   U       (input) SuperMatrix*
kusano 7d535a
 *           The factor U from the factorization Pr*A*Pc=L*U as computed by
kusano 7d535a
 *           sgstrf(). Use column-wise storage scheme, 
kusano 7d535a
 *           i.e., U has types: Stype = SLU_NC, Dtype = SLU_S, Mtype = SLU_TRU.
kusano 7d535a
 *
kusano 7d535a
 *   perm_c  (input) int*, dimension (A->ncol)
kusano 7d535a
 *	     Column permutation vector, which defines the 
kusano 7d535a
 *           permutation matrix Pc; perm_c[i] = j means column i of A is 
kusano 7d535a
 *           in position j in A*Pc.
kusano 7d535a
 *
kusano 7d535a
 *   perm_r  (input) int*, dimension (A->nrow)
kusano 7d535a
 *           Row permutation vector, which defines the permutation matrix Pr;
kusano 7d535a
 *           perm_r[i] = j means row i of A is in position j in Pr*A.
kusano 7d535a
 *
kusano 7d535a
 *   equed   (input) Specifies the form of equilibration that was done.
kusano 7d535a
 *           = 'N': No equilibration.
kusano 7d535a
 *           = 'R': Row equilibration, i.e., A was premultiplied by diag(R).
kusano 7d535a
 *           = 'C': Column equilibration, i.e., A was postmultiplied by
kusano 7d535a
 *                  diag(C).
kusano 7d535a
 *           = 'B': Both row and column equilibration, i.e., A was replaced 
kusano 7d535a
 *                  by diag(R)*A*diag(C).
kusano 7d535a
 *
kusano 7d535a
 *   R       (input) float*, dimension (A->nrow)
kusano 7d535a
 *           The row scale factors for A.
kusano 7d535a
 *           If equed = 'R' or 'B', A is premultiplied by diag(R).
kusano 7d535a
 *           If equed = 'N' or 'C', R is not accessed.
kusano 7d535a
 * 
kusano 7d535a
 *   C       (input) float*, dimension (A->ncol)
kusano 7d535a
 *           The column scale factors for A.
kusano 7d535a
 *           If equed = 'C' or 'B', A is postmultiplied by diag(C).
kusano 7d535a
 *           If equed = 'N' or 'R', C is not accessed.
kusano 7d535a
 *
kusano 7d535a
 *   B       (input) SuperMatrix*
kusano 7d535a
 *           B has types: Stype = SLU_DN, Dtype = SLU_S, Mtype = SLU_GE.
kusano 7d535a
 *           The right hand side matrix B.
kusano 7d535a
 *           if equed = 'R' or 'B', B is premultiplied by diag(R).
kusano 7d535a
 *
kusano 7d535a
 *   X       (input/output) SuperMatrix*
kusano 7d535a
 *           X has types: Stype = SLU_DN, Dtype = SLU_S, Mtype = SLU_GE.
kusano 7d535a
 *           On entry, the solution matrix X, as computed by sgstrs().
kusano 7d535a
 *           On exit, the improved solution matrix X.
kusano 7d535a
 *           if *equed = 'C' or 'B', X should be premultiplied by diag(C)
kusano 7d535a
 *               in order to obtain the solution to the original system.
kusano 7d535a
 *
kusano 7d535a
 *   FERR    (output) float*, dimension (B->ncol)   
kusano 7d535a
 *           The estimated forward error bound for each solution vector   
kusano 7d535a
 *           X(j) (the j-th column of the solution matrix X).   
kusano 7d535a
 *           If XTRUE is the true solution corresponding to X(j), FERR(j) 
kusano 7d535a
 *           is an estimated upper bound for the magnitude of the largest 
kusano 7d535a
 *           element in (X(j) - XTRUE) divided by the magnitude of the   
kusano 7d535a
 *           largest element in X(j).  The estimate is as reliable as   
kusano 7d535a
 *           the estimate for RCOND, and is almost always a slight   
kusano 7d535a
 *           overestimate of the true error.
kusano 7d535a
 *
kusano 7d535a
 *   BERR    (output) float*, dimension (B->ncol)   
kusano 7d535a
 *           The componentwise relative backward error of each solution   
kusano 7d535a
 *           vector X(j) (i.e., the smallest relative change in   
kusano 7d535a
 *           any element of A or B that makes X(j) an exact solution).
kusano 7d535a
 *
kusano 7d535a
 *   stat     (output) SuperLUStat_t*
kusano 7d535a
 *            Record the statistics on runtime and floating-point operation count.
kusano 7d535a
 *            See util.h for the definition of 'SuperLUStat_t'.
kusano 7d535a
 *
kusano 7d535a
 *   info    (output) int*   
kusano 7d535a
 *           = 0:  successful exit   
kusano 7d535a
 *            < 0:  if INFO = -i, the i-th argument had an illegal value   
kusano 7d535a
 *
kusano 7d535a
 *    Internal Parameters   
kusano 7d535a
 *    ===================   
kusano 7d535a
 *
kusano 7d535a
 *    ITMAX is the maximum number of steps of iterative refinement.   
kusano 7d535a
 *
kusano 7d535a
 * 
kusano 7d535a
 */
kusano 7d535a
void
kusano 7d535a
sgsrfs(trans_t trans, SuperMatrix *A, SuperMatrix *L, SuperMatrix *U,
kusano 7d535a
       int *perm_c, int *perm_r, char *equed, float *R, float *C,
kusano 7d535a
       SuperMatrix *B, SuperMatrix *X, float *ferr, float *berr,
kusano 7d535a
       SuperLUStat_t *stat, int *info)
kusano 7d535a
{
kusano 7d535a
kusano 7d535a
kusano 7d535a
#define ITMAX 5
kusano 7d535a
    
kusano 7d535a
    /* Table of constant values */
kusano 7d535a
    int    ione = 1;
kusano 7d535a
    float ndone = -1.;
kusano 7d535a
    float done = 1.;
kusano 7d535a
    
kusano 7d535a
    /* Local variables */
kusano 7d535a
    NCformat *Astore;
kusano 7d535a
    float   *Aval;
kusano 7d535a
    SuperMatrix Bjcol;
kusano 7d535a
    DNformat *Bstore, *Xstore, *Bjcol_store;
kusano 7d535a
    float   *Bmat, *Xmat, *Bptr, *Xptr;
kusano 7d535a
    int      kase;
kusano 7d535a
    float   safe1, safe2;
kusano 7d535a
    int      i, j, k, irow, nz, count, notran, rowequ, colequ;
kusano 7d535a
    int      ldb, ldx, nrhs;
kusano 7d535a
    float   s, xk, lstres, eps, safmin;
kusano 7d535a
    char     transc[1];
kusano 7d535a
    trans_t  transt;
kusano 7d535a
    float   *work;
kusano 7d535a
    float   *rwork;
kusano 7d535a
    int      *iwork;
kusano 7d535a
kusano 7d535a
    extern int slacon_(int *, float *, float *, int *, float *, int *);
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
    extern int SCOPY(int *, float *, int *, float *, int *);
kusano 7d535a
    extern int SSAXPY(int *, float *, float *, int *, float *, int *);
kusano 7d535a
#else
kusano 7d535a
    extern int scopy_(int *, float *, int *, float *, int *);
kusano 7d535a
    extern int saxpy_(int *, float *, float *, int *, float *, int *);
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    Astore = A->Store;
kusano 7d535a
    Aval   = Astore->nzval;
kusano 7d535a
    Bstore = B->Store;
kusano 7d535a
    Xstore = X->Store;
kusano 7d535a
    Bmat   = Bstore->nzval;
kusano 7d535a
    Xmat   = Xstore->nzval;
kusano 7d535a
    ldb    = Bstore->lda;
kusano 7d535a
    ldx    = Xstore->lda;
kusano 7d535a
    nrhs   = B->ncol;
kusano 7d535a
    
kusano 7d535a
    /* Test the input parameters */
kusano 7d535a
    *info = 0;
kusano 7d535a
    notran = (trans == NOTRANS);
kusano 7d535a
    if ( !notran && trans != TRANS && trans != CONJ ) *info = -1;
kusano 7d535a
    else if ( A->nrow != A->ncol || A->nrow < 0 ||
kusano 7d535a
	      A->Stype != SLU_NC || A->Dtype != SLU_S || A->Mtype != SLU_GE )
kusano 7d535a
	*info = -2;
kusano 7d535a
    else if ( L->nrow != L->ncol || L->nrow < 0 ||
kusano 7d535a
 	      L->Stype != SLU_SC || L->Dtype != SLU_S || L->Mtype != SLU_TRLU )
kusano 7d535a
	*info = -3;
kusano 7d535a
    else if ( U->nrow != U->ncol || U->nrow < 0 ||
kusano 7d535a
 	      U->Stype != SLU_NC || U->Dtype != SLU_S || U->Mtype != SLU_TRU )
kusano 7d535a
	*info = -4;
kusano 7d535a
    else if ( ldb < SUPERLU_MAX(0, A->nrow) ||
kusano 7d535a
 	      B->Stype != SLU_DN || B->Dtype != SLU_S || B->Mtype != SLU_GE )
kusano 7d535a
        *info = -10;
kusano 7d535a
    else if ( ldx < SUPERLU_MAX(0, A->nrow) ||
kusano 7d535a
 	      X->Stype != SLU_DN || X->Dtype != SLU_S || X->Mtype != SLU_GE )
kusano 7d535a
	*info = -11;
kusano 7d535a
    if (*info != 0) {
kusano 7d535a
	i = -(*info);
kusano 7d535a
	xerbla_("sgsrfs", &i);
kusano 7d535a
	return;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Quick return if possible */
kusano 7d535a
    if ( A->nrow == 0 || nrhs == 0) {
kusano 7d535a
	for (j = 0; j < nrhs; ++j) {
kusano 7d535a
	    ferr[j] = 0.;
kusano 7d535a
	    berr[j] = 0.;
kusano 7d535a
	}
kusano 7d535a
	return;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    rowequ = lsame_(equed, "R") || lsame_(equed, "B");
kusano 7d535a
    colequ = lsame_(equed, "C") || lsame_(equed, "B");
kusano 7d535a
    
kusano 7d535a
    /* Allocate working space */
kusano 7d535a
    work = floatMalloc(2*A->nrow);
kusano 7d535a
    rwork = (float *) SUPERLU_MALLOC( A->nrow * sizeof(float) );
kusano 7d535a
    iwork = intMalloc(2*A->nrow);
kusano 7d535a
    if ( !work || !rwork || !iwork ) 
kusano 7d535a
        ABORT("Malloc fails for work/rwork/iwork.");
kusano 7d535a
    
kusano 7d535a
    if ( notran ) {
kusano 7d535a
	*(unsigned char *)transc = 'N';
kusano 7d535a
        transt = TRANS;
kusano 7d535a
    } else {
kusano 7d535a
	*(unsigned char *)transc = 'T';
kusano 7d535a
	transt = NOTRANS;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* NZ = maximum number of nonzero elements in each row of A, plus 1 */
kusano 7d535a
    nz     = A->ncol + 1;
kusano 7d535a
    eps    = slamch_("Epsilon");
kusano 7d535a
    safmin = slamch_("Safe minimum");
kusano 7d535a
    /* Set SAFE1 essentially to be the underflow threshold times the
kusano 7d535a
       number of additions in each row. */
kusano 7d535a
    safe1  = nz * safmin;
kusano 7d535a
    safe2  = safe1 / eps;
kusano 7d535a
kusano 7d535a
    /* Compute the number of nonzeros in each row (or column) of A */
kusano 7d535a
    for (i = 0; i < A->nrow; ++i) iwork[i] = 0;
kusano 7d535a
    if ( notran ) {
kusano 7d535a
	for (k = 0; k < A->ncol; ++k)
kusano 7d535a
	    for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) 
kusano 7d535a
		++iwork[Astore->rowind[i]];
kusano 7d535a
    } else {
kusano 7d535a
	for (k = 0; k < A->ncol; ++k)
kusano 7d535a
	    iwork[k] = Astore->colptr[k+1] - Astore->colptr[k];
kusano 7d535a
    }	
kusano 7d535a
kusano 7d535a
    /* Copy one column of RHS B into Bjcol. */
kusano 7d535a
    Bjcol.Stype = B->Stype;
kusano 7d535a
    Bjcol.Dtype = B->Dtype;
kusano 7d535a
    Bjcol.Mtype = B->Mtype;
kusano 7d535a
    Bjcol.nrow  = B->nrow;
kusano 7d535a
    Bjcol.ncol  = 1;
kusano 7d535a
    Bjcol.Store = (void *) SUPERLU_MALLOC( sizeof(DNformat) );
kusano 7d535a
    if ( !Bjcol.Store ) ABORT("SUPERLU_MALLOC fails for Bjcol.Store");
kusano 7d535a
    Bjcol_store = Bjcol.Store;
kusano 7d535a
    Bjcol_store->lda = ldb;
kusano 7d535a
    Bjcol_store->nzval = work; /* address aliasing */
kusano 7d535a
	
kusano 7d535a
    /* Do for each right hand side ... */
kusano 7d535a
    for (j = 0; j < nrhs; ++j) {
kusano 7d535a
	count = 0;
kusano 7d535a
	lstres = 3.;
kusano 7d535a
	Bptr = &Bmat[j*ldb];
kusano 7d535a
	Xptr = &Xmat[j*ldx];
kusano 7d535a
kusano 7d535a
	while (1) { /* Loop until stopping criterion is satisfied. */
kusano 7d535a
kusano 7d535a
	    /* Compute residual R = B - op(A) * X,   
kusano 7d535a
	       where op(A) = A, A**T, or A**H, depending on TRANS. */
kusano 7d535a
	    
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
	    SCOPY(&A->nrow, Bptr, &ione, work, &ione);
kusano 7d535a
#else
kusano 7d535a
	    scopy_(&A->nrow, Bptr, &ione, work, &ione);
kusano 7d535a
#endif
kusano 7d535a
	    sp_sgemv(transc, ndone, A, Xptr, ione, done, work, ione);
kusano 7d535a
kusano 7d535a
	    /* Compute componentwise relative backward error from formula 
kusano 7d535a
	       max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )   
kusano 7d535a
	       where abs(Z) is the componentwise absolute value of the matrix
kusano 7d535a
	       or vector Z.  If the i-th component of the denominator is less
kusano 7d535a
	       than SAFE2, then SAFE1 is added to the i-th component of the   
kusano 7d535a
	       numerator before dividing. */
kusano 7d535a
kusano 7d535a
	    for (i = 0; i < A->nrow; ++i) rwork[i] = fabs( Bptr[i] );
kusano 7d535a
	    
kusano 7d535a
	    /* Compute abs(op(A))*abs(X) + abs(B). */
kusano 7d535a
	    if (notran) {
kusano 7d535a
		for (k = 0; k < A->ncol; ++k) {
kusano 7d535a
		    xk = fabs( Xptr[k] );
kusano 7d535a
		    for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
kusano 7d535a
			rwork[Astore->rowind[i]] += fabs(Aval[i]) * xk;
kusano 7d535a
		}
kusano 7d535a
	    } else {
kusano 7d535a
		for (k = 0; k < A->ncol; ++k) {
kusano 7d535a
		    s = 0.;
kusano 7d535a
		    for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) {
kusano 7d535a
			irow = Astore->rowind[i];
kusano 7d535a
			s += fabs(Aval[i]) * fabs(Xptr[irow]);
kusano 7d535a
		    }
kusano 7d535a
		    rwork[k] += s;
kusano 7d535a
		}
kusano 7d535a
	    }
kusano 7d535a
	    s = 0.;
kusano 7d535a
	    for (i = 0; i < A->nrow; ++i) {
kusano 7d535a
		if (rwork[i] > safe2) {
kusano 7d535a
		    s = SUPERLU_MAX( s, fabs(work[i]) / rwork[i] );
kusano 7d535a
		} else if ( rwork[i] != 0.0 ) {
kusano 7d535a
                    /* Adding SAFE1 to the numerator guards against
kusano 7d535a
                       spuriously zero residuals (underflow). */
kusano 7d535a
		    s = SUPERLU_MAX( s, (safe1 + fabs(work[i])) / rwork[i] );
kusano 7d535a
                }
kusano 7d535a
                /* If rwork[i] is exactly 0.0, then we know the true 
kusano 7d535a
                   residual also must be exactly 0.0. */
kusano 7d535a
	    }
kusano 7d535a
	    berr[j] = s;
kusano 7d535a
kusano 7d535a
	    /* Test stopping criterion. Continue iterating if   
kusano 7d535a
	       1) The residual BERR(J) is larger than machine epsilon, and   
kusano 7d535a
	       2) BERR(J) decreased by at least a factor of 2 during the   
kusano 7d535a
	          last iteration, and   
kusano 7d535a
	       3) At most ITMAX iterations tried. */
kusano 7d535a
kusano 7d535a
	    if (berr[j] > eps && berr[j] * 2. <= lstres && count < ITMAX) {
kusano 7d535a
		/* Update solution and try again. */
kusano 7d535a
		sgstrs (trans, L, U, perm_c, perm_r, &Bjcol, stat, info);
kusano 7d535a
		
kusano 7d535a
#ifdef _CRAY
kusano 7d535a
		SAXPY(&A->nrow, &done, work, &ione,
kusano 7d535a
		       &Xmat[j*ldx], &ione);
kusano 7d535a
#else
kusano 7d535a
		saxpy_(&A->nrow, &done, work, &ione,
kusano 7d535a
		       &Xmat[j*ldx], &ione);
kusano 7d535a
#endif
kusano 7d535a
		lstres = berr[j];
kusano 7d535a
		++count;
kusano 7d535a
	    } else {
kusano 7d535a
		break;
kusano 7d535a
	    }
kusano 7d535a
        
kusano 7d535a
	} /* end while */
kusano 7d535a
kusano 7d535a
	stat->RefineSteps = count;
kusano 7d535a
kusano 7d535a
	/* Bound error from formula:
kusano 7d535a
	   norm(X - XTRUE) / norm(X) .le. FERR = norm( abs(inv(op(A)))*   
kusano 7d535a
	   ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)   
kusano 7d535a
          where   
kusano 7d535a
            norm(Z) is the magnitude of the largest component of Z   
kusano 7d535a
            inv(op(A)) is the inverse of op(A)   
kusano 7d535a
            abs(Z) is the componentwise absolute value of the matrix or
kusano 7d535a
	       vector Z   
kusano 7d535a
            NZ is the maximum number of nonzeros in any row of A, plus 1   
kusano 7d535a
            EPS is machine epsilon   
kusano 7d535a
kusano 7d535a
          The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))   
kusano 7d535a
          is incremented by SAFE1 if the i-th component of   
kusano 7d535a
          abs(op(A))*abs(X) + abs(B) is less than SAFE2.   
kusano 7d535a
kusano 7d535a
          Use SLACON to estimate the infinity-norm of the matrix   
kusano 7d535a
             inv(op(A)) * diag(W),   
kusano 7d535a
          where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) */
kusano 7d535a
	
kusano 7d535a
	for (i = 0; i < A->nrow; ++i) rwork[i] = fabs( Bptr[i] );
kusano 7d535a
	
kusano 7d535a
	/* Compute abs(op(A))*abs(X) + abs(B). */
kusano 7d535a
	if ( notran ) {
kusano 7d535a
	    for (k = 0; k < A->ncol; ++k) {
kusano 7d535a
		xk = fabs( Xptr[k] );
kusano 7d535a
		for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i)
kusano 7d535a
		    rwork[Astore->rowind[i]] += fabs(Aval[i]) * xk;
kusano 7d535a
	    }
kusano 7d535a
	} else {
kusano 7d535a
	    for (k = 0; k < A->ncol; ++k) {
kusano 7d535a
		s = 0.;
kusano 7d535a
		for (i = Astore->colptr[k]; i < Astore->colptr[k+1]; ++i) {
kusano 7d535a
		    irow = Astore->rowind[i];
kusano 7d535a
		    xk = fabs( Xptr[irow] );
kusano 7d535a
		    s += fabs(Aval[i]) * xk;
kusano 7d535a
		}
kusano 7d535a
		rwork[k] += s;
kusano 7d535a
	    }
kusano 7d535a
	}
kusano 7d535a
	
kusano 7d535a
	for (i = 0; i < A->nrow; ++i)
kusano 7d535a
	    if (rwork[i] > safe2)
kusano 7d535a
		rwork[i] = fabs(work[i]) + (iwork[i]+1)*eps*rwork[i];
kusano 7d535a
	    else
kusano 7d535a
		rwork[i] = fabs(work[i])+(iwork[i]+1)*eps*rwork[i]+safe1;
kusano 7d535a
kusano 7d535a
	kase = 0;
kusano 7d535a
kusano 7d535a
	do {
kusano 7d535a
	    slacon_(&A->nrow, &work[A->nrow], work,
kusano 7d535a
		    &iwork[A->nrow], &ferr[j], &kase);
kusano 7d535a
	    if (kase == 0) break;
kusano 7d535a
kusano 7d535a
	    if (kase == 1) {
kusano 7d535a
		/* Multiply by diag(W)*inv(op(A)**T)*(diag(C) or diag(R)). */
kusano 7d535a
		if ( notran && colequ )
kusano 7d535a
		    for (i = 0; i < A->ncol; ++i) work[i] *= C[i];
kusano 7d535a
		else if ( !notran && rowequ )
kusano 7d535a
		    for (i = 0; i < A->nrow; ++i) work[i] *= R[i];
kusano 7d535a
		
kusano 7d535a
		sgstrs (transt, L, U, perm_c, perm_r, &Bjcol, stat, info);
kusano 7d535a
		
kusano 7d535a
		for (i = 0; i < A->nrow; ++i) work[i] *= rwork[i];
kusano 7d535a
	    } else {
kusano 7d535a
		/* Multiply by (diag(C) or diag(R))*inv(op(A))*diag(W). */
kusano 7d535a
		for (i = 0; i < A->nrow; ++i) work[i] *= rwork[i];
kusano 7d535a
		
kusano 7d535a
		sgstrs (trans, L, U, perm_c, perm_r, &Bjcol, stat, info);
kusano 7d535a
		
kusano 7d535a
		if ( notran && colequ )
kusano 7d535a
		    for (i = 0; i < A->ncol; ++i) work[i] *= C[i];
kusano 7d535a
		else if ( !notran && rowequ )
kusano 7d535a
		    for (i = 0; i < A->ncol; ++i) work[i] *= R[i];  
kusano 7d535a
	    }
kusano 7d535a
	    
kusano 7d535a
	} while ( kase != 0 );
kusano 7d535a
kusano 7d535a
kusano 7d535a
	/* Normalize error. */
kusano 7d535a
	lstres = 0.;
kusano 7d535a
 	if ( notran && colequ ) {
kusano 7d535a
	    for (i = 0; i < A->nrow; ++i)
kusano 7d535a
	    	lstres = SUPERLU_MAX( lstres, C[i] * fabs( Xptr[i]) );
kusano 7d535a
  	} else if ( !notran && rowequ ) {
kusano 7d535a
	    for (i = 0; i < A->nrow; ++i)
kusano 7d535a
	    	lstres = SUPERLU_MAX( lstres, R[i] * fabs( Xptr[i]) );
kusano 7d535a
	} else {
kusano 7d535a
	    for (i = 0; i < A->nrow; ++i)
kusano 7d535a
	    	lstres = SUPERLU_MAX( lstres, fabs( Xptr[i]) );
kusano 7d535a
	}
kusano 7d535a
	if ( lstres != 0. )
kusano 7d535a
	    ferr[j] /= lstres;
kusano 7d535a
kusano 7d535a
    } /* for each RHS j ... */
kusano 7d535a
    
kusano 7d535a
    SUPERLU_FREE(work);
kusano 7d535a
    SUPERLU_FREE(rwork);
kusano 7d535a
    SUPERLU_FREE(iwork);
kusano 7d535a
    SUPERLU_FREE(Bjcol.Store);
kusano 7d535a
kusano 7d535a
    return;
kusano 7d535a
kusano 7d535a
} /* sgsrfs */