SPIVec.cpp

Go to the documentation of this file.

#include "SPIVec.hpp"
// use 1 to use the GPU and 0 or anything else to only use CPU and MPI.  Be sure this matches what is in SPIMat.cpp
#define USE_GPU 0
 
namespace SPI{
 
    // constructors
    SPIVec::SPIVec(
            std::string _name 
            ){name=_name; }
    SPIVec::SPIVec(
            const SPIVec &A, 
            std::string _name 
            ){
        name=_name; 
        (*this) = A;
    }
    SPIVec::SPIVec(
            const PetscInt _rows,  
            const std::string _name
            ){
        //this->Init(_rows,_name);
        this->name=_name;
        this->rows=_rows;
        //PetscInt _rows2 = this->rows;
        //SPI::printf("_rows at initialization %D",_rows);
        //SPI::printf(" rows at initialization %D",_rows2);
        //std::cout<<" rows at initialization "<<this->rows<<std::endl;
        ierr = VecCreate(PETSC_COMM_WORLD,&vec);CHKERRXX(ierr);
        ierr = VecSetSizes(vec,PETSC_DECIDE,_rows);CHKERRXX(ierr);
#if USE_GPU == 1
        ierr = VecSetType(vec,VECMPICUDA);CHKERRXX(ierr);
#else
        ierr = VecSetType(vec,VECMPI);CHKERRXX(ierr);
#endif
        flag_init=PETSC_TRUE;
    }
 
    // Initialize vector
    PetscInt SPIVec::Init(
            PetscInt _rows,     
            const std::string _name   
            ){
        this->name=_name;
        rows=_rows;
        SPI::printf("_rows at initialization %D",_rows);
        SPI::printf("rows at initialization %D",this->rows);
        ierr = VecCreate(PETSC_COMM_WORLD,&vec);CHKERRQ(ierr);
        ierr = VecSetSizes(vec,PETSC_DECIDE,_rows);CHKERRQ(ierr);
#if USE_GPU == 1
        ierr = VecSetType(vec,VECMPICUDA);CHKERRQ(ierr);
#else
        ierr = VecSetType(vec,VECMPI);CHKERRQ(ierr);
#endif
        flag_init=PETSC_TRUE;
        return 0;
    }
 
    PetscInt SPIVec::set(
            const PetscInt _row,  
            const PetscScalar v   
            ){
        PetscInt low,high;
        VecGetOwnershipRange(vec,&low,&high);
        if ((low <= _row) && (_row < high)){
            ierr = VecSetValue(vec,_row,v,INSERT_VALUES);CHKERRQ(ierr);
        }
        return 0;
    }
    PetscInt SPIVec::set(
            const PetscScalar v   
            ){
        ierr = VecSet(vec,v); CHKERRQ(ierr);
        return 0;
    }
    PetscInt SPIVec::add(
            PetscInt _row,      
            const PetscScalar v 
            ){
        PetscInt low,high;
        VecGetOwnershipRange(vec,&low,&high);
        if ((low <= _row) && (_row < high)){
            ierr = VecSetValue(vec,_row,v,ADD_VALUES);CHKERRQ(ierr);
        }
        return 0;
    }
    // get size of vector
    PetscInt SPIVec::size(){
        PetscInt n;
        ierr = VecGetSize(vec,&n);CHKERRXX(ierr);
        rows=n; // update rows
        return n;
    }
 
