Changeset 25995


Ignore:
Timestamp:
02/17/21 09:25:47 (4 years ago)
Author:
bulthuis
Message:

CGH: adding elementwise multiplication and power for vectors

Location:
issm/trunk-jpl/src/c/toolkits
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h

    r23643 r25995  
    66 *        IssmSeqVec and IssmMpiVec
    77 *
    8  */ 
     8 */
    99
    1010#ifndef _ISSM_ABS_VEC_H_
    … …  
    1616
    1717/*We need to template this class, in case we want to create Vectors that hold
    18   IssmDouble* vector or IssmPDouble* vector. 
     18  IssmDouble* vector or IssmPDouble* vector.
    1919  Such vectors are useful for use without or with the matlab or python
    2020  interface (which do not care for IssmDouble types, but only rely on
    2121  IssmPDouble types)
    2222*/
    23 template <class doubletype> 
     23template <class doubletype>
    2424class IssmAbsVec{
    2525
    … …  
    5050                virtual doubletype Dot(IssmAbsVec* input)=0;
    5151                virtual void PointwiseDivide(IssmAbsVec* x,IssmAbsVec* y)=0;
     52                virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0;
     53                virtual void Pow(doubletype scale_factor)=0;
    5254};
    5355

  • issm/trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h

    r25336 r25995  
    11/*!\file:  IssmMpiVec.h
    2  * \brief implementation of parallel dense ISSM vector. Internally, the parallel dense vector is 
    3  * split in rows across each cpu. Each vector (representing a subset of rows) on each cpu is fully 
    4  * dense, and is represented by a linear buffer of type doubletype. 
    5  * This object needs to answer the API defined by the virtual functions in IssmAbsVec, 
     2 * \brief implementation of parallel dense ISSM vector. Internally, the parallel dense vector is
     3 * split in rows across each cpu. Each vector (representing a subset of rows) on each cpu is fully
     4 * dense, and is represented by a linear buffer of type doubletype.
     5 * This object needs to answer the API defined by the virtual functions in IssmAbsVec,
    66 * and the contructors required by IssmVec (see IssmVec.h)
    7  */ 
     7 */
    88
    99#ifndef _ISSM_MPI_VEC_H_
    … …  
    2626/*}}}*/
    2727
    28 /*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. 
    29   Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, 
     28/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector.
     29  Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types,
    3030  but only rely on IssmPDouble types)*/
    3131template <class doubletype> class IssmAbsVec;
    3232
    33 template <class doubletype> 
     33template <class doubletype>
    3434class IssmMpiVec:public IssmAbsVec<doubletype>{
    3535
    … …  
    206206
    207207                        /*Recap, each cpu has num_procs datasets of buckets. For a certain cpu j, for a given dataset i, the buckets this  {{{
    208                          * 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 
     208                         * 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
    209209                         * vectors that will be shipped around the cluster: */
    210210                        this->BucketsBuildScatterBuffers(&numvalues_forcpu,&row_indices_forcpu,&values_forcpu,&modes_forcpu,bucketsforcpu,num_procs);
    … …  
    244244
    245245                        /*Scatter values around: {{{*/
    246                         /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given 
    247                          * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the ISSM_MPI_Scatterv prototype: 
     246                        /*Now, to scatter values across the cluster, we need sendcnts and displs. Our sendbufs have been built by BucketsBuildScatterBuffers, with a stride given
     247                         * by numvalues_forcpu. Get this ready to go before starting the scatter itslef. For reference, here is the ISSM_MPI_Scatterv prototype:
    248248                         * 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) :*/
    249249                        sendcnts=xNew<int>(num_procs);
    … …  
    308308                        xDelete<int>(sendcnts);
    309309                        xDelete<int>(displs);
    310                        
     310
    311311                        /*Get rid of all buckets, as we have already added them to the matrix!: */
    312312                        delete this->buckets;
    … …  
    318318                void SetValues(int ssize, int* list, doubletype* values, InsMode mode){/*{{{*/
    319319
    320                         /*we need to store all the values we collect here in order to Assemble later. 
    321                          * Indeed, the values we are collecting here most of the time will not belong 
     320                        /*we need to store all the values we collect here in order to Assemble later.
     321                         * Indeed, the values we are collecting here most of the time will not belong
    322322                         * to us, but to another part of the vector on another cpu: */
    323323                        _assert_(buckets);
    … …  
    329329                void SetValue(int dof, doubletype value, InsMode mode){/*{{{*/
    330330
    331                         /*we need to store the value we collect here in order to Assemble later. 
    332                          * Indeed, the value we are collecting here most of the time will not belong 
     331                        /*we need to store the value we collect here in order to Assemble later.
     332                         * Indeed, the value we are collecting here most of the time will not belong
    333333                         * to us, but to another part of the vector on another cpu: */
    334334                        _assert_(buckets);
    … …  
    365365
    366366                        /*Get Ownership range*/
    367                         int lower_row,upper_row; 
     367                        int lower_row,upper_row;
    368368                        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    369                         int range=upper_row-lower_row;   
     369                        int range=upper_row-lower_row;
    370370
    371371                        /*return NULL if no range*/
    … …  
    377377
    378378                        /*Build indices*/
    379                         int* indices=xNew<int>(range); 
     379                        int* indices=xNew<int>(range);
    380380                        for(int i=0;i<range;i++) indices[i]=lower_row+i;
    381381
    … …  
    518518                                        break;
    519519                                case NORM_TWO:
    520                                         local_norm=0.; 
     520                                        local_norm=0.;
    521521                                        for(i=0;i<this->m;i++)local_norm+=this->vector[i]*this->vector[i];
    522522                                        ISSM_MPI_Reduce(&local_norm, &norm, 1, ISSM_MPI_DOUBLE, ISSM_MPI_SUM, 0, IssmComm::GetComm());
    … …  
    570570                        /*pointwise w=x/y where this->vector is w: */
    571571                        for(i=0;i<this->m;i++)this->vector[i]=x->vector[i]/y->vector[i];
     572                }
     573                /*}}}*/
     574                void PointwiseMult(IssmAbsVec<doubletype>* xin,IssmAbsVec<doubletype>* yin){/*{{{*/
     575
     576                        int i;
     577
     578                        /*Assume xin and yin are of the correct type, and downcast: */
     579                        IssmMpiVec* x=NULL;
     580                        IssmMpiVec* y=NULL;
     581
     582                        x=(IssmMpiVec<doubletype>*)xin;
     583                        y=(IssmMpiVec<doubletype>*)yin;
     584
     585                        /*pointwise w=x*y where this->vector is w: */
     586                        for(i=0;i<this->m;i++)this->vector[i]=x->vector[i]*y->vector[i];
     587                }
     588                /*}}}*/
     589                void Pow(doubletype scale_factor){/*{{{*/
     590
     591                        int i;
     592                        for(i=0;i<this->M;i++)this->vector[i]=pow(this->vector[i],scale_factor);
     593
    572594                }
    573595                /*}}}*/
    … …  
    623645
    624646                        /*we are going to march through the buffers, and marshall data onto them, so in order to not
    625                          *lose track of where these buffers are located in memory, we are going to work using copies 
     647                         *lose track of where these buffers are located in memory, we are going to work using copies
    626648                         of them: */
    627649                        temp_row_indices_forcpu=row_indices_forcpu;
    … …  
    649671                        *pmodes_forcpu       = modes_forcpu;
    650672                }
    651                 /*}}}*/         
     673                /*}}}*/
    652674};
    653 #endif //#ifndef _ISSM_MPI_VEC_H_       
     675#endif //#ifndef _ISSM_MPI_VEC_H_

  • issm/trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h

    r24682 r25995  
    11/*!\file:  IssmSeqVec.h
    2  * \brief implementation of an ISSM vector which run serially (1 cpu only), which is made of a fully dense 
    3  * vector. Internally, this dense vector is just a linear buffer of type doubletype. 
    4  * This object needs to answer the API defined by the virtual functions in IssmAbsVec, 
     2 * \brief implementation of an ISSM vector which run serially (1 cpu only), which is made of a fully dense
     3 * vector. Internally, this dense vector is just a linear buffer of type doubletype.
     4 * This object needs to answer the API defined by the virtual functions in IssmAbsVec,
    55 * and the contructors required by IssmVec (see IssmVec.h)
    6  */ 
     6 */
    77
    88#ifndef _ISSM_SEQ_VEC_H_
    … …  
    2222/*}}}*/
    2323
    24 /*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector. 
    25   Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types, 
     24/*We need to template this class, in case we want to create vectors that hold IssmDouble* vector or IssmPDouble* vector.
     25  Such vectors would be useful for use without or with the matlab or python interface (which do not care for IssmDouble types,
    2626  but only rely on IssmPDouble types)*/
    2727
    2828template <class doubletype> class IssmAbsVec;
    2929
    30 template <class doubletype> 
     30template <class doubletype>
    3131class IssmSeqVec: public IssmAbsVec<doubletype>{
    3232
    … …  
    155155
    156156                        /*Build indices*/
    157                         int* indices=xNew<int>(vector_size); 
     157                        int* indices=xNew<int>(vector_size);
    158158                        for(int i=0;i<vector_size;i++) indices[i]=i;
    159159
    … …  
    255255                                        break;
    256256                                case NORM_TWO:
    257                                         norm=0.; 
     257                                        norm=0.;
    258258                                        for(i=0;i<this->M;i++)norm+=this->vector[i]*this->vector[i];
    259259                                        return sqrt(norm);
    … …  
    302302                }
    303303                /*}}}*/
     304                void PointwiseMult(IssmAbsVec<doubletype>* xin,IssmAbsVec<doubletype>* yin){/*{{{*/
     305
     306                        int i;
     307
     308                        /*Assume xin and yin are of the correct type, and downcast: */
     309                        IssmSeqVec* x=NULL;
     310                        IssmSeqVec* y=NULL;
     311
     312                        x=(IssmSeqVec<doubletype>*)xin;
     313                        y=(IssmSeqVec<doubletype>*)yin;
     314
     315                        /*pointwise w=x*y where this->vector is w: */
     316                        for(i=0;i<this->M;i++)this->vector[i]=x->vector[i]*y->vector[i];
     317                }
     318                /*}}}*/
     319                void Pow(doubletype scale_factor){/*{{{*/
     320
     321                        int i;
     322                        for(i=0;i<this->M;i++)this->vector[i]=pow(this->vector[i],scale_factor);
     323
     324                }
     325                /*}}}*/
     326               
    304327};
    305328#endif //#ifndef _ISSM_SEQ_VEC_H_

  • issm/trunk-jpl/src/c/toolkits/issm/IssmVec.h

    r25390 r25995  
    11/*!\file:  IssmVec.h
    2  * \brief Main Vector class for the Issm toolkit. 
    3  */ 
     2 * \brief Main Vector class for the Issm toolkit.
     3 */
    44
    55#ifndef _ISSMVEC_H_
    … …  
    2222
    2323/*We need to template this class, in case we want to create Vectors that hold
    24   IssmDouble* matrix or IssmPDouble* matrix. 
     24  IssmDouble* matrix or IssmPDouble* matrix.
    2525  Such vectors are useful for use without or with the matlab or python
    2626  interface (which do not care for IssmDouble types, but only rely on
    … …  
    3131template <class doubletype> class IssmMpiVec;
    3232
    33 template <class doubletype> 
     33template <class doubletype>
    3434class IssmVec{
    3535
    … …  
    4747                        switch(IssmVecTypeFromToolkitOptions()){
    4848
    49                                 case SeqEnum: 
     49                                case SeqEnum:
    5050                                        this->vector=new IssmSeqVec<doubletype>(M);
    5151                                        break;
    … …  
    6666                        switch(IssmVecTypeFromToolkitOptions()){
    6767
    68                                 case SeqEnum: 
     68                                case SeqEnum:
    6969                                        this->vector=new IssmSeqVec<doubletype>(m,M);
    7070                                        break;
    … …  
    8585                        switch(IssmVecTypeFromToolkitOptions()){
    8686
    87                                 case SeqEnum: 
     87                                case SeqEnum:
    8888                                        this->vector=new IssmSeqVec<doubletype>(M,fromlocalsize);
    8989                                        break;
    … …  
    104104                        switch(IssmVecTypeFromToolkitOptions()){
    105105
    106                                 case SeqEnum: 
     106                                case SeqEnum:
    107107                                        this->vector=new IssmSeqVec<doubletype>(buffer,M);
    108108                                        break;
    … …  
    212212                }
    213213                /*}}}*/
     214                void PointwiseMult(IssmVec* x,IssmVec* y){/*{{{*/
     215                        vector->PointwiseMult(x->vector,y->vector);
     216                }
     217                /*}}}*/
     218                void Pow(doubletype scale_factor){/*{{{*/
     219                        vector->Pow(scale_factor);
     220                }
     221                /*}}}*/
    214222};
    215223

  • issm/trunk-jpl/src/c/toolkits/objects/Vector.h

    r23643 r25995  
    11/*!\file:  Vector.h
    2  * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?) 
     2 * \brief wrapper to vector objects. The goal is to control which API (PETSc,Scalpack, Plapack?)
    33 * implements our underlying vector format.
    4  */ 
     4 */
    55
    66#ifndef _VECTOR_H_
    … …  
    2222enum vectortype { PetscVecType, IssmVecType };
    2323
    24 template <class doubletype> 
     24template <class doubletype>
    2525class Vector{
    2626
    … …  
    3131                PetscVec* pvector;
    3232                #endif
    33                 IssmVec<doubletype>* ivector; 
     33                IssmVec<doubletype>* ivector;
    3434
    3535                /*Vector constructors, destructors */
    … …  
    105105                        /*retrieve toolkittype: */
    106106                        char* toolkittype=ToolkitOptions::GetToolkitType();
    107                         _assert_(toolkittype); 
     107                        _assert_(toolkittype);
    108108
    109109                        /*set vector type: */
    110110                        if(strcmp(toolkittype,"petsc")==0){
    111111                                #ifdef _HAVE_PETSC_
    112                                 type=PetscVecType; 
     112                                type=PetscVecType;
    113113                                #else
    114114                                _error_("cannot create petsc vector without PETSC compiled!");
    … …  
    200200                        this->GetSize(&M);
    201201
    202                         if(M==0) 
     202                        if(M==0)
    203203                                return true;
    204204                        else
    … …  
    386386                }
    387387                /*}}}*/
     388                void PointwiseMult(Vector* x,Vector* y){_assert_(this);/*{{{*/
     389
     390                        if(type==PetscVecType){
     391                                #ifdef _HAVE_PETSC_
     392                                this->pvector->PointwiseMult(x->pvector,y->pvector);
     393                                #endif
     394                        }
     395                        else this->ivector->PointwiseMult(x->ivector,y->ivector);
     396                }
     397                /*}}}*/
     398                void Pow(doubletype scale_factor){_assert_(this);/*{{{*/
     399
     400                        if(type==PetscVecType){
     401                                #ifdef _HAVE_PETSC_
     402                                this->pvector->Pow(scale_factor);
     403                                #endif
     404                        }
     405                        else this->ivector->Pow(scale_factor);
     406                }
     407                /*}}}*/
    388408};
    389409#endif //#ifndef _VECTOR_H_
Note: See TracChangeset for help on using the changeset viewer.