Changeset 23044


Ignore:
Timestamp:
08/02/18 16:17:04 (7 years ago)
Author:
erobo
Message:

CHG: allow for MUMPS serial

Location:
issm/trunk-jpl/src/c/toolkits
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/toolkits/issm/IssmDenseMat.h

    r22737 r23044  
    272272                        pf=(IssmSeqVec<IssmDouble>*)pfin;
    273273
     274                        #ifdef _HAVE_MUMPS_
     275                        /*Assume we have a sequential vec, downcast*/
     276                        uf=((IssmSeqVec<IssmDouble>*)pfin)->Duplicate();
     277                        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);
     278                        return uf;
     279                        #endif
     280
    274281                        #ifdef _HAVE_GSL_
    275282                        DenseGslSolve(/*output*/ &x,/*stiffness matrix:*/ this->matrix,this->M,this->N, /*right hand side load vector: */ pf->vector,pf->M,parameters);
     
    277284                        uf=new IssmSeqVec<IssmDouble>(x,this->N); xDelete(x);
    278285                        return uf;
    279                         #else
    280                                 _error_("GSL support not compiled in!");
    281286                        #endif
     287
     288                        _error_("No solver available");
     289                        return NULL;
    282290
    283291                }/*}}}*/
  • issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h

    r16415 r23044  
    1717#include "../ToolkitOptions.h"
    1818#include "./IssmToolkitUtils.h"
    19 
     19#include "../mumps/mumpsincludes.h"
    2020/*}}}*/
    2121
  • issm/trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

    r21615 r23044  
    2727        theMumpsStruc.par          = 1; 
    2828        theMumpsStruc.sym          = 0;
     29        #ifdef _HAVE_MPI_
    2930        theMumpsStruc.comm_fortran = MPI_Comm_c2f(IssmComm::GetComm());
     31        #endif
    3032        theMumpsStruc.job          = -1;
    3133        dmumps_c(&theMumpsStruc);
     
    3537void MumpsSettings(DMUMPS_STRUC_C &theMumpsStruc) {
    3638        /*Control statements:{{{ */
    37         theMumpsStruc.icntl[1-1] = 6; //error verbose: default 6, 0 or negative -> suppressed
     39        theMumpsStruc.icntl[1-1] = 0; //error verbose: default 6, 0 or negative -> suppressed
    3840        theMumpsStruc.icntl[2-1] = 0; //std verbose: default 1, 0 or negative -> suppressed
    3941        theMumpsStruc.icntl[3-1] = 0; //global information verbose: default 6, 0 or negative -> suppressed
     
    108110#endif
    109111
     112#ifdef _HAVE_MPI_
    110113void 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){ /*{{{*/
    111114
     
    201204        xDelete<int>(displs);
    202205} /*}}}*/
    203 
    204206void 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){ /*{{{*/
    205207
     
    287289        xDelete<int>(displs);
    288290} /*}}}*/
     291#endif
     292void 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){/*{{{*/
     293
     294        /*Variables*/
     295        int        *irn_loc = NULL;
     296        int        *jcn_loc = NULL;
     297        IssmDouble *a_loc   = NULL;
     298        IssmDouble *rhs     = NULL;
     299
     300        /*First, some checks*/
     301        if (Kff_M!=Kff_N)_error_("stiffness matrix Kff should be square");
     302        if (Kff_M!=Kff_m)_error_("stiffness matrix Kff is not serial");
     303        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");
     304        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");
     305
     306        /*Initialize matrix */
     307        /*figure out number of non-zero entries: */
     308        int nnz = 0;
     309        for(int i=0;i<Kff_M;i++){
     310                for(int j=0;j<Kff_N;j++){
     311                        if(Kff[i*Kff_N+j]!=0)nnz++;
     312                }
     313        }
     314
     315        /*Allocate: */
     316        if(nnz){
     317                irn_loc = xNew<int>(nnz);
     318                jcn_loc = xNew<int>(nnz);
     319                a_loc   = xNew<IssmDouble>(nnz);
     320        }
     321
     322        /*Populate the triplets: */
     323        int count=0;
     324        for(int i=0;i<Kff_M;i++){
     325                for(int j=0;j<Kff_N;j++){
     326                        if(Kff[i*Kff_N+j]!=0){
     327                                irn_loc[count] = i+1; //fortran indexing
     328                                jcn_loc[count] = j+1; //fortran indexing
     329                                a_loc[count]   = Kff[i*Kff_N+j];
     330                                count++;
     331                        }
     332                }
     333        }
     334
     335        /*Deal with right hand side. We need to ISSM_MPI_Gather it onto cpu 0: */
     336// AD performance is sensitive to calls to ensurecontiguous.
     337// Providing "t" will cause ensurecontiguous to be called.
     338#ifdef _HAVE_AD_
     339        rhs=xNew<IssmDouble>(pf_M,"t");
     340#else
     341        rhs=xNew<IssmDouble>(pf_M);
     342#endif
     343        xMemCpy<IssmDouble>(rhs,pf,Kff_M);
     344
     345        MumpsSolve(Kff_M,nnz,nnz,irn_loc,jcn_loc,a_loc,rhs,parameters);
     346
     347        /*Now scatter from cpu 0 to all other cpus*/
     348        xMemCpy<IssmDouble>(uf,rhs,Kff_M);
     349
     350        /*Cleanup*/
     351        xDelete<int>(irn_loc);
     352        xDelete<int>(jcn_loc);
     353        xDelete<IssmDouble>(a_loc);
     354        xDelete<IssmDouble>(rhs);
     355}/*}}}*/
    289356
    290357#ifdef _HAVE_ADOLC_
  • issm/trunk-jpl/src/c/toolkits/mumps/mumpsincludes.h

    r16415 r23044  
    1919template <class doubletype> class SparseRow;
    2020
     21#ifdef _HAVE_MPI_
    2122void 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);
    2223void 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);
     24#endif
     25void 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);
    2326
    2427#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_)
Note: See TracChangeset for help on using the changeset viewer.