Changeset 23044 for issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

Timestamp:

08/02/18 16:17:04 (7 years ago)

Author:

erobo

Message:

CHG: allow for MUMPS serial

File:

• issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp (modified) (5 diffs)

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

@@ -27 +27 @@
         theMumpsStruc.par          = 1;  
         theMumpsStruc.sym          = 0;
+        #ifdef _HAVE_MPI_
         theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm());
+        #endif
         theMumpsStruc.job          = -1;
         dmumps_c(&theMumpsStruc);
@@ -35 +37 @@
 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
@@ -108 +110 @@
 #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){ /*{{{*/
 
@@ -201 +204 @@
         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){ /*{{{*/
 
@@ -287 +289 @@
         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_