void tlin::factorize(SuperMatrix *A, SuperFactors *&F, superlu_options_t *opt) {
  assert(A->nrow == A->ncol);
  int n = A->nrow;

  if (!F) F = (SuperFactors *)SUPERLU_MALLOC(sizeof(SuperFactors));

  if (!opt) opt = &defaultOpt;

  F->perm_c = intMalloc(n);

  get_perm_c(3, A, F->perm_c);

  SuperMatrix AC;
  int *etree = intMalloc(n);

  sp_preorder(opt, A, F->perm_c, etree, &AC);

  F->L      = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));
  F->U      = (SuperMatrix *)SUPERLU_MALLOC(sizeof(SuperMatrix));
  F->perm_r = intMalloc(n);

  SuperLUStat_t stat;
  StatInit(&stat);

  int result;
  dgstrf(opt, &AC, sp_ienv(1), sp_ienv(2), etree, NULL, 0, F->perm_c, F->perm_r,
         F->L, F->U, &stat, &result);

  StatFree(&stat);

  Destroy_CompCol_Permuted(&AC);
  SUPERLU_FREE(etree);

  if (result != 0) freeF(F), F = 0;
}

/*! \brief Set up pointers for real working arrays.
 */
void
sSetRWork(int m, int panel_size, float *dworkptr,
	 float **dense, float **tempv)
{
    float zero = 0.0;

    int maxsuper = SUPERLU_MAX( sp_ienv(3), sp_ienv(7) ),
        rowblk   = sp_ienv(4);
    *dense = dworkptr;
    *tempv = *dense + panel_size*m;
    sfill (*dense, m * panel_size, zero);
    sfill (*tempv, NUM_TEMPV(m,panel_size,maxsuper,rowblk), zero);     
}

/*! \brief Set up pointers for real working arrays.
 */
void
zSetRWork(int m, int panel_size, doublecomplex *dworkptr,
	 doublecomplex **dense, doublecomplex **tempv)
{
    doublecomplex zero = {0.0, 0.0};

    int maxsuper = SUPERLU_MAX( sp_ienv(3), sp_ienv(7) ),
        rowblk   = sp_ienv(4);
    *dense = dworkptr;
    *tempv = *dense + panel_size*m;
    zfill (*dense, m * panel_size, zero);
    zfill (*tempv, NUM_TEMPV(m,panel_size,maxsuper,rowblk), zero);     
}

/*! \brief 
 *
 * <pre>
 * mem_usage consists of the following fields:
 *    - for_lu (float)
 *      The amount of space used in bytes for the L\U data structures.
 *    - total_needed (float)
 *      The amount of space needed in bytes to perform factorization.
 * </pre>
 */
int sQuerySpace(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage)
{
    SCformat *Lstore;
    NCformat *Ustore;
    register int n, iword, dword, panel_size = sp_ienv(1);

    Lstore = L->Store;
    Ustore = U->Store;
    n = L->ncol;
    iword = sizeof(int);
    dword = sizeof(float);

    /* For LU factors */
    mem_usage->for_lu = (float)( (4.0*n + 3.0) * iword +
                                 Lstore->nzval_colptr[n] * dword +
                                 Lstore->rowind_colptr[n] * iword );
    mem_usage->for_lu += (float)( (n + 1.0) * iword +
				 Ustore->colptr[n] * (dword + iword) );

    /* Working storage to support factorization */
    mem_usage->total_needed = mem_usage->for_lu +
	(float)( (2.0 * panel_size + 4.0 + NO_MARKER) * n * iword +
		(panel_size + 1.0) * n * dword );

    return 0;
} /* sQuerySpace */

/*! \brief
 *
 * <pre>
 * mem_usage consists of the following fields:
 *    - for_lu (float)
 *      The amount of space used in bytes for the L\U data structures.
 *    - total_needed (float)
 *      The amount of space needed in bytes to perform factorization.
 * </pre>
 */
int ilu_sQuerySpace(SuperMatrix *L, SuperMatrix *U, mem_usage_t *mem_usage)
{
    SCformat *Lstore;
    NCformat *Ustore;
    register int n, panel_size = sp_ienv(1);
    register float iword, dword;

    Lstore = L->Store;
    Ustore = U->Store;
    n = L->ncol;
    iword = sizeof(int);
    dword = sizeof(double);

    /* For LU factors */
    mem_usage->for_lu = (float)( (4.0f * n + 3.0f) * iword +
				 Lstore->nzval_colptr[n] * dword +
				 Lstore->rowind_colptr[n] * iword );
    mem_usage->for_lu += (float)( (n + 1.0f) * iword +
				 Ustore->colptr[n] * (dword + iword) );

    /* Working storage to support factorization.
       ILU needs 5*n more integers than LU */
    mem_usage->total_needed = mem_usage->for_lu +
	(float)( (2.0f * panel_size + 9.0f + NO_MARKER) * n * iword +
		(panel_size + 1.0f) * n * dword );

    return 0;
} /* ilu_sQuerySpace */

/*! \brief Allocate known working storage. Returns 0 if success, otherwise
   returns the number of bytes allocated so far when failure occurred. */
