Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementMatrix.cpp (revision 11734) @@ -259,47 +259,43 @@ /*In debugging mode, check consistency (no NaN, and values not too big)*/ this->CheckConsistency(); - #ifdef _HAVE_PETSC_ - if(this->dofsymmetrical){ - /*only use row dofs to add values into global matrices: */ - - if(this->row_fsize){ - /*first, retrieve values that are in the f-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_fsize;j++){ - *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); - } - } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Kff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); + if(this->dofsymmetrical){ + /*only use row dofs to add values into global matrices: */ - /*Free ressources:*/ - xfree((void**)&localvalues); + if(this->row_fsize){ + /*first, retrieve values that are in the f-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_fsize;j++){ + *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); + } } + /*add local values into global matrix, using the fglobaldoflist: */ + Kff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL); + /*Free ressources:*/ + xfree((void**)&localvalues); + } - if((this->row_ssize!=0) && (this->row_fsize!=0)){ - /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_ssize;j++){ - *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]); - } - } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Kfs->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,(const double*)localvalues,ADD_VALUES); - /*Free ressources:*/ - xfree((void**)&localvalues); + if((this->row_ssize!=0) && (this->row_fsize!=0)){ + /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_ssize;j++){ + *(localvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]); + } } + /*add local values into global matrix, using the fglobaldoflist: */ + Kfs->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,localvalues,ADD_VAL); + + /*Free ressources:*/ + xfree((void**)&localvalues); } - else{ - _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); - } - #else - _error_("not supported yet!"); - #endif + } + else{ + _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); + } } /*}}}*/ @@ -315,33 +311,29 @@ /*In debugging mode, check consistency (no NaN, and values not too big)*/ this->CheckConsistency(); - #ifdef _HAVE_PETSC_ - if(this->dofsymmetrical){ - /*only use row dofs to add values into global matrices: */ + if(this->dofsymmetrical){ + /*only use row dofs to add values into global matrices: */ - if(this->row_fsize){ - /*first, retrieve values that are in the f-set from the g-set values matrix: */ - localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); - for(i=0;irow_fsize;i++){ - for(j=0;jrow_fsize;j++){ - *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); - } + if(this->row_fsize){ + /*first, retrieve values that are in the f-set from the g-set values matrix: */ + localvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_fsize;j++){ + *(localvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]); } - /*add local values into global matrix, using the fglobaldoflist: */ - MatSetValues(Jff->matrix,this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); } + /*add local values into global matrix, using the fglobaldoflist: */ + Jff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL); + /*Free ressources:*/ + xfree((void**)&localvalues); } - else{ - _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); - } - #else - _error_("not supported yet!"); - #endif + } + else{ + _error_(" non dofsymmetrical matrix AddToGlobal routine not support yet!"); + } + } /*}}}*/ /*FUNCTION ElementMatrix::CheckConsistency{{{1*/ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.cpp (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.cpp (revision 11734) @@ -24,7 +24,6 @@ /*FUNCTION Vector::Vector(){{{1*/ Vector::Vector(){ - this->M=0; #ifdef _HAVE_PETSC_ this->vector=NULL; #else @@ -41,8 +40,7 @@ #ifdef _HAVE_PETSC_ this->vector=NewVec(pM,fromlocalsize); #else - this->M=pM; - this->vector=(double*)xcalloc(pM,sizeof(double)); + this->vector=new SeqVec(pM,fromlocalsize); #endif #ifdef _HAVE_ADOLC_ this->avector=(adouble*)xmalloc(pM*sizeof(adouble)); @@ -50,14 +48,14 @@ } /*}}}*/ /*FUNCTION Vector::Vector(double* serial_vec,int M){{{1*/ -Vector::Vector(double* serial_vec,int pM){ +Vector::Vector(double* serial_vec,int M){ int i,j; #ifdef _HAVE_PETSC_ int* idxm=NULL; - this->vector=NewVec(pM); + this->vector=NewVec(M); idxm=(int*)xmalloc(M*sizeof(int)); for(i=0;iM=pM; - this->vector=(double*)xcalloc(pM,sizeof(double)); + this->vector=new SeqVec(serial_vec,M); #endif #ifdef _HAVE_ADOLC_ - this->avector=(adouble*)xmalloc(pM*sizeof(adouble)); + this->avector=(adouble*)xmalloc(M*sizeof(adouble)); #endif } /*}}}*/ @@ -81,9 +78,6 @@ /*FUNCTION Vector::Vector(Vec petsc_vec){{{1*/ Vector::Vector(Vec petsc_vec){ - /*Get Vector size*/ - VecGetSize(petsc_vec,&this->M); - /*copy vector*/ VecDuplicate(petsc_vec,&this->vector); VecCopy(petsc_vec,this->vector); @@ -97,7 +91,7 @@ #ifdef _HAVE_PETSC_ VecFree(&this->vector); #else - xfree((void**)&this->vector); + delete this->vector; #endif #ifdef _HAVE_ADOLC_ xfree((void**)&this->avector); @@ -114,18 +108,11 @@ #ifdef _HAVE_PETSC_ VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD); #else - printf("Vector size: %i\n",M); - for(i=0;ivector->Echo(); #endif #ifdef _HAVE_ADOLC_ - /*Not sure about that one. Should we use the overloaded operator >>?*/ - printf("ADOLC Vector equivalent:" ); - for(i=0;ivector); #else - _error_("not implemented yet!"); + dataref=this->vector->ToMatlabVector(); #endif return dataref; } /*}}}*/ +/*FUNCTION MatlabVectorToVector{{{1*/ +Vector* MatlabVectorToVector(const mxArray* mxvector){ + + int dummy; + Vector* vector=NULL; + + /*allocate vector object: */ + vector=new Vector(); + + #ifdef _HAVE_PETSC_ + MatlabVectorToPetscVector(&vector->vector,&dummy,mxvector); + #else + vector->vector=MatlabVectorToSeqVec(mxvector); + #endif + + return vector; +} +/*}}}*/ #endif /*FUNCTION Vector::Assemble{{{1*/ void Vector::Assemble(void){ @@ -153,7 +158,7 @@ VecAssemblyBegin(this->vector); VecAssemblyEnd(this->vector); #else - /*do nothing*/ + this->vector->Assemble(); #endif } @@ -166,7 +171,7 @@ _assert_(this->vector); VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode)); #else - _error_("not implemented yet!"); + this->vector->SetValues(ssize,list,values,mode); #endif } @@ -178,7 +183,7 @@ _assert_(this->vector); VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode)); #else - _error_("not implemented yet!"); + this->vector->SetValue(dof,value,mode); #endif } @@ -190,7 +195,7 @@ _assert_(this->vector); VecGetValues(this->vector,1,&dof,pvalue); #else - _error_("not implemented yet!"); + this->vector->GetValue(pvalue,dof); #endif } /*}}}*/ @@ -201,7 +206,7 @@ _assert_(this->vector); VecGetSize(this->vector,pM); #else - _error_("not implemented yet!"); + this->vector->GetSize(pM); #endif } @@ -213,7 +218,7 @@ _assert_(this->vector); VecGetLocalSize(this->vector,pM); #else - _error_("not implemented yet!"); + this->vector->GetLocalSize(pM); #endif } @@ -228,9 +233,8 @@ _assert_(this->vector); VecDuplicate(this->vector,&vec_output); output->vector=vec_output; - VecGetSize(output->vector,&output->M); #else - _error_("not implemented yet!"); + output->vector=this->vector->Duplicate(); #endif return output; @@ -243,7 +247,7 @@ _assert_(this->vector); VecSet(this->vector,value); #else - _error_("not implemented yet!"); + this->vector->Set(value); #endif } @@ -255,7 +259,7 @@ _assert_(this->vector); VecAXPY(this->vector,a,X->vector); #else - _error_("not implemented yet!"); + this->vector->AXPY(X->vector,a); #endif } /*}}}*/ @@ -266,7 +270,7 @@ _assert_(this->vector); VecAYPX(this->vector,a,X->vector); #else - _error_("not implemented yet!"); + this->vector->AYPX(X->vector,a); #endif } @@ -279,7 +283,7 @@ #ifdef _HAVE_PETSC_ VecToMPISerial(&vec_serial, this->vector); #else - _error_("not implemented yet!"); + vec_serial=this->vector->ToMPISerial(); #endif return vec_serial; @@ -292,7 +296,7 @@ #ifdef _HAVE_PETSC_ VecCopy(this->vector,to->vector); #else - _error_("not implemented yet!"); + this->vector->Copy(to->vector); #endif } @@ -305,7 +309,7 @@ _assert_(this->vector); VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm); #else - _error_("not implemented yet!"); + norm=this->vector->Norm(norm_type); #endif return norm; } @@ -317,7 +321,7 @@ _assert_(this->vector); VecScale(this->vector,scale_factor); #else - _error_("not implemented yet!"); + this->vector->Scale(scale_factor); #endif } /*}}}*/ @@ -329,7 +333,7 @@ _assert_(this->vector); VecDot(this->vector,vector->vector,&dot); #else - _error_("not implemented yet!"); + dot=this->vector->Dot(vector->vector); #endif return dot; } @@ -341,7 +345,7 @@ _assert_(this->vector); VecPointwiseDivide(this->vector,x->vector,y->vector); #else - _error_("not implemented yet!"); + this->vector->PointwiseDivide(x->vector,y->vector); #endif } /*}}}*/ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.cpp (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.cpp (revision 11734) @@ -47,7 +47,7 @@ #ifdef _HAVE_PETSC_ this->matrix=NewMat(pM,pN); #else - this->matrix=(double*)xcalloc(pM*pN,sizeof(double)); + this->matrix=new SeqMat(pM,pN); #endif #ifdef _HAVE_ADOLC_ this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble)); @@ -63,7 +63,7 @@ #ifdef _HAVE_PETSC_ this->matrix=NewMat(pM,pN,sparsity); #else - this->matrix=(double*)xcalloc(pM*pN,sizeof(double)); + this->matrix=new SeqMat(pM,pN,sparsity); #endif #ifdef _HAVE_ADOLC_ this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble)); @@ -96,7 +96,7 @@ xfree((void**)&idxm); xfree((void**)&idxn); #else - this->matrix=(double*)xcalloc(pM*pN,sizeof(double)); + this->matrix=new SeqMat(serial_mat,pM,pN,sparsity); #endif #ifdef _HAVE_ADOLC_ this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble)); @@ -112,7 +112,7 @@ #ifdef _HAVE_PETSC_ this->matrix=NewMat(pM,pN,connectivity,numberofdofspernode); #else - this->matrix=(double*)xcalloc(pM*pN,sizeof(double)); + this->matrix=new SeqMat(pM,pN,connectivity,numberofdofspernode); #endif #ifdef _HAVE_ADOLC_ this->amatrix=(adouble*)xmalloc(pM*pN*sizeof(adouble)); @@ -125,7 +125,7 @@ #ifdef _HAVE_PETSC_ MatFree(&this->matrix); #else - xfree((void**)&this->matrix); + delete this->matrix; #endif #ifdef _HAVE_ADOLC_ xfree((void**)&this->amatrix); @@ -142,13 +142,7 @@ #ifdef _HAVE_PETSC_ MatView(this->matrix,PETSC_VIEWER_STDOUT_WORLD); #else - printf("Matrix size: %i-%i\n",M,N); - for(i=0;imatrix->Echo(); #endif #ifdef _HAVE_ADOLC_ @@ -172,12 +166,30 @@ #ifdef _HAVE_PETSC_ PetscMatrixToMatlabMatrix(&dataref,this->matrix); #else - _error_("not implemented yet!"); + dataref=this->matrix->ToMatlabMatrix(); #endif return dataref; } /*}}}*/ +/*FUNCTION MatlabMatrixToMatrix{{{1*/ +Matrix* MatlabMatrixToMatrix(const mxArray* mxmatrix){ + + int dummy; + Matrix* matrix=NULL; + + /*allocate matrix object: */ + matrix=new Matrix(); + + #ifdef _HAVE_PETSC_ + MatlabMatrixToPetscMatrix(&matrix->matrix,&matrix->M,&matrix->N,mxmatrix); + #else + matrix->matrix=MatlabMatrixToSeqMat(mxmatrix); + #endif + + return matrix; +} +/*}}}*/ #endif /*FUNCTION Matrix::Assemble{{{1*/ void Matrix::Assemble(void){ @@ -189,7 +201,7 @@ MatCompress(this->matrix); #endif #else - /*do nothing:*/ + this->matrix->Assemble(); #endif } @@ -202,7 +214,7 @@ _assert_(this->matrix); MatNorm(this->matrix,ISSMToPetscNormMode(norm_type),&norm); #else - _error_("not implemented yet!"); + norm=this->matrix->Norm(norm_type); #endif return norm; } @@ -214,7 +226,7 @@ _assert_(this->matrix); MatGetSize(this->matrix,pM,pN); #else - _error_("not implemented yet!"); + this->matrix->GetSize(pM,pN); #endif } /*}}}*/ @@ -225,7 +237,7 @@ _assert_(this->matrix); MatGetLocalSize(this->matrix,pM,pN); #else - _error_("not implemented yet!"); + this->matrix->GetLocalSize(pM,pN); #endif } /*}}}*/ @@ -237,7 +249,7 @@ _assert_(X->vector); MatMultPatch(this->matrix,X->vector,AX->vector); #else - _error_("not implemented yet!"); + this->matrix->MatMult(X->vector,AX->vector); #endif } /*}}}*/ @@ -248,11 +260,11 @@ output=new Matrix(); - #ifdef _HAVE_PETSC + #ifdef _HAVE_PETSC_ _assert_(this->matrix); MatDuplicate(this->matrix,MAT_COPY_VALUES,&output->matrix); #else - _error_("not implemented yet!"); + output->matrix=this->matrix->Duplicate(); #endif } /*}}}*/ @@ -261,11 +273,21 @@ double* output=NULL; - #ifdef _HAVE_PETSC + #ifdef _HAVE_PETSC_ MatToSerial(&output,this->matrix); #else - _error_("not implemented yet!"); + output=this->matrix->ToSerial(); #endif return output; } /*}}}*/ +/*FUNCTION Matrix::SetValues{{{1*/ +void Matrix::SetValues(int m,int* idxm,int n,int* idxn,double* 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 +} +/*}}}*/ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.h (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Vector.h (revision 11734) @@ -30,12 +30,10 @@ public: - int M; - #ifdef _HAVE_PETSC_ Vec vector; #else - double* vector; + SeqVec* vector; #endif #ifdef _HAVE_ADOLC_ adouble* avector; @@ -73,4 +71,10 @@ double Dot(Vector* vector); /*}}}*/ }; + +/*API: */ +#ifdef _SERIAL_ +Vector* MatlabVectorToVector(const mxArray* mxvector); +#endif + #endif //#ifndef _VECTOR_H_ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.h =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.h (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/Matrix.h (revision 11734) @@ -36,7 +36,7 @@ #ifdef _HAVE_PETSC_ Mat matrix; #else - double* matrix; + SeqMat* matrix; #endif #ifdef _HAVE_ADOLC_ adouble* amatrix; @@ -62,7 +62,13 @@ void MatMult(Vector* X,Vector* AX); Matrix* Duplicate(void); double* ToSerial(void); + void SetValues(int m,int* idxm,int n,int* idxn,double* values,InsMode mode); /*}}}*/ }; +/*API: */ +#ifdef _SERIAL_ +Matrix* MatlabMatrixToMatrix(const mxArray* mxmatrix); +#endif + #endif //#ifndef _MATRIX_H_ Index: /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp =================================================================== --- /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp (revision 11733) +++ /proj/ice/larour/issm-uci-clean/trunk-jpl/src/c/objects/Numerics/ElementVector.cpp (revision 11734) @@ -165,22 +165,18 @@ int i; double* localvalues=NULL; - #ifdef _HAVE_PETSC_ - if(this->fsize){ - /*first, retrieve values that are in the f-set from the g-set values vector: */ - localvalues=(double*)xmalloc(this->fsize*sizeof(double)); - for(i=0;ifsize;i++){ - localvalues[i]=this->values[this->flocaldoflist[i]]; - } - /*add local values into global vector, using the fglobaldoflist: */ - VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,ADD_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); + if(this->fsize){ + /*first, retrieve values that are in the f-set from the g-set values vector: */ + localvalues=(double*)xmalloc(this->fsize*sizeof(double)); + for(i=0;ifsize;i++){ + localvalues[i]=this->values[this->flocaldoflist[i]]; } - #else - _error_("not supported yet!"); - #endif + /*add local values into global vector, using the fglobaldoflist: */ + pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL); + + /*Free ressources:*/ + xfree((void**)&localvalues); + } } /*}}}*/ @@ -190,23 +186,19 @@ int i; double* localvalues=NULL; - #ifdef _HAVE_PETSC_ - if(this->fsize){ - /*first, retrieve values that are in the f-set from the g-set values vector: */ - localvalues=(double*)xmalloc(this->fsize*sizeof(double)); - for(i=0;ifsize;i++){ - localvalues[i]=this->values[this->flocaldoflist[i]]; - } - /*add local values into global vector, using the fglobaldoflist: */ - VecSetValues(pf->vector,this->fsize,this->fglobaldoflist,(const double*)localvalues,INSERT_VALUES); - - /*Free ressources:*/ - xfree((void**)&localvalues); + if(this->fsize){ + /*first, retrieve values that are in the f-set from the g-set values vector: */ + localvalues=(double*)xmalloc(this->fsize*sizeof(double)); + for(i=0;ifsize;i++){ + localvalues[i]=this->values[this->flocaldoflist[i]]; } - #else - _error_("not supported yet!"); - #endif + /*add local values into global vector, using the fglobaldoflist: */ + pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,INS_VAL); + /*Free ressources:*/ + xfree((void**)&localvalues); + } + } /*}}}*/ /*FUNCTION ElementVector::Echo{{{1*/