Implementation of SuperLU_MT?

[SalmanMaths]

I am the beginner for SuperLU. I am trying to solve my system A*X= B through C++ coding. In my code A is SK and B is f. I trying to solve this by
X = SK.Solve_superlu(f); where
FEM_Matrix_2<> SK(mesh);
template
std::vector FEM_Matrix_2::Solve_superlu(std::vector const &f) const
{
SuperMatrix A, AC, L, U, B;
NCformat *Astore;
SCPformat *Lstore;
NCPformat *Ustore;
superlumt_options_t superlumt_options;
pxgstrf_shared_t pxgstrf_shared;
pdgstrf_threadarg_t *pdgstrf_threadarg;
int_t nprocs;
fact_t fact;
trans_t trans;
yes_no_t refact, usepr;
double u, drop_tol;
double *a;
int_t *asub, *xa;
int_t perm_c; / column permutation vector */
int_t perm_r; / row permutations from partial pivoting */
void *work;
int_t info, lwork, nrhs, ldx;
int_t m, n, nnz, permc_spec, panel_size, relax;
int_t i, firstfact;
double *rhsb, *xact;
Gstat_t Gstat;
flops_t flopcnt;
int argc; char *argv; //int_t nprocs;
void parse_command_line(int argc, char *argv[], int_t *nprocs);

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

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

//      Initialize matrix A.
int m1 = f.size();

m = n = f.size();
nnz = _sk.size();

if ( !(a = doubleMalloc(nnz)) ) SUPERLU_ABORT("Malloc fails for a[].");
if ( !(asub = intMalloc(nnz)) ) SUPERLU_ABORT("Malloc fails for asub[].");
if ( !(xa = intMalloc(n+1)) ) SUPERLU_ABORT("Malloc fails for xa[].");

dCompRow_to_CompCol(m, n, nnz,
                    const_cast<double*>(_sk.data()), const_cast<int_t*>(_ik.data()), const_cast<int_t*>(_id.data()),
                    &a, &asub, &xa);
           
/* Read the input matrix stored in Harwell-Boeing format. */                  
//dreadhb(&m, &n, &nnz, &a, &asub, &xa);

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


// Create right-hand side matrix B.
//nrhs = 1;
if (!(rhsb = doubleMalloc(m * nrhs))) SUPERLU_ABORT("Malloc fails for rhsb[].");

for (int_t j = 0; j < m; ++j)
{
rhsb[j] = f[j];
}
dCreate_Dense_Matrix(&B, m, nrhs, rhsb, m, SLU_DN, SLU_D, SLU_GE);

//xact = doubleMalloc(n * nrhs);

// ldx = n;
//dGenXtrue(n, nrhs, xact, ldx);
//dFillRHS(trans, nrhs, xact, ldx, &A, &B);

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

SUPERLU_ABORT("Malloc fails for etree[].");
if ( !(superlumt_options.colcnt_h = intMalloc(n)) )
SUPERLU_ABORT("Malloc fails for colcnt_h[].");
if ( !(superlumt_options.part_super_h = intMalloc(n)) )
SUPERLU_ABORT("Malloc fails for colcnt_h[].");

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

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


/* ------------------------------------------------------------
   Get column permutation vector perm_c[], according to permc_spec:
   permc_spec = 0: natural ordering
   permc_spec = 1: minimum degree ordering on structure of A'*A
   permc_spec = 2: minimum degree ordering on structure of A'+A
   permc_spec = 3: approximate minimum degree for unsymmetric matrices
   ------------------------------------------------------------*/
permc_spec = 1;
get_perm_c(permc_spec, &A, perm_c);
/* ------------------------------------------------------------
   Initialize the option structure superlumt_options using the
   user-input parameters;
   Apply perm_c to the columns of original A to form AC.
   ------------------------------------------------------------*/
refact= NO;
pdgstrf_init(nprocs, fact, trans, refact, panel_size, relax,
             u, usepr, drop_tol, perm_c, perm_r,
             work, lwork, &A, &AC, &superlumt_options, &Gstat);
/* ------------------------------------------------------------
   Compute the LU factorization of A.
   The following routine will create nprocs threads.
   ------------------------------------------------------------*/
pdgstrf(&superlumt_options, &AC, perm_r, &L, &U, &Gstat, &info);
flopcnt = 0;
for (i = 0; i < nprocs; ++i) flopcnt += Gstat.procstat[i].fcops;
Gstat.ops[FACT] = flopcnt;
/* ------------------------------------------------------------
   Solve the system A*X=B, overwriting B with X.
   ------------------------------------------------------------*/
dgstrs(trans, &L, &U, perm_r, perm_c, &B, &Gstat, &info);

printf("\n** Result of sparse LU **\n");
std::vector<double> sol(f.size());
for(size_t j=0; i<sol.size(); ++j)
{
    NCformat* p = static_cast<NCformat*>(B.Store);
    double* val = static_cast<double*>(p->nzval);
    sol[j]=val[j];
   
}
//dinf_norm_error(nrhs, &B, xact); /* Check inf. norm of the error */

Destroy_CompCol_Permuted(&AC); /* Free extra arrays in AC. */

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

/* ------------------------------------------------------------
   Re-initialize statistics variables and options used by the
   factorization routine pdgstrf().
   ------------------------------------------------------------*/
StatInit(n, nprocs, &Gstat);
refact= YES;
pdgstrf_init(nprocs, fact, trans, refact, panel_size, relax,

u, usepr, drop_tol, perm_c, perm_r,
work, lwork, &A, &AC, &superlumt_options, &Gstat);

/* ------------------------------------------------------------
   Compute the LU factorization of A.
   The following routine will create nprocs threads.
   ------------------------------------------------------------*/
pdgstrf(&superlumt_options, &AC, perm_r, &L, &U, &Gstat, &info);

flopcnt = 0;
for (i = 0; i < nprocs; ++i) flopcnt += Gstat.procstat[i].fcops;
Gstat.ops[FACT] = flopcnt;

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


 /* ------------------------------------------------------------
   Deallocate storage after factorization.
   ------------------------------------------------------------*/
pxgstrf_finalize(&superlumt_options, &AC);

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

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

SUPERLU_FREE (rhsb);
SUPERLU_FREE (xact);
SUPERLU_FREE (perm_r);
SUPERLU_FREE (perm_c);
SUPERLU_FREE (superlumt_options.etree);
SUPERLU_FREE (superlumt_options.colcnt_h);
SUPERLU_FREE (superlumt_options.part_super_h);
Destroy_CompCol_Matrix(&A);
Destroy_SuperMatrix_Store(&B);
if ( lwork == 0 ) {
    Destroy_SuperNode_SCP(&L);
    Destroy_CompCol_NCP(&U);
} else if ( lwork > 0 ) {
    SUPERLU_FREE(work);
}
StatFree(&Gstat);

return sol;

}

