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
/*
kusano 7d535a
 * File name:		sdrive.c
kusano 7d535a
 * Purpose:             MAIN test program
kusano 7d535a
 */
kusano 7d535a
#include <string.h></string.h>
kusano 7d535a
#include "slu_sdefs.h"
kusano 7d535a
kusano 7d535a
#define NTESTS    5      /* Number of test types */
kusano 7d535a
#define NTYPES    11     /* Number of matrix types */
kusano 7d535a
#define NTRAN     2    
kusano 7d535a
#define THRESH    20.0
kusano 7d535a
#define FMT1      "%10s:n=%d, test(%d)=%12.5g\n"
kusano 7d535a
#define	FMT2      "%10s:fact=%4d, trans=%4d, equed=%c, n=%d, imat=%d, test(%d)=%12.5g\n"
kusano 7d535a
#define FMT3      "%10s:info=%d, izero=%d, n=%d, nrhs=%d, imat=%d, nfail=%d\n"
kusano 7d535a
kusano 7d535a
static void
kusano 7d535a
parse_command_line(int argc, char *argv[], char *matrix_type,
kusano 7d535a
		   int *n, int *w, int *relax, int *nrhs, int *maxsuper,
kusano 7d535a
		   int *rowblk, int *colblk, int *lwork, double *u);
