pnslac.c
/*
Use PARPACK + SuperLU_DIST to solve generalized eigenvalue
arising from Accelerator Design (SLAC problem)
This version uses the nonsymmetric solver of ARPACK: pdnaupd.
NOTE: Aztec is used for sparse matrix-vector multiplication.
*/
#include <math.h>
#include "superlu_ddefs.h"
#define workd(i) workd[(i)-1]
#define z(i,j) z[nloc*( (j) -1 )+ (i) - 1]
#define zglb(i,j) zglb[na*( (j) -1 )+ (i) - 1]
#define asub(i) asub[(i) - 1]
#define amtx(i) amtx[(i) - 1]
#define USECRS
#define USEAZTEC
#ifdef USEAZTEC
#include "az_aztec.h"
#endif
#ifdef USECRS
/* Convert csr format to msr format */
void csr2msr(int n,int nnz,int *row_start, int *col_index, double *val,
int *bindex,double *newval,int *update)
{
int row,j,count,offdiags;
/* Copy val and col_index to second part of bindex and newval */
bindex[0]=n+1;
count = n+1;
offdiags=0;
for(row=0; row < n;row++) {
if(row > 0) {
bindex[row]=bindex[row-1]+offdiags;
offdiags=0;
}
for(j=row_start[row];j < row_start[row+1];j++) {
if(col_index[j] == update[row]) {
newval[row]=val[j]; /* Diagonal Element */
} else {
newval[count]=val[j];
bindex[count]=col_index[j];
count++;
offdiags++;
}
}
}
bindex[n]=count;
}
/*
Convert a global compressed column matrix to
local compressed row matrix
Input matrix is in ncols, col_start, row_index, col_val
Local row range is given by my_starting_row <= row < my_finishing_row
Matrix is output in p_row_start, p_col_index, and p_val
Returns number of local entries
*/
int gccs2lcrs(int ncols,int *col_start, int *row_index, double *col_val,
int my_starting_row, int my_finishing_row,
int **p_row_start, int **p_col_index, double **p_val)
{
int *row_start;
int *col_index;
double *val;
int n_local_rows = my_finishing_row - my_starting_row;
int *row_length;
int *seen_in_row;
int total_entries, gcol,i,current_row,lrow;
total_entries = 0;
row_length = (int *) calloc(n_local_rows,sizeof(int));
for(gcol=0;gcol < ncols;gcol++) {
for(i=col_start[gcol];i < col_start[gcol+1];i++) {
current_row = row_index[i];
if((current_row >= my_starting_row) &&
(current_row < my_finishing_row)) {
row_length[current_row-my_starting_row]++;
total_entries++;
}
}
}
printf("ncols: %d\n",ncols);
printf("n_local_rows: %d\n",n_local_rows);
printf("Total entries: %d\n",total_entries);
*p_row_start = row_start = (int *) calloc(n_local_rows+1,sizeof(int));
*p_col_index = col_index = (int *) calloc(total_entries,sizeof(int));
*p_val = val = (double *) calloc(total_entries,sizeof(double));
row_start[0]=0;
for(lrow=1;lrow <= n_local_rows;lrow++) {
row_start[lrow]=row_start[lrow-1]+row_length[lrow-1];
}
seen_in_row = (int *) calloc(n_local_rows,sizeof(int));
for(lrow=0;lrow < n_local_rows;lrow++) {
seen_in_row[lrow]=row_start[lrow];
}
for(gcol=0;gcol < ncols;gcol++) {
for(i=col_start[gcol];i < col_start[gcol+1];i++) {
current_row = row_index[i];
if((current_row >= my_starting_row) &&
(current_row < my_finishing_row)) {
lrow = current_row - my_starting_row;
col_index[seen_in_row[lrow]]=gcol;
val[seen_in_row[lrow]]=col_val[i];
seen_in_row[lrow]++;
}
}
}
free(row_length);
free(seen_in_row);
return(total_entries);
}
void crsmvm(int n,int *row_start, int *col_index, double *val,
double *x,double *y)
{
int i,j;
for(i=0;i < n;i++) {
y[i]=0.0;
for(j=row_start[i];j < row_start[i+1];j++) {
y[i]+=val[j]*x[col_index[j]];
}
}
}
#endif
main(int argc, char *argv[])
/*
* Purpose
* =======
*
* The driver program PEIGEN
*
* This example illustrates how to use SuperLU_DIST and PARPACK
* to solve a large-scale sparse standard eigenvalue problem
* in shift-invert mode.
* In this case, we factorize A only once in the first call to
* pdgssvx_ABglobal, and reuse the following data structures
* in the subsequent call to pdgssvx_ABglobal:
* ScalePermstruct : DiagScale, R, C, perm_r, perm_c
* LUstruct : Glu_persist, Llu
*
* On the Cray T3E, the program may be run by typing:
* mpprun -n <procs> pddrive1 -r <proc rows> -c <proc columns> <input_matrix>
*
*/
{
superlu_options_t options;
SuperLUStat_t stat;
SuperMatrix A;
ScalePermstruct_t ScalePermstruct;
LUstruct_t LUstruct;
gridinfo_t grid;
double *berr;
double *amtx, *bmtx, *b, *b1, *b2;
int_t *asub, *xa, *bsub, *xb;
int_t i, j, ma, na, nnza, mb, nb, nnzb;
int_t nprow, npcol;
int iam, info, ldb, ldx, nrhs;
char trans[1];
char **cpp, c;
FILE *fp1, *fp2, *fopen();
double t0, t1;
/* ----------------
ARPACK variables
----------------*/
MPI_Comm arcomm;
MPI_Status mpistat;
int nprocs, ibeg;
int ido, nloc, nrem, nev , ncv , lworkl, mypid,
sid, arflg, pin, pout, maxitr, rvec , nconv,
iopt;
double tol, sigmar =2000.0, sigmai=0.0;
char bmat[1] = "G", which[2]="LR", howmny[1]="A";
int ipntr[14], iparam[11], *select, *psize, *nbase;
double *v , *z , *workd, *workl, *workev, *resid, *recvbuf,
*dr, *di, *zglb , *res;
int ione = 1, stiff = 1, irow, ibegin, iend;
double done = 1.0, dzero = 0.0;
/* make a copy of the matrix for residual checking */
int *colptr, *rowind;
double *nzvals;
#ifdef USECRS
int *row_start;
int *col_index;
double *val;
#ifdef USEAZTEC
int proc_config[AZ_PROC_SIZE];
AZ_MATRIX *Bmat;
double *dsrc;
double *ddest,*ddest2;
int N_update;
int *external, *update_index, *update, *extern_index, *data_org;
int aztotal;
int *bindx;
double *msrval;
int local_nnz;
#endif
#endif
nprow = 1; /* Default process rows. */
npcol = 1; /* Default process columns. */
nrhs = 1; /* Number of right-hand side. */
/* Parse command line argv[]. */
if (argc != 7) ABORT("Number of arguements incorrect");
for (cpp = argv+1; *cpp; ++cpp) {
if ( **cpp == '-' ) {
c = *(*cpp+1);
++cpp;
switch (c) {
case 'h':
printf("Options:\n");
printf("\t-r <int>: process rows (default %d)\n", nprow);
printf("\t-c <int>: process columns (default %d)\n", npcol);
exit(0);
break;
case 'r': nprow = atoi(*cpp);
break;
case 'c': npcol = atoi(*cpp);
break;
}
} else { /* Last two args are considered a filename */
if (stiff) {
/* open the stiffness matrix */
if ( !(fp1 = fopen(*cpp, "r")) ) {
ABORT("File does not exist");
}
stiff = 0;
}
else {
if ( !(fp2 = fopen(*cpp, "r")) ) {
ABORT("File does not exist");
}
break;
}
}
}
/* ------------------------------------------------------------
INITIALIZE MPI ENVIRONMENT.
------------------------------------------------------------*/
MPI_Init( &argc, &argv );
t0 = MPI_Wtime();
#ifdef USEAZTEC
AZ_set_proc_config(proc_config,MPI_COMM_WORLD);
#endif
/* ------------------------------------------------------------
INITIALIZE THE SUPERLU PROCESS GRID.
------------------------------------------------------------*/
superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, &grid);
/* Bail out if I do not belong in the grid. */
iam = grid.iam;
if ( iam >= nprow * npcol )
goto out;
#if ( DEBUGlevel>=1 )
CHECK_MALLOC(iam, "Enter main()");
#endif
/* ------------------------------------------------------------
PROCESS 0 READS THE MATRIX A, AND THEN BROADCASTS IT TO ALL
THE OTHER PROCESSES.
------------------------------------------------------------*/
if ( !iam ) {
/* Print the CPP definitions. */
cpp_defs();
/* Read the stiffness matrix stored on disk in binary format. */
fread(&ma , 4, 1, fp1);
na = ma;
fread(&nnza, 4, 1, fp1);
xa = (int*)calloc(ma+1, sizeof(int));
asub = (int*)calloc(nnza, sizeof(int));
amtx = (double*)calloc(nnza, sizeof(double));
fread((int*)xa , 4, na+1, fp1);
fread((int*)asub , 4, nnza, fp1);
fread((double*)amtx, 8, nnza, fp1);
fclose(fp1);
printf("\tA: Dimension\t%dx%d\t # nonzeros %d\n", ma, na, nnza);
/* Change to 0-based index */
for (i=0;i<na+1;i++) xa[i]--;
for (i=0;i<nnza;i++) asub[i]--;
/* Read the mass matrix */
fread(&mb , 4, 1, fp2);
nb = mb;
fread(&nnzb, 4, 1, fp2);
xb = (int*)calloc(mb+1, sizeof(int));
bsub = (int*)calloc(nnzb, sizeof(int));
bmtx = (double*)calloc(nnzb, sizeof(double));
fread((int*)xb , 4, nb+1, fp2);
fread((int*)bsub , 4, nnzb, fp2);
fread((double*)bmtx, 8, nnzb, fp2);
fclose(fp2);
/* Change to 0-based index */
for (i=0;i<nb+1;i++) xb[i]--;
for (i=0;i<nnzb;i++) bsub[i]--;
printf("\tB: Dimension\t%dx%d\t # nonzeros %d\n", mb, nb, nnzb);
printf("\tProcess grid\t%d X %d\n", grid.nprow, grid.npcol);
/* Broadcast matrix A to the other PEs. */
MPI_Bcast( &ma, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &na, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nnza, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( amtx, nnza, MPI_DOUBLE, 0, grid.comm );
MPI_Bcast( asub, nnza, mpi_int_t, 0, grid.comm );
MPI_Bcast( xa, na+1, mpi_int_t, 0, grid.comm );
/* Broadcast matrix B to the other PEs. */
MPI_Bcast( &mb, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nb, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nnzb, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( bmtx, nnzb, MPI_DOUBLE, 0, grid.comm );
MPI_Bcast( bsub, nnzb, mpi_int_t, 0, grid.comm );
MPI_Bcast( xb, nb+1, mpi_int_t, 0, grid.comm );
} else {
/* Receive matrix A from PE 0. */
MPI_Bcast( &ma, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &na, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nnza, 1, mpi_int_t, 0, grid.comm );
/* Allocate storage for compressed column representation. */
dallocateA_dist(na, nnza, &amtx, &asub, &xa);
MPI_Bcast( amtx, nnza, MPI_DOUBLE, 0, grid.comm );
MPI_Bcast( asub, nnza, mpi_int_t, 0, grid.comm );
MPI_Bcast( xa, na+1, mpi_int_t, 0, grid.comm );
/* Receive matrix B from PE 0. */
MPI_Bcast( &mb, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nb, 1, mpi_int_t, 0, grid.comm );
MPI_Bcast( &nnzb, 1, mpi_int_t, 0, grid.comm );
/* Allocate storage for compressed column representation. */
dallocateA_dist(nb, nnzb, &bmtx, &bsub, &xb);
MPI_Bcast( bmtx, nnzb, MPI_DOUBLE, 0, grid.comm );
MPI_Bcast( bsub, nnzb, mpi_int_t, 0, grid.comm );
MPI_Bcast( xb, nb+1, mpi_int_t, 0, grid.comm );
}
/* Make a copy of the stiffness matrix A (amtx will be corrupted by fact)*/
colptr = (int*)calloc((na+1),sizeof(int));
rowind = (int*)calloc(nnza,sizeof(int));
nzvals = (double*)calloc(nnza,sizeof(double));
for (i=0;i<=na;i++) colptr[i]=xa[i]+1;
for (i=0;i<nnza;i++) {
rowind[i]=asub[i]+1;
nzvals[i]=amtx[i];
}
/* A <--- A - sigmar B */
for (i=0;i<nnza;i++) {
amtx[i] = amtx[i] - sigmar*bmtx[i];
}
/* Create compressed column matrix for A. */
dCreate_CompCol_Matrix_dist(&A, ma, na, nnza, amtx, asub, xa, NC, D, GE);
/* Generate the exact solution and compute the right-hand side. */
if ( !(b = doubleMalloc_dist(ma * nrhs)) ) ABORT("Malloc fails for b[]");
if ( !(b1 = doubleMalloc_dist(ma * nrhs)) ) ABORT("Malloc fails for b1[]");
if ( !(b2 = doubleMalloc_dist(ma * nrhs)) ) ABORT("Malloc fails for b2[]");
*trans = 'N';
ldx = na;
ldb = ma;
if ( !(berr = doubleMalloc_dist(nrhs)) )
ABORT("Malloc fails for berr[].");
/* ------------------------------------------------------------
WE SOLVE THE LINEAR SYSTEM FOR THE FIRST TIME.
------------------------------------------------------------*/
/* Set the default input options. */
set_default_options(&options);
options.RowPerm = NO; /* XSL 4-11-02 */
/* Initialize ScalePermstruct and LUstruct. */
ScalePermstructInit(ma, na, &ScalePermstruct);
LUstructInit(ma, na, &LUstruct);
/* Initialize the statistics variables. */
PStatInit(&stat);
/* Call the linear equation solver: just to factorize when nrhs = 0. */
nrhs = 0;
pdgssvx_ABglobal(&options, &A, &ScalePermstruct, b, ldb, nrhs, &grid,
&LUstruct, berr, &stat, &info);
PStatPrint(&stat, &grid); /* Print the statistics. */
PStatFree(&stat);
/* ------------------------------------------------------------
NOW WE CALL pdnaupd to solve the eigenvalue problem
WE WILL USE THE EXISTING L AND U FACTORS IN
LUSTRUCT OBTAINED FROM A PREVIOUS FATORIZATION.
------------------------------------------------------------*/
options.Fact = FACTORED; /* Indicate the factored form of A is supplied. */
options.IterRefine = NO; /* XSL 4-9-02 */
PStatInit(&stat); /* Initialize the statistics variables. */
if (mypid == 0) {
printf(" nnza = %d\n", nnza);
printf(" xa[na] = %d\n", xa[na]);
printf(" asub[0] = %d\n", asub[0]);
printf(" asub[nnza-1] = %d\n", asub[nnza-1]);
}
arcomm = MPI_COMM_WORLD;
MPI_Comm_rank(arcomm, &mypid);
nev = 20;
ncv = 40;
lworkl = 3*ncv*ncv+6*ncv;
/* distribute ARPACK work space as uniformly as possible */
nprocs = nprow*npcol;
psize = (int*)malloc(nprocs*sizeof(int));
nbase = (int*)malloc(nprocs*sizeof(int));
if (mypid == 0) {
nloc = na/nprocs;
nrem = na - nloc*nprocs;
for (i = 0; i<nprocs; i++) {
psize[i] = nloc;
if (i<nrem) psize[i]++;
}
nbase[0] = 0;
for (i = 1; i<nprocs; i++)
nbase[i] = nbase[i-1]+psize[i-1];
recvbuf = (double*)malloc(nloc*sizeof(double));
}
if (mypid == 0) printf(" nprocs = %d\n", nprocs);
MPI_Bcast( psize, nprocs, MPI_INT, 0, arcomm );
MPI_Bcast( nbase, nprocs, MPI_INT, 0, arcomm );
nloc = psize[mypid];
#ifdef USECRS
#ifdef USEAZTEC
local_nnz = gccs2lcrs(nb,xb,bsub,bmtx,nbase[mypid],nbase[mypid]+nloc,
&row_start,&col_index,&val);
#else
gccs2lcrs(nb,xb,bsub,bmtx,0,nb,&row_start,&col_index,&val);
#endif
#endif
/* Change the indices of the mass matrix to fortran 1-based indices */
for (i=0;i<=nb;i++) xb[i]++;
for (i=0;i<nnzb;i++) bsub[i]++;
v = (double*) malloc( nloc*ncv*sizeof(double) );
z = (double*) malloc( nloc*ncv*sizeof(double) );
workl = (double*) malloc( lworkl*sizeof(double) );
workev = (double*) malloc( 3*nloc*sizeof(double) );
resid = (double*) malloc( nloc*sizeof(double) );
dr = (double*) malloc( ncv*sizeof(double) );
di = (double*) malloc( ncv*sizeof(double) );
select = (int*) malloc( ncv*sizeof(int));
for (i=0;i<ncv;i++) {
dr[i] = 0.0;
di[i] = 0.0;
}
workd = (double*) malloc( 3*nloc*sizeof(double) );
tol = 1.0e-6;
maxitr = 100;
iparam[0] = 1;
iparam[2] = maxitr;
iparam[6] = 3;
printf(" mypid = %d, nloc = %d\n", mypid, nloc);
printf(" mypid = %d, nbase = %d\n", mypid, nbase[mypid]);
#ifdef USEAZTEC
update = calloc(nloc,sizeof(int));
N_update = nloc;
for(i=0;i < nloc;i++) {
update[i]=nbase[mypid]+i;
}
bindx = calloc(local_nnz+N_update+1,sizeof(int));
msrval = calloc(local_nnz+N_update+1,sizeof(double));
csr2msr(N_update,local_nnz,row_start,col_index,val,bindx,msrval,update);
free(row_start);
free(col_index);
free(val);
AZ_transform(proc_config, &external, bindx, msrval, update, &update_index,
&extern_index, &data_org, N_update, NULL, NULL, NULL, NULL,
AZ_MSR_MATRIX);
Bmat = AZ_matrix_create(data_org[AZ_N_internal]+data_org[AZ_N_border]);
AZ_set_MSR(Bmat,bindx,msrval,data_org,0,NULL,AZ_LOCAL);
aztotal= data_org[AZ_N_internal]+data_org[AZ_N_border]+
data_org[AZ_N_external] + 1;
dsrc = calloc(aztotal,sizeof(double));
ddest = calloc(aztotal,sizeof(double));
ddest2 = calloc(aztotal,sizeof(double));
#endif
ido = 0;
arflg = 0;
fprintf(stderr, " calling pdnaupd ...\n");
while (ido != 99) {
pdnaupd(&arcomm, &ido , bmat, &nloc, which, &nev , &tol ,
resid , &ncv , v , &nloc, iparam, ipntr, workd,
workl , &lworkl, &arflg);
if (ido == -1) {
/* collect the local pieces from all processors
and put it in b1 */
pin = ipntr[0];
pout = ipntr[1];
MPI_Allgatherv(&workd(pin),nloc,MPI_DOUBLE,b1,psize,nbase,
MPI_DOUBLE,arcomm);
/* do b2 <--- B*b1 */
#ifdef USECRS
#ifdef USEAZTEC
/* Copy source */
bzero(dsrc,aztotal*sizeof(double));
dcopy(&nloc, &workd(pin), &ione, dsrc, &ione);
AZ_reorder_vec(dsrc,data_org,update_index,NULL);
Bmat->matvec(dsrc,ddest,Bmat,proc_config);
AZ_invorder_vec(ddest,data_org,update_index,NULL,ddest2);
MPI_Allgatherv(ddest2,nloc,MPI_DOUBLE,b2,psize,nbase,
MPI_DOUBLE,arcomm);
#else
crsmvm(nb,row_start,col_index,val,b1,b2);
#endif
#else
iopt = 1;
dmvm( &nb, bmtx, bsub, xb, b1, b2, &iopt );
#endif
nrhs = 1;
pdgssvx_ABglobal(&options, &A , &ScalePermstruct,
b2 , ldb , nrhs ,
&grid , &LUstruct, berr ,
&stat , &info);
/* put the solution into workd */
ibeg = nbase[mypid];
dcopy(&nloc, &b2[ibeg], &ione, &workd(pout), &ione);
}
else if (ido == 1) {
/* collect the local pieces from all processors
and put it in b1 */
pin = ipntr[2];
pout = ipntr[1];
MPI_Allgatherv(&workd(pin),nloc,MPI_DOUBLE,b1,psize,nbase,
MPI_DOUBLE,arcomm);
nrhs = 1;
pdgssvx_ABglobal(&options, &A , &ScalePermstruct,
b1 , ldb , nrhs ,
&grid , &LUstruct, berr ,
&stat , &info);
/* put the solution into workd */
ibeg = nbase[mypid];
dcopy(&nloc, &b1[ibeg], &ione, &workd(pout), &ione);
}
else if (ido == 2) {
/* collect the local pieces from all processors
and put it in b1 */
pin = ipntr[0];
pout = ipntr[1];
MPI_Allgatherv(&workd(pin),nloc,MPI_DOUBLE,b1,psize,nbase,
MPI_DOUBLE,arcomm);
#ifdef USECRS
#ifdef USEAZTEC
bzero(dsrc,aztotal*sizeof(double));
dcopy(&nloc, &workd(pin), &ione, dsrc, &ione);
AZ_reorder_vec(dsrc,data_org,update_index,NULL);
Bmat->matvec(dsrc,ddest,Bmat,proc_config);
AZ_invorder_vec(ddest,data_org,update_index,NULL,ddest2);
MPI_Allgatherv(ddest2,nloc,MPI_DOUBLE,b2,psize,nbase,
MPI_DOUBLE,arcomm);
#else
crsmvm(nb,row_start,col_index,val,b1,b2);
#endif
#else
iopt = 1;
dmvm( &nb, bmtx, bsub, xb, b1, b2, &iopt);
#endif
/* put the product into workd */
ibeg = nbase[mypid];
dcopy(&nloc, &b2[ibeg], &ione, &workd(pout), &ione);
}
}
/* DEALLOCATE SOME STORAGE */
Destroy_CompCol_Matrix_dist(&A);
Destroy_LU(na, &grid, &LUstruct);
ScalePermstructFree(&ScalePermstruct);
LUstructFree(&LUstruct);
if (arflg == 0) {
rvec = 1;
nconv = iparam[4];
pdneupd(&arcomm , &rvec , howmny , select , dr , di ,
z , &nloc , &sigmar, &sigmai, workev, bmat ,
&nloc , which , &nev , &tol , resid , &ncv ,
v , &nloc , iparam , ipntr , workd , workl,
&lworkl, &arflg);
if (arflg==0) {
/* compute residual */
if (!iam ) printf(" nconv = %d\n", nconv);
/* Check the accuracy of the solution. */
if (mypid == 0) {
zglb = (double*) malloc(na*nconv*sizeof(double));
res = (double*) malloc(nconv*sizeof(double));
}
for (i=0;i<nconv;i++) {
if (mypid == 0) {
dcopy(&nloc, &z(1,i+1), &ione, &zglb(1,i+1), &ione);
for (j=1;j<nprocs;j++) {
MPI_Recv(recvbuf, psize[j] , MPI_DOUBLE,
j , MPI_ANY_TAG, arcomm ,
&mpistat);
sid = mpistat.MPI_SOURCE;
ibeg = nbase[sid];
dcopy(&psize[j], recvbuf, &ione, &zglb(ibeg+1,i+1), &ione);
}
}
else {
MPI_Send(&z(1,i+1), nloc , MPI_DOUBLE, 0,
mypid , arcomm);
}
} /* end for */
if (mypid == 0) {
dngeres(&na , &nconv, colptr, rowind, nzvals, xb,
bsub, bmtx , dr , di , zglb , &na,
res);
for (i = 0; i<nconv; i++)
printf( "dr = %12.5e, di = %12.5e, res = %11.4e\n",
dr[i], di[i], res[i]);
free(zglb);
free(res);
}
} /* end if (arflg) */
}
/* Print the statistics. */
PStatPrint(&stat, &grid);
t1 = MPI_Wtime();
t1 = t1 - t0;
if (!iam) printf("Total run time = %11.3e\n", t1);
/* ------------------------------------------------------------
DEALLOCATE STORAGE.
------------------------------------------------------------*/
PStatFree(&stat);
SUPERLU_FREE(b);
SUPERLU_FREE(b1);
SUPERLU_FREE(b2);
SUPERLU_FREE(berr);
free(v);
free(z);
free(workl);
free(workev);
free(workd);
free(resid);
free(psize);
free(nbase);
free(dr);
free(di);
free(external);
free(update);
free(update_index);
free(extern_index);
free(data_org);
free(dsrc);
free(ddest);
free(ddest2);
free(bindx);
free(msrval);
AZ_free(Bmat);
if (mypid == 0) free(recvbuf);
/* ------------------------------------------------------------
RELEASE THE SUPERLU PROCESS GRID.
------------------------------------------------------------*/
out:
superlu_gridexit(&grid);
/* ------------------------------------------------------------
TERMINATES THE MPI EXECUTION ENVIRONMENT.
------------------------------------------------------------*/
MPI_Finalize();
#if ( DEBUGlevel>=1 )
CHECK_MALLOC(iam, "Exit main()");
#endif
}
int cpp_defs()
{
printf(".. CPP definitions:\n");
#if ( PRNTlevel>=1 )
printf("\tPRNTlevel = %d\n", PRNTlevel);
#endif
#if ( DEBUGlevel>=1 )
printf("\tDEBUGlevel = %d\n", DEBUGlevel);
#endif
#if ( PROFlevel>=1 )
printf("\tPROFlevel = %d\n", PROFlevel);
#endif
#if ( StaticPivot>=1 )
printf("\tStaticPivot = %d\n", StaticPivot);
#endif
#ifdef USECRS
#ifdef USEAZTEC
printf("\tAztec Matrix Multiplication\n");
#else
printf("\tCRS Matrix Multiplication\n");
#endif
#endif
printf("....\n");
}
syntax highlighted by Code2HTML, v. 0.9.1