Index: /issm/trunk-jpl/src/c/classes/objects/Bucket.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Bucket.h (revision 14749) +++ /issm/trunk-jpl/src/c/classes/objects/Bucket.h (revision 14750) @@ -10,13 +10,17 @@ #include "./Object.h" #include "../../shared/Alloc/alloc.h" +#include "../../include/macros.h" #include "../../Container/DataSet.h" #include "../../toolkits/toolkitsenums.h" /*}}}*/ -#define BUCKETSIZEOFREQUESTS 6 /*how many MPI_Isend requests does it take to transfer the contents of a bucket to another cpu?*/ - +/*how many MPI_Isend requests does it take to transfer the contents of a bucket to another cpu?*/ +#define MATRIXBUCKETSIZEOFREQUESTS 7 +#define VECTORBUCKETSIZEOFREQUESTS 5 +typedef enum {VECTOR_BUCKET, MATRIX_BUCKET} BucketType; template class Bucket: public Object{ private: + int type; //either a VECTOR_BUCKET or MATRIX_BUCKET int m,n; /*size of local matrix we are storing*/ /*row and column indices of the matrix we are storing*/ @@ -30,4 +34,5 @@ /*constructors, destructors: */ Bucket(){ /*{{{*/ + this->type=0; this->m=0; this->n=0; @@ -38,4 +43,5 @@ } /*}}}*/ Bucket(int min,int* idxmin,int nin,int* idxnin,doubletype* valuesin,InsMode modein){ /*{{{*/ + this->type=MATRIX_BUCKET; this->m=min; this->n=nin; @@ -55,4 +61,5 @@ } /*}}}*/ Bucket(int min,int* idxmin,doubletype* valuesin,InsMode modein){ /*{{{*/ + this->type=VECTOR_BUCKET; this->m=min; this->n=1; @@ -75,6 +82,7 @@ /*object virtual functions definitions:*/ void Echo(){ /*{{{*/ - printf("Bucket echo (cpu #: %i): \n",IssmComm::GetRank()); - printf("# rows: %i, #cols: %i\n",this->m,this->n); + _printLine_("Bucket echo (cpu #: "<m<<" num cols: "<n); } /*}}}*/ void DeepEcho(){ /*{{{*/ @@ -82,15 +90,24 @@ _printLine_("Bucket echo (cpu #: "<m<<" num cols: "<n); - for (i=0;im;i++){ - _printLine_("row "<idxm[i]<<", column indices: "); - for (j=0;jn;j++){ - _printLine_(" "<idxn[j]); - } - _printLine_("values: "); - for (j=0;jn;j++){ - _printLine_(" "<values[m*i+j]); - } - } + if(type=MATRIX_BUCKET){ + for (i=0;im;i++){ + _printLine_("row "<idxm[i]<<", column indices: "); + for (j=0;jn;j++){ + _printLine_(" "<idxn[j]); + } + _printLine_("values: "); + for (j=0;jn;j++){ + _printLine_(" "<values[m*i+j]); + } + } + } + else if(type=VECTOR_BUCKET){ + for (i=0;im;i++){ + _printLine_("row "<idxm[i]<<", value " << this->values[i]); + } + } + else _error_("unknown type of bucket!"); } /*}}}*/ @@ -115,5 +132,10 @@ if (rowranks[idxm[i]]==rank_i){ /*This row belongs to cpu rank_i, so spawn a bucket with this row, and add it to the bucketsofcpu_i dataset: */ - bucketsofcpu_i->AddObject(new Bucket(1,idxm+i,n,idxn,values+n*i,mode)); + if(type==MATRIX_BUCKET){ + bucketsofcpu_i->AddObject(new Bucket(1,idxm+i,n,idxn,values+n*i,mode)); + } + else{ + bucketsofcpu_i->AddObject(new Bucket(1,idxm+i,values+i,mode)); + } } } @@ -132,4 +154,12 @@ }; /*}}}*/ + void SetLocalVectorValues(double* local_vector,int lower_row){ /*{{{*/ + + int i; + for(i=0;iSize()*BUCKETSIZEOFREQUESTS; //this is to take into account all the MPI_ISend calls in each bucket. + num_requests+=bucketspercpu[i]->Size()*MATRIXBUCKETSIZEOFREQUESTS; //this is to take into account all the MPI_ISend calls in each bucket. num_requests++; //this is to take into account on MPI_ISend in BucketsSend. } Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14749) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14750) @@ -124,7 +124,122 @@ /*}}}*/ /*FUNCTION Assemble{{{*/ - void Assemble(void){ - printf("not imlemented yet!"); - exit(1); + void Assemble(){ + + + int i; + int j; + int k; + int my_rank; + int num_procs; + int *RowRank = NULL; + + DataSet **bucketspercpu = NULL; + int *bucketspercpu_sizes = NULL; + MPI_Request *requests = NULL; + MPI_Status *statuses = NULL; + MPI_Status status; + int num_requests = 0; + DataSet *mybuckets = NULL; + int lower_row; + int upper_row; + int count = 0; + + int size; + + + + /*some communicator info: */ + num_procs=IssmComm::GetSize(); + my_rank=IssmComm::GetRank(); + MPI_Comm comm=IssmComm::GetComm(); + + /*First, make a vector of size M, which for each row between 0 and M-1, tells which cpu this row belongs to: */ + RowRank=DetermineRowRankFromLocalSize(M,m,comm); + + /*Now, sort out our dataset of buckets according to cpu ownership of rows: */ + bucketspercpu=xNew(num_procs); + bucketspercpu_sizes=xNew(num_procs); + + for(i=0;iSize();j++){ + Bucket* bucket=(Bucket*)buckets->GetObjectByOffset(j); + bucket->SpawnBucketsPerCpu(bucketsofcpu_i,i,RowRank); + } + bucketspercpu[i]=bucketsofcpu_i; + bucketspercpu_sizes[i]=bucketsofcpu_i->Size(); + } + + /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this + * dataset owns correspond to rows that are owned by cpu i, not j!:*/ + + /*First, figure out how many requests are going to be sent by MPI_Isend. Do this a little bit better? */ + for(i=0;iSize()*VECTORBUCKETSIZEOFREQUESTS; //this is to take into account all the MPI_ISend calls in each bucket. + num_requests++; //this is to take into account on MPI_ISend in BucketsSend. + } + } + + /*Initialize array to track requests and statuses: */ + requests=new MPI_Request[num_requests]; + statuses=new MPI_Status[num_requests]; + + /*Now, go through all our bucketspercpu datasets, and send them to the corresponding cpus. Do not send our own buckets though!: */ + count=0; //count requests + for(i=0;iSize();k++){ + Bucket* bucket=(Bucket*)buckets->GetObjectByOffset(k); + bucket->Isend(j,requests,&count,comm); + } + } + } + } + else{ + + /*Receive buckets from cpu i, and add them to my own my_rank bucket list: */ + /*First, are we receiving anything from sender_rank? :*/ + MPI_Recv(&size,1, MPI_INT,i,1, comm, &status); + + /*If so, started receiving extra buckets and plug them into out buckets: */ + if(size){ + for(j=0;j* bucket=new Bucket(); + bucket->Recv(i,comm); + bucketspercpu[my_rank]->AddObject(bucket); + } + } + } + } + /*Wait for all requests to complete: */ + MPI_Waitall(num_requests,requests,statuses); + + /*Every cpu now has a dataset of buckets in bucketspercpu[my_rank], which holds all the values + *local to this cpu that should be added to the global matrix. Just do that: */ + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,comm); + mybuckets=bucketspercpu[my_rank]; + + for(i=0;iSize();i++){ + Bucket* bucket=(Bucket*)mybuckets->GetObjectByOffset(i); + bucket->SetLocalVectorValues(this->vector,lower_row); + } + + /*Free ressources:{{{*/ + xDelete(RowRank); + for(i=0;i(bucketspercpu); + xDelete(bucketspercpu_sizes); + xDelete(requests); + /*}}}*/ } /*}}}*/