Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 12849) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 12850) @@ -331,4 +331,5 @@ ./modules/Solverx/Solverx.cpp\ ./modules/Solverx/Solverx.h\ + ./modules/Solverx/SolverxSeq.cpp\ ./modules/VecMergex/VecMergex.cpp\ ./modules/VecMergex/VecMergex.h\ @@ -784,4 +785,9 @@ ./toolkits/petsc/patches/ISSMToPetscInsertMode.cpp\ ./toolkits/petsc/patches/ISSMToPetscNormMode.cpp\ + ./toolkits/petsc/objects/petscobjects.h\ + ./toolkits/petsc/objects/PetscMat.h\ + ./toolkits/petsc/objects/PetscMat.cpp\ + ./toolkits/petsc/objects/PetscVec.h\ + ./toolkits/petsc/objects/PetscVec.cpp\ ./toolkits/petsc/petscincludes.h\ ./shared/Numerics/PetscOptionsFromAnalysis.cpp\ @@ -789,7 +795,4 @@ ./modules/Solverx/DofTypesToIndexSet.cpp -#}}} -#Gsl sources {{{ -gsl_sources= ./modules/Solverx/SolverxGsl.cpp #}}} #Mpi sources {{{ @@ -831,6 +834,6 @@ #}}} #Matlab and Petsc sources {{{ -matlabpetsc_sources= ./matlab/io/MatlabMatrixToPetscMatrix.cpp\ - ./matlab/io/MatlabVectorToPetscVector.cpp +matlabpetsc_sources= ./matlab/io/MatlabMatrixToPetscMat.cpp\ + ./matlab/io/MatlabVectorToPetscVec.cpp #}}} @@ -925,8 +928,4 @@ endif -if GSL -issm_sources += $(gsl_sources) -endif - if TRANSIENT issm_sources += $(transient_sources) Index: /issm/trunk-jpl/src/c/classes/matrix/Matrix.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/Matrix.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/classes/matrix/Matrix.cpp (revision 12850) @@ -26,78 +26,95 @@ Matrix::Matrix(){ - #ifdef _HAVE_PETSC_ - this->matrix=NULL; - #else - this->matrix=NULL; + pmatrix=NULL; + smatrix=NULL; + + type=PetscMatType; //default + #ifndef _HAVE_PETSC_ + type=SeqMatType; #endif - #ifdef _HAVE_ADOLC_ - this->amatrix=NULL; - #endif -} -/*}}}*/ -/*FUNCTION Matrix::Matrix(int M,int N){{{*/ -Matrix::Matrix(int M,int N){ - - #ifdef _HAVE_PETSC_ - this->matrix=NewMat(M,N); - #else - this->matrix=new SeqMat(M,N); - #endif - #ifdef _HAVE_ADOLC_ - this->amatrix=xNew(M*N); - #endif -} -/*}}}*/ -/*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity){{{*/ -Matrix::Matrix(int M,int N,IssmDouble sparsity){ - - #ifdef _HAVE_PETSC_ - this->matrix=NewMat(M,N,sparsity); - #else - this->matrix=new SeqMat(M,N,sparsity); - #endif - #ifdef _HAVE_ADOLC_ - this->amatrix=xNew(M*N); - #endif -} -/*}}}*/ -/*FUNCTION Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){{{*/ -Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity){ - - #ifdef _HAVE_PETSC_ - int i; - - - int* idxm=xNew(M); - int* idxn=xNew(N); - for(i=0;imatrix=NewMat(M,N,sparsity); - MatSetValues(this->matrix,M,idxm,N,idxn,serial_mat,INSERT_VALUES); - MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY); - - xDelete(idxm); - xDelete(idxn); - #else - this->matrix=new SeqMat(serial_mat,M,N,sparsity); - #endif - #ifdef _HAVE_ADOLC_ - this->amatrix=xNew(M*N); - #endif -} -/*}}}*/ -/*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){{{*/ -Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode){ - - #ifdef _HAVE_PETSC_ - this->matrix=NewMat(M,N,connectivity,numberofdofspernode); - #else - this->matrix=new SeqMat(M,N,connectivity,numberofdofspernode); - #endif - #ifdef _HAVE_ADOLC_ - this->amatrix=xNew(M*N); - #endif + +} +/*}}}*/ +/*FUNCTION Matrix::Matrix(int M,int N,int type){{{*/ +Matrix::Matrix(int M,int N,int in_type){ + + pmatrix=NULL; + smatrix=NULL; + type=in_type; + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix=new PetscMat(M,N); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix=new SeqMat(M,N); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + +} +/*}}}*/ +/*FUNCTION Matrix::Matrix(int M,int N,IssmDouble sparsity,int type){{{*/ +Matrix::Matrix(int M,int N,IssmDouble sparsity,int in_type){ + + pmatrix=NULL; + smatrix=NULL; + type=in_type; + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix=new PetscMat(M,N,sparsity); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix=new SeqMat(M,N,sparsity); + } + else _error2_("Matrix type: " << type << " not supported yet!"); +} +/*}}}*/ +/*FUNCTION Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int type){{{*/ +Matrix::Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int in_type){ + + pmatrix=NULL; + smatrix=NULL; + type=in_type; + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix=new PetscMat(serial_mat,M,N,sparsity); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix=new SeqMat(serial_mat,M,N,sparsity); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + +} +/*}}}*/ +/*FUNCTION Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode,int type){{{*/ +Matrix::Matrix(int M,int N,int connectivity,int numberofdofspernode,int in_type){ + + pmatrix=NULL; + smatrix=NULL; + type=in_type; + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix=new PetscMat(M,N,connectivity,numberofdofspernode); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix=new SeqMat(M,N,connectivity,numberofdofspernode); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ @@ -105,12 +122,16 @@ Matrix::~Matrix(){ - #ifdef _HAVE_PETSC_ - MatFree(&this->matrix); - #else - delete this->matrix; - #endif - #ifdef _HAVE_ADOLC_ - xDelete(this->amatrix); - #endif + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + delete this->pmatrix; + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + delete this->smatrix; + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ @@ -120,37 +141,32 @@ void Matrix::Echo(void){ - int i,j; - - #ifdef _HAVE_PETSC_ - MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD); - #else - this->matrix->Echo(); - #endif - - #ifdef _HAVE_ADOLC_ - /*Not sure about that one. Should we use the overloaded operator >>?*/ - _printString_("ADOLC Matrix equivalent:" ); -// for(i=0;ipmatrix->Echo(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->Echo(); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ /*FUNCTION Matrix::Assemble{{{*/ void Matrix::Assemble(void){ - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY); - #if _PETSC_MAJOR_ == 2 - MatCompress(this->matrix); - #endif - #else - this->matrix->Assemble(); - #endif - + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->Assemble(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->Assemble(); + } + else _error2_("Matrix type: " << type << " not supported yet!"); } /*}}}*/ @@ -159,10 +175,17 @@ IssmDouble norm=0; - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - MatNorm(this->matrix,ISSMToPetscNormMode(norm_type),&norm); - #else - norm=this->matrix->Norm(norm_type); - #endif + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + norm=this->pmatrix->Norm(norm_type); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + norm=this->smatrix->Norm(norm_type); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + return norm; } @@ -170,11 +193,17 @@ /*FUNCTION Matrix::GetSize{{{*/ void Matrix::GetSize(int* pM,int* pN){ - - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - MatGetSize(this->matrix,pM,pN); - #else - this->matrix->GetSize(pM,pN); - #endif + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->GetSize(pM,pN); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->GetSize(pM,pN); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ @@ -182,10 +211,16 @@ void Matrix::GetLocalSize(int* pM,int* pN){ - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - MatGetLocalSize(this->matrix,pM,pN); - #else - this->matrix->GetLocalSize(pM,pN); - #endif + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->GetLocalSize(pM,pN); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->GetLocalSize(pM,pN); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ @@ -193,11 +228,16 @@ void Matrix::MatMult(Vector* X,Vector* AX){ - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - _assert_(X->vector); - MatMultPatch(this->matrix,X->vector,AX->vector); - #else - this->matrix->MatMult(X->vector,AX->vector); - #endif + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->MatMult(X->pvector,AX->pvector); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->MatMult(X->svector,AX->svector); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + } /*}}}*/ @@ -209,10 +249,15 @@ output=new Matrix(); - #ifdef _HAVE_PETSC_ - _assert_(this->matrix); - MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix); - #else - output->matrix=this->matrix->Duplicate(); - #endif + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + output->pmatrix=this->pmatrix->Duplicate(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + output->smatrix=this->smatrix->Duplicate(); + } + else _error2_("Matrix type: " << type << " not supported yet!"); return output; @@ -223,10 +268,18 @@ IssmDouble* output=NULL; - - #ifdef _HAVE_PETSC_ - MatToSerial(&output,this->matrix); - #else - output=this->matrix->ToSerial(); - #endif + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + output=this->pmatrix->ToSerial(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + output=this->smatrix->ToSerial(); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + + return output; } @@ -234,10 +287,16 @@ /*FUNCTION Matrix::SetValues{{{*/ void Matrix::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){ - - #ifdef _HAVE_PETSC_ - MatSetValues(this->matrix,m,idxm,n,idxn,values,ISSMToPetscInsertMode(mode)); - #else - this->matrix->SetValues(m,idxm,n,idxn,values,mode); - #endif + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->SetValues(m,idxm,n,idxn,values,mode); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->SetValues(m,idxm,n,idxn,values,mode); + } + else _error2_("Matrix type: " << type << " not supported yet!"); } /*}}}*/ @@ -245,9 +304,16 @@ void Matrix::Convert(MatrixType type){ - #ifdef _HAVE_PETSC_ - MatConvert(this->matrix,ISSMToPetscMatrixType(type),MAT_REUSE_MATRIX,&this->matrix); - #else - this->matrix->Convert(type); - #endif -} -/*}}}*/ + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + this->pmatrix->Convert(type); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqMatType){ + this->smatrix->Convert(type); + } + else _error2_("Matrix type: " << type << " not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/matrix/Matrix.h =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/Matrix.h (revision 12849) +++ /issm/trunk-jpl/src/c/classes/matrix/Matrix.h (revision 12850) @@ -18,26 +18,24 @@ class Vector; +/*}}}*/ -/*}}}*/ +enum matrixtype { PetscMatType, SeqMatType }; class Matrix{ public: - + #ifdef _HAVE_PETSC_ - Mat matrix; - #else - SeqMat* matrix; + PetscMat* pmatrix; #endif - #ifdef _HAVE_ADOLC_ - IssmDouble* amatrix; - #endif + SeqMat* smatrix; + int type; /*Matrix constructors, destructors {{{*/ Matrix(); - Matrix(int M,int N); - Matrix(int M,int N,IssmDouble sparsity); - Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity); - Matrix(int M,int N,int connectivity,int numberofdofspernode); + Matrix(int M,int N,int type=PetscMatType); + Matrix(int M,int N,IssmDouble sparsity,int type=PetscMatType); + Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int type=PetscMatType); + Matrix(int M,int N,int connectivity,int numberofdofspernode,int type=PetscMatType); ~Matrix(); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/classes/matrix/Vector.cpp (revision 12850) @@ -25,90 +25,71 @@ Vector::Vector(){ - #ifdef _HAVE_PETSC_ - this->vector=NULL; - #else - this->vector=NULL; + this->pvector=NULL; + this->svector=NULL; + + type=PetscVecType; //default + #ifndef _HAVE_PETSC_ + type=SeqVecType; #endif - #ifdef _HAVE_ADOLC_ - this->avector=NULL; - #endif -} -/*}}}*/ -/*FUNCTION Vector::Vector(int M,bool fromlocalsize){{{*/ -Vector::Vector(int pM,bool fromlocalsize){ - - #ifdef _HAVE_PETSC_ - this->vector=NewVec(pM,fromlocalsize); - #else - this->vector=new SeqVec(pM,fromlocalsize); - #endif - #ifdef _HAVE_ADOLC_ - this->avector=xNew(pM); - #endif -} -/*}}}*/ -/*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M){{{*/ -Vector::Vector(IssmDouble* serial_vec,int M){ - - #ifdef _HAVE_PETSC_ - int* idxm=xNew(M); - for(int i=0;ivector=NewVec(M); - VecSetValues(this->vector,M,idxm,serial_vec,INSERT_VALUES); - VecAssemblyBegin(this->vector); - VecAssemblyEnd(this->vector); - - xDelete(idxm); - #else - this->vector=new SeqVec(serial_vec,M); - #endif - #ifdef _HAVE_ADOLC_ - this->avector=xNew(M); - #endif -} -/*}}}*/ -#ifdef _HAVE_PETSC_ -/*FUNCTION Vector::Vector(Vec petsc_vec){{{*/ -Vector::Vector(Vec petsc_vec){ - - if(petsc_vec==NULL){ - this->vector=NewVec(0); - } - else{ - /*copy vector*/ - VecDuplicate(petsc_vec,&this->vector); - VecCopy(petsc_vec,this->vector); - } - -} -/*}}}*/ -#endif -#if defined(_HAVE_GSL_) && !defined(_HAVE_PETSC_) -/*FUNCTION Vector::Vector(SeqVec* seq_vec){{{*/ -Vector::Vector(SeqVec* seq_vec){ - - if(seq_vec==NULL){ - this->vector=NULL; - } - else{ - /*copy vector*/ - this->vector=seq_vec->Duplicate(); - } -} -/*}}}*/ -#endif - - /*FUNCTION Vector::~Vector(){{{*/ + +} +/*}}}*/ +/*FUNCTION Vector::Vector(int M,bool fromlocalsize,int type){{{*/ +Vector::Vector(int M,bool fromlocalsize,int in_type){ + + pvector=NULL; + svector=NULL; + type=in_type; + + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector=new PetscVec(M,fromlocalsize); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector=new SeqVec(M,fromlocalsize); + } + else _error2_("Vector type: " << type << " not supported yet!"); + +} +/*}}}*/ +/*FUNCTION Vector::Vector(IssmDouble* serial_vec,int M,int type){{{*/ +Vector::Vector(IssmDouble* serial_vec,int M,int in_type){ + + pvector=NULL; + svector=NULL; + type=in_type; + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector=new PetscVec(serial_vec,M); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector=new SeqVec(serial_vec,M); + } + else _error2_("Vector type: " << type << " not supported yet!"); + +} +/*}}}*/ +/*FUNCTION Vector::~Vector(){{{*/ Vector::~Vector(){ - #ifdef _HAVE_PETSC_ - VecFree(&this->vector); - #else - delete this->vector; - #endif - #ifdef _HAVE_ADOLC_ - xDelete(this->avector); - #endif + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + delete this->pvector; + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + delete this->svector; + } + else _error2_("Vector type: " << type << " not supported yet!"); } /*}}}*/ @@ -118,28 +99,32 @@ void Vector::Echo(void){ - int i; - - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD); - #else - this->vector->Echo(); - #endif - - #ifdef _HAVE_ADOLC_ - /*do nothing for now: */ - #endif + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Echo(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->Echo(); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } /*}}}*/ /*FUNCTION Vector::Assemble{{{*/ void Vector::Assemble(void){ - - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecAssemblyBegin(this->vector); - VecAssemblyEnd(this->vector); - #else - this->vector->Assemble(); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Assemble(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->Assemble(); + } + else _error2_("Vector type: " << type << " not supported yet!"); } @@ -148,11 +133,18 @@ void Vector::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->SetValues(ssize,list,values,mode); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->SetValues(ssize,list,values,mode); + } + else _error2_("Vector type: " << type << " not supported yet!"); + - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode)); - #else - this->vector->SetValues(ssize,list,values,mode); - #endif } @@ -160,11 +152,16 @@ /*FUNCTION Vector::SetValue{{{*/ void Vector::SetValue(int dof, IssmDouble value, InsMode mode){ - - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode)); - #else - this->vector->SetValue(dof,value,mode); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->SetValue(dof,value,mode); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->SetValue(dof,value,mode); + } + else _error2_("Vector type: " << type << " not supported yet!"); } @@ -173,21 +170,33 @@ void Vector::GetValue(IssmDouble* pvalue,int dof){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecGetValues(this->vector,1,&dof,pvalue); - #else - this->vector->GetValue(pvalue,dof); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->GetValue(pvalue,dof); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->GetValue(pvalue,dof); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } /*}}}*/ /*FUNCTION Vector::GetSize{{{*/ void Vector::GetSize(int* pM){ - - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecGetSize(this->vector,pM); - #else - this->vector->GetSize(pM); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->GetSize(pM); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->GetSize(pM); + } + else _error2_("Vector type: " << type << " not supported yet!"); } @@ -198,5 +207,4 @@ int M; - _assert_(this->vector); this->GetSize(&M); @@ -210,10 +218,16 @@ void Vector::GetLocalSize(int* pM){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecGetLocalSize(this->vector,pM); - #else - this->vector->GetLocalSize(pM); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->GetLocalSize(pM); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->GetLocalSize(pM); + } + else _error2_("Vector type: " << type << " not supported yet!"); } @@ -223,14 +237,20 @@ Vector* output=NULL; - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - Vec vec_output=NULL; - VecDuplicate(this->vector,&vec_output); - output=new Vector(vec_output); - VecFree(&vec_output); - #else + + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ output=new Vector(); - output->vector=this->vector->Duplicate(); - #endif + output->pvector=this->pvector->Duplicate(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + output=new Vector(); + output->svector=this->svector->Duplicate(); + } + else _error2_("Vector type: " << type << " not supported yet!"); + return output; @@ -240,10 +260,16 @@ void Vector::Set(IssmDouble value){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecSet(this->vector,value); - #else - this->vector->Set(value); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Set(value); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->Set(value); + } + else _error2_("Vector type: " << type << " not supported yet!"); } @@ -252,10 +278,17 @@ void Vector::AXPY(Vector* X, IssmDouble a){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecAXPY(this->vector,a,X->vector); - #else - this->vector->AXPY(X->vector,a); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->AXPY(X->pvector,a); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->AXPY(X->svector,a); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } /*}}}*/ @@ -263,10 +296,17 @@ void Vector::AYPX(Vector* X, IssmDouble a){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecAYPX(this->vector,a,X->vector); - #else - this->vector->AYPX(X->vector,a); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->AYPX(X->pvector,a); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->AYPX(X->svector,a); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } @@ -277,9 +317,15 @@ IssmDouble* vec_serial=NULL; - #ifdef _HAVE_PETSC_ - VecToMPISerial(&vec_serial, this->vector); - #else - vec_serial=this->vector->ToMPISerial(); - #endif + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + vec_serial=this->pvector->ToMPISerial(); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + vec_serial=this->svector->ToMPISerial(); + } + else _error2_("Vector type: " << type << " not supported yet!"); return vec_serial; @@ -290,9 +336,17 @@ void Vector::Copy(Vector* to){ - #ifdef _HAVE_PETSC_ - if(this->vector) VecCopy(this->vector,to->vector); - #else - this->vector->Copy(to->vector); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Copy(to->pvector); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->Copy(to->svector); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } @@ -302,10 +356,17 @@ IssmDouble norm=0; - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm); - #else - norm=this->vector->Norm(norm_type); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + norm=this->pvector->Norm(norm_type); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + norm=this->svector->Norm(norm_type); + } + else _error2_("Vector type: " << type << " not supported yet!"); + return norm; } @@ -314,10 +375,17 @@ void Vector::Scale(IssmDouble scale_factor){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecScale(this->vector,scale_factor); - #else - this->vector->Scale(scale_factor); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Scale(scale_factor); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->Scale(scale_factor); + } + else _error2_("Vector type: " << type << " not supported yet!"); + } /*}}}*/ @@ -326,10 +394,17 @@ IssmDouble dot; - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecDot(this->vector,vector->vector,&dot); - #else - dot=this->vector->Dot(vector->vector); - #endif + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + dot=this->pvector->Dot(vector->pvector); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + dot=this->svector->Dot(vector->svector); + } + else _error2_("Vector type: " << type << " not supported yet!"); + return dot; } @@ -338,10 +413,17 @@ void Vector::PointwiseDivide(Vector* x,Vector* y){ - #ifdef _HAVE_PETSC_ - _assert_(this->vector); - VecPointwiseDivide(this->vector,x->vector,y->vector); - #else - this->vector->PointwiseDivide(x->vector,y->vector); - #endif -} -/*}}}*/ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->PointwiseDivide(x->pvector,y->pvector); + #else + _error2_("Petsc matrix format not usable, as Petsc has not been compiled!"); + #endif + } + else if(type==SeqVecType){ + this->svector->PointwiseDivide(x->svector,y->svector); + } + else _error2_("Vector type: " << type << " not supported yet!"); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/matrix/Vector.h =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/Vector.h (revision 12849) +++ /issm/trunk-jpl/src/c/classes/matrix/Vector.h (revision 12850) @@ -19,31 +19,25 @@ /*}}}*/ +enum vectortype { PetscVecType, SeqVecType }; + class Vector{ public: + + #ifdef _HAVE_PETSC_ + PetscVec* pvector; + #endif + SeqVec* svector; + int type; - #ifdef _HAVE_PETSC_ - Vec vector; - #else - SeqVec* vector; - #endif - #ifdef _HAVE_ADOLC_ - IssmDouble* avector; - #endif /*Vector constructors, destructors {{{*/ Vector(); - Vector(int M,bool fromlocalsize=false); - Vector(IssmDouble* serial_vec,int pM); - #ifdef _HAVE_PETSC_ - Vector(Vec petsc_vec); - #endif - #if defined(_HAVE_GSL_) && !defined(_HAVE_PETSC_) - Vector(SeqVec* seq_vec); - #endif + Vector(int M,bool fromlocalsize=false,int type=PetscVecType); + Vector(IssmDouble* serial_vec,int pM,int type=PetscVecType); ~Vector(); /*}}}*/ /*Vector specific routines {{{*/ - void Echo(void); + void Echo(void); void AXPY(Vector *X, IssmDouble a); void AYPX(Vector *X, IssmDouble a); Index: /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToMatrix.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToMatrix.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToMatrix.cpp (revision 12850) @@ -31,7 +31,7 @@ #ifdef _HAVE_PETSC_ - MatlabMatrixToPetscMatrix(&matrix->matrix,NULL,NULL,mxmatrix); + matrix->pmatrix=MatlabMatrixToPetscMat(mxmatrix); #else - matrix->matrix=MatlabMatrixToSeqMat(mxmatrix); + matrix->smatrix=MatlabMatrixToSeqMat(mxmatrix); #endif Index: /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMat.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMat.cpp (revision 12850) +++ /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMat.cpp (revision 12850) @@ -0,0 +1,124 @@ +/* \file MatlabMatrixToPetscMatrix.cpp + * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix: + */ + + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "../../shared/shared.h" + +/*Petsc includes: */ +#include "petscmat.h" +#include "petscvec.h" +#include "petscksp.h" + +/*Matlab includes: */ +#include "mex.h" +#include "matlabio.h" + +PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix){ + + int dummy; + PetscMat* matrix=new PetscMat(); + + MatlabMatrixToPetscMat(&matrix->matrix, &dummy, &dummy, mxmatrix); + + return matrix; +} +int MatlabMatrixToPetscMat(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){ + + int rows, cols; + double *mxmatrix_ptr = NULL; + double *tmatrix = NULL; + int ierr; + int i,j; + + /*output: */ + Mat matrix = NULL; + + /*matlab indices: */ + mwIndex *ir = NULL; + mwIndex *jc = NULL; + double *pr = NULL; + int count; + int nnz; + int nz; + + /*petsc indices: */ + int *idxm = NULL; + int *idxn = NULL; + + /*Ok, first check if we are dealing with a sparse or full matrix: */ + if (mxIsSparse(mxmatrix)){ + + /*Dealing with sparse matrix: recover size first: */ + mxmatrix_ptr=(double*)mxGetPr(mxmatrix); + rows=mxGetM(mxmatrix); + cols=mxGetN(mxmatrix); + nnz=mxGetNzmax(mxmatrix); + if(rows){ + nz=(int)((double)nnz/(double)rows); + } + else{ + nz=0; + } + + ierr=MatCreateSeqAIJ(PETSC_COMM_SELF,rows,cols,nz,PETSC_NULL,&matrix);CHKERRQ(ierr); + + /*Now, get ir,jc and pr: */ + pr=mxGetPr(mxmatrix); + ir=mxGetIr(mxmatrix); + jc=mxGetJc(mxmatrix); + + /*Now, start inserting data into sparse matrix: */ + count=0; + for(i=0;i(rows*cols); + for(i=0;i(rows); + idxn=xNew(cols); + + for(i=0;i(tmatrix); + } + + /*Assemble matrix: */ + MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY); + + /*Assign output pointer: */ + *pmatrix=matrix; + if(pmatrix_rows) *pmatrix_rows=rows; + if(pmatrix_cols) *pmatrix_cols=cols; + + return 1; +} Index: sm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMatrix.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/MatlabMatrixToPetscMatrix.cpp (revision 12849) +++ (revision ) @@ -1,114 +1,0 @@ -/* \file MatlabMatrixToPetscMatrix.cpp - * \brief: convert a sparse or dense matlab matrix to a serial Petsc matrix: - */ - - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif -#include "../../shared/shared.h" - -/*Petsc includes: */ -#include "petscmat.h" -#include "petscvec.h" -#include "petscksp.h" - -/*Matlab includes: */ -#include "mex.h" - -int MatlabMatrixToPetscMatrix(Mat* pmatrix,int* pmatrix_rows,int* pmatrix_cols,const mxArray* mxmatrix){ - - int rows, cols; - double *mxmatrix_ptr = NULL; - double *tmatrix = NULL; - int ierr; - int i,j; - - /*output: */ - Mat matrix = NULL; - - /*matlab indices: */ - mwIndex *ir = NULL; - mwIndex *jc = NULL; - double *pr = NULL; - int count; - int nnz; - int nz; - - /*petsc indices: */ - int *idxm = NULL; - int *idxn = NULL; - - /*Ok, first check if we are dealing with a sparse or full matrix: */ - if (mxIsSparse(mxmatrix)){ - - /*Dealing with sparse matrix: recover size first: */ - mxmatrix_ptr=(double*)mxGetPr(mxmatrix); - rows=mxGetM(mxmatrix); - cols=mxGetN(mxmatrix); - nnz=mxGetNzmax(mxmatrix); - if(rows){ - nz=(int)((double)nnz/(double)rows); - } - else{ - nz=0; - } - - ierr=MatCreateSeqAIJ(PETSC_COMM_SELF,rows,cols,nz,PETSC_NULL,&matrix);CHKERRQ(ierr); - - /*Now, get ir,jc and pr: */ - pr=mxGetPr(mxmatrix); - ir=mxGetIr(mxmatrix); - jc=mxGetJc(mxmatrix); - - /*Now, start inserting data into sparse matrix: */ - count=0; - for(i=0;i(rows*cols); - for(i=0;i(rows); - idxn=xNew(cols); - - for(i=0;i(tmatrix); - } - - /*Assemble matrix: */ - MatAssemblyBegin(matrix,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(matrix,MAT_FINAL_ASSEMBLY); - - /*Assign output pointer: */ - *pmatrix=matrix; - if(pmatrix_rows) *pmatrix_rows=rows; - if(pmatrix_cols) *pmatrix_cols=cols; - - return 1; -} Index: /issm/trunk-jpl/src/c/matlab/io/MatlabVectorToPetscVec.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/MatlabVectorToPetscVec.cpp (revision 12850) +++ /issm/trunk-jpl/src/c/matlab/io/MatlabVectorToPetscVec.cpp (revision 12850) @@ -0,0 +1,109 @@ +/* \file MatlabVectorToPetscVector.cpp + * \brief: convert a sparse or dense matlab vector to a serial Petsc vector: + */ + + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +/*Petsc includes: */ +#include "petscmat.h" +#include "petscvec.h" +#include "petscksp.h" + +/*Matlab includes: */ +#include "mex.h" + +#include "../../shared/shared.h" +#include "matlabio.h" + +PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector){ + + int dummy; + PetscVec* vector=new PetscVec(); + + MatlabVectorToPetscVec(&vector->vector,&dummy, mxvector); + + return vector; +} + +int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector){ + + int rows, cols; + double* mxvector_ptr=NULL; + int ierr; + int i,j; + + /*output: */ + Vec vector=NULL; + + /*matlab indices: */ + mwIndex* ir=NULL; + mwIndex* jc=NULL; + double* pr=NULL; + int count; + int nnz; + int nz; + + /*petsc indices: */ + int* idxm=NULL; + + /*Ok, first check if we are dealing with a sparse or full vector: */ + if (mxIsSparse(mxvector)){ + + /*Dealing with sparse vector: recover size first: */ + mxvector_ptr=(double*)mxGetPr(mxvector); + rows=mxGetM(mxvector); + cols=mxGetN(mxvector); + nnz=mxGetNzmax(mxvector); + nz=(int)((double)nnz/(double)rows); + + ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr); + + /*Now, get ir,jc and pr: */ + pr=mxGetPr(mxvector); + ir=mxGetIr(mxvector); + jc=mxGetJc(mxvector); + + /*Now, start inserting data into sparse vector: */ + count=0; + for(i=0;i(rows); + + for(i=0;i -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -/*Petsc includes: */ -#include "petscmat.h" -#include "petscvec.h" -#include "petscksp.h" - -/*Matlab includes: */ -#include "mex.h" - -#include "../../shared/shared.h" - -int MatlabVectorToPetscVector(Vec* pvector,int* pvector_rows,const mxArray* mxvector){ - - int rows, cols; - double* mxvector_ptr=NULL; - int ierr; - int i,j; - - /*output: */ - Vec vector=NULL; - - /*matlab indices: */ - mwIndex* ir=NULL; - mwIndex* jc=NULL; - double* pr=NULL; - int count; - int nnz; - int nz; - - /*petsc indices: */ - int* idxm=NULL; - - /*Ok, first check if we are dealing with a sparse or full vector: */ - if (mxIsSparse(mxvector)){ - - /*Dealing with sparse vector: recover size first: */ - mxvector_ptr=(double*)mxGetPr(mxvector); - rows=mxGetM(mxvector); - cols=mxGetN(mxvector); - nnz=mxGetNzmax(mxvector); - nz=(int)((double)nnz/(double)rows); - - ierr=VecCreateSeq(PETSC_COMM_SELF,rows,&vector);CHKERRQ(ierr); - - /*Now, get ir,jc and pr: */ - pr=mxGetPr(mxvector); - ir=mxGetIr(mxvector); - jc=mxGetJc(mxvector); - - /*Now, start inserting data into sparse vector: */ - count=0; - for(i=0;i(rows); - - for(i=0;ivector,&dummy,mxvector); + vector->pvector=MatlabVectorToPetscVec(mxvector); #else - vector->vector=MatlabVectorToSeqVec(mxvector); + vector->svector=MatlabVectorToSeqVec(mxvector); #endif Index: /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/matlab/io/WriteMatlabData.cpp (revision 12850) @@ -227,10 +227,6 @@ if(matrix){ - #ifdef _HAVE_PETSC_ - PetscMatrixToDoubleMatrix(&matrix_ptr,&rows,&cols,matrix->matrix); - #else - matrix_ptr=matrix->matrix->ToSerial(); - matrix->matrix->GetSize(&rows,&cols); - #endif + matrix_ptr=matrix->ToSerial(); + matrix->GetSize(&rows,&cols); /*Now transpose the matrix and allocate with Matlab's memory manager: */ @@ -269,10 +265,6 @@ if(vector){ /*call toolkit routine: */ - #ifdef _HAVE_PETSC_ - PetscVectorToDoubleVector(&vector_ptr,&rows,vector->vector); - #else - vector_ptr=vector->vector->ToMPISerial(); - vector->vector->GetSize(&rows); - #endif + vector_ptr=vector->ToMPISerial(); + vector->GetSize(&rows); /*now create the matlab vector with Matlab's memory manager */ Index: /issm/trunk-jpl/src/c/matlab/io/matlabio.h =================================================================== --- /issm/trunk-jpl/src/c/matlab/io/matlabio.h (revision 12849) +++ /issm/trunk-jpl/src/c/matlab/io/matlabio.h (revision 12850) @@ -79,6 +79,8 @@ /*Matlab to Petsc routines: */ #ifdef _HAVE_PETSC_ -int MatlabMatrixToPetscMatrix(Mat* matrix,int* prows,int* pcols, const mxArray* mxmatrix); -int MatlabVectorToPetscVector(Vec* pvector,int* pvector_rows,const mxArray* mxvector); +int MatlabMatrixToPetscMat(Mat* matrix,int* prows,int* pcols, const mxArray* mxmatrix); +PetscMat* MatlabMatrixToPetscMat(const mxArray* mxmatrix); +int MatlabVectorToPetscVec(Vec* pvector,int* pvector_rows,const mxArray* mxvector); +PetscVec* MatlabVectorToPetscVec(const mxArray* mxvector); #endif Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.cpp (revision 12850) @@ -26,36 +26,20 @@ _assert_(pf); - #ifdef _HAVE_PETSC_ - Vec uf0_vector = NULL; - Vec df_vector = NULL; - Vec uf_vector = NULL; - if(uf0) uf0_vector = uf0->vector; - if(df) df_vector = df->vector; + /*Initialize vector: */ + uf=new Vector(); - /*In serial mode, the Petsc Options database has not been initialized properly: */ + /*According to matrix type, use specific solvers: */ + if(Kff->type==PetscMatType){ + PetscVec* uf0_vector = NULL; + PetscVec* df_vector = NULL; + if(uf0) uf0_vector = uf0->pvector; + if(df) df_vector = df->pvector; - SolverxPetsc(&uf_vector,Kff->matrix,pf->vector,uf0_vector,df_vector,parameters); - - /*Create vector out of petsc vector: */ - uf=new Vector(uf_vector); - - /*Free ressources: */ - VecFree(&uf_vector); - #else - #ifdef _HAVE_GSL_ - SeqVec* uf_vector=NULL; - - SolverxGsl(&uf_vector,Kff->matrix,pf->vector); - - /*Create vector out of SeqVec vector: */ - uf=new Vector(uf_vector); - - /*Free ressources: */ - delete uf_vector; - - #else - _error2_("GSL support not compiled in!"); - #endif - #endif + SolverxPetsc(&uf->pvector,Kff->pmatrix,pf->pvector,uf0_vector,df_vector,parameters); + } + else if(Kff->type==SeqMatType){ + SolverxSeq(&uf->svector,Kff->smatrix,pf->svector); + } + else _error2_("Matrix type: " << Kff->type << " not supported yet!"); /*Assign output pointers: */ Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 12849) +++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h (revision 12850) @@ -18,12 +18,11 @@ #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 -#ifdef _HAVE_GSL_ -void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf); -void SolverxGsl(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n); -#endif +void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf); +void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n); #endif /* _SOLVERX_H */ Index: sm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/SolverxGsl.cpp (revision 12849) +++ (revision ) @@ -1,78 +1,0 @@ -/*!\file SolverxGsl - * \brief Gsl implementation of solver - */ - -#ifdef HAVE_CONFIG_H - #include -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif -#include -#include - -#include "./Solverx.h" -#include "../../shared/shared.h" -#include "../../include/include.h" -#include "../../io/io.h" - -void SolverxGsl(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/ - - /*Intermediary: */ - int M,N,N2,s; - SeqVec *uf = NULL; - IssmDouble *x = NULL; - - Kff->GetSize(&M,&N); - pf->GetSize(&N2); - - if(N!=N2)_error2_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !"); - if(M!=N)_error2_("Stiffness matrix should be square!"); - - SolverxGsl(&x,Kff->matrix,pf->vector,N); - uf=new SeqVec(x,N); - - /*Assign output pointers:*/ - *puf=uf; -}/*}}}*/ -void SolverxGsl(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/ - - /*GSL Matrices and vectors: */ - int s; - gsl_matrix_view a; - gsl_vector_view b; - gsl_vector *x = NULL; - gsl_permutation *p = NULL; -#ifdef _HAVE_ADOLC_ - // if we use Adol-C then the IssmDouble will be an adouble - // and the calls to gsl_... will not work - // and we should call a suitable wrapped solve instead - _error2_("SolverxGsl: should not be here with Adol-C"); -#else - /*A will be modified by LU decomposition. Use copy*/ - IssmDouble* Acopy = xNew(n*n); - xMemCpy(Acopy,A,n*n); - - /*Initialize gsl matrices and vectors: */ - a = gsl_matrix_view_array (Acopy,n,n); - b = gsl_vector_view_array (B,n); - x = gsl_vector_alloc (n); - - /*Run LU and solve: */ - p = gsl_permutation_alloc (n); - gsl_linalg_LU_decomp (&a.matrix, p, &s); - gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x); - - //printf ("x = \n"); - //gsl_vector_fprintf (stdout, x, "%g"); - - /*Copy result*/ - IssmDouble* X = xNew(n); - memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble)); - - /*Clean up and assign output pointer*/ - xDelete(Acopy); - gsl_permutation_free(p); - gsl_vector_free(x); - *pX=X; -#endif -}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxPetsc.cpp (revision 12850) @@ -14,4 +14,11 @@ #endif +void SolverxPetsc(PetscVec** puf, PetscMat* Kff, PetscVec* pf, PetscVec* uf0,PetscVec* df, Parameters* parameters){ + + PetscVec* uf=new PetscVec(); + SolverxPetsc(&uf->vector, Kff->matrix, pf->vector, uf0->vector, df->vector, parameters); + *puf=uf; + +} void SolverxPetsc(Vec* puf, Mat Kff, Vec pf, Vec uf0,Vec df, Parameters* parameters){ Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 12850) +++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp (revision 12850) @@ -0,0 +1,89 @@ +/*!\file SolverxSeq + * \brief implementation of sequential solver using the GSL librarie + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include + +#include "./Solverx.h" +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../io/io.h" + +#ifdef _HAVE_GSL_ +#include +#endif + +void SolverxSeq(SeqVec** puf,SeqMat* Kff, SeqVec* pf){/*{{{*/ + + #ifdef _HAVE_GSL_ + /*Intermediary: */ + int M,N,N2,s; + SeqVec *uf = NULL; + IssmDouble *x = NULL; + + Kff->GetSize(&M,&N); + pf->GetSize(&N2); + + if(N!=N2)_error2_("Right hand side vector of size " << N2 << ", when matrix is of size " << M << "-" << N << " !"); + if(M!=N)_error2_("Stiffness matrix should be square!"); + + SolverxSeq(&x,Kff->matrix,pf->vector,N); + uf=new SeqVec(x,N); + + /*Assign output pointers:*/ + *puf=uf; + + #else + _error2_("GSL support not compiled in!"); + #endif + +}/*}}}*/ +void SolverxSeq(IssmDouble** pX,IssmDouble* A,IssmDouble* B,int n){/*{{{*/ + + #ifdef _HAVE_GSL_ + /*GSL Matrices and vectors: */ + int s; + gsl_matrix_view a; + gsl_vector_view b; + gsl_vector *x = NULL; + gsl_permutation *p = NULL; + #ifdef _HAVE_ADOLC_ + // if we use Adol-C then the IssmDouble will be an adouble + // and the calls to gsl_... will not work + // and we should call a suitable wrapped solve instead + _error2_("SolverxSeq: should not be here with Adol-C"); + #else + /*A will be modified by LU decomposition. Use copy*/ + IssmDouble* Acopy = xNew(n*n); + xMemCpy(Acopy,A,n*n); + + /*Initialize gsl matrices and vectors: */ + a = gsl_matrix_view_array (Acopy,n,n); + b = gsl_vector_view_array (B,n); + x = gsl_vector_alloc (n); + + /*Run LU and solve: */ + p = gsl_permutation_alloc (n); + gsl_linalg_LU_decomp (&a.matrix, p, &s); + gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x); + + //printf ("x = \n"); + //gsl_vector_fprintf (stdout, x, "%g"); + + /*Copy result*/ + IssmDouble* X = xNew(n); + memcpy(X,gsl_vector_ptr(x,0),n*sizeof(IssmDouble)); + + /*Clean up and assign output pointer*/ + xDelete(Acopy); + gsl_permutation_free(p); + gsl_vector_free(x); + *pX=X; + #endif + #endif +}/*}}}*/ Index: /issm/trunk-jpl/src/c/shared/Alloc/alloc_module.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Alloc/alloc_module.cpp (revision 12849) +++ /issm/trunk-jpl/src/c/shared/Alloc/alloc_module.cpp (revision 12850) @@ -22,7 +22,9 @@ /*Actually, still get rid of internal Petsc matrix. Quick fix until Matlab handles C++ * correctly: */ - #ifdef _HAVE_PETSC_ + /*Keeping this for legacy reasons: */ + /*#ifdef _HAVE_PETSC_ MatFree(&(*pv)->matrix); - #endif + #endif*/ + delete *pv; *pv=NULL; } @@ -34,7 +36,9 @@ /*Actually, still get rid of internal Petsc vector. Quick fix until Matlab handles C++ * correctly: */ - #ifdef _HAVE_PETSC_ + /*Keeping this for legacy reasons: shoudl be done by the vector itself!*/ + /*#ifdef _HAVE_PETSC_ VecFree(&(*pv)->vector); - #endif + #endif*/ + delete *pv; *pv=NULL; } Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp (revision 12850) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp (revision 12850) @@ -0,0 +1,174 @@ +/*!\file PetscMat.cpp + * \brief: implementation of the PetscMat object + */ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../petscincludes.h" +#include "../../../shared/shared.h" + +/*}}}*/ + +/*PetscMat constructors and destructor*/ +/*FUNCTION PetscMat::PetscMat(){{{*/ +PetscMat::PetscMat(){ + this->matrix=NULL; + #ifdef _HAVE_ADOLC_ + this->amatrix=NULL; + #endif + +} +/*}}}*/ +/*FUNCTION PetscMat::PetscMat(int M,int N){{{*/ +PetscMat::PetscMat(int M,int N){ + + this->matrix=NULL; + if(M*N)this->matrix=NewMat(M,N); +} +/*}}}*/ +/*FUNCTION PetscMat::PetscMat(int M,int N, IssmDouble sparsity){{{*/ +PetscMat::PetscMat(int M,int N, IssmDouble sparsity){ + + this->matrix=NULL; + if(M*N) this->matrix=NewMat(M,N,sparsity); +} +/*}}}*/ +/*FUNCTION PetscMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity){{{*/ +PetscMat::PetscMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity){ + + this->matrix=NULL; + + if(M*N){ + int i; + int* idxm=xNew(M); + int* idxn=xNew(N); + + for(i=0;imatrix=NewMat(M,N,sparsity); + MatSetValues(this->matrix,M,idxm,N,idxn,serial_mat,INSERT_VALUES); + MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY); + + xDelete(idxm); + xDelete(idxn); + } + +} +/*}}}*/ +/*FUNCTION PetscMat::PetscMat(int M,int N, int connectivity, int numberofdofspernode){{{*/ +PetscMat::PetscMat(int M,int N, int connectivity,int numberofdofspernode){ + + this->matrix=NULL; + if(M*N) this->matrix=NewMat(M,N,connectivity,numberofdofspernode); + +} +/*}}}*/ +/*FUNCTION PetscMat::~PetscMat(){{{*/ +PetscMat::~PetscMat(){ + MatFree(&this->matrix); +} +/*}}}*/ + +/*PetscMat specific routines: */ +/*FUNCTION PetscMat::Echo{{{*/ +void PetscMat::Echo(void){ + + MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD); +} +/*}}}*/ +/*FUNCTION PetscMat::Assemble{{{*/ +void PetscMat::Assemble(void){ + + _assert_(this->matrix); + MatAssemblyBegin(this->matrix,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(this->matrix,MAT_FINAL_ASSEMBLY); + #if _PETSC_MAJOR_ == 2 + MatCompress(this->matrix); + #endif + +} +/*}}}*/ +/*FUNCTION PetscMat::Norm{{{*/ +IssmDouble PetscMat::Norm(NormMode mode){ + + + IssmDouble norm=0; + _assert_(this->matrix); + MatNorm(this->matrix,ISSMToPetscNormMode(mode),&norm); + + return norm; + +} +/*}}}*/ +/*FUNCTION PetscMat::GetSize{{{*/ +void PetscMat::GetSize(int* pM,int* pN){ + + _assert_(this->matrix); + MatGetSize(this->matrix,pM,pN); +} +/*}}}*/ +/*FUNCTION PetscMat::GetLocalSize{{{*/ +void PetscMat::GetLocalSize(int* pM,int* pN){ + + _assert_(this->matrix); + MatGetLocalSize(this->matrix,pM,pN); + +} +/*}}}*/ +/*FUNCTION PetscMat::MatMult{{{*/ +void PetscMat::MatMult(PetscVec* X,PetscVec* AX){ + + _assert_(this->matrix); + _assert_(X->vector); + MatMultPatch(this->matrix,X->vector,AX->vector); + +} +/*}}}*/ +/*FUNCTION PetscMat::Duplicate{{{*/ +PetscMat* PetscMat::Duplicate(void){ + + PetscMat* output=NULL; + + output=new PetscMat(); + _assert_(this->matrix); + MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix); + + return output; + +} +/*}}}*/ +/*FUNCTION PetscMat::ToSerial{{{*/ +IssmDouble* PetscMat::ToSerial(void){ + + IssmDouble* output=NULL; + + MatToSerial(&output,this->matrix); + return output; + +} +/*}}}*/ +/*FUNCTION PetscMat::SetValues{{{*/ +void PetscMat::SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){ + + MatSetValues(this->matrix,m,idxm,n,idxn,values,ISSMToPetscInsertMode(mode)); + +} +/*}}}*/ +/*FUNCTION PetscMat::Convert{{{*/ +void PetscMat::Convert(MatrixType type){ + + MatConvert(this->matrix,ISSMToPetscMatrixType(type),MAT_REUSE_MATRIX,&this->matrix); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h (revision 12850) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h (revision 12850) @@ -0,0 +1,56 @@ +/*!\file: PetscMat.h + * \brief wrapper to our own PetscMat object, which is needed to add AD capabilities (using ADOLC) + * to a C-coded Petsc API. We are just wrapping the Petsc objects into C++ equivalent, so that + * later, we can map all of the Petsc routines into Adolc equivalents. + */ + +#ifndef _PETSCMAT_H_ +#define _PETSCMAT_H_ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../petscincludes.h" +#include "../../../include/include.h" + +/*}}}*/ +class PetscVec; + +class PetscMat{ + + public: + Mat matrix; + + #ifdef _HAVE_ADOLC_ + IssmDouble* amatrix; + #endif + + /*PetscMat constructors, destructors {{{*/ + PetscMat(); + PetscMat(int M,int N); + PetscMat(int M,int N,IssmDouble sparsity); + PetscMat(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity); + PetscMat(int M,int N,int connectivity,int numberofdofspernode); + ~PetscMat(); + /*}}}*/ + /*PetscMat specific routines {{{*/ + void Echo(void); + void Assemble(void); + IssmDouble Norm(NormMode norm_type); + void GetSize(int* pM,int* pN); + void GetLocalSize(int* pM,int* pN); + void MatMult(PetscVec* X,PetscVec* AX); + PetscMat* Duplicate(void); + IssmDouble* ToSerial(void); + void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode); + void Convert(MatrixType type); + /*}}}*/ + +}; + +#endif //#ifndef _PETSCMAT_H_ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 12850) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 12850) @@ -0,0 +1,217 @@ +/*!\file PetscVec.cpp + * \brief: implementation of the PetscVec object + */ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include +#include "../petscincludes.h" +#include "../../../shared/shared.h" + +/*}}}*/ + +/*PetscVec constructors and destructor*/ +/*FUNCTION PetscVec::PetscVec(){{{*/ +PetscVec::PetscVec(){ + this->vector=NULL; + #ifdef _HAVE_ADOLC_ + this->avector=NULL; + #endif +} +/*}}}*/ +/*FUNCTION PetscVec::PetscVec(int M,bool fromlocalsize){{{*/ +PetscVec::PetscVec(int M,bool fromlocalsize){ + this->vector=NULL; + if(M) this->vector=NewVec(M,fromlocalsize); +} +/*}}}*/ +/*FUNCTION PetscVec::PetscVec(IssmDouble* serial_vec,int M){{{*/ +PetscVec::PetscVec(IssmDouble* serial_vec,int M){ + + this->vector=NULL; + if(M){ + int* idxm=xNew(M); + for(int i=0;ivector=NewVec(M); + VecSetValues(this->vector,M,idxm,serial_vec,INSERT_VALUES); + VecAssemblyBegin(this->vector); + VecAssemblyEnd(this->vector); + + xDelete(idxm); + } +} +/*}}}*/ +/*FUNCTION PetscVec::~PetscVec(){{{*/ +PetscVec::~PetscVec(){ + VecFree(&this->vector); +} +/*}}}*/ +/*FUNCTION Vector::Vector(Vec petsc_vec){{{*/ +PetscVec::PetscVec(Vec petsc_vec){ + + if(petsc_vec==NULL){ + this->vector=NewVec(0); + } + else{ + /*copy vector*/ + VecDuplicate(petsc_vec,&this->vector); + VecCopy(petsc_vec,this->vector); + } + +} +/*}}}*/ + +/*PetscVec specific routines: */ +/*FUNCTION PetscVec::Echo{{{*/ +void PetscVec::Echo(void){ + + _assert_(this->vector); + VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD); +} +/*}}}*/ +/*FUNCTION PetscVec::Assemble{{{*/ +void PetscVec::Assemble(void){ + + _assert_(this->vector); + VecAssemblyBegin(this->vector); + VecAssemblyEnd(this->vector); + +} +/*}}}*/ +/*FUNCTION PetscVec::SetValues{{{*/ +void PetscVec::SetValues(int ssize, int* list, IssmDouble* values, InsMode mode){ + + _assert_(this->vector); + VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode)); + +} +/*}}}*/ +/*FUNCTION PetscVec::SetValue{{{*/ +void PetscVec::SetValue(int dof, IssmDouble value, InsMode mode){ + + _assert_(this->vector); + VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode)); + +} +/*}}}*/ +/*FUNCTION PetscVec::GetValue{{{*/ +void PetscVec::GetValue(IssmDouble* pvalue,int dof){ + + _assert_(this->vector); + VecGetValues(this->vector,1,&dof,pvalue); + +} +/*}}}*/ +/*FUNCTION PetscVec::GetSize{{{*/ +void PetscVec::GetSize(int* pM){ + + _assert_(this->vector); + VecGetSize(this->vector,pM); + +} +/*}}}*/ +/*FUNCTION PetscVec::GetLocalSize{{{*/ +void PetscVec::GetLocalSize(int* pM){ + + _assert_(this->vector); + VecGetLocalSize(this->vector,pM); + +} +/*}}}*/ +/*FUNCTION PetscVec::Duplicate{{{*/ +PetscVec* PetscVec::Duplicate(void){ + + PetscVec* output=NULL; + _assert_(this->vector); + Vec vec_output=NULL; + VecDuplicate(this->vector,&vec_output); + output=new PetscVec(vec_output); + VecFree(&vec_output); + + return output; +} +/*}}}*/ +/*FUNCTION PetscVec::Set{{{*/ +void PetscVec::Set(IssmDouble value){ + + _assert_(this->vector); + VecSet(this->vector,value); + +} +/*}}}*/ +/*FUNCTION PetscVec::AXPY{{{*/ +void PetscVec::AXPY(PetscVec* X, IssmDouble a){ + + _assert_(this->vector); + VecAXPY(this->vector,a,X->vector); + +} +/*}}}*/ +/*FUNCTION PetscVec::AYPX{{{*/ +void PetscVec::AYPX(PetscVec* X, IssmDouble a){ + + _assert_(this->vector); + VecAYPX(this->vector,a,X->vector); + +} +/*}}}*/ +/*FUNCTION PetscVec::ToMPISerial{{{*/ +IssmDouble* PetscVec::ToMPISerial(void){ + + IssmDouble* vec_serial=NULL; + VecToMPISerial(&vec_serial, this->vector); + return vec_serial; + +} +/*}}}*/ +/*FUNCTION PetscVec::Copy{{{*/ +void PetscVec::Copy(PetscVec* to){ + + if(this->vector) VecCopy(this->vector,to->vector); + +} +/*}}}*/ +/*FUNCTION PetscVec::Norm{{{*/ +IssmDouble PetscVec::Norm(NormMode mode){ + + IssmDouble norm=0; + _assert_(this->vector); + VecNorm(this->vector,ISSMToPetscNormMode(mode),&norm); + return norm; + +} +/*}}}*/ +/*FUNCTION PetscVec::Scale{{{*/ +void PetscVec::Scale(IssmDouble scale_factor){ + + _assert_(this->vector); + VecScale(this->vector,scale_factor); + +} +/*}}}*/ +/*FUNCTION PetscVec::Dot{{{*/ +IssmDouble PetscVec::Dot(PetscVec* input){ + + IssmDouble dot; + _assert_(this->vector); + VecDot(this->vector,input->vector,&dot); + return dot; + +} +/*}}}*/ +/*FUNCTION PetscVec::PointwiseDivide{{{*/ +void PetscVec::PointwiseDivide(PetscVec* x,PetscVec* y){ + + _assert_(this->vector); + VecPointwiseDivide(this->vector,x->vector,y->vector); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 12850) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 12850) @@ -0,0 +1,62 @@ +/*!\file: PetscVec.h + * \brief wrapper to our own PetscVec object, which is needed to add AD capabilities (using ADOLC) + * to a C-coded Petsc API. We are just wrapping the Petsc objects into C++ equivalent, so that + * later, we can map all of the Petsc routines into Adolc equivalents. + */ + +#ifndef _PETSCVEC_H_ +#define _PETSCVEC_H_ + +/*Headers:*/ +/*{{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "../petscincludes.h" +#include "../../../include/include.h" + + +/*}}}*/ + +class PetscVec{ + + public: + Vec vector; + + #ifdef _HAVE_ADOLC_ + IssmDouble* avector; + #endif + + + /*PetscVec constructors, destructors {{{*/ + PetscVec(); + PetscVec(int M,bool fromlocalsize=false); + PetscVec(IssmDouble* buffer, int M); + PetscVec(Vec petsc_vec); + ~PetscVec(); + /*}}}*/ + /*PetscVec specific routines {{{*/ + void Echo(void); + void Assemble(void); + void SetValues(int ssize, int* list, IssmDouble* values, InsMode mode); + void SetValue(int dof, IssmDouble value, InsMode mode); + void GetValue(IssmDouble* pvalue, int dof); + void GetSize(int* pM); + void GetLocalSize(int* pM); + PetscVec* Duplicate(void); + void Set(IssmDouble value); + void AXPY(PetscVec* X, IssmDouble a); + void AYPX(PetscVec* X, IssmDouble a); + IssmDouble* ToMPISerial(void); + void Copy(PetscVec* to); + IssmDouble Norm(NormMode norm_type); + void Scale(IssmDouble scale_factor); + void PointwiseDivide(PetscVec* x,PetscVec* y); + IssmDouble Dot(PetscVec* vector); + /*}}}*/ +}; + +#endif //#ifndef _PETSCVEC_H_ Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h (revision 12850) +++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/petscobjects.h (revision 12850) @@ -0,0 +1,11 @@ +/* \file petscobjects.h + * \brief all includes for our own petsc object implementation + */ + +#ifndef _PETSC_OBJECTS_H_ +#define _PETSC_OBJECTS_H_ + +#include "./PetscMat.h" +#include "./PetscVec.h" + +#endif Index: /issm/trunk-jpl/src/c/toolkits/petsc/petscincludes.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/petsc/petscincludes.h (revision 12849) +++ /issm/trunk-jpl/src/c/toolkits/petsc/petscincludes.h (revision 12850) @@ -16,4 +16,5 @@ /*our own patches: */ #include "patches/petscpatches.h" +#include "objects/petscobjects.h" #endif