diff --git a/opm/simulators/linalg/mixed/README.md b/opm/simulators/linalg/mixed/README.md index 988ff58e91..1e27c55ea8 100644 --- a/opm/simulators/linalg/mixed/README.md +++ b/opm/simulators/linalg/mixed/README.md @@ -1,17 +1,19 @@ # Mixed-precision linear solvers This folder contains mixed-precision building blocks for Krylov subspace methods -and a highly optimized mixed-precision implementation of ILU0 preconditioned bicgstab. -Hopefully, this will inspire the exploration of mixed-precision algorithms in OPM. +and a highly optimized mixed-precision implementation of ILU0 and DILU preconditioned +bicgstab. Hopefully, this will inspire the exploration of mixed-precision algorithms +in OPM. The initial implementation is specialized for 3x3 block-sparse matrices due to their importance in reservoir simulation and restricted to serial runs. Extending the work to parallel runs is a matter of extending OPM's parallel infrastructure and should be relatively straight-forward. -The mixed-precision solver is selected by the command-line option `--linear-solver=mixed-ilu0`. -The command-line option `--matrix-add-well-contributions=true` must also be set as the -mixed-precision solver operates directly on block-sparse matrices, not on linear operators as -other OPM solvers do. For convenience, a wrapper similar to the one below can be used. +The mixed-precision solver is selected by the command-line options `--linear-solver=mixed-ilu0` +or `--linear-solver=mixed-dilu`. The command-line option `--matrix-add-well-contributions=true` +must also be set as the mixed-precision solver operates directly on block-sparse matrices, not +on linear operators as other OPM solvers do. For convenience, a wrapper similar to the one below +can be used. ``` bash OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ diff --git a/opm/simulators/linalg/mixed/bslv.h b/opm/simulators/linalg/mixed/bslv.h index c2b83f8849..6f31d8d010 100644 --- a/opm/simulators/linalg/mixed/bslv.h +++ b/opm/simulators/linalg/mixed/bslv.h @@ -10,19 +10,28 @@ extern "C" { #include -// Solver memory struct +/*! + * @brief Linear solver memory. + */ typedef struct bslv_memory { + // perform dilu (ilu0) of true (False) bool use_dilu; + // tolerance double tol; + // maximum number of iterations int max_iter; + // convergencey history double *e; + // size of linear system int n; + // pointer to temporary arrays double **dtmp; + // pointer to preconditioner prec_t *P; } bslv_memory; diff --git a/opm/simulators/linalg/mixed/bsr.h b/opm/simulators/linalg/mixed/bsr.h index f58a022179..11e72787ea 100644 --- a/opm/simulators/linalg/mixed/bsr.h +++ b/opm/simulators/linalg/mixed/bsr.h @@ -4,18 +4,28 @@ extern "C" { #endif -// bsr matrix struct +/*! + * @brief Mixed-precision bsr matrix. + */ typedef struct bsr_matrix { + // number or rows int nrows; + // number or columns int ncols; + // number or non-zero blocks int nnz; + // block-size int b; + // pointer ot row offsets int *rowptr; + // pointer to column indices int *colidx; + // pointer to double-precision values double *dbl; + // pointer to single-precision values float *flt; } bsr_matrix; diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index f0b9ea6d01..23f11f57c5 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -122,6 +122,13 @@ void prec_init(prec_t *P, bsr_matrix const *A) P->noffsets=count; } +/** + * @brief Matrix-matrix multiplication of 3x3 matrices. + * + * @param C Pointer to output matrix. + * @param A Pointer to left input matrix. + * @param B Pointer to right input matrix. + */ void mat3_matmul(double *C, const double *A, const double *B) { // assume 3x3 column-major matrices @@ -141,6 +148,13 @@ void mat3_matmul(double *C, const double *A, const double *B) for(int k=0;k<9;k++) C[k]=M[k]; } +/** + * @brief Fused multiply-subtract for 3x3 matrices. + * + * @param C Pointer to output matrix. + * @param A Pointer to left input matrix. + * @param B Pointer to right input matrix. + */ void mat3_matfms(double *C, const double *A, const double *B) { // assume 3x3 column-major matrices @@ -160,7 +174,12 @@ void mat3_matfms(double *C, const double *A, const double *B) for(int k=0;k<9;k++) C[k]-=M[k]; } - +/** + * @brief Matrix invese for 3x3 matrix. + * + * @param invA Pointer to inverse matrix. + * @param A Pointer to input matrix. + */ void mat3_inv(double *invA, const double *A) { double M[9]; @@ -178,11 +197,25 @@ void mat3_inv(double *invA, const double *A) for(int k=0;k<9;k++) invA[k]=M[k]/detA; } +/** + * @brief vector copy of 9-element vectors. + * + * @param y Pointer to output vector. + * @param x Pointer to input vector. + */ inline void vec_copy9(double *y, double const *x) { for(int i=0;i<9;i++) y[i]=x[i]; } +/** + * @brief In-place right matrix-matrix multiplication for 3x3 matrices. + * + * @note Function is optimized with avx2 intrinsics + * + * @param A Pointer to left input and output matrix. + * @param B Pointer to right input matrix. + */ void mat3_rmul(double *A, double const *B) { // load left hand matrix @@ -211,6 +244,14 @@ void mat3_rmul(double *A, double const *B) } } +/** + * @brief In-place left matrix-matrix multiplication for 3x3 matrices. + * + * @note Function is optimized with avx2 intrinsics + * + * @param A Pointer to left input matrix. + * @param B Pointer to right input and output matrix. + */ void mat3_lmul(double const *A, double *B) { // load left hand matrix @@ -239,6 +280,15 @@ void mat3_lmul(double const *A, double *B) } } +/** + * @brief Fused multiply-subtract for 3x3 matrices. + * + * @note Function is optimized with avx2 intrinsics + * + * @param C Pointer to output matrix. + * @param A Pointer to left input matrix. + * @param B Pointer to right input matrix. + */ void mat3_vfms(double *C, double const *A, double const *B) { // load left hand matrix @@ -408,7 +458,6 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) next=P->offsets[++idx][0]; } - for(int k=L->rowptr[i];krowptr[i+1];k++) { //scale row i of U @@ -484,8 +533,6 @@ void prec_mapply3c(prec_t *restrict P, double *x) __m256d vxj = _mm256_loadu_pd(xj); __m256d vz = (vxj - vA[0]) - (vA[1] + vA[2]); - //vz =_mm256_blend_pd(vxj,vz,0x7); // 4th element unchanged - //_mm256_storeu_pd(xj,vz); double z[4]; _mm256_store_pd(z,vz); for(int n=0;n<3;n++) xj[n]=z[n]; @@ -498,8 +545,6 @@ void prec_mapply3c(prec_t *restrict P, double *x) vA[2] = _mm256_cvtps_pd(_mm_loadu_ps(A+6))*vx[2]; //0b01010101 __m256d vz = vA[0] + vA[1] + vA[2]; - //vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged - //_mm256_storeu_pd(xi,vz); double z[4]; _mm256_store_pd(z,vz); for(int k=0;k<3;k++) xi[k]=z[k]; @@ -524,10 +569,6 @@ void prec_mapply3c(prec_t *restrict P, double *x) __m256d vz = (vxi - vA[0]) - (vA[1] + vA[2]); vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged _mm256_storeu_pd(xi,vz); - - //double z[4]; - //_mm256_store_pd(z,vz); - //for(int k=0;k<3;k++) xi[k]=z[k]; } } @@ -565,8 +606,6 @@ void prec_dapply3c(prec_t *restrict P, double *x) __m256d vxj = _mm256_loadu_pd(xj); __m256d vz = (vxj - vA[0]) - (vA[1] + vA[2]); - //vz =_mm256_blend_pd(vxj,vz,0x7); // 4th element unchanged - //_mm256_storeu_pd(xj,vz); double z[4]; _mm256_store_pd(z,vz); for(int kk=0;kk<3;kk++) xj[kk]=z[kk]; @@ -579,8 +618,6 @@ void prec_dapply3c(prec_t *restrict P, double *x) vA[2] = _mm256_loadu_pd(A+6)*vx[2]; //0b01010101 __m256d vz = vA[0] + vA[1] + vA[2]; - //vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged - //_mm256_storeu_pd(xi,vz); double z[4]; _mm256_store_pd(z,vz); for(int k=0;k<3;k++) xi[k]=z[k]; @@ -605,10 +642,6 @@ void prec_dapply3c(prec_t *restrict P, double *x) __m256d vz = (vxi - vA[0]) - (vA[1] + vA[2]); vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged _mm256_storeu_pd(xi,vz); - - //double z[4]; - //_mm256_store_pd(z,vz); - //for(int k=0;k<3;k++) xi[k]=z[k]; } } diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index b442c7f7af..7d47480192 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -6,14 +6,21 @@ extern "C" { #include "bsr.h" -// preconditioner struct +/*! + * @brief Preconditioner struct. + */ typedef struct prec_t { + // lower triangular factor bsr_matrix *L; + // diagonal factor bsr_matrix *D; + // upper triangular factor bsr_matrix *U; + // number of off-diagonal ilu0 updates int noffsets; + // triplets uniquely identifying off-diagonal ilu0 updates int(*offsets)[3]; } prec_t; diff --git a/opm/simulators/linalg/mixed/vec.h b/opm/simulators/linalg/mixed/vec.h index 0c5a4bdad3..6ea9e52563 100644 --- a/opm/simulators/linalg/mixed/vec.h +++ b/opm/simulators/linalg/mixed/vec.h @@ -5,11 +5,25 @@ extern "C" { #endif +/** + * @brief Copy vector. + * + * @param y Pointer to output vector. + * @param x Pointer to inout vector. + * @param n Length of vectors. + */ void vec_copy(double *y, double const * x, int n) { for(int i=0;i