source: issm/oecreview/Archive/22819-23185/ISSM-23043-23044.diff

../trunk-jpl/src/c/toolkits/issm/IssmMat.h

@@ -16 +16 @@
 #include "../../shared/shared.h"
 #include "../ToolkitOptions.h"
 #include "./IssmToolkitUtils.h"
-
+#include "../mumps/mumpsincludes.h"
 /*}}}*/
 
 /*We need to template this class, in case we want to create Matrices that hold

../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h

@@ -18 +18 @@
 class Parameters;
 template <class doubletype> class SparseRow;
 
+#ifdef _HAVE_MPI_
 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<IssmDouble>** Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters);
+#endif
+void SeqDenseMumpsSolve(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);
 
 #if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
 // call back functions:

../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h

@@ -271 +271 @@
                         /*Assume we are getting an IssmMpiVec in input, downcast: */
                         pf=(IssmSeqVec<IssmDouble>*)pfin;
 
+                        #ifdef _HAVE_MUMPS_
+                        /*Assume we have a sequential vec, downcast*/
+                        uf=((IssmSeqVec<IssmDouble>*)pfin)->Duplicate();
+                        SeqDenseMumpsSolve(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);
+                        return uf;
+                        #endif
+
                         #ifdef _HAVE_GSL_
                         DenseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters);
 
                         uf=new IssmSeqVec<IssmDouble>(x,this->N); xDelete(x);
                         return uf;
-                        #else
-                                _error_("GSL support not compiled in!");
                         #endif
 
+                        _error_("No solver available");
+                        return NULL;
+
                 }/*}}}*/
                 #endif
 };

../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

@@ -26 +26 @@
 void MumpsInit(DMUMPS_STRUC_C &theMumpsStruc){ 
         theMumpsStruc.par          = 1;  
         theMumpsStruc.sym          = 0;
+        #ifdef _HAVE_MPI_
         theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm());
+        #endif
         theMumpsStruc.job          = -1;
         dmumps_c(&theMumpsStruc);
 }
@@ -34 +36 @@
 // 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[1-1] = 0; //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
@@ -107 +109 @@
                 Parameters* parameters);
 #endif 
 
+#ifdef _HAVE_MPI_
 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*/
@@ -200 +203 @@
         xDelete<int>(recvcounts);
         xDelete<int>(displs);
 } /*}}}*/
-
 void MpiSparseMumpsSolve( /*output: */ IssmDouble* uf, int uf_M, int uf_m, /*matrix input: */ SparseRow<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*/
@@ -286 +288 @@
         xDelete<int>(recvcounts);
         xDelete<int>(displs);
 } /*}}}*/
+#endif
+void SeqDenseMumpsSolve(IssmDouble* uf,int uf_M,int uf_m, IssmDouble* Kff,int Kff_M, int Kff_N, int Kff_m, IssmDouble* pf, int pf_M, int pf_m, Parameters* parameters){/*{{{*/
 
+        /*Variables*/
+        int        *irn_loc = NULL;
+        int        *jcn_loc = NULL;
+        IssmDouble *a_loc   = NULL;
+        IssmDouble *rhs     = NULL;
+
+        /*First, some checks*/
+        if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
+        if (Kff_M!=Kff_m)_error_("stiffness matrix Kff is not serial");
+        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: */
+        int nnz = 0;
+        for(int i=0;i<Kff_M;i++){
+                for(int j=0;j<Kff_N;j++){
+                        if(Kff[i*Kff_N+j]!=0)nnz++;
+                }
+        }
+
+        /*Allocate: */
+        if(nnz){
+                irn_loc = xNew<int>(nnz);
+                jcn_loc = xNew<int>(nnz);
+                a_loc   = xNew<IssmDouble>(nnz);
+        }
+
+        /*Populate the triplets: */
+        int count=0;
+        for(int i=0;i<Kff_M;i++){
+                for(int j=0;j<Kff_N;j++){
+                        if(Kff[i*Kff_N+j]!=0){
+                                irn_loc[count] = i+1; //fortran indexing
+                                jcn_loc[count] = j+1; //fortran indexing
+                                a_loc[count]   = Kff[i*Kff_N+j];
+                                count++;
+                        }
+                }
+        }
+
+        /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
+// AD performance is sensitive to calls to ensurecontiguous.
+// Providing "t" will cause ensurecontiguous to be called.
+#ifdef _HAVE_AD_
+        rhs=xNew<IssmDouble>(pf_M,"t");
+#else
+        rhs=xNew<IssmDouble>(pf_M);
+#endif
+        xMemCpy<IssmDouble>(rhs,pf,Kff_M);
+
+        MumpsSolve(Kff_M,nnz,nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
+
+        /*Now scatter from cpu 0 to all other cpus*/
+        xMemCpy<IssmDouble>(uf,rhs,Kff_M);
+
+        /*Cleanup*/
+        xDelete<int>(irn_loc);
+        xDelete<int>(jcn_loc);
+        xDelete<IssmDouble>(a_loc);
+        xDelete<IssmDouble>(rhs);
+}/*}}}*/
+
 #ifdef _HAVE_ADOLC_
 
 int mumpsSolveEDF(int iArrLength, int* iArr, int /* ignored */, IssmPDouble* dp_x, int /* ignored */, IssmPDouble* dp_y) {