Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 25994) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 25995) @@ -6,5 +6,5 @@ * IssmSeqVec and IssmMpiVec * - */ + */ #ifndef _ISSM_ABS_VEC_H_ @@ -16,10 +16,10 @@ /*We need to template this class, in case we want to create Vectors that hold - IssmDouble* vector or IssmPDouble* vector. + IssmDouble* vector or IssmPDouble* vector. Such vectors are useful for use without or with the matlab or python interface (which do not care for IssmDouble types, but only rely on IssmPDouble types) */ -template +template class IssmAbsVec{ @@ -50,4 +50,6 @@ virtual doubletype Dot(IssmAbsVec* input)=0; virtual void PointwiseDivide(IssmAbsVec* x,IssmAbsVec* y)=0; + virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0; + virtual void Pow(doubletype scale_factor)=0; }; Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 25994) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 25995) @@ -1,9 +1,9 @@ /*!\file: IssmMpiVec.h - * \brief implementation of parallel dense ISSM vector. Internally, the parallel dense vector is - * split in rows across each cpu. Each vector (representing a subset of rows) on each cpu is fully - * dense, and is represented by a linear buffer of type doubletype. - * This object needs to answer the API defined by the virtual functions in IssmAbsVec, + * \brief implementation of parallel dense ISSM vector. Internally, the parallel dense vector is + * split in rows across each cpu. Each vector (representing a subset of rows) on each cpu is fully + * dense, and is represented by a linear buffer of type doubletype. + * This object needs to answer the API defined by the virtual functions in IssmAbsVec, * and the contructors required by IssmVec (see IssmVec.h) - */ + */ #ifndef _ISSM_MPI_VEC_H_ @@ -26,10 +26,10 @@ /*}}}*/ -/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. - Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, +/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. + Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, but only rely on IssmPDouble types)*/ template class IssmAbsVec; -template +template class IssmMpiVec:public IssmAbsVec{ @@ -206,5 +206,5 @@ /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this {{{ - * dataset owns correspond to rows that are owned by cpu i, not j!. Out of all the buckets we own, make row,col,value,insert_mode + * dataset owns correspond to rows that are owned by cpu i, not j!. Out of all the buckets we own, make row,col,value,insert_mode * vectors that will be shipped around the cluster: */ this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs); @@ -244,6 +244,6 @@ /*Scatter values around: {{{*/ - /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given - * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the ISSM_MPI_Scatterv prototype: + /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given + * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the ISSM_MPI_Scatterv prototype: * int ISSM_MPI_Scatterv( void *sendbuf, int *sendcnts, int *displs, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcnt, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm) :*/ sendcnts=xNew(num_procs); @@ -308,5 +308,5 @@ xDelete(sendcnts); xDelete(displs); - + /*Get rid of all buckets, as we have already added them to the matrix!: */ delete this->buckets; @@ -318,6 +318,6 @@ void SetValues(int ssize, int* list, doubletype* values, InsMode mode){/*{{{*/ - /*we need to store all the values we collect here in order to Assemble later. - * Indeed, the values we are collecting here most of the time will not belong + /*we need to store all the values we collect here in order to Assemble later. + * Indeed, the values we are collecting here most of the time will not belong * to us, but to another part of the vector on another cpu: */ _assert_(buckets); @@ -329,6 +329,6 @@ void SetValue(int dof, doubletype value, InsMode mode){/*{{{*/ - /*we need to store the value we collect here in order to Assemble later. - * Indeed, the value we are collecting here most of the time will not belong + /*we need to store the value we collect here in order to Assemble later. + * Indeed, the value we are collecting here most of the time will not belong * to us, but to another part of the vector on another cpu: */ _assert_(buckets); @@ -365,7 +365,7 @@ /*Get Ownership range*/ - int lower_row,upper_row; + int lower_row,upper_row; GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - int range=upper_row-lower_row; + int range=upper_row-lower_row; /*return NULL if no range*/ @@ -377,5 +377,5 @@ /*Build indices*/ - int* indices=xNew(range); + int* indices=xNew(range); for(int i=0;im;i++)local_norm+=this->vector[i]*this->vector[i]; ISSM_MPI_Reduce(&local_norm, &norm, 1, ISSM_MPI_DOUBLE, ISSM_MPI_SUM, 0, IssmComm::GetComm()); @@ -570,4 +570,26 @@ /*pointwise w=x/y where this->vector is w: */ for(i=0;im;i++)this->vector[i]=x->vector[i]/y->vector[i]; + } + /*}}}*/ + void PointwiseMult(IssmAbsVec* xin,IssmAbsVec* yin){/*{{{*/ + + int i; + + /*Assume xin and yin are of the correct type, and downcast: */ + IssmMpiVec* x=NULL; + IssmMpiVec* y=NULL; + + x=(IssmMpiVec*)xin; + y=(IssmMpiVec*)yin; + + /*pointwise w=x*y where this->vector is w: */ + for(i=0;im;i++)this->vector[i]=x->vector[i]*y->vector[i]; + } + /*}}}*/ + void Pow(doubletype scale_factor){/*{{{*/ + + int i; + for(i=0;iM;i++)this->vector[i]=pow(this->vector[i],scale_factor); + } /*}}}*/ @@ -623,5 +645,5 @@ /*we are going to march through the buffers, and marshall data onto them, so in order to not - *lose track of where these buffers are located in memory, we are going to work using copies + *lose track of where these buffers are located in memory, we are going to work using copies of them: */ temp_row_indices_forcpu=row_indices_forcpu; @@ -649,5 +671,5 @@ *pmodes_forcpu = modes_forcpu; } - /*}}}*/ + /*}}}*/ }; -#endif //#ifndef _ISSM_MPI_VEC_H_ +#endif //#ifndef _ISSM_MPI_VEC_H_ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 25994) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 25995) @@ -1,8 +1,8 @@ /*!\file: IssmSeqVec.h - * \brief implementation of an ISSM vector which run serially (1 cpu only), which is made of a fully dense - * vector. Internally, this dense vector is just a linear buffer of type doubletype. - * This object needs to answer the API defined by the virtual functions in IssmAbsVec, + * \brief implementation of an ISSM vector which run serially (1 cpu only), which is made of a fully dense + * vector. Internally, this dense vector is just a linear buffer of type doubletype. + * This object needs to answer the API defined by the virtual functions in IssmAbsVec, * and the contructors required by IssmVec (see IssmVec.h) - */ + */ #ifndef _ISSM_SEQ_VEC_H_ @@ -22,11 +22,11 @@ /*}}}*/ -/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. - Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, +/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. + Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, but only rely on IssmPDouble types)*/ template class IssmAbsVec; -template +template class IssmSeqVec: public IssmAbsVec{ @@ -155,5 +155,5 @@ /*Build indices*/ - int* indices=xNew(vector_size); + int* indices=xNew(vector_size); for(int i=0;iM;i++)norm+=this->vector[i]*this->vector[i]; return sqrt(norm); @@ -302,4 +302,27 @@ } /*}}}*/ + void PointwiseMult(IssmAbsVec* xin,IssmAbsVec* yin){/*{{{*/ + + int i; + + /*Assume xin and yin are of the correct type, and downcast: */ + IssmSeqVec* x=NULL; + IssmSeqVec* y=NULL; + + x=(IssmSeqVec*)xin; + y=(IssmSeqVec*)yin; + + /*pointwise w=x*y where this->vector is w: */ + for(i=0;iM;i++)this->vector[i]=x->vector[i]*y->vector[i]; + } + /*}}}*/ + void Pow(doubletype scale_factor){/*{{{*/ + + int i; + for(i=0;iM;i++)this->vector[i]=pow(this->vector[i],scale_factor); + + } + /*}}}*/ + }; #endif //#ifndef _ISSM_SEQ_VEC_H_ Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 25994) +++ /issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 25995) @@ -1,5 +1,5 @@ /*!\file: IssmVec.h - * \brief Main Vector class for the Issm toolkit. - */ + * \brief Main Vector class for the Issm toolkit. + */ #ifndef _ISSMVEC_H_ @@ -22,5 +22,5 @@ /*We need to template this class, in case we want to create Vectors that hold - IssmDouble* matrix or IssmPDouble* matrix. + IssmDouble* matrix or IssmPDouble* matrix. Such vectors are useful for use without or with the matlab or python interface (which do not care for IssmDouble types, but only rely on @@ -31,5 +31,5 @@ template class IssmMpiVec; -template +template class IssmVec{ @@ -47,5 +47,5 @@ switch(IssmVecTypeFromToolkitOptions()){ - case SeqEnum: + case SeqEnum: this->vector=new IssmSeqVec(M); break; @@ -66,5 +66,5 @@ switch(IssmVecTypeFromToolkitOptions()){ - case SeqEnum: + case SeqEnum: this->vector=new IssmSeqVec(m,M); break; @@ -85,5 +85,5 @@ switch(IssmVecTypeFromToolkitOptions()){ - case SeqEnum: + case SeqEnum: this->vector=new IssmSeqVec(M,fromlocalsize); break; @@ -104,5 +104,5 @@ switch(IssmVecTypeFromToolkitOptions()){ - case SeqEnum: + case SeqEnum: this->vector=new IssmSeqVec(buffer,M); break; @@ -212,4 +212,12 @@ } /*}}}*/ + void PointwiseMult(IssmVec* x,IssmVec* y){/*{{{*/ + vector->PointwiseMult(x->vector,y->vector); + } + /*}}}*/ + void Pow(doubletype scale_factor){/*{{{*/ + vector->Pow(scale_factor); + } + /*}}}*/ }; Index: /issm/trunk-jpl/src/c/toolkits/objects/Vector.h =================================================================== --- /issm/trunk-jpl/src/c/toolkits/objects/Vector.h (revision 25994) +++ /issm/trunk-jpl/src/c/toolkits/objects/Vector.h (revision 25995) @@ -1,6 +1,6 @@ /*!\file: Vector.h - * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) + * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) * implements our underlying vector format. - */ + */ #ifndef _VECTOR_H_ @@ -22,5 +22,5 @@ enum vectortype { PetscVecType, IssmVecType }; -template +template class Vector{ @@ -31,5 +31,5 @@ PetscVec* pvector; #endif - IssmVec* ivector; + IssmVec* ivector; /*Vector constructors, destructors */ @@ -105,10 +105,10 @@ /*retrieve toolkittype: */ char* toolkittype=ToolkitOptions::GetToolkitType(); - _assert_(toolkittype); + _assert_(toolkittype); /*set vector type: */ if(strcmp(toolkittype,"petsc")==0){ #ifdef _HAVE_PETSC_ - type=PetscVecType; + type=PetscVecType; #else _error_("cannot create petsc vector without PETSC compiled!"); @@ -200,5 +200,5 @@ this->GetSize(&M); - if(M==0) + if(M==0) return true; else @@ -386,4 +386,24 @@ } /*}}}*/ + void PointwiseMult(Vector* x,Vector* y){_assert_(this);/*{{{*/ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->PointwiseMult(x->pvector,y->pvector); + #endif + } + else this->ivector->PointwiseMult(x->ivector,y->ivector); + } + /*}}}*/ + void Pow(doubletype scale_factor){_assert_(this);/*{{{*/ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Pow(scale_factor); + #endif + } + else this->ivector->Pow(scale_factor); + } + /*}}}*/ }; #endif //#ifndef _VECTOR_H_