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

Index: ../trunk-jpl/src/c/toolkits/issm/IssmMat.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmMat.h  (revision 23043)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmMat.h  (revision 23044)
@@ -16,7 +16,7 @@
 #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
Index: ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h   (revision 23043)
+++ ../trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h   (revision 23044)
@@ -18,8 +18,11 @@
 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:
Index: ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h
===================================================================
--- ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h     (revision 23043)
+++ ../trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h     (revision 23044)
@@ -271,15 +271,23 @@
                        /*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
 };
Index: ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp
===================================================================
--- ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp    (revision 23043)
+++ ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp    (revision 23044)
@@ -26,7 +26,9 @@
 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,7 +36,7 @@
 // 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,6 +109,7 @@
                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,7 +203,6 @@
        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,7 +288,72 @@
        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) {