diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f56331543ca..2d14e3593b3 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -281,6 +281,7 @@ list(APPEND PRIVATE_HEADER_FILES if (HAVE_AVX2_EXTENSION) set (AVX2_SOURCE_FILES + opm/simulators/linalg/mixed/dot.c opm/simulators/linalg/mixed/bsr.c opm/simulators/linalg/mixed/prec.c opm/simulators/linalg/mixed/bslv.c) diff --git a/opm/simulators/linalg/FlexibleSolver_impl.hpp b/opm/simulators/linalg/FlexibleSolver_impl.hpp index 4b0f58f8957..8209df90d46 100644 --- a/opm/simulators/linalg/FlexibleSolver_impl.hpp +++ b/opm/simulators/linalg/FlexibleSolver_impl.hpp @@ -219,13 +219,13 @@ namespace Dune maxiter, // maximum number of iterations verbosity); #if HAVE_AVX2_EXTENSION - } else if (solver_type == "mixed-bicgstab") { + } else if (solver_type == "mixed-legacy") { if constexpr (Opm::is_gpu_operator_v) { - OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operators"); + OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for GPU operators"); } else if constexpr (Opm::detail::is_multi_type_block_vector_v) { - OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for multi-type block vectors."); + OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for multi-type block vectors."); } else if constexpr (std::is_same_v){ - OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for single precision."); + OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for single precision."); } else { const std::string prec_type = prm.get("preconditioner.type", "error"); bool use_mixed_dilu= (prec_type=="legacy-mixed-dilu"); @@ -238,13 +238,13 @@ namespace Dune ); } // MixedBiCGSTABSolver starts here - } else if (solver_type == "mixed-precision") { + } else if (solver_type == "mixed-bicgstab") { if constexpr (Opm::is_gpu_operator_v) { - OPM_THROW(std::invalid_argument, "mixed-precision solver not supported for GPU operators"); + OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operators"); } else if constexpr (Opm::detail::is_multi_type_block_vector_v) { OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for multi-type block vectors."); } else if constexpr (std::is_same_v){ - OPM_THROW(std::invalid_argument, "mixed-precision solver not supported for single precision."); + OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for single precision."); } else { linsolver_ = std::make_shared>(linearoperator_for_solver_, scalarproduct_, diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp new file mode 100644 index 00000000000..4d26a866717 --- /dev/null +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -0,0 +1,163 @@ +#ifndef OPM_SCALAR_PRODUCTS_HEADER_INCLUDED +#define OPM_SCALAR_PRODUCTS_HEADER_INCLUDED + +#include + +namespace Dune +{ + +/// A parallel scalar product that takes advantage of the fact that all +/// elements associated with ghost cells are located at the end of the +/// vector. This allows us to ignore the block structure of the vector +/// and eliminate the use of a mask to exclude ghost entries from being +/// included in the scalar product +template +class GhostLastScalarProduct : public ScalarProduct +{ + public: + + /*! \brief constructor + * \param com The communication object for syncing overlap and copy + * data points. + * \param cat parallel solver category (nonoverlapping or overlapping) + */ + GhostLastScalarProduct (std::shared_ptr com, SolverCategory::Category cat) + : _communication(com), _category(cat) + { + count_ = getLocalCount(); // number or local cells + int verify = verifyLocalCount(); // redundant check on numbef of local cells + if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n"); + } + + /*! \brief constructor + * \param com The communication object for syncing overlap and copy + * data points. + * \param cat parallel solver category (nonoverlapping or overlapping) + * \note if you use this constructor you have to make sure com stays alive + */ + GhostLastScalarProduct (const Comm& com, SolverCategory::Category cat) + : GhostLastScalarProduct(stackobject_to_shared_ptr(com), cat) + {} + + /*! \brief Dot product of two vectors. + * \param vx first input vector + * \param vy second input vector + */ + double dot (const Vector& vx, const Vector& vy) const override + { + + // access underlying data + double const *x = &vx[0][0]; + double const *y = &vy[0][0]; + + // total array length + int NN = block_size*count_; + + auto cc = _communication->communicator(); + return cc.sum(vec_dot(x,y,NN)); + } + + /*! \brief Vector L2-norm. + * \param vx input vector + */ + double norm (const Vector& vx) const override + { + return sqrt(dot(vx,vx)); + } + + //! Category of the scalar product (see SolverCategory::Category) + virtual SolverCategory::Category category() const override + { + return _category; + } + + private: + + ///Exctract block size from vector type + static constexpr auto block_size = Vector::block_type::dimension; + + std::shared_ptr _communication; + SolverCategory::Category _category; + int count_; + + /*! \brief Count number of local cells. + */ + int getLocalCount() const + { + int count = 0; + // Loop over index set + auto indexSet = _communication->indexSet(); + for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) { + if (idx->local().attribute()==1) count++; // count non-local indices + } + return count; + } + + /*! \brief Infer number of local cells from largest local index. + */ + int verifyLocalCount() const + { + auto indexSet = _communication->indexSet(); + + size_t is = 0; + // Loop over index set + for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) { + //Only take "owner" indices + if (idx->local().attribute()==1) { + //get local index + auto loc = idx->local().local(); + // if loc is higher than "old interior size", update it + if (loc > is) { + is = loc; + } + } + } + return is + 1; //size is plus 1 since we start at 0 + } + +}; + + + +/// A sequential scalar product that ignores block structure of the vector +/// to facilitate well-known optimization techniques +template +class SeqOptmizedProduct : public Dune::SeqScalarProduct +{ +public: + + /*! \brief Dot product of two vectors. + * \param vx first input vector + * \param vy second input vector + */ + double dot(const Vector& vx, const Vector& vy) const override + { + // access underlying data + double const *x = &vx[0][0]; + double const *y = &vy[0][0]; + + // total array length + int NN = block_size*vx.N(); + + return vec_dot(x,y,NN); + } + + /*! \brief Vector L2-norm. + * \param vx input vector + */ + double norm(const Vector& vx) const override { + return std::sqrt(this->dot(vx, vx)); + } + + private: + + // extract block size + static constexpr auto block_size = Vector::block_type::dimension; + + +}; + +} + +#endif //OPM_SCALAR_PRODUCTS_HEADER_INCLUDED + diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index f11c080b72f..c145e44f0bb 100644 --- a/opm/simulators/linalg/WellOperators.hpp +++ b/opm/simulators/linalg/WellOperators.hpp @@ -262,6 +262,8 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator const matrix_type& getmat() const override { return A_; } + const LinearOperatorExtra& getwellOper() const { return wellOper_; } + void addWellPressureEquations(PressureMatrix& jacobian, const X& weights, const bool use_well_weights) const @@ -357,6 +359,8 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator& getwellOper() const { return wellOper_; } + void addWellPressureEquations(PressureMatrix& jacobian, const X& weights, const bool use_well_weights) const diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 1364fb3a6fd..8ca328332f1 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -13,7 +13,6 @@ namespace Opm //! operations are performed in double-precision //! //! @tparam Vector the block-vector used by linear operator -//! @tparam b block size template class MixedMatrixWrapper { @@ -22,18 +21,12 @@ class MixedMatrixWrapper // extract block size static constexpr auto block_size = Vector::block_type::dimension; - virtual void mv(const Vector& x, Vector& y) const; - virtual void umv(const Vector& x, Vector& y) const; - virtual void usmv(double alpha, const Vector& x, Vector& y) const; - //! @brief constructor //! //! @param nrows number of block rows //! @param nnz number of nonzero blocks MixedMatrixWrapper(int nrows, int nnz) { - if constexpr(block_size!=3) OPM_THROW(std::invalid_argument, "MixedMatrixWrapper only supports block size == 3! \n"); - nnz_=nnz; M_ = bsr_alloc(); bsr_init(M_, nrows, nnz, block_size); @@ -42,40 +35,177 @@ class MixedMatrixWrapper //! @brief destructor ~MixedMatrixWrapper() {bsr_free(M_);} + //! @brief update matrix entries + //! + //! @note downcasts from double precision and transposes + //! each non-zero block entry + //! + //! @param data pointer to double precision data void update(double const *data); + //! @brief block-sparse matrix-vector multiplication (y = M.x) + //! + //! @param x input vector + //! @param y output vector + void mv(const Vector& x, Vector& y) const; + + //! @brief block-sparse matrix-vector multiplication with + //! update (y += M.x) + //! + //! @param x input vector + //! @param y output vector + void umv(const Vector& x, Vector& y) const; + + //! @brief block-sparse matrix-vector multiplication with + //! scaled update (y += alpha * M.x) + //! + //! @param alpha scaling factor + //! @param x input vector + //! @param y output vector + void usmv(double alpha, const Vector& x, Vector& y) const; + + //! @brief access row offset pointer int *rowptr(){return M_->rowptr;} + + //! @brief access column index pointer int *colidx(){return M_->colidx;} + int nrows() const {return M_->nrows;} + private: int nnz_; bsr_matrix *M_; }; +//! @brief mixed-precision block-sparse matrix-vector multiplication +//! (y = M.x) +//! +//! @note hand-optimized versions are provided for block-sizes +//! 2,3, and 4. A generic implementation is provided for block- +//! sizes > 4 +//! +//! @param x input vector +//! @param y output vector template void MixedMatrixWrapper:: mv(const Vector& x, Vector& y) const { - // mixed-precision block spmv (y = M.x) - bsr_vmspmv3(M_, &x[0][0], &y[0][0]); + int const 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 bb = b*b; + double yy[bb]; + for(int i=0;i 4 is +//! NOT provided +//! +//! @param x input vector +//! @param y output vector template void MixedMatrixWrapper:: umv(const Vector& x, Vector& y) const { - // mixed-precision block spmv with update (y += M.x) - bsr_vmspumv3(M_, &x[0][0], &y[0][0], 1.0); + int const 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");} } +//! @brief mixed-precision block-sparse matrix-vector multiplication +//! with scaled update (y += alpha * M.x) +//! +//! @note hand-optimized versions are provided for block-sizes +//! 2,3, and 4. A generic implementation is provided for block- +//! sizes > 4 +//! +//! @param alpha scaling factor +//! @param x input vector +//! @param y output vector template void MixedMatrixWrapper:: usmv(double alpha, const Vector& x, Vector& y) const { - // scaled mixed-precision block spmv with update (y += alpha * M.x) - bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha); + int const 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 bb = b*b; + double yy[bb]; + for(int i=0;i void MixedMatrixWrapper:: update(double const *data) @@ -87,11 +217,10 @@ update(double const *data) for(int k=0;kdbl[bb*k + i] = B[i]; + for(int i=0;iflt[bb*k + i] = B[i]; } - - // downcast to single precision - bsr_downcast(M_); } + + } // namespace Opm #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED diff --git a/opm/simulators/linalg/mixed/Operators.hpp b/opm/simulators/linalg/mixed/Operators.hpp new file mode 100644 index 00000000000..62a5d2de804 --- /dev/null +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -0,0 +1,156 @@ +#ifndef OPM_MIXED_OPERATORS_HEADER_INCLUDED +#define OPM_MIXED_OPERATORS_HEADER_INCLUDED + + +namespace Opm +{ + +//! @brief Adapter to take advantage of the fact that all matrix rows +//! associated with ghost cells are located at the end of the matrix +//! +//! @note The underlying mixed-matrix already ignores ghost rows. +//! +//! @param M matrix class +//! @param V vector class +//! @param C communicator class +template +class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + //! constructor: just store a reference to matrix and communicator + MixedGhostLastMatrixAdapter (const M& A, const C& comm) : A_( A ), comm_(comm) {} + + // y = A * x + void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + ghostLast_project(y); + } + + // y += \alpha * A * x + void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + ghostLast_project(y); + } + + // accessor to matix object + const M& getmat() const override { return A_; } + + // solver category + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::overlapping; + } + +private: + // extract block size + static constexpr auto block_size = V::block_type::dimension; + + void ghostLast_project( V& y ) const + { + double *yy = &y[0][0]; + int n = block_size*A_.nrows(); + int N = block_size*y.N(); + for (int i=n;i +class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + using field_type = typename V::field_type; + using PressureMatrix = Dune::BCRSMatrix>; + + // extract block size + static constexpr auto block_size = V::block_type::dimension; + + //! constructor: just store a reference to a matrix + WellModelMixedGhostLastMatrixAdapter (const M& A, + const LinearOperatorExtra& wellOper, + const C& comm + ) + : A_( A ), wellOper_( wellOper ), comm_ ( comm ) + {} + + // y = A * x + void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + wellOper_.apply(x, y); + //comm_.project(y); + ghostLast_project(y); + } + + // y += \alpha * A * x + void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + wellOper_.applyscaleadd(alpha, x, y); + //comm_.project(y); + ghostLast_project(y); + } + + // accessor to matix object + const M& getmat() const override { return A_; } + + void addWellPressureEquations(PressureMatrix& jacobian, + const V& weights, + const bool use_well_weights) const + { + OPM_TIMEBLOCK(addWellPressureEquations); + wellOper_.addWellPressureEquations(jacobian, weights, use_well_weights); + } + + void addWellPressureEquationsStruct(PressureMatrix& jacobian) const + { + OPM_TIMEBLOCK(addWellPressureEquationsStruct); + wellOper_.addWellPressureEquationsStruct(jacobian); + } + + int getNumberOfExtraEquations() const + { + return wellOper_.getNumberOfExtraEquations(); + } + + // solver category + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::overlapping; + } + + +protected: + + const M& A_ ; + const LinearOperatorExtra& wellOper_; + const C& comm_ ; + +private: + + void ghostLast_project( V& y ) const + { + double *yy = &y[0][0]; + int n = block_size*A_.nrows(); + int N = block_size*y.N(); + for (int i=n;i mixed_matrix_ = bsr_alloc(); bsr_init(mixed_matrix_, nrows, nnz, block_size); - // copy sparsity pattern from double preccision matrix + // copy sparsity pattern from double-precision matrix int *rows = mixed_matrix_->rowptr; int *cols = mixed_matrix_->colidx; @@ -70,12 +70,25 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate prec_free(prec_); } + //! @brief Update ilu0/dilu factorization + //! + //! Transposes double-precision blocks before factorization. + //! Demotes factors after factorization void update() override; + + //! @brief Mixed-precision ilu0/dilu application + //! + //! @param y input vector + //! @param x output vector + void apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) override; + + //! @brief Solver category + Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; }; + bool hasPerfectUpdate() const override {return true;} + void pre ([[maybe_unused]] X& x, [[maybe_unused]] Y& y) override {}; void post ([[maybe_unused]] X& x) override {}; - void apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) override; - Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; }; private: bool use_dilu_; @@ -83,33 +96,382 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate bsr_matrix *mixed_matrix_; prec_t *prec_; int nnz_; + + + //! @brief Dense mixed-precision matrix-vector multiplication + //! (y = A.x) + //! + //! @param y output vector + //! @param A column-major matrix + //! @param x input vector + void matvec_mul(double *y, float const *A, double const * x); + + //! @brief Dense mixed-precision matrix-vector multiply-subtract + //! (y -= A.x) + //! + //! @param y output vector + //! @param A column-major matrix + //! @param x input vector + void matvec_mulsub(double *y, float const *A, double const * x); + + //! @brief Dense matrix copy (C = A) + //! + //! @param C output matrix + //! @param A input matrix + void mat_copy(double *C, double const * A); + + //! @brief Dense matrix inverse + //! (invA = A^{-1}) + //! + //! @param invA output matrix + //! @param A input matrix + void mat_inv(double *invA, const double *A); + + //! @brief Dense matrix-matrix multiply-subtract (C -= A.B) + //! + //! @param C column-major output matrix + //! @param A left column-major input matrix + //! @param B right column-major input matrix + void mat_mulsub(double *C, double const *A, double const * B); + + //! @brief In-place matrix-matrix multiplication (C = C.A) + //! + //! @param C column-major input/output matrix + //! @param A left column-major input matrix + void mat_rmul(double *C, double const *A); + + //! @brief In-place matrix-matrix multiplication (C = A.C) + //! + //! @param C column-major input/output matrix + //! @param A right column-major input matrix + void mat_lmul(double const *A, double *C); }; +//! @brief Update ilu0/dilu factorization +//! +//! Transposes double-precision blocks before factorization. +//! Demotes factors after factorization +//! +//! @note hand-optimized versions are provided for block-sizes +//! 2,3, and 4. A generic implementation is provided for block- +//! sizes > 4 template void MixedPreconditioner:: update () { // transpose each dense block to make them column-major - constexpr int b = block_size; - constexpr int bb=b*b; - double B[bb]; + constexpr int N = block_size; + constexpr int NN=N*N; + + double B[NN]; for(int k=0;kdbl[bb*k + i] = B[i]; + for(int i=0;idbl[NN*k + i] = B[i]; + } + + if constexpr(N==1){OPM_THROW(std::invalid_argument, "MixedMatrixPreconditioner::update does not support block size == 1!\n");} + else if constexpr(N==2) prec_ilu0_factorize2(prec_, mixed_matrix_, use_dilu_); + else if constexpr(N==3) prec_ilu0_factorize3(prec_, mixed_matrix_, use_dilu_); + else if constexpr(N==4) prec_ilu0_factorize4(prec_, mixed_matrix_, use_dilu_); + else + { + bsr_matrix const *A = mixed_matrix_; + bsr_matrix *L=prec_->L; + bsr_matrix *D=prec_->D; + bsr_matrix *U=prec_->U; + + int const nrows = A->nrows; + + // Splitting values of A into L, D, and U, respectively + int kU=0; + for(int i=0;irowptr[i];krowptr[i+1];k++) + { + int j=A->colidx[k]; + if(jrowptr[j]; + mat_copy(L->dbl + NN*kL, A->dbl + NN*k); + L->rowptr[j]++; + } + else if(j==i) // struct-copy of D + { + mat_copy(D->dbl + NN*i, A->dbl + NN*k); + } + else if(j>i) // struct-copy of U + { + mat_copy(U->dbl + NN*kU, A->dbl + NN*k); + kU++; + } + } + } + // reset rowptr of L + for(int i=nrows;i>0;i--) L->rowptr[i]=L->rowptr[i-1]; + L->rowptr[0]=0; + + // Factorizing + int idx=0; + int next = prec_->offsets[idx][0]; + double scale[NN]; + for(int i=0;inrows;i++) + { + mat_inv(scale,D->dbl+i*NN); + mat_copy(D->dbl+NN*i, scale); //store inverse instead to simplify application + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + //scale column i of L + mat_rmul(L->dbl+k*NN,scale); + + //update diagonal D + int j=L->colidx[k]; + mat_mulsub(D->dbl+j*NN,L->dbl+k*NN,U->dbl+k*NN); + } + + if (!use_dilu_) + while(nextrowptr[i+1]) + { + int ij = prec_->offsets[idx][0]; + int ik = prec_->offsets[idx][1]; + int jk = prec_->offsets[idx][2]; + + //update off-diagonals L and U + mat_mulsub(U->dbl+jk*NN,L->dbl+ij*NN,U->dbl+ik*NN); + mat_mulsub(L->dbl+jk*NN,L->dbl+ik*NN,U->dbl+ij*NN); + + //update marker + next=prec_->offsets[++idx][0]; + } + + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + //scale row i of U + mat_lmul(scale,U->dbl+k*NN); + } + } + //prec_test(); getchar(); } - use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0 prec_downcast(prec_); } +//! @brief Mixed-precision ilu0/dilu application +//! +//! @param y input vector +//! @param x output vector +//! +//! @note hand-optimized versions are provided for block-sizes +//! 2,3, and 4. A generic implementation is provided for block- +//! sizes > 4 template void MixedPreconditioner:: apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) { x=y; - prec_mapply3c(prec_,&x[0][0]); + + int const b = block_size; + if constexpr(b==1){OPM_THROW(std::invalid_argument, "MixedMatrixPreconditioner::apply does not support block size == 1!\n");} + else if constexpr(b==2) prec_mapply2c(prec_,&x[0][0]); + else if constexpr(b==3) prec_mapply3c(prec_,&x[0][0]); + else if constexpr(b==4) prec_mapply4c(prec_,&x[0][0]); + else //if constexpr(b==4) + { + bsr_matrix const *L = prec_->L; + bsr_matrix const *D = prec_->D; + bsr_matrix const *U = prec_->U; + + int const N = block_size; + int const NN = N*N; + + // Lower triangular solve assuming ones on diagonal + for(int i=0;incols;i++) + { + double *xi = &x[0][0]+N*i; + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + const float *A = L->flt+k*NN; + int j=U->colidx[k]; // should be L + double *xj = &x[0][0]+N*j; + matvec_mulsub(xj,A,xi); + } + + // Muliply by (inverse) diagonal block + const float *A = D->flt+i*NN; + matvec_mul(xi,A,xi); + } + + // Upper triangular solve assuming ones on diagonal` + for(int i=U->ncols;i>0;i--) + { + double *xi = &x[0][0]+N*(i-1); + for(int k=U->rowptr[i]-1;k>U->rowptr[i-1]-1;k--) + { + const float *A = U->flt+k*NN; + int j=U->colidx[k]; + double const *xj =&x[0][0]+N*j; + matvec_mulsub(xi,A,xj); + } + } + + } +} + +//! @brief Dense mixed-precision matrix-vector multiplication +//! (y = A.x) +//! +//! @param y output vector +//! @param A column-major matrix +//! @param x input vector +template +void MixedPreconditioner:: +matvec_mul(double *y, float const *A, double const * x) +{ + int const N = block_size; + double z[N]; + for(int i=0;i +void MixedPreconditioner:: +matvec_mulsub(double *y, float const *A, double const * x) +{ + int const N = block_size; + double z[N]; + for(int i=0;i +void MixedPreconditioner:: +mat_copy(double *C, double const * A) +{ + int const N = block_size; + int const NN =N*N; + for(int i=0;i +void MixedPreconditioner:: +mat_inv(double *invA, const double *A) +{ + int const N = block_size; + int const NN =N*N; + double T[NN]; + mat_copy(T,A); + + for(int k=0;k +void MixedPreconditioner:: +mat_mulsub(double *C, double const *A, double const * B) +{ + int const N = block_size; + double z[N]; + for(int j=0;j +void MixedPreconditioner:: +mat_rmul(double *C, double const *A) +{ + int const N = block_size; + int const NN =N*N; + double T[NN]; + for(int j=0;j +void MixedPreconditioner:: +mat_lmul(double const *A, double *C) +{ + int const N = block_size; + double z[N]; + for(int j=0;j 1. However, only block-sizes 2,3, and 4 benefit from hand-optimized implementations, +and suboptimal performance is expected for other block sizes. To run the simulator with +mixed-precision ILU0+BiCGSTAB, you can modify the wrapper script below to your liking ``` bash -OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ +OMP_NUM_THREADS=1 mpirun -np 8 --map-by l3cache --bind-to core build/bin/flow \ --matrix-add-well-contributions=true \ --linear-solver=mixed-ilu0 \ --linear-solver-reduction=1e-3 \ --linear-solver-max-iter=1024 \ $@ ``` +Similarly, a sample wrapper for running the simulator with mixed-precision CPR+AMG++BiCGSTAB +is given by +``` bash +OMP_NUM_THREADS=1 mpirun -np 8 --map-by l3cache --bind-to core build/bin/flow \ + --matrix-add-well-contributions=false \ + --linear-solver=mixed-cprw \ + --linear-solver-reduction=1e-3 \ + --linear-solver-max-iter=1024 \ + $@ +``` +Fine-tuning CPR+AMG can be done via a JSON specification file, e.g. by using the wrapper script +below +``` +OMP_NUM_THREADS=1 mpirun -np 8 --map-by l3cache --bind-to core build/bin/flow \ + --linear-solver=../mixed-cprw.json \ + $@ +``` +and modifying the following `mixed-cprw.json` file to your liking +``` +{ + "maxiter": "1024", + "tol": "0.001", + "verbosity": "0", + "solver": "mixed-bicgstab", + "preconditioner": { + "type": "cprw", + "use_well_weights": "false", + "add_wells": "true", + "weight_type": "trueimpes", + "pre_smooth": "0", + "post_smooth": "1", + "finesmoother": { + "type": "mixed-ilu0", + "relaxation": "1" + }, + "verbosity": "0", + "coarsesolver": { + "maxiter": "1", + "tol": "0.10000000000000001", + "solver": "loopsolver", + "verbosity": "0", + "preconditioner": { + "type": "amg", + "alpha": "0.33333333333300003", + "relaxation": "1", + "iterations": "1", + "coarsenTarget": "1200", + "pre_smooth": "1", + "post_smooth": "1", + "beta": "0", + "smoother": "ilu0", + "verbosity": "0", + "maxlevel": "15", + "skip_isolated": "0", + "accumulate": "1", + "prolongationdamping": "1", + "maxdistance": "2", + "maxconnectivity": "15", + "maxaggsize": "6", + "minaggsize": "4" + } + } + } +} +``` + +The legacy mixed-precision implementation is still available. Unlike the current +implementation, it does not leverage the ISTL-framework and consequently only works in serial. +It is considered a developer option and is a little faster than the serial version of the +ISTL-based algorithms documented above. Note that only block sizes 2, 3 and 4 are currently +supported. To run the simulator with legacy mixed-precision ILU0+BiCGSTAB, you can modify the +wrapper script below to your liking +``` bash +OMP_NUM_THREADS=1 mpirun -np 1 --map-by l3cache --bind-to core build/bin/flow \ + --matrix-add-well-contributions=true \ + --linear-solver=legacy-mixed-ilu0 \ + --linear-solver-reduction=1e-3 \ + --linear-solver-max-iter=1024 \ + $@ +``` -To invoke the original serial implementation add the `legacy-` prefix to the mixed-precision -linear solver options. +Have fun! diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 295889b9393..83a8c6b4a5a 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -9,84 +9,17 @@ #include #include #include - #include #include - +#include +#include namespace Dune { -//#include - - -//! @brief Optimized sequential scalar product. -//! -//! @tparam Vector block-vector class with data stored as contiguous double array -template -class SeqOptmizedProduct : public Dune::SeqScalarProduct -{ -public: - - // extract block size - static constexpr auto block_size = Vector::block_type::dimension; - - // Compute the dot product - double dot(const Vector& vx, const Vector& vy) const override - { - // access underlying data - double const *x = &vx[0][0]; - double const *y = &vy[0][0]; - - // total array length - int NN = block_size*vx.N(); - - // unroll loop in multiples of 8 - int n=NN/8; - int N=8*n; - double agg[8]; - for(int i=0;i<8;i++) agg[i]=0.0; - for(int i=0;idot(x, x)); - } -}; - -//! @brief Generalized mixed precision operator interface -//! -//! @tparam Matrix the block-matrix used by linear operator -//! @tparam Vector the block-vector used by linear operator -//! @tparam Comm the communicator used by linear operator -template -struct MixedOperator -{ - using type = Dune::OverlappingSchwarzOperator; -}; -//! @brief Generalized mixed precision operator interface -//! -//! @tparam Matrix the block-matrix used by linear operator -//! @tparam Vector the block-vector used by linear operator -template -struct MixedOperator -{ - using type = Dune::MatrixAdapter; -}; - - -//! @brief Wraps mixed precision +//! @brief Adapts BiCGSTAB to mixed precision //! //! @tparam Comm the communicator passed to FlexibleLinearSolver //! @tparam Operator the linear operator passed to FlexibleLinearSolver @@ -101,8 +34,6 @@ class MixedBiCGSTABSolver:public InverseOperator static constexpr auto block_size = Vector::block_type::dimension; using MixedMatrixType = Opm::MixedMatrixWrapper; - using MixedOperatorType = MixedOperator::type; - using OptimizedProductType = SeqOptmizedProduct; //! @brief constructor //! @@ -121,11 +52,38 @@ class MixedBiCGSTABSolver:public InverseOperator const int& verbosity, const Comm &comm) { + int halo; + size_t nrows; + int nnz=0; - // Access matrix data from double precision operator auto &A = op->getmat(); - int nrows = A.N(); - int nnz = A.nonzeroes(); + // trivially determine size of halo==0 for serial linear operators + if constexpr (std::is_same_v) + { + halo = 0; + nrows = A.N(); + nnz = A.nonzeroes(); + } + // Determine size of halo for parallel linear operators + else + { + local_ = new int[A.N()]; + + // number of ghost cells + halo = getHaloCount(comm); + + // number of local cells + nrows = A.N() - halo; + + // number of nonzeros for local cells + int irow=0; + for(auto row=A.begin(); row.index() < nrows; row++) + { + nnz += local_[irow++] ? std::distance(row->begin(), row->end()) : 0; + } + } + + // Access matrix data from double precision operator double_data_ = &A[0][0][0][0]; //allocate mixed matrix @@ -138,7 +96,7 @@ class MixedBiCGSTABSolver:public InverseOperator int irow = 0; int icol = 0; rows[0] = 0; - for(auto row=A.begin(); row!=A.end(); row++) + for(auto row=A.begin(); row.index() < nrows; row++) { for(auto col = row->begin(); col != row->end(); ++col) { @@ -148,32 +106,59 @@ class MixedBiCGSTABSolver:public InverseOperator irow++; } - // The following skeleton is in preparation for better support for various MatrixAdapter. For now, it simply throws an error if - // the operator provided is not supported. + // initialize mixed operator and scalar product depending on the linear operator type provided to the constructor + double_operator_ = op; using MatrixType = std::remove_const_tgetmat())>>; + + // serial runs with plain block-sparse matrices, i.e. Dune::MatrixAdapter if constexpr (std::is_same_v, Dune::MatrixAdapter>) { + using MixedOperatorType = Dune::MatrixAdapter; + mixed_operator_ = std::make_shared(*mixed_matrix_); + + using ScalarProductType = SeqOptmizedProduct; + scalar_product_ = std::make_shared(); } - else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) + // serial runs with separate linear operator for wells, i.e. Opm::WellModelMatrixAdapter + else if constexpr (std::is_same_v, Opm::WellModelMatrixAdapter>) { + using MixedOperatorType = Opm::WellModelMatrixAdapter; + using WellOperatorType = Opm::LinearOperatorExtra; + const WellOperatorType &wellOper = op->getwellOper(); + mixed_operator_ = std::make_shared(*mixed_matrix_, wellOper); + + using ScalarProductType = SeqOptmizedProduct; + scalar_product_ = std::make_shared(); } - else + // parallel runs with plain block-sparse matrices and all ghost cells sorted after local cells, i.e. Opm::GhostLastMatrixAdapter + else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) { - OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver only supports Dune::MatrixAdapter and Opm::GhostLastMatrixAdapter\n"); - } + using MixedOperatorType = Opm::MixedGhostLastMatrixAdapter; + mixed_operator_ = std::make_shared(*mixed_matrix_,comm); - //initialize mixed operator and optimized scalar product - double_operator_ = op; - if constexpr (std::is_same_v) - { - mixed_operator_ = std::make_shared(*mixed_matrix_); - scalar_product_ = std::make_shared(); + using ScalarProductType = GhostLastScalarProduct; + scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); } - else + // parallel runs with separate linear operators for wells and all ghost cells sorted after local cells, i.e. Opm::WellModelGhostLastMatrixAdapter + else if constexpr (std::is_same_v, Opm::WellModelGhostLastMatrixAdapter>) { - mixed_operator_ = std::make_shared(*mixed_matrix_,comm); - scalar_product_ = sp; + using MixedOperatorType = Opm::WellModelMixedGhostLastMatrixAdapter; + using WellOperatorType = Opm::LinearOperatorExtra; + const WellOperatorType &wellOper = op->getwellOper(); + mixed_operator_ = std::make_shared(*mixed_matrix_, wellOper, comm); + + if constexpr (std::is_same_v) + { + scalar_product_ = sp; + } + else + { + using ScalarProductType = GhostLastScalarProduct; + scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); + } } + // throw an exception for all other linear operator types + else { OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver: Unsupported linear operator type!!\n");} //initialize bicgstab solver from Dune solver_ = std::make_shared>( @@ -183,8 +168,17 @@ class MixedBiCGSTABSolver:public InverseOperator tol, // desired residual reduction factor maxiter, // maximum number of iterations verbosity); + + } + + //! @brief destructor + ~MixedBiCGSTABSolver() + { + if constexpr (std::is_same_v) return; + delete [] local_; } + //! @brief Solver application void apply(Vector &x, Vector &b, InverseOperatorResult &res) override { //transpose dense blocks and demote to single precision @@ -194,6 +188,7 @@ class MixedBiCGSTABSolver:public InverseOperator solver_->apply(x,b,res); } + //! @brief Unused variant of solver application void apply(Vector &x, Vector &b, double reduction, InverseOperatorResult &res) override { x=0; @@ -202,19 +197,41 @@ class MixedBiCGSTABSolver:public InverseOperator OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver::apply(...) not implemented yet."); } + //! @brief Solver category Dune::SolverCategory::Category category() const override{return Dune::SolverCategory::overlapping;}; private: - using AbstractSolverType = Dune::InverseOperator; + //! @brief Count number of ghost cells + //! + //! @param comm communicator object + int getHaloCount(const Comm& comm) const + { + int count = 0; + // Loop over index set + auto indexSet = comm.indexSet(); + for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) + { + if (idx->local().attribute()!=1) count++; // count ghost indices + + int i=idx->local().local(); // tag local indices + local_[i] = (idx->local().attribute()==1) ? 1 : 0; + } + + return count; + } + + using AbstractSolverType = Dune::InverseOperator; + using AbstractOperatorType = Dune::AssembledLinearOperator; Operator *double_operator_; std::shared_ptr solver_; - std::shared_ptr mixed_operator_; + std::shared_ptr mixed_operator_; std::shared_ptr mixed_matrix_; std::shared_ptr scalar_product_; double const *double_data_; + int *local_; }; } diff --git a/opm/simulators/linalg/mixed/bslv.c b/opm/simulators/linalg/mixed/bslv.c index 77ba71b6b52..9d1dd7854ec 100644 --- a/opm/simulators/linalg/mixed/bslv.c +++ b/opm/simulators/linalg/mixed/bslv.c @@ -70,6 +70,30 @@ void bslv_init(bslv_memory *mem, double tol, int max_iter, bsr_matrix const *A, mem->P = prec_alloc(); prec_init(mem->P, A); // initialize structure of L,D,U components of P + + // pick spmv, factorization, and apply functions according to block size + switch (A->b) + { + case 2: + mem->bsr_spmv = bsr_vmspmv2; + mem->prec_factorize = prec_ilu0_factorize2; + mem->prec_apply = prec_mapply2c; + break; + case 3: + mem->bsr_spmv = bsr_vmspmv3; + mem->prec_factorize = prec_ilu0_factorize3; + mem->prec_apply = prec_mapply3c; + break; + case 4: + mem->bsr_spmv = bsr_vmspmv4; + mem->prec_factorize = prec_ilu0_factorize4; + mem->prec_apply = prec_mapply4c; + break; + default: + mem->bsr_spmv = NULL; + mem->prec_factorize = NULL; + mem->prec_apply = NULL; + } } /** @@ -101,9 +125,8 @@ double __attribute__((noinline)) vec_inner2(const double *a, const double *b, in return agg[0]; } -int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x) +int bslv_pbicgstabm(bslv_memory *mem, bsr_matrix *A, const double *b, double *x) { - double tol = mem->tol; int max_iter = mem->max_iter; int n = mem->n; @@ -120,7 +143,7 @@ int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x double * restrict x_j = x; prec_t * restrict P = mem->P; - mem->use_dilu ? prec_dilu_factorize(P,A) : prec_ilu0_factorize(P,A); // choose dilu or ilu0 + mem->prec_factorize(P,A,mem->use_dilu); // choose dilu or ilu0 prec_downcast(P); vec_fill(x_j,0.0,n); @@ -128,54 +151,45 @@ int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x vec_copy(p_j,b,n); vec_copy(q_j,p_j,n); - //double norm_0 = sqrt(vec_inner(r_j,r_j,n)); double norm_0 = sqrt(vec_inner2(r_j,r_j,n)); - //double rho_j = vec_inner(r0,r_j,n); double rho_j = vec_inner2(r0,r_j,n); int j; for(j=0;jprec_apply(P,q_j); //q_j=P.q_j; + mem->bsr_spmv(A,q_j,v_j); //v_j= A.q_j - //double alpha_j = rho_j/vec_inner(r0,v_j,n); double alpha_j = rho_j/vec_inner2(r0,v_j,n); - for (int k=0;kprec_apply(P,q_j); //q_j=P.q_j; + mem->bsr_spmv(A,q_j,t_j); //t_j= A.q_j - //double w_j = vec_inner(s_j,t_j,n)/vec_inner(t_j,t_j,n); double w_j = vec_inner2(s_j,t_j,n)/vec_inner2(t_j,t_j,n); for (int k=0;kprec_apply(P,x_j); //x_j=P.x_j; return j == max_iter ? j : ++j; } - int bslv_pbicgstab3d(bslv_memory *mem, bsr_matrix *A, const double *b, double *x) { - double tol = mem->tol; int max_iter = mem->max_iter; int n = mem->n; @@ -192,7 +206,7 @@ int bslv_pbicgstab3d(bslv_memory *mem, bsr_matrix *A, const double *b, double *x double * restrict x_j = x; prec_t * restrict P = mem->P; - mem->use_dilu ? prec_dilu_factorize(P,A) : prec_ilu0_factorize(P,A); // choose dilu or ilu0 + prec_ilu0_factorize3(P,A,mem->use_dilu); // choose dilu or ilu0 prec_downcast(P); vec_fill(x_j,0.0,n); @@ -200,52 +214,41 @@ int bslv_pbicgstab3d(bslv_memory *mem, bsr_matrix *A, const double *b, double *x vec_copy(p_j,b,n); vec_copy(q_j,p_j,n); - //double norm_0 = sqrt(vec_inner(r_j,r_j,n)); double norm_0 = sqrt(vec_inner2(r_j,r_j,n)); - //double rho_j = vec_inner(r0,r_j,n); double rho_j = vec_inner2(r0,r_j,n); int j; for(j=0;jnrows=nrows; @@ -47,8 +55,9 @@ void bsr_init(bsr_matrix *A, int nrows, int nnz, int b) A->rowptr = malloc((nrows+1)*sizeof(int)); A->colidx = malloc(nnz*sizeof(int)); - A->dbl = malloc(b*b*nnz*sizeof(double)); - A->flt = malloc(b*b*nnz*sizeof(float)); + + A->dbl = buffered_alloc(64,b*b*nnz*sizeof(double)); + A->flt = buffered_alloc(64,b*b*nnz*sizeof(float)); assert(A->rowptr); assert(A->colidx); @@ -185,6 +194,141 @@ void bsr_vdspmv3(bsr_matrix *A, const double *x, double *y) } +void bsr_vmspmv4(bsr_matrix *A, const double *x, double *y) +{ + int nrows = A->nrows; + int *rowptr=A->rowptr; + int *colidx=A->colidx; + const float *data=A->flt; + + const int b=4; + + __m256d mm_zeros =_mm256_setzero_pd(); + for(int i=0;inrows; + int *rowptr=A->rowptr; + int *colidx=A->colidx; + const float *data=A->flt; + + const int b=4; + + __m256d valpha = _mm256_set1_pd(alpha); + + __m256d mm_zeros =_mm256_setzero_pd(); + for(int i=0;inrows; + int *rowptr=A->rowptr; + int *colidx=A->colidx; + const float *data=A->flt; + + const int b=2; + + __m256d mm_zeros =_mm256_setzero_pd(); + for(int i=0;inrows; + int *rowptr=A->rowptr; + int *colidx=A->colidx; + const float *data=A->flt; + + const int b=2; + + __m128d valpha = _mm_set1_pd(alpha); + + __m256d mm_zeros =_mm256_setzero_pd(); + for(int i=0;i + +double vec_dot (double const *x, double const *y, int NN) +{ + // unroll loop in multiples of 8 + int n=NN/8; + int N=8*n; + double agg[8]; + for(int i=0;i<8;i++) agg[i]=0.0; + for(int i=0;i #include +#include "matvec.h" prec_t *prec_alloc() { @@ -194,6 +195,11 @@ void mat3_inv(double *invA, const double *A) for(int k=0;k<9;k++) invA[k]=M[k]/detA; } +static inline void vec_copy4(double *y, double const *x) +{ + for(int i=0;i<4;i++) y[i]=x[i]; +} + /** * @brief vector copy of 9-element vectors. * @@ -205,6 +211,11 @@ static inline void vec_copy9(double *y, double const *x) for(int i=0;i<9;i++) y[i]=x[i]; } +static inline void vec_copy16(double *y, double const *x) +{ + for(int i=0;i<16;i++) y[i]=x[i]; +} + /** * @brief In-place right matrix-matrix multiplication for 3x3 matrices. * @@ -314,8 +325,7 @@ void mat3_vfms(double *C, double const *A, double const *B) } } - -void prec_dilu_factorize(prec_t *P, bsr_matrix *A) +void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A, bool use_dilu) { int nrows = A->nrows; int b = A->b; @@ -335,16 +345,16 @@ void prec_dilu_factorize(prec_t *P, bsr_matrix *A) if(jrowptr[j]; - vec_copy9(L->dbl + bb*kL, A->dbl + bb*k); + vec_copy4(L->dbl + bb*kL, A->dbl + bb*k); L->rowptr[j]++; } else if(j==i) // struct-copy of D { - vec_copy9(D->dbl + bb*i, A->dbl + bb*k); + vec_copy4(D->dbl + bb*i, A->dbl + bb*k); } else if(j>i) // struct-copy of U { - vec_copy9(U->dbl + bb*kU, A->dbl + bb*k); + vec_copy4(U->dbl + bb*kU, A->dbl + bb*k); kU++; } } @@ -354,38 +364,47 @@ void prec_dilu_factorize(prec_t *P, bsr_matrix *A) L->rowptr[0]=0; // Factorizing - double scale[9]; //hard-coded to 3x3 blocks for now + int idx=0; + int next = use_dilu ? A->nnz : P->offsets[idx][0]; + double scale[4]; //hard-coded to 2x2 blocks for(int i=0;inrows;i++) { - mat3_inv(scale,D->dbl+i*bb); - vec_copy9(D->dbl+bb*i, scale); //store inverse instead to simplify application + mat2_inv(scale,D->dbl+i*bb); + vec_copy4(D->dbl+bb*i, scale); //store inverse instead to simplify application for(int k=L->rowptr[i];krowptr[i+1];k++) { //scale column i of L - mat3_rmul(L->dbl+k*bb,scale); + mat2_rmul(L->dbl+k*bb,scale); - //update diagonal of U + //update diagonal D int j=L->colidx[k]; - mat3_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); + mat2_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); + } - //scale row i of U - mat3_lmul(scale,U->dbl+k*bb); + while(nextrowptr[i+1]) + { + int ij = P->offsets[idx][0]; + int ik = P->offsets[idx][1]; + int jk = P->offsets[idx][2]; - //NOT IMPLEMENTED! - for(int m=L->rowptr[j];mrowptr[j+1];m++) - { - if(L->colidx[m]==j) - { - printf("ILU OFF_DIAGONALS NOT IMPLEMENTED!\n"); - printf("(%d,%d)",m,j); - getchar(); - } - } + //update off-diagonals L and U + mat2_vfms(U->dbl+jk*bb,L->dbl+ij*bb,U->dbl+ik*bb); + mat2_vfms(L->dbl+jk*bb,L->dbl+ik*bb,U->dbl+ij*bb); + + //update marker + next=P->offsets[++idx][0]; } + + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + //scale row i of U + mat2_lmul(scale,U->dbl+k*bb); + } + } } -void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) +void prec_ilu0_factorize3(prec_t *P, bsr_matrix *A, bool use_dilu) { int nrows = A->nrows; int b = A->b; @@ -425,7 +444,7 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) // Factorizing int idx=0; - int next = P->offsets[idx][0]; + int next = use_dilu ? A->nnz : P->offsets[idx][0]; double scale[9]; //hard-coded to 3x3 blocks for now for(int i=0;inrows;i++) { @@ -464,51 +483,142 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) } } -#if 0 -/** - * @brief In-place matrix-vector multiplication for 3x3 matrices. - * - * @param A Pointer to input matrix. - * @param x Pointer to input/output vector. - */ -static inline void mat3_vecmul(const double *A, double *x) +void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A, bool use_dilu) { - const int b=3; - double z[3]; - for(int k=0;k<3;k++) z[k]=0; - for(int c=0;cnrows; + //int const b = 4; + int const bb =16; + + bsr_matrix *L=P->L; + bsr_matrix *D=P->D; + bsr_matrix *U=P->U; + + // Splitting values of A into L, D, and U, respectively + int kU=0; + for(int i=0;irowptr[i];krowptr[i+1];k++) { - z[r]+=A[c*b+r]*x[c]; + int j=A->colidx[k]; + if(jrowptr[j]; + vec_copy16(L->dbl + bb*kL, A->dbl + bb*k); + L->rowptr[j]++; + } + else if(j==i) // struct-copy of D + { + vec_copy16(D->dbl + bb*i, A->dbl + bb*k); + } + else if(j>i) // struct-copy of U + { + vec_copy16(U->dbl + bb*kU, A->dbl + bb*k); + kU++; + } } } - for(int k=0;k<3;k++) x[k]=z[k]; -} + // reset rowptr of L + for(int i=nrows;i>0;i--) L->rowptr[i]=L->rowptr[i-1]; + L->rowptr[0]=0; -/** - * @brief In-place fused matrix-vector multiply-subtract for 3x3 matrices. - * - * @param y Pointer to input/output vector. - * @param A Pointer to input matrix. - * @param x Pointer to input vector. - */ + // Factorizing + int idx=0; + int next = use_dilu ? A->nnz : P->offsets[idx][0]; + double scale[16] __attribute__((aligned(64))); //hard-coded to 4x4 blocks + for(int i=0;inrows;i++) + { + mat4_vinv(scale,D->dbl+i*bb); + vec_copy16(D->dbl+bb*i, scale); //store inverse instead to simplify application + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + //scale column i of L + mat4_rmul(L->dbl+k*bb,scale); -static inline void mat3_vecfms(double *y, const double *A, const double *x) + //update diagonal D + int j=L->colidx[k]; + mat4_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); + } + + while(nextrowptr[i+1]) + { + int ij = P->offsets[idx][0]; + int ik = P->offsets[idx][1]; + int jk = P->offsets[idx][2]; + + //update off-diagonals L and U + mat4_vfms(U->dbl+jk*bb,L->dbl+ij*bb,U->dbl+ik*bb); + mat4_vfms(L->dbl+jk*bb,L->dbl+ik*bb,U->dbl+ij*bb); + + //update marker + next=P->offsets[++idx][0]; + } + + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + //scale row i of U + mat4_lmul(scale,U->dbl+k*bb); + } + } +} + +void prec_mapply2c(prec_t *restrict P, double *x) { - const int b=3; - double z[3]; - for(int k=0;k<3;k++) z[k]=0; - for(int c=0;cL; + bsr_matrix *D = P->D; + bsr_matrix *U = P->U; + + int b=L->b; + int bb=b*b; + + __m256d mm256_zero_pd =_mm256_setzero_pd(); + + // Lower triangular solve assuming ones on diagonal + for(int i=0;incols;i++) { - for(int r=0;rrowptr[i];krowptr[i+1];k++) { - z[r]+=A[b*c+r]*x[c]; + const float *A = L->flt+k*bb; + int j=U->colidx[k]; // should be L, but does not matter due to structural symmetry? + vA = _mm256_cvtps_pd(_mm_loadu_ps(A))*vx; + + double *xj = x+b*j; + __m128d vxj = _mm_loadu_pd(xj) - (_mm256_extractf128_pd(vA,0) +_mm256_extractf128_pd(vA,1)); + _mm_storeu_pd(xj,vxj); } + + // Muliply by (inverse) diagonal block + const float *A = D->flt+i*bb; + vA = _mm256_cvtps_pd(_mm_loadu_ps(A))*vx; + __m128d vz = _mm256_extractf128_pd(vA,0) +_mm256_extractf128_pd(vA,1); + + _mm_storeu_pd(xi,vz); + } + + // Upper triangular solve assuming nonzeros stored in original order + for(int i=U->ncols;i>0;i--) + { + __m256d vA; + vA=mm256_zero_pd; + for(int k=U->rowptr[i]-1;k>U->rowptr[i-1]-1;k--) + { + const float *A = U->flt+k*bb; + int j=U->colidx[k]; + __m256d vxj = _mm256_loadu_pd(x+b*j); + vA += _mm256_cvtps_pd(_mm_loadu_ps(A))*_mm256_permute4x64_pd(vxj,0x50); + } + + double *xi = x+b*(i-1); + __m128d vxi = _mm_loadu_pd(xi) - (_mm256_extractf128_pd(vA,0) +_mm256_extractf128_pd(vA,1)); + _mm_storeu_pd(xi,vxi); } - for(int k=0;k<3;k++) y[k]-=z[k]; } -#endif void prec_mapply3c(prec_t *restrict P, double *x) { @@ -583,6 +693,78 @@ void prec_mapply3c(prec_t *restrict P, double *x) } } +void prec_mapply4c(prec_t *restrict P, double *x) +{ + bsr_matrix *L = P->L; + bsr_matrix *D = P->D; + bsr_matrix *U = P->U; + + int const b=4; + int const bb=16; + + __m256d mm256_zero_pd =_mm256_setzero_pd(); + + // Lower triangular solve assuming ones on diagonal + for(int i=0;incols;i++) + { + __m256d vA[4], vx[4]; + + double *xi = x+b*i; + __m256d vxi = _mm256_loadu_pd(xi); + + vx[0] = _mm256_permute4x64_pd(vxi,0x00); // 0b00000000 + vx[1] = _mm256_permute4x64_pd(vxi,0x55); // 0b01010101 + vx[2] = _mm256_permute4x64_pd(vxi,0xAA); // 0b10101010 + vx[3] = _mm256_permute4x64_pd(vxi,0xFF); // 0b11111111 + for(int k=L->rowptr[i];krowptr[i+1];k++) + { + const float *A = L->flt+k*bb; + int j=U->colidx[k]; // should be L, but does not matter due to structural + vA[0] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 0))*vx[0]; + vA[1] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 4))*vx[1]; + vA[2] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 8))*vx[2]; + vA[3] = _mm256_cvtps_pd(_mm_loadu_ps(A+12))*vx[3]; + + double *xj = x+b*j; + __m256d vxj = _mm256_loadu_pd(xj); + __m256d vz = vxj - (vA[0]+vA[1]) - (vA[2]+vA[3]); + _mm256_storeu_pd(xj,vz); + } + + // Muliply by (inverse) diagonal block + const float *A = D->flt+i*bb; + vA[0] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 0))*vx[0]; //0b00000000 + vA[1] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 4))*vx[1]; //0b01010101 + vA[2] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 8))*vx[2]; //0b10101010 + vA[3] = _mm256_cvtps_pd(_mm_loadu_ps(A+12))*vx[3]; //0b11111111 + + __m256d vz = vA[0] + vA[1] + vA[2] + vA[3]; + _mm256_storeu_pd(xi,vz); + } + + // Upper triangular solve assuming ones on diagonal` + for(int i=U->ncols;i>0;i--) + { + __m256d vA[4]; + for(int k=0;k<4;k++) vA[k]=mm256_zero_pd; + for(int k=U->rowptr[i]-1;k>U->rowptr[i-1]-1;k--) + { + const float *A = U->flt+k*bb; + int j=U->colidx[k]; + __m256d vxj = _mm256_loadu_pd(x+b*j); + vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(A+ 0))*_mm256_permute4x64_pd(vxj,0x00); // 0b00000000 + vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(A+ 4))*_mm256_permute4x64_pd(vxj,0x55); // 0b01010101 + vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(A+ 8))*_mm256_permute4x64_pd(vxj,0xAA); // 0b10101010 + vA[3] += _mm256_cvtps_pd(_mm_loadu_ps(A+12))*_mm256_permute4x64_pd(vxj,0xFF); // 0b11111111 + } + + double *xi = x+b*(i-1); + __m256d vxi = _mm256_loadu_pd(xi); + __m256d vz = vxi - (vA[0]+vA[1]) - (vA[2]+vA[3]); + _mm256_storeu_pd(xi,vz); + } +} + void prec_dapply3c(prec_t *restrict P, double *x) { bsr_matrix *L = P->L; @@ -669,3 +851,4 @@ void prec_info(prec_t *P) bsr_info(P->D); bsr_info(P->U); } + diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index 7d474801925..4c58a528fe5 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -5,6 +5,7 @@ extern "C" { #endif #include "bsr.h" +#include /*! * @brief Preconditioner struct. @@ -58,20 +59,14 @@ void prec_init(prec_t *P, bsr_matrix const *A); int prec_analyze(bsr_matrix *M, int (*offsets)[3]); /** - * @brief DILU factorization. + * @brief ILU0/DILU factorization. * * @param P Pointer preconditioner object. * @param A Pointer to bsr matrix. */ -void prec_dilu_factorize(prec_t *P, bsr_matrix *A); - -/** - * @brief ILU0 factorization. - * - * @param P Pointer preconditioner object. - * @param A Pointer to bsr matrix. - */ -void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); +void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A, bool use_dilu); +void prec_ilu0_factorize3(prec_t *P, bsr_matrix *A, bool use_dilu); +void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A, bool use_dilu); /** * @brief Preconditioner application in mixed-precision. @@ -81,7 +76,9 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); * @param P Pointer to preconditioner object. * @apram x Pointer to input/output vector */ +void prec_mapply2c(prec_t *P, double *x); void prec_mapply3c(prec_t *P, double *x); +void prec_mapply4c(prec_t *P, double *x); /** * @brief Preconditioner applicationin double-precision. diff --git a/opm/simulators/linalg/mixed/wrapper.hpp b/opm/simulators/linalg/mixed/wrapper.hpp index cd4bc2253e8..14d60c9ade9 100644 --- a/opm/simulators/linalg/mixed/wrapper.hpp +++ b/opm/simulators/linalg/mixed/wrapper.hpp @@ -25,8 +25,8 @@ class MixedSolver : public InverseOperator int nnz = A.nonzeroes(); int b = A[0][0].N(); - // verify that block size is 3x3 - if (b!=3) {OPM_THROW(std::logic_error, "Block sizes other than 3x3 are not supported by mixed precision.");} + // verify that block size is 3x3 or 4x4 + if (b<2 || b>4) {OPM_THROW(std::logic_error, "Legacy mixed precision only supports 3x3 and 4x4 blocks.");} // create jacobian matrix object and allocate various arrays jacobian_ = bsr_alloc(); @@ -61,25 +61,25 @@ class MixedSolver : public InverseOperator { bsr_free(jacobian_); bslv_free(mem_); - } void apply (X& x, X& b, InverseOperatorResult& res) override { // transpose each dense block to make them column-major - double B[9]; + int const N = block_size; + int const NN = N*N; + double B[NN]; for(int k=0;knnz;k++) { - for(int i=0;i<3;i++) for(int j=0;j<3;j++) B[3*j+i] = data_[9*k + 3*i + j]; - for(int i=0;i<9;i++) jacobian_->dbl[9*k + i] = B[i]; + for(int i=0;idbl[NN*k + i] = B[i]; } // downcast to allow mixed precision bsr_downcast(jacobian_); // solve linear system - int count = bslv_pbicgstab3m(mem_, jacobian_, &b[0][0], &x[0][0]); - //int count = bslv_pbicgstab3d(mem_, jacobian_, &b[0][0], &x[0][0]); + int count = bslv_pbicgstabm(mem_, jacobian_, &b[0][0], &x[0][0]); // return convergence information res.converged = (mem_->e[count] < mem_->tol); @@ -99,6 +99,10 @@ class MixedSolver : public InverseOperator Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::sequential; }; private: + + // extract block size + static constexpr auto block_size = X::block_type::dimension; + bsr_matrix *jacobian_; bslv_memory *mem_; double const *data_; diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 29047733943..41fd56c2433 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -284,7 +284,7 @@ setupPropertyTree(FlowLinearSolverParameters p, // Note: copying the parameters return setupCPR(conf, p); } - if ((conf == "cpr") || (conf == "cprw")) { + if ((conf == "cpr") || (conf == "cprw") || (conf == "mixed-cprw")) { if (!linearSolverMaxIterSet) { // Use our own default unless it was explicitly overridden by user. p.linear_solver_maxiter_ = cprDefaultMaxIter; @@ -379,14 +379,16 @@ setupCPRW(const std::string& /*conf*/, const FlowLinearSolverParameters& p) prm.put("maxiter", p.linear_solver_maxiter_); prm.put("tol", p.linear_solver_reduction_); prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", getSolverString(p)); + //prm.put("solver", getSolverString(p)); + prm.put("solver", (p.linsolver_ == "mixed-cprw")?"mixed-bicgstab":getSolverString(p)); prm.put("preconditioner.type", "cprw"s); prm.put("preconditioner.use_well_weights", "false"s); prm.put("preconditioner.add_wells", "true"s); prm.put("preconditioner.weight_type", "trueimpes"s); prm.put("preconditioner.pre_smooth", 0); prm.put("preconditioner.post_smooth", 1); - prm.put("preconditioner.finesmoother.type", "paroverilu0"s); + //prm.put("preconditioner.finesmoother.type", "paroverilu0"s); + prm.put("preconditioner.finesmoother.type", (p.linsolver_ == "mixed-cprw")?"mixed-ilu0":"paroverilu0"s); prm.put("preconditioner.finesmoother.relaxation", 1.0); prm.put("preconditioner.verbosity", 0); prm.put("preconditioner.coarsesolver.maxiter", 1); @@ -533,7 +535,7 @@ setupMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverPa prm.put("tol", p.linear_solver_reduction_); prm.put("maxiter", p.linear_solver_maxiter_); prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", "mixed-precision"s); + prm.put("solver", "mixed-bicgstab"s); prm.put("preconditioner.type", "mixed-ilu0"s); return prm; } @@ -546,7 +548,7 @@ setupMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearSolverP prm.put("tol", p.linear_solver_reduction_); prm.put("maxiter", p.linear_solver_maxiter_); prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", "mixed-precision"s); + prm.put("solver", "mixed-bicgstab"s); prm.put("preconditioner.type", "mixed-dilu"s); return prm; } @@ -560,7 +562,7 @@ setupLegacyMixedILU([[maybe_unused]] const std::string& conf, const FlowLinearSo prm.put("tol", p.linear_solver_reduction_); prm.put("maxiter", p.linear_solver_maxiter_); prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", "mixed-bicgstab"s); + prm.put("solver", "mixed-legacy"s); prm.put("preconditioner.type", "legacy-mixed-ilu0"s); return prm; } @@ -573,7 +575,7 @@ setupLegacyMixedDILU([[maybe_unused]] const std::string& conf, const FlowLinearS prm.put("tol", p.linear_solver_reduction_); prm.put("maxiter", p.linear_solver_maxiter_); prm.put("verbosity", p.linear_solver_verbosity_); - prm.put("solver", "mixed-bicgstab"s); + prm.put("solver", "mixed-legacy"s); prm.put("preconditioner.type", "legacy-mixed-dilu"s); return prm; }