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
 */
kusano 7d535a
#include "slu_sdefs.h"
kusano 7d535a
kusano 7d535a
int sgst02(trans_t trans, int m, int n, int nrhs, SuperMatrix *A,
kusano 7d535a
	      float *x, int ldx, float *b, int ldb, float *resid)
kusano 7d535a
{
kusano 7d535a
/*  
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    SGST02 computes the residual for a solution of a system of linear   
kusano 7d535a
    equations  A*x = b  or  A'*x = b:   
kusano 7d535a
       RESID = norm(B - A*X) / ( norm(A) * norm(X) * EPS ),   
kusano 7d535a
    where EPS is the machine epsilon.   
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   
kusano 7d535a
            = TRANS  :  A'*x = b, where A' is the transpose of A   
kusano 7d535a
            = CONJ   :  A'*x = b, where A' is the transpose of A   
kusano 7d535a
kusano 7d535a
    M       (input) INTEGER   
kusano 7d535a
            The number of rows of the matrix A.  M >= 0.   
kusano 7d535a
kusano 7d535a
    N       (input) INTEGER   
kusano 7d535a
            The number of columns of the matrix A.  N >= 0.   
kusano 7d535a
kusano 7d535a
    NRHS    (input) INTEGER   
kusano 7d535a
            The number of columns of B, the matrix of right hand sides.   
kusano 7d535a
            NRHS >= 0.
kusano 7d535a
	    
kusano 7d535a
    A       (input) SuperMatrix*, dimension (LDA,N)   
kusano 7d535a
            The original M x N sparse matrix A.   
kusano 7d535a
kusano 7d535a
    X       (input) FLOAT PRECISION array, dimension (LDX,NRHS)   
kusano 7d535a
            The computed solution vectors for the system of linear   
kusano 7d535a
            equations.   
kusano 7d535a
kusano 7d535a
    LDX     (input) INTEGER   
kusano 7d535a
            The leading dimension of the array X.  If TRANS = NOTRANS,   
kusano 7d535a
            LDX >= max(1,N); if TRANS = TRANS or CONJ, LDX >= max(1,M).   
kusano 7d535a
kusano 7d535a
    B       (input/output) FLOAT PRECISION array, dimension (LDB,NRHS)   
kusano 7d535a
            On entry, the right hand side vectors for the system of   
kusano 7d535a
            linear equations.   
kusano 7d535a
            On exit, B is overwritten with the difference B - A*X.   
kusano 7d535a
kusano 7d535a
    LDB     (input) INTEGER   
kusano 7d535a
            The leading dimension of the array B.  IF TRANS = NOTRANS,
kusano 7d535a
            LDB >= max(1,M); if TRANS = TRANS or CONJ, LDB >= max(1,N).
kusano 7d535a
	    
kusano 7d535a
    RESID   (output) FLOAT PRECISION   
kusano 7d535a
            The maximum over the number of right hand sides of   
kusano 7d535a
            norm(B - A*X) / ( norm(A) * norm(X) * EPS ).   
kusano 7d535a
kusano 7d535a
    =====================================================================
kusano 7d535a
*/
kusano 7d535a
kusano 7d535a
    /* Table of constant values */
kusano 7d535a
    float alpha = -1.;
kusano 7d535a
    float beta  = 1.;
kusano 7d535a
    int    c__1  = 1;
kusano 7d535a
    
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    float d__1, d__2;
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    int j;
kusano 7d535a
    int n1, n2;
kusano 7d535a
    float anorm, bnorm;
kusano 7d535a
    float xnorm;
kusano 7d535a
    float eps;
kusano 7d535a
    char transc[1];
kusano 7d535a
kusano 7d535a
    /* Function prototypes */
kusano 7d535a
    extern int lsame_(char *, char *);
kusano 7d535a
    extern float slangs(char *, SuperMatrix *);
kusano 7d535a
    extern float sasum_(int *, float *, int *);
kusano 7d535a
    
kusano 7d535a
    /* Function Body */
kusano 7d535a
    if ( m <= 0 || n <= 0 || nrhs == 0) {
kusano 7d535a
	*resid = 0.;
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    if ( (trans == TRANS) || (trans == CONJ) ) {
kusano 7d535a
	n1 = n;
kusano 7d535a
	n2 = m;
kusano 7d535a
        *transc = 'T';
kusano 7d535a
    } else {
kusano 7d535a
	n1 = m;
kusano 7d535a
	n2 = n;
kusano 7d535a
	*transc = 'N';
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Exit with RESID = 1/EPS if ANORM = 0. */
kusano 7d535a
kusano 7d535a
    eps = slamch_("Epsilon");
kusano 7d535a
    anorm = slangs("1", A);
kusano 7d535a
    if (anorm <= 0.) {
kusano 7d535a
	*resid = 1. / eps;
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Compute  B - A*X  (or  B - A'*X ) and store in B. */
kusano 7d535a
kusano 7d535a
    sp_sgemm(transc, "N", n1, nrhs, n2, alpha, A, x, ldx, beta, b, ldb);
kusano 7d535a
kusano 7d535a
    /* Compute the maximum over the number of right hand sides of   
kusano 7d535a
       norm(B - A*X) / ( norm(A) * norm(X) * EPS ) . */
kusano 7d535a
kusano 7d535a
    *resid = 0.;
kusano 7d535a
    for (j = 0; j < nrhs; ++j) {
kusano 7d535a
	bnorm = sasum_(&n1, &b[j*ldb], &c__1);
kusano 7d535a
	xnorm = sasum_(&n2, &x[j*ldx], &c__1);
kusano 7d535a
	if (xnorm <= 0.) {
kusano 7d535a
	    *resid = 1. / eps;
kusano 7d535a
	} else {
kusano 7d535a
	    /* Computing MAX */
kusano 7d535a
	    d__1 = *resid, d__2 = bnorm / anorm / xnorm / eps;
kusano 7d535a
	    *resid = SUPERLU_MAX(d__1, d__2);
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
} /* sgst02 */
kusano 7d535a