In this way, the SuperLU is working but SuperLU_MT (the current code is not working). The following error gives me:
** On entry to sp_ienv, parameter number 1 had an illegal value
Storage for L subscripts exceeded; Current column 0; Need at least 72;
You may set it by the 8-th parameter in routine sp_ienv().
Memory allocation failed at line 222 in file pmemory.c

Please help me with this, what is the problem?
Best Regards,
Salman Ahmad

[xiaoyeli]

It seems your set up is similar to the example EXAMPLE/pdrrepeat.c.
Are you able to run that example?

[SalmanMaths]

Dear Sherry Li, Thank you so much for your reply. Yes I took this idea from the EXAMPLE/pdrepeat.c. The only difference is just defining the matrix and right hand side. That program runs but when I change the program to my problem then they give me an error. My problem is time dependent and nonlinear, so I repeat the iterations (time iteration by backward euler method as well as newton iteration for nonlinearity). So I need this approach to speed up my program. Best Regards, Salman AHmad

…

On Sun, Feb 26, 2023 at 5:14 PM X. Sherry Li ***@***.***> wrote: It seems your set up is similar to the example EXAMPLE/pdrrepeat.c. Are you able to run that example? — Reply to this email directly, view it on GitHub <#133 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AXKALTRYRMKGVQJP5ULZLR3WZMNH3ANCNFSM6AAAAAAVIKHWTU> . You are receiving this because you authored the thread.Message ID: ***@***.***>

[SalmanMaths]

Dear Sherry Li,

When, I run the dreadhb(&m, &n, &nnz, &a, &asub, &xa); function in SuperLU_MT, they gives me the error: Segmentation fault (core dumped).
I am using the C++17. Is SuperLU_MT compatible with C++17?

Best Regards,
Salman Ahmad