Index: /issm/trunk-jpl/src/c/classes/objects/Bucket.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Bucket.h (revision 14833) +++ /issm/trunk-jpl/src/c/classes/objects/Bucket.h (revision 14834) @@ -187,4 +187,5 @@ int *modes_forcpu = NULL; + /*initialize buffers: */ row_indices_forcpu=*prow_indices_forcpu; @@ -216,7 +217,42 @@ }; /*}}}*/ + void Marshall(int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu){ /*{{{*/ + + /*intermediary: */ + int i; + int j; + + /*buffers: */ + int *row_indices_forcpu = NULL; + doubletype *values_forcpu = NULL; + int *modes_forcpu = NULL; + + + /*initialize buffers: */ + row_indices_forcpu=*prow_indices_forcpu; + values_forcpu=*pvalues_forcpu; + modes_forcpu=*pmodes_forcpu; + + /*fill buffers with out values and indices and modes: */ + for(i=0;i class IssmAbsVec; +class Parameters; + template class IssmAbsMat{ @@ -53,4 +56,5 @@ virtual void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode)=0; virtual void Convert(MatrixType type)=0; + virtual IssmAbsVec* Solve(IssmAbsVec* pf, Parameters* parameters)=0; }; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h (revision 14833) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h (revision 14834) @@ -44,5 +44,5 @@ /*IssmDenseMat constructors, destructors*/ - /*FUNCTION IssmDenseMat(){{{*/ + /*IssmDenseMat(){{{*/ IssmDenseMat(){ @@ -52,5 +52,5 @@ } /*}}}*/ - /*FUNCTION IssmDenseMat(int M,int N){{{*/ + /*IssmDenseMat(int M,int N){{{*/ IssmDenseMat(int pM,int pN){ @@ -61,5 +61,5 @@ } /*}}}*/ - /*FUNCTION IssmDenseMat(int M,int N, doubletype sparsity){{{*/ + /*IssmDenseMat(int M,int N, doubletype sparsity){{{*/ IssmDenseMat(int pM,int pN, doubletype sparsity){ @@ -70,5 +70,5 @@ } /*}}}*/ - /*FUNCTION IssmDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/ + /*IssmDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/ IssmDenseMat(int m,int n,int pM,int pN,int* d_nnz,int* o_nnz){ @@ -79,5 +79,5 @@ } /*}}}*/ - /*FUNCTION IssmDenseMat(doubletype* serial_mat,int M,int N,doubletype sparsity){{{*/ + /*IssmDenseMat(doubletype* serial_mat,int M,int N,doubletype sparsity){{{*/ IssmDenseMat(doubletype* serial_mat,int pM,int pN,doubletype sparsity){ @@ -92,5 +92,5 @@ } /*}}}*/ - /*FUNCTION IssmDenseMat(int M,int N, int connectivity, int numberofdofspernode){{{*/ + /*IssmDenseMat(int M,int N, int connectivity, int numberofdofspernode){{{*/ IssmDenseMat(int pM,int pN, int connectivity,int numberofdofspernode){ @@ -101,5 +101,5 @@ } /*}}}*/ - /*FUNCTION ~IssmDenseMat(){{{*/ + /*~IssmDenseMat(){{{*/ ~IssmDenseMat(){ @@ -110,6 +110,6 @@ /*}}}*/ - /*IssmDenseMat specific routines */ - /*FUNCTION Echo{{{*/ + /*IssmAbsMat virtual functions*/ + /*Echo{{{*/ void Echo(void){ @@ -124,5 +124,5 @@ } /*}}}*/ - /*FUNCTION Assemble{{{*/ + /*Assemble{{{*/ void Assemble(void){ @@ -131,5 +131,5 @@ } /*}}}*/ - /*FUNCTION Norm{{{*/ + /*Norm{{{*/ doubletype Norm(NormMode mode){ @@ -156,5 +156,5 @@ } /*}}}*/ - /*FUNCTION GetSize{{{*/ + /*GetSize{{{*/ void GetSize(int* pM,int* pN){ @@ -164,5 +164,5 @@ } /*}}}*/ - /*FUNCTION GetLocalSize{{{*/ + /*GetLocalSize{{{*/ void GetLocalSize(int* pM,int* pN){ @@ -172,5 +172,5 @@ } /*}}}*/ - /*FUNCTION MatMult{{{*/ + /*MatMult{{{*/ void MatMult(IssmAbsVec* Xin,IssmAbsVec* AXin){ @@ -203,5 +203,5 @@ } /*}}}*/ - /*FUNCTION Duplicate{{{*/ + /*Duplicate{{{*/ IssmDenseMat* Duplicate(void){ @@ -212,5 +212,5 @@ } /*}}}*/ - /*FUNCTION ToSerial{{{*/ + /*ToSerial{{{*/ doubletype* ToSerial(void){ @@ -225,5 +225,5 @@ } /*}}}*/ - /*FUNCTION SetValues{{{*/ + /*SetValues{{{*/ void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode){ @@ -243,5 +243,5 @@ } /*}}}*/ - /*FUNCTION Convert{{{*/ + /*Convert{{{*/ void Convert(MatrixType type){ @@ -250,4 +250,43 @@ } /*}}}*/ + /*Solve{{{*/ + IssmAbsVec* Solve(IssmAbsVec* pfin, Parameters* parameters){ + + /*First off, we assume that the type of IssmAbsVec is IssmSeqVec. So downcast: */ + IssmSeqVec* pf=(IssmSeqVec*)pfin; + + #ifdef _HAVE_GSL_ + /*Intermediary: */ + int M,N,N2; + IssmSeqVec *uf = NULL; + + Kff->GetSize(&M,&N); + pf->GetSize(&N2); + + if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !"); + if(M!=N)_error_("Stiffness matrix should be square!"); + #ifdef _HAVE_ADOLC_ + ensureContiguousLocations(N); + #endif + IssmDouble *x = xNew(N); + #ifdef _HAVE_ADOLC_ + SolverxSeq(x,Kff->matrix,pf->vector,N,parameters); + #els + SolverxSeq(x,Kff->matrix,pf->vector,N); + #endif + + uf=new IssmSeqVec(x,N); + xDelete(x); + + /*Assign output pointers:*/ + IssmVec* out=new IssmVec(); + out->vector=(IssmAbsVec*)uf; + *pout=out; + + #else + _error_("GSL support not compiled in!"); + #endif + + }/*}}}*/ }; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h (revision 14833) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h (revision 14834) @@ -34,4 +34,5 @@ template class IssmDenseMat; template class IssmMpiDenseMat; +class Parameters; template @@ -199,4 +200,13 @@ matrix->convert(type); }/*}}}*/ + IssmVec* Solve(IssmVec* pf, Parameters* parameters){ /*{{{*/ + + IssmVec* outvector=NULL; + + outvector=new IssmVec(); + + outvector->vector=this->matrix->Solve(pf->vector,parameters); + + }/*}}}*/ }; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 14833) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 14834) @@ -142,4 +142,5 @@ printf("%g ",this->matrix[j*this->N+k]); } + printf("\n"); } } @@ -302,5 +303,5 @@ RowRank=DetermineRowRankFromLocalSize(M,m,comm); - /*Now, sort out our dataset of buckets according to cpu ownership of rows: */ + /*Now, sort out our dataset of buckets according to cpu ownership of rows: {{{*/ bucketsforcpu=xNew(num_procs); @@ -313,23 +314,41 @@ bucketsforcpu[i]=bucketsofcpu_i; } - - /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this + /*}}}*/ + + /*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!. Out of all the buckets we own, make row,col,value,insert_mode * vectors that will be shipped around the cluster: */ this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&col_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs); - - /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need + /*}}}*/ + + /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need {{{ *some scatter calls: */ numvalues_fromcpu = xNew(num_procs); for(i=0;i(numvalues_fromcpu[i]); - col_indices_fromcpu[i]=xNew(numvalues_fromcpu[i]); - values_fromcpu[i]=xNew(numvalues_fromcpu[i]); - modes_fromcpu[i]=xNew(numvalues_fromcpu[i]); - } - + MPI_Scatter(numvalues_forcpu,1,MPI_INT,numvalues_fromcpu+i,1,MPI_INT,i,comm); + } + + row_indices_fromcpu=xNew(num_procs); + col_indices_fromcpu=xNew(num_procs); + values_fromcpu=xNew(num_procs); + modes_fromcpu=xNew(num_procs); + for(i=0;i(size); + col_indices_fromcpu[i]=xNew(size); + values_fromcpu[i]=xNew(size); + modes_fromcpu[i]=xNew(size); + } + else{ + row_indices_fromcpu[i]=NULL; + col_indices_fromcpu[i]=NULL; + values_fromcpu[i]=NULL; + modes_fromcpu[i]=NULL; + } + } + /*}}}*/ + + /*Scatter values around: {{{*/ /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the MPI_Scatterv prototype: @@ -344,5 +363,4 @@ } - /*Start the scattering: */ for(i=0;i(RowRank); @@ -523,4 +543,6 @@ /*figure out size of buffers per cpu: */ + + numvalues_forcpu=xNew(num_procs); for(i=0;i* Solve(IssmAbsVec* pfin, Parameters* parameters){ + + printf("IssmMpiDenseMat solve not implemented yet!"); + exit(1); + + }/*}}}*/ }; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14833) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 14834) @@ -126,4 +126,154 @@ /*FUNCTION Assemble{{{*/ void Assemble(){ + + + int i,j; + + int *RowRank = NULL; + int num_procs; + + int *row_indices_forcpu = NULL; + int *col_indices_forcpu = NULL; + int *modes_forcpu = NULL; + doubletype *values_forcpu = NULL; + int *numvalues_forcpu = NULL; + DataSet **bucketsforcpu = NULL; + + int **row_indices_fromcpu = NULL; + int **col_indices_fromcpu = NULL; + int **modes_fromcpu = NULL; + doubletype **values_fromcpu = NULL; + int *numvalues_fromcpu = NULL; + + int lower_row; + int upper_row; + int* sendcnts = NULL; + int* displs = NULL; + int count = 0; + + /*some communicator info: */ + num_procs=IssmComm::GetSize(); + 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: {{{*/ + bucketsforcpu=xNew(num_procs); + + for(i=0;iSize();j++){ + Bucket* bucket=(Bucket*)buckets->GetObjectByOffset(j); + bucket->SpawnBucketsPerCpu(bucketsofcpu_i,i,RowRank); + } + bucketsforcpu[i]=bucketsofcpu_i; + } + /*}}}*/ + + /*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!. Out of all the buckets we own, make row,col,value,insert_mode + * vectors that will be shipped around the cluster: */ + this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs); + /*}}}*/ + + /*Now, we need to allocate on each cpu arrays to receive data from all the other cpus. To know what we need to allocate, we need {{{ + *some scatter calls: */ + numvalues_fromcpu = xNew(num_procs); + for(i=0;i(num_procs); + values_fromcpu=xNew(num_procs); + modes_fromcpu=xNew(num_procs); + for(i=0;i(size); + values_fromcpu[i]=xNew(size); + modes_fromcpu[i]=xNew(size); + } + else{ + row_indices_fromcpu[i]=NULL; + values_fromcpu[i]=NULL; + modes_fromcpu[i]=NULL; + } + } + /*}}}*/ + + /*Scatter values around: {{{*/ + /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given + * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the MPI_Scatterv prototype: + * int MPI_Scatterv( void *sendbuf, int *sendcnts, int *displs, MPI_Datatype sendtype, void *recvbuf, int recvcnt, MPI_Datatype recvtype, int root, MPI_Comm comm) :*/ + sendcnts=xNew(num_procs); + displs=xNew(num_procs); + count=0; + for(i=0;i(RowRank); + xDelete(row_indices_forcpu); + xDelete(modes_forcpu); + xDelete(values_forcpu); + xDelete(numvalues_forcpu); + + for(i=0;i(bucketsforcpu); + + for(i=0;i(rows); + xDelete(modes); + xDelete(values); + } + xDelete(row_indices_fromcpu); + xDelete(modes_fromcpu); + xDelete(values_fromcpu); + xDelete(numvalues_fromcpu); + + xDelete(sendcnts); + xDelete(displs); + /*}}}*/ + + + } + /*}}}*/ + /*FUNCTION Assemble2{{{*/ + void Assemble2(){ int i; @@ -461,4 +611,77 @@ } /*}}}*/ + /*FUNCTION BucketsBuildScatterBuffers{{{*/ + void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){ + + + /*intermediary: */ + int i,j; + int count = 0; + int total_size = 0; + int *temp_row_indices_forcpu = NULL; + doubletype *temp_values_forcpu = NULL; + int *temp_modes_forcpu = NULL; + + /*output: */ + int *numvalues_forcpu = NULL; + int *row_indices_forcpu = NULL; + doubletype *values_forcpu = NULL; + int *modes_forcpu = NULL; + + /*figure out size of buffers per cpu: */ + + numvalues_forcpu=xNew(num_procs); + for(i=0;iSize();j++){ + Bucket* bucket =(Bucket*)buckets->GetObjectByOffset(j); + count+=bucket->MarshallSize(); + } + + numvalues_forcpu[i]=count; + } + + /*now, figure out size of total buffers (for all cpus!): */ + count=0; + for(i=0;i(total_size); + values_forcpu = xNew(total_size); + modes_forcpu = xNew(total_size); + + /*we are going to march through the buffers, and marshall data onto them, so in order to not + *lose track of where these buffers are located in memory, we are going to work using copies + of them: */ + temp_row_indices_forcpu=row_indices_forcpu; + temp_values_forcpu=values_forcpu; + temp_modes_forcpu=modes_forcpu; + + /*Fill buffers: */ + for(i=0;iSize();j++){ + Bucket* bucket =(Bucket*)buckets->GetObjectByOffset(j); + bucket->Marshall(&temp_row_indices_forcpu,&temp_values_forcpu,&temp_modes_forcpu); //pass in the address of the buffers, so as to have the Marshall routine increment them. + } + } + + /*sanity check: */ + if (temp_row_indices_forcpu!=row_indices_forcpu+total_size)_error_("problem with marshalling of buckets"); + if (temp_values_forcpu!=values_forcpu+total_size)_error_("problem with marshalling of buckets"); + if (temp_modes_forcpu!=modes_forcpu+total_size)_error_("problem with marshalling of buckets"); + + /*output buffers: */ + *pnumvalues_forcpu = numvalues_forcpu; + *prow_indices_forcpu = row_indices_forcpu; + *pvalues_forcpu = values_forcpu; + *pmodes_forcpu = modes_forcpu; + } + /*}}}*/ }; #endif //#ifndef _ISSM_MPI_VEC_H_ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp (revision 14833) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmSolver.cpp (revision 14834) @@ -1,4 +1,4 @@ -/*!\file SolverxSeq - * \brief implementation of sequential solver using the GSL librarie +/*!\file IssmSolver + * \brief implementation of solver */ @@ -20,43 +20,14 @@ #endif -void IssmSolve(IssmVec** pout,IssmMat* Kffin, IssmVec* pfin, Parameters* parameters){/*{{{*/ - - /*First off, we assume that the type of IssmVec is IssmSeqVec and the type of IssmMat is IssmDenseMat. So downcast: */ - IssmSeqVec* pf=(IssmSeqVec*)pfin->vector; - IssmDenseMat* Kff=(IssmDenseMat*)Kffin->matrix; - -#ifdef _HAVE_GSL_ - /*Intermediary: */ - int M,N,N2; - IssmSeqVec *uf = NULL; - - Kff->GetSize(&M,&N); - pf->GetSize(&N2); - - if(N!=N2)_error_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !"); - if(M!=N)_error_("Stiffness matrix should be square!"); -#ifdef _HAVE_ADOLC_ - ensureContiguousLocations(N); -#endif - IssmDouble *x = xNew(N); -#ifdef _HAVE_ADOLC_ - SolverxSeq(x,Kff->matrix,pf->vector,N,parameters); -#else - SolverxSeq(x,Kff->matrix,pf->vector,N); -#endif - - uf=new IssmSeqVec(x,N); - xDelete(x); - - /*Assign output pointers:*/ - IssmVec* out=new IssmVec(); - out->vector=(IssmAbsVec*)uf; - *pout=out; - -#else - _error_("GSL support not compiled in!"); -#endif - -}/*}}}*/ +void IssmSolve(IssmVec** pout,IssmMat* Kff, IssmVec* pf, Parameters* parameters){/*{{{*/ + + /*Let matrix decide, to retain object orientation: */ + IssmVec* outvector=NULL; + + outvector=Kff->Solve(pf,parameters); + + *pout=outvector; +} +/*}}}*/ void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/