Changeset 14792

Timestamp:

04/29/13 15:36:48 (12 years ago)

Author:

Eric.Larour

Message:

CHG: finished implementation of IssmMpiVec and IssmMpiDenseMat objects. Still need to V&V.
CHG: moved the src/c/classes/matrix/Matrix and Vector objects to src/c/classes/toolkits.
CHG and NEW: offloaded the solver code to the toolkits. To do so, created a Solver class
in the src/c/classes/toolkits, which connects the toolkits to the abstract classes Matrix
and Vector.

Location:

issm/trunk-jpl/src/c

Files:

• Makefile.am (modified) (4 diffs)
• classes/classes.h (modified) (1 diff)
• classes/matrix/matrixobjects.h (modified) (1 diff)
• classes/toolkits (added)
• classes/toolkits/Matrix.h (moved) (moved from issm/trunk-jpl/src/c/classes/matrix/Matrix.h ) (1 diff)
• classes/toolkits/Solver.h (added)
• classes/toolkits/Vector.h (moved) (moved from issm/trunk-jpl/src/c/classes/matrix/Vector.h )
• classes/toolkits/toolkitobjects.h (added)
• modules/Solverx/CMakeLists.txt (modified) (1 diff)
• modules/Solverx/DofTypesToIndexSet.cpp (deleted)
• modules/Solverx/Solverx.cpp (modified) (2 diffs)
• modules/Solverx/Solverx.h (modified) (1 diff)
• modules/Solverx/SolverxPetsc.cpp (deleted)
• modules/Solverx/SolverxSeq.cpp (deleted)
• toolkits/issm/IssmMpiDenseMat.h (modified) (2 diffs)
• toolkits/issm/IssmMpiVec.h (modified) (4 diffs)
• toolkits/issm/IssmSolver.cpp (added)
• toolkits/issm/IssmSolver.h (added)
• toolkits/issm/issmtoolkit.h (modified) (1 diff)
• toolkits/petsc/objects/PetscSolver.cpp (added)
• toolkits/petsc/objects/PetscSolver.h (added)
• toolkits/petsc/objects/petscobjects.h (modified) (1 diff)

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

@@ -119 +119 @@
                                         ./classes/matrix/ElementVector.h\
                                         ./classes/matrix/ElementVector.cpp\
-                                        ./classes/matrix/Matrix.h\
-                                        ./classes/matrix/Vector.h\
+                                        ./classes/toolkits/toolkitobjects.h\
+                                        ./classes/toolkits/Matrix.h\
+                                        ./classes/toolkits/Vector.h\
+                                        ./classes/toolkits/Solver.h\
                                         ./classes/objects/Params/Param.h\
                                         ./classes/objects/Params/GenericParam.h\
@@ -229 +231 @@
                                         ./toolkits/issm/IssmSeqVec.h\
                                         ./toolkits/issm/IssmVec.h\
+                                        ./toolkits/issm/IssmSolver.h\
+                                        ./toolkits/issm/IssmSolver.cpp\
                                         ./toolkits/triangle/triangleincludes.h\
                                         ./toolkitsenums.h\
@@ -320 +324 @@
                                         ./modules/Solverx/Solverx.cpp\
                                         ./modules/Solverx/Solverx.h\
-                                        ./modules/Solverx/SolverxSeq.cpp\
                                         ./modules/VecMergex/VecMergex.cpp\
                                         ./modules/VecMergex/VecMergex.h\
@@ -760 +763 @@
                                         ./toolkits/petsc/objects/PetscVec.h\
                                         ./toolkits/petsc/objects/PetscVec.cpp\
-                                        ./toolkits/petsc/petscincludes.h\
-                                        ./modules/Solverx/SolverxPetsc.cpp\
-                                        ./modules/Solverx/DofTypesToIndexSet.cpp
+                                        ./toolkits/petsc/objects/PetscSolver.cpp\
+                                        ./toolkits/petsc/objects/PetscSolver.h\
+                                        ./toolkits/petsc/petscincludes.h
 
 #}}}

