Changeset 14864

Timestamp:

05/02/13 18:39:09 (12 years ago)

Author:

Eric.Larour

Message:

NEW: implemented a new hook up for the ISSM toolkit to the MUMPS solver, witout recourse to the Petsc
library. Needs V&V still, but I can get test101 to run. No success yet.

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (2 diffs)
• toolkits/issm/IssmMat.h (modified) (1 diff)
• toolkits/issm/IssmMpiDenseMat.h (modified) (3 diffs)
• toolkits/issm/IssmMpiVec.h (modified) (4 diffs)
• toolkits/issm/IssmVec.h (modified) (1 diff)
• toolkits/mumps (added)
• toolkits/mumps/MpiDenseMumpsSolve.cpp (added)
• toolkits/mumps/mumpsincludes.h (added)

issm/trunk-jpl/src/c/Makefile.am

@@ -768 +768 @@
 
 #}}}
+#Mumps sources  {{{
+mumps_sources=      ./toolkits/mumps\
+                                        ./toolkits/mumps/mumpsincludes.h\
+                                        ./toolkits/mumps/MpiDenseMumpsSolve.cpp
+#}}}
 #Mpi sources  {{{
 mpi_sources= ./toolkits/mpi/mpiincludes.h\
@@ -866 +871 @@
 endif
 
+if MUMPS
+issm_sources  +=  $(mumps_sources)
+endif
+
 if TRANSIENT
 issm_sources  +=  $(transient_sources)

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

@@ -208 +208 @@
                         outvector->vector=this->matrix->Solve(pf->vector,parameters);
 
+                        return outvector;
+
                 }/*}}}*/
 

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

@@ -25 +25 @@
 #include "../../classes/IssmComm.h"
 #include "../../classes/objects/Bucket.h"
-#include "../../toolkits/toolkits.h"
+#include "../mumps/mumpsincludes.h"
+#include "./IssmMpiVec.h"
 #include <math.h>
 
@@ -69 +70 @@
                 /*}}}*/
                 /*FUNCTION IssmMpiDenseMat(int m,int n,int M,int N,int* d_nnz,int* o_nnz){{{*/
-                IssmMpiDenseMat(int m,int n,int pM,int pN,int* d_nnz,int* o_nnz){
+                IssmMpiDenseMat(int min,int nin,int Min,int Nin,int* d_nnz,int* o_nnz){
                         /*not needed, we are fully dense!: */
-                        this->Init(pM,pN);
+                        
+                        this->buckets=new DataSet();
+
+                        this->M=Min;
+                        this->N=Nin;
+                        this->m=min;
+
+                        /*Initialize pointer: */
+                        this->matrix=NULL;
+
+                        /*Allocate: */
+                        if (m*N)this->matrix=xNewZeroInit<doubletype>(this->m*N);
                 }
                 /*}}}*/
@@ -604 +616 @@
                 IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pfin, Parameters* parameters){
 
+                        int i;
+                        
+                        /*output: */
+                        IssmMpiVec<IssmDouble>* uf=NULL;
+                        IssmMpiVec<IssmDouble>* pf=NULL;
+
                         /*Assume we are getting an IssmMpiVec in input, downcast: */
-                        IssmMpiVec<IssmDouble>* pf=(IssmMpiVec<IssmDouble>*)pfin;
+                        pf=(IssmMpiVec<IssmDouble>*)pfin;
+
+                        /*Initialize output: */
+                        uf=pf->Duplicate();
 
                         /*Let's try and use the MUMPS solver here: */
                         #ifdef _HAVE_MUMPS_
-
+                        MpiDenseMumpsSolve(/*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);
                         #else
                         _error_("IssmMpiDenseMat solver requires MUMPS solver");
                         #endif
+                        return (IssmAbsVec<IssmDouble>*)uf;
 
                 }/*}}}*/

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

@@ -60 +60 @@
                 }
                 /*}}}*/
+                /*FUNCTION IssmMpiVec(int m,int M){{{*/
+                IssmMpiVec(int min, int Min){
+                        this->Init(min,true);
+                }
+                /*}}}*/
                 /*FUNCTION IssmMpiVec(int M,bool fromlocalsize){{{*/
                 IssmMpiVec(int Min, bool fromlocalsize){
@@ -71 +76 @@
 
                         if(this->M){
+                                this->vector=xNew<doubletype>(this->m);
+                                xMemCpy<doubletype>(this->vector,buffer,this->m);
+                        }
+                }
+                /*}}}*/
+                /*FUNCTION IssmMpiVec(doubletype* serial_vec,int M,int m){{{*/
+                IssmMpiVec(doubletype* buffer,int Min,int min){
+
+                        this->vector=NULL;
+                        this->buckets=new DataSet();
+                        this->M=Min;
+                        this->m=min;
+
+                        if(this->m){
                                 this->vector=xNew<doubletype>(this->m);
                                 xMemCpy<doubletype>(this->vector,buffer,this->m);
@@ -119 +138 @@
                                                 printf("row %i %g",j,this->vector[j]);
                                         }
+                                        printf("\n");
                                 }
                                 MPI_Barrier(IssmComm::GetComm());
@@ -437 +457 @@
                 IssmMpiVec<doubletype>* Duplicate(void){
 
-                        return new IssmMpiVec<doubletype>(this->vector,this->M);
+                        return new IssmMpiVec<doubletype>(this->vector,this->M,this->m);
 
                 }

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

@@ -42 +42 @@
                 /*IssmVec constructors, destructors*/
                 IssmVec(){ /*{{{*/
-
-                        switch(IssmVecTypeFromToolkitOptions()){
-
-                                case SeqEnum: 
-                                        this->vector=new IssmSeqVec<doubletype>();
-                                        break;
-                                case MpiEnum:
-                                        #ifdef _HAVE_MPI_
-                                        this->vector=new IssmMpiVec<doubletype>();
-                                        #else
-                                        _error_("Mpi vector requires compilation of MPI!");
-                                        #endif
-                                        break;
-                                default:
-                                        _error_("vector type not supported yet!");
-                        }
+                        this->vector=NULL;
                 }
                 /*}}}*/