自学教程：C++ sp_ienv函数代码示例

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. */voidsSetRWork(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. */voidzSetRWork(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. */intsLUWorkInit(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;}

voidStatInit(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);//.........这里部分代码省略.........

voidc_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);}

voidcpanel_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 *///.........这里部分代码省略.........

voiddgssv(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.

voidsgstrf (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);

int HYPRE_ParCSR_SuperLUSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix A_csr,                              HYPRE_ParVector b, HYPRE_ParVector x ){#ifdef HAVE_SUPERLU   int    startRow, endRow, nrows, *partition, *AdiagI, *AdiagJ, nnz;   int    irow, colNum, index, *cscI, *cscJ, jcol, *colLengs;   int    *etree, permcSpec, lwork, panelSize, relax, info;   double *AdiagA, *cscA, diagPivotThresh, dropTol;   char              refact[1];   hypre_CSRMatrix   *Adiag;   HYPRE_SuperLU     *sluPtr;   SuperMatrix       sluAmat, auxAmat;   superlu_options_t slu_options;   SuperLUStat_t     slu_stat;   /* ---------------------------------------------------------------- */   /* get matrix information                                           */   /* ---------------------------------------------------------------- */   sluPtr = (HYPRE_SuperLU *) solver;   assert ( sluPtr != NULL );   HYPRE_ParCSRMatrixGetRowPartitioning( A_csr, &partition );   startRow = partition[0];   endRow   = partition[1] - 1;   nrows    = endRow - startRow + 1;   free( partition );   if ( startRow != 0 )   {      printf("HYPRE_ParCSR_SuperLUSetup ERROR - start row != 0./n");      return -1;   }   /* ---------------------------------------------------------------- */   /* get hypre matrix                                                 */   /* ---------------------------------------------------------------- */   Adiag  = hypre_ParCSRMatrixDiag((hypre_ParCSRMatrix *) A_csr);   AdiagI = hypre_CSRMatrixI(Adiag);   AdiagJ = hypre_CSRMatrixJ(Adiag);   AdiagA = hypre_CSRMatrixData(Adiag);   nnz    = AdiagI[nrows];   /* ---------------------------------------------------------------- */   /* convert the csr matrix into csc matrix                           */   /* ---------------------------------------------------------------- */   colLengs = (int *) malloc(nrows * sizeof(int));   for ( irow = 0; irow < nrows; irow++ ) colLengs[irow] = 0;   for ( irow = 0; irow < nrows; irow++ )      for ( jcol = AdiagI[irow]; jcol < AdiagI[irow+1]; jcol++ )         colLengs[AdiagJ[jcol]]++;   cscJ = (int *)    malloc( (nrows+1) * sizeof(int) );   cscI = (int *)    malloc( nnz * sizeof(int) );   cscA = (double *) malloc( nnz * sizeof(double) );   cscJ[0] = 0;   nnz = 0;   for ( jcol = 1; jcol <= nrows; jcol++ )   {      nnz += colLengs[jcol-1];      cscJ[jcol] = nnz;   }   for ( irow = 0; irow < nrows; irow++ )   {      for ( jcol = AdiagI[irow]; jcol < AdiagI[irow+1]; jcol++ )      {         colNum = AdiagJ[jcol];         index  = cscJ[colNum]++;         cscI[index] = irow;         cscA[index] = AdiagA[jcol];      }   }   cscJ[0] = 0;   nnz = 0;   for ( jcol = 1; jcol <= nrows; jcol++ )   {      nnz += colLengs[jcol-1];      cscJ[jcol] = nnz;   }   free(colLengs);   /* ---------------------------------------------------------------- */   /* create SuperMatrix                                                */   /* ---------------------------------------------------------------- */                                                                                   dCreate_CompCol_Matrix(&sluAmat,nrows,nrows,cscJ[nrows],cscA,cscI,                          cscJ, SLU_NC, SLU_D, SLU_GE);   etree   = (int *) malloc(nrows * sizeof(int));   sluPtr->permC_  = (int *) malloc(nrows * sizeof(int));   sluPtr->permR_  = (int *) malloc(nrows * sizeof(int));   permcSpec = 0;   get_perm_c(permcSpec, &sluAmat, sluPtr->permC_);   slu_options.Fact = DOFACT;   slu_options.SymmetricMode = NO;   sp_preorder(&slu_options, &sluAmat, sluPtr->permC_, etree, &auxAmat);   diagPivotThresh = 1.0;   dropTol = 0.0;   panelSize = sp_ienv(1);   relax = sp_ienv(2);   StatInit(&slu_stat);   lwork = 0;//.........这里部分代码省略.........

//.........这里部分代码省略.........	    else rowcnd = 1.;	}	if (colequ && *info == 0) {	    rcmin = bignum;	    rcmax = 0.;	    for (j = 0; j < A->nrow; ++j) {		rcmin = SUPERLU_MIN(rcmin, C[j]);		rcmax = SUPERLU_MAX(rcmax, C[j]);	    }	    if (rcmin <= 0.) *info = -11;	    else if (A->nrow > 0)		colcnd = SUPERLU_MAX(rcmin,smlnum) / SUPERLU_MIN(rcmax,bignum);	    else colcnd = 1.;	}	if (*info == 0) {	    if ( lwork < -1 ) *info = -15;	    else if ( B->ncol < 0 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||		      B->Stype != DN || B->Dtype != D_ || 		      B->Mtype != GE )		*info = -16;	    else if ( X->ncol < 0 || Xstore->lda < SUPERLU_MAX(0, A->nrow) ||		      B->ncol != X->ncol || X->Stype != DN ||		      X->Dtype != D_ || X->Mtype != GE )		*info = -17;	}    }    if (*info != 0) {	i = -(*info);	xerbla_("dgssvx", &i);	return;    }        /* Default values for factor_params */    panel_size = sp_ienv(1);    relax      = sp_ienv(2);    diag_pivot_thresh = 1.0;    drop_tol   = 0.0;    if ( factor_params != NULL ) {	if ( factor_params->panel_size != -1 )	    panel_size = factor_params->panel_size;	if ( factor_params->relax != -1 ) relax = factor_params->relax;	if ( factor_params->diag_pivot_thresh != -1 )	    diag_pivot_thresh = factor_params->diag_pivot_thresh;	if ( factor_params->drop_tol != -1 )	    drop_tol = factor_params->drop_tol;    }    StatInit(panel_size, relax);    utime = SuperLUStat.utime;        /* Convert A to NC format when necessary. */    if ( A->Stype == 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,			       NC, A->Dtype, A->Mtype);	if ( notran ) { /* Reverse the transpose argument. */	    *trant = 'T';	    notran = 0;	} else {	    *trant = 'N';	    notran = 1;	}    } else { /* A->Stype == NC */	*trant = *trans;

//.........这里部分代码省略.........    }    /* ------------------------------------------------------------       Scale the right hand side.       ------------------------------------------------------------*/    if ( notran ) {        if ( rowequ ) {            for (j = 0; j < nrhs; ++j)                for (i = 0; i < A->nrow; ++i) {                    Bmat[i + j*ldb] *= R[i];                }        }    } else if ( colequ ) {        for (j = 0; j < nrhs; ++j)            for (i = 0; i < A->nrow; ++i) {                Bmat[i + j*ldb] *= C[i];            }    }    /* ------------------------------------------------------------       Perform the LU factorization.       ------------------------------------------------------------*/    if ( dofact || equil ) {        /* Obtain column etree, the column count (colcnt_h) and supernode        partition (part_super_h) for the Householder matrix. */        t0 = SuperLU_timer_();        sp_colorder(AA, perm_c, superlumt_options, &AC);        utime[ETREE] = SuperLU_timer_() - t0;#if ( PRNTlevel >= 2 )        printf("Factor PA = LU ... relax %d/tw %d/tmaxsuper %d/trowblk %d/n",               relax, panel_size, sp_ienv(3), sp_ienv(4));        fflush(stdout);#endif        /* Compute the LU factorization of A*Pc. */        t0 = SuperLU_timer_();        psgstrf(superlumt_options, &AC, perm_r, L, U, &Gstat, info);        utime[FACT] = SuperLU_timer_() - t0;        flopcnt = 0;        for (i = 0; i < nprocs; ++i) flopcnt += Gstat.procstat[i].fcops;        ops[FACT] = flopcnt;        if ( superlumt_options->lwork == -1 ) {            superlu_memusage->total_needed = *info - A->ncol;            return;        }    }    if ( *info > 0 ) {        if ( *info <= A->ncol ) {            /* Compute the reciprocal pivot growth factor of the leading               rank-deficient *info columns of A. */            *recip_pivot_growth = sPivotGrowth(*info, AA, perm_c, L, U);        }    } else {        /* ------------------------------------------------------------           Compute the reciprocal pivot growth factor *recip_pivot_growth.           ------------------------------------------------------------*/        *recip_pivot_growth = sPivotGrowth(A->ncol, AA, perm_c, L, U);        /* ------------------------------------------------------------

//.........这里部分代码省略......... * info    (output) int* *	   = 0: successful exit *         > 0: if info = i, and i is *             <= A->ncol: U(i,i) is exactly zero. The factorization has *                been completed, but the factor U is exactly singular, *                so the solution could not be computed. *             > A->ncol: number of bytes allocated when memory allocation *                failure occurred, plus A->ncol. *    */    DNformat *Bstore;    SuperMatrix *AA = NULL;/* 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 */    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_S || 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_S || B->Mtype != SLU_GE )	*info = -7;    if ( *info != 0 ) {	i = -(*info);	xerbla_("sgssv", &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) );	sCreate_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. */    sgstrf(options, &AC, 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. */        sgstrs (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);    }}

//.........这里部分代码省略.........    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, &marker_relax, &marker);    cSetRWork(m, panel_size, cwork, &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;    swap = (int *)intMalloc(n);    for (k = 0; k < n; k++) swap[k] = iperm_c[k];    iswap = (int *)intMalloc(n);    for (k = 0; k < n; k++) iswap[k] = perm_c[k];    amax = (float *) floatMalloc(panel_size);    if (drop_rule & DROP_SECONDARY)        swork2 = (float *)floatMalloc(n);    else        swork2 = NULL;    nnzAj = 0;    nnzLj = 0;    nnzUj = 0;    last_drop = SUPERLU_MAX(min_mn - 2 * sp_ienv(7), (int)(min_mn * 0.95));    alpha = pow((double)n, -1.0 / options->ILU_MILU_Dim);    /* Identify relaxed snodes */    relax_end = (int *) intMalloc(n);    relax_fsupc = (int *) intMalloc(n);    if ( options->SymmetricMode == YES )        ilu_heap_relax_snode(n, etree, relax, marker, relax_end, relax_fsupc);    else        ilu_relax_snode(n, etree, relax, marker, relax_end, relax_fsupc);    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;    /* Mark the rows used by relaxed supernodes */    ifill (marker_relax, m, EMPTY);    i = mark_relax(m, relax_end, relax_fsupc, xa_begin, xa_end,                 asub, marker_relax);#if ( PRNTlevel >= 1)    printf("%d relaxed supernodes./n", i);#endif    /*     * 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; ) {        if ( relax_end[jcol] != EMPTY ) { /* start of a relaxed snode */

intpzgstrf_column_dfs(		   const int  pnum,    /* process number */		   const int  m,       /* number of rows in the matrix */		   const int  jcol,    /* current column in the panel */		   const int  fstcol,  /* first column in the panel */		   int *perm_r,   /* row pivotings that are done so far */		   int *ispruned, /* in */		   int *col_lsub, /* the RHS vector to start the dfs */		   int lsub_end,  /* size of col_lsub[] */		   int *super_bnd,/* supernode partition by upper bound */		   int *nseg,     /* modified - with new segments appended */		   int *segrep,   /* modified - with new segments appended */		   int *repfnz,   /* modified */		   int *xprune,   /* modified */		   int *marker2,  /* modified */		   int *parent,   /* working array */		   int *xplore,   /* working array */		   pxgstrf_shared_t *pxgstrf_shared /* modified */		   ){/* * -- SuperLU MT routine (version 2.0) -- * Lawrence Berkeley National Lab, Univ. of California Berkeley, * and Xerox Palo Alto Research Center. * September 10, 2007 * * Purpose * ======= *   pzgstrf_column_dfs() performs a symbolic factorization on column jcol,  *   and detects whether column jcol belongs in the same supernode as jcol-1. * * Local parameters * ================ *   A supernode representative is the last column of a supernode. *   The nonzeros in U[*,j] are segments that end at supernodal *   representatives. The routine returns a list of such supernodal  *   representatives in topological order of the dfs that generates them. *   The location of the first nonzero in each such supernodal segment *   (supernodal entry location) is also returned. * *   nseg: no of segments in current U[*,j] *   samesuper: samesuper=NO if column j does not belong in the same *	        supernode as j-1. Otherwise, samesuper=YES. * *   marker2: A-row --> A-row/col (0/1) *   repfnz: SuperA-col --> PA-row *   parent: SuperA-col --> SuperA-col *   xplore: SuperA-col --> index to L-structure * * Return value * ============ *     0  success; *   > 0  number of bytes allocated when run out of space. * */    GlobalLU_t *Glu = pxgstrf_shared->Glu; /* modified */    Gstat_t *Gstat = pxgstrf_shared->Gstat; /* modified */    register int jcolm1, jcolm1size, nextl, ifrom;    register int k, krep, krow, kperm, samesuper, nsuper;    register int no_lsub;    int	    fsupc;		/* first column in a supernode */    int     myfnz;		/* first nonz column in a U-segment */    int	    chperm, chmark, chrep, kchild;    int     xdfs, maxdfs, kpar;    int     ito;	        /* Used to compress row subscripts */    int     mem_error;    int     *xsup, *xsup_end, *supno, *lsub, *xlsub, *xlsub_end;    static  int  first = 1, maxsuper;    if ( first ) {	maxsuper = sp_ienv(3);	first = 0;    }    /* Initialize pointers */    xsup      = Glu->xsup;    xsup_end  = Glu->xsup_end;    supno     = Glu->supno;    lsub      = Glu->lsub;    xlsub     = Glu->xlsub;    xlsub_end = Glu->xlsub_end;    jcolm1    = jcol - 1;    nextl     = lsub_end;    no_lsub   = 0;    samesuper = YES;    /* Test whether the row structure of column jcol is contained       in that of column jcol-1. */    for (k = 0; k < lsub_end; ++k) {	krow = col_lsub[k];	if ( perm_r[krow] == EMPTY ) { /* krow is in L */	    ++no_lsub;	    if (marker2[krow] != jcolm1) 	        samesuper = NO; /* row subset test */	    marker2[krow] = jcol;	}    }#if ( DEBUGlevel>=2 )//.........这里部分代码省略.........

/*! /brief * * <pre> * Purpose * ======= *   ILU_SCOLUMN_DFS performs a symbolic factorization on column jcol, and *   decide the supernode boundary. * *   This routine does not use numeric values, but only use the RHS *   row indices to start the dfs. * *   A supernode representative is the last column of a supernode. *   The nonzeros in U[*,j] are segments that end at supernodal *   representatives. The routine returns a list of such supernodal *   representatives in topological order of the dfs that generates them. *   The location of the first nonzero in each such supernodal segment *   (supernodal entry location) is also returned. * * Local parameters * ================ *   nseg: no of segments in current U[*,j] *   jsuper: jsuper=EMPTY if column j does not belong to the same *	supernode as j-1. Otherwise, jsuper=nsuper. * *   marker2: A-row --> A-row/col (0/1) *   repfnz: SuperA-col --> PA-row *   parent: SuperA-col --> SuperA-col *   xplore: SuperA-col --> index to L-structure * * Return value * ============ *     0  success; *   > 0  number of bytes allocated when run out of space. * </pre> */intilu_scolumn_dfs(	   const int  m,	 /* in - number of rows in the matrix */	   const int  jcol,	 /* in */	   int	      *perm_r,	 /* in */	   int	      *nseg,	 /* modified - with new segments appended */	   int	      *lsub_col, /* in - defines the RHS vector to start the				    dfs */	   int	      *segrep,	 /* modified - with new segments appended */	   int	      *repfnz,	 /* modified */	   int	      *marker,	 /* modified */	   int	      *parent,	 /* working array */	   int	      *xplore,	 /* working array */	   GlobalLU_t *Glu	 /* modified */	   ){    int     jcolp1, jcolm1, jsuper, nsuper, nextl;    int     k, krep, krow, kmark, kperm;    int     *marker2;		/* Used for small panel LU */    int     fsupc;		/* First column of a snode */    int     myfnz;		/* First nonz column of a U-segment */    int     chperm, chmark, chrep, kchild;    int     xdfs, maxdfs, kpar, oldrep;    int     jptr, jm1ptr;    int     ito, ifrom; 	/* Used to compress row subscripts */    int     mem_error;    int     *xsup, *supno, *lsub, *xlsub;    int     nzlmax;    static  int  first = 1, maxsuper;    xsup    = Glu->xsup;    supno   = Glu->supno;    lsub    = Glu->lsub;    xlsub   = Glu->xlsub;    nzlmax  = Glu->nzlmax;    if ( first ) {	maxsuper = sp_ienv(3);	first = 0;    }    jcolp1  = jcol + 1;    jcolm1  = jcol - 1;    nsuper  = supno[jcol];    jsuper  = nsuper;    nextl   = xlsub[jcol];    marker2 = &marker[2*m];    /* For each nonzero in A[*,jcol] do dfs */    for (k = 0; lsub_col[k] != EMPTY; k++) {	krow = lsub_col[k];	lsub_col[k] = EMPTY;	kmark = marker2[krow];	/* krow was visited before, go to the next nonzero */	if ( kmark == jcol ) continue;	/* For each unmarked nbr krow of jcol	 *	krow is in L: place it in structure of L[*,jcol]	 */	marker2[krow] = jcol;	kperm = perm_r[krow];//.........这里部分代码省略.........

//.........这里部分代码省略.........	}	if (colequ && *info == 0) {	    rcmin = bignum;	    rcmax = 0.;	    for (j = 0; j < A->nrow; ++j) {		rcmin = SUPERLU_MIN(rcmin, C[j]);		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 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||		      B->Stype != SLU_DN || B->Dtype != SLU_D || 		      B->Mtype != SLU_GE )		*info = -13;	    else if ( X->ncol < 0 || Xstore->lda < SUPERLU_MAX(0, A->nrow) ||		      (B->ncol != 0 && B->ncol != X->ncol) ||                      X->Stype != SLU_DN ||		      X->Dtype != SLU_D || X->Mtype != SLU_GE )		*info = -14;	}    }    if (*info != 0) {	i = -(*info);	superlu_xerbla("dgssvx", &i);	return;    }        /* Initialization for factor parameters */    panel_size = sp_ienv(1);    relax      = sp_ienv(2);    drop_tol   = 0.0;    utime = stat->utime;        /* Convert A to SLU_NC format when necessary. */    if ( A->Stype == SLU_NR ) {	NRformat *Astore = (NRformat*) A->Store;	AA = (SuperMatrix *) SUPERLU_MALLOC( sizeof(SuperMatrix) );	dCreate_CompCol_Matrix(AA, A->ncol, A->nrow, Astore->nnz, 			       (double*) 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. */	dgsequ(AA, R, C, &rowcnd, &colcnd, &amax, &info1);		if ( info1 == 0 ) {	    /* Equilibrate matrix A. */	    dlaqgs(AA, R, C, rowcnd, colcnd, amax, equed);

voidpsgstrf_panel_bmod(		   const int  pnum, /* process number */		   const int  m,    /* number of rows in the matrix */		   const int  w,    /* current panel width */		   const int  jcol, /* leading column of the current panel */		   const int  bcol, /* first column of the farthest busy snode*/		   int   *inv_perm_r,/* in; inverse of the row pivoting */		   int   *etree,     /* in */		   int   *nseg,      /* modified */		   int   *segrep,    /* modified */		   int   *repfnz,    /* modified, size n-by-w */		   int   *panel_lsub,/* modified */		   int   *w_lsub_end,/* modified */		   int   *spa_marker,/* modified; size n-by-w */		   float *dense, /* modified, size n-by-w */		   float *tempv, /* working array - zeros on input/output */		   pxgstrf_shared_t *pxgstrf_shared /* modified */		   ){/* * -- SuperLU MT routine (version 2.0) -- * Lawrence Berkeley National Lab, Univ. of California Berkeley, * and Xerox Palo Alto Research Center. * September 10, 2007 * * Purpose * ======= * *    Performs numeric block updates (sup-panel) in topological order. *    It features combined 1D and 2D blocking of the source updating s-node. *    It consists of two steps: *       (1) accumulates updates from "done" s-nodes. *       (2) accumulates updates from "busy" s-nodes. * *    Before entering this routine, the nonzeros of the original A in *    this panel were already copied into the SPA dense[n,w]. * * Updated/Output arguments * ======================== *    L[*,j:j+w-1] and U[*,j:j+w-1] are returned collectively in the *    m-by-w vector dense[*,w]. The locations of nonzeros in L[*,j:j+w-1] *    are given by lsub[*] and U[*,j:j+w-1] by (nseg,segrep,repfnz). * */    GlobalLU_t *Glu = pxgstrf_shared->Glu;  /* modified */    Gstat_t *Gstat = pxgstrf_shared->Gstat; /* modified */    register int j, k, ksub;    register int fsupc, nsupc, nsupr, nrow;    register int kcol, krep, ksupno, dadsupno;    register int jj;	      /* index through each column in the panel */    int          *xsup, *xsup_end, *supno;    int          *lsub, *xlsub, *xlsub_end;    int          *repfnz_col; /* repfnz[] for a column in the panel */    float       *dense_col;  /* dense[] for a column in the panel */    int          *col_marker; /* each column of the spa_marker[*,w] */    int          *col_lsub;   /* each column of the panel_lsub[*,w] */    static   int first = 1, rowblk, colblk;#ifdef PROFILE    double   t1, t2; /* temporary time */#endif    #ifdef PREDICT_OPT        register float pmod, max_child_eft = 0, sum_pmod = 0, min_desc_eft = 0;    register float pmod_eft;    register int   kid, ndesc = 0;#endif    #if ( DEBUGlevel>=2 )    int dbg_addr = 0*m;#endif        if ( first ) {	rowblk   = sp_ienv(4);	colblk   = sp_ienv(5);	first = 0;    }        xsup      = Glu->xsup;    xsup_end  = Glu->xsup_end;    supno     = Glu->supno;    lsub      = Glu->lsub;    xlsub     = Glu->xlsub;    xlsub_end = Glu->xlsub_end;#if ( DEBUGlevel>=2 )    /*if (jcol >= LOCOL && jcol <= HICOL)    check_panel_dfs_list(pnum, "begin", jcol, *nseg, segrep);*/if (jcol == BADPAN)    printf("(%d) Enter psgstrf_panel_bmod() jcol %d,BADCOL %d,dense_col[%d] %.10f/n",	   pnum, jcol, BADCOL, BADROW, dense[dbg_addr+BADROW]);#endif        /* --------------------------------------------------------------------       For each non-busy supernode segment of U[*,jcol] in topological order,       perform sup-panel update.       -------------------------------------------------------------------- */    k = *nseg - 1;    for (ksub = 0; ksub < *nseg; ++ksub) {//.........这里部分代码省略.........

//.........这里部分代码省略.........	}	if (colequ && *info == 0) {	    rcmin = bignum;	    rcmax = 0.;	    for (j = 0; j < A->nrow; ++j) {		rcmin = SUPERLU_MIN(rcmin, C[j]);		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 || Bstore->lda < SUPERLU_MAX(0, A->nrow) ||		      B->Stype != SLU_DN || B->Dtype != SLU_Z || 		      B->Mtype != SLU_GE )		*info = -13;	    else if ( X->ncol < 0 || Xstore->lda < SUPERLU_MAX(0, A->nrow) ||		      (B->ncol != 0 && B->ncol != X->ncol) ||		      X->Stype != SLU_DN ||		      X->Dtype != SLU_Z || X->Mtype != SLU_GE )		*info = -14;	}    }    if (*info != 0) {	i = -(*info);	xerbla_("zgsisx", &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) );	zCreate_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 ) {	register int i, j;	NCformat *Astore = AA->Store;	int nnz = Astore->nnz;	int *colptr = Astore->colptr;	int *rowind = Astore->rowind;	doublecomplex *nzval = (doublecomplex *)Astore->nzval;	int n = AA->nrow;

intdcolumn_dfs(	   const int  m,         /* in - number of rows in the matrix */	   const int  jcol,      /* in */	   int        *perm_r,   /* in */	   int        *nseg,     /* modified - with new segments appended */	   int        *lsub_col, /* in - defines the RHS vector to start the dfs */	   int        *segrep,   /* modified - with new segments appended */	   int        *repfnz,   /* modified */	   int        *xprune,   /* modified */	   int        *marker,   /* modified */	   int        *parent,	 /* working array */	   int        *xplore,   /* working array */	   GlobalLU_t *Glu       /* modified */	   ){/*  * Purpose * ======= *   "column_dfs" performs a symbolic factorization on column jcol, and *   decide the supernode boundary. * *   This routine does not use numeric values, but only use the RHS  *   row indices to start the dfs. * *   A supernode representative is the last column of a supernode. *   The nonzeros in U[*,j] are segments that end at supernodal *   representatives. The routine returns a list of such supernodal  *   representatives in topological order of the dfs that generates them. *   The location of the first nonzero in each such supernodal segment *   (supernodal entry location) is also returned. * * Local parameters * ================ *   nseg: no of segments in current U[*,j] *   jsuper: jsuper=NO if column j does not belong to the same *	supernode as j-1. Otherwise, jsuper=nsuper. * *   marker2: A-row --> A-row/col (0/1) *   repfnz: SuperA-col --> PA-row *   parent: SuperA-col --> SuperA-col *   xplore: SuperA-col --> index to L-structure * * Return value * ============ *     0  success; *   > 0  number of bytes allocated when run out of space. * */    int     jcolp1, jcolm1, jsuper, nsuper, nextl;    int     k, krep, krow, kmark, kperm;    int     *marker2;           /* Used for small panel LU */    int	    fsupc;		/* First column of a snode */    int     myfnz;		/* First nonz column of a U-segment */    int	    chperm, chmark, chrep, kchild;    int     xdfs, maxdfs, kpar, oldrep;    int     jptr, jm1ptr;    int     ito, ifrom, istop;	/* Used to compress row subscripts */    int     mem_error;    int     *xsup, *supno, *lsub, *xlsub;    int     nzlmax;    static  int  first = 1, maxsuper;        xsup    = Glu->xsup;    supno   = Glu->supno;    lsub    = Glu->lsub;    xlsub   = Glu->xlsub;    nzlmax  = Glu->nzlmax;    if ( first ) {	maxsuper = sp_ienv(3);	first = 0;    }    jcolp1  = jcol + 1;    jcolm1  = jcol - 1;    nsuper  = supno[jcol];    jsuper  = nsuper;    nextl   = xlsub[jcol];    marker2 = &marker[2*m];    /* For each nonzero in A[*,jcol] do dfs */    for (k = 0; lsub_col[k] != EMPTY; k++) {	krow = lsub_col[k];    	lsub_col[k] = EMPTY;	kmark = marker2[krow];    		/* krow was visited before, go to the next nonz */        if ( kmark == jcol ) continue; 	/* For each unmarked nbr krow of jcol	 *	krow is in L: place it in structure of L[*,jcol]	 */	marker2[krow] = jcol;	kperm = perm_r[krow];   	if ( kperm == EMPTY ) {	    lsub[nextl++] = krow; 	/* krow is indexed into A */	    if ( nextl >= nzlmax ) {//.........这里部分代码省略.........

//.........这里部分代码省略......... *             Pr*transpose(A)*Pc=L*U   (if A->Stype = SLU_NR). *         Uses column-wise storage scheme, i.e., U has types: *         Stype = SLU_NC, Dtype = SLU_D, Mtype = TRU. * * B       (input/output) SuperMatrix* *         B has types: Stype = SLU_DN, Dtype = SLU_D, Mtype = SLU_GE. *         On entry, the right hand side matrix. *         On exit, the solution matrix if info = 0; * * info    (output) int* *	   = 0: successful exit *         > 0: if info = i, and i is *             <= A->ncol: U(i,i) is exactly zero. The factorization has *                been completed, but the factor U is exactly singular, *                so the solution could not be computed. *             > A->ncol: number of bytes allocated when memory allocation *                failure occurred, plus A->ncol. *    */    double   t1;	/* Temporary time */    char     refact[1], trans[1];    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   diag_pivot_thresh = 1.0;    double   drop_tol = 0;    int      panel_size;     /* panel size */    int      relax;          /* no of columns in a relaxed snodes */    double   *utime;    extern SuperLUStat_t SuperLUStat;    /* Test the input parameters ... */    *info = 0;    Bstore = B->Store;    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 = -1;    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 = -6;    if ( *info != 0 ) {	i = -(*info);	xerbla_("dgssv", &i);	return;    }        *refact = 'N';    *trans = 'N';    panel_size = sp_ienv(1);    relax = sp_ienv(2);    StatInit(panel_size, relax);    utime = SuperLUStat.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 = 'T';    } else if ( A->Stype == SLU_NC ) AA = A;    etree = intMalloc(A->ncol);    t1 = SuperLU_timer_();    sp_preorder(refact, AA, perm_c, etree, &AC);    utime[ETREE] = SuperLU_timer_() - t1;    /*printf("Factor PA = LU ... relax %d/tw %d/tmaxsuper %d/trowblk %d/n", 	  relax, panel_size, sp_ienv(3), sp_ienv(4));*/    t1 = SuperLU_timer_();     /* Compute the LU factorization of A. */    dgstrf(refact, &AC, diag_pivot_thresh, drop_tol, relax, panel_size,	   etree, NULL, lwork, perm_r, perm_c, L, U, info);    utime[FACT] = SuperLU_timer_() - t1;    t1 = SuperLU_timer_();    if ( *info == 0 ) {        /* Solve the system A*X=B, overwriting B with X. */        dgstrs (trans, L, U, perm_r, perm_c, B, info);    }    utime[SOLVE] = SuperLU_timer_() - t1;    SUPERLU_FREE (etree);    Destroy_CompCol_Permuted(&AC);    if ( A->Stype == SLU_NR ) {	Destroy_SuperMatrix_Store(AA);	SUPERLU_FREE(AA);    }    /*PrintStat( &SuperLUStat );*/    StatFree();}