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_zdefs.h"
kusano 7d535a
kusano 7d535a
main(int argc, char *argv[])
kusano 7d535a
{
kusano 7d535a
    SuperMatrix A;
kusano 7d535a
    NCformat *Astore;
kusano 7d535a
    doublecomplex   *a;
kusano 7d535a
    int      *asub, *xa;
kusano 7d535a
    int      *perm_c; /* column permutation vector */
kusano 7d535a
    int      *perm_r; /* row permutations from partial pivoting */
kusano 7d535a
    SuperMatrix L;      /* factor L */
kusano 7d535a
    SCformat *Lstore;
kusano 7d535a
    SuperMatrix U;      /* factor U */
kusano 7d535a
    NCformat *Ustore;
kusano 7d535a
    SuperMatrix B;
kusano 7d535a
    int      nrhs, ldx, info, m, n, nnz;
kusano 7d535a
    doublecomplex   *xact, *rhs;
kusano 7d535a
    mem_usage_t   mem_usage;
kusano 7d535a
    superlu_options_t options;
kusano 7d535a
    SuperLUStat_t stat;
kusano 7d535a
    
kusano 7d535a
#if ( DEBUGlevel>=1 )
kusano 7d535a
    CHECK_MALLOC("Enter main()");
kusano 7d535a
#endif
kusano 7d535a
kusano 7d535a
    /* Set the default input options:
kusano 7d535a
	options.Fact = DOFACT;
kusano 7d535a
        options.Equil = YES;
kusano 7d535a
    	options.ColPerm = COLAMD;
kusano 7d535a
	options.DiagPivotThresh = 1.0;
kusano 7d535a
    	options.Trans = NOTRANS;
kusano 7d535a
    	options.IterRefine = NOREFINE;
kusano 7d535a
    	options.SymmetricMode = NO;
kusano 7d535a
    	options.PivotGrowth = NO;
kusano 7d535a
    	options.ConditionNumber = NO;
kusano 7d535a
    	options.PrintStat = YES;
kusano 7d535a
     */
kusano 7d535a
    set_default_options(&options);
kusano 7d535a
kusano 7d535a
    /* Now we modify the default options to use the symmetric mode. */
kusano 7d535a
    options.SymmetricMode = YES;
kusano 7d535a
    options.ColPerm = MMD_AT_PLUS_A;
kusano 7d535a
    options.DiagPivotThresh = 0.001;
kusano 7d535a
kusano 7d535a
    /* Read the matrix in Harwell-Boeing format. */
kusano 7d535a
    zreadhb(&m, &n, &nnz, &a, &asub, &xa);
kusano 7d535a
kusano 7d535a
    zCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_Z, SLU_GE);
kusano 7d535a
    Astore = A.Store;
kusano 7d535a
    printf("Dimension %dx%d; # nonzeros %d\n", A.nrow, A.ncol, Astore->nnz);
kusano 7d535a
    
kusano 7d535a
    nrhs   = 1;
kusano 7d535a
    if ( !(rhs = doublecomplexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhs[].");
kusano 7d535a
    zCreate_Dense_Matrix(&B, m, nrhs, rhs, m, SLU_DN, SLU_Z, SLU_GE);
kusano 7d535a
    xact = doublecomplexMalloc(n * nrhs);
kusano 7d535a
    ldx = n;
kusano 7d535a
    zGenXtrue(n, nrhs, xact, ldx);
kusano 7d535a
    zFillRHS(options.Trans, nrhs, xact, ldx, &A, &B);
kusano 7d535a
kusano 7d535a
    if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");
kusano 7d535a
    if ( !(perm_r = intMalloc(m)) ) ABORT("Malloc fails for perm_r[].");
kusano 7d535a
kusano 7d535a
    /* Initialize the statistics variables. */
kusano 7d535a
    StatInit(&stat);
kusano 7d535a
    
kusano 7d535a
    zgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);
kusano 7d535a
    
kusano 7d535a
    if ( info == 0 ) {
kusano 7d535a
kusano 7d535a
	/* This is how you could access the solution matrix. */
kusano 7d535a
        doublecomplex *sol = (doublecomplex*) ((DNformat*) B.Store)->nzval; 
kusano 7d535a
kusano 7d535a
	 /* Compute the infinity norm of the error. */
kusano 7d535a
	zinf_norm_error(nrhs, &B, xact);
kusano 7d535a
kusano 7d535a
	Lstore = (SCformat *) L.Store;
kusano 7d535a
	Ustore = (NCformat *) U.Store;
kusano 7d535a
    	printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
kusano 7d535a
    	printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
kusano 7d535a
    	printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz - n);
kusano 7d535a
    	printf("FILL ratio = %.1f\n", (float)(Lstore->nnz + Ustore->nnz - n)/nnz);
kusano 7d535a
	
kusano 7d535a
	zQuerySpace(&L, &U, &mem_usage);
kusano 7d535a
	printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
kusano 7d535a
	       mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
kusano 7d535a
	
kusano 7d535a
    } else {
kusano 7d535a
	printf("zgssv() error returns INFO= %d\n", info);
kusano 7d535a
	if ( info <= n ) { /* factorization completes */
kusano 7d535a
	    zQuerySpace(&L, &U, &mem_usage);
kusano 7d535a
	    printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
kusano 7d535a
		   mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
kusano 7d535a
	}
kusano 7d535a
    }
kusano 7d535a
kusano 7d535a
    if ( options.PrintStat ) StatPrint(&stat);
kusano 7d535a
    StatFree(&stat);
kusano 7d535a
kusano 7d535a
    SUPERLU_FREE (rhs);
kusano 7d535a
    SUPERLU_FREE (xact);
kusano 7d535a
    SUPERLU_FREE (perm_r);
kusano 7d535a
    SUPERLU_FREE (perm_c);
kusano 7d535a
    Destroy_CompCol_Matrix(&A);
kusano 7d535a
    Destroy_SuperMatrix_Store(&B);
kusano 7d535a
    Destroy_SuperNode_Matrix(&L);
kusano 7d535a
    Destroy_CompCol_Matrix(&U);
kusano 7d535a
kusano 7d535a
#if ( DEBUGlevel>=1 )
kusano 7d535a
    CHECK_MALLOC("Exit main()");
kusano 7d535a
#endif
kusano 7d535a
}
kusano 7d535a