IssmSeqVec.h

Go to the documentation of this file.

 
#ifndef _ISSM_SEQ_VEC_H_
#define _ISSM_SEQ_VEC_H_
 
/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
   #include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif
 
#include "../../shared/shared.h"
#include <math.h>
 
/*}}}*/
 
/*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 doubletype> class IssmAbsVec;
 
template <class doubletype> 
class IssmSeqVec: public IssmAbsVec<doubletype>{
 
   public:
 
      doubletype* vector;
      int M;
 
      /*IssmSeqVec constructors, destructors*/
      IssmSeqVec(){/*{{{*/
 
         this->M=0;
         this->vector=NULL;
      }
      /*}}}*/
      IssmSeqVec(int pM){/*{{{*/
 
         this->M=pM;
         this->vector=NULL;
         if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
      }
      /*}}}*/
      IssmSeqVec(int pm,int pM){/*{{{*/
 
         this->M=pM;
         this->vector=NULL;
         if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
      }
      /*}}}*/
      IssmSeqVec(int pM,bool fromlocalsize){/*{{{*/
 
         this->M=pM;
         this->vector=NULL;
         if(this->M) this->vector=xNewZeroInit<doubletype>(pM);
      }
      /*}}}*/
      IssmSeqVec(doubletype* buffer,int pM){/*{{{*/
 
         this->M=pM;
         this->vector=NULL;
         if(this->M){
            this->vector=xNew<doubletype>(pM);
            xMemCpy<doubletype>(this->vector,buffer,pM);
         }
      }
      /*}}}*/
      ~IssmSeqVec(){/*{{{*/
         if(this->M)xDelete<doubletype>(this->vector);
         M=0;
      }
      /*}}}*/
 
      /*IssmSeqVec specific routines*/
      void Echo(void){/*{{{*/
 
         int i;
         _printf_("IssmSeqVec size " << this->M << "\n");
         for(i=0;i<M;i++){
            _printf_(vector[i] << "\n ");
         }
      }
      /*}}}*/
      void Assemble(void){/*{{{*/
 
         /*do nothing*/
 
      }
      /*}}}*/
      void SetValues(int ssize, int* list, doubletype* values, InsMode mode){/*{{{*/
 
         int i;
         switch(mode){
            case ADD_VAL:
               for(i=0;i<ssize;i++) this->vector[list[i]]+=values[i];
               break;
            case INS_VAL:
               for(i=0;i<ssize;i++) this->vector[list[i]]=values[i];
               break;
            default:
               _error_("unknown insert mode!");
               break;
         }
 
      }
      /*}}}*/
      void SetValue(int dof, doubletype value, InsMode mode){/*{{{*/
 
         switch(mode){
            case ADD_VAL:
               this->vector[dof]+=value;
               break;
            case INS_VAL:
               this->vector[dof]=value;
               break;
            default:
               _error_("unknown insert mode!");
               break;
         }
      }
      /*}}}*/
      void GetValue(doubletype* pvalue,int dof){/*{{{*/
 
         *pvalue=this->vector[dof];
 
      }
      /*}}}*/
      void GetSize(int* pM){/*{{{*/
 
         *pM=this->M;
 
      }
      /*}}}*/
      void GetLocalSize(int* pM){/*{{{*/
         *pM=this->M;
      }/*}}}*/
      void GetLocalVector(doubletype** pvector,int** pindices){/*{{{*/
 
         /*First, check that vector size is not 0*/
         int vector_size;
         this->GetSize(&vector_size);
         if(vector_size==0){
            *pvector=NULL;
            *pindices=NULL;
            return;
         }
 
         /*Build indices*/
         int* indices=xNew<int>(vector_size); 
         for(int i=0;i<vector_size;i++) indices[i]=i;
 
         /*Get vector*/
         doubletype* values =xNew<doubletype>(vector_size);
         xMemCpy<doubletype>(values,this->vector,vector_size);
 
         *pvector  = values;
         *pindices = indices;
      } /*}}}*/
      IssmSeqVec<doubletype>* Duplicate(void){/*{{{*/
 
         return new IssmSeqVec<doubletype>(this->vector,this->M);
 
      }
      /*}}}*/
      void Set(doubletype value){/*{{{*/
 
         for(int i=0;i<this->M;i++)this->vector[i]=value;
 
      }
      /*}}}*/
      void AXPY(IssmAbsVec<doubletype>* Xin, doubletype a){/*{{{*/
 
         int i;
 
         /*Assume X is of the correct type, and downcast: */
         IssmSeqVec* X=NULL;
 
         X=(IssmSeqVec<doubletype>*)Xin;
 
         /*y=a*x+y where this->vector is y*/
         for(i=0;i<this->M;i++)this->vector[i]=a*X->vector[i]+this->vector[i];
 
      }
      /*}}}*/
      void AYPX(IssmAbsVec<doubletype>* Xin, doubletype a){/*{{{*/
 
         int i;
 
         /*Assume X is of the correct type, and downcast: */
         IssmSeqVec* X=NULL;
 
         X=(IssmSeqVec<doubletype>*)Xin;
 
         /*y=x+a*y where this->vector is y*/
         for(i=0;i<this->M;i++)this->vector[i]=X->vector[i]+a*this->vector[i];
 
      }
      /*}}}*/
      doubletype* ToMPISerial(void){/*{{{*/
 
         doubletype* buffer=NULL;
 
         if(this->M){
            buffer=xNew<doubletype>(this->M);
            xMemCpy<doubletype>(buffer,this->vector,this->M);
         }
         return buffer;
 
      }
      /*}}}*/
      doubletype* ToMPISerial0(void){/*{{{*/
 
         return this->ToMPISerial();
 
      }
      /*}}}*/
      void Shift(doubletype shift){/*{{{*/
 
         for(int i=0;i<this->M;i++)this->vector[i]+=shift;
 
      }
      /*}}}*/
      void Copy(IssmAbsVec<doubletype>* toin){/*{{{*/
 
         int i;
 
         /*Assume toin is of the correct type, and downcast: */
         IssmSeqVec* to=NULL;
 
         to=(IssmSeqVec<doubletype>*)toin;
 
         to->M=this->M;
         for(i=0;i<this->M;i++)to->vector[i]=this->vector[i];
 
      }
      /*}}}*/
      doubletype Norm(NormMode mode){/*{{{*/
 
         doubletype norm;
         int i;
 
         switch(mode){
            case NORM_INF:
               norm=0.; for(i=0;i<this->M;i++)norm=max(norm,fabs(this->vector[i]));
               //norm=0.; for(i=0;i<this->M;i++)norm=max(norm,this->vector[i]);
               return norm;
               break;
            case NORM_TWO:
               norm=0.; 
               for(i=0;i<this->M;i++)norm+=this->vector[i]*this->vector[i];
               return sqrt(norm);
               break;
            default:
               _error_("unknown norm !");
               break;
         }
      }
      /*}}}*/
      void Scale(doubletype scale_factor){/*{{{*/
 
         int i;
         for(i=0;i<this->M;i++)this->vector[i]=scale_factor*this->vector[i];
 
      }
      /*}}}*/
      doubletype Dot(IssmAbsVec<doubletype>* inputin){/*{{{*/
 
         int i;
 
         /*Assume inputin is of the correct type, and downcast: */
         IssmSeqVec* input=NULL;
 
         input=(IssmSeqVec<doubletype>*)inputin;
 
         doubletype dot=0;
         for(i=0;i<this->M;i++)dot+=this->vector[i]*input->vector[i];
         return dot;
 
      }
      /*}}}*/
      void PointwiseDivide(IssmAbsVec<doubletype>* xin,IssmAbsVec<doubletype>* yin){/*{{{*/
 
         int i;
 
         /*Assume xin and yin are of the correct type, and downcast: */
         IssmSeqVec* x=NULL;
         IssmSeqVec* y=NULL;
 
         x=(IssmSeqVec<doubletype>*)xin;
         y=(IssmSeqVec<doubletype>*)yin;
 
         /*pointwise w=x/y where this->vector is w: */
         for(i=0;i<this->M;i++)this->vector[i]=x->vector[i]/y->vector[i];
      }
      /*}}}*/
};
#endif //#ifndef _ISSM_SEQ_VEC_H_