|
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_cdefs.h"
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
main(int argc, char *argv[])
|
|
kusano |
7d535a |
{
|
|
kusano |
7d535a |
SuperMatrix A;
|
|
kusano |
7d535a |
NCformat *Astore;
|
|
kusano |
7d535a |
complex *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 |
complex *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 |
creadhb(&m, &n, &nnz, &a, &asub, &xa);
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
cCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_C, 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 = complexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhs[].");
|
|
kusano |
7d535a |
cCreate_Dense_Matrix(&B, m, nrhs, rhs, m, SLU_DN, SLU_C, SLU_GE);
|
|
kusano |
7d535a |
xact = complexMalloc(n * nrhs);
|
|
kusano |
7d535a |
ldx = n;
|
|
kusano |
7d535a |
cGenXtrue(n, nrhs, xact, ldx);
|
|
kusano |
7d535a |
cFillRHS(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 |
cgssv(&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 |
complex *sol = (complex*) ((DNformat*) B.Store)->nzval;
|
|
kusano |
7d535a |
|
|
kusano |
7d535a |
/* Compute the infinity norm of the error. */
|
|
kusano |
7d535a |
cinf_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 |
cQuerySpace(&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("cgssv() error returns INFO= %d\n", info);
|
|
kusano |
7d535a |
if ( info <= n ) { /* factorization completes */
|
|
kusano |
7d535a |
cQuerySpace(&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 |
|