Using the Dense matrix computation package wrong. about strumpack HOT 4 CLOSED

Hi jerome92,

Are you using the STRUMPACK-Dense-1.1.1 code, as found here:
http://portal.nersc.gov/project/sparse/strumpack/STRUMPACK-Dense-1.1.1.tar.gz

Or are you using the newer STRUMPACK v2.x.y?

Perhaps send us a piece of code so we can help you figure out what is wrong.

Pieter

from strumpack.

I used the STRUMPACK-Dense-1.1.1 code, and modify some in solve.cpp, the code as follow:
`#include "StrumpackDensePackage.hpp"

#define myscalar double
#define myreal double

int main (int argc, char *argv[]) {
myscalar *A=NULL, *X=NULL, *B=NULL;
int descA[BLACSCTXTSIZE], descXB[BLACSCTXTSIZE];
int n;
int nrhs;
int nb;
int locr, locc;
int i, j, ii, jj;
int ierr;
int dummy;
int myid, np;
int myrow, mycol, nprow, npcol;
int ctxt;
FILE *cp;

n=3; /* Size of the problem /
nrhs=1; / Number of RHS /
nb=16; / Blocksize for the 2D block-cyclic distribution */

/* Initialize MPI */
if((ierr=MPI_Init(&argc,&argv)))
return 1;
myid=-1;
if((ierr=MPI_Comm_rank(MPI_COMM_WORLD,&myid)))
return 1;
np=-1;
if((ierr=MPI_Comm_size(MPI_COMM_WORLD,&np)))
return 1;

/* Initialize the BLACS grid */
nprow=floor(sqrt((float)np));
npcol=np/nprow;
blacs_get_(&IZERO,&IZERO,&ctxt);
blacs_gridinit_(&ctxt,"R",&nprow,&npcol);
blacs_gridinfo_(&ctxt,&nprow,&npcol,&myrow,&mycol);

/* A is a dense n x n distributed Toeplitz matrix /
if(myid<nprownpcol) {
locr=numroc_(&n,&nb,&myrow,&IZERO,&nprow);
locc=numroc_(&n,&nb,&mycol,&IZERO,&npcol);
A=new myscalar[locr*locc];
char file[256];
sprintf(file,"a_%d.txt",myid);
cp=fopen(file,"w+");
dummy=std::max(1,locr);
descinit_(descA,&n,&n,&nb,&nb,&IZERO,&IZERO,&ctxt,&dummy,&ierr);

for(i=1;i<=locr;i++)
{
  for(j=1;j<=locc;j++) {
    ii=indxl2g_(&i,&nb,&myrow,&IZERO,&nprow);
    jj=indxl2g_(&j,&nb,&mycol,&IZERO,&npcol);
    // Toeplitz matrix from Quantum Chemistry.
    // myreal pi=3.1416, d=0.1;
    A[locr*(j-1)+(i-1)]=(i-1)*3+j;
	fprintf(cp,"%lf ",A[locr*(j-1)+(i-1)]);
  }
  fprintf(cp,"\n");
}

} else {
descset_(descA,&n,&n,&nb,&nb,&IZERO,&IZERO,&INONE,&IONE);
}
fclose(cp);

/* Initialize the solver and set parameters */
StrumpackDensePackage<myscalar,myreal> sdp(MPI_COMM_WORLD);
sdp.use_HSS=true;
sdp.levels_HSS=4;
sdp.min_rand_HSS=10;
sdp.lim_rand_HSS=5;
sdp.inc_rand_HSS=10;
sdp.max_rand_HSS=100;
sdp.tol_HSS=1e-12;
sdp.steps_IR=10;
sdp.tol_IR=1e-10;

/* Compression */
sdp.compress(A,descA);

/* Accuracy checking */
sdp.check_compression(A,descA);

/* Factorization */
sdp.factor(A,descA);

/* Set the RHS (random vector) and the solution space /
if(myid<nprownpcol) {
locr=numroc_(&n,&nb,&myrow,&IZERO,&nprow);
locc=numroc_(&nrhs,&nb,&mycol,&IZERO,&npcol);
B=new myscalar[locrlocc];
char file[256];
sprintf(file,"B_%d.txt",myid);
cp=fopen(file,"w+");
dummy=std::max(1,locr);
descinit_(descXB,&n,&nrhs,&nb,&nb,&IZERO,&IZERO,&ctxt,&dummy,&ierr);
for(i=0;i<locrlocc;i++) {
// B[i]=static_cast(rand())/(static_cast(RAND_MAX));
B[i]=1;
fprintf(cp,"%lf\n",B[i]);
}
X=new myscalar[locr*locc];
} else {
descset_(descXB,&n,&nrhs,&nb,&nb,&IZERO,&IZERO,&INONE,&IONE);
}

fclose(cp);
char file[256];
sprintf(file,"X_%d.txt",myid);
cp=fopen(file,"w+");
/* Triangular solution */
sdp.solve(X,descXB,B,descXB);

for(i=0;i<locrlocc;i++) {
fprintf(cp,"%lf\n",X[i]);
}
fclose(cp);
/ Accuracy checking */
sdp.check_solution(A,descA,X,descXB,B,descXB);

/* Iterative refinement */
sdp.refine(A,descA,X,descXB,B,descXB);

/* Statistics */
sdp.print_statistics();

/* Clean-up */
delete[] A;
delete[] B;
delete[] X;

/* The end */
MPI_Finalize();
return 0;

}
`,the right result is [-2.5000
4.0000
-1.5000],but I obtain [-1.0000
1.00000
0.00000].Thanks you.

from strumpack.

I think something is wrong with the matrix you generate. The matrix always has rank 2. Do you agree?
So, I'm not sure how you got to the 'right result'.
If I run this, I get different results depending on how many MPI processes are used. Ideally the code should return some kind of a warning when encountering a singular matrix. I'm not sure why it doen't do that in this case. I will investigate.

from strumpack.

OK,THANKS.

from strumpack.