    // overloaded operators, get
    PetscScalar SPIVec::operator()(
            PetscInt _row, 
            PetscBool global 
            )const{
        PetscScalar v,v_global=0.;
        PetscInt low,high;
        PetscErrorCode ierr2;
        ierr2 = VecGetOwnershipRange(vec,&low, &high);CHKERRXX(ierr2);
        if ((low<=_row) && (_row<high)){
            ierr2 = VecGetValues(vec,1,&_row,&v);CHKERRXX(ierr2);
        }
        if (global){
            MPIU_Allreduce(&v,&v_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
        }
        else{
            v_global=v; // return local value
        }
        return v_global;
    }
    PetscScalar SPIVec::operator()(
            PetscInt _row, 
            PetscBool global 
            ){
        PetscScalar v,v_global=0.;
        PetscInt low,high;
        ierr = VecGetOwnershipRange(vec,&low, &high);CHKERRXX(ierr);
        if ((low<=_row) && (_row<high)){
            ierr = VecGetValues(vec,1,&_row,&v);CHKERRXX(ierr);
        }
        if (global){
            MPIU_Allreduce(&v,&v_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
        }
        else{
            v_global=v; // return local value
        }
        return v_global;
    }
    // overloaded operator, set
    PetscInt SPIVec::operator()(
            PetscInt _row,      
            const PetscScalar v 
            ){
        PetscInt low,high;
        VecGetOwnershipRange(vec,&low,&high);
        if ((low <= _row) && (_row < high)){
            ierr = VecSetValue(vec,_row,v,INSERT_VALUES);CHKERRQ(ierr);
        }
        //(*this)(); // assemble after every insertion
        return 0;
    }
    PetscInt SPIVec::operator()(PetscInt _row, const double v){
        ierr = (*this)(_row,(PetscScalar)(v+0.0*PETSC_i));CHKERRQ(ierr);
        return 0;
    }
    PetscInt SPIVec::operator()(PetscInt _row, const int v){
        ierr = (*this)(_row,(PetscScalar)((double)v+0.0*PETSC_i));CHKERRQ(ierr);
        return 0;
    }
 
    // overloaded operator, assemble
    SPIVec& SPIVec::operator()(){
        ierr = VecAssemblyBegin(vec);CHKERRXX(ierr);
        ierr = VecAssemblyEnd(vec);CHKERRXX(ierr);
        return (*this);
    }
    // overloaded operator, VecAXPY
    SPIVec& SPIVec::operator+=(
            const SPIVec &X 
            ){
        ierr = VecAXPY(this->vec,1.,X.vec);CHKERRXX(ierr);
        return *this;
    }
    SPIVec& SPIVec::axpy(
            const PetscScalar a,    
            const SPIVec &X         
            ){
        ierr = VecAXPY(this->vec,a,X.vec);CHKERRXX(ierr);
        return (*this);
    }
    // overloaded operator, VecAXPY
    SPIVec SPIVec::operator+(
            const SPIVec &X 
            ){
        SPIVec A;
        A=*this;
        ierr = VecAXPY(A.vec,1.,X.vec);CHKERRXX(ierr);
#if USE_GPU == 1
        ierr = VecSetType(A.vec,VECMPICUDA);CHKERRXX(ierr);
#else
        ierr = VecSetType(A.vec,VECMPI);CHKERRXX(ierr);
#endif
        return A;
    }
    SPIVec SPIVec::operator+(
            const PetscScalar a 
            ){ // Y + a operation
        SPIVec A;
        A=(*this);
        A += a*ones(rows);
        return A;
    }
    SPIVec SPIVec::operator+(
            const double a 
            ){ // Y + a operation
        SPIVec A;
        A=(*this);
        A += a*ones(rows);
        return A;
    }
    SPIVec SPIVec::operator-(
            const PetscScalar a 
            ){ // Y - a operation
        SPIVec A;
        A=(*this);
        A -= a*ones(rows);
        return A;
    }
    SPIVec SPIVec::operator-(
            const PetscInt a 
            ){ // Y - a operation
        SPIVec A;
        A=(*this);
        A -= a*ones(rows);
        return A;
    }
    // overloaded operator, VecAXPY
    SPIVec& SPIVec::operator-=(
            const SPIVec &X 
            ){
        ierr = VecAXPY(this->vec,-1.,X.vec);CHKERRXX(ierr);
        return *this;
    }
    // overloaded operator, VecAXPY
    SPIVec SPIVec::operator-(
            const SPIVec &X 
            ){
        SPIVec A;
        A=*this;
        ierr = VecAXPY(A.vec,-1.,X.vec);CHKERRXX(ierr);
#if USE_GPU == 1
        ierr = VecSetType(A.vec,VECMPICUDA);CHKERRXX(ierr);
#else
        ierr = VecSetType(A.vec,VECMPI);CHKERRXX(ierr);
#endif
        return A;
    }
    SPIVec SPIVec::operator-() const{
        return -1.*(*this);
    }
    // overload operator, scale with scalar
    SPIVec SPIVec::operator*(
            const PetscScalar a 
            ){
        SPIVec A;
        A=(*this);
        ierr = VecScale(A.vec,a);CHKERRXX(ierr);
        return A;
    }
    SPIVec SPIVec::operator*(const double a){
        PetscScalar as=a;
        SPIVec A;
        A=*this;
        ierr = VecScale(A.vec,as);CHKERRXX(ierr);
        return A;
    }
    // overload operator, scale with scalar
    SPIVec& SPIVec::operator*=(
            const PetscScalar a 
            ){
        ierr = VecScale(this->vec,a);CHKERRXX(ierr);
        return *this;
    }
    SPIVec& SPIVec::operator*=(
            const double a 
            ){
        ierr = VecScale(this->vec,(PetscScalar)(a+0.*PETSC_i));CHKERRXX(ierr);
        return *this;
    }
    SPIVec& SPIVec::operator*=(
            const SPIVec &a 
            ){
        ierr = VecPointwiseMult(vec,a.vec,(*this).vec);CHKERRXX(ierr);
        return *this;
    }
    // overload operator, pointwise multiply
    SPIVec SPIVec::operator*(
            const SPIVec& X 
            ){
        SPIVec C;
        C.rows=rows;
        C=(*this);
        ierr = VecPointwiseMult(C.vec,X.vec,(*this).vec);CHKERRXX(ierr);
        // ierr = VecSetType(C.vec,VECMPI);CHKERRXX(ierr);
        return C;
    }
    // overload operator, pointwise divide
    SPIVec SPIVec::operator/(
            const PetscScalar a 
            ){
        SPIVec A;
        A=(*this);
        ierr = VecScale(A.vec,1./a);CHKERRXX(ierr);
        return A;
    }
    SPIVec SPIVec::operator/(const double a){
        PetscScalar as=a;
        SPIVec A;
        A=*this;
        ierr = VecScale(A.vec,1./as);CHKERRXX(ierr);
        return A;
    }
    SPIVec SPIVec::operator/(
            const SPIVec X 
            ){
        SPIVec Z(X.rows);
        ierr = VecPointwiseDivide(Z.vec,this->vec,X.vec);CHKERRXX(ierr);
        return Z;
    }
    // overload operator, scale with scalar
    SPIVec& SPIVec::operator/=(
            const PetscScalar a 
            ){
        ierr = VecScale(this->vec,1./a);CHKERRXX(ierr);
        return *this;
    }
    // ^ operator
    SPIVec SPIVec::operator^(
            const PetscScalar p 
            ){
        return pow(*this,p);
    }
    SPIVec SPIVec::operator^(
            const double p 
            ){
        return pow(*this,(PetscScalar)p);
    }
    SPIVec SPIVec::operator^(
            const int p 
            ){
        return pow(*this,(PetscScalar)p);
    }
    SPIVec SPIVec::operator^(
            SPIVec p 
            ){
        return pow(*this,p);
    }
    // overload operator, copy and initialize
    SPIVec& SPIVec::operator=(const SPIVec &X){
        if(flag_init){
            if(X.rows==this->rows){
                ierr = VecCopy(X.vec,vec);CHKERRXX(ierr); // use copy if size matches
            }
            else{
                this->~SPIVec(); // destroy and recreate from scratch
                this->rows=X.rows;
                ierr = VecDuplicate(X.vec,&vec); CHKERRXX(ierr);
                ierr = VecCopy(X.vec,vec); CHKERRXX(ierr);
            }
        }
        else{
            this->rows=X.rows;
            ierr = VecDuplicate(X.vec,&vec); CHKERRXX(ierr);
            ierr = VecCopy(X.vec,vec); CHKERRXX(ierr);
            flag_init=PETSC_TRUE;
        }
        //ierr = VecSetType(vec,VECMPI);CHKERRXX(ierr);
        return (*this);
    }
    PetscBool SPIVec::operator==(
            const SPIVec &x2    
            ){
        PetscBool iftrue;
        ierr = VecEqual(vec,x2.vec,&iftrue); CHKERRXX(ierr);
        return iftrue;
 
    }
 
 
    // overload % for inner product
    //SPIVec operator%(SPIVec A){
    //return *this;
    //}     
    SPIVec& SPIVec::conj(){
        ierr = VecConjugate(vec);CHKERRXX(ierr);
        return (*this);
    }
 
    PetscScalar SPIVec::max(){
        PetscInt argmax;
        PetscReal max=0.;
        PetscScalar maxscalar;
        ierr = VecMax(this->vec,&argmax,&max);CHKERRXX(ierr);
        maxscalar = (*this)(argmax,PETSC_TRUE);
        //PetscScalar maxscalar_global = 0.;
        //MPIU_Allreduce(&maxscalar,&maxscalar_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
 
        return maxscalar;
    }
    PetscInt SPIVec::argmax(){
        PetscInt argmax;
        PetscReal max=0.;
        //PetscScalar maxscalar;
        ierr = VecMax(this->vec,&argmax,&max);CHKERRXX(ierr);
        //maxscalar = (*this)(argmax,PETSC_TRUE);
        //PetscScalar maxscalar_global = 0.;
        //MPIU_Allreduce(&maxscalar,&maxscalar_global,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
 
        return argmax;
    }
    SPIVec& SPIVec::real(){
        ierr = VecRealPart(vec); CHKERRXX(ierr);
        return (*this);
    }
    SPIVec& SPIVec::imag(){
        ierr = VecImaginaryPart(vec); CHKERRXX(ierr);
        return (*this);
    }
    PetscScalar SPIVec::dot(
            SPIVec y    
            ){
        PetscScalar val;
        ierr = VecDot(vec,y.vec,&val); CHKERRXX(ierr);
        return val;
    }
 
 
    // print vector to screen
    PetscInt SPIVec::print(){
        (*this)();// assemble
        printf("\n---------------- "+this->name+"---start------");
        //PetscPrintf(PETSC_COMM_WORLD,("\n---------------- "+name+"---start------\n").c_str());
        //SPI::printf("shape = %d x 1",this->rows);
        SPI::printf("shape = "+std::to_string(this->rows)+" x 1");
        ierr = VecView(vec,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
        //PetscPrintf(PETSC_COMM_WORLD,("---------------- "+name+"---done-------\n\n").c_str());
        printf("---------------- "+this->name+"---done-------\n");
        return 0;
    }
 
    SPIVec::~SPIVec(){
        if(flag_init){
            flag_init=PETSC_FALSE;
            ierr = VecDestroy(&vec);CHKERRXX(ierr);
        }
        //else{
            //(*this)=zeros(1);
            //ierr = VecDestroy(&vec);CHKERRXX(ierr);
 
        //}
    }
 
    // overload operator, scale with scalar
    SPIVec operator/(
            const PetscScalar a, 
            const SPIVec &Y     
            ){
        SPIVec B(Y.rows);
        SPIVec A(ones(Y.rows)*a);
        //B=Y;
        //B.ierr = VecScale(B.vec,a);CHKERRXX(B.ierr);
        B.ierr = VecPointwiseDivide(B.vec,A.vec,Y.vec);CHKERRXX(B.ierr);
        return B;
    }
 
    SPIVec operator*(
            const PetscScalar a, 
            const SPIVec &Y     
            ){
        SPIVec B;
        B=Y;
        B.ierr = VecScale(B.vec,a);CHKERRXX(B.ierr);
        return B;
    }
 
    SPIVec operator+(
            const PetscScalar a,  
            const SPIVec &Y         
            ){
        SPIVec B;
        B = Y;
        B += a*ones(B.rows);
        return B;
    }
 
    SPIVec operator-(
            const PetscScalar a, 
            const SPIVec &Y     
            ){
        SPIVec B;
        B = -1.*Y;
        B += a*ones(B.rows);
        return B;
    }
 
    PetscInt save(
            const SPIVec &A,            
            const std::string filename  
            ){ // save A to hdf5 to filename as variable A.name
        PetscErrorCode ierr;
        ierr = PetscObjectSetName((PetscObject)A.vec, A.name.c_str());CHKERRQ(ierr);
        PetscViewer viewer;
        std::ifstream f(filename.c_str());
        if(f.good()){// if filename exists, then append
            ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_APPEND, &viewer); CHKERRQ(ierr);
        }
        else{ // if not, then write a new file
            ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_WRITE, &viewer); CHKERRQ(ierr);
        }
        ierr = VecView(A.vec,viewer); CHKERRQ(ierr);
        ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
        return 0;
    }
    PetscInt save(
            std::vector<PetscScalar> A, 
            std::string variablename,   
            std::string filename        
            ){
        PetscInt n=A.size();
        SPIVec Avec(n,variablename);
        for(PetscInt i=0; i<n; i++){
            Avec(i,A[i]);
        }
        save(Avec,filename);
        return 0;
    }
    PetscInt save(
            std::vector<SPIVec> A,      
            std::string variablename,   
            std::string filename        
            ){
        PetscInt n=A.size();
        std::string sep="_";
        for(PetscInt i=0; i<n; i++){
            A[i].name = variablename+sep+std::to_string(i);
            save(A[i],filename);
        }
        return 0;
    }
 
    PetscInt load( 
            SPIVec &A,                  
            const std::string filename  
            ){
        A.ierr = PetscObjectSetName((PetscObject)A.vec, A.name.c_str());CHKERRQ(A.ierr);
        PetscViewer viewer;
        //std::ifstream f(filename.c_str());
        A.ierr = PetscViewerHDF5Open(PETSC_COMM_WORLD, filename.c_str(), FILE_MODE_READ, &viewer); CHKERRQ(A.ierr);
        A.ierr = VecLoad(A.vec,viewer); CHKERRQ(A.ierr);
        A.ierr = PetscViewerDestroy(&viewer); CHKERRQ(A.ierr);
        return 0;
    }
 
    SPIVec ones(
            const PetscInt rows 
            ){
        SPIVec A(rows);
        A.ierr = VecSet(A.vec,1.);CHKERRXX(A.ierr);
        return A;
    }
 
    SPIVec zeros(
            const PetscInt rows 
            ){
        SPIVec A(rows);
        A.ierr = VecSet(A.vec,0.);CHKERRXX(A.ierr);
        return A;
    }
    SPIVec conj(
            const SPIVec &A 
            ){
        SPIVec B;
        B=A;
        B.ierr = VecConjugate(B.vec);CHKERRXX(B.ierr);
        return B;
    }
    SPIVec real(
            const SPIVec &A     
            ){
        SPIVec B(A);
        return B.real();
    }
    SPIVec imag(
            const SPIVec &A     
            ){
        SPIVec B(A);
        return B.imag();
    }
    SPIVec linspace(
            const PetscScalar begin,    
            const PetscScalar end,      
            const PetscInt _rows         
            ){ // return linspace of number of rows equally spaced points between begin and end
        SPIVec y(_rows);
        //PetscInt _rows2 = y.rows;
        //SPI::printf("y.rows in linspace = %D",_rows2);
        PetscScalar step = (end-begin)/((PetscScalar)(_rows-1));
        PetscScalar value=begin;
        //PetscInt i=0;
        for (PetscInt i=0; i<_rows; ++i){
            y(i,value);
            value += step;
        }
        y();
        //_rows2 = y.rows;
        //SPI::printf("y.rows in linspace = %D",_rows2);
        //SPI::printf("_rows in linspace = %D",_rows);
        return y;
    }
    SPIVec arange(
            const PetscScalar begin,    
            const PetscScalar end,      
            const PetscScalar stepsize  
            ){ // return linspace of number of rows equally spaced points between begin and end
        PetscInt _rows=ceil(PetscRealPart((end-begin)/stepsize));
        SPIVec y(_rows);
        //PetscScalar step = (end-begin)/((PetscScalar)(_rows-1));
        PetscScalar value=begin;
        //PetscInt i=0;
        for (PetscInt i=0; i<_rows; ++i){
            y(i,value);
            value += stepsize;
        }
        y();
        return y;
    }
    SPIVec arange(
            const PetscScalar end      
            ){ // return linspace of number of rows equally spaced points between begin and end
        PetscScalar begin=0.;
        return arange(begin,end);
    }
 
    template <class T>
        SPIVec _Function_on_each_element(
                T (*f)(T const&),   
                const SPIVec &A     
                ){
            SPIVec out(A);
            for (PetscInt i=0; i<out.rows; ++i){
                out(i,(*f)(out(i)));                // TODO speed up by getting all values at once on local processor and looping through those
            }
            out();
            return out;
        }
 
    SPIVec sin(const SPIVec &A){ return _Function_on_each_element(&std::sin<PetscReal>, A); }
    SPIVec cos(const SPIVec &A){ return _Function_on_each_element(&std::cos<PetscReal>, A); }
    SPIVec tan(const SPIVec &A){ return _Function_on_each_element(&std::tan<PetscReal>, A); }
    SPIVec exp(const SPIVec &A){ 
        //return _Function_on_each_element(&std::exp<PetscReal>, A); 
        SPIVec B(A);
        B.ierr = VecExp(B.vec); CHKERRXX(B.ierr);
        return B;
    }
    SPIVec log(const SPIVec &A){ return _Function_on_each_element(&std::log<PetscReal>, A); }
    SPIVec log10(const SPIVec &A){ return _Function_on_each_element(&std::log10<PetscReal>, A); }
    SPIVec sinh(const SPIVec &A){ return _Function_on_each_element(&std::sinh<PetscReal>, A); }
    SPIVec cosh(const SPIVec &A){ return _Function_on_each_element(&std::cosh<PetscReal>, A); }
    SPIVec tanh(const SPIVec &A){ return _Function_on_each_element(&std::tanh<PetscReal>, A); }
    SPIVec asin(const SPIVec &A){ return _Function_on_each_element(&std::asin<PetscReal>, A); }
    SPIVec acos(const SPIVec &A){ return _Function_on_each_element(&std::acos<PetscReal>, A); }
    SPIVec atan(const SPIVec &A){ return _Function_on_each_element(&std::atan<PetscReal>, A); }
    SPIVec asinh(const SPIVec &A){ return _Function_on_each_element(&std::asinh<PetscReal>, A); }
    SPIVec acosh(const SPIVec &A){ return _Function_on_each_element(&std::acosh<PetscReal>, A); }
    SPIVec atanh(const SPIVec &A){ return _Function_on_each_element(&std::atanh<PetscReal>, A); }
    SPIVec sqrt(const SPIVec &A){ return _Function_on_each_element(&std::sqrt<PetscReal>, A); }
    SPIVec erf(const SPIVec &A){
            SPIVec out(A);
            for (PetscInt i=0; i<out.rows; ++i){
                out(i,std::erf((double)(PetscRealPart(out(i)))));                // TODO speed up by getting all values at once on local processor and looping through those
            }
            out();
            return out;
    }
    SPIVec erfc(const SPIVec &A){
            SPIVec out(A);
            for (PetscInt i=0; i<out.rows; ++i){
                out(i,std::erfc((double)(PetscRealPart(out(i)))));                // TODO speed up by getting all values at once on local processor and looping through those
            }
            out();
            return out;
    }
    template <class T>
        SPIVec _Function_on_each_element(
                T (*f)(T const&, T const&),     
                const SPIVec &A,                
                SPIVec &B                 
                ){
            SPIVec out(A);
            for (PetscInt i=0; i<out.rows; ++i){
                out(i,(*f)(out(i),B(i)));                // TODO speed up by getting all values at once on local processor and looping through those
            }
            out();
            return out;
        }
    SPIVec pow(const SPIVec &A,SPIVec &B){ return _Function_on_each_element(&std::pow<PetscReal>, A,B); }
    SPIVec pow(
            const SPIVec &A, 
            PetscScalar b       
            ){
        SPIVec B(A);
        B.ierr = VecPow(B.vec,b);CHKERRXX(B.ierr);
        return B;
 
    }
 
    PetscScalar dot(
            SPIVec x,   
            SPIVec y    
            ){
        PetscScalar innerproduct;
        x.ierr = VecDot(x.vec,y.vec,&innerproduct); CHKERRXX(x.ierr);
        return innerproduct;
    }
 
    SPIVec abs(const SPIVec &A){ 
        SPIVec B(A);
        VecAbs(B.vec);
        return B;
    }
 
    PetscScalar sum(
            SPIVec x1   
            ){
        PetscScalar sum;
        x1.ierr = VecSum(x1.vec,&sum); CHKERRXX(x1.ierr);
        return sum;
    }
    PetscScalar integrate_coeffs(
            const SPIVec &a     
            ){
        PetscInt n=a.rows;
        PetscInt N=n-1;
        SPIVec d(n+1,"d");
        PetscScalar value=0.0,d0=0.0;
        for(PetscInt i=1; i<=N+1; ++i){
            if(i==1){
                value=0.5*(2.0*a(i-1,PETSC_TRUE)-a(i+1,PETSC_TRUE));
            }
            else if(i==N+1){
                value = 0.5/((PetscScalar)i)*(2.0*a(i-1,PETSC_TRUE));
            }
            else if(i==N){
                value = 0.5/((PetscScalar)i)*(1.0*a(i-1,PETSC_TRUE));
            }
            else{
                value = 0.5/((PetscScalar)i)*(1.0*a(i-1,PETSC_TRUE) - a(i+1,PETSC_TRUE));
            }
            d(i,value);
            if((i%2)==0){ // even
                d0 -= value;
            }
            else{ // odd
                d0 += value;
            }
        }
        d(0,d0);
        d(); // assemble
        return sum(d);
    }
    PetscScalar integrate_coeffs(
            const SPIVec &a,      
            const SPIVec &y       
            ){
        PetscScalar ai=y(0,PETSC_TRUE);
        PetscScalar bi=y(y.rows-1,PETSC_TRUE);
        PetscScalar mul = (bi-ai)/2.0;
        return mul*integrate_coeffs(a);
    }
 
 
 
    PetscReal L2(
            SPIVec x1,      
            const SPIVec x2,      
            NormType type   
            ){
        PetscReal error;
        VecNorm((x1-x2).vec,type,&error);
        return error;
    }
 
 
    PetscReal L2(
            const SPIVec x1,        
            NormType type           
            ){
        PetscReal error;
        VecNorm(x1.vec,type,&error);
        return error;
    }
 
    SPIVec diff( 
            SPIVec x       
            ){ 
        SPIVec x0(x.rows-1), x1(x.rows-1);
        // set x0=x[i] and x1=x[i+1]
        for (PetscInt i=0; i<x.rows-1; ++i){
            x0(i,x(i,PETSC_TRUE));
            x1(i,x(i+1,PETSC_TRUE));
        }
        // assemble
        x0();
        x1();
        // return difference x[i+1]-x[i]
        return x1-x0;
    }
 
    PetscScalar trapz(
            SPIVec y      
            ){
        SPIVec y0(y.rows-1), y1(y.rows-1);
        // set y0=y[i] and y1=y[i+1]
        for (PetscInt i=0; i<y.rows-1; ++i){
            y0(i,y(i,PETSC_TRUE));
            y1(i,y(i+1,PETSC_TRUE));
        }
        // assemble
        y0();
        y1();
        // return trapezoidal rule sum((y[i+1]+y[i])/2.
        return sum((y1+y0)/2.);
    }
 
    PetscScalar trapz(
            SPIVec y,     
            SPIVec x      
            ){
        SPIVec y0(y.rows-1), y1(y.rows-1);
        // set y0=y[i] and y1=y[i+1]
        for (PetscInt i=0; i<y.rows-1; ++i){
            y0(i,y(i,PETSC_TRUE));
            y1(i,y(i+1,PETSC_TRUE));
        }
        // assemble
        y0();
        y1();
        // return trapezoidal rule sum((y[i+1]+y[i])/2. * diff(x))
        return sum((y1+y0)/2. * diff(x));
    }
    PetscInt draw(
            const SPIVec &x     
            ){
        PetscViewer viewer;
        PetscErrorCode ierr;
        ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,NULL,x.name.c_str(),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,&viewer);CHKERRQ(ierr);
        ierr = VecView(x.vec,viewer);CHKERRQ(ierr);
        SPI::printf("  draw(SPIVec) with title=%s, hit Enter to continue",x.name.c_str());
        std::cin.ignore();
 
        ierr = PetscViewerDestroy(&viewer); CHKERRQ(ierr);
        return 0;
    }
    
}