Changeset 12850


Ignore:
Timestamp:
08/01/12 13:27:31 (13 years ago)
Author:
Eric.Larour
Message:

CHG: homoegeization of the matrix class. Is now split into a pector (petsc vector) and svector (sequential vector).
New objects PetscMat and PetscVec in the toolkits allow for this homoegeization. Things are not perfect yet,
as we reference the matrix and vector internals across the code, which could create some serious issues.

Location:
issm/trunk-jpl/src/c
Files:
6 added
14 edited
3 moved

Legend:

Unmodified
Added
Removed

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

    r12835 r12850  
    331331                                        ./modules/Solverx/Solverx.cpp\
    332332                                        ./modules/Solverx/Solverx.h\
     333                                        ./modules/Solverx/SolverxSeq.cpp\
    333334                                        ./modules/VecMergex/VecMergex.cpp\
    334335                                        ./modules/VecMergex/VecMergex.h\
    … …  
    784785                                        ./toolkits/petsc/patches/ISSMToPetscInsertMode.cpp\
    785786                                        ./toolkits/petsc/patches/ISSMToPetscNormMode.cpp\
     787                                        ./toolkits/petsc/objects/petscobjects.h\
     788                                        ./toolkits/petsc/objects/PetscMat.h\
     789                                        ./toolkits/petsc/objects/PetscMat.cpp\
     790                                        ./toolkits/petsc/objects/PetscVec.h\
     791                                        ./toolkits/petsc/objects/PetscVec.cpp\
    786792                                        ./toolkits/petsc/petscincludes.h\
    787793                                        ./shared/Numerics/PetscOptionsFromAnalysis.cpp\
    … …  
    789795                                        ./modules/Solverx/DofTypesToIndexSet.cpp
    790796
    791 #}}}
    792 #Gsl sources  {{{
    793 gsl_sources= ./modules/Solverx/SolverxGsl.cpp
    794797#}}}
    795798#Mpi sources  {{{
    … …  
    831834#}}}
    832835#Matlab and Petsc sources  {{{
    833 matlabpetsc_sources= ./matlab/io/MatlabMatrixToPetscMatrix.cpp\
    834                                          ./matlab/io/MatlabVectorToPetscVector.cpp
     836matlabpetsc_sources= ./matlab/io/MatlabMatrixToPetscMat.cpp\
     837                                         ./matlab/io/MatlabVectorToPetscVec.cpp
    835838
    836839#}}}
    … …  
    925928endif
    926929
    927 if GSL
    928 issm_sources  +=  $(gsl_sources)
    929 endif
    930 
    931930if TRANSIENT
    932931issm_sources  +=  $(transient_sources)

  • issm/trunk-jpl/src/c/classes/matrix/Matrix.cpp

    r12832 r12850  
    2626Matrix::Matrix(){
    2727
    28         #ifdef _HAVE_PETSC_
    29         this->matrix=NULL;
    30         #else
    31         this->matrix=NULL;
     28        pmatrix=NULL;
     29        smatrix=NULL;
     30
     31        type=PetscMatType; //default
     32        #ifndef _HAVE_PETSC_
     33        type=SeqMatType;
    3234        #endif
    33         #ifdef _HAVE_ADOLC_
    34         this->amatrix=NULL;
    35         #endif
    36 }
    37 /*}}}*/
    38 /*FUNCTION Matrix::Matrix(int M,int N){{{*/
    39 Matrix::Matrix(int M,int N){
    40 
    41         #ifdef _HAVE_PETSC_
    42         this->matrix=NewMat(M,N);
    43         #else
    44         this->matrix=new SeqMat(M,N);
    45         #endif
    46         #ifdef _HAVE_ADOLC_
    47         this->amatrix=xNew<IssmDouble>(M*N);
    48         #endif
    49 }
    50 /*}}}*/
    51 /*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity){{{*/
    52 Matrix::Matrix(int M,int N,IssmDouble sparsity){
    53 
    54         #ifdef _HAVE_PETSC_
    55         this->matrix=NewMat(M,N,sparsity);
    56         #else
    57         this->matrix=new SeqMat(M,N,sparsity);
    58         #endif
    59         #ifdef _HAVE_ADOLC_
    60         this->amatrix=xNew<IssmDouble>(M*N);
    61         #endif
    62 }
    63 /*}}}*/
    64 /*FUNCTION Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/
    65 Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){
    66 
    67         #ifdef _HAVE_PETSC_
    68         int     i;
    69 
    70 
    71         int* idxm=xNew<int>(M);
    72         int* idxn=xNew<int>(N);
    73         for(i=0;i<M;i++)idxm[i]=i;
    74         for(i=0;i<N;i++)idxn[i]=i;
    75 
    76         this->matrix=NewMat(M,N,sparsity);
    77         MatSetValues(this->matrix,M,idxm,N,idxn,serial_mat,INSERT_VALUES);
    78         MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY);
    79         MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY);
    80 
    81         xDelete<int>(idxm);
    82         xDelete<int>(idxn);
    83         #else
    84         this->matrix=new SeqMat(serial_mat,M,N,sparsity);
    85         #endif
    86         #ifdef _HAVE_ADOLC_
    87         this->amatrix=xNew<IssmDouble>(M*N);
    88         #endif
    89 }
    90 /*}}}*/
    91 /*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/
    92 Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){
    93        
    94         #ifdef _HAVE_PETSC_
    95         this->matrix=NewMat(M,N,connectivity,numberofdofspernode);
    96         #else
    97         this->matrix=new SeqMat(M,N,connectivity,numberofdofspernode);
    98         #endif
    99         #ifdef _HAVE_ADOLC_
    100         this->amatrix=xNew<IssmDouble>(M*N);
    101         #endif
     35       
     36}
     37/*}}}*/
     38/*FUNCTION Matrix::Matrix(int M,int N,int type){{{*/
     39Matrix::Matrix(int M,int N,int in_type){
     40
     41        pmatrix=NULL;
     42        smatrix=NULL;
     43        type=in_type;
     44
     45        if(type==PetscMatType){
     46                #ifdef _HAVE_PETSC_
     47                this->pmatrix=new PetscMat(M,N);
     48                #else
     49                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     50                #endif
     51        }
     52        else if(type==SeqMatType){
     53                this->smatrix=new SeqMat(M,N);
     54        }
     55        else _error2_("Matrix type: " << type << " not supported yet!");
     56
     57}
     58/*}}}*/
     59/*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity,int type){{{*/
     60Matrix::Matrix(int M,int N,IssmDouble sparsity,int in_type){
     61
     62        pmatrix=NULL;
     63        smatrix=NULL;
     64        type=in_type;
     65
     66        if(type==PetscMatType){
     67                #ifdef _HAVE_PETSC_
     68                this->pmatrix=new PetscMat(M,N,sparsity);
     69                #else
     70                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     71                #endif
     72        }
     73        else if(type==SeqMatType){
     74                this->smatrix=new SeqMat(M,N,sparsity);
     75        }
     76        else _error2_("Matrix type: " << type << " not supported yet!");
     77}
     78/*}}}*/
     79/*FUNCTION Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int type){{{*/
     80Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int in_type){
     81
     82        pmatrix=NULL;
     83        smatrix=NULL;
     84        type=in_type;
     85
     86        if(type==PetscMatType){
     87                #ifdef _HAVE_PETSC_
     88                this->pmatrix=new PetscMat(serial_mat,M,N,sparsity);
     89                #else
     90                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     91                #endif
     92        }
     93        else if(type==SeqMatType){
     94                this->smatrix=new SeqMat(serial_mat,M,N,sparsity);
     95        }
     96        else _error2_("Matrix type: " << type << " not supported yet!");
     97       
     98}
     99/*}}}*/
     100/*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode,int type){{{*/
     101Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode,int in_type){
     102
     103        pmatrix=NULL;
     104        smatrix=NULL;
     105        type=in_type;
     106
     107        if(type==PetscMatType){
     108                #ifdef _HAVE_PETSC_
     109                this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode);
     110                #else
     111                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     112                #endif
     113        }
     114        else if(type==SeqMatType){
     115                this->smatrix=new SeqMat(M,N,connectivity,numberofdofspernode);
     116        }
     117        else _error2_("Matrix type: " << type << " not supported yet!");
     118       
    102119}
    103120/*}}}*/
    … …  
    105122Matrix::~Matrix(){
    106123
    107         #ifdef _HAVE_PETSC_
    108         MatFree(&this->matrix);
    109         #else
    110         delete this->matrix;
    111         #endif
    112         #ifdef _HAVE_ADOLC_
    113         xDelete<IssmDouble>(this->amatrix);
    114         #endif
     124        if(type==PetscMatType){
     125                #ifdef _HAVE_PETSC_
     126                delete this->pmatrix;
     127                #else
     128                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     129                #endif
     130        }
     131        else if(type==SeqMatType){
     132                delete this->smatrix;
     133        }
     134        else _error2_("Matrix type: " << type << " not supported yet!");
     135
    115136}
    116137/*}}}*/
    … …  
    120141void Matrix::Echo(void){
    121142
    122         int i,j;
    123 
    124         #ifdef _HAVE_PETSC_
    125         MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD);
    126         #else
    127         this->matrix->Echo();
    128         #endif
    129 
    130         #ifdef _HAVE_ADOLC_
    131         /*Not sure about that one. Should we use the overloaded operator >>?*/
    132         _printString_("ADOLC Matrix equivalent:" );
    133 //      for(i=0;i<M;i++){
    134 //              for(j=0;j<N;j++){
    135 //                      _printString_(*(amatrix+N*i+j) << " ");
    136 //              }
    137 //              _printLine_("");
    138 //      }
    139         #endif
     143        if(type==PetscMatType){
     144                #ifdef _HAVE_PETSC_
     145                this->pmatrix->Echo();
     146                #else
     147                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     148                #endif
     149        }
     150        else if(type==SeqMatType){
     151                this->smatrix->Echo();
     152        }
     153        else _error2_("Matrix type: " << type << " not supported yet!");
     154
    140155}
    141156/*}}}*/
    142157/*FUNCTION Matrix::Assemble{{{*/
    143158void Matrix::Assemble(void){
    144         #ifdef _HAVE_PETSC_
    145                 _assert_(this->matrix);
    146                 MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY);
    147                 MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY);
    148                 #if _PETSC_MAJOR_ == 2
    149                         MatCompress(this->matrix);
    150                 #endif
    151         #else
    152                 this->matrix->Assemble();
    153         #endif
    154 
     159
     160        if(type==PetscMatType){
     161                #ifdef _HAVE_PETSC_
     162                this->pmatrix->Assemble();
     163                #else
     164                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     165                #endif
     166        }
     167        else if(type==SeqMatType){
     168                this->smatrix->Assemble();
     169        }
     170        else _error2_("Matrix type: " << type << " not supported yet!");
    155171}
    156172/*}}}*/
    … …  
    159175       
    160176        IssmDouble norm=0;
    161         #ifdef _HAVE_PETSC_
    162                 _assert_(this->matrix);
    163                 MatNorm(this->matrix,ISSMToPetscNormMode(norm_type),&norm);
    164         #else
    165                 norm=this->matrix->Norm(norm_type);
    166         #endif
     177
     178        if(type==PetscMatType){
     179                #ifdef _HAVE_PETSC_
     180                norm=this->pmatrix->Norm(norm_type);
     181                #else
     182                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     183                #endif
     184        }
     185        else if(type==SeqMatType){
     186                norm=this->smatrix->Norm(norm_type);
     187        }
     188        else _error2_("Matrix type: " << type << " not supported yet!");
     189
    167190        return norm;
    168191}
    … …  
    170193/*FUNCTION Matrix::GetSize{{{*/
    171194void Matrix::GetSize(int* pM,int* pN){
    172        
    173         #ifdef _HAVE_PETSC_
    174                 _assert_(this->matrix);
    175                 MatGetSize(this->matrix,pM,pN);
    176         #else
    177                 this->matrix->GetSize(pM,pN);
    178         #endif
     195
     196        if(type==PetscMatType){
     197                #ifdef _HAVE_PETSC_
     198                this->pmatrix->GetSize(pM,pN);
     199                #else
     200                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     201                #endif
     202        }
     203        else if(type==SeqMatType){
     204                this->smatrix->GetSize(pM,pN);
     205        }
     206        else _error2_("Matrix type: " << type << " not supported yet!");
     207       
    179208}
    180209/*}}}*/
    … …  
    182211void Matrix::GetLocalSize(int* pM,int* pN){
    183212       
    184         #ifdef _HAVE_PETSC_
    185                 _assert_(this->matrix);
    186                 MatGetLocalSize(this->matrix,pM,pN);
    187         #else
    188                 this->matrix->GetLocalSize(pM,pN);
    189         #endif
     213        if(type==PetscMatType){
     214                #ifdef _HAVE_PETSC_
     215                this->pmatrix->GetLocalSize(pM,pN);
     216                #else
     217                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     218                #endif
     219        }
     220        else if(type==SeqMatType){
     221                this->smatrix->GetLocalSize(pM,pN);
     222        }
     223        else _error2_("Matrix type: " << type << " not supported yet!");
     224
    190225}
    191226/*}}}*/
    … …  
    193228void Matrix::MatMult(Vector* X,Vector* AX){
    194229
    195         #ifdef _HAVE_PETSC_
    196                 _assert_(this->matrix);
    197                 _assert_(X->vector);
    198                 MatMultPatch(this->matrix,X->vector,AX->vector);
    199         #else
    200                 this->matrix->MatMult(X->vector,AX->vector);
    201         #endif
     230        if(type==PetscMatType){
     231                #ifdef _HAVE_PETSC_
     232                this->pmatrix->MatMult(X->pvector,AX->pvector);
     233                #else
     234                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     235                #endif
     236        }
     237        else if(type==SeqMatType){
     238                this->smatrix->MatMult(X->svector,AX->svector);
     239        }
     240        else _error2_("Matrix type: " << type << " not supported yet!");
     241
    202242}
    203243/*}}}*/
    … …  
    209249        output=new Matrix();
    210250
    211         #ifdef _HAVE_PETSC_
    212                 _assert_(this->matrix);
    213                 MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix);
    214         #else
    215                 output->matrix=this->matrix->Duplicate();
    216         #endif
     251        if(type==PetscMatType){
     252                #ifdef _HAVE_PETSC_
     253                output->pmatrix=this->pmatrix->Duplicate();
     254                #else
     255                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     256                #endif
     257        }
     258        else if(type==SeqMatType){
     259                output->smatrix=this->smatrix->Duplicate();
     260        }
     261        else _error2_("Matrix type: " << type << " not supported yet!");
    217262       
    218263        return output;
    … …  
    223268
    224269        IssmDouble* output=NULL;
    225 
    226         #ifdef _HAVE_PETSC_
    227                 MatToSerial(&output,this->matrix);
    228         #else
    229                 output=this->matrix->ToSerial();
    230         #endif
     270       
     271        if(type==PetscMatType){
     272                #ifdef _HAVE_PETSC_
     273                output=this->pmatrix->ToSerial();
     274                #else
     275                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     276                #endif
     277        }
     278        else if(type==SeqMatType){
     279                output=this->smatrix->ToSerial();
     280        }
     281        else _error2_("Matrix type: " << type << " not supported yet!");
     282
     283
    231284        return output;
    232285}
    … …  
    234287/*FUNCTION Matrix::SetValues{{{*/
    235288void Matrix::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){
    236 
    237         #ifdef _HAVE_PETSC_
    238                 MatSetValues(this->matrix,m,idxm,n,idxn,values,ISSMToPetscInsertMode(mode));
    239         #else
    240                 this->matrix->SetValues(m,idxm,n,idxn,values,mode);
    241         #endif
     289               
     290        if(type==PetscMatType){
     291                #ifdef _HAVE_PETSC_
     292                this->pmatrix->SetValues(m,idxm,n,idxn,values,mode);
     293                #else
     294                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     295                #endif
     296        }
     297        else if(type==SeqMatType){
     298                this->smatrix->SetValues(m,idxm,n,idxn,values,mode);
     299        }
     300        else _error2_("Matrix type: " << type << " not supported yet!");
    242301}
    243302/*}}}*/
    … …  
    245304void Matrix::Convert(MatrixType type){
    246305
    247         #ifdef _HAVE_PETSC_
    248                 MatConvert(this->matrix,ISSMToPetscMatrixType(type),MAT_REUSE_MATRIX,&this->matrix);
    249         #else
    250                 this->matrix->Convert(type);
    251         #endif
    252 }
    253 /*}}}*/
     306        if(type==PetscMatType){
     307                #ifdef _HAVE_PETSC_
     308                this->pmatrix->Convert(type);
     309                #else
     310                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     311                #endif
     312        }
     313        else if(type==SeqMatType){
     314                this->smatrix->Convert(type);
     315        }
     316        else _error2_("Matrix type: " << type << " not supported yet!");
     317
     318}
     319/*}}}*/

  • issm/trunk-jpl/src/c/classes/matrix/Matrix.h

    r12821 r12850  
    1818
    1919class Vector;
     20/*}}}*/
    2021
    21 /*}}}*/
     22enum matrixtype { PetscMatType, SeqMatType };
    2223
    2324class Matrix{
    2425
    2526        public:
    26        
     27
    2728                #ifdef _HAVE_PETSC_
    28                 Mat matrix;
    29                 #else
    30                 SeqMat* matrix;
     29                PetscMat* pmatrix;
    3130                #endif
    32                 #ifdef _HAVE_ADOLC_
    33                 IssmDouble* amatrix;
    34                 #endif
     31                SeqMat* smatrix;
     32                int     type;
    3533
    3634                /*Matrix constructors, destructors {{{*/
    3735                Matrix();
    38                 Matrix(int M,int N);
    39                 Matrix(int M,int N,IssmDouble sparsity);
    40                 Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity);
    41                 Matrix(int M,int N,int connectivity,int numberofdofspernode);
     36                Matrix(int M,int N,int type=PetscMatType);
     37                Matrix(int M,int N,IssmDouble sparsity,int type=PetscMatType);
     38                Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int type=PetscMatType);
     39                Matrix(int M,int N,int connectivity,int numberofdofspernode,int type=PetscMatType);
    4240                ~Matrix();
    4341                /*}}}*/

  • issm/trunk-jpl/src/c/classes/matrix/Vector.cpp

    r12832 r12850  
    2525Vector::Vector(){
    2626
    27         #ifdef _HAVE_PETSC_
    28         this->vector=NULL;
    29         #else
    30         this->vector=NULL;
     27        this->pvector=NULL;
     28        this->svector=NULL;
     29       
     30        type=PetscVecType; //default
     31        #ifndef _HAVE_PETSC_
     32        type=SeqVecType;
    3133        #endif
    32         #ifdef _HAVE_ADOLC_
    33         this->avector=NULL;
    34         #endif
    35 }
    36 /*}}}*/
    37 /*FUNCTION Vector::Vector(int M,bool fromlocalsize){{{*/
    38 Vector::Vector(int pM,bool fromlocalsize){
    39 
    40         #ifdef _HAVE_PETSC_
    41         this->vector=NewVec(pM,fromlocalsize);
    42         #else
    43         this->vector=new SeqVec(pM,fromlocalsize);
    44         #endif
    45         #ifdef _HAVE_ADOLC_
    46         this->avector=xNew<IssmDouble>(pM);
    47         #endif
    48 }
    49 /*}}}*/
    50 /*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M){{{*/
    51 Vector::Vector(IssmDouble* serial_vec,int M){
    52 
    53         #ifdef _HAVE_PETSC_
    54                 int* idxm=xNew<int>(M);
    55                 for(int i=0;i<M;i++) idxm[i]=i;
    56 
    57                 this->vector=NewVec(M);
    58                 VecSetValues(this->vector,M,idxm,serial_vec,INSERT_VALUES);
    59                 VecAssemblyBegin(this->vector);
    60                 VecAssemblyEnd(this->vector);
    61 
    62                 xDelete<int>(idxm);
    63         #else
    64                 this->vector=new SeqVec(serial_vec,M);
    65         #endif
    66         #ifdef _HAVE_ADOLC_
    67                 this->avector=xNew<IssmDouble>(M);
    68         #endif
    69 }
    70 /*}}}*/
    71 #ifdef _HAVE_PETSC_
    72 /*FUNCTION Vector::Vector(Vec petsc_vec){{{*/
    73 Vector::Vector(Vec petsc_vec){
    74 
    75         if(petsc_vec==NULL){
    76                 this->vector=NewVec(0);
    77         }
    78         else{
    79                 /*copy vector*/
    80                 VecDuplicate(petsc_vec,&this->vector);
    81                 VecCopy(petsc_vec,this->vector);
    82         }
    83 
    84 }
    85 /*}}}*/
    86 #endif
    87 #if defined(_HAVE_GSL_) && !defined(_HAVE_PETSC_)
    88 /*FUNCTION Vector::Vector(SeqVec* seq_vec){{{*/
    89 Vector::Vector(SeqVec*  seq_vec){
    90 
    91         if(seq_vec==NULL){
    92                 this->vector=NULL;
    93         }
    94         else{
    95                 /*copy vector*/
    96                 this->vector=seq_vec->Duplicate();
    97         }
    98 }
    99 /*}}}*/
    100 #endif
    101 
    102                 /*FUNCTION Vector::~Vector(){{{*/
     34
     35}
     36/*}}}*/
     37/*FUNCTION Vector::Vector(int M,bool fromlocalsize,int type){{{*/
     38Vector::Vector(int M,bool fromlocalsize,int in_type){
     39       
     40        pvector=NULL;
     41        svector=NULL;
     42        type=in_type;
     43
     44
     45        if(type==PetscVecType){
     46                #ifdef _HAVE_PETSC_
     47                this->pvector=new PetscVec(M,fromlocalsize);
     48                #else
     49                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     50                #endif
     51        }
     52        else if(type==SeqVecType){
     53                this->svector=new SeqVec(M,fromlocalsize);
     54        }
     55        else _error2_("Vector type: " << type << " not supported yet!");
     56
     57}
     58/*}}}*/
     59/*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M,int type){{{*/
     60Vector::Vector(IssmDouble* serial_vec,int M,int in_type){
     61
     62        pvector=NULL;
     63        svector=NULL;
     64        type=in_type;
     65
     66        if(type==PetscVecType){
     67                #ifdef _HAVE_PETSC_
     68                this->pvector=new PetscVec(serial_vec,M);
     69                #else
     70                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     71                #endif
     72        }
     73        else if(type==SeqVecType){
     74                this->svector=new SeqVec(serial_vec,M);
     75        }
     76        else _error2_("Vector type: " << type << " not supported yet!");
     77
     78}
     79/*}}}*/
     80/*FUNCTION Vector::~Vector(){{{*/
    10381Vector::~Vector(){
    10482
    105         #ifdef _HAVE_PETSC_
    106         VecFree(&this->vector);
    107         #else
    108         delete this->vector;
    109         #endif
    110         #ifdef _HAVE_ADOLC_
    111         xDelete<IssmDouble>(this->avector);
    112         #endif
     83        if(type==PetscVecType){
     84                #ifdef _HAVE_PETSC_
     85                delete this->pvector;
     86                #else
     87                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     88                #endif
     89        }
     90        else if(type==SeqVecType){
     91                delete this->svector;
     92        }
     93        else _error2_("Vector type: " << type << " not supported yet!");
    11394}
    11495/*}}}*/
    … …  
    11899void Vector::Echo(void){
    119100
    120         int i;
    121 
    122         #ifdef _HAVE_PETSC_
    123         _assert_(this->vector);
    124         VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
    125         #else
    126         this->vector->Echo();
    127         #endif
    128 
    129         #ifdef _HAVE_ADOLC_
    130         /*do nothing for now: */
    131         #endif
     101        if(type==PetscVecType){
     102                #ifdef _HAVE_PETSC_
     103                this->pvector->Echo();
     104                #else
     105                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     106                #endif
     107        }
     108        else if(type==SeqVecType){
     109                this->svector->Echo();
     110        }
     111        else _error2_("Vector type: " << type << " not supported yet!");
     112
    132113}
    133114/*}}}*/
    134115/*FUNCTION Vector::Assemble{{{*/
    135116void Vector::Assemble(void){
    136                
    137         #ifdef _HAVE_PETSC_
    138                 _assert_(this->vector);
    139                 VecAssemblyBegin(this->vector);
    140                 VecAssemblyEnd(this->vector);
    141         #else
    142                 this->vector->Assemble();
    143         #endif
     117
     118        if(type==PetscVecType){
     119                #ifdef _HAVE_PETSC_
     120                this->pvector->Assemble();
     121                #else
     122                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     123                #endif
     124        }
     125        else if(type==SeqVecType){
     126                this->svector->Assemble();
     127        }
     128        else _error2_("Vector type: " << type << " not supported yet!");
    144129
    145130}
    … …  
    148133void Vector::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){
    149134               
     135
     136        if(type==PetscVecType){
     137                #ifdef _HAVE_PETSC_
     138                this->pvector->SetValues(ssize,list,values,mode);
     139                #else
     140                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     141                #endif
     142        }
     143        else if(type==SeqVecType){
     144                this->svector->SetValues(ssize,list,values,mode);
     145        }
     146        else _error2_("Vector type: " << type << " not supported yet!");
     147
    150148               
    151         #ifdef _HAVE_PETSC_
    152                 _assert_(this->vector);
    153                 VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode));
    154         #else
    155                 this->vector->SetValues(ssize,list,values,mode);
    156         #endif
    157149
    158150}
    … …  
    160152/*FUNCTION Vector::SetValue{{{*/
    161153void Vector::SetValue(int dof, IssmDouble value, InsMode mode){
    162                
    163         #ifdef _HAVE_PETSC_
    164                 _assert_(this->vector);
    165                 VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode));
    166         #else
    167                 this->vector->SetValue(dof,value,mode);
    168         #endif
     154       
     155        if(type==PetscVecType){
     156                #ifdef _HAVE_PETSC_
     157                this->pvector->SetValue(dof,value,mode);
     158                #else
     159                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     160                #endif
     161        }
     162        else if(type==SeqVecType){
     163                this->svector->SetValue(dof,value,mode);
     164        }
     165        else _error2_("Vector type: " << type << " not supported yet!");
    169166
    170167}
    … …  
    173170void Vector::GetValue(IssmDouble* pvalue,int dof){
    174171               
    175         #ifdef _HAVE_PETSC_
    176                 _assert_(this->vector);
    177                 VecGetValues(this->vector,1,&dof,pvalue);
    178         #else
    179         this->vector->GetValue(pvalue,dof);
    180         #endif
     172
     173        if(type==PetscVecType){
     174                #ifdef _HAVE_PETSC_
     175                this->pvector->GetValue(pvalue,dof);
     176                #else
     177                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     178                #endif
     179        }
     180        else if(type==SeqVecType){
     181                this->svector->GetValue(pvalue,dof);
     182        }
     183        else _error2_("Vector type: " << type << " not supported yet!");
     184
    181185}
    182186/*}}}*/
    183187/*FUNCTION Vector::GetSize{{{*/
    184188void Vector::GetSize(int* pM){
    185                
    186         #ifdef _HAVE_PETSC_
    187                 _assert_(this->vector);
    188                 VecGetSize(this->vector,pM);
    189         #else
    190                 this->vector->GetSize(pM);
    191         #endif
     189       
     190        if(type==PetscVecType){
     191                #ifdef _HAVE_PETSC_
     192                this->pvector->GetSize(pM);
     193                #else
     194                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     195                #endif
     196        }
     197        else if(type==SeqVecType){
     198                this->svector->GetSize(pM);
     199        }
     200        else _error2_("Vector type: " << type << " not supported yet!");
    192201
    193202}
    … …  
    198207        int M;
    199208
    200         _assert_(this->vector);
    201209        this->GetSize(&M);
    202210
    … …  
    210218void Vector::GetLocalSize(int* pM){
    211219               
    212         #ifdef _HAVE_PETSC_
    213                 _assert_(this->vector);
    214                 VecGetLocalSize(this->vector,pM);
    215         #else
    216                 this->vector->GetLocalSize(pM);
    217         #endif
     220
     221        if(type==PetscVecType){
     222                #ifdef _HAVE_PETSC_
     223                this->pvector->GetLocalSize(pM);
     224                #else
     225                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     226                #endif
     227        }
     228        else if(type==SeqVecType){
     229                this->svector->GetLocalSize(pM);
     230        }
     231        else _error2_("Vector type: " << type << " not supported yet!");
    218232
    219233}
    … …  
    223237       
    224238        Vector* output=NULL;
    225         #ifdef _HAVE_PETSC_
    226                 _assert_(this->vector);
    227                 Vec vec_output=NULL;
    228                 VecDuplicate(this->vector,&vec_output);
    229                 output=new Vector(vec_output);
    230                 VecFree(&vec_output);
    231         #else
     239
     240
     241        if(type==PetscVecType){
     242                #ifdef _HAVE_PETSC_
    232243                output=new Vector();
    233                 output->vector=this->vector->Duplicate();
    234         #endif
     244                output->pvector=this->pvector->Duplicate();
     245                #else
     246                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     247                #endif
     248        }
     249        else if(type==SeqVecType){
     250                output=new Vector();
     251                output->svector=this->svector->Duplicate();
     252        }
     253        else _error2_("Vector type: " << type << " not supported yet!");
     254
    235255        return output;
    236256
    … …  
    240260void Vector::Set(IssmDouble value){
    241261       
    242         #ifdef _HAVE_PETSC_
    243                 _assert_(this->vector);
    244                 VecSet(this->vector,value);
    245         #else
    246                 this->vector->Set(value);
    247         #endif
     262       
     263        if(type==PetscVecType){
     264                #ifdef _HAVE_PETSC_
     265                this->pvector->Set(value);
     266                #else
     267                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     268                #endif
     269        }
     270        else if(type==SeqVecType){
     271                this->svector->Set(value);
     272        }
     273        else _error2_("Vector type: " << type << " not supported yet!");
    248274
    249275}
    … …  
    252278void Vector::AXPY(Vector* X, IssmDouble a){
    253279       
    254         #ifdef _HAVE_PETSC_
    255                 _assert_(this->vector);
    256                 VecAXPY(this->vector,a,X->vector);
    257         #else
    258                 this->vector->AXPY(X->vector,a);
    259         #endif
     280
     281        if(type==PetscVecType){
     282                #ifdef _HAVE_PETSC_
     283                this->pvector->AXPY(X->pvector,a);
     284                #else
     285                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     286                #endif
     287        }
     288        else if(type==SeqVecType){
     289                this->svector->AXPY(X->svector,a);
     290        }
     291        else _error2_("Vector type: " << type << " not supported yet!");
     292
    260293}
    261294/*}}}*/
    … …  
    263296void Vector::AYPX(Vector* X, IssmDouble a){
    264297       
    265         #ifdef _HAVE_PETSC_
    266                 _assert_(this->vector);
    267                 VecAYPX(this->vector,a,X->vector);
    268         #else
    269                 this->vector->AYPX(X->vector,a);
    270         #endif
     298
     299        if(type==PetscVecType){
     300                #ifdef _HAVE_PETSC_
     301                this->pvector->AYPX(X->pvector,a);
     302                #else
     303                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     304                #endif
     305        }
     306        else if(type==SeqVecType){
     307                this->svector->AYPX(X->svector,a);
     308        }
     309        else _error2_("Vector type: " << type << " not supported yet!");
     310
    271311
    272312}
    … …  
    277317        IssmDouble* vec_serial=NULL;
    278318
    279         #ifdef _HAVE_PETSC_
    280                 VecToMPISerial(&vec_serial, this->vector);
    281         #else
    282                 vec_serial=this->vector->ToMPISerial();
    283         #endif
     319        if(type==PetscVecType){
     320                #ifdef _HAVE_PETSC_
     321                vec_serial=this->pvector->ToMPISerial();
     322                #else
     323                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     324                #endif
     325        }
     326        else if(type==SeqVecType){
     327                vec_serial=this->svector->ToMPISerial();
     328        }
     329        else _error2_("Vector type: " << type << " not supported yet!");
    284330
    285331        return vec_serial;
    … …  
    290336void Vector::Copy(Vector* to){
    291337
    292         #ifdef _HAVE_PETSC_
    293                 if(this->vector) VecCopy(this->vector,to->vector);
    294         #else
    295                 this->vector->Copy(to->vector);
    296         #endif
     338       
     339        if(type==PetscVecType){
     340                #ifdef _HAVE_PETSC_
     341                this->pvector->Copy(to->pvector);
     342                #else
     343                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     344                #endif
     345        }
     346        else if(type==SeqVecType){
     347                this->svector->Copy(to->svector);
     348        }
     349        else _error2_("Vector type: " << type << " not supported yet!");
     350
    297351
    298352}
    … …  
    302356       
    303357        IssmDouble norm=0;
    304         #ifdef _HAVE_PETSC_
    305                 _assert_(this->vector);
    306                 VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm);
    307         #else
    308                 norm=this->vector->Norm(norm_type);
    309         #endif
     358
     359        if(type==PetscVecType){
     360                #ifdef _HAVE_PETSC_
     361                norm=this->pvector->Norm(norm_type);
     362                #else
     363                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     364                #endif
     365        }
     366        else if(type==SeqVecType){
     367                norm=this->svector->Norm(norm_type);
     368        }
     369        else _error2_("Vector type: " << type << " not supported yet!");
     370
    310371        return norm;
    311372}
    … …  
    314375void Vector::Scale(IssmDouble scale_factor){
    315376       
    316         #ifdef _HAVE_PETSC_
    317                 _assert_(this->vector);
    318                 VecScale(this->vector,scale_factor);
    319         #else
    320                 this->vector->Scale(scale_factor);
    321         #endif
     377
     378        if(type==PetscVecType){
     379                #ifdef _HAVE_PETSC_
     380                this->pvector->Scale(scale_factor);
     381                #else
     382                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     383                #endif
     384        }
     385        else if(type==SeqVecType){
     386                this->svector->Scale(scale_factor);
     387        }
     388        else _error2_("Vector type: " << type << " not supported yet!");
     389
    322390}
    323391/*}}}*/
    … …  
    326394
    327395        IssmDouble dot;
    328         #ifdef _HAVE_PETSC_
    329                 _assert_(this->vector);
    330                 VecDot(this->vector,vector->vector,&dot);
    331         #else
    332                 dot=this->vector->Dot(vector->vector);
    333         #endif
     396       
     397        if(type==PetscVecType){
     398                #ifdef _HAVE_PETSC_
     399                dot=this->pvector->Dot(vector->pvector);
     400                #else
     401                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     402                #endif
     403        }
     404        else if(type==SeqVecType){
     405                dot=this->svector->Dot(vector->svector);
     406        }
     407        else _error2_("Vector type: " << type << " not supported yet!");
     408
    334409        return dot;
    335410}
    … …  
    338413void Vector::PointwiseDivide(Vector* x,Vector* y){
    339414
    340         #ifdef _HAVE_PETSC_
    341                 _assert_(this->vector);
    342                 VecPointwiseDivide(this->vector,x->vector,y->vector);
    343         #else
    344                 this->vector->PointwiseDivide(x->vector,y->vector);
    345         #endif
    346 }
    347 /*}}}*/
     415
     416        if(type==PetscVecType){
     417                #ifdef _HAVE_PETSC_
     418                this->pvector->PointwiseDivide(x->pvector,y->pvector);
     419                #else
     420                _error2_("Petsc matrix format not usable, as Petsc has not been compiled!");
     421                #endif
     422        }
     423        else if(type==SeqVecType){
     424                this->svector->PointwiseDivide(x->svector,y->svector);
     425        }
     426        else _error2_("Vector type: " << type << " not supported yet!");
     427
     428}
     429/*}}}*/

  • issm/trunk-jpl/src/c/classes/matrix/Vector.h

    r12821 r12850  
    1919/*}}}*/
    2020
     21enum vectortype { PetscVecType, SeqVecType };
     22
    2123class Vector{
    2224
    2325        public:
     26
     27                #ifdef _HAVE_PETSC_
     28                PetscVec* pvector;
     29                #endif
     30                SeqVec* svector;
     31                int     type;
    2432       
    25                 #ifdef _HAVE_PETSC_
    26                         Vec vector;
    27                 #else
    28                         SeqVec* vector;
    29                 #endif
    30                 #ifdef _HAVE_ADOLC_
    31                         IssmDouble* avector;
    32                 #endif
    3333
    3434                /*Vector constructors, destructors {{{*/
    3535                Vector();
    36                 Vector(int M,bool fromlocalsize=false);
    37                 Vector(IssmDouble* serial_vec,int pM);
    38                 #ifdef _HAVE_PETSC_
    39                 Vector(Vec petsc_vec);
    40                 #endif
    41                 #if defined(_HAVE_GSL_) && !defined(_HAVE_PETSC_)
    42                 Vector(SeqVec*  seq_vec);
    43                 #endif
     36                Vector(int M,bool fromlocalsize=false,int type=PetscVecType);
     37                Vector(IssmDouble* serial_vec,int pM,int type=PetscVecType);
    4438                ~Vector();
    4539                /*}}}*/
    4640                /*Vector specific routines {{{*/
    47                 void Echo(void);
     41                void    Echo(void);
    4842                void    AXPY(Vector *X, IssmDouble a);
    4943                void    AYPX(Vector *X, IssmDouble a);

  • issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToMatrix.cpp

    r12835 r12850  
    3131
    3232        #ifdef _HAVE_PETSC_
    33         MatlabMatrixToPetscMatrix(&matrix->matrix,NULL,NULL,mxmatrix);
     33        matrix->pmatrix=MatlabMatrixToPetscMat(mxmatrix);
    3434        #else
    35         matrix->matrix=MatlabMatrixToSeqMat(mxmatrix);
     35        matrix->smatrix=MatlabMatrixToSeqMat(mxmatrix);
    3636        #endif
    3737       

  • issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMat.cpp

    r12849 r12850  
    1818/*Matlab includes: */
    1919#include "mex.h"
     20#include "matlabio.h"
    2021
    21 int MatlabMatrixToPetscMatrix(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){
     22PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix){
     23
     24        int dummy;
     25        PetscMat* matrix=new PetscMat();
     26
     27        MatlabMatrixToPetscMat(&matrix->matrix, &dummy, &dummy, mxmatrix);
     28
     29        return matrix;
     30}
     31int MatlabMatrixToPetscMat(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){
    2232
    2333        int rows, cols;

  • issm/trunk-jpl/src/c/matlab/io/MatlabVectorToPetscVec.cpp

    r12849 r12850  
    1919
    2020#include "../../shared/shared.h"
     21#include "matlabio.h"
    2122
    22 int MatlabVectorToPetscVector(Vec* pvector,int* pvector_rows,const mxArray* mxvector){
     23PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector){
     24
     25        int dummy;
     26        PetscVec* vector=new PetscVec();
     27       
     28        MatlabVectorToPetscVec(&vector->vector,&dummy, mxvector);
     29
     30        return vector;
     31}
     32
     33int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector){
    2334
    2435        int rows, cols;

  • issm/trunk-jpl/src/c/matlab/io/MatlabVectorToVector.cpp

    r12835 r12850  
    3131
    3232        #ifdef _HAVE_PETSC_
    33         MatlabVectorToPetscVector(&vector->vector,&dummy,mxvector);
     33        vector->pvector=MatlabVectorToPetscVec(mxvector);
    3434        #else
    35         vector->vector=MatlabVectorToSeqVec(mxvector);
     35        vector->svector=MatlabVectorToSeqVec(mxvector);
    3636        #endif
    3737       

  • issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp

    r12434 r12850  
    227227        if(matrix){
    228228               
    229                 #ifdef _HAVE_PETSC_
    230                 PetscMatrixToDoubleMatrix(&matrix_ptr,&rows,&cols,matrix->matrix);
    231                 #else
    232                 matrix_ptr=matrix->matrix->ToSerial();
    233                 matrix->matrix->GetSize(&rows,&cols);
    234                 #endif
     229                matrix_ptr=matrix->ToSerial();
     230                matrix->GetSize(&rows,&cols);
    235231
    236232                /*Now transpose the matrix and allocate with Matlab's memory manager: */
    … …  
    269265        if(vector){
    270266                /*call toolkit routine: */
    271                 #ifdef _HAVE_PETSC_
    272                 PetscVectorToDoubleVector(&vector_ptr,&rows,vector->vector);
    273                 #else
    274                 vector_ptr=vector->vector->ToMPISerial();
    275                 vector->vector->GetSize(&rows);
    276                 #endif
     267                vector_ptr=vector->ToMPISerial();
     268                vector->GetSize(&rows);
    277269               
    278270                /*now create the matlab vector with Matlab's memory manager */

  • issm/trunk-jpl/src/c/matlab/io/matlabio.h

    r12835 r12850  
    7979/*Matlab to Petsc routines: */
    8080#ifdef _HAVE_PETSC_
    81 int MatlabMatrixToPetscMatrix(Mat* matrix,int* prows,int* pcols, const mxArray* mxmatrix);
    82 int MatlabVectorToPetscVector(Vec* pvector,int* pvector_rows,const mxArray* mxvector);
     81int MatlabMatrixToPetscMat(Mat* matrix,int* prows,int* pcols, const mxArray* mxmatrix);
     82PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix);
     83int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector);
     84PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector);
    8385#endif
    8486

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

    r12493 r12850  
    2626        _assert_(pf);
    2727
    28         #ifdef _HAVE_PETSC_
    29                 Vec uf0_vector = NULL;
    30                 Vec df_vector  = NULL;
    31                 Vec uf_vector  = NULL;
    32                 if(uf0) uf0_vector = uf0->vector;
    33                 if(df)  df_vector  = df->vector;
     28        /*Initialize vector: */
     29        uf=new Vector();
    3430
    35                 /*In serial mode, the Petsc Options database has not been initialized properly: */
     31        /*According to matrix type, use specific solvers: */
     32        if(Kff->type==PetscMatType){
     33                PetscVec* uf0_vector = NULL;
     34                PetscVec* df_vector  = NULL;
     35                if(uf0) uf0_vector = uf0->pvector;
     36                if(df)  df_vector  = df->pvector;
    3637
    37                 SolverxPetsc(&uf_vector,Kff->matrix,pf->vector,uf0_vector,df_vector,parameters);
    38 
    39                 /*Create vector out of petsc vector: */
    40                 uf=new Vector(uf_vector);
    41 
    42                 /*Free ressources: */
    43                 VecFree(&uf_vector);
    44         #else
    45                 #ifdef _HAVE_GSL_
    46                 SeqVec* uf_vector=NULL;
    47 
    48                 SolverxGsl(&uf_vector,Kff->matrix,pf->vector);
    49 
    50                 /*Create vector out of SeqVec vector: */
    51                 uf=new Vector(uf_vector);
    52 
    53                 /*Free ressources: */
    54                 delete uf_vector;
    55 
    56                 #else
    57                         _error2_("GSL support not compiled in!");
    58                 #endif
    59         #endif
     38                SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters);
     39        }
     40        else if(Kff->type==SeqMatType){
     41                SolverxSeq(&uf->svector,Kff->smatrix,pf->svector);
     42        }
     43        else _error2_("Matrix type: " << Kff->type << " not supported yet!");
    6044
    6145        /*Assign output pointers: */

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

    r12832 r12850  
    1818
    1919#ifdef _HAVE_PETSC_
     20void    SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters);
    2021void    SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters);
    2122void  DofTypesToIndexSet(IS* pisv, IS* pisp, Vec df,int typeenum);
    2223#endif
    2324
    24 #ifdef _HAVE_GSL_
    25 void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
    26 void SolverxGsl(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n);
    27 #endif
     25void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf);
     26void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n);
    2827
    2928#endif  /* _SOLVERX_H */

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

    r12515 r12850  
    1414#endif
    1515
     16void    SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){
     17
     18        PetscVec* uf=new PetscVec();
     19        SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0->vector, df->vector, parameters);
     20        *puf=uf;
     21
     22}
    1623void    SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){
    1724

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

    r12849 r12850  
    1 /*!\file SolverxGsl
    2  * \brief Gsl implementation of solver
     1/*!\file SolverxSeq
     2 * \brief implementation of sequential solver using the GSL librarie
    33 */
    44
    … …  
    99#endif
    1010#include <cstring>
    11 #include <gsl/gsl_linalg.h>
    1211
    1312#include "./Solverx.h"
    … …  
    1615#include "../../io/io.h"
    1716
    18 void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/
     17#ifdef _HAVE_GSL_
     18#include <gsl/gsl_linalg.h>
     19#endif
    1920
     21void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/
     22
     23        #ifdef _HAVE_GSL_
    2024        /*Intermediary: */
    2125        int M,N,N2,s;
    … …  
    2933        if(M!=N)_error2_("Stiffness matrix should be square!");
    3034
    31         SolverxGsl(&x,Kff->matrix,pf->vector,N);
     35        SolverxSeq(&x,Kff->matrix,pf->vector,N);
    3236        uf=new SeqVec(x,N);
    3337
    3438        /*Assign output pointers:*/
    3539        *puf=uf;
     40
     41        #else
     42                _error2_("GSL support not compiled in!");
     43        #endif
     44
    3645}/*}}}*/
    37 void SolverxGsl(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/
     46void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/
    3847
    39         /*GSL Matrices and vectors: */
    40         int              s;
    41         gsl_matrix_view  a;
    42         gsl_vector_view  b;
    43         gsl_vector      *x = NULL;
    44         gsl_permutation *p = NULL;
    45 #ifdef _HAVE_ADOLC_
    46         // if we use Adol-C then the IssmDouble will be an adouble
    47         // and the calls to gsl_... will not work
    48         // and we should call  a suitable wrapped solve instead
    49         _error2_("SolverxGsl: should not be here with Adol-C");
    50 #else
    51         /*A will be modified by LU decomposition. Use copy*/
    52         IssmDouble* Acopy = xNew<IssmDouble>(n*n);
    53         xMemCpy<IssmDouble>(Acopy,A,n*n);
     48        #ifdef _HAVE_GSL_
     49                /*GSL Matrices and vectors: */
     50                int              s;
     51                gsl_matrix_view  a;
     52                gsl_vector_view  b;
     53                gsl_vector      *x = NULL;
     54                gsl_permutation *p = NULL;
     55                #ifdef _HAVE_ADOLC_
     56                        // if we use Adol-C then the IssmDouble will be an adouble
     57                        // and the calls to gsl_... will not work
     58                        // and we should call  a suitable wrapped solve instead
     59                        _error2_("SolverxSeq: should not be here with Adol-C");
     60                #else
     61                        /*A will be modified by LU decomposition. Use copy*/
     62                        IssmDouble* Acopy = xNew<IssmDouble>(n*n);
     63                        xMemCpy<IssmDouble>(Acopy,A,n*n);
    5464
    55         /*Initialize gsl matrices and vectors: */
    56         a = gsl_matrix_view_array (Acopy,n,n);
    57         b = gsl_vector_view_array (B,n);
    58         x = gsl_vector_alloc (n);
     65                        /*Initialize gsl matrices and vectors: */
     66                        a = gsl_matrix_view_array (Acopy,n,n);
     67                        b = gsl_vector_view_array (B,n);
     68                        x = gsl_vector_alloc (n);
    5969
    60         /*Run LU and solve: */
    61         p = gsl_permutation_alloc (n);
    62         gsl_linalg_LU_decomp (&a.matrix, p, &s);
    63         gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
     70                        /*Run LU and solve: */
     71                        p = gsl_permutation_alloc (n);
     72                        gsl_linalg_LU_decomp (&a.matrix, p, &s);
     73                        gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
    6474
    65         //printf ("x = \n");
    66         //gsl_vector_fprintf (stdout, x, "%g");
     75                        //printf ("x = \n");
     76                        //gsl_vector_fprintf (stdout, x, "%g");
    6777
    68         /*Copy result*/
    69         IssmDouble* X = xNew<IssmDouble>(n);
    70         memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble));
     78                        /*Copy result*/
     79                        IssmDouble* X = xNew<IssmDouble>(n);
     80                        memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble));
    7181
    72         /*Clean up and assign output pointer*/
    73         xDelete<IssmDouble>(Acopy);
    74         gsl_permutation_free(p);
    75         gsl_vector_free(x);
    76         *pX=X;
    77 #endif
     82                        /*Clean up and assign output pointer*/
     83                        xDelete<IssmDouble>(Acopy);
     84                        gsl_permutation_free(p);
     85                        gsl_vector_free(x);
     86                        *pX=X;
     87                #endif
     88        #endif
    7889}/*}}}*/

  • issm/trunk-jpl/src/c/shared/Alloc/alloc_module.cpp

    r12832 r12850  
    2222                /*Actually, still get rid of internal Petsc matrix. Quick fix until Matlab handles C++
    2323                 * correctly: */
    24                 #ifdef _HAVE_PETSC_
     24                /*Keeping this for legacy reasons: */
     25                /*#ifdef _HAVE_PETSC_
    2526                        MatFree(&(*pv)->matrix);
    26                 #endif
     27                #endif*/
     28                delete *pv;
    2729                *pv=NULL;
    2830        }
    … …  
    3436                /*Actually, still get rid of internal Petsc vector. Quick fix until Matlab handles C++
    3537                 * correctly: */
    36                 #ifdef _HAVE_PETSC_
     38                /*Keeping this for legacy reasons: shoudl be done by the vector itself!*/
     39                /*#ifdef _HAVE_PETSC_
    3740                        VecFree(&(*pv)->vector);
    38                 #endif
     41                #endif*/
     42                delete *pv;
    3943                *pv=NULL;
    4044        }

  • issm/trunk-jpl/src/c/toolkits/petsc/petscincludes.h

    r11695 r12850  
    1616/*our own patches: */
    1717#include "patches/petscpatches.h"
     18#include "objects/petscobjects.h"
    1819
    1920#endif
Note: See TracChangeset for help on using the changeset viewer.