scalapack pdgesv c++

Open discussion regarding features, bugs, issues, vendors, etc.

scalapack pdgesv c++

Postby savasagduk » Thu Mar 29, 2007 5:23 am

hi !

I'm new in c++ i 'm trying to learn parallel programming but most of the example codes in the web is written in fortran ..

I want to use scalapack routine pdgesv and i find an example code written in c++..

but that program not include matrix distribution ,just it generates partial matrices and operates on them .. sorry i'm not well in english .. yes i found an example program that generates a matrix and a vector but is written in fortran..i have tried to convert that code to c++ but i confuse..

can you help me for adding "a matrix and vector distribution " to the A and B respectively , in the following program?

// STL
#include <iostream>

// MPI (must be before #define REAL)
#include "mpi.h"

using std::cout;
using std::endl;

extern "C" // {{{
{
void Cblacs_pinfo (int* mypnum, int* nprocs);
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);
void Cblacs_gridexit (int context);
void Cblacs_exit (int error_code);
int numroc_ (int *n, int *nb, int *iproc, int *isrcproc, int *nprocs);
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 pdgesv_ (int *n, int *nrhs, double *A, int *ia, int *ja, int *desca, int* ipiv, double *B, int *ib, int *jb, int *descb, int *info);
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);
} // }}}

int main(int argc, char **argv)
{
// Initialize
MPI::Init(argc, argv);
int rank = MPI::COMM_WORLD.Get_rank();
int size = MPI::COMM_WORLD.Get_size();

// Info
cout << "Process #" << rank << " of " << size << " processes just started ...\n";
cout << "" << endl;

// Data
int n = 500; // nrows=ncols of a squared matrix
int nrhs = 1; // number of Right-Hand-Side columns
int nprow = 2; // grid size
int npcol = 2; // grid size
int nb = 64; // number of blocks

// Check
if (nb>n) nb=n;
if (nprow*npcol>size)
{
if (rank==0)
{
cout << "ERROR: There is not enough processes (" << size << "<" << nprow*npcol << ") to make a p-by-q process grid" << endl;
MPI::Finalize();
return 1;
}
}

// Initialize BLACS grid
int iam, nprocs, ictxt;
int myrow, mycol;
Cblacs_pinfo (&iam, &nprocs);
Cblacs_get (-1, 0, &ictxt);
Cblacs_gridinit (&ictxt, /*order*/"Col", nprow, npcol);
Cblacs_gridinfo (ictxt, &nprow, &npcol, &myrow, &mycol);

// Work only the process in the process grid
if ((myrow>-1)&&(mycol>-1)&&(myrow<nprow)&&(mycol<npcol))
{
// Compute the size of the local matrices
int izero = 0;
int np = numroc_(&n , &nb, &myrow, &izero, &nprow);
int nq = numroc_(&n , &nb, &mycol, &izero, &npcol);
int nqrhs = numroc_(&nrhs, &nb, &mycol, &izero, &npcol);

// Allocate and fill the matrices A and B
double * A = new double [np*nq];
double * Acpy = new double [np*nq];
double * B = new double [np*nqrhs];
double * X = new double [np*nqrhs];
double * R = new double [np*nqrhs];
int * ippiv = new int [np+nb];

// Generate random A matrix
int k=0;
for (int i=0; i<np; i++)
for (int j=0; j<nq; j++)
{
A[k] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
k++;
}

// Generate random B vector
k=0;
for (int i=0; i<np; i++)
for (int j=0; j<nqrhs; j++)
{
B[k] = ((double) rand()) / ((double) RAND_MAX) - 0.5;
k++;
}

// Initialize the array descriptor for the matrix A and B
int itemp = (np<1 ? 1 : np);
int descA[9];
int descB[9];
int info;
descinit_(descA, &n, &n , &nb, &nb, &izero, &izero, &ictxt, &itemp, &info);
descinit_(descB, &n, &nrhs, &nb, &nb, &izero, &izero, &ictxt, &itemp, &info);

// Make a copy of A and the rhs for checking purposes
int ione = 1;
pdlacpy_("All", &n, &n , A, &ione, &ione, descA, Acpy, &ione, &ione, descA);
pdlacpy_("All", &n, &nrhs, B, &ione, &ione, descB, X , &ione, &ione, descB);

// Call ScaLAPACK PDGESV routine
double MPIt1 = MPI::Wtime();
pdgesv_(&n, &nrhs, A, &ione, &ione, descA, ippiv, X, &ione, &ione, descB, &info);
double MPIt2 = MPI::Wtime();
double MPIelapsed = MPIt2-MPIt1;

// Froebenius norm = ||A * X - B|| / ( ||X|| * ||A|| * eps * N )
double mone=-1.0;
double pone= 1.0;
pdlacpy_("All", &n, &nrhs, B, &ione, &ione, descB, R, &ione, &ione, descB);
pdgemm_ ("N", "N", &n, &nrhs, &n, &pone, Acpy, &ione, &ione, descA, X, &ione, &ione, descB, &mone, R, &ione, &ione, descB);
double work = 0.0;
double eps = pdlamch_(&ictxt, "Epsilon");
double AnormF = pdlange_("F", &n, &n , A, &ione, &ione, descA, &work);
double XnormF = pdlange_("F", &n, &nrhs, X, &ione, &ione, descB, &work);
double RnormF = pdlange_("F", &n, &nrhs, R, &ione, &ione, descB, &work);
double residF = RnormF/(AnormF*XnormF*eps*static_cast<double>(n));

// Information
if (iam==0)
{
cout << endl;
cout << "***********************************************" << endl;
cout << " Example of ScaLAPACK routine call: (PDGESV) " << endl;
cout << "***********************************************" << endl;
cout << endl;
cout << " n = " << n << endl;
cout << " nrhs = " << nrhs << endl;
cout << " grid = " << nprow << ", " << npcol << endl;
cout << " blocks = " << nb << "x" << nb << endl;
cout << " np = " << np << endl;
cout << " nq = " << nq << endl;
cout << " nqrhs = " << nqrhs << endl;
cout << endl;
cout << "MPI elapsed time during (pdgesv) = " << MPIelapsed << endl;
cout << "Info returned by (pdgesv) = " << info << endl;
cout << endl;
cout << "Froebenius norm (residual) = " << residF << endl;
cout << endl;
}

// Clean up
delete [] A;
delete [] Acpy;
delete [] B;
delete [] X;
delete [] R;
delete [] ippiv;

// Grid exit
Cblacs_gridexit(0);
}

// Information
if(iam==0)
{
cout << "***********************************************\n";
cout << " END \n";
cout << "***********************************************\n";
}

// Finalize
MPI::Finalize();

return 0;
}




[code][/code]
savasagduk
 
Posts: 2
Joined: Thu Mar 29, 2007 4:22 am

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 4 guests

cron