Using the Dense matrix computation package wrong.

[jerome92]

Hello, I used the Dense matrix computation package to solve a 3*3 matrix,A=[1,2,3;4,5,6;7,8,9],b=[1,1,1];and the result is wrong. Please tell me how can i do? Thanks.

[pghysels]

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

[jerome92]

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; blacsget(&IZERO,&IZERO,&ctxt); blacsgridinit(&ctxt,"R",&nprow,&npcol); blacsgridinfo(&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[locrlocc]; 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[locrlocc]; } 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.

[pghysels]

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.

[jerome92]

OK,THANKS.