int
sLUWorkInit(int m, int n, int panel_size, int **iworkptr, 
            float **dworkptr, GlobalLU_t *Glu)
{
    int    isize, dsize, extra;
    float *old_ptr;
    int    maxsuper = SUPERLU_MAX( sp_ienv(3), sp_ienv(7) ),
           rowblk   = sp_ienv(4);

    isize = ( (2 * panel_size + 3 + NO_MARKER ) * m + n ) * sizeof(int);
    dsize = (m * panel_size +
	     NUM_TEMPV(m,panel_size,maxsuper,rowblk)) * sizeof(float);
    
    if ( Glu->MemModel == SYSTEM ) 
	*iworkptr = (int *) intCalloc(isize/sizeof(int));
    else
	*iworkptr = (int *) suser_malloc(isize, TAIL, Glu);
    if ( ! *iworkptr ) {
	fprintf(stderr, "sLUWorkInit: malloc fails for local iworkptr[]\n");
	return (isize + n);
    }

    if ( Glu->MemModel == SYSTEM )
	*dworkptr = (float *) SUPERLU_MALLOC(dsize);
    else {
	*dworkptr = (float *) suser_malloc(dsize, TAIL, Glu);
	if ( NotDoubleAlign(*dworkptr) ) {
	    old_ptr = *dworkptr;
	    *dworkptr = (float*) DoubleAlign(*dworkptr);
	    *dworkptr = (float*) ((double*)*dworkptr - 1);
	    extra = (char*)old_ptr - (char*)*dworkptr;
#ifdef DEBUG	    
	    printf("sLUWorkInit: not aligned, extra %d\n", extra);
#endif	    
	    Glu->stack.top2 -= extra;
	    Glu->stack.used += extra;
	}
    }
    if ( ! *dworkptr ) {
	fprintf(stderr, "malloc fails for local dworkptr[].");
	return (isize + dsize + n);
    }
	
    return 0;
}

void
StatInit(SuperLUStat_t *stat)
{
    register int i, w, panel_size, relax;

    panel_size = sp_ienv(1);
    relax = sp_ienv(2);
    w = SUPERLU_MAX(panel_size, relax);
    stat->panel_histo = intCalloc(w+1);
    stat->utime = (double *) SUPERLU_MALLOC(NPHASES * sizeof(double));
    if (!stat->utime) ABORT("SUPERLU_MALLOC fails for stat->utime");
    stat->ops = (flops_t *) SUPERLU_MALLOC(NPHASES * sizeof(flops_t));
    if (!stat->ops) ABORT("SUPERLU_MALLOC fails for stat->ops");
    for (i = 0; i < NPHASES; ++i) {
        stat->utime[i] = 0.;
        stat->ops[i] = 0.;
    }
}

  //! \brief The constructor initializes the default input options.
  explicit SuperLUdata(int numThreads = 0) :
    A{}, L{}, U{}
  {
    equed[0] = 0;
    R = C = 0;
    perm_r = perm_c = etree = 0;
    rcond = rpg = 0.0;
    if (numThreads > 0)
    {
      opts = new sluop_t;
#ifdef HAS_SUPERLU_MT
      opts->nprocs = numThreads;
      opts->fact = DOFACT;
      opts->trans = NOTRANS;
      opts->refact = NO;
      opts->panel_size = sp_ienv(1);
      opts->relax = sp_ienv(2);
      opts->diag_pivot_thresh = 1.0;
      opts->drop_tol = 0.0;
      opts->ColPerm = MMD_ATA;
      opts->usepr = NO;
      opts->SymmetricMode = NO;
      opts->PrintStat = NO;
      opts->perm_c = 0;
      opts->perm_r = 0;
      opts->work = 0;
      opts->lwork = 0;
      opts->etree = 0;
      opts->colcnt_h = 0;
      opts->part_super_h = 0;
#else
      set_default_options(opts);
      opts->SymmetricMode = YES;
      opts->ColPerm = MMD_AT_PLUS_A;
      opts->DiagPivotThresh = 0.001;
#endif
    }
    else
      opts = 0;
  }