issm/trunk-jpl/src/c/classes/classes.h

@@ -11 +11 @@
 /*matrix: */
 #include "./matrix/matrixobjects.h"
+
+/*toolkit objects: */
+#include "./toolkits/toolkitobjects.h"
 
 /*bamg: */

issm/trunk-jpl/src/c/classes/matrix/matrixobjects.h

@@ -9 +9 @@
 #include "./ElementMatrix.h"
 #include "./ElementVector.h"
-#include "./Vector.h"
-#include "./Matrix.h"
 
 #endif

issm/trunk-jpl/src/c/classes/toolkits/Matrix.h

@@ -324 +324 @@
                 }
                 /*}}}*/
+
 };
 

issm/trunk-jpl/src/c/modules/Solverx/CMakeLists.txt

@@ -5 +5 @@
 # }}}
 # CORE_SOURCES {{{
-set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp
-              $ENV{ISSM_DIR}/src/c/modules/Solverx/SolverxSeq.cpp PARENT_SCOPE)
+set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/modules/Solverx/Solverx.cpp PARENT_SCOPE)
 # }}}

issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp

@@ -2 +2 @@
  * \brief solver
  */
-
-#include "./Solverx.h"
-#include "../../shared/shared.h"
-#include "../../include/include.h"
-#include "../../io/io.h"
 
 #ifdef HAVE_CONFIG_H
@@ -14 +9 @@
 #endif
 
