Changeset 11734 for issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp

Timestamp:

03/19/12 20:11:39 (13 years ago)

Author:

Eric.Larour

Message:

Hook up new issm toolkit

File:

• issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp (modified) (25 diffs)

issm/trunk-jpl/src/c/objects/Numerics/Vector.cpp

@@ -25 +25 @@
 Vector::Vector(){
 
-        this->M=0;
         #ifdef _HAVE_PETSC_
         this->vector=NULL;
@@ -42 +41 @@
         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_
@@ -51 +49 @@
 /*}}}*/
 /*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;
@@ -58 +56 @@
                 int* idxm=NULL;
 
-                this->vector=NewVec(pM);
+                this->vector=NewVec(M);
                 
                 idxm=(int*)xmalloc(M*sizeof(int));
@@ -70 +68 @@
 
         #else
-                this->M=pM;
-                this->vector=(double*)xcalloc(pM,sizeof(double));
-        #endif
-        #ifdef _HAVE_ADOLC_
-                this->avector=(adouble*)xmalloc(pM*sizeof(adouble));
+                this->vector=new SeqVec(serial_vec,M);
+        #endif
+        #ifdef _HAVE_ADOLC_
+                this->avector=(adouble*)xmalloc(M*sizeof(adouble));
         #endif
 }
@@ -82 +79 @@
 Vector::Vector(Vec petsc_vec){
 
-        /*Get Vector size*/
-        VecGetSize(petsc_vec,&this->M);
-
         /*copy vector*/
         VecDuplicate(petsc_vec,&this->vector);
@@ -98 +92 @@
         VecFree(&this->vector);
         #else
-        xfree((void**)&this->vector);
+        delete this->vector;
         #endif
         #ifdef _HAVE_ADOLC_
@@ -115 +109 @@
         VecView(this->vector,PETSC_VIEWER_STDOUT_WORLD);
         #else
-        printf("Vector size: %i\n",M);
-        for(i=0;i<M;i++){
-                printf("%g\n ",*(vector+i));
-        }
-        #endif
-
-        #ifdef _HAVE_ADOLC_
-        /*Not sure about that one. Should we use the overloaded operator >>?*/
-        printf("ADOLC Vector equivalent:" );
-        for(i=0;i<M;i++){
-                printf("%g\n ",*(avector+i));
-        }
+        this->vector->Echo();
+        #endif
+
+        #ifdef _HAVE_ADOLC_
+        /*do nothing for now: */
         #endif
 }
@@ -139 +126 @@
         PetscVectorToMatlabVector(&dataref,this->vector);
         #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;
 }
 /*}}}*/
@@ -154 +159 @@
                 VecAssemblyEnd(this->vector);
         #else
-                /*do nothing*/
+                this->vector->Assemble();
         #endif
 
@@ -167 +172 @@
                 VecSetValues(this->vector,ssize,list,values,ISSMToPetscInsertMode(mode));
         #else
-                _error_("not implemented yet!");
+                this->vector->SetValues(ssize,list,values,mode);
         #endif
 
@@ -179 +184 @@
                 VecSetValues(this->vector,1,&dof,&value,ISSMToPetscInsertMode(mode));
         #else
-                _error_("not implemented yet!");
+                this->vector->SetValue(dof,value,mode);
         #endif
 
@@ -191 +196 @@
                 VecGetValues(this->vector,1,&dof,pvalue);
         #else
-                _error_("not implemented yet!");
+        this->vector->GetValue(pvalue,dof);
         #endif
 }
@@ -202 +207 @@
                 VecGetSize(this->vector,pM);
         #else
-                _error_("not implemented yet!");
+                this->vector->GetSize(pM);
         #endif
 
@@ -214 +219 @@
                 VecGetLocalSize(this->vector,pM);
         #else
-                _error_("not implemented yet!");
+                this->vector->GetLocalSize(pM);
         #endif
 
@@ -229 +234 @@
                 VecDuplicate(this->vector,&vec_output);
                 output->vector=vec_output;
-                VecGetSize(output->vector,&output->M);
-        #else
-                _error_("not implemented yet!");
+        #else
+                output->vector=this->vector->Duplicate();
         #endif
         return output;
@@ -244 +248 @@
                 VecSet(this->vector,value);
         #else
-                _error_("not implemented yet!");
+                this->vector->Set(value);
         #endif
 
@@ -256 +260 @@
                 VecAXPY(this->vector,a,X->vector);
         #else
-                _error_("not implemented yet!");
+                this->vector->AXPY(X->vector,a);
         #endif
 }
@@ -267 +271 @@
                 VecAYPX(this->vector,a,X->vector);
         #else
-                _error_("not implemented yet!");
+                this->vector->AYPX(X->vector,a);
         #endif
 
@@ -280 +284 @@
                 VecToMPISerial(&vec_serial, this->vector);
         #else
-                _error_("not implemented yet!");
+                vec_serial=this->vector->ToMPISerial();
         #endif
 
@@ -293 +297 @@
                 VecCopy(this->vector,to->vector);
         #else
-                _error_("not implemented yet!");
+                this->vector->Copy(to->vector);
         #endif
 
@@ -306 +310 @@
                 VecNorm(this->vector,ISSMToPetscNormMode(norm_type),&norm);
         #else
-                _error_("not implemented yet!");
+                norm=this->vector->Norm(norm_type);
         #endif
         return norm;
@@ -318 +322 @@
                 VecScale(this->vector,scale_factor); 
         #else
-                _error_("not implemented yet!");
+                this->vector->Scale(scale_factor);
         #endif
 }
@@ -330 +334 @@
                 VecDot(this->vector,vector->vector,&dot);
         #else
-                _error_("not implemented yet!");
+                dot=this->vector->Dot(vector->vector);
         #endif
         return dot;
@@ -342 +346 @@
                 VecPointwiseDivide(this->vector,x->vector,y->vector);
         #else
-                _error_("not implemented yet!");
-        #endif
-}
-/*}}}*/
+                this->vector->PointwiseDivide(x->vector,y->vector);
+        #endif
+}
+/*}}}*/