MatrixWrapper.hpp

Go to the documentation of this file.

#ifndef OPM_MIXED_MATRIX_HEADER_INCLUDED
#define OPM_MIXED_MATRIX_HEADER_INCLUDED
 
 
#include <opm/simulators/linalg/mixed/bsr.h>
 
#include <opm/common/ErrorMacros.hpp>
 
#include <stdexcept>
 
namespace Opm
{
 
template <class Vector>
class MixedMatrixWrapper
{
    public:
 
    // extract block size
    static constexpr auto block_size = Vector::block_type::dimension;
 
    MixedMatrixWrapper(int nrows, int nnz)
    {
        nnz_=nnz;
        M_ = bsr_alloc();
        bsr_init(M_, nrows, nnz, block_size);
    }
 
    ~MixedMatrixWrapper() {bsr_free(M_);}
 
    void update(double const *data);
 
    void mv(const Vector& x, Vector& y) const;
 
    void umv(const Vector& x, Vector& y) const;
 
    void usmv(double alpha, const Vector& x, Vector& y) const;
 
    int *rowptr(){return M_->rowptr;}
 
    int *colidx(){return M_->colidx;}
 
    int nrows() const {return M_->nrows;}
 
    private:
    int nnz_;
    bsr_matrix  *M_;
};
 
template <class Vector>
void MixedMatrixWrapper<Vector>::
mv(const Vector& x, Vector& y) const
{
    int constexpr b = block_size;
    if      constexpr(b==1){OPM_THROW(std::invalid_argument, "MixedMatrixWrapper::mv does not support block size == 1!\n");}
    else if constexpr(b==2) bsr_vmspmv2(M_, &x[0][0], &y[0][0]);
    else if constexpr(b==3) bsr_vmspmv3(M_, &x[0][0], &y[0][0]);
    else if constexpr(b==4) bsr_vmspmv4(M_, &x[0][0], &y[0][0]);
    else
    {
        int nrows = M_->nrows;
        int *rowptr=M_->rowptr;
        int *colidx=M_->colidx;
        const float *data=M_->flt;
 
        int constexpr bb = b*b;
        double yy[bb];
        for(int i=0;i<nrows;i++)
        {
            for(int k=0;k<bb;k++) yy[k] = 0.0;
            for(int k=rowptr[i];k<rowptr[i+1];k++)
            {
                const float *AA=data+bb*k;
 
                int j = colidx[k];
                double const *xx = &x[0][0]+b*j;
                for(int n=0;n<b;n++) for(int m=0;m<b;m++) yy[m+b*n] += AA[m+b*n]*xx[n];
            }
 
            // sum over columns
            double z[b];
            for(int m=0;m<b;m++) z[m] = yy[m];
            for(int n=1;n<b;n++) for(int m=0;m<b;m++) z[m] += yy[m+b*n];
 
            double *y_i = &y[0][0] + b*i;
            for(int m=0;m<b;m++) y_i[m] = z[m];
        }
    }
}
 
template <class Vector>
void MixedMatrixWrapper<Vector>::
umv(const Vector& x, Vector& y) const
{
    int constexpr b = block_size;
    if      constexpr(b==1){OPM_THROW(std::invalid_argument, "MixedMatrixWrapper::umv does not support block size == 1!\n");}
    else if constexpr(b==2) bsr_vmspumv2(M_, &x[0][0], &y[0][0], 1.0);
    else if constexpr(b==3) bsr_vmspumv3(M_, &x[0][0], &y[0][0], 1.0);
    else if constexpr(b==4) bsr_vmspumv4(M_, &x[0][0], &y[0][0], 1.0);
    else {OPM_THROW(std::invalid_argument, "MixedMatrixWrapper::umv does not support block size == 1!\n");}
}
 
template <class Vector>
void MixedMatrixWrapper<Vector>::
usmv(double alpha, const Vector& x, Vector& y) const
{
    int constexpr b = block_size;
    if      constexpr(b==1){OPM_THROW(std::invalid_argument, "MixedMatrixWrapper::usmv does not support block size == 1!\n");}
    else if constexpr(b==2) bsr_vmspumv2(M_, &x[0][0], &y[0][0], alpha);
    else if constexpr(b==3) bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha);
    else if constexpr(b==4) bsr_vmspumv4(M_, &x[0][0], &y[0][0], alpha);
    else
    {
        int nrows = M_->nrows;
        int *rowptr=M_->rowptr;
        int *colidx=M_->colidx;
        const float *data=M_->flt;
 
        int constexpr bb = b*b;
        double yy[bb];
        for(int i=0;i<nrows;i++)
        {
            for(int k=0;k<bb;k++) yy[k] = 0.0;
            for(int k=rowptr[i];k<rowptr[i+1];k++)
            {
                const float *AA=data+bb*k;
 
                int j = colidx[k];
                double const *xx = &x[0][0]+b*j;
                for(int n=0;n<b;n++) for(int m=0;m<b;m++) yy[m+b*n] += AA[m+b*n]*xx[n];
            }
 
            // sum over columns
            double z[b];
            for(int m=0;m<b;m++) z[m] = yy[m];
            for(int n=1;n<b;n++) for(int m=0;m<b;m++) z[m] += yy[m+b*n];
 
            double *y_i = &y[0][0] + b*i;
            for(int m=0;m<b;m++) y_i[m] += alpha*z[m];
        }
 
    }
}
 
template <class Vector>
void MixedMatrixWrapper<Vector>::
update(double const *data)
{
    // transpose each dense block to make them column-major
    int constexpr b = block_size;
    int constexpr bb=b*b;
    double B[bb];
    for(int k=0;k<nnz_;k++)
    {
        for(int i=0;i<b;i++) for(int j=0;j<b;j++) B[b*j+i] = data[bb*k + b*i + j];
        for(int i=0;i<bb;i++) M_->flt[bb*k + i] = B[i];
    }
}
 
 
} // namespace Opm
#endif // OPM_MIXED_MATRIX_HEADER_INCLUDED