Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16414) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16415) @@ -229,4 +229,5 @@ ./toolkits/issm/IssmSolver.h\ ./toolkits/issm/IssmSolver.cpp\ + ./toolkits/issm/SparseRow.h\ ./toolkits/issm/Bucket.h\ ./toolkits/mpi/issmmpi.h\ @@ -804,5 +805,5 @@ mumps_sources= ./toolkits/mumps\ ./toolkits/mumps/mumpsincludes.h\ - ./toolkits/mumps/MpiDenseMumpsSolve.cpp + ./toolkits/mumps/MumpsSolve.cpp #}}} #Gsl sources {{{ Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16414) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 16415) @@ -627,4 +627,5 @@ DenseEnum, MpiDenseEnum, + MpiSparseEnum, SeqEnum, MpiEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16414) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 16415) @@ -599,4 +599,5 @@ case DenseEnum : return "Dense"; case MpiDenseEnum : return "MpiDense"; + case MpiSparseEnum : return "MpiSparse"; case SeqEnum : return "Seq"; case MpiEnum : return "Mpi"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16414) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 16415) @@ -611,4 +611,5 @@ else if (strcmp(name,"Dense")==0) return DenseEnum; else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; + else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Seq")==0) return SeqEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h (revision 16414) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h (revision 16415) @@ -29,4 +29,5 @@ template class IssmDenseMat; template class IssmMpiDenseMat; +template class IssmMpiSparseMat; class Parameters; @@ -53,4 +54,11 @@ #endif break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); + #endif + break; default: _error_("matrix type not supported yet!"); @@ -72,4 +80,11 @@ #endif break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(M,N); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); + #endif + break; default: _error_("matrix type not supported yet!"); @@ -91,4 +106,11 @@ #endif break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(M,N,sparsity); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); + #endif + break; default: _error_("matrix type not supported yet!"); @@ -110,4 +132,11 @@ #endif break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(m,n,M,N,d_nnz,o_nnz); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); + #endif + break; default: _error_("matrix type not supported yet!"); @@ -129,4 +158,11 @@ #endif break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(serial_mat,M,N,sparsity); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); + #endif + break; default: _error_("matrix type not supported yet!"); @@ -147,4 +183,11 @@ #else _error_("MpiDense matrix requires compilation of MPI!"); + #endif + break; + case MpiSparseEnum: + #ifdef _HAVE_MPI_ + this->matrix=new IssmMpiSparseMat(M,N,connectivity,numberofdofspernode); + #else + _error_("MpiSparse matrix requires compilation of MPI!"); #endif break; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 16414) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiDenseMat.h (revision 16415) @@ -146,7 +146,7 @@ _printf_("cpu " << i << " #rows: " << this->m << "\n"); for (j=0;jm;j++){ - _printf_("row " << j << "\n"); + _printf_("row " << j << ":"); for (k=0;kN;k++){ - _printf_(this->matrix[j*this->N+k] << " \n"); + if(this->matrix[j*this->N+k]!=0)_printf_("(" << k << "," << this->matrix[j*this->N+k] << ") "); } _printf_("\n"); Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiSparseMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiSparseMat.h (revision 16415) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiSparseMat.h (revision 16415) @@ -0,0 +1,581 @@ +/*!\file: IssmMpiSparseMat.h + * \brief implementation of parallel sparse ISSM matrix. Internally, the parallel sparse matrix is + * split in rows across each cpu. Locally, on each cpu, the local matrix is represented by a vector of sparse rows. + * This object needs to answer the API defined by the virtual functions in IssmAbsMat, + * and the contructors required by IssmMat (see IssmMat.h) + */ + +#ifndef _ISSM_MPI_SPARSE_MAT_H_ +#define _ISSM_MPI_SPARSE_MAT_H_ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../../datastructures/datastructures.h" +#include "../../shared/shared.h" +#include "../mumps/mumpsincludes.h" +#include "./Bucket.h" +#include "./IssmMpiVec.h" +#include "./SparseRow.h" +#include + +/*}}}*/ + +/*We need to template this class, in case we want to create Matrices that hold + IssmDouble* matrix or IssmPDouble* matrix. + Such matrices would be useful for use without or with the matlab or python + interface (which do not care for IssmDouble types, but only rely on + IssmPDouble types)*/ + +template class IssmAbsMat; + +template +class IssmMpiSparseMat:public IssmAbsMat{ + + public: + + int M,N; //global size + int m; //local number of rows + SparseRow** matrix; /*here, doubletype is either IssmDouble or IssmPDouble*/ + DataSet* buckets; /*here, we store buckets of values that we will Assemble into a global matrix.*/ + /*IssmMpiSparseMat constructors, destructors*/ + /*FUNCTION IssmMpiSparseMat(){{{*/ + IssmMpiSparseMat(){ + this->M=0; + this->N=0; + this->m=0; + this->matrix=NULL; + this->buckets=new DataSet(); + } + /*}}}*/ + /*FUNCTION IssmMpiSparseMat(int M,int N){{{*/ + IssmMpiSparseMat(int Min,int Nin){ + this->Init(Min,Nin); + } + /*}}}*/ + /*FUNCTION IssmMpiSparseMat(int M,int N, doubletype sparsity){{{*/ + IssmMpiSparseMat(int pM,int pN, doubletype sparsity){ + /*no sparsity involved here, the sparsity pattern is resolve during the assemble phase: */ + this->Init(pM,pN); + } + /*}}}*/ + /*FUNCTION IssmMpiSparseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/ + IssmMpiSparseMat(int min,int nin,int Min,int Nin,int* d_nnz,int* o_nnz){ + + int i; + + /*no sparsity involved here, the sparsity pattern is resolved at the assemble phase: */ + this->buckets=new DataSet(); + + this->M=Min; + this->N=Nin; + this->m=min; + + /*Initialize pointer: */ + this->matrix=NULL; + + /*Allocate: */ + if (m*N){ + this->matrix=xNew*>(m); + for(i=0;imatrix[i]=new SparseRow(N); + } + } + } + /*}}}*/ + /*FUNCTION IssmMpiSparseMat(int M,int N, int connectivity, int numberofdofspernode){{{*/ + IssmMpiSparseMat(int pM,int pN, int connectivity,int numberofdofspernode){ + /*this is not needed, sparsity pattern is resolved at assemble phase: */ + this->Init(pM,pN); + } + /*}}}*/ + /*FUNCTION IssmMpiSparseMat::Init(int Min,int Nin){{{*/ + void Init(int Min,int Nin){ + + int i; + + this->buckets=new DataSet(); + + this->M=Min; + this->N=Nin; + + /*Figure out local number of rows: */ + this->m=DetermineLocalSize(this->M,IssmComm::GetComm()); + + /*Initialize pointer: */ + this->matrix=NULL; + + /*Allocate: */ + if (m*N){ + this->matrix=xNew*>(m); + for(i=0;imatrix[i]=new SparseRow(N); + } + } + } + /*}}}*/ + /*FUNCTION ~IssmMpiSparseMat(){{{*/ + ~IssmMpiSparseMat(){ + int i; + + if(m){ + for(i=0;i* sparserow=this->matrix[i]; + delete sparserow; this->matrix[i]=NULL; + } + xDelete*>(this->matrix); + } + M=0; + N=0; + m=0; + delete this->buckets; + } + /*}}}*/ + + /*IssmMpiSparseMat specific routines */ + /*FUNCTION Echo{{{*/ + void Echo(void){ + + int my_rank; + int i,j,k; + + /*Do a synchronized dump across all the rows: */ + my_rank=IssmComm::GetRank(); + for(i=0;im << "\n"); + for (j=0;jm;j++){ + _printf_("row " << j << ":"); + this->matrix[j]->Echo(); + _printf_("\n"); + } + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + } + + } + /*}}}*/ + /*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; + + int this_row_numvalues; + int* this_row_cols = NULL; + int* this_row_mods = NULL; + doubletype* this_row_values = NULL; + int row; + + /*some communicator info: */ + num_procs=IssmComm::GetSize(); + ISSM_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,&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 {{{ + *some scatter calls: */ + numvalues_fromcpu = xNew(num_procs); + for(i=0;i(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 ISSM_MPI_Scatterv prototype: + * int ISSM_MPI_Scatterv( void *sendbuf, int *sendcnts, int *displs, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcnt, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm) :*/ + sendcnts=xNew(num_procs); + displs=xNew(num_procs); + count=0; + for(i=0;im;row++){ + + /*figure out how many values this row is going to get: */ + count=0; + for(i=0;i(this_row_numvalues); + this_row_values=xNew(this_row_numvalues); + this_row_mods=xNew(this_row_numvalues); + + /*Collect: */ + count=0; + for(i=0;imatrix[row]->SetValues(this_row_numvalues,this_row_cols,this_row_values,this_row_mods); + + /*Free ressources: */ + xDelete(this_row_cols); + xDelete(this_row_values); + xDelete(this_row_mods); + } + /*}}}*/ + + /*Free ressources:{{{*/ + xDelete(RowRank); + xDelete(row_indices_forcpu); + xDelete(col_indices_forcpu); + xDelete(modes_forcpu); + xDelete(values_forcpu); + xDelete(numvalues_forcpu); + + for(i=0;i(bucketsforcpu); + + for(i=0;i(rows); + xDelete(cols); + xDelete(modes); + xDelete(values); + } + xDelete(row_indices_fromcpu); + xDelete(col_indices_fromcpu); + xDelete(modes_fromcpu); + xDelete(values_fromcpu); + xDelete(numvalues_fromcpu); + + xDelete(sendcnts); + xDelete(displs); + /*}}}*/ + + } + /*}}}*/ + /*FUNCTION Norm{{{*/ + doubletype Norm(NormMode mode){ + + doubletype norm,local_norm; + doubletype absolute; + int i,j; + + switch(mode){ + case NORM_INF: + local_norm=0; + for(i=0;im;i++){ + local_norm=max(local_norm,this->matrix[i]->Norm(mode)); + } + ISSM_MPI_Reduce(&local_norm, &norm, 1, ISSM_MPI_DOUBLE, ISSM_MPI_MAX, 0, IssmComm::GetComm()); + ISSM_MPI_Bcast(&norm,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + return norm; + break; + case NORM_FROB: + local_norm=0; + for(i=0;im;i++){ + local_norm+=this->matrix[i]->Norm(mode); + } + ISSM_MPI_Reduce(&local_norm, &norm, 1, ISSM_MPI_DOUBLE, ISSM_MPI_SUM, 0, IssmComm::GetComm()); + ISSM_MPI_Bcast(&norm,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + return sqrt(norm); + break; + + default: + _error_("unknown norm !"); + break; + } + } + /*}}}*/ + /*FUNCTION GetSize{{{*/ + void GetSize(int* pM,int* pN){ + *pM=M; + *pN=N; + } + /*}}}*/ + /*FUNCTION GetLocalSize{{{*/ + void GetLocalSize(int* pM,int* pN){ + *pM=m; + *pN=N; + } + /*}}}*/ + /*FUNCTION MatMult{{{*/ + void MatMult(IssmAbsVec* Xin,IssmAbsVec* AXin){ + + int i,j; + doubletype *X_serial = NULL; + + /*A check on the types: */ + if(IssmVecTypeFromToolkitOptions()!=MpiEnum)_error_("MatMult operation only possible with 'mpi' vectors"); + + /*Now that we are sure, cast vectors: */ + IssmMpiVec* X=(IssmMpiVec*)Xin; + IssmMpiVec* AX=(IssmMpiVec*)AXin; + + /*Serialize input Xin: */ + X_serial=X->ToMPISerial(); + + /*Every cpu has a serial version of the input vector. Use it to do the Matrix-Vector multiply + *locally and plug it into AXin: */ + for(i=0;im;i++){ + AX->vector[i]=this->matrix[i]->Mult(X_serial); + } + + /*Free ressources: */ + xDelete(X_serial); + } + /*}}}*/ + /*FUNCTION Duplicate{{{*/ + IssmMpiSparseMat* Duplicate(void){ + + _error_("not supported yet!"); + + } + /*}}}*/ + /*FUNCTION ToSerial{{{*/ + doubletype* ToSerial(void){ + _error_("not supported yet!"); + } + /*}}}*/ + /*FUNCTION SetValues{{{*/ + void SetValues(int min,int* idxm,int nin,int* idxn,doubletype* values,InsMode mode){ + + /*we need to store all the values we collect here in order to Assemble later. + * Indeed, the values we are collecting here most of the time will not belong + * to us, but to another part of the matrix on another cpu: */ + _assert_(buckets); + + buckets->AddObject(new Bucket(min,idxm,nin,idxn,values,mode)); + + } + /*}}}*/ + /*FUNCTION Convert{{{*/ + void Convert(MatrixType type){ + _error_("not supported yet!"); + } + /*}}}*/ + /*FUNCTION BucketsBuildScatterBuffers{{{*/ + void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,int** pcol_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; + int *temp_col_indices_forcpu = NULL; + doubletype *temp_values_forcpu = NULL; + int *temp_modes_forcpu = NULL; + + /*output: */ + int *numvalues_forcpu = NULL; + int *row_indices_forcpu = NULL; + int *col_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); + col_indices_forcpu = xNew(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_col_indices_forcpu=col_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_col_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_col_indices_forcpu!=col_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; + *pcol_indices_forcpu = col_indices_forcpu; + *pvalues_forcpu = values_forcpu; + *pmodes_forcpu = modes_forcpu; + } + /*}}}*/ + #ifndef _HAVE_WRAPPERS_ + /*Solve{{{*/ + IssmAbsVec* Solve(IssmAbsVec* pfin, Parameters* parameters){ + + /*output: */ + IssmMpiVec* uf=NULL; + IssmMpiVec* pf=NULL; + + /*Assume we are getting an IssmMpiVec in input, downcast: */ + pf=(IssmMpiVec*)pfin; + + switch(IssmSolverTypeFromToolkitOptions()){ + case MumpsEnum: { + /*Initialize output: */ + uf=pf->Duplicate(); + #ifdef _HAVE_MUMPS_ + MpiSparseMumpsSolve(/*output*/ uf->vector,uf->M,uf->m, /*stiffness matrix:*/ this->matrix,this->M,this->N,this->m, /*right hand side load vector: */ pf->vector,pf->M,pf->m,parameters); + #else + _error_("IssmMpiSparseMat solver requires MUMPS solver"); + #endif + return (IssmAbsVec*)uf; + } + case GslEnum: { + + IssmDouble* x=NULL; + #ifdef _HAVE_GSL_ + + _error_("not implemented yet!"); + //SparseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters); + + uf=new IssmMpiVec(x,this->N); xDelete(x); + + return (IssmAbsVec*)uf; + #else + _error_("GSL support not compiled in!"); + #endif + } + default: + _error_("solver type not supported yet!"); + } + + }/*}}}*/ + #endif +}; + +#endif //#ifndef _ISSM_MPI_SPARSE_MAT_H_ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp (revision 16414) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmToolkitUtils.cpp (revision 16415) @@ -38,4 +38,7 @@ if (strcmp(mat_type,"mpidense")==0){ mat_type_enum=MpiDenseEnum; + } + else if (strcmp(mat_type,"mpisparse")==0){ + mat_type_enum=MpiSparseEnum; } else if (strcmp(mat_type,"dense")==0){ Index: /issm/trunk-jpl/src/c/toolkits/issm/SparseRow.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/SparseRow.h (revision 16415) +++ /issm/trunk-jpl/src/c/toolkits/issm/SparseRow.h (revision 16415) @@ -0,0 +1,169 @@ +/*!\file: SparseRow: + * \brief implementation of a sparse row, which can then be used to make a sparse matrix. + */ + +#ifndef _SPARSE_ROW_H_ +#define _SPARSE_ROW_H_ + +/*Headers:*/ +#include "../toolkitsenums.h" +#include "../../shared/shared.h" + +template +class SparseRow{ + + public: + + int M; //real size + int ncols; //number of non-zeros + int* indices; + doubletype* values; + + /*SparseRow constructors, destructors*/ + SparseRow(){ /*{{{*/ + M=0; + ncols=0; + indices=NULL; + values=NULL; + } /*}}}*/ + SparseRow(int in_M){/*{{{*/ + + M=in_M; + ncols=0; + indices=NULL; + values=NULL; + } /*}}}*/ + ~SparseRow(){/*{{{*/ + if(ncols){ + xDelete(indices); + xDelete(values); + } + } /*}}}*/ + + /*SparseRow specific routines*/ + void Echo(){ /*{{{*/ + int i; + + for(i=0;i(indices); + xDelete(values); + } + + /*check numvalues: */ + if(!numvalues)return; + + /*Figure out number of non-zero values ncols: */ + buffer=xNewZeroInit(M); + for(i=0;i(ncols); _assert_(indices); + values=xNewZeroInit(ncols); _assert_(values); + + /*Go through buffer and add values: */ + count=0; + for(i=0;i(buffer); + + } /*}}}*/ + doubletype Norm(NormMode mode){ /*{{{*/ + + int i; + doubletype norm=0; + + switch(mode){ + case NORM_INF: + for(i=0;i +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../../shared/Numerics/types.h" +#include "../../shared/MemOps/MemOps.h" +#include "../../shared/Exceptions/exceptions.h" +#include "../../shared/io/Comm/IssmComm.h" +#include "../../classes/Params/GenericParam.h" +#include "../../classes/Params/Parameters.h" +#include "../mpi/issmmpi.h" +#include "../adolc/adolcincludes.h" +#include "../issm/SparseRow.h" +#include "./mumpsincludes.h" + +/*Mumps header files: */ +#include + +void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc){ + theMumpsStruc.par = 1; + theMumpsStruc.sym = 0; + theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm()); + theMumpsStruc.job = -1; + dmumps_c(&theMumpsStruc); +} + +// must be preceded by a call to MumpsInit +void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc) { + /*Control statements:{{{ */ + theMumpsStruc.icntl[1-1] = 6; //error verbose: default 6, 0 or negative -> suppressed + theMumpsStruc.icntl[2-1] = 0; //std verbose: default 1, 0 or negative -> suppressed + theMumpsStruc.icntl[3-1] = 0; //global information verbose: default 6, 0 or negative -> suppressed + theMumpsStruc.icntl[4-1] = 0; //verbose everything: default is 4 + theMumpsStruc.icntl[5-1] = 0; + theMumpsStruc.icntl[18-1] = 3; + + theMumpsStruc.icntl[20-1] = 0; + theMumpsStruc.icntl[21-1] = 0; + theMumpsStruc.icntl[30-1] = 0; + /*}}}*/ +} + +// must be preceded by a call to MumpsInit +void MumpsAnalyze(DMUMPS_STRUC_C &theMumpsStruc) { + theMumpsStruc.job = 1; + dmumps_c(&theMumpsStruc); +} + +// must be preceded by a call to MumpsAnalyze +void MumpsFactorize(DMUMPS_STRUC_C &theMumpsStruc) { + theMumpsStruc.job = 2; + dmumps_c(&theMumpsStruc); +} + +// must be preceded by a call to MumpsFactorize +void MumpsBacksubstitute(DMUMPS_STRUC_C &theMumpsStruc) { + theMumpsStruc.job = 3; + dmumps_c(&theMumpsStruc); +} + +// must be preceded at least by a call to MumpsInit +void MumpsFinalize(DMUMPS_STRUC_C &theMumpsStruc) { + theMumpsStruc.job = -2; + dmumps_c(&theMumpsStruc); +} + +void MumpsSolve(int n,int nnz,int local_nnz,int* irn_loc,int* jcn_loc, IssmPDouble *a_loc, IssmPDouble *rhs, Parameters* parameters=0 /*unused here*/){ + /*Initialize mumps*/ + DMUMPS_STRUC_C theMumpsStruc; + MumpsInit(theMumpsStruc); + MumpsSettings(theMumpsStruc); + + /*now setup the rest of theMumpsStruc */ + theMumpsStruc.n=n; + theMumpsStruc.nz=nnz; + theMumpsStruc.nz_loc=local_nnz; + theMumpsStruc.irn_loc=irn_loc; + theMumpsStruc.jcn_loc=jcn_loc; + theMumpsStruc.a_loc=a_loc; + theMumpsStruc.rhs=rhs; + theMumpsStruc.nrhs=1; + theMumpsStruc.lrhs=1; + + /*Solve system*/ + MumpsAnalyze(theMumpsStruc); + MumpsFactorize(theMumpsStruc); + MumpsBacksubstitute(theMumpsStruc); + MumpsFinalize(theMumpsStruc); +} + +#ifdef _HAVE_ADOLC_ +// prototype for active variant +void MumpsSolve(int n, + int nnz, + int local_nnz, + int* irn_loc, + int* jcn_loc, + IssmDouble *a_loc, + IssmDouble *rhs, + Parameters* parameters); +#endif + +void MpiDenseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ IssmDouble* Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){ /*{{{*/ + + /*Variables*/ + ISSM_MPI_Comm comm; + int num_procs; + int i,j; + int nnz,local_nnz; + int *irn_loc = NULL; + int *jcn_loc = NULL; + IssmDouble *a_loc = NULL; + int count; + int lower_row; + int upper_row; + IssmDouble *rhs = NULL; + int *recvcounts = NULL; + int *displs = NULL; + + /*Communicator info */ + num_procs = IssmComm::GetSize(); + comm = IssmComm::GetComm(); + + /*First, some checks*/ + if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square"); + if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf"); + if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf"); + + /*Initialize matrix */ + /*figure out number of non-zero entries: */ + local_nnz=0; + for(i=0;i(local_nnz); + jcn_loc=xNew(local_nnz); + a_loc=xNew(local_nnz); + } + + /*Populate the triplets: */ + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,Kff_m,comm); + count=0; + for(i=0;i(pf_M); + + recvcounts=xNew(num_procs); + displs=xNew(num_procs); + + /*recvcounts:*/ + ISSM_MPI_Allgather(&pf_m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,comm); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,comm); + + /*Gather:*/ + ISSM_MPI_Gatherv(pf, pf_m, ISSM_MPI_DOUBLE, rhs, recvcounts, displs, ISSM_MPI_DOUBLE,0,comm); + + MumpsSolve(Kff_M,nnz,local_nnz,irn_loc,jcn_loc,a_loc,rhs,parameters); + + /*Now scatter from cpu 0 to all other cpus*/ + ISSM_MPI_Scatterv( rhs, recvcounts, displs, ISSM_MPI_DOUBLE, uf, uf_m, ISSM_MPI_DOUBLE, 0, comm); + + /*Cleanup*/ + xDelete(irn_loc); + xDelete(jcn_loc); + xDelete(a_loc); + xDelete(rhs); + xDelete(recvcounts); + xDelete(displs); +} /*}}}*/ + +void MpiSparseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ SparseRow** Kff, int Kff_M, int Kff_N, int Kff_m, /*right hand side vector: */ IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){ /*{{{*/ + + /*Variables*/ + ISSM_MPI_Comm comm; + int num_procs; + int i,j; + int nnz,local_nnz; + int *irn_loc = NULL; + int *jcn_loc = NULL; + IssmDouble *a_loc = NULL; + int count; + int lower_row; + int upper_row; + IssmDouble *rhs = NULL; + int *recvcounts = NULL; + int *displs = NULL; + + /*Communicator info */ + num_procs = IssmComm::GetSize(); + comm = IssmComm::GetComm(); + + /*First, some checks*/ + if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square"); + if (uf_M!=Kff_M || uf_M!=pf_M)_error_("solution vector should be the same size as stiffness matrix Kff and load vector pf"); + if (uf_m!=Kff_m || uf_m!=pf_m)_error_("solution vector should be locally the same size as stiffness matrix Kff and load vector pf"); + + /*Initialize matrix */ + /*figure out number of non-zero entries: */ + local_nnz=0; + for(i=0;iNnz(); + } + + ISSM_MPI_Reduce(&local_nnz,&nnz,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,comm); + ISSM_MPI_Bcast(&nnz,1,ISSM_MPI_INT,0,comm); + + /*Allocate: */ + if(local_nnz){ + irn_loc=xNew(local_nnz); + jcn_loc=xNew(local_nnz); + a_loc=xNew(local_nnz); + } + + /*Populate the triplets: */ + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,Kff_m,comm); + count=0; + for(i=0;iSetIrnJcnA(irn_loc,jcn_loc,a_loc,lower_row+i+1,count); + count+=Kff[i]->Nnz(); + } + /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */ + rhs=xNew(pf_M); + + recvcounts=xNew(num_procs); + displs=xNew(num_procs); + + /*recvcounts:*/ + ISSM_MPI_Allgather(&pf_m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,comm); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,comm); + + /*Gather:*/ + ISSM_MPI_Gatherv(pf, pf_m, ISSM_MPI_DOUBLE, rhs, recvcounts, displs, ISSM_MPI_DOUBLE,0,comm); + + MumpsSolve(Kff_M,nnz,local_nnz,irn_loc,jcn_loc,a_loc,rhs,parameters); + + /*Now scatter from cpu 0 to all other cpus*/ + ISSM_MPI_Scatterv( rhs, recvcounts, displs, ISSM_MPI_DOUBLE, uf, uf_m, ISSM_MPI_DOUBLE, 0, comm); + + /*Cleanup*/ + xDelete(irn_loc); + xDelete(jcn_loc); + xDelete(a_loc); + xDelete(rhs); + xDelete(recvcounts); + xDelete(displs); +} /*}}}*/ + +#ifdef _HAVE_ADOLC_ + +int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) { + // unpack parameters + int n=iArr[0]; + int nnz=iArr[1]; + int local_nnz=iArr[2]; + int *local_irn=xNew(local_nnz); + int *local_jcn=xNew(local_nnz); + IssmPDouble *A=xNew(local_nnz); + for (int i=0;i(n); + for (int i=0;i(packedDimsSparseArrLength); + packedDimsSparseArr[0]=n; + packedDimsSparseArr[1]=nnz; + packedDimsSparseArr[2]=local_nnz; + for (int i=0;i(local_nnz+n); + for (int i=0;i(local_nnz+n); + IssmPDouble *passiveSol=xNew(n); + IssmDouble *sol=xNew(n); + call_ext_fct(dynamic_cast * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p, + packedDimsSparseArrLength, packedDimsSparseArr, + local_nnz+n, passivePack_A_rhs, pack_A_rhs, + n, passiveSol,sol); + for (int i=0;i(local_nnz); + int *local_jcn=xNew(local_nnz); + IssmPDouble *a_loc=xNew(local_nnz); + for (int i=0;i(n); + for (int i=0;i class SparseRow; void MpiDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_n, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters); +void MpiSparseMumpsSolve(IssmDouble* uf,int uf_M,int uf_n, SparseRow** Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters); #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16414) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 16415) @@ -591,4 +591,5 @@ def DenseEnum(): return StringToEnum("Dense")[0] def MpiDenseEnum(): return StringToEnum("MpiDense")[0] +def MpiSparseEnum(): return StringToEnum("MpiSparse")[0] def SeqEnum(): return StringToEnum("Seq")[0] def MpiEnum(): return StringToEnum("Mpi")[0] Index: /issm/trunk-jpl/src/m/enum/MpiSparseEnum.m =================================================================== --- /issm/trunk-jpl/src/m/enum/MpiSparseEnum.m (revision 16415) +++ /issm/trunk-jpl/src/m/enum/MpiSparseEnum.m (revision 16415) @@ -0,0 +1,11 @@ +function macro=MpiSparseEnum() +%MPISPARSEENUM - Enum of MpiSparse +% +% WARNING: DO NOT MODIFY THIS FILE +% this file has been automatically generated by src/c/shared/Enum/Synchronize.sh +% Please read src/c/shared/Enum/README for more information +% +% Usage: +% macro=MpiSparseEnum() + +macro=StringToEnum('MpiSparse'); Index: /issm/trunk-jpl/src/m/solvers/issmmumpssolver.m =================================================================== --- /issm/trunk-jpl/src/m/solvers/issmmumpssolver.m (revision 16414) +++ /issm/trunk-jpl/src/m/solvers/issmmumpssolver.m (revision 16415) @@ -11,5 +11,5 @@ %default issmoptions options issmoptions.toolkit='issm'; -issmoptions.mat_type=getfieldvalue(options,'mat_type','mpidense'); +issmoptions.mat_type=getfieldvalue(options,'mat_type','mpisparse'); issmoptions.vec_type=getfieldvalue(options,'vec_type','mpi'); issmoptions.solver_type=getfieldvalue(options,'solver_type','mumps');