/*
 * -- SuperLU MT routine (version 1.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * August 15, 1997
 *
 * Purpose: This program illustrates how to perform the multiple 
 *          factorizations for the matrix with same sparsity pattern yet
 *          different numerical values. Here, computing a fill-reducing
 *          ordering and performing column permutation are done only once.
 *          In addition, the memory for the factors L and U is allocated
 *          once, and re-used in the subsequent factorizations.
 * 
 */
#include "pdsp_defs.h"
#include "util.h"

main(int argc, char *argv[])
{
    SuperMatrix A, AC, L, U, B;
    NCformat    *Astore;
    SCPformat   *Lstore;
    NCPformat   *Ustore;
    pdgstrf_options_t pdgstrf_options;
    pxgstrf_shared_t pxgstrf_shared;
    pdgstrf_threadarg_t *pdgstrf_threadarg;
    int         nprocs;
    fact_t      fact;
    trans_t     trans;
    yes_no_t    refact, usepr;
    double      u, drop_tol;
    double      *a;
    int         *asub, *xa;
    int         *perm_c; /* column permutation vector */
    int         *perm_r; /* row permutations from partial pivoting */
    void        *work;
    int         info, lwork, nrhs, ldx; 
    int         m, n, nnz, permc_spec, panel_size, relax;
    int         i, firstfact;
    double      *rhsb, *xact;
    Gstat_t Gstat;
    flops_t     flopcnt;
    void parse_command_line();

    /* Default parameters to control factorization. */
    nprocs = 1;
    fact  = DOFACT;
    trans = NOTRANS;
    panel_size = sp_ienv(1);
    relax = sp_ienv(2);
    u     = 1.0;
    usepr = NO;
    drop_tol = 0.0;
    work = NULL;
    lwork = 0;
    nrhs  = 1;

    /* Get the number of processes from command line. */
    parse_command_line(argc, argv, &nprocs);

    /* Read the input matrix stored in Harwell-Boeing format. */
    dreadhb(&m, &n, &nnz, &a, &asub, &xa);

    /* Set up the sparse matrix data structure for A. */
    dCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_D, SLU_GE);

    if ( !(rhsb = doubleMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsb[].");
    dCreate_Dense_Matrix(&B, m, nrhs, rhsb, m, SLU_DN, SLU_D, SLU_GE);
    xact = doubleMalloc(n * nrhs);
    ldx = n;
    dGenXtrue(n, nrhs, xact, ldx);
    dFillRHS(trans, nrhs, xact, ldx, &A, &B);
    
    if ( !(perm_r = intMalloc(m)) ) ABORT("Malloc fails for perm_r[].");
    if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");


    /********************************
     * THE FIRST TIME FACTORIZATION *
     ********************************/

    /* ------------------------------------------------------------
       Allocate storage and initialize statistics variables. 
       ------------------------------------------------------------*/
    StatAlloc(n, nprocs, panel_size, relax, &Gstat);
    StatInit(n, nprocs, &Gstat);

    /* ------------------------------------------------------------
       Get column permutation vector perm_c[], according to permc_spec:
       permc_spec = 0: natural ordering 
       permc_spec = 1: minimum degree ordering on structure of A'*A
       permc_spec = 2: minimum degree ordering on structure of A'+A
       permc_spec = 3: approximate minimum degree for unsymmetric matrices
       ------------------------------------------------------------*/ 	
    permc_spec = 1;
    get_perm_c(permc_spec, &A, perm_c);

    /* ------------------------------------------------------------
       Initialize the option structure pdgstrf_options using the
       user-input parameters;
       Apply perm_c to the columns of original A to form AC.
       ------------------------------------------------------------*/
    refact= NO;
    pdgstrf_init(nprocs, refact, panel_size, relax,
		 u, usepr, drop_tol, perm_c, perm_r,
		 work, lwork, &A, &AC, &pdgstrf_options, &Gstat);

    /* ------------------------------------------------------------
       Compute the LU factorization of A.
       The following routine will create nprocs threads.
       ------------------------------------------------------------*/
    pdgstrf(&pdgstrf_options, &AC, perm_r, &L, &U, &Gstat, &info);
    
    flopcnt = 0;
    for (i = 0; i < nprocs; ++i) flopcnt += Gstat.procstat[i].fcops;
    Gstat.ops[FACT] = flopcnt;

    /* ------------------------------------------------------------
       Solve the system A*X=B, overwriting B with X.
       ------------------------------------------------------------*/
    dgstrs(trans, &L, &U, perm_r, perm_c, &B, &Gstat, &info);
    
    printf("\n** Result of sparse LU **\n");
    dinf_norm_error(nrhs, &B, xact); /* Check inf. norm of the error */


    /*********************************
     * THE SUBSEQUENT FACTORIZATIONS *
     *********************************/

    /* ------------------------------------------------------------
       Re-initialize statistics variables and options used by the
       factorization routine pdgstrf().
       ------------------------------------------------------------*/
    StatInit(n, nprocs, &Gstat);
    refact= YES;
    pdgstrf_init(nprocs, refact, panel_size, relax,
		 u, usepr, drop_tol, perm_c, perm_r,
		 work, lwork, &A, &AC, &pdgstrf_options, &Gstat);

    /* ------------------------------------------------------------
       Compute the LU factorization of A.
       The following routine will create nprocs threads.
       ------------------------------------------------------------*/
    pdgstrf(&pdgstrf_options, &AC, perm_r, &L, &U, &Gstat, &info);
    
    flopcnt = 0;
    for (i = 0; i < nprocs; ++i) flopcnt += Gstat.procstat[i].fcops;
    Gstat.ops[FACT] = flopcnt;

    /* ------------------------------------------------------------
       Re-generate right-hand side B, then solve A*X= B.
       ------------------------------------------------------------*/
    dFillRHS(trans, nrhs, xact, ldx, &A, &B);
    dgstrs(trans, &L, &U, perm_r, perm_c, &B, &Gstat, &info);

    
     /* ------------------------------------------------------------
       Deallocate storage after factorization.
       ------------------------------------------------------------*/
    pdgstrf_finalize(&pdgstrf_options, &AC);

    printf("\n** Result of sparse LU **\n");
    dinf_norm_error(nrhs, &B, xact); /* Check inf. norm of the error */

    Lstore = (SCPformat *) L.Store;
    Ustore = (NCPformat *) U.Store;
    printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
    printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
    printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz - n);
    fflush(stdout);

    SUPERLU_FREE (rhsb);
    SUPERLU_FREE (xact);
    SUPERLU_FREE (perm_r);
    SUPERLU_FREE (perm_c);
    Destroy_CompCol_Matrix(&A);
    Destroy_SuperMatrix_Store(&B);
    if ( lwork >= 0 ) {
        Destroy_SuperNode_SCP(&L);
        Destroy_CompCol_NCP(&U);
    }
    StatFree(&Gstat);
}


/*  
 * Parse command line to get nprocs, the number of processes.
 */
void
parse_command_line(int argc, char *argv[], int *nprocs)
{
    register int c;
    extern char *optarg;

    while ( (c = getopt(argc, argv, "hp:")) != EOF ) {
	switch (c) {
	  case 'h':
	    printf("Options: (default values are in parenthesis)\n");
	    printf("\t-p <int> - number of processes     ( %d )\n", *nprocs);
	    exit(1);
	    break;
	  case 'p': *nprocs = atoi(optarg); 
	            break;
  	}
    }
}