kusano 7d535a
kusano 7d535a
main(int argc, char *argv[])
kusano 7d535a
{
kusano 7d535a
/* 
kusano 7d535a
 * Purpose
kusano 7d535a
 * =======
kusano 7d535a
 *
kusano 7d535a
 * SDRIVE is the main test program for the FLOAT linear 
kusano 7d535a
 * equation driver routines SGSSV and SGSSVX.
kusano 7d535a
 * 
kusano 7d535a
 * The program is invoked by a shell script file -- stest.csh.
kusano 7d535a
 * The output from the tests are written into a file -- stest.out.
kusano 7d535a
 *
kusano 7d535a
 * =====================================================================
kusano 7d535a
 */
kusano 7d535a
    float         *a, *a_save;
kusano 7d535a
    int            *asub, *asub_save;
kusano 7d535a
    int            *xa, *xa_save;
kusano 7d535a
    SuperMatrix  A, B, X, L, U;
kusano 7d535a
    SuperMatrix  ASAV, AC;
kusano 7d535a
    mem_usage_t    mem_usage;
kusano 7d535a
    int            *perm_r; /* row permutation from partial pivoting */
kusano 7d535a
    int            *perm_c, *pc_save; /* column permutation */
kusano 7d535a
    int            *etree;
kusano 7d535a
    float  zero = 0.0;
kusano 7d535a
    float         *R, *C;
kusano 7d535a
    float         *ferr, *berr;
kusano 7d535a
    float         *rwork;
kusano 7d535a
    float	   *wwork;
kusano 7d535a
    void           *work;
kusano 7d535a
    int            info, lwork, nrhs, panel_size, relax;
kusano 7d535a
    int            m, n, nnz;
kusano 7d535a
    float         *xact;
kusano 7d535a
    float         *rhsb, *solx, *bsav;
kusano 7d535a
    int            ldb, ldx;
kusano 7d535a
    float         rpg, rcond;
kusano 7d535a
    int            i, j, k1;
kusano 7d535a
    float         rowcnd, colcnd, amax;
kusano 7d535a
    int            maxsuper, rowblk, colblk;
kusano 7d535a
    int            prefact, nofact, equil, iequed;
kusano 7d535a
    int            nt, nrun, nfail, nerrs, imat, fimat, nimat;
kusano 7d535a
    int            nfact, ifact, itran;
kusano 7d535a
    int            kl, ku, mode, lda;
kusano 7d535a
    int            zerot, izero, ioff;
kusano 7d535a
    double         u;
kusano 7d535a
    float         anorm, cndnum;
kusano 7d535a
    float         *Afull;
kusano 7d535a
    float         result[NTESTS];
kusano 7d535a
    superlu_options_t options;
kusano 7d535a
    fact_t         fact;
kusano 7d535a
    trans_t        trans;
kusano 7d535a
    SuperLUStat_t  stat;
kusano 7d535a
    static char    matrix_type[8];
kusano 7d535a
    static char    equed[1], path[4], sym[1], dist[1];
kusano 7d535a
kusano 7d535a
    /* Fixed set of parameters */
kusano 7d535a
    int            iseed[]  = {1988, 1989, 1990, 1991};
kusano 7d535a
    static char    equeds[]  = {'N', 'R', 'C', 'B'};
kusano 7d535a
    static fact_t  facts[] = {FACTORED, DOFACT, SamePattern,
kusano 7d535a
			      SamePattern_SameRowPerm};
kusano 7d535a
    static trans_t transs[]  = {NOTRANS, TRANS, CONJ};
kusano 7d535a
kusano 7d535a
    /* Some function prototypes */ 
kusano 7d535a
    extern int sgst01(int, int, SuperMatrix *, SuperMatrix *, 
kusano 7d535a
		      SuperMatrix *, int *, int *, float *);
kusano 7d535a
    extern int sgst02(trans_t, int, int, int, SuperMatrix *, float *,
kusano 7d535a
                      int, float *, int, float *resid);
kusano 7d535a
    extern int sgst04(int, int, float *, int, 
kusano 7d535a
                      float *, int, float rcond, float *resid);
kusano 7d535a
    extern int sgst07(trans_t, int, int, SuperMatrix *, float *, int,
kusano 7d535a
                         float *, int, float *, int, 
kusano 7d535a
                         float *, float *, float *);
kusano 7d535a
    extern int slatb4_(char *, int *, int *, int *, char *, int *, int *, 
kusano 7d535a
	               float *, int *, float *, char *);
kusano 7d535a
    extern int slatms_(int *, int *, char *, int *, char *, float *d,
kusano 7d535a
                       int *, float *, float *, int *, int *,
kusano 7d535a
                       char *, float *, int *, float *, int *);
kusano 7d535a
    extern int sp_sconvert(int, int, float *, int, int, int,
kusano 7d535a
	                   float *a, int *, int *, int *);
kusano 7d535a
kusano 7d535a
kusano 7d535a
    /* Executable statements */
kusano 7d535a
kusano 7d535a
    strcpy(path, "SGE");
kusano 7d535a
    nrun  = 0;
kusano 7d535a
    nfail = 0;
kusano 7d535a
    nerrs = 0;
kusano 7d535a
kusano 7d535a
    /* Defaults */
kusano 7d535a
    lwork      = 0;
kusano 7d535a
    n          = 1;
kusano 7d535a
    nrhs       = 1;
kusano 7d535a
    panel_size = sp_ienv(1);
kusano 7d535a
    relax      = sp_ienv(2);
kusano 7d535a
    u          = 1.0;
kusano 7d535a
    strcpy(matrix_type, "LA");
kusano 7d535a
    parse_command_line(argc, argv, matrix_type, &n,
kusano 7d535a
		       &panel_size, &relax, &nrhs, &maxsuper,
kusano 7d535a
		       &rowblk, &colblk, &lwork, &u);
kusano 7d535a
    if ( lwork > 0 ) {
kusano 7d535a
	work = SUPERLU_MALLOC(lwork);
kusano 7d535a
	if ( !work ) {
kusano 7d535a
	    fprintf(stderr, "expert: cannot allocate %d bytes\n", lwork);
kusano 7d535a
	    exit (-1);
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    /* Set the default input options. */
kusano 7d535a
    set_default_options(&options);
kusano 7d535a
    options.DiagPivotThresh = u;
kusano 7d535a
    options.PrintStat = NO;
kusano 7d535a
    options.PivotGrowth = YES;
kusano 7d535a
    options.ConditionNumber = YES;
kusano 7d535a
    options.IterRefine = SINGLE;
kusano 7d535a
    
kusano 7d535a
    if ( strcmp(matrix_type, "LA") == 0 ) {
kusano 7d535a
	/* Test LAPACK matrix suite. */
kusano 7d535a
	m = n;
kusano 7d535a
	lda = SUPERLU_MAX(n, 1);
kusano 7d535a
	nnz = n * n;        /* upper bound */
kusano 7d535a
	fimat = 1;
kusano 7d535a
	nimat = NTYPES;
kusano 7d535a
	Afull = floatCalloc(lda * n);
kusano 7d535a
	sallocateA(n, nnz, &a, &asub, &xa);
kusano 7d535a
    } else {
kusano 7d535a
	/* Read a sparse matrix */
kusano 7d535a
	fimat = nimat = 0;
kusano 7d535a
	sreadhb(&m, &n, &nnz, &a, &asub, &xa);
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    sallocateA(n, nnz, &a_save, &asub_save, &xa_save);
kusano 7d535a
    rhsb = floatMalloc(m * nrhs);
kusano 7d535a
    bsav = floatMalloc(m * nrhs);
kusano 7d535a
    solx = floatMalloc(n * nrhs);
kusano 7d535a
    ldb  = m;
kusano 7d535a
    ldx  = n;
kusano 7d535a
    sCreate_Dense_Matrix(&B, m, nrhs, rhsb, ldb, SLU_DN, SLU_S, SLU_GE);
kusano 7d535a
    sCreate_Dense_Matrix(&X, n, nrhs, solx, ldx, SLU_DN, SLU_S, SLU_GE);
kusano 7d535a
    xact = floatMalloc(n * nrhs);
kusano 7d535a
    etree   = intMalloc(n);
kusano 7d535a
    perm_r  = intMalloc(n);
kusano 7d535a
    perm_c  = intMalloc(n);
kusano 7d535a
    pc_save = intMalloc(n);
kusano 7d535a
    R       = (float *) SUPERLU_MALLOC(m*sizeof(float));
kusano 7d535a
    C       = (float *) SUPERLU_MALLOC(n*sizeof(float));
kusano 7d535a
    ferr    = (float *) SUPERLU_MALLOC(nrhs*sizeof(float));
kusano 7d535a
    berr    = (float *) SUPERLU_MALLOC(nrhs*sizeof(float));
kusano 7d535a
    j = SUPERLU_MAX(m,n) * SUPERLU_MAX(4,nrhs);    
kusano 7d535a
    rwork   = (float *) SUPERLU_MALLOC(j*sizeof(float));
kusano 7d535a
    for (i = 0; i < j; ++i) rwork[i] = 0.;
kusano 7d535a
    if ( !R ) ABORT("SUPERLU_MALLOC fails for R");
kusano 7d535a
    if ( !C ) ABORT("SUPERLU_MALLOC fails for C");
kusano 7d535a
    if ( !ferr ) ABORT("SUPERLU_MALLOC fails for ferr");
kusano 7d535a
    if ( !berr ) ABORT("SUPERLU_MALLOC fails for berr");
kusano 7d535a
    if ( !rwork ) ABORT("SUPERLU_MALLOC fails for rwork");
kusano 7d535a
    wwork   = floatCalloc( SUPERLU_MAX(m,n) * SUPERLU_MAX(4,nrhs) );
kusano 7d535a
kusano 7d535a
    for (i = 0; i < n; ++i) perm_c[i] = pc_save[i] = i;
kusano 7d535a
    options.ColPerm = MY_PERMC;
kusano 7d535a
kusano 7d535a
    for (imat = fimat; imat <= nimat; ++imat) { /* All matrix types */
kusano 7d535a
	
kusano 7d535a
	if ( imat ) {
kusano 7d535a
kusano 7d535a
	    /* Skip types 5, 6, or 7 if the matrix size is too small. */
kusano 7d535a
	    zerot = (imat >= 5 && imat <= 7);
kusano 7d535a
	    if ( zerot && n < imat-4 )
kusano 7d535a
		continue;
kusano 7d535a
	    
kusano 7d535a
	    /* Set up parameters with SLATB4 and generate a test matrix
kusano 7d535a
	       with SLATMS.  */
kusano 7d535a
	    slatb4_(path, &imat, &n, &n, sym, &kl, &ku, &anorm, &mode,
kusano 7d535a
		    &cndnum, dist);
kusano 7d535a
kusano 7d535a
	    slatms_(&n, &n, dist, iseed, sym, &rwork[0], &mode, &cndnum,
kusano 7d535a
		    &anorm, &kl, &ku, "No packing", Afull, &lda,
kusano 7d535a
		    &wwork[0], &info);
kusano 7d535a
kusano 7d535a
	    if ( info ) {
kusano 7d535a
		printf(FMT3, "SLATMS", info, izero, n, nrhs, imat, nfail);
kusano 7d535a
		continue;
kusano 7d535a
	    }
kusano 7d535a
kusano 7d535a
	    /* For types 5-7, zero one or more columns of the matrix
kusano 7d535a
	       to test that INFO is returned correctly.   */
kusano 7d535a
	    if ( zerot ) {
kusano 7d535a
		if ( imat == 5 ) izero = 1;
kusano 7d535a
		else if ( imat == 6 ) izero = n;
kusano 7d535a
		else izero = n / 2 + 1;
kusano 7d535a
		ioff = (izero - 1) * lda;
kusano 7d535a
		if ( imat < 7 ) {
kusano 7d535a
		    for (i = 0; i < n; ++i) Afull[ioff + i] = zero;
kusano 7d535a
		} else {
kusano 7d535a
		    for (j = 0; j < n - izero + 1; ++j)
kusano 7d535a
			for (i = 0; i < n; ++i)
kusano 7d535a
			    Afull[ioff + i + j*lda] = zero;
kusano 7d535a
		}
kusano 7d535a
	    } else {
kusano 7d535a
		izero = 0;
kusano 7d535a
	    }
kusano 7d535a
kusano 7d535a
	    /* Convert to sparse representation. */
kusano 7d535a
	    sp_sconvert(n, n, Afull, lda, kl, ku, a, asub, xa, &nnz);
kusano 7d535a
kusano 7d535a
	} else {
kusano 7d535a
	    izero = 0;
kusano 7d535a
	    zerot = 0;
kusano 7d535a
	}
kusano 7d535a
	
kusano 7d535a
	sCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_S, SLU_GE);
kusano 7d535a
kusano 7d535a
	/* Save a copy of matrix A in ASAV */
kusano 7d535a
	sCreate_CompCol_Matrix(&ASAV, m, n, nnz, a_save, asub_save, xa_save,
kusano 7d535a
			      SLU_NC, SLU_S, SLU_GE);
kusano 7d535a
	sCopy_CompCol_Matrix(&A, &ASAV);
kusano 7d535a
	
kusano 7d535a
	/* Form exact solution. */
kusano 7d535a
	sGenXtrue(n, nrhs, xact, ldx);
kusano 7d535a
	
kusano 7d535a
	StatInit(&stat);
kusano 7d535a
kusano 7d535a
	for (iequed = 0; iequed < 4; ++iequed) {
kusano 7d535a
	    *equed = equeds[iequed];
kusano 7d535a
	    if (iequed == 0) nfact = 4;
kusano 7d535a
	    else nfact = 1; /* Only test factored, pre-equilibrated matrix */
kusano 7d535a
kusano 7d535a
	    for (ifact = 0; ifact < nfact; ++ifact) {
kusano 7d535a
		fact = facts[ifact];
kusano 7d535a
		options.Fact = fact;
kusano 7d535a
kusano 7d535a
		for (equil = 0; equil < 2; ++equil) {
kusano 7d535a
		    options.Equil = equil;
kusano 7d535a
		    prefact   = ( options.Fact == FACTORED ||
kusano 7d535a
				  options.Fact == SamePattern_SameRowPerm );
kusano 7d535a
                                /* Need a first factor */
kusano 7d535a
		    nofact    = (options.Fact != FACTORED);  /* Not factored */
kusano 7d535a
kusano 7d535a
		    /* Restore the matrix A. */
kusano 7d535a
		    sCopy_CompCol_Matrix(&ASAV, &A);
kusano 7d535a
			
kusano 7d535a
		    if ( zerot ) {
kusano 7d535a
                        if ( prefact ) continue;
kusano 7d535a
		    } else if ( options.Fact == FACTORED ) {
kusano 7d535a
                        if ( equil || iequed ) {
kusano 7d535a
			    /* Compute row and column scale factors to
kusano 7d535a
			       equilibrate matrix A.    */
kusano 7d535a
			    sgsequ(&A, R, C, &rowcnd, &colcnd, &amax, &info);
kusano 7d535a
kusano 7d535a
			    /* Force equilibration. */
kusano 7d535a
			    if ( !info && n > 0 ) {
kusano 7d535a
				if ( lsame_(equed, "R") ) {
kusano 7d535a
				    rowcnd = 0.;
kusano 7d535a
				    colcnd = 1.;
kusano 7d535a
				} else if ( lsame_(equed, "C") ) {
kusano 7d535a
				    rowcnd = 1.;
kusano 7d535a
				    colcnd = 0.;
kusano 7d535a
				} else if ( lsame_(equed, "B") ) {
kusano 7d535a
				    rowcnd = 0.;
kusano 7d535a
				    colcnd = 0.;
kusano 7d535a
				}
kusano 7d535a
			    }
kusano 7d535a
			
kusano 7d535a
			    /* Equilibrate the matrix. */
kusano 7d535a
			    slaqgs(&A, R, C, rowcnd, colcnd, amax, equed);
kusano 7d535a
			}
kusano 7d535a
		    }
kusano 7d535a
		    
kusano 7d535a
		    if ( prefact ) { /* Need a factor for the first time */
kusano 7d535a
			
kusano 7d535a
		        /* Save Fact option. */
kusano 7d535a
		        fact = options.Fact;
kusano 7d535a
			options.Fact = DOFACT;
kusano 7d535a
kusano 7d535a
			/* Preorder the matrix, obtain the column etree. */
kusano 7d535a
			sp_preorder(&options, &A, perm_c, etree, &AC);
kusano 7d535a
kusano 7d535a
			/* Factor the matrix AC. */
kusano 7d535a
			sgstrf(&options, &AC, relax, panel_size,
kusano 7d535a
                               etree, work, lwork, perm_c, perm_r, &L, &U,
kusano 7d535a
                               &stat, &info);
kusano 7d535a
kusano 7d535a
			if ( info ) { 
kusano 7d535a
                            printf("** First factor: info %d, equed %c\n",
kusano 7d535a
				   info, *equed);
kusano 7d535a
                            if ( lwork == -1 ) {
kusano 7d535a
                                printf("** Estimated memory: %d bytes\n",
kusano 7d535a
                                        info - n);
kusano 7d535a
                                exit(0);
kusano 7d535a
                            }
kusano 7d535a
                        }
kusano 7d535a
	
kusano 7d535a
                        Destroy_CompCol_Permuted(&AC);
kusano 7d535a
			
kusano 7d535a
		        /* Restore Fact option. */
kusano 7d535a
			options.Fact = fact;
kusano 7d535a
		    } /* if .. first time factor */
kusano 7d535a
		    
kusano 7d535a
		    for (itran = 0; itran < NTRAN; ++itran) {
kusano 7d535a
			trans = transs[itran];
kusano 7d535a
                        options.Trans = trans;
kusano 7d535a
kusano 7d535a
			/* Restore the matrix A. */
kusano 7d535a
			sCopy_CompCol_Matrix(&ASAV, &A);
kusano 7d535a
			
kusano 7d535a
 			/* Set the right hand side. */
kusano 7d535a
			sFillRHS(trans, nrhs, xact, ldx, &A, &B);
kusano 7d535a
			sCopy_Dense_Matrix(m, nrhs, rhsb, ldb, bsav, ldb);
kusano 7d535a
kusano 7d535a
			/*----------------
kusano 7d535a
			 * Test sgssv
kusano 7d535a
			 *----------------*/
kusano 7d535a
			if ( options.Fact == DOFACT && itran == 0) {
kusano 7d535a
                            /* Not yet factored, and untransposed */
kusano 7d535a
	
kusano 7d535a
			    sCopy_Dense_Matrix(m, nrhs, rhsb, ldb, solx, ldx);
kusano 7d535a
			    sgssv(&options, &A, perm_c, perm_r, &L, &U, &X,
kusano 7d535a
                                  &stat, &info);
kusano 7d535a
			    
kusano 7d535a
			    if ( info && info != izero ) {
kusano 7d535a
                                printf(FMT3, "sgssv",
kusano 7d535a
				       info, izero, n, nrhs, imat, nfail);
kusano 7d535a
			    } else {
kusano 7d535a
                                /* Reconstruct matrix from factors and
kusano 7d535a
	                           compute residual. */
kusano 7d535a
                                sgst01(m, n, &A, &L, &U, perm_c, perm_r,
kusano 7d535a
                                         &result[0]);
kusano 7d535a
				nt = 1;
kusano 7d535a
				if ( izero == 0 ) {
kusano 7d535a
				    /* Compute residual of the computed
kusano 7d535a
				       solution. */
kusano 7d535a
				    sCopy_Dense_Matrix(m, nrhs, rhsb, ldb,
kusano 7d535a
						       wwork, ldb);
kusano 7d535a
				    sgst02(trans, m, n, nrhs, &A, solx,
kusano 7d535a
                                              ldx, wwork,ldb, &result[1]);
kusano 7d535a
				    nt = 2;
kusano 7d535a
				}
kusano 7d535a
				
kusano 7d535a
				/* Print information about the tests that
kusano 7d535a
				   did not pass the threshold.      */
kusano 7d535a
				for (i = 0; i < nt; ++i) {
kusano 7d535a
				    if ( result[i] >= THRESH ) {
kusano 7d535a
					printf(FMT1, "sgssv", n, i,
kusano 7d535a
					       result[i]);
kusano 7d535a
					++nfail;
kusano 7d535a
				    }
kusano 7d535a
				}
kusano 7d535a
				nrun += nt;
kusano 7d535a
			    } /* else .. info == 0 */
kusano 7d535a
kusano 7d535a
			    /* Restore perm_c. */
kusano 7d535a
			    for (i = 0; i < n; ++i) perm_c[i] = pc_save[i];
kusano 7d535a
kusano 7d535a
		            if (lwork == 0) {
kusano 7d535a
			        Destroy_SuperNode_Matrix(&L);
kusano 7d535a
			        Destroy_CompCol_Matrix(&U);
kusano 7d535a
			    }
kusano 7d535a
			} /* if .. end of testing sgssv */
kusano 7d535a
    
kusano 7d535a
			/*----------------
kusano 7d535a
			 * Test sgssvx
kusano 7d535a
			 *----------------*/
kusano 7d535a
    
kusano 7d535a
			/* Equilibrate the matrix if fact = FACTORED and
kusano 7d535a
			   equed = 'R', 'C', or 'B'.   */
kusano 7d535a
			if ( options.Fact == FACTORED &&
kusano 7d535a
			     (equil || iequed) && n > 0 ) {
kusano 7d535a
			    slaqgs(&A, R, C, rowcnd, colcnd, amax, equed);
kusano 7d535a
			}
kusano 7d535a
			
kusano 7d535a
			/* Solve the system and compute the condition number
kusano 7d535a
			   and error bounds using sgssvx.      */
kusano 7d535a
			sgssvx(&options, &A, perm_c, perm_r, etree,
kusano 7d535a
                               equed, R, C, &L, &U, work, lwork, &B, &X, &rpg,
kusano 7d535a
                               &rcond, ferr, berr, &mem_usage, &stat, &info);
kusano 7d535a
kusano 7d535a
			if ( info && info != izero ) {
kusano 7d535a
			    printf(FMT3, "sgssvx",
kusano 7d535a
				   info, izero, n, nrhs, imat, nfail);
kusano 7d535a
                            if ( lwork == -1 ) {
kusano 7d535a
                                printf("** Estimated memory: %.0f bytes\n",
kusano 7d535a
                                        mem_usage.total_needed);
kusano 7d535a
                                exit(0);
kusano 7d535a
                            }
kusano 7d535a
			} else {
kusano 7d535a
			    if ( !prefact ) {
kusano 7d535a
			    	/* Reconstruct matrix from factors and
kusano 7d535a
	 			   compute residual. */
kusano 7d535a
                                sgst01(m, n, &A, &L, &U, perm_c, perm_r,
kusano 7d535a
                                         &result[0]);
kusano 7d535a
				k1 = 0;
kusano 7d535a
			    } else {
kusano 7d535a
			   	k1 = 1;
kusano 7d535a
			    }
kusano 7d535a
kusano 7d535a
			    if ( !info ) {
kusano 7d535a
				/* Compute residual of the computed solution.*/
kusano 7d535a
				sCopy_Dense_Matrix(m, nrhs, bsav, ldb,
kusano 7d535a
						  wwork, ldb);
kusano 7d535a
				sgst02(trans, m, n, nrhs, &ASAV, solx, ldx,
kusano 7d535a
					  wwork, ldb, &result[1]);
kusano 7d535a
kusano 7d535a
				/* Check solution from generated exact
kusano 7d535a
				   solution. */
kusano 7d535a
				sgst04(n, nrhs, solx, ldx, xact, ldx, rcond,
kusano 7d535a
					  &result[2]);
kusano 7d535a
kusano 7d535a
				/* Check the error bounds from iterative
kusano 7d535a
				   refinement. */
kusano 7d535a
				sgst07(trans, n, nrhs, &ASAV, bsav, ldb,
kusano 7d535a
					  solx, ldx, xact, ldx, ferr, berr,
kusano 7d535a
					  &result[3]);
kusano 7d535a
kusano 7d535a
				/* Print information about the tests that did
kusano 7d535a
				   not pass the threshold.    */
kusano 7d535a
				for (i = k1; i < NTESTS; ++i) {
kusano 7d535a
				    if ( result[i] >= THRESH ) {
kusano 7d535a
					printf(FMT2, "sgssvx",
kusano 7d535a
					       options.Fact, trans, *equed,
kusano 7d535a
					       n, imat, i, result[i]);
kusano 7d535a
					++nfail;
kusano 7d535a
				    }
kusano 7d535a
				}
kusano 7d535a
				nrun += NTESTS;
kusano 7d535a
			    } /* if .. info == 0 */
kusano 7d535a
			} /* else .. end of testing sgssvx */
kusano 7d535a
kusano 7d535a
		    } /* for itran ... */
kusano 7d535a
kusano 7d535a
		    if ( lwork == 0 ) {
kusano 7d535a
			Destroy_SuperNode_Matrix(&L);
kusano 7d535a
			Destroy_CompCol_Matrix(&U);
kusano 7d535a
		    }
kusano 7d535a
kusano 7d535a
		} /* for equil ... */
kusano 7d535a
	    } /* for ifact ... */
kusano 7d535a
	} /* for iequed ... */
kusano 7d535a
#if 0    
kusano 7d535a
    if ( !info ) {
kusano 7d535a
	PrintPerf(&L, &U, &mem_usage, rpg, rcond, ferr, berr, equed);
kusano 7d535a
    }
kusano 7d535a
#endif    
kusano 7d535a
kusano 7d535a
    } /* for imat ... */
kusano 7d535a
kusano 7d535a
    /* Print a summary of the results. */
kusano 7d535a
    PrintSumm("SGE", nfail, nrun, nerrs);
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE (rhsb);
kusano 7d535a
    SUPERLU_FREE (bsav);
kusano 7d535a
    SUPERLU_FREE (solx);    
kusano 7d535a
    SUPERLU_FREE (xact);
kusano 7d535a
    SUPERLU_FREE (etree);
kusano 7d535a
    SUPERLU_FREE (perm_r);
kusano 7d535a
    SUPERLU_FREE (perm_c);
kusano 7d535a
    SUPERLU_FREE (pc_save);
kusano 7d535a
    SUPERLU_FREE (R);
kusano 7d535a
    SUPERLU_FREE (C);
kusano 7d535a
    SUPERLU_FREE (ferr);
kusano 7d535a
    SUPERLU_FREE (berr);
kusano 7d535a
    SUPERLU_FREE (rwork);
kusano 7d535a
    SUPERLU_FREE (wwork);
kusano 7d535a
    Destroy_SuperMatrix_Store(&B);
kusano 7d535a
    Destroy_SuperMatrix_Store(&X);
kusano 7d535a
    Destroy_CompCol_Matrix(&A);
kusano 7d535a
    Destroy_CompCol_Matrix(&ASAV);
kusano 7d535a
    if ( lwork > 0 ) {
kusano 7d535a
	SUPERLU_FREE (work);
kusano 7d535a
	Destroy_SuperMatrix_Store(&L);
kusano 7d535a
	Destroy_SuperMatrix_Store(&U);
kusano 7d535a
    }
kusano 7d535a
    StatFree(&stat);
kusano 7d535a
kusano 7d535a
    return 0;
kusano 7d535a
}
kusano 7d535a
kusano 7d535a
/*  
kusano 7d535a
 * Parse command line options to get relaxed snode size, panel size, etc.
kusano 7d535a
 */
kusano 7d535a
static void
kusano 7d535a
parse_command_line(int argc, char *argv[], char *matrix_type,
kusano 7d535a
		   int *n, int *w, int *relax, int *nrhs, int *maxsuper,
kusano 7d535a
		   int *rowblk, int *colblk, int *lwork, double *u)
kusano 7d535a
{
kusano 7d535a
    int c;
kusano 7d535a
    extern char *optarg;
kusano 7d535a
kusano 7d535a
    while ( (c = getopt(argc, argv, "ht:n:w:r:s:m:b:c:l:")) != EOF ) {
kusano 7d535a
	switch (c) {
kusano 7d535a
	  case 'h':
kusano 7d535a
	    printf("Options:\n");
kusano 7d535a
	    printf("\t-w <int> - panel size\n");</int>
kusano 7d535a
	    printf("\t-r <int> - granularity of relaxed supernodes\n");</int>
kusano 7d535a
	    exit(1);
kusano 7d535a
	    break;
kusano 7d535a
	  case 't': strcpy(matrix_type, optarg);
kusano 7d535a
	            break;
kusano 7d535a
	  case 'n': *n = atoi(optarg);
kusano 7d535a
	            break;
kusano 7d535a
	  case 'w': *w = atoi(optarg);
kusano 7d535a
	            break;
kusano 7d535a
	  case 'r': *relax = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 's': *nrhs = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 'm': *maxsuper = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 'b': *rowblk = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 'c': *colblk = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 'l': *lwork = atoi(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
	  case 'u': *u = atof(optarg); 
kusano 7d535a
	            break;
kusano 7d535a
  	}
kusano 7d535a
    }
kusano 7d535a
}