main(int argc, char *argv[])
{
/* 
 * Purpose
 * =======
 *
 * SDRIVE is the main test program for the FLOAT linear 
 * equation driver routines SGSSV and SGSSVX.
 * 
 * The program is invoked by a shell script file -- stest.csh.
 * The output from the tests are written into a file -- stest.out.
 *
 * =====================================================================
 */
    float         *a, *a_save;
    int            *asub, *asub_save;
    int            *xa, *xa_save;
    SuperMatrix  A, B, X, L, U;
    SuperMatrix  ASAV, AC;
    GlobalLU_t   Glu; /* Not needed on return. */
    mem_usage_t    mem_usage;
    int            *perm_r; /* row permutation from partial pivoting */
    int            *perm_c, *pc_save; /* column permutation */
    int            *etree;
    float  zero = 0.0;
    float         *R, *C;
    float         *ferr, *berr;
    float         *rwork;
    float	   *wwork;
    void           *work;
    int            info, lwork, nrhs, panel_size, relax;
    int            m, n, nnz;
    float         *xact;
    float         *rhsb, *solx, *bsav;
    int            ldb, ldx;
    float         rpg, rcond;
    int            i, j, k1;
    float         rowcnd, colcnd, amax;
    int            maxsuper, rowblk, colblk;
    int            prefact, nofact, equil, iequed;
    int            nt, nrun, nfail, nerrs, imat, fimat, nimat;
    int            nfact, ifact, itran;
    int            kl, ku, mode, lda;
    int            zerot, izero, ioff;
    double         u;
    float         anorm, cndnum;
    float         *Afull;
    float         result[NTESTS];
    superlu_options_t options;
    fact_t         fact;
    trans_t        trans;
    SuperLUStat_t  stat;
    static char    matrix_type[8];
    static char    equed[1], path[4], sym[1], dist[1];
    FILE           *fp;

    /* Fixed set of parameters */
    int            iseed[]  = {1988, 1989, 1990, 1991};
    static char    equeds[]  = {'N', 'R', 'C', 'B'};
    static fact_t  facts[] = {FACTORED, DOFACT, SamePattern,
			      SamePattern_SameRowPerm};
    static trans_t transs[]  = {NOTRANS, TRANS, CONJ};

    /* Some function prototypes */ 
    extern int sgst01(int, int, SuperMatrix *, SuperMatrix *, 
		      SuperMatrix *, int *, int *, float *);
    extern int sgst02(trans_t, int, int, int, SuperMatrix *, float *,
                      int, float *, int, float *resid);
    extern int sgst04(int, int, float *, int, 
                      float *, int, float rcond, float *resid);
    extern int sgst07(trans_t, int, int, SuperMatrix *, float *, int,
                         float *, int, float *, int, 
                         float *, float *, float *);
    extern int slatb4_slu(char *, int *, int *, int *, char *, int *, int *, 
	               float *, int *, float *, char *);
    extern int slatms_slu(int *, int *, char *, int *, char *, float *d,
                       int *, float *, float *, int *, int *,
                       char *, float *, int *, float *, int *);
    extern int sp_sconvert(int, int, float *, int, int, int,
	                   float *a, int *, int *, int *);


    /* Executable statements */

    strcpy(path, "SGE");
    nrun  = 0;
    nfail = 0;
    nerrs = 0;

    /* Defaults */
    lwork      = 0;
    n          = 1;
    nrhs       = 1;
    panel_size = sp_ienv(1);
    relax      = sp_ienv(2);
    u          = 1.0;
    strcpy(matrix_type, "LA");
    parse_command_line(argc, argv, matrix_type, &n,
		       &panel_size, &relax, &nrhs, &maxsuper,
		       &rowblk, &colblk, &lwork, &u, &fp);
//.........这里部分代码省略.........

void
c_fortran_zgssv_(int *iopt, int *n, int *nnz, int *nrhs, 
                 doublecomplex *values, int *rowind, int *colptr,
                 doublecomplex *b, int *ldb,
		 fptr *f_factors, /* a handle containing the address
				     pointing to the factored matrices */
		 int *info)

{
/* 
 * This routine can be called from Fortran.
 *
 * iopt (input) int
 *      Specifies the operation:
 *      = 1, performs LU decomposition for the first time
 *      = 2, performs triangular solve
 *      = 3, free all the storage in the end
 *
 * f_factors (input/output) fptr* 
 *      If iopt == 1, it is an output and contains the pointer pointing to
 *                    the structure of the factored matrices.
 *      Otherwise, it it an input.
 *
 */
 
    SuperMatrix A, AC, B;
    SuperMatrix *L, *U;
    int *perm_r; /* row permutations from partial pivoting */
    int *perm_c; /* column permutation vector */
    int *etree;  /* column elimination tree */
    SCformat *Lstore;
    NCformat *Ustore;
    int      i, panel_size, permc_spec, relax;
    trans_t  trans;
    mem_usage_t   mem_usage;
    superlu_options_t options;
    SuperLUStat_t stat;
    factors_t *LUfactors;

    trans = TRANS;
    
    if ( *iopt == 1 ) { /* LU decomposition */

        /* Set the default input options. */
        set_default_options(&options);

	/* Initialize the statistics variables. */
	StatInit(&stat);

	/* Adjust to 0-based indexing */
	for (i = 0; i < *nnz; ++i) --rowind[i];
	for (i = 0; i <= *n; ++i) --colptr[i];

	zCreate_CompCol_Matrix(&A, *n, *n, *nnz, values, rowind, colptr,
			       SLU_NC, SLU_Z, SLU_GE);
	L = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
	U = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
	if ( !(perm_r = intMalloc(*n)) ) ABORT("Malloc fails for perm_r[].");
	if ( !(perm_c = intMalloc(*n)) ) ABORT("Malloc fails for perm_c[].");
	if ( !(etree = intMalloc(*n)) ) ABORT("Malloc fails for etree[].");

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

	panel_size = sp_ienv(1);
	relax = sp_ienv(2);

	zgstrf(&options, &AC, relax, panel_size, etree,
                NULL, 0, perm_c, perm_r, L, U, &stat, info);

	if ( *info == 0 ) {
	    Lstore = (SCformat *) L->Store;
	    Ustore = (NCformat *) 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);
	    zQuerySpace(L, U, &mem_usage);
	    printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
		   mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
	} else {
	    printf("zgstrf() error returns INFO= %d\n", *info);
	    if ( *info <= *n ) { /* factorization completes */
		zQuerySpace(L, U, &mem_usage);
		printf("L\\U MB %.3f\ttotal MB needed %.3f\n",
		       mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
	    }
	}
	
	/* Restore to 1-based indexing */
	for (i = 0; i < *nnz; ++i) ++rowind[i];
	for (i = 0; i <= *n; ++i) ++colptr[i];
//.........这里部分代码省略.........

main(int argc, char *argv[])
{
    char           fact[1], equed[1], trans[1], refact[1];
    SuperMatrix  A, L, U;
    SuperMatrix  B, X;
    NCformat       *Astore;
    NCformat       *Ustore;
    SCformat       *Lstore;
    complex         *a;
    int            *asub, *xa;
    int            *perm_r; /* row permutations from partial pivoting */
    int            *perm_c; /* column permutation vector */
    int            *etree;
    void           *work;
    factor_param_t iparam;
    int            info, lwork, nrhs, ldx, panel_size, relax;
    int            m, n, nnz, permc_spec;
    complex         *rhsb, *rhsx, *xact;
    float         *R, *C;
    float         *ferr, *berr;
    float         u, rpg, rcond;
    int            i, firstfact;
    mem_usage_t    mem_usage;
    void    parse_command_line();

    /* Defaults */
    lwork = 0;
    *fact      = 'E';
    *equed     = 'N';
    *trans     = 'N';
    *refact    = 'N';
    nrhs       = 1;
    panel_size = sp_ienv(1);
    relax      = sp_ienv(2);
    u          = 1.0;
    parse_command_line(argc, argv, &lwork, &panel_size, &relax, &u,
		       fact, trans, refact);
    firstfact = lsame_(fact, "F") || lsame_(refact, "Y");

    iparam.panel_size        = panel_size;
    iparam.relax             = relax;
    iparam.diag_pivot_thresh = u;
    iparam.drop_tol          = -1;
    
    if ( lwork > 0 ) {
	work = SUPERLU_MALLOC(lwork);
	if ( !work ) {
	    ABORT("CLINSOLX: cannot allocate work[]");
	}
    }

    
    creadhb(&m, &n, &nnz, &a, &asub, &xa);
    
    cCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_C, SLU_GE);
    Astore = A.Store;
    printf("Dimension %dx%d; # nonzeros %d\n", A.nrow, A.ncol, Astore->nnz);
    
    if ( !(rhsb = complexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsb[].");
    if ( !(rhsx = complexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhsx[].");
    cCreate_Dense_Matrix(&B, m, nrhs, rhsb, m, SLU_DN, SLU_C, SLU_GE);
    cCreate_Dense_Matrix(&X, m, nrhs, rhsx, m, SLU_DN, SLU_C, SLU_GE);
    xact = complexMalloc(n * nrhs);
    ldx = n;
    cGenXtrue(n, nrhs, xact, ldx);
    cFillRHS(trans, nrhs, xact, ldx, &A, &B);
    
    if ( !(etree = intMalloc(n)) ) ABORT("Malloc fails for etree[].");
    if ( !(perm_r = intMalloc(m)) ) ABORT("Malloc fails for perm_r[].");
    if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");

    /*
     * Get column permutation vector perm_c[], according to permc_spec:
     *   permc_spec = 0: natural ordering 
     *   permc_spec = 1: minimum degree on structure of A'*A
     *   permc_spec = 2: minimum degree 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);

    if ( !(R = (float *) SUPERLU_MALLOC(A.nrow * sizeof(float))) ) 
        ABORT("SUPERLU_MALLOC fails for R[].");
    if ( !(C = (float *) SUPERLU_MALLOC(A.ncol * sizeof(float))) )
        ABORT("SUPERLU_MALLOC fails for C[].");
    if ( !(ferr = (float *) SUPERLU_MALLOC(nrhs * sizeof(float))) )
        ABORT("SUPERLU_MALLOC fails for ferr[].");
    if ( !(berr = (float *) SUPERLU_MALLOC(nrhs * sizeof(float))) ) 
        ABORT("SUPERLU_MALLOC fails for berr[].");

    
    /* Solve the system and compute the condition number
       and error bounds using dgssvx.      */
    
    cgssvx(fact, trans, refact, &A, &iparam, perm_c, perm_r, etree,
	   equed, R, C, &L, &U, work, lwork, &B, &X, &rpg, &rcond,
	   ferr, berr, &mem_usage, &info);

    printf("cgssvx(): info %d\n", info);

//.........这里部分代码省略.........

int main ( int argc, char *argv[] )

/******************************************************************************/
/*
  Purpose:

    MAIN is the main program for PSLINSOL.

  Licensing:

    This code is distributed under the GNU LGPL license. 

  Modified:

    10 February 2014

  Author:

    Xiaoye Li
*/
{
  SuperMatrix   A;
  NCformat *Astore;
  float   *a;
  int      *asub, *xa;
  int      *perm_r; /* row permutations from partial pivoting */
  int      *perm_c; /* column permutation vector */
  SuperMatrix   L;       /* factor L */
  SCPformat *Lstore;
  SuperMatrix   U;       /* factor U */
  NCPformat *Ustore;
  SuperMatrix   B;
  int      nrhs, ldx, info, m, n, nnz, b;
  int      nprocs; /* maximum number of processors to use. */
  int      panel_size, relax, maxsup;
  int      permc_spec;
  trans_t  trans;
  float   *xact, *rhs;
  superlu_memusage_t   superlu_memusage;
  void   parse_command_line();

  timestamp ( );
  printf ( "\n" );
  printf ( "PSLINSOL:\n" );
  printf ( "  C/OpenMP version\n" );
  printf ( "  Call the OpenMP version of SuperLU to solve a linear system.\n" );

  nrhs              = 1;
  trans             = NOTRANS;
  nprocs             = 1;
  n                 = 1000;
  b                 = 1;
  panel_size        = sp_ienv(1);
  relax             = sp_ienv(2);
  maxsup            = sp_ienv(3);
/*
  Check for any commandline input.
*/  
  parse_command_line ( argc, argv, &nprocs, &n, &b, &panel_size, 
    &relax, &maxsup );

#if ( PRNTlevel>=1 || DEBUGlevel>=1 )
    cpp_defs();
#endif

#define HB
#if defined( DEN )
    m = n;
    nnz = n * n;
    sband(n, n, nnz, &a, &asub, &xa);
#elif defined( BAND )
    m = n;
    nnz = (2*b+1) * n;
    sband(n, b, nnz, &a, &asub, &xa);
#elif defined( BD )
    nb = 5;
    bs = 200;
    m = n = bs * nb;
    nnz = bs * bs * nb;
    sblockdiag(nb, bs, nnz, &a, &asub, &xa);
#elif defined( HB )
    sreadhb(&m, &n, &nnz, &a, &asub, &xa);
#else    
    sreadmt(&m, &n, &nnz, &a, &asub, &xa);
#endif

    sCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_S, SLU_GE);
    Astore = A.Store;
    printf("Dimension %dx%d; # nonzeros %d\n", A.nrow, A.ncol, Astore->nnz);
    
    if (!(rhs = floatMalloc(m * nrhs))) SUPERLU_ABORT("Malloc fails for rhs[].");
    sCreate_Dense_Matrix(&B, m, nrhs, rhs, m, SLU_DN, SLU_S, SLU_GE);
    xact = floatMalloc(n * nrhs);
    ldx = n;
    sGenXtrue(n, nrhs, xact, ldx);
    sFillRHS(trans, nrhs, xact, ldx, &A, &B);

    if (!(perm_r = intMalloc(m))) SUPERLU_ABORT("Malloc fails for perm_r[].");
    if (!(perm_c = intMalloc(n))) SUPERLU_ABORT("Malloc fails for perm_c[].");

//.........这里部分代码省略.........

/* Here is a driver inspired by A. Sheffer's "cow flattener". */
static NLboolean __nlFactorize_SUPERLU(__NLContext *context, NLint *permutation) {

	/* OpenNL Context */
	__NLSparseMatrix* M = (context->least_squares)? &context->MtM: &context->M;
	NLuint n = context->n;
	NLuint nnz = __nlSparseMatrixNNZ(M); /* number of non-zero coeffs */

	/* Compressed Row Storage matrix representation */
	NLint	*xa		= __NL_NEW_ARRAY(NLint, n+1);
	NLfloat	*rhs	= __NL_NEW_ARRAY(NLfloat, n);
	NLfloat	*a		= __NL_NEW_ARRAY(NLfloat, nnz);
	NLint	*asub	= __NL_NEW_ARRAY(NLint, nnz);
	NLint	*etree	= __NL_NEW_ARRAY(NLint, n);

	/* SuperLU variables */
	SuperMatrix At, AtP;
	NLint info, panel_size, relax;
	superlu_options_t options;

	/* Temporary variables */
	NLuint i, jj, count;
	
	__nl_assert(!(M->storage & __NL_SYMMETRIC));
	__nl_assert(M->storage & __NL_ROWS);
	__nl_assert(M->m == M->n);
	
	/* Convert M to compressed column format */
	for(i=0, count=0; i<n; i++) {
		__NLRowColumn *Ri = M->row + i;
		xa[i] = count;

		for(jj=0; jj<Ri->size; jj++, count++) {
			a[count] = Ri->coeff[jj].value;
			asub[count] = Ri->coeff[jj].index;
		}
	}
	xa[n] = nnz;

	/* Free M, don't need it anymore at this point */
	__nlSparseMatrixClear(M);

	/* Create superlu A matrix transposed */
	sCreate_CompCol_Matrix(
		&At, n, n, nnz, a, asub, xa, 
		SLU_NC,		/* Colum wise, no supernode */
		SLU_S,		/* floats */ 
		SLU_GE		/* general storage */
	);

	/* Set superlu options */
	set_default_options(&options);
	options.ColPerm = MY_PERMC;
	options.Fact = DOFACT;

	StatInit(&(context->slu.stat));

	panel_size = sp_ienv(1); /* sp_ienv give us the defaults */
	relax = sp_ienv(2);

	/* Compute permutation and permuted matrix */
	context->slu.perm_r = __NL_NEW_ARRAY(NLint, n);
	context->slu.perm_c = __NL_NEW_ARRAY(NLint, n);

	if ((permutation == NULL) || (*permutation == -1)) {
		get_perm_c(3, &At, context->slu.perm_c);

		if (permutation)
			memcpy(permutation, context->slu.perm_c, sizeof(NLint)*n);
	}
	else
		memcpy(context->slu.perm_c, permutation, sizeof(NLint)*n);

	sp_preorder(&options, &At, context->slu.perm_c, etree, &AtP);

	/* Decompose into L and U */
	sgstrf(&options, &AtP, relax, panel_size,
		etree, NULL, 0, context->slu.perm_c, context->slu.perm_r,
		&(context->slu.L), &(context->slu.U), &(context->slu.stat), &info);

	/* Cleanup */

	Destroy_SuperMatrix_Store(&At);
	Destroy_CompCol_Permuted(&AtP);

	__NL_DELETE_ARRAY(etree);
	__NL_DELETE_ARRAY(xa);
	__NL_DELETE_ARRAY(rhs);
	__NL_DELETE_ARRAY(a);
	__NL_DELETE_ARRAY(asub);

	context->slu.alloc_slu = NL_TRUE;

	return (info == 0);
}

void
cpanel_bmod (
            const int  m,          /* in - number of rows in the matrix */
            const int  w,          /* in */
            const int  jcol,       /* in */
            const int  nseg,       /* in */
            complex     *dense,     /* out, of size n by w */
            complex     *tempv,     /* working array */
            int        *segrep,    /* in */
            int        *repfnz,    /* in, of size n by w */
            GlobalLU_t *Glu,       /* modified */
            SuperLUStat_t *stat    /* output */
            )
{


#ifdef USE_VENDOR_BLAS
#ifdef _CRAY
    _fcd ftcs1 = _cptofcd("L", strlen("L")),
         ftcs2 = _cptofcd("N", strlen("N")),
         ftcs3 = _cptofcd("U", strlen("U"));
#endif
    int          incx = 1, incy = 1;
    complex       alpha, beta;
#endif

    register int k, ksub;
    int          fsupc, nsupc, nsupr, nrow;
    int          krep, krep_ind;
    complex       ukj, ukj1, ukj2;
    int          luptr, luptr1, luptr2;
    int          segsze;
    int          block_nrow;  /* no of rows in a block row */
    register int lptr;        /* Points to the row subscripts of a supernode */
    int          kfnz, irow, no_zeros;
    register int isub, isub1, i;
    register int jj;          /* Index through each column in the panel */
    int          *xsup, *supno;
    int          *lsub, *xlsub;
    complex       *lusup;
    int          *xlusup;
    int          *repfnz_col; /* repfnz[] for a column in the panel */
    complex       *dense_col;  /* dense[] for a column in the panel */
    complex       *tempv1;             /* Used in 1-D update */
    complex       *TriTmp, *MatvecTmp; /* used in 2-D update */
    complex      zero = {0.0, 0.0};
    complex      one = {1.0, 0.0};
    complex      comp_temp, comp_temp1;
    register int ldaTmp;
    register int r_ind, r_hi;
    static   int first = 1, maxsuper, rowblk, colblk;
    flops_t  *ops = stat->ops;

    xsup    = Glu->xsup;
    supno   = Glu->supno;
    lsub    = Glu->lsub;
    xlsub   = Glu->xlsub;
    lusup   = Glu->lusup;
    xlusup  = Glu->xlusup;

    if ( first ) {
        maxsuper = SUPERLU_MAX( sp_ienv(3), sp_ienv(7) );
        rowblk   = sp_ienv(4);
        colblk   = sp_ienv(5);
        first = 0;
    }
    ldaTmp = maxsuper + rowblk;

    /*
     * For each nonz supernode segment of U[*,j] in topological order
     */
    k = nseg - 1;
    for (ksub = 0; ksub < nseg; ksub++) { /* for each updating supernode */

        /* krep = representative of current k-th supernode
         * fsupc = first supernodal column
         * nsupc = no of columns in a supernode
         * nsupr = no of rows in a supernode
         */
        krep = segrep[k--];
        fsupc = xsup[supno[krep]];
        nsupc = krep - fsupc + 1;
        nsupr = xlsub[fsupc+1] - xlsub[fsupc];
        nrow = nsupr - nsupc;
        lptr = xlsub[fsupc];
        krep_ind = lptr + nsupc - 1;

        repfnz_col = repfnz;
        dense_col = dense;

        if ( nsupc >= colblk && nrow > rowblk ) { /* 2-D block update */

            TriTmp = tempv;

            /* Sequence through each column in panel -- triangular solves */
            for (jj = jcol; jj < jcol + w; jj++,
                 repfnz_col += m, dense_col += m, TriTmp += ldaTmp ) {

                kfnz = repfnz_col[krep];
                if ( kfnz == EMPTY ) continue;  /* Skip any zero segment */
//.........这里部分代码省略.........

void
dgssv(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r,
      SuperMatrix *L, SuperMatrix *U, SuperMatrix *B,
      SuperLUStat_t *stat, int *info )
{

    DNformat *Bstore;
    SuperMatrix *AA;/* A in SLU_NC format used by the factorization routine.*/
    SuperMatrix AC; /* Matrix postmultiplied by Pc */
    int      lwork = 0, *etree, i;
    
    /* Set default values for some parameters */
    double   drop_tol = 0.;
    int      panel_size;     /* panel size */
    int      relax;          /* no of columns in a relaxed snodes */
    int      permc_spec;
    trans_t  trans = NOTRANS;
    double   *utime;
    double   t;	/* Temporary time */

    /* Test the input parameters ... */
    *info = 0;
    Bstore = B->Store;
    if ( options->Fact != DOFACT ) *info = -1;
    else if ( A->nrow != A->ncol || A->nrow < 0 ||
	 (A->Stype != SLU_NC && A->Stype != SLU_NR) ||
	 A->Dtype != SLU_D || A->Mtype != SLU_GE )
	*info = -2;
    else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||
	B->Stype != SLU_DN || B->Dtype != SLU_D || B->Mtype != SLU_GE )
	*info = -7;
    if ( *info != 0 ) {
	i = -(*info);
	xerbla_("dgssv", &i);
	return;
    }

    utime = stat->utime;

    /* Convert A to SLU_NC format when necessary. */
    if ( A->Stype == SLU_NR ) {
	NRformat *Astore = A->Store;
	AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
	dCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 
			       Astore->nzval, Astore->colind, Astore->rowptr,
			       SLU_NC, A->Dtype, A->Mtype);
	trans = TRANS;
    } else {
        if ( A->Stype == SLU_NC ) AA = A;
    }

    t = SuperLU_timer_();
    /*
     * Get column permutation vector perm_c[], according to permc_spec:
     *   permc_spec = NATURAL:  natural ordering 
     *   permc_spec = MMD_AT_PLUS_A: minimum degree on structure of A'+A
     *   permc_spec = MMD_ATA:  minimum degree on structure of A'*A
     *   permc_spec = COLAMD:   approximate minimum degree column ordering
     *   permc_spec = MY_PERMC: the ordering already supplied in perm_c[]
     */
    permc_spec = options->ColPerm;
    if ( permc_spec != MY_PERMC && options->Fact == DOFACT )
      get_perm_c(permc_spec, AA, perm_c);
    utime[COLPERM] = SuperLU_timer_() - t;

    etree = intMalloc(A->ncol);

    t = SuperLU_timer_();
    sp_preorder(options, AA, perm_c, etree, &AC);
    utime[ETREE] = SuperLU_timer_() - t;

    panel_size = sp_ienv(1);
    relax = sp_ienv(2);

    /*printf("Factor PA = LU ... relax %d\tw %d\tmaxsuper %d\trowblk %d\n", 
	  relax, panel_size, sp_ienv(3), sp_ienv(4));*/
    t = SuperLU_timer_(); 
    /* Compute the LU factorization of A. */
    dgstrf(options, &AC, drop_tol, relax, panel_size,
	   etree, NULL, lwork, perm_c, perm_r, L, U, stat, info);
    utime[FACT] = SuperLU_timer_() - t;

    t = SuperLU_timer_();
    if ( *info == 0 ) {
        /* Solve the system A*X=B, overwriting B with X. */
        dgstrs (trans, L, U, perm_c, perm_r, B, stat, info);
    }
    utime[SOLVE] = SuperLU_timer_() - t;

    SUPERLU_FREE (etree);
    Destroy_CompCol_Permuted(&AC);
    if ( A->Stype == SLU_NR ) {
	Destroy_SuperMatrix_Store(AA);
	SUPERLU_FREE(AA);
    }

}

//.........这里部分代码省略.........
    double   *utime;
    flops_t  *ops, flopcnt;

    /* ------------------------------------------------------------
       Test the input parameters.
       ------------------------------------------------------------*/
    Astore = A->Store;
    Bstore = B->Store;
    *info = 0;
    if ( nprocs <= 0 ) *info = -1;
    else if ( A->nrow != A->ncol || A->nrow < 0 || 
	      (A->Stype != SLU_NC && A->Stype != SLU_NR) ||
	      A->Dtype != SLU_Z || A->Mtype != SLU_GE )
	*info = -2;
    else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(1, A->nrow) )*info = -7;
    if ( *info != 0 ) {
        i = -(*info);
	xerbla_("pzgssv", &i);
	return;
    }

#if 0
    /* Use the best sequential code. 
       if this part is commented out, we will use the parallel code 
       run on one processor. */
    if ( nprocs == 1 ) {
        return;
    }
#endif

    fact               = EQUILIBRATE;
    refact             = NO;
    trans              = NOTRANS;
    panel_size         = sp_ienv(1);
    relax              = sp_ienv(2);
    diag_pivot_thresh  = 1.0;
    usepr              = NO;
    drop_tol           = 0.0;
    work               = NULL;
    lwork              = 0;

    /* ------------------------------------------------------------
       Allocate storage and initialize statistics variables. 
       ------------------------------------------------------------*/
    n = A->ncol;
    StatAlloc(n, nprocs, panel_size, relax, &Gstat);
    StatInit(n, nprocs, &Gstat);
    utime = Gstat.utime;
    ops = Gstat.ops;

    /* ------------------------------------------------------------
       Convert A to NC format when necessary.
       ------------------------------------------------------------*/
    if ( A->Stype == SLU_NR ) {
	NRformat *Astore = A->Store;
	AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
	zCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 
			       Astore->nzval, Astore->colind, Astore->rowptr,
			       SLU_NC, A->Dtype, A->Mtype);
	trans = TRANS;
    } else if ( A->Stype == SLU_NC ) AA = A;

    /* ------------------------------------------------------------
       Initialize the option structure superlumt_options using the
       user-input parameters;
       Apply perm_c to the columns of original A to form AC.

void
sgstrf (superlu_options_t *options, SuperMatrix *A,
        int relax, int panel_size, int *etree, void *work, int lwork,
        int *perm_c, int *perm_r, SuperMatrix *L, SuperMatrix *U,
    	GlobalLU_t *Glu, /* persistent to facilitate multiple factorizations */
        SuperLUStat_t *stat, int *info)
{
    /* Local working arrays */
    NCPformat *Astore;
    int       *iperm_r = NULL; /* inverse of perm_r; used when 
                                  options->Fact == SamePattern_SameRowPerm */
    int       *iperm_c; /* inverse of perm_c */
    int       *iwork;
    float    *swork;
    int	      *segrep, *repfnz, *parent, *xplore;
    int	      *panel_lsub; /* dense[]/panel_lsub[] pair forms a w-wide SPA */
    int	      *xprune;
    int	      *marker;
    float    *dense, *tempv;
    int       *relax_end;
    float    *a;
    int       *asub;
    int       *xa_begin, *xa_end;
    int       *xsup, *supno;
    int       *xlsub, *xlusup, *xusub;
    int       nzlumax;
    float fill_ratio = sp_ienv(6);  /* estimated fill ratio */

    /* Local scalars */
    fact_t    fact = options->Fact;
    double    diag_pivot_thresh = options->DiagPivotThresh;
    int       pivrow;   /* pivotal row number in the original matrix A */
    int       nseg1;	/* no of segments in U-column above panel row jcol */
    int       nseg;	/* no of segments in each U-column */
    register int jcol;	
    register int kcol;	/* end column of a relaxed snode */
    register int icol;
    register int i, k, jj, new_next, iinfo;
    int       m, n, min_mn, jsupno, fsupc, nextlu, nextu;
    int       w_def;	/* upper bound on panel width */
    int       usepr, iperm_r_allocated = 0;
    int       nnzL, nnzU;
    int       *panel_histo = stat->panel_histo;
    flops_t   *ops = stat->ops;

    iinfo    = 0;
    m        = A->nrow;
    n        = A->ncol;
    min_mn   = SUPERLU_MIN(m, n);
    Astore   = A->Store;
    a        = Astore->nzval;
    asub     = Astore->rowind;
    xa_begin = Astore->colbeg;
    xa_end   = Astore->colend;

    /* Allocate storage common to the factor routines */
    *info = sLUMemInit(fact, work, lwork, m, n, Astore->nnz,
                       panel_size, fill_ratio, L, U, Glu, &iwork, &swork);
    if ( *info ) return;
    
    xsup    = Glu->xsup;
    supno   = Glu->supno;
    xlsub   = Glu->xlsub;
    xlusup  = Glu->xlusup;
    xusub   = Glu->xusub;
    
    SetIWork(m, n, panel_size, iwork, &segrep, &parent, &xplore,
	     &repfnz, &panel_lsub, &xprune, &marker);
    sSetRWork(m, panel_size, swork, &dense, &tempv);
    
    usepr = (fact == SamePattern_SameRowPerm);
    if ( usepr ) {
	/* Compute the inverse of perm_r */
	iperm_r = (int *) intMalloc(m);
	for (k = 0; k < m; ++k) iperm_r[perm_r[k]] = k;
	iperm_r_allocated = 1;
    }
    iperm_c = (int *) intMalloc(n);
    for (k = 0; k < n; ++k) iperm_c[perm_c[k]] = k;

    /* Identify relaxed snodes */
    relax_end = (int *) intMalloc(n);
    if ( options->SymmetricMode == YES ) {
        heap_relax_snode(n, etree, relax, marker, relax_end); 
    } else {
        relax_snode(n, etree, relax, marker, relax_end); 
    }
    
    ifill (perm_r, m, EMPTY);
    ifill (marker, m * NO_MARKER, EMPTY);
    supno[0] = -1;
    xsup[0]  = xlsub[0] = xusub[0] = xlusup[0] = 0;
    w_def    = panel_size;

    /* 
     * Work on one "panel" at a time. A panel is one of the following: 
     *	   (a) a relaxed supernode at the bottom of the etree, or
     *	   (b) panel_size contiguous columns, defined by the user
     */
    for (jcol = 0; jcol < min_mn; ) {
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		rcmax = SUPERLU_MAX(rcmax, C[j]);
	    }
	    if (rcmin <= 0.) *info = -8;
	    else if (A->nrow > 0)
		colcnd = SUPERLU_MAX(rcmin,smlnum) / SUPERLU_MIN(rcmax,bignum);
	    else colcnd = 1.;
	}
	if (*info == 0) {
	    if ( lwork < -1 ) *info = -12;
	    else if ( B->ncol < 0 ) *info = -13;
	    else if ( B->ncol > 0 ) { /* no checking if B->ncol=0 */
	         if ( Bstore->lda < SUPERLU_MAX(0, A->nrow) ||
		      B->Stype != SLU_DN || B->Dtype != SLU_S || 
		      B->Mtype != SLU_GE )
		*info = -13;
            }
	    if ( X->ncol < 0 ) *info = -14;
            else if ( X->ncol > 0 ) { /* no checking if X->ncol=0 */
                 if ( Xstore->lda < SUPERLU_MAX(0, A->nrow) ||
		      (B->ncol != 0 && B->ncol != X->ncol) ||
                      X->Stype != SLU_DN ||
		      X->Dtype != SLU_S || X->Mtype != SLU_GE )
		*info = -14;
            }
	}
    }
    if (*info != 0) {
	i = -(*info);
	xerbla_("sgssvx", &i);
	return;
    }
    
    /* Initialization for factor parameters */
    panel_size = sp_ienv(1);
    relax      = sp_ienv(2);
    diag_pivot_thresh = options->DiagPivotThresh;

    utime = stat->utime;
    
    /* Convert A to SLU_NC format when necessary. */
    if ( A->Stype == SLU_NR ) {
	NRformat *Astore = A->Store;
	AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );
	sCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 
			       Astore->nzval, Astore->colind, Astore->rowptr,
			       SLU_NC, A->Dtype, A->Mtype);
	if ( notran ) { /* Reverse the transpose argument. */
	    trant = TRANS;
	    notran = 0;
	} else {
	    trant = NOTRANS;
	    notran = 1;
	}
    } else { /* A->Stype == SLU_NC */
	trant = options->Trans;
	AA = A;
    }

    if ( nofact && equil ) {
	t0 = SuperLU_timer_();
	/* Compute row and column scalings to equilibrate the matrix A. */
	sgsequ(AA, R, C, &rowcnd, &colcnd, &amax, &info1);
	
	if ( info1 == 0 ) {
	    /* Equilibrate matrix A. */
	    slaqgs(AA, R, C, rowcnd, colcnd, amax, equed);