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_cdefs.h"
kusano 7d535a
kusano 7d535a
int cgst04(int n, int nrhs, complex *x, int ldx, complex *xact,
kusano 7d535a
	      int ldxact, float rcond, float *resid)
kusano 7d535a
{
kusano 7d535a
/*
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    CGST04 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) 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) 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) 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) FLOAT 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
    float 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
    float xnorm;
kusano 7d535a
    float eps;
kusano 7d535a
    float diffnm;
kusano 7d535a
kusano 7d535a
    /* Function prototypes */
kusano 7d535a
    extern int icamax_(int *, complex *, 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 = slamch_("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 = icamax_(&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
} /* cgst04_ */