kusano 7d535a
kusano 7d535a
/*
kusano 7d535a
 * -- SuperLU routine (version 2.0) --
kusano 7d535a
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
kusano 7d535a
 * and Lawrence Berkeley National Lab.
kusano 7d535a
 * November 15, 1997
kusano 7d535a
 *
kusano 7d535a
 */
kusano 7d535a
#include <math.h></math.h>
kusano 7d535a
#include "slu_zdefs.h"
kusano 7d535a
kusano 7d535a
int zgst04(int n, int nrhs, doublecomplex *x, int ldx, doublecomplex *xact,
kusano 7d535a
	      int ldxact, double rcond, double *resid)
kusano 7d535a
{
kusano 7d535a
/*
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    ZGST04 computes the difference between a computed solution and the   
kusano 7d535a
    true solution to a system of linear equations.   
kusano 7d535a
    RESID =  ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS ),   
kusano 7d535a
    where RCOND is the reciprocal of the condition number and EPS is the 
kusano 7d535a
    machine epsilon.   
kusano 7d535a
kusano 7d535a
    Arguments   
kusano 7d535a
    =========   
kusano 7d535a
kusano 7d535a
    N       (input) INT   
kusano 7d535a
            The number of rows of the matrices X and XACT.  N >= 0.   
kusano 7d535a
kusano 7d535a
    NRHS    (input) INT   
kusano 7d535a
            The number of columns of the matrices X and XACT.  NRHS >= 0. 
kusano 7d535a
kusano 7d535a
    X       (input) DOUBLE COMPLEX PRECISION array, dimension (LDX,NRHS)   
kusano 7d535a
            The computed solution vectors.  Each vector is stored as a   
kusano 7d535a
            column of the matrix X.   
kusano 7d535a
kusano 7d535a
    LDX     (input) INT   
kusano 7d535a
            The leading dimension of the array X.  LDX >= max(1,N).   
kusano 7d535a
kusano 7d535a
    XACT    (input) DOUBLE COMPLEX PRECISION array, dimension( LDX, NRHS )   
kusano 7d535a
            The exact solution vectors.  Each vector is stored as a   
kusano 7d535a
            column of the matrix XACT.   
kusano 7d535a
kusano 7d535a
    LDXACT  (input) INT   
kusano 7d535a
            The leading dimension of the array XACT.  LDXACT >= max(1,N). 
kusano 7d535a
kusano 7d535a
    RCOND   (input) DOUBLE COMPLEX PRECISION   
kusano 7d535a
            The reciprocal of the condition number of the coefficient   
kusano 7d535a
            matrix in the system of equations.   
kusano 7d535a
kusano 7d535a
    RESID   (output) DOUBLE PRECISION   
kusano 7d535a
            The maximum over the NRHS solution vectors of   
kusano 7d535a
            ( norm(X-XACT) * RCOND ) / ( norm(XACT) * EPS )   
kusano 7d535a
kusano 7d535a
    ===================================================================== 
kusano 7d535a
*/
kusano 7d535a
    /* Table of constant values */
kusano 7d535a
    int c__1 = 1;
kusano 7d535a
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    double d__1, d__2, d__3, d__4;
kusano 7d535a
kusano 7d535a
    /* Local variables */
kusano 7d535a
    int    i, j, n__1;
kusano 7d535a
    int    ix;
kusano 7d535a
    double xnorm;
kusano 7d535a
    double eps;
kusano 7d535a
    double diffnm;
kusano 7d535a
kusano 7d535a
    /* Function prototypes */
kusano 7d535a
    extern int izamax_(int *, doublecomplex *, int *);
kusano 7d535a
kusano 7d535a
    /* Quick exit if N = 0 or NRHS = 0. */
kusano 7d535a
   if ( n <= 0 || nrhs <= 0 ) {
kusano 7d535a
	*resid = 0.;
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Exit with RESID = 1/EPS if RCOND is invalid. */
kusano 7d535a
kusano 7d535a
    eps = dlamch_("Epsilon");
kusano 7d535a
    if ( rcond < 0. ) {
kusano 7d535a
	*resid = 1. / eps;
kusano 7d535a
	return 0;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Compute the maximum of norm(X - XACT) / ( norm(XACT) * EPS )   
kusano 7d535a
       over all the vectors X and XACT . */
kusano 7d535a
kusano 7d535a
    *resid = 0.;
kusano 7d535a
    for (j = 0; j < nrhs; ++j) {
kusano 7d535a
	n__1 = n;
kusano 7d535a
	ix = izamax_(&n__1, &xact[j*ldxact], &c__1);
kusano 7d535a
	xnorm = (d__1 = xact[ix-1 + j*ldxact].r, fabs(d__1)) +
kusano 7d535a
		(d__2 = xact[ix-1 + j*ldxact].i, fabs(d__2));
kusano 7d535a
kusano 7d535a
	diffnm = 0.;
kusano 7d535a
	for (i = 0; i < n; ++i) {
kusano 7d535a
	    /* Computing MAX */
kusano 7d535a
	    d__3 = diffnm;
kusano 7d535a
	    d__4 = (d__1 = x[i+j*ldx].r-xact[i+j*ldxact].r, fabs(d__1)) +
kusano 7d535a
                   (d__2 = x[i+j*ldx].i-xact[i+j*ldxact].i, fabs(d__2));
kusano 7d535a
	    diffnm = SUPERLU_MAX(d__3,d__4);
kusano 7d535a
	}
kusano 7d535a
	if (xnorm <= 0.) {
kusano 7d535a
	    if (diffnm > 0.) {
kusano 7d535a
		*resid = 1. / eps;
kusano 7d535a
	    }
kusano 7d535a
	} else {
kusano 7d535a
	    /* Computing MAX */
kusano 7d535a
	    d__1 = *resid, d__2 = diffnm / xnorm * rcond;
kusano 7d535a
	    *resid = SUPERLU_MAX(d__1,d__2);
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
    if (*resid * eps < 1.) {
kusano 7d535a
	*resid /= eps;
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
} /* zgst04_ */