+#include "./Solverx.h"
+#include "../../shared/shared.h"
+
 void    Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters){
 
+        /*intermediary: */
+        Solver<IssmDouble> *solver=NULL;
+        
         /*output: */
         Vector<IssmDouble> *uf=NULL;
 
-        /*In debugging mode, check that stiffness matrix and load vectors are not NULL (they can be empty)*/
-        _assert_(Kff);
-        _assert_(pf);
-        
         if(VerboseModule()) _pprintLine_("   Solving matrix system");
 
-        /*Initialize vector: */
-        uf=new Vector<IssmDouble>();
+        /*Initialize solver: */
+        solver=new Solver<IssmDouble>(Kff,pf,uf0,df,parameters);
 
-        /*According to matrix type, use specific solvers: */
-        switch(Kff->type){
-                #ifdef _HAVE_PETSC_
-                case PetscMatType:{
-                        PetscVec* uf0_vector = NULL;
-                        PetscVec* df_vector  = NULL;
-                        if(uf0) uf0_vector = uf0->pvector;
-                        if(df)  df_vector  = df->pvector;
-                        SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
-                        break;}
-                #endif
-                case IssmMatType:{
-                        SolverxSeq(&uf->ivector,Kff->imatrix,pf->ivector,parameters);
-                        break;}
-                default:
-                          _error_("Matrix type: " << Kff->type << " not supported yet!");
-        }
+        /*Solve:*/
+        uf=solver->Solve();
 
         /*Assign output pointers:*/

issm/trunk-jpl/src/c/modules/Solverx/Solverx.h

@@ -12 +12 @@
 #endif
 
-#include "../../classes/objects/objects.h"
-#include "../../toolkits/toolkits.h"
+#include "../../classes/toolkits/toolkitobjects.h"
 
 /* local prototypes: */
 void    Solverx(Vector<IssmDouble>** puf, Matrix<IssmDouble>* Kff, Vector<IssmDouble>* pf, Vector<IssmDouble>* uf0,Vector<IssmDouble>* df, Parameters* parameters);
 
-#ifdef _HAVE_PETSC_
-void    SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
-void    SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
-void  DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
-#endif
-
-void SolverxSeq(IssmVec<IssmDouble>** puf,IssmMat<IssmDouble>* Kff, IssmVec<IssmDouble>* pf,Parameters* parameters);
-void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
-void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
-
-#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters);
-// call back functions:
-ADOLC_ext_fct EDF_for_solverx;
-ADOLC_ext_fct_fos_forward EDF_fos_forward_for_solverx;
-ADOLC_ext_fct_fos_reverse EDF_fos_reverse_for_solverx;
-ADOLC_ext_fct_fov_forward EDF_fov_forward_for_solverx;
-ADOLC_ext_fct_fov_reverse EDF_fov_reverse_for_solverx;
-#endif
-
 #endif  /* _SOLVERX_H */

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

@@ -272 +272 @@
                 /*FUNCTION Norm{{{*/
                 doubletype Norm(NormMode mode){
-                        _error_("not supported yet!");
+                        
+                        
+                        doubletype norm,local_norm;
+                        doubletype absolute;
+                        int i,j;
+
+                        switch(mode){
+                                case NORM_INF:
+                                        local_norm=0;
+                                        for(i=0;i<this->M;i++){
+                                                absolute=0;
+                                                for(j=0;j<this->N;j++){
+                                                        absolute+=fabs(this->matrix[N*i+j]);
+                                                }
+                                                local_norm=max(local_norm,absolute);
+                                        }
+                                        MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
+                                        MPI_Bcast(&norm,1,MPI_DOUBLE,0,IssmComm::GetComm());
+                                        return norm;
+                                        break; 
+                                default:
+                                        _error_("unknown norm !");
+                                        break;
+                        }
                 }
                 /*}}}*/
                 /*FUNCTION GetSize{{{*/
                 void GetSize(int* pM,int* pN){
-                        _error_("not supported yet!");
+                        *pM=M;
+                        *pN=N;
                 }
                 /*}}}*/
                 /*FUNCTION GetLocalSize{{{*/
                 void GetLocalSize(int* pM,int* pN){
-                        _error_("not supported yet!");
+                        *pM=m;
+                        *pN=N;
                 }
                 /*}}}*/
                 /*FUNCTION MatMult{{{*/
                 void MatMult(IssmAbsVec<doubletype>* Xin,IssmAbsVec<doubletype>* AXin){
-                        _error_("not supported yet!");
+
+
+                        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<doubletype>* X=(IssmMpiVec<doubletype>*)Xin;
+                        IssmMpiVec<doubletype>* AX=(IssmMpiVec<doubletype>*)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;i<this->m;i++){
+                                for(j=0;j<this->N;j++){
+                                        AX->vector[i]+=this->matrix[i*N+j]*X_serial[j];
+                                }
+                        }
+
+                        /*Free ressources: */
+                        xDelete<doubletype>(X_serial);
                 }
                 /*}}}*/
                 /*FUNCTION Duplicate{{{*/
                 IssmMpiDenseMat<doubletype>* Duplicate(void){
-                        _error_("not supported yet!");
+
+                        IssmMpiDenseMat<doubletype>* dup=new IssmMpiDenseMat<doubletype>(this->matrix,this->M,this->N,0);
+                        return dup;
+
                 }
                 /*}}}*/
@@ -314 +367 @@
                 /*FUNCTION Convert{{{*/
                 void Convert(MatrixType type){
-                        /*do nothing: */
+                        _error_("not supported yet!");
                 }
                 /*}}}*/         

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

@@ -22 +22 @@
 #include "../../shared/Alloc/alloc.h"
 #include "../../include/macros.h"
+#include "../../io/io.h"
 #ifdef _HAVE_MPI_
 #include "../mpi/mpiincludes.h"
