Skip to Content.
Sympa Menu

charm - [charm] Cholesky decomposition -- please help and advice

charm AT lists.cs.illinois.edu

Subject: Charm++ parallel programming system

List archive

[charm] Cholesky decomposition -- please help and advice


Chronological Thread 
  • From: Rabab Al-omairy <rabab.omairy AT kaust.edu.sa>
  • To: charm AT lists.cs.illinois.edu
  • Subject: [charm] Cholesky decomposition -- please help and advice
  • Date: Sat, 29 Apr 2017 11:39:46 +0300


Dear Sir/ Madam,

I am trying to implement cholesky factorization using charm++. I have problem in which program is blocked and I haven't received any respond from program.

I initialise matrix and decompose them into tiles and distributed them to chars. However, when I am trying to invoke cholesky (run function) the program is blocked.

Could you please help me and guide me on solving this issue.? 

In the attachment, you can find the codes  that I have developed

--

Best Regards,

Rabab Alomairy


King Abdullah University of Science and Technology

Kingdom of Saudi Arabia

Mobile KSA: +966 (0) 548884792

E-mail: rabab.alomairy AT kaust.edu.sa



 

 



This message and its contents, including attachments are intended solely for the original recipient. If you are not the intended recipient or have received this message in error, please notify me immediately and delete this message from your computer system. Any unauthorized use or distribution is prohibited. Please consider the environment before printing this email.

PNG image

#include "cholesky2.decl.h"
#include "rand48_replacement.h"
#include <cmath>
#include <limits>
#include <map>
#include <string.h>
#include <mkl.h>

CProxy_Main mainProxy;

//void dgemm_ (const char *transa, const char *transb, int *l, int *n, int *m, double *alpha,
//             const void *a, int *lda, void *b, int *ldb, double *beta, void *c, int *ldc);
//void dtrsm_ (char *side, char *uplo, char *transa, char *diag, int *m, int *n, double *alpha,
//             double *a, int *lda, double *b, int *ldb);
//void dtrmm_ (char *side, char *uplo, char *transa, char *diag, int *m, int *n, double *alpha,
//             double *a, int *lda, double *b, int *ldb);
//void dsyrk_ (char *uplo, char *trans, int *n, int *k, double *alpha, double *a, int *lda,
//             double *beta, double *c, int *ldc);

void charm_potrf(double * const A, int ts, int ld)
{
   static int INFO;
   static const char L = 'L';
   dpotrf(&L, &ts, A, &ld, &INFO);
}

void charm_trsm(double *A, double *B, int ts, int ld)
{
   static char LO = 'L', TR = 'T', NU = 'N', RI = 'R';
   static double DONE = 1.0;
   dtrsm(&RI, &LO, &TR, &NU, &ts, &ts, &DONE, A, &ld, B, &ld );
}

void charm_syrk(double *A, double *B, int ts, int ld)
{
   static char LO = 'L', NT = 'N';
   static double DONE = 1.0, DMONE = -1.0;
   dsyrk(&LO, &NT, &ts, &ts, &DMONE, A, &ld, &DONE, B, &ld );
}

void charm_gemm(double *A, double *B, double *C, int ts, int ld)
{
   static const char TR = 'T', NT = 'N';
   static double DONE = 1.0, DMONE = -1.0;
   dgemm(&NT, &TR, &ts, &ts, &ts, &DMONE, A, &ld, B, &ld, &DONE, C, &ld);
}

class Main : public CBase_Main {
  double startTime;
    //
    unsigned int blockSize, numBlocks, matrix_size;
    double * matrixA; double * original_matrix; double **sub_matrix;
CProxy_Block a;
public:
  Main(CkArgMsg* m) {
    if (m->argc > 2) {
        //
      matrix_size= atoi(m->argv[1]);
      blockSize = atoi(m->argv[2]);
    } else {
      CkAbort("Usage: cholesky2 matrixSize blockSize ");
    }
numBlocks=matrix_size/blockSize;
    mainProxy = thisProxy;
//      a = CProxy_Block::ckNew(blockSize, numBlocks, numBlocks, numBlocks); //for computing
  
    startTime = CkWallTimer();
	Matrix_ini();
//       a = CProxy_Block::ckNew(blockSize, numBlocks, numBlocks, numBlocks);
//print();  
a = CProxy_Block::ckNew(blockSize, numBlocks, numBlocks, numBlocks);     
decompose_matrix();
//a.run();
ckout<<"I am in the main 1\n";
//a.run(CkCallback(CkReductionTarget(Main, done), thisProxy));
//ckout<<"Fun run";
}

void Matrix_ini()
    {
       int nn=matrix_size*matrix_size;
	const int n=matrix_size;
        matrixA=new double[nn];
	original_matrix= new double [nn];
        int ISEED[4] = {0,0,0,1};
        int intONE=1;

        for (int i = 0; i < nn; i+=matrix_size) {
            dlarnv_(&intONE, &ISEED[0], &n, &matrixA[i]);
        }
    
        for (int i=0; i<matrix_size; i++) {
            for (int j=0; j<matrix_size; j++) {
                matrixA[j*matrix_size + i] = matrixA[j*matrix_size + i] + matrixA[i*matrix_size + j];
                matrixA[i*matrix_size + j] = matrixA[j*matrix_size + i];
            }
        }

        double alpha= (double) matrix_size;
        for (int i = 0; i < matrix_size; i++)
            matrixA[ i + i * matrix_size ] += alpha;
 	
	for (int i = 0; i < n * n; i++ ) {
      	original_matrix[i] = matrixA[i];
   	}
    }
    
    void print()
    {
        int i, j, m, n, k, r;
        for (i=0;i<matrix_size;i++){
            for ( j=0; j<matrix_size; j++){
               ckout<<matrixA[i*matrix_size + j] <<" ";
            }
            ckout<< endl;
        }
	ckout<< endl;
       for (m=0;m<matrix_size;m+=blockSize){
	for (n=0;n<matrix_size;n+=blockSize){
	 for (i=m, k=0; k<blockSize;i++, k++){
	  for (j=n, r=0; r<blockSize;j++, r++)
		{
             ckout<< matrixA[i*matrix_size+j]<<" ";
		} //j, r
		ckout<< endl;
        }//i, k
		 ckout<< endl<<endl;
	} //n
         } //m
    }
void decompose_matrix()
{
int total_number_blocks=(matrix_size/blockSize) * (matrix_size/blockSize);
int dimx=0, dimy=0; int i, j, n, k, r, m;
int total_number_wn_block=blockSize*blockSize;

sub_matrix= new double* [total_number_blocks];
for(int i = 0; i < matrix_size; ++i)
    sub_matrix[i] = new double[total_number_wn_block];

 for (m=0;m<matrix_size;m+=blockSize){
        for (n=0;n<matrix_size;n+=blockSize){
         for (i=m, k=0; k<blockSize;i++, k++){
          for (j=n, r=0; r<blockSize;j++, r++)
                {
        //     ckout<< matrixA[i*matrix_size+j]<<" ";
		sub_matrix[dimx][dimy]=matrixA[i*matrix_size+j];
		//ckout<<sub_matrix[dimx][dimy]<<" ";
		dimy++;
                } //j, r
               //ckout<< endl;
        }//i, k
	 dimx++;
    	dimy=0;
        } //n
         } //m
int blocks=0;
int x, y;
for (y=0;y<numBlocks;y++){
   for (x=0;x<numBlocks;x++){
a(x,y).initlize_block(sub_matrix[blocks], total_number_wn_block, blocks, mainProxy);
blocks++;
} 
}
}
  void done() {
    double endTime = CkWallTimer();
    CkPrintf("Matrix multiply of %u blocks with %u elements each (%u^2) finished in %f seconds\n", 
             numBlocks, blockSize, numBlocks*blockSize, endTime - startTime);
   CkExit();
  }
//It is used just for syn
void done_decomp(int value){
CkPrintf("The the total value is: %d \n", value); 
a.run();

}
   
};

class Block : public CBase_Block {
  unsigned int blockSize, numBlocks, block, matrix_size, elems_within_block;
  double *matrix;
  int k, n, m;
  Block_SDAG_CODE
  public:
  Block(unsigned int blockSize_, unsigned int numBlocks_)
    : blockSize(blockSize_), numBlocks(numBlocks_)
  {
    elems_within_block = blockSize * blockSize;
    matrix = new double[elems_within_block];
    k=0;
  }
   void initlize_block(double *matrix2, int blockSize, int total,  CProxy_Main mainProxy)  
{
    for (int i=0; i<elems_within_block ;i++){
   matrix[i]=matrix2[i];
//ckout<<matrix[i] <<" ";
}
//ckout<<endl;


//CkCallback cb(CkReductionTarget(Main, done), mainProxy);
//contribute(cb);
CkCallback cb(CkReductionTarget(Main, done_decomp), mainProxy);
contribute(sizeof (int), &total, CkReduction::sum_int, cb);
} 
void print_blocks()
{

for (int i=0;i<elems_within_block;i++){
        CkPrintf(" %d ", matrix[i]);
}
}
void send_dpotrf(int iter){
    if(thisIndex.x==iter && thisIndex.y==iter)
      {
	for (int i=iter+1;i<numBlocks;i++)
	thisProxy(iter, i).dtrsm_receive(iter, matrix, elems_within_block);
      }

	}//send_dpotrf
void send_dtrsm(int iter_k, int iter_n){

 if (thisIndex.x==iter_k && thisIndex.y==iter_n)
   {
	   for (int j=iter_k+1; j<=iter_n;j++)
		{
		if (j==iter_n) thisProxy(iter_n, iter_n).dsyrk_receive(iter_k, matrix, elems_within_block);
		else thisProxy(j, iter_n).dgemm_receive_horizontal(iter_k, matrix, elems_within_block);
		}//for
          for (int i=iter_n+1; i<numBlocks;i++)
		{
		thisProxy(iter_n, i).dgemm_receive_vertical(iter_k, matrix, elems_within_block);
		}
	}
} //send_dtrsm

void send_dsyrk(int n, int k){
     if ((thisIndex.x==n) && (thisIndex.y==k))
	thisProxy(n, k).dpotrf_receive(n, matrix, elems_within_block);
}
  Block(CkMigrateMessage*) {}

};

#include "cholesky2.def.h"

Attachment: cholesky2.ci
Description: Binary data

Attachment: Makefile
Description: Binary data



  • [charm] Cholesky decomposition -- please help and advice, Rabab Al-omairy, 04/29/2017

Archive powered by MHonArc 2.6.19.

Top of Page