MKL Scalapack PDGESV

[mfactor]

Hello all,

I'm trying to do the block cyclic matrix distribution using MKL and Scalapack's PDGESV. The matrix is transposed, as in Lapack, but the results are not correct. In the PDGESV call, the index IA and JA is the starting index in the local matrix, and so always one, or should be the index that this specific block has in the total matrix?

Thank you!

#include
#include
#include
#include
#include
#include
#include
#include
#define mat(matriz,coluna,i,j) (matriz[i*coluna+j])


#define p_of_i(i,bs,p) ( int((i-1)/bs)%p)
#define l_of_i(i,bs,p) (int((i-1)/(p*bs)))
#define x_of_i(i,bs,p) (((i-1)%bs)+1)


using namespace std;


extern "C" {
/* BLACS C interface */
void Cblacs_get( int context, int request, int* value);
int Cblacs_gridinit( int* context, char * order, int np_row, int np_col);
void Cblacs_gridinfo( int context, int* np_row, int* np_col, int* my_row, int* my_col);
int numroc_( int *n, int *nb, int *iproc, int *isrcproc, int *nprocs);
void Cblacs_gridexit(int ictxt);
void pdgesv_(MKL_INT *n, MKL_INT *nrhs, double *a, MKL_INT *ia,
MKL_INT *ja, MKL_INT *desca, MKL_INT *ipiv, double *b, MKL_INT *ib, MKL_INT *jb, MKL_INT *descb, MKL_INT *info);
void dgesv_( MKL_INT* n, MKL_INT* nrhs, double* a, MKL_INT* lda, MKL_INT* ipiv, double* b, MKL_INT* ldb, MKL_INT* info );

void descinit_( int *desc, int *m, int *n, int *mb, int *nb, int *irsrc, int *icsrc,
int *ictxt, int *lld, int *info);
double pdlamch_( int *ictxt , char *cmach);
double pdlange_( char *norm, int *m, int *n, double *A, int *ia, int *ja, int *desca, double *work);

void pdlacpy_( char *uplo, int *m, int *n, double *a, int *ia, int *ja, int *desca,
double *b, int *ib, int *jb, int *descb);

void pdgemm_( char *TRANSA, char *TRANSB, int * M, int * N, int * K, double * ALPHA,
double * A, int * IA, int * JA, int * DESCA, double * B, int * IB, int * JB, int * DESCB,
double * BETA, double * C, int * IC, int * JC, int * DESCC );
int indxg2p_( int *indxglob, int *nb, int *iproc, int *isrcproc, int *nprocs);

//void descinit_(int *descA, int *n2, int *n , int *nb, int *nb2, int *izero2, int *izero, int * ictxt, int * itemp, int *info);
}

void find_nps(int np, int &nprow, int & npcol);
int getIndex(int row, int col,int NCOLS) {return row*NCOLS+col;}




int main(int argv, char ** argc) {

MPI::Init();
int rank = MPI::COMM_WORLD.Get_rank();
int nprocs = MPI::COMM_WORLD.Get_size();

// std::cout<<"Returned: "<<" ";
// std::cout << "Hello World! I am " << rank << " of " << nprocs <<
// std::endl;
//srand(1);
int myrow=0;
int mycol=0;
int ictxt;
int nprow=0,npcol=0;
int locR=0, locC=0;
int block = 2;
int myzero = 0;
int matrix_size = 9;
int myone = 1;
int nrhs = 1;
int info=0;
int i=0,j=0;
double mone=(-1.0e0),pone=(1.0e0);
double AnormF, XnormF, RnormF, BnormF, residF,eps;

find_nps(nprocs,nprow,npcol);

if(rank == 0) {
cout <<"nprocs: "<<<" nprow "<<<" npcol "<<<<"rank "<<< " nprow " <<<" "<<" npcol "<<< " myrow "<<< " mycol "<<< matrix_size; j++) {
// A[getIndex(i,j,matrix_size)] = 0;//(double)(rand()%100)/100.0;

Acpy[getIndex(i,j,matrix_size)] = A[getIndex(i,j,matrix_size)];
}

//B = 0;//(double)(rand()%100)/100.0;
Bcpy = B;

}


// cout<< "Global Matrix" << matrix_size; i++) {
for(j=0; j< matrix_size; j++) {
A[getIndex(i,j,matrix_size)] = Acpy[getIndex(j,i,matrix_size)];

//cout<<" "<<<<"++++++++++"<< matrix_size; i++) {
for(j=0; j< matrix_size; j++) {
Acpy[getIndex(i,j,matrix_size)] = A[getIndex(i,j,matrix_size)];
}
}

// if(rank ==0) {
// cout <<"Global Vector" << matrix_size; i++) {
// cout<<"* "<<<<"++++++++"<<<"rank "<<< " locR " <<<" "<<" locC "<<< matrix_size; i++) {
cout <<"dgesv "<<< matrix_size+1; i++) {
for(j=1; j< matrix_size+1; j++) {

int pi = p_of_i(i,block,nprow);
int li = l_of_i(i,block,nprow);
int xi = x_of_i(i,block,nprow);

int pj = p_of_i(j,block,npcol);
int lj = l_of_i(j,block,npcol);
int xj = x_of_i(j,block,npcol);

if( (pi == myrow) && (pj == mycol)) {
il = li*block+xi;
jl = lj*block+xj;
local_matrix[getIndex(il-1, jl-1, locC)] = A[getIndex(i-1,j-1,matrix_size)];
}

if( (pi == myrow) &&(mycol==0) ){
il = li*block+xi;
local_know_vector[il-1] = B[i-1];
}
}

}

//if(rank ==0) {
cout <<"Local Matrix myrow " <<<" mycol "<<< locR; i++) {
for(j=0; j< locC; j++) {
cout<<" "<<<<"++++++++"<<<"Local Vector rank " <<< locR; i++) {
cout<<"* "<<<<"++++++++"<< locR; i++)
{
cout<<"**"<<"rank "<<<" "<<<<"++++++++"<< min_nprow) min_nprow = nprow;
}

if(nprow ==1 ) break;

}

nprow = min_nprow;
npcol = min_npcol;

}

[Artem_V_Intel]

Hello,

IA and JA should represent specific block in the global matrix A.

Thanks,
Art

[mfactor]

Quoting - Artem Vorobiev (Intel)

Hello,

IA and JA should represent specific block in the global matrix A.

Thanks,
Art

Not quite. I managed to do an solver using local distributed matrices with both IA and JA equal to one, independent of the block size or processor row and column. But thank you for your help.

[Artem_V_Intel]

Quoting - mfactor

Quoting - Artem Vorobiev (Intel)

Hello,

IA and JA should represent specific block in the global matrix A.

Thanks,
Art

Not quite. I managed to do an solver using local distributed matrices with both IA and JA equal to one, independent of the block size or processor row and column. But thank you for your help.

Hello,

It is very interesting, because according to MKL Manual p2202:

ia, ja (global) INTEGER. The row and column indices in the global
array A indicating the first row and the first column of
sub(A), respectively.

Could you please submit an example of code? (Looks like first example is broken:

[cpp]cout<<" "<<<<"++++++++"<<<"Local Vector rank " <<< locR; i++) {[/cpp]

) Could you also specify what MKL version do you use?

Thanks,
Art