@@ -303 +304 @@
                 /*FUNCTION AXPY{{{*/
                 void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){
-                        printf("AXPY not implemented yet!");
-                        exit(1);
+                                
+                        int i;
+
+                        /*Assume X is of the correct type, and downcast: */
+                        IssmMpiVec* X=NULL;
+
+                        X=(IssmMpiVec<doubletype>*)Xin;
+
+                        /*y=a*x+y where this->vector is y*/
+                        for(i=0;i<this->m;i++)this->vector[i]=a*X->vector[i]+this->vector[i];
+
+
                 }
                 /*}}}*/
                 /*FUNCTION AYPX{{{*/
                 void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){
-                        printf("AYPX not implemented yet!");
-                        exit(1);
+                        int i;
+
+                        /*Assume X is of the correct type, and downcast: */
+                        IssmMpiVec* X=NULL;
+
+                        X=(IssmMpiVec<doubletype>*)Xin;
+
+                        /*y=x+a*y where this->vector is y*/
+                        for(i=0;i<this->m;i++)this->vector[i]=X->vector[i]+a*this->vector[i];
+
                 }
                 /*}}}*/
                 /*FUNCTION ToMPISerial{{{*/
                 doubletype* ToMPISerial(void){
-                        
-                        printf("IssmMpiVec ToMPISerial not implemented yet!");
-                        exit(1);
+
+                        /*communicator info: */
+                        MPI_Comm comm;
+                        int num_procs;
+
+                        /*MPI_Allgatherv info: */
+                        int  lower_row,upper_row;
+                        int* recvcounts=NULL;
+                        int* displs=NULL;
+                        
+                        /*output: */
+                        doubletype* buffer=NULL;
+
+                        /*initialize comm info: */
+                        comm=IssmComm::GetComm();
+                        num_procs=IssmComm::GetSize();
+
+                        /*Allocate: */
+                        buffer=xNew<doubletype>(M);
+                        recvcounts=xNew<int>(num_procs);
+                        displs=xNew<int>(num_procs);
+
+                        /*recvcounts:*/
+                        MPI_Allgather(&this->m,1,MPI_INT,recvcounts,1,MPI_INT,comm);
+                        
+                        /*get lower_row: */
+                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,this->m,comm);
+
+                        /*displs: */
+                        MPI_Allgather(&lower_row,1,MPI_INT,displs,1,MPI_INT,comm);
+
+                        /*All gather:*/
+                        MPI_Allgatherv(this->vector, this->m, MPI_DOUBLE, buffer, recvcounts, displs, MPI_DOUBLE,comm);
+
+                        /*free ressources: */
+                        xDelete<int>(recvcounts);
+                        xDelete<int>(displs);
+
+                        /*return: */
+                        return buffer;
 
                 }
@@ -338 +394 @@
                 /*FUNCTION Norm{{{*/
                 doubletype Norm(NormMode mode){
-
-                        printf("IssmMpiVec Norm not implemented yet!");
-                        exit(1);
+                
+                        doubletype local_norm;
+                        doubletype norm;
+                        int i;
+
+                        switch(mode){
+                                case NORM_INF:
+                                        //local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,fabs(this->vector[i]));
+                                        local_norm=0; for(i=0;i<this->m;i++)local_norm=max(local_norm,this->vector[i]);
+                                        MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_MAX, 0, IssmComm::GetComm());
+                                        return norm;
+                                        break;
+                                case NORM_TWO:
+                                        local_norm=0; 
+                                        for(i=0;i<this->m;i++)local_norm+=pow(this->vector[i],2);
+                                        MPI_Reduce(&local_norm, &norm, 1, MPI_DOUBLE, MPI_SUM, 0, IssmComm::GetComm());
+                                        return sqrt(norm);
+                                        break;
+                                default:
+                                        _error_("unknown norm !");
+                                        break;
+                        }
                 }
                 /*}}}*/
@@ -392 +467 @@
                 /*}}}*/
 };
-#endif //#ifndef _ISSM_MPI_VEC_H_
+#endif //#ifndef _ISSM_MPI_VEC_H_       

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

@@ -18 +18 @@
 #include "./IssmSeqVec.h"
 #include "./IssmVec.h"
+#include "./IssmSolver.h"
 
 #ifdef _HAVE_MPI_

issm/trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h

@@ -8 +8 @@
 #include "./PetscMat.h"
 #include "./PetscVec.h"
+#include "./PetscSolver.h"
 
 #endif