From d36c6979e3bdafee11aa645371ebc74d7fea0197 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Mon, 1 Jun 2026 17:21:56 -0500 Subject: [PATCH 01/39] mixed: matrix now supports 4x4 blocks --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 32 +++++++- opm/simulators/linalg/mixed/bsr.c | 73 +++++++++++++++++++ opm/simulators/linalg/mixed/bsr.h | 2 + opm/simulators/linalg/setupPropertyTree.cpp | 3 +- 4 files changed, 106 insertions(+), 4 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 1364fb3a6fd..81bf7e30bbe 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -57,7 +57,15 @@ 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]); + if constexpr(b==1){printf("MixedMatrixWrapper::mv does not support block size == 1!\n");getchar();} + else if constexpr(b==2){printf("MixedMatrixWrapper::mv does not support block size == 2!\n");getchar();} + 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 + { + printf("MixedMatrixWrapper::mv only supports block sizes < 5!\n"); + getchar(); + } } template @@ -65,7 +73,15 @@ 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); + if constexpr(b==1){printf("MixedMatrixWrapper::umv does not support block size == 1!\n");getchar();} + else if constexpr(b==2){printf("MixedMatrixWrapper::umv does not support block size == 2!\n");getchar();} + 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 + { + printf("MixedMatrixWrapper::umv only supports block sizes < 4!\n"); + getchar(); + } } template @@ -73,7 +89,17 @@ 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); + //bsr_vmspumv3(M_, &x[0][0], &y[0][0], alpha); + if constexpr(b==1){printf("MixedMatrixWrapper::usmv does not support block size == 1!\n");getchar();} + else if constexpr(b==2){printf("MixedMatrixWrapper::usmv does not support block size == 2!\n");getchar();} + 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 if constexpr(b==4){} + else + { + printf("MixedMatrixWrapper::usmv only supports block sizes < 5!\n"); + getchar(); + } } template diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c index ba23e68e01f..65032790e6e 100644 --- a/opm/simulators/linalg/mixed/bsr.c +++ b/opm/simulators/linalg/mixed/bsr.c @@ -185,6 +185,79 @@ 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;i Date: Mon, 1 Jun 2026 21:30:10 -0500 Subject: [PATCH 02/39] mixed: matrix now supports 2x2 blocks --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 7 +-- opm/simulators/linalg/mixed/bsr.c | 63 +++++++++++++++++++ opm/simulators/linalg/mixed/bsr.h | 2 + 3 files changed, 68 insertions(+), 4 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 81bf7e30bbe..67432a56f14 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -58,7 +58,7 @@ mv(const Vector& x, Vector& y) const { // mixed-precision block spmv (y = M.x) if constexpr(b==1){printf("MixedMatrixWrapper::mv does not support block size == 1!\n");getchar();} - else if constexpr(b==2){printf("MixedMatrixWrapper::mv does not support block size == 2!\n");getchar();} + 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 @@ -74,7 +74,7 @@ umv(const Vector& x, Vector& y) const { // mixed-precision block spmv with update (y += M.x) if constexpr(b==1){printf("MixedMatrixWrapper::umv does not support block size == 1!\n");getchar();} - else if constexpr(b==2){printf("MixedMatrixWrapper::umv does not support block size == 2!\n");getchar();} + 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 @@ -91,10 +91,9 @@ 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); if constexpr(b==1){printf("MixedMatrixWrapper::usmv does not support block size == 1!\n");getchar();} - else if constexpr(b==2){printf("MixedMatrixWrapper::usmv does not support block size == 2!\n");getchar();} + 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 if constexpr(b==4){} else { printf("MixedMatrixWrapper::usmv only supports block sizes < 5!\n"); diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c index 65032790e6e..a458bf41683 100644 --- a/opm/simulators/linalg/mixed/bsr.c +++ b/opm/simulators/linalg/mixed/bsr.c @@ -259,6 +259,69 @@ void bsr_vmspumv4(bsr_matrix *A, const double *x, double *y, double alpha) } } +void bsr_vmspmv2(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=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;innz; diff --git a/opm/simulators/linalg/mixed/bsr.h b/opm/simulators/linalg/mixed/bsr.h index 0e3743d42cd..5384dba8c20 100644 --- a/opm/simulators/linalg/mixed/bsr.h +++ b/opm/simulators/linalg/mixed/bsr.h @@ -64,6 +64,7 @@ void bsr_init(bsr_matrix *A, int nrows, int nnz, int b); * @param x Pointer to input vector. * @param y Pointer to output vector. */ +void bsr_vmspumv2(bsr_matrix *A, const double *x, double *y, double alpha); void bsr_vmspumv3(bsr_matrix *A, const double *x, double *y, double alpha); void bsr_vmspumv4(bsr_matrix *A, const double *x, double *y, double alpha); @@ -79,6 +80,7 @@ void bsr_vmspumv4(bsr_matrix *A, const double *x, double *y, double alpha); * @param y Pointer to output vector. */ +void bsr_vmspmv2(bsr_matrix *A, const double *x, double *y); void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y); void bsr_vmspmv4(bsr_matrix *A, const double *x, double *y); From 01b50907b1c308845ced4431989b3867a6eef1d1 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 2 Jun 2026 10:48:12 -0500 Subject: [PATCH 03/39] mixed: matrix now has fallback for blocks larger than 4x4 --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 60 +++++++++++++++++-- 1 file changed, 55 insertions(+), 5 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 67432a56f14..d89a3aa2316 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -63,8 +63,33 @@ mv(const Vector& x, Vector& y) const else if constexpr(b==4) bsr_vmspmv4(M_, &x[0][0], &y[0][0]); else { - printf("MixedMatrixWrapper::mv only supports block sizes < 5!\n"); - getchar(); + 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:: 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); if constexpr(b==1){printf("MixedMatrixWrapper::usmv does not support block size == 1!\n");getchar();} 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 { - printf("MixedMatrixWrapper::usmv only supports block sizes < 5!\n"); - getchar(); + 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 Date: Tue, 2 Jun 2026 17:48:39 -0500 Subject: [PATCH 04/39] mixed: replace binary literals with hex counterparts --- opm/simulators/linalg/mixed/bsr.c | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c index a458bf41683..afadadf2af4 100644 --- a/opm/simulators/linalg/mixed/bsr.c +++ b/opm/simulators/linalg/mixed/bsr.c @@ -206,10 +206,10 @@ void bsr_vmspmv4(bsr_matrix *A, const double *x, double *y) int j = colidx[k]; __m256d vx = _mm256_loadu_pd(x+b*j); - vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 0))*_mm256_permute4x64_pd(vx,0b00000000); - vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 4))*_mm256_permute4x64_pd(vx,0b01010101); - vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 8))*_mm256_permute4x64_pd(vx,0b10101010); - vA[3] += _mm256_cvtps_pd(_mm_loadu_ps(AA+12))*_mm256_permute4x64_pd(vx,0b11111111); + vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 0))*_mm256_permute4x64_pd(vx,0x00); // 0b00000000 + vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 4))*_mm256_permute4x64_pd(vx,0x55); // 0b01010101 + vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 8))*_mm256_permute4x64_pd(vx,0xAA); // 0b10101010 + vA[3] += _mm256_cvtps_pd(_mm_loadu_ps(AA+12))*_mm256_permute4x64_pd(vx,0xFF); // 0b11111111 } // sum over columns @@ -244,10 +244,10 @@ void bsr_vmspumv4(bsr_matrix *A, const double *x, double *y, double alpha) int j = colidx[k]; __m256d vx = _mm256_loadu_pd(x+b*j); - vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 0))*_mm256_permute4x64_pd(vx,0b00000000); - vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 4))*_mm256_permute4x64_pd(vx,0b01010101); - vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 8))*_mm256_permute4x64_pd(vx,0b10101010); - vA[3] += _mm256_cvtps_pd(_mm_loadu_ps(AA+12))*_mm256_permute4x64_pd(vx,0b11111111); + vA[0] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 0))*_mm256_permute4x64_pd(vx,0x00); // 0b00000000 + vA[1] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 4))*_mm256_permute4x64_pd(vx,0x55); // 0b01010101 + vA[2] += _mm256_cvtps_pd(_mm_loadu_ps(AA+ 8))*_mm256_permute4x64_pd(vx,0xAA); // 0b10101010 + vA[3] += _mm256_cvtps_pd(_mm_loadu_ps(AA+12))*_mm256_permute4x64_pd(vx,0xFF); // 0b11111111 } // sum over columns @@ -277,7 +277,7 @@ void bsr_vmspmv2(bsr_matrix *A, const double *x, double *y) const float *AA=data+4*k; int j = colidx[k]; __m256d vx = _mm256_loadu_pd(x+b*j); - vA += _mm256_cvtps_pd(_mm_loadu_ps(AA))*_mm256_permute4x64_pd(vx,0b01010000); + vA += _mm256_cvtps_pd(_mm_loadu_ps(AA))*_mm256_permute4x64_pd(vx,0x50); // 0b01010000 } // sum over columns @@ -309,7 +309,7 @@ void bsr_vmspumv2(bsr_matrix *A, const double *x, double *y, double alpha) int j = colidx[k]; __m256d vx = _mm256_loadu_pd(x+b*j); - vA += _mm256_cvtps_pd(_mm_loadu_ps(AA))*_mm256_permute4x64_pd(vx,0b01010000); + vA += _mm256_cvtps_pd(_mm_loadu_ps(AA))*_mm256_permute4x64_pd(vx,0x50); // 0b01010000 } // sum over columns From 15795ba844d510cbddf89ec4faa8cd5718bf757e Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Wed, 3 Jun 2026 12:28:10 -0500 Subject: [PATCH 05/39] mixed: preconditioner checks block size --- .../linalg/mixed/PreconditionerWrapper.hpp | 31 +++++++++++++++++-- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 779bf9ff608..b34067ff53b 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -35,7 +35,7 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate 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; @@ -92,6 +92,7 @@ update () // transpose each dense block to make them column-major constexpr int b = block_size; constexpr int bb=b*b; + double B[bb]; for(int k=0;kdbl[bb*k + i] = B[i]; } - use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0 + //use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0 + + if constexpr(b==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} + else if constexpr(b==2){printf("MixedPreconditioner::update does not support block size == 2!\n");getchar();} + else if constexpr(b==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); + else if constexpr(b==4){printf("MixedPreconditioner::update does not support block size == 4!\n");getchar();} + else + { + printf("MixedPreconditioner::update only supports block sizes < 5!\n"); + getchar(); + } + prec_downcast(prec_); } @@ -108,7 +120,20 @@ void MixedPreconditioner:: apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) { x=y; - prec_mapply3c(prec_,&x[0][0]); + + //prec_mapply3c(prec_,&x[0][0]); + + int const b = block_size; + if constexpr(b==1){printf("MixedPreconditioner::apply does not support block size == 1!\n");getchar();} + else if constexpr(b==2){printf("MixedPreconditioner::apply does not support block size == 2!\n");getchar();} + else if constexpr(b==3) prec_mapply3c(prec_,&x[0][0]); + else if constexpr(b==4){printf("MixedPreconditioner::apply does not support block size == 4!\n");getchar();} + else + { + printf("MixedPreconditioner::apply only supports block sizes < 5!\n"); + getchar(); + } + } } From f3a46c14d00fb37107b15f1e56b5c00c10360c76 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 4 Jun 2026 10:29:09 -0500 Subject: [PATCH 06/39] mixed: inverse and matrix multiplications for 4x4 blocks --- .../linalg/mixed/PreconditionerWrapper.hpp | 1 + opm/simulators/linalg/mixed/matvec.h | 194 ++++++++++++++++++ opm/simulators/linalg/mixed/prec.c | 46 +++++ opm/simulators/linalg/mixed/prec.h | 2 + 4 files changed, 243 insertions(+) create mode 100644 opm/simulators/linalg/mixed/matvec.h diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index b34067ff53b..d3e5ed2261c 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -108,6 +108,7 @@ update () else if constexpr(b==4){printf("MixedPreconditioner::update does not support block size == 4!\n");getchar();} else { + prec_test(); printf("MixedPreconditioner::update only supports block sizes < 5!\n"); getchar(); } diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h new file mode 100644 index 00000000000..8b2299f0b27 --- /dev/null +++ b/opm/simulators/linalg/mixed/matvec.h @@ -0,0 +1,194 @@ +#pragma once + +#ifdef __cplusplus +extern "C" { +#endif + +void mat_show(double const *A, int n, char const *name) +{ + printf("%s = [\n",name); + for(int i = 0;i #include +#include "matvec.h" prec_t *prec_alloc() { @@ -205,6 +206,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. * @@ -669,3 +675,43 @@ void prec_info(prec_t *P) bsr_info(P->D); bsr_info(P->U); } + +void prec_test() +{ +#if 0 + // verify 2x2 inverse and matrix-matrix multiplications + double A[4] = {1,0.2,0.3,4}; + double B[4] = {1,0.2,0.3,4}; + double C[4] = {1,0.2,0.3,4}; + + mat2_inv(A,A); + mat_show(A,2,"A"); + + mat2_rmul(B,A); + mat_show(B,2,"B"); + + mat2_lmul(A,C); + mat_show(C,2,"C"); +#endif + + + // verify 4x4 inverse and matrix-matrix multiplications + double AA[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double BB[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double CC[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double II[16] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; + mat_show(AA,4,"AA"); + + mat4_inv(AA,AA); + mat_show(AA,4,"AA"); + + mat4_vfms(II,AA,BB); + mat_show(II,4,"II"); + + mat4_rmul(BB,AA); + mat_show(BB,4,"BB"); + + mat4_lmul(AA,CC); + mat_show(CC,4,"CC"); +} + diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index 7d474801925..cfcae6b79e5 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -107,6 +107,8 @@ void prec_downcast(prec_t *P); */ void prec_info(prec_t *P); + +void prec_test(); #ifdef __cplusplus } #endif From 594e9154f4c02242dacb47bd42093f97f480363d Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 4 Jun 2026 11:44:41 -0500 Subject: [PATCH 07/39] mixed: ilu0 and dilu now supports 4x4 blocks --- .../linalg/mixed/PreconditionerWrapper.hpp | 4 +- opm/simulators/linalg/mixed/matvec.h | 2 +- opm/simulators/linalg/mixed/prec.c | 228 ++++++++++++++++++ opm/simulators/linalg/mixed/prec.h | 3 + 4 files changed, 234 insertions(+), 3 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index d3e5ed2261c..7dd02f6984d 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -105,7 +105,7 @@ update () if constexpr(b==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} else if constexpr(b==2){printf("MixedPreconditioner::update does not support block size == 2!\n");getchar();} else if constexpr(b==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); - else if constexpr(b==4){printf("MixedPreconditioner::update does not support block size == 4!\n");getchar();} + else if constexpr(b==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); else { prec_test(); @@ -128,7 +128,7 @@ apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) if constexpr(b==1){printf("MixedPreconditioner::apply does not support block size == 1!\n");getchar();} else if constexpr(b==2){printf("MixedPreconditioner::apply does not support block size == 2!\n");getchar();} else if constexpr(b==3) prec_mapply3c(prec_,&x[0][0]); - else if constexpr(b==4){printf("MixedPreconditioner::apply does not support block size == 4!\n");getchar();} + else if constexpr(b==4) prec_mapply4c(prec_,&x[0][0]); else { printf("MixedPreconditioner::apply only supports block sizes < 5!\n"); diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h index 8b2299f0b27..da5c3515180 100644 --- a/opm/simulators/linalg/mixed/matvec.h +++ b/opm/simulators/linalg/mixed/matvec.h @@ -184,7 +184,7 @@ void mat4_vfms(double *C, double const *A, double const *B) vz = _mm256_loadu_pd(C+4*j) - vz; // Store result in column j of matrix C - _mm256_store_pd(C+4*j,vz); + _mm256_storeu_pd(C+4*j,vz); } } diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index bc48c7278d8..e9b08f7b0d4 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -391,6 +391,76 @@ void prec_dilu_factorize(prec_t *P, bsr_matrix *A) } } +void prec_dilu_factorize4(prec_t *P, bsr_matrix *A) +{ + int nrows = A->nrows; + int b = A->b; + int bb = b*b; + + 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++) + { + 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++; + } + } + } + // reset rowptr of L + for(int i=nrows;i>0;i--) L->rowptr[i]=L->rowptr[i-1]; + L->rowptr[0]=0; + + // Factorizing + double scale[16]; //hard-coded to 4x4 blocks + for(int i=0;inrows;i++) + { + mat4_inv(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); + + //update diagonal of U + int j=L->colidx[k]; + mat4_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); + + //scale row i of U + mat4_lmul(scale,U->dbl+k*bb); + + //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(); + } + } + } + } +} + void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) { int nrows = A->nrows; @@ -470,6 +540,84 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) } } +void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A) +{ + + int nrows = A->nrows; + int b = A->b; + int bb = b*b; + + 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++) + { + 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++; + } + } + } + // 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 = P->offsets[idx][0]; + double scale[16]; //hard-coded to 4x4 blocks + for(int i=0;inrows;i++) + { + mat4_inv(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); + + //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); + } + } +} #if 0 /** * @brief In-place matrix-vector multiplication for 3x3 matrices. @@ -589,6 +737,86 @@ 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 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++) + { + __m256d vA[4], vx[4]; + + double *xi = x+b*i; + __m256d vxi = _mm256_loadu_pd(xi); + + vx[0] = _mm256_permute4x64_pd(vxi,0x00); + vx[1] = _mm256_permute4x64_pd(vxi,0x55); // 0b01010101 + vx[2] = _mm256_permute4x64_pd(vxi,0xAA); // 0b10101010 + vx[3] = _mm256_permute4x64_pd(vxi,0xFF); // 0b10101010 + 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 die to structural symmetry? + 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]); + __m256d vz = vxj - (vA[0]+vA[1]) - (vA[2]+vA[3]); + + //double z[4]; + _mm256_storeu_pd(xj,vz); + //for(int n=0;n<3;n++) xj[n]=z[n]; + } + + // Muliply by (inverse) diagonal block + const float *A = D->flt+i*bb; + vA[0] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 0))*vx[0]; //0b01010101 + 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]; //0b01010101 + vA[3] = _mm256_cvtps_pd(_mm_loadu_ps(A+12))*vx[3]; //0b01010101 + __m256d vz = vA[0] + vA[1] + vA[2] + vA[3]; + + //double z[4]; + //_mm256_store_pd(z,vz); + //for(int k=0;k<3;k++) xi[k]=z[k]; + _mm256_storeu_pd(xi,vz); + } + + // Upper triangular solve assuming nonzeros stored in original order + 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); + 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); // 0b10101010 + } + double *xi = x+b*(i-1); + __m256d vxi = _mm256_loadu_pd(xi); + //__m256d vz = (vxi - vA[0]) - (vA[1] + vA[2]); + __m256d vz = vxi - (vA[0]+vA[1]) - (vA[2]+vA[3]); + //vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged + _mm256_storeu_pd(xi,vz); + } +} + void prec_dapply3c(prec_t *restrict P, double *x) { bsr_matrix *L = P->L; diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index cfcae6b79e5..84d11fe7b49 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -64,6 +64,7 @@ int prec_analyze(bsr_matrix *M, int (*offsets)[3]); * @param A Pointer to bsr matrix. */ void prec_dilu_factorize(prec_t *P, bsr_matrix *A); +void prec_dilu_factorize4(prec_t *P, bsr_matrix *A); /** * @brief ILU0 factorization. @@ -72,6 +73,7 @@ void prec_dilu_factorize(prec_t *P, bsr_matrix *A); * @param A Pointer to bsr matrix. */ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); +void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A); /** * @brief Preconditioner application in mixed-precision. @@ -82,6 +84,7 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); * @apram x Pointer to input/output vector */ void prec_mapply3c(prec_t *P, double *x); +void prec_mapply4c(prec_t *P, double *x); /** * @brief Preconditioner applicationin double-precision. From c9a3d2161b9c400f74b1565e2626b840b58e1b3f Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 4 Jun 2026 14:32:09 -0500 Subject: [PATCH 08/39] mixed: ilu0 and dilu now supports 2x2 blocks --- .../linalg/mixed/PreconditionerWrapper.hpp | 4 +- opm/simulators/linalg/mixed/matvec.h | 13 +- opm/simulators/linalg/mixed/prec.c | 220 +++++++++++++++++- opm/simulators/linalg/mixed/prec.h | 3 + 4 files changed, 233 insertions(+), 7 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 7dd02f6984d..b9adec3fd7a 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -103,7 +103,7 @@ update () //use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0 if constexpr(b==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} - else if constexpr(b==2){printf("MixedPreconditioner::update does not support block size == 2!\n");getchar();} + else if constexpr(b==2) use_dilu_ ? prec_dilu_factorize2(prec_, mixed_matrix_) : prec_ilu0_factorize2(prec_, mixed_matrix_); else if constexpr(b==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); else if constexpr(b==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); else @@ -126,7 +126,7 @@ apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) int const b = block_size; if constexpr(b==1){printf("MixedPreconditioner::apply does not support block size == 1!\n");getchar();} - else if constexpr(b==2){printf("MixedPreconditioner::apply does not support block size == 2!\n");getchar();} + 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 diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h index da5c3515180..c2fb78b21ef 100644 --- a/opm/simulators/linalg/mixed/matvec.h +++ b/opm/simulators/linalg/mixed/matvec.h @@ -14,7 +14,7 @@ void mat_show(double const *A, int n, char const *name) } printf("]\n\n"); } -#if 0 + void mat2_rmul(double *A, double const *B) { double M[4]; @@ -46,7 +46,16 @@ void mat2_inv(double *A, double const *B) double inv_det = 1.0/(M[0]*M[3]-M[1]*M[2]); for(int k=0;k<4;k++) A[k]=inv_det*M[k]; } -#endif + +void mat2_vfms(double *C, double const *A, double const *B) +{ + double M[4]; + M[0] = A[0]*B[0] + A[2]*B[1]; + M[1] = A[1]*B[0] + A[3]*B[1]; + M[2] = A[0]*B[2] + A[2]*B[3]; + M[3] = A[1]*B[2] + A[3]*B[3]; + for(int k=0;k<4;k++) C[k]-=M[k]; +} /** * @brief Matrix inverse for 4x4 matrix. diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index e9b08f7b0d4..65ba90bf571 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -195,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. * @@ -320,6 +325,76 @@ void mat3_vfms(double *C, double const *A, double const *B) } } +void prec_dilu_factorize2(prec_t *P, bsr_matrix *A) +{ + int nrows = A->nrows; + int b = A->b; + int bb = b*b; + + 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++) + { + int j=A->colidx[k]; + if(jrowptr[j]; + vec_copy4(L->dbl + bb*kL, A->dbl + bb*k); + L->rowptr[j]++; + } + else if(j==i) // struct-copy of D + { + vec_copy4(D->dbl + bb*i, A->dbl + bb*k); + } + else if(j>i) // struct-copy of U + { + vec_copy4(U->dbl + bb*kU, A->dbl + bb*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 + double scale[4]; //hard-coded to 2x2 blocks + for(int i=0;inrows;i++) + { + 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 + mat2_rmul(L->dbl+k*bb,scale); + + //update diagonal of U + int j=L->colidx[k]; + mat2_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); + + //scale row i of U + mat2_lmul(scale,U->dbl+k*bb); + + //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(); + } + } + } + } +} + void prec_dilu_factorize(prec_t *P, bsr_matrix *A) { @@ -461,6 +536,85 @@ void prec_dilu_factorize4(prec_t *P, bsr_matrix *A) } } +void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A) +{ + int nrows = A->nrows; + int b = A->b; + int bb = b*b; + + 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++) + { + int j=A->colidx[k]; + if(jrowptr[j]; + vec_copy4(L->dbl + bb*kL, A->dbl + bb*k); + L->rowptr[j]++; + } + else if(j==i) // struct-copy of D + { + vec_copy4(D->dbl + bb*i, A->dbl + bb*k); + } + else if(j>i) // struct-copy of U + { + vec_copy4(U->dbl + bb*kU, A->dbl + bb*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 = P->offsets[idx][0]; + double scale[4]; //hard-coded to 2x2 blocks + for(int i=0;inrows;i++) + { + 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 + mat2_rmul(L->dbl+k*bb,scale); + + //update diagonal D + int j=L->colidx[k]; + mat2_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 + 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) { int nrows = A->nrows; @@ -663,6 +817,63 @@ static inline void mat3_vecfms(double *y, const double *A, const double *x) for(int k=0;k<3;k++) y[k]-=z[k]; } #endif +void prec_mapply2c(prec_t *restrict P, double *x) +{ + bsr_matrix *L = P->L; + 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++) + { + __m256d vA, vx; + + double *xi = x+b*i; + __m256d vxi = _mm256_loadu_pd(xi); + + vx = _mm256_permute4x64_pd(vxi,0x50); // 0b01010000 + 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 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); + } +} void prec_mapply3c(prec_t *restrict P, double *x) { @@ -906,23 +1117,25 @@ void prec_info(prec_t *P) void prec_test() { -#if 0 // verify 2x2 inverse and matrix-matrix multiplications double A[4] = {1,0.2,0.3,4}; double B[4] = {1,0.2,0.3,4}; double C[4] = {1,0.2,0.3,4}; + double I[4] = {1,0,0,1}; mat2_inv(A,A); mat_show(A,2,"A"); + mat2_vfms(I,A,B); + mat_show(I,2,"I"); + mat2_rmul(B,A); mat_show(B,2,"B"); mat2_lmul(A,C); mat_show(C,2,"C"); -#endif - +#if 0 // verify 4x4 inverse and matrix-matrix multiplications double AA[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; double BB[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; @@ -941,5 +1154,6 @@ void prec_test() mat4_lmul(AA,CC); mat_show(CC,4,"CC"); +#endif } diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index 84d11fe7b49..cca6c114804 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -63,6 +63,7 @@ int prec_analyze(bsr_matrix *M, int (*offsets)[3]); * @param P Pointer preconditioner object. * @param A Pointer to bsr matrix. */ +void prec_dilu_factorize2(prec_t *P, bsr_matrix *A); void prec_dilu_factorize(prec_t *P, bsr_matrix *A); void prec_dilu_factorize4(prec_t *P, bsr_matrix *A); @@ -72,6 +73,7 @@ void prec_dilu_factorize4(prec_t *P, bsr_matrix *A); * @param P Pointer preconditioner object. * @param A Pointer to bsr matrix. */ +void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A); void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A); @@ -83,6 +85,7 @@ void prec_ilu0_factorize4(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); From 23847c71376e9fa15116eb03794035730799c2b9 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Mon, 8 Jun 2026 09:22:14 -0500 Subject: [PATCH 09/39] mixed: avx2 version of 4x4 inverse --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 4 +- opm/simulators/linalg/mixed/bsr.c | 5 +- opm/simulators/linalg/mixed/matvec.h | 281 ++++++++++++++++++ opm/simulators/linalg/mixed/prec.c | 71 +++-- 4 files changed, 329 insertions(+), 32 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index d89a3aa2316..3b8ab1ebe28 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -162,11 +162,11 @@ 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_); + //bsr_downcast(M_); } } // namespace Opm #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c index afadadf2af4..90b9fa8eb0b 100644 --- a/opm/simulators/linalg/mixed/bsr.c +++ b/opm/simulators/linalg/mixed/bsr.c @@ -47,8 +47,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 = aligned_alloc(64,b*b*nnz*sizeof(double)); + A->flt = aligned_alloc(64,b*b*nnz*sizeof(float)); assert(A->rowptr); assert(A->colidx); diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h index c2fb78b21ef..97edd0cd4bc 100644 --- a/opm/simulators/linalg/mixed/matvec.h +++ b/opm/simulators/linalg/mixed/matvec.h @@ -4,6 +4,17 @@ extern "C" { #endif +void mat_fshow(float const *A, int n, char const *name) +{ + printf("%s = [\n",name); + for(int i = 0;i [0 1 4 5] + row1 = _mm_loadh_pi(_mm_loadl_pi(row1, (__m64*)(A+8)), (__m64*)(A +12)); // [8 9 x x] -> [8 9 C D] + + row0 = _mm_shuffle_ps(tmp1, row1, 0x88); // [0 4 8 C] 0b10 00 10 00 + row1 = _mm_shuffle_ps(row1, tmp1, 0xDD); // [9 D 1 5] 0b11 01 11 01 +*/ + tmp1 = _mm_unpacklo_ps(col0,col1); // [0 4 1 5] + row1 = _mm_unpacklo_ps(col2,col3); // [8 C 9 D] + + row0 = _mm_shuffle_ps(tmp1, row1, 0x44); // [0 4 8 C] 0b 01 00 01 00 + row1 = _mm_shuffle_ps(row1, tmp1, 0xEE); // [9 D 1 5] 0b 11 10 11 10 // notice flipped order +/* + tmp1 = _mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(A+ 2)), (__m64*)(A+ 6)); // [2 3 x x] -> [2 3 6 7] + row3 = _mm_loadh_pi(_mm_loadl_pi(row3, (__m64*)(A+10)), (__m64*)(A+14)); // [A B x x] -> [A B E F] + + row2 = _mm_shuffle_ps(tmp1, row3, 0x88); // [2 6 A E] 0b 10 00 10 00 + row3 = _mm_shuffle_ps(row3, tmp1, 0xDD); // [B F 3 7] 0b 11 01 11 01 +*/ + tmp1 = _mm_unpackhi_ps(col0,col1); // [2 6 3 7] 0b 01 00 01 00 + row3 = _mm_unpackhi_ps(col2,col3); // [A E B F] 0b 11 10 11 10 + + row2 = _mm_shuffle_ps(tmp1, row3, 0x44); // [2 6 A E] 0b 01 00 01 00 + row3 = _mm_shuffle_ps(row3, tmp1, 0xEE); // [B F 3 7] 0b 11 10 11 10 // notice flipped order + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row2, row3); // [2B 6F A3 E7] + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); // 1011 0001 [2B 6F A3 E7] -> [6F 2B E7 A3] + + col0 = _mm_mul_ps(row1, tmp1); // [96F D2B 1E7 5A3] + col1 = _mm_mul_ps(row0, tmp1); // [06F 42B 8E7 CA3] + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); // 0100 1110 [2B 6F A3 E7] -> [A3 E7 2B 6F] + + col0 = _mm_sub_ps(_mm_mul_ps(row1, tmp1), col0);// [9A3 DE7 12B 56F] - [96F D2B 1E7 5A3] = [9(A3-6F) D(E7-2B) 1(2B-E7) 5(6F-A3)] + col1 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col1);// [0A3 4E7 82B C6F] - [06F 42B 8E7 CA3] = [0(A3-6F) 4(E7-2B) 8(2B-E7) C(6F-A3)] + col1 = _mm_shuffle_ps(col1, col1, 0x4E); // [8(2B-E7) C(6F-A3) 0(A3-6F) 4(E7-2B)] + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row1, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col0 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col0); + col3 = _mm_mul_ps(row0, tmp1); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col0 = _mm_sub_ps(col0, _mm_mul_ps(row3, tmp1)); + col3 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col3); + col3 = _mm_shuffle_ps(col3, col3, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + row2 = _mm_shuffle_ps(row2, row2, 0x4E); + + col0 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col0); + col2 = _mm_mul_ps(row0, tmp1); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col0 = _mm_sub_ps(col0, _mm_mul_ps(row2, tmp1)); + col2 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col2); + col2 = _mm_shuffle_ps(col2, col2, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col2); + col3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), col3); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), col2); + col3 = _mm_sub_ps(col3, _mm_mul_ps(row2, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col1 = _mm_sub_ps(col1, _mm_mul_ps(row2, tmp1)); + col2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col2); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col1); + col2 = _mm_sub_ps(col2, _mm_mul_ps(row1, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col1); + col3 = _mm_sub_ps(col3, _mm_mul_ps(row1, tmp1)); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col1 = _mm_sub_ps(col1, _mm_mul_ps(row3, tmp1)); + col3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col3); +// ----------------------------------------------- + det = _mm_mul_ps(row0, col0); + det = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det); + det = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det); + tmp1 = _mm_rcp_ss(det); + + det = _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1))); + det = _mm_shuffle_ps(det, det, 0x00); + +/* + _mm_store_ps(A+ 0, col0); + _mm_store_ps(A+ 4, col1); + _mm_store_ps(A+ 8, col2); + _mm_store_ps(A+12, col3); + return; +*/ +/* + col0 = _mm_mul_ps(det, col0); + _mm_storel_pi((__m64*)(A+0), col0); + _mm_storeh_pi((__m64*)(A+2), col0); + + col1 = _mm_mul_ps(det, col1); + _mm_storel_pi((__m64*)(A+4), col1); + _mm_storeh_pi((__m64*)(A+6), col1); + + col2 = _mm_mul_ps(det, col2); + _mm_storel_pi((__m64*)(A+ 8), col2); + _mm_storeh_pi((__m64*)(A+10), col2); + + col3 = _mm_mul_ps(det, col3); + _mm_storel_pi((__m64*)(A+12), col3); + _mm_storeh_pi((__m64*)(A+14), col3); +*/ + + _mm_store_ps(A+ 0,_mm_mul_ps(det, col0)); + _mm_store_ps(A+ 4,_mm_mul_ps(det, col1)); + _mm_store_ps(A+ 8,_mm_mul_ps(det, col2)); + _mm_store_ps(A+12,_mm_mul_ps(det, col3)); + +} + +// AVX2 double-precision translation of Intel's SSE single-precision 4x4 matrix inverse using Cramer's Rule +void mat4_vinv(double *invA, double const *A) +{ + __m256d col0, col1, col2, col3; + __m256d row0, row1, row2, row3; + __m256d tmp1, det; +/* + 0 4 8 C + 1 5 9 D + 2 6 A E + 3 7 B F +*/ + // ----------------------------------------------- + // extract rows from columns + // ----------------------------------------------- + col0 = _mm256_load_pd(A+ 0); + col1 = _mm256_load_pd(A+ 4); + col2 = _mm256_load_pd(A+ 8); + col3 = _mm256_load_pd(A+12); + + tmp1 = _mm256_unpacklo_pd(col0,col1); // [0 4 1 5] //all maps have to be updated due to per-lane unpack + row2 = _mm256_unpacklo_pd(col2,col3); // [8 C 9 D] + + row0 = _mm256_permute2f128_pd(tmp1, row2, 0x20); // [0 4 8 C] 0b 00 10 00 00 extract lower pairs + row2 = _mm256_permute2f128_pd(row2, tmp1, 0x13); // [9 D 1 5] 0b 00 01 00 11 extract and swap upper pairs + + tmp1 = _mm256_unpackhi_pd(col0,col1); // [2 6 3 7] 0b 01 00 01 00 + row1 = _mm256_unpackhi_pd(col2,col3); // [A E B F] 0b 11 10 11 10 + + row3 = _mm256_permute2f128_pd(tmp1, row1, 0x13); // [2 6 A E] 0b 01 00 01 00 + row1 = _mm256_permute2f128_pd(row1, tmp1, 0x20); // [B F 3 7] 0b 11 10 11 10 // notice flipped order + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(row2, row3); // [2B 6F A3 E7] + tmp1 = _mm256_permute_pd(tmp1, 0x05); // 1011 0001 [2B 6F A3 E7] -> [6F 2B E7 A3] reverse upper and lower pairs + + col0 = _mm256_mul_pd(row1, tmp1); // [96F D2B 1E7 5A3] + col1 = _mm256_mul_pd(row0, tmp1); // [06F 42B 8E7 CA3] + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col0 = _mm256_fmsub_pd(row1,tmp1,col0); // [9A3 DE7 12B 56F] - [96F D2B 1E7 5A3] = [9(A3-6F) D(E7-2B) 1(2B-E7) 5(6F-A3)] + col1 = _mm256_fmsub_pd(row0,tmp1,col1); // [0A3 4E7 82B C6F] - [06F 42B 8E7 CA3] = [0(A3-6F) 4(E7-2B) 8(2B-E7) C(6F-A3)] + col1 = _mm256_permute4x64_pd(col1,0x4E); + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(row1, row2); + tmp1 = _mm256_permute_pd(tmp1, 0x05); + + col0 = _mm256_fmadd_pd(row3,tmp1,col0); + col3 = _mm256_mul_pd(row0, tmp1); + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col0 = _mm256_fnmadd_pd(row3,tmp1,col0); + col3 = _mm256_fmsub_pd(row0,tmp1,col3); + col3 = _mm256_permute4x64_pd(col3,0x4E); + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(_mm256_permute4x64_pd(row1, 0x4E), row3); + tmp1 = _mm256_permute_pd(tmp1, 0x05); + row2 = _mm256_permute4x64_pd(row2, 0x4E); + + col0 = _mm256_fmadd_pd(row2,tmp1,col0); + col2 = _mm256_mul_pd(row0, tmp1); + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col0 = _mm256_fnmadd_pd(row2,tmp1,col0); + col2 = _mm256_fmsub_pd(row0,tmp1,col2); + col2 = _mm256_permute4x64_pd(col2,0x4E); + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(row0, row1); + tmp1 = _mm256_permute_pd(tmp1, 0x05); + + col2 = _mm256_fmadd_pd(row3,tmp1,col2); + col3 = _mm256_fmsub_pd(row2,tmp1,col3); + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col2 = _mm256_fmsub_pd(row3,tmp1,col2); + col3 = _mm256_fnmadd_pd(row2,tmp1,col3); + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(row0, row3); + tmp1 = _mm256_permute_pd(tmp1, 0x05); + + col1 = _mm256_fnmadd_pd(row2,tmp1,col1); + col2 = _mm256_fmadd_pd(row1,tmp1,col2); + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col1 = _mm256_fmadd_pd(row2,tmp1,col1); + col2 = _mm256_fnmadd_pd(row1,tmp1,col2); + // ----------------------------------------------- + tmp1 = _mm256_mul_pd(row0, row2); + tmp1 = _mm256_permute_pd(tmp1, 0x05); + + col1 = _mm256_fmadd_pd(row3,tmp1,col1); + col3 = _mm256_fnmadd_pd(row1,tmp1,col3); + + tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); + + col1 = _mm256_fnmadd_pd(row3,tmp1,col1); + col3 = _mm256_fmadd_pd(row1,tmp1,col3); +// ----------------------------------------------- + det = _mm256_mul_pd(row0, col0); + det = _mm256_add_pd(_mm256_permute4x64_pd(det, 0x4E), det); + det = _mm256_add_pd(_mm256_permute_pd(det, 0x05), det); + det = _mm256_permute4x64_pd(det,0x00); + det = _mm256_div_pd(_mm256_set1_pd(1.0),det); + + _mm256_store_pd(invA+ 0,_mm256_mul_pd(det, col0)); + _mm256_store_pd(invA+ 4,_mm256_mul_pd(det, col1)); + _mm256_store_pd(invA+ 8,_mm256_mul_pd(det, col2)); + _mm256_store_pd(invA+12,_mm256_mul_pd(det, col3)); + +} #ifdef __cplusplus } #endif diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index 65ba90bf571..1bcfabc472c 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -508,7 +508,7 @@ void prec_dilu_factorize4(prec_t *P, bsr_matrix *A) double scale[16]; //hard-coded to 4x4 blocks for(int i=0;inrows;i++) { - mat4_inv(scale,D->dbl+i*bb); + 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++) { @@ -698,8 +698,8 @@ void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A) { int nrows = A->nrows; - int b = A->b; - int bb = b*b; + int const b = 4; + int const bb =16; bsr_matrix *L=P->L; bsr_matrix *D=P->D; @@ -736,10 +736,10 @@ void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A) // Factorizing int idx=0; int next = P->offsets[idx][0]; - double scale[16]; //hard-coded to 4x4 blocks + double scale[16] __attribute__((aligned(64))); //hard-coded to 4x4 blocks for(int i=0;inrows;i++) { - mat4_inv(scale,D->dbl+i*bb); + 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++) { @@ -954,8 +954,8 @@ void prec_mapply4c(prec_t *restrict P, double *x) bsr_matrix *D = P->D; bsr_matrix *U = P->U; - int b=L->b; - int bb=b*b; + int const b=4; + int const bb=16; __m256d mm256_zero_pd =_mm256_setzero_pd(); @@ -967,14 +967,14 @@ void prec_mapply4c(prec_t *restrict P, double *x) double *xi = x+b*i; __m256d vxi = _mm256_loadu_pd(xi); - vx[0] = _mm256_permute4x64_pd(vxi,0x00); + 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); // 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 die to structural symmetry? + 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]; @@ -982,29 +982,22 @@ void prec_mapply4c(prec_t *restrict P, double *x) double *xj = x+b*j; __m256d vxj = _mm256_loadu_pd(xj); - //__m256d vz = (vxj - vA[0]) - (vA[1] + vA[2]); __m256d vz = vxj - (vA[0]+vA[1]) - (vA[2]+vA[3]); - - //double z[4]; _mm256_storeu_pd(xj,vz); - //for(int n=0;n<3;n++) xj[n]=z[n]; } // Muliply by (inverse) diagonal block const float *A = D->flt+i*bb; - vA[0] = _mm256_cvtps_pd(_mm_loadu_ps(A+ 0))*vx[0]; //0b01010101 + 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]; //0b01010101 - vA[3] = _mm256_cvtps_pd(_mm_loadu_ps(A+12))*vx[3]; //0b01010101 - __m256d vz = vA[0] + vA[1] + vA[2] + vA[3]; + 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 - //double z[4]; - //_mm256_store_pd(z,vz); - //for(int k=0;k<3;k++) xi[k]=z[k]; + __m256d vz = vA[0] + vA[1] + vA[2] + vA[3]; _mm256_storeu_pd(xi,vz); } - // Upper triangular solve assuming nonzeros stored in original order + // Upper triangular solve assuming ones on diagonal` for(int i=U->ncols;i>0;i--) { __m256d vA[4]; @@ -1014,16 +1007,15 @@ void prec_mapply4c(prec_t *restrict P, double *x) 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); + 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); // 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]); __m256d vz = vxi - (vA[0]+vA[1]) - (vA[2]+vA[3]); - //vz =_mm256_blend_pd(vxi,vz,0x7); // 4th element unchanged _mm256_storeu_pd(xi,vz); } } @@ -1117,6 +1109,7 @@ void prec_info(prec_t *P) void prec_test() { +#if 0 // verify 2x2 inverse and matrix-matrix multiplications double A[4] = {1,0.2,0.3,4}; double B[4] = {1,0.2,0.3,4}; @@ -1134,7 +1127,7 @@ void prec_test() mat2_lmul(A,C); mat_show(C,2,"C"); - +#endif #if 0 // verify 4x4 inverse and matrix-matrix multiplications double AA[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; @@ -1143,7 +1136,7 @@ void prec_test() double II[16] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; mat_show(AA,4,"AA"); - mat4_inv(AA,AA); + mat4_inv2(AA,AA); mat_show(AA,4,"AA"); mat4_vfms(II,AA,BB); @@ -1155,5 +1148,27 @@ void prec_test() mat4_lmul(AA,CC); mat_show(CC,4,"CC"); #endif + // verify 4x4 inverse and matrix-matrix multiplications + double A[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double B[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double C[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; + double I[16] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; +/* + mat_fshow(A,4,"A"); + mat4_inv2(A); + mat_fshow(A,4,"A"); +*/ + mat_show(A,4,"A"); + mat4_vinv(A,A); + mat_show(A,4,"A"); + + mat4_vfms(I,A,B); + mat_show(I,4,"I"); + + mat4_rmul(B,A); + mat_show(B,4,"B"); + + mat4_lmul(A,C); + mat_show(C,4,"C"); } From 2ab0e13a30d0e7a532da374c6134d3c2f44c80e8 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Mon, 8 Jun 2026 10:32:13 -0500 Subject: [PATCH 10/39] mixed: legacy implementation now supports 4x4 blocks --- opm/simulators/linalg/mixed/bslv.c | 71 +++++++++++++++++++++ opm/simulators/linalg/mixed/bslv.h | 1 + opm/simulators/linalg/mixed/prec.c | 2 +- opm/simulators/linalg/mixed/wrapper.hpp | 21 +++--- opm/simulators/linalg/setupPropertyTree.cpp | 4 +- 5 files changed, 88 insertions(+), 11 deletions(-) diff --git a/opm/simulators/linalg/mixed/bslv.c b/opm/simulators/linalg/mixed/bslv.c index 77ba71b6b52..2fa2fef4254 100644 --- a/opm/simulators/linalg/mixed/bslv.c +++ b/opm/simulators/linalg/mixed/bslv.c @@ -172,6 +172,77 @@ int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x return j == max_iter ? j : ++j; } +int bslv_pbicgstab4m(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; + + double * restrict e = mem->e; + const double *r0 = b; + //const double * restrict r0 = mem->dtmp[0]; //access randomly initialized one-dimensional shadow space + double * restrict p_j = mem->dtmp[1]; + double * restrict q_j = mem->dtmp[2]; + double * restrict r_j = mem->dtmp[3]; + double * restrict s_j = mem->dtmp[4]; + double * restrict t_j = mem->dtmp[5]; + double * restrict v_j = mem->dtmp[6]; + double * restrict x_j = x; + + prec_t * restrict P = mem->P; + mem->use_dilu ? prec_dilu_factorize4(P,A) : prec_ilu0_factorize4(P,A); // choose dilu or ilu0 + prec_downcast(P); + + vec_fill(x_j,0.0,n); + vec_copy(r_j,b,n); + 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; - int const b = 4; + //int const b = 4; int const bb =16; bsr_matrix *L=P->L; diff --git a/opm/simulators/linalg/mixed/wrapper.hpp b/opm/simulators/linalg/mixed/wrapper.hpp index cd4bc2253e8..c8e0ee37077 100644 --- a/opm/simulators/linalg/mixed/wrapper.hpp +++ b/opm/simulators/linalg/mixed/wrapper.hpp @@ -14,6 +14,9 @@ class MixedSolver : public InverseOperator { public: + // extract block size + static constexpr auto block_size = X::block_type::dimension; + MixedSolver(const M &A, double tol, int maxiter, bool use_dilu) { // verify that well contributions are added to the matrix @@ -25,8 +28,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<3 || 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 +64,27 @@ 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 = 0; + if constexpr(N==3) count = bslv_pbicgstab3m(mem_, jacobian_, &b[0][0], &x[0][0]); + else if constexpr(N==4) count = bslv_pbicgstab4m(mem_, jacobian_, &b[0][0], &x[0][0]); // return convergence information res.converged = (mem_->e[count] < mem_->tol); diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 03cf4d0f28b..76ea086a6d3 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -367,8 +367,8 @@ getSolverString(const FlowLinearSolverParameters& p) } else { - //return {"bicgstab"}; - return {"mixed-precision"}; + return {"bicgstab"}; + //return {"mixed-precision"}; } } From aa9564e600ae1931d125a658c7c4132c2d4978d2 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 9 Jun 2026 10:51:04 -0500 Subject: [PATCH 11/39] mixed: ilu0/dilu apply supports blocks larger than 4x4 --- .../linalg/mixed/PreconditionerWrapper.hpp | 77 ++++++++++++++++++- 1 file changed, 75 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index b9adec3fd7a..37c8175e2d2 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -83,6 +83,10 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate bsr_matrix *mixed_matrix_; prec_t *prec_; int nnz_; + + + void matvec_mul(double *y, float const *A, double const * x); + void matvec_mulsub(double *y, float const *A, double const * x); }; template @@ -122,19 +126,88 @@ 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){printf("MixedPreconditioner::apply does not support block size == 1!\n");getchar();} 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); + } + } + + } +/* else { printf("MixedPreconditioner::apply only supports block sizes < 5!\n"); getchar(); } +*/ +} +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 Date: Tue, 9 Jun 2026 13:22:12 -0500 Subject: [PATCH 12/39] mixed: ilu0/dilu update supports blocks larger than 4x4 --- .../linalg/mixed/PreconditionerWrapper.hpp | 200 ++++++++++++++++-- 1 file changed, 179 insertions(+), 21 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 37c8175e2d2..35621624d8d 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -87,6 +87,11 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate void matvec_mul(double *y, float const *A, double const * x); void matvec_mulsub(double *y, float const *A, double const * x); + void mat_copy(double *C, double const * A); + void mat_inv(double *invA, const double *A); + void mat_mulsub(double *C, double const *A, double const * B); + void mat_rmul(double *C, double const *A); + void mat_lmul(double const *A, double *C); }; template @@ -94,27 +99,97 @@ void MixedPreconditioner:: update () { // transpose each dense block to make them column-major - constexpr int b = block_size; - constexpr int bb=b*b; + constexpr int N = block_size; + constexpr int NN=N*N; - double B[bb]; + double B[NN]; for(int k=0;kdbl[bb*k + i] = B[i]; + for(int i=0;idbl[NN*k + i] = B[i]; } - //use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); // choose dilu or ilu0 - - if constexpr(b==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} - else if constexpr(b==2) use_dilu_ ? prec_dilu_factorize2(prec_, mixed_matrix_) : prec_ilu0_factorize2(prec_, mixed_matrix_); - else if constexpr(b==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); - else if constexpr(b==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); + if constexpr(N==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} + else if constexpr(N==2) use_dilu_ ? prec_dilu_factorize2(prec_, mixed_matrix_) : prec_ilu0_factorize2(prec_, mixed_matrix_); + else if constexpr(N==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); + else if constexpr(N==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); else { - prec_test(); - printf("MixedPreconditioner::update only supports block sizes < 5!\n"); - getchar(); + 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(); } prec_downcast(prec_); @@ -171,13 +246,6 @@ apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) } } -/* - else - { - printf("MixedPreconditioner::apply only supports block sizes < 5!\n"); - getchar(); - } -*/ } template @@ -210,5 +278,95 @@ matvec_mulsub(double *y, float const *A, double const * x) 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 Date: Wed, 17 Jun 2026 10:44:13 -0500 Subject: [PATCH 13/39] mixed: support WellModelMatrixAdapter --- opm/simulators/linalg/WellOperators.hpp | 2 + .../linalg/mixed/PreconditionerWrapper.hpp | 2 +- opm/simulators/linalg/mixed/SolverAdapter.hpp | 55 ++++++++++++++++--- opm/simulators/linalg/mixed/matvec.h | 21 ++++++- opm/simulators/linalg/mixed/prec.c | 8 ++- 5 files changed, 75 insertions(+), 13 deletions(-) diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index f11c080b72f..95510203565 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 diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 35621624d8d..73d89baf311 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -189,7 +189,7 @@ update () mat_lmul(scale,U->dbl+k*NN); } } - //prec_test(); + //prec_test(); getchar(); } prec_downcast(prec_); diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 295889b9393..3fec74b0c27 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -63,7 +63,7 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct return std::sqrt(this->dot(x, x)); } }; - +/* //! @brief Generalized mixed precision operator interface //! //! @tparam Matrix the block-matrix used by linear operator @@ -83,8 +83,10 @@ template struct MixedOperator { using type = Dune::MatrixAdapter; + //using type = Dune::AssembledLinearOperator; + //using type = Opm::WellModelMatrixAdapter; }; - +*/ //! @brief Wraps mixed precision //! @@ -100,8 +102,8 @@ class MixedBiCGSTABSolver:public InverseOperator using AbstractScalarProductType = Dune::ScalarProduct; static constexpr auto block_size = Vector::block_type::dimension; - using MixedMatrixType = Opm::MixedMatrixWrapper; - using MixedOperatorType = MixedOperator::type; + using MixedMatrixType = Opm::MixedMatrixWrapper; + //using MixedOperatorType = MixedOperator::type; using OptimizedProductType = SeqOptmizedProduct; //! @brief constructor @@ -164,17 +166,50 @@ class MixedBiCGSTABSolver:public InverseOperator //initialize mixed operator and optimized scalar product double_operator_ = op; - if constexpr (std::is_same_v) + using MatrixType = std::remove_const_tgetmat())>>; + if constexpr (std::is_same_v, Dune::MatrixAdapter>) { + //OPM_THROW(std::invalid_argument, "Dune::MatrixAdapter\n"); + using MixedOperatorType = Dune::MatrixAdapter; + mixed_operator_ = std::make_shared(*mixed_matrix_); + scalar_product_ = std::make_shared(); + } + else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) + { + //OPM_THROW(std::invalid_argument, "Opm::GhostLastMatrixAdapter\n"); + using MixedOperatorType = Dune::OverlappingSchwarzOperator; + mixed_operator_ = std::make_shared(*mixed_matrix_,comm); + scalar_product_ = sp; + } + else if constexpr (std::is_same_v, Opm::WellModelMatrixAdapter>) + { + //OPM_THROW(std::invalid_argument, "Opm::WellModelMatrixAdapter\n"); + using MixedOperatorType = Opm::WellModelMatrixAdapter; + using WellOperatorType = Opm::LinearOperatorExtra; + const WellOperatorType &wellOper = op->getwellOper(); + mixed_operator_ = std::make_shared(*mixed_matrix_, wellOper); + scalar_product_ = std::make_shared(); + //scalar_product_ = sp; + } + else if constexpr (std::is_same_v, Opm::WellModelGhostLastMatrixAdapter>) + { + OPM_THROW(std::invalid_argument, "Opm::WellModelGhostLastMatrixAdapter\n"); + } +/* + else if constexpr (std::is_same_v) + { + //mixed_operator_ = std::make_shared(*mixed_matrix_,op->wellOper_); + using MixedOperatorType = Dune::MatrixAdapter; mixed_operator_ = std::make_shared(*mixed_matrix_); scalar_product_ = std::make_shared(); } else { + using MixedOperatorType = Dune::OverlappingSchwarzOperator; mixed_operator_ = std::make_shared(*mixed_matrix_,comm); scalar_product_ = sp; } - +*/ //initialize bicgstab solver from Dune solver_ = std::make_shared>( *mixed_operator_, @@ -183,6 +218,7 @@ class MixedBiCGSTABSolver:public InverseOperator tol, // desired residual reduction factor maxiter, // maximum number of iterations verbosity); + } void apply(Vector &x, Vector &b, InverseOperatorResult &res) override @@ -205,12 +241,13 @@ class MixedBiCGSTABSolver:public InverseOperator Dune::SolverCategory::Category category() const override{return Dune::SolverCategory::overlapping;}; private: - using AbstractSolverType = Dune::InverseOperator; - + 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_operator_; std::shared_ptr mixed_matrix_; std::shared_ptr scalar_product_; double const *double_data_; diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h index 97edd0cd4bc..b1d07a14454 100644 --- a/opm/simulators/linalg/mixed/matvec.h +++ b/opm/simulators/linalg/mixed/matvec.h @@ -36,6 +36,25 @@ void mat2_rmul(double *A, double const *B) for(int k=0;k<4;k++) A[k]=M[k]; } + +/* + 0 2 | A C 0A 2B + 1 3 | B D 1A 3B 0A 1A 2B 3B 00 00 01 01 10 10 11 11 + +*/ +/* +void mat2_rmul(double *A, double const *B) +{ + // load matrices + __m256d vA, vB, vC; + vA = _mm256_loadu_pd(A); + vB = _mm256_loadu_pd(B); + vC = vA*_mm256_permute4x64_pd(vB,0x50);// + vA*_mm256_permute4x64_pd(vB,0x50); // 0b01010000 + //vC = _mm256_permute4x64_pd(vB,0xFA);// + vA*_mm256_permute4x64_pd(vB,0x50); // 0b01010000 + __m128d va = _mm256_extractf128_pd(vC,0) + _mm256_extractf128_pd(vC,1); + _mm_storeu_pd(A,va); +} +*/ void mat2_lmul(double const *A, double *B) { double M[4]; @@ -81,7 +100,7 @@ void mat4_inv(double *invA, const double *A) for(int k=0;k<4;k++) { - double scale=-1.0/M[5*k]; + double scale=-1.0/M[5*k]; for(int i=0;i<4;i++) M[i+4*k] *= i==k?0:scale; // scale column k for(int j=0;j<4;j++) { diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index 29073d866ec..d2ad01ce650 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -1109,13 +1109,15 @@ void prec_info(prec_t *P) void prec_test() { -#if 0 +//#if 0 // verify 2x2 inverse and matrix-matrix multiplications double A[4] = {1,0.2,0.3,4}; double B[4] = {1,0.2,0.3,4}; double C[4] = {1,0.2,0.3,4}; double I[4] = {1,0,0,1}; + mat_show(A,2,"A"); + mat2_inv(A,A); mat_show(A,2,"A"); @@ -1127,7 +1129,7 @@ void prec_test() mat2_lmul(A,C); mat_show(C,2,"C"); -#endif +//#endif #if 0 // verify 4x4 inverse and matrix-matrix multiplications double AA[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; @@ -1148,6 +1150,7 @@ void prec_test() mat4_lmul(AA,CC); mat_show(CC,4,"CC"); #endif +#if 0 // verify 4x4 inverse and matrix-matrix multiplications double A[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; double B[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; @@ -1170,5 +1173,6 @@ void prec_test() mat4_lmul(A,C); mat_show(C,4,"C"); +#endif } From e3367e286216cd9490bf551f1f5494c355a06735 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Mon, 22 Jun 2026 10:03:11 -0500 Subject: [PATCH 14/39] mixed: custom MixedGhostLastMatrixAdapter --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 37 ++++ opm/simulators/linalg/mixed/SolverAdapter.hpp | 206 ++++++++++++++++-- 2 files changed, 223 insertions(+), 20 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 3b8ab1ebe28..24d28db0ffa 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -168,5 +168,42 @@ update(double const *data) // downcast to single precision //bsr_downcast(M_); } + + +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 + virtual void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + comm_.project(y); + } + + // y += \alpha * A * x + virtual void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + comm_.project(y); + } + + virtual const M& getmat() const override { return A_; } + + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::overlapping; + } + +private: + const M& A_; + const C& comm_; +}; + + } // namespace Opm #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 3fec74b0c27..dfc0950ec35 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -63,6 +63,113 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct return std::sqrt(this->dot(x, x)); } }; + + +template +class GhostLastScalarProduct : public ScalarProduct +{ + public: + + static constexpr auto block_size = Vector::block_type::dimension; + + //! \brief The type of the vector to compute the scalar product on. + //! + //! E.g. BlockVector or another type fulfilling the ISTL + //! vector interface. + + /*! + * \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(); + } + + /*! + * \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. + It is assumed that the vectors are consistent on the interior+border + partition. + */ + void setLocalCount(int count){ count_ = count; } + + virtual 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_; + + //double result(0); + + // 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;icommunicator(); + double result = cc.sum(agg[0]); + return result; + } + + /*! \brief Norm of a right-hand side vector. + The vector must be consistent on the interior+border partition + */ + virtual double norm (const Vector& x) const override + { + return sqrt(dot(x,x)); + } + + //! Category of the scalar product (see SolverCategory::Category) + virtual SolverCategory::Category category() const override + { + return _category; + } + + private: + std::shared_ptr _communication; + SolverCategory::Category _category; + int count_; + + //int getLocalCount(const Comm& comm) const + 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 Generalized mixed precision operator interface //! @@ -104,7 +211,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 //! @@ -123,11 +229,53 @@ class MixedBiCGSTABSolver:public InverseOperator const int& verbosity, const Comm &comm) { +#if 1 + int halo; + int nrows; + int nnz=0; + + auto &A = op->getmat(); + if constexpr (std::is_same_v) + { + halo = 0; + nrows = A.N(); + nnz = A.nonzeroes(); + } + 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!=A.end(); row++) + for(auto row=A.begin(); row.index() < nrows; row++) + { + if(local_[irow++]==1) for(auto col = row->begin(); col != row->end(); col++) nnz++; + //else break; // This line is used to verify that all local nodes have indices lower than all ghost nodes. This appears to be true! + //irow++; + } + } + + printf("nnz = %d\n",nnz); + printf("A.nonzeroes() = %ld\n",A.nonzeroes()); + + printf("local = %d\n",nrows); + printf("halo = %d\n",halo); + printf("total = %ld\n",A.N()); + //getchar(); +#endif // Access matrix data from double precision operator +#if 0 auto &A = op->getmat(); int nrows = A.N(); int nnz = A.nonzeroes(); +#endif double_data_ = &A[0][0][0][0]; //allocate mixed matrix @@ -140,7 +288,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) { @@ -172,14 +320,19 @@ class MixedBiCGSTABSolver:public InverseOperator //OPM_THROW(std::invalid_argument, "Dune::MatrixAdapter\n"); using MixedOperatorType = Dune::MatrixAdapter; mixed_operator_ = std::make_shared(*mixed_matrix_); + using OptimizedProductType = SeqOptmizedProduct; scalar_product_ = std::make_shared(); } else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) { //OPM_THROW(std::invalid_argument, "Opm::GhostLastMatrixAdapter\n"); - using MixedOperatorType = Dune::OverlappingSchwarzOperator; + //using MixedOperatorType = Dune::OverlappingSchwarzOperator; + using MixedOperatorType = Opm::MixedGhostLastMatrixAdapter; mixed_operator_ = std::make_shared(*mixed_matrix_,comm); - scalar_product_ = sp; + + using OptimizedScalarProductType = GhostLastScalarProduct; + scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); + //scalar_product_ = sp; } else if constexpr (std::is_same_v, Opm::WellModelMatrixAdapter>) { @@ -188,6 +341,7 @@ class MixedBiCGSTABSolver:public InverseOperator using WellOperatorType = Opm::LinearOperatorExtra; const WellOperatorType &wellOper = op->getwellOper(); mixed_operator_ = std::make_shared(*mixed_matrix_, wellOper); + using OptimizedProductType = SeqOptmizedProduct; scalar_product_ = std::make_shared(); //scalar_product_ = sp; } @@ -195,21 +349,7 @@ class MixedBiCGSTABSolver:public InverseOperator { OPM_THROW(std::invalid_argument, "Opm::WellModelGhostLastMatrixAdapter\n"); } -/* - else if constexpr (std::is_same_v) - { - //mixed_operator_ = std::make_shared(*mixed_matrix_,op->wellOper_); - using MixedOperatorType = Dune::MatrixAdapter; - mixed_operator_ = std::make_shared(*mixed_matrix_); - scalar_product_ = std::make_shared(); - } - else - { - using MixedOperatorType = Dune::OverlappingSchwarzOperator; - mixed_operator_ = std::make_shared(*mixed_matrix_,comm); - scalar_product_ = sp; - } -*/ + //initialize bicgstab solver from Dune solver_ = std::make_shared>( *mixed_operator_, @@ -241,17 +381,43 @@ class MixedBiCGSTABSolver:public InverseOperator Dune::SolverCategory::Category category() const override{return Dune::SolverCategory::overlapping;}; private: + + 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 +/* + // This code snippet is used to check wheter or not all ghost indices occur last in the index set + // In general, that is NOT the case + if (idx->local().attribute()==1) count++; // count local indices + else + { + printf("debug: %d: %ld %d\n",count, idx->local().local(), idx->local().attribute()); + break; + } +*/ + 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_; }; } From 2f2432a01f440173529cf6dae19ac3c073528210 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 23 Jun 2026 09:54:52 -0500 Subject: [PATCH 15/39] mixed: verify local cell count in GhostLastScalarProduct --- opm/simulators/linalg/mixed/SolverAdapter.hpp | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index dfc0950ec35..bdadc38d6f7 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -86,6 +86,8 @@ class GhostLastScalarProduct : public ScalarProduct : _communication(com), _category(cat) { count_ = getLocalCount(); + int verify = verifyLocalCount(); + if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n"); } /*! @@ -166,6 +168,26 @@ class GhostLastScalarProduct : public ScalarProduct return count; } + 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 + } + }; From 89761b6b3eac82c4922af986deac06bcb734ee7b Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Wed, 8 Jul 2026 15:42:13 -0500 Subject: [PATCH 16/39] mixed: custom WellModelMixedGhostLastMatrixAdapter --- opm/simulators/linalg/WellOperators.hpp | 2 + opm/simulators/linalg/mixed/MatrixWrapper.hpp | 77 +++++++++++++++++++ opm/simulators/linalg/mixed/SolverAdapter.hpp | 16 +++- 3 files changed, 94 insertions(+), 1 deletion(-) diff --git a/opm/simulators/linalg/WellOperators.hpp b/opm/simulators/linalg/WellOperators.hpp index 95510203565..c145e44f0bb 100644 --- a/opm/simulators/linalg/WellOperators.hpp +++ b/opm/simulators/linalg/WellOperators.hpp @@ -359,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 24d28db0ffa..ac69e33bc5a 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -205,5 +205,82 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator }; +/*! + \brief Adapter to combine a matrix and another linear operator into + a combined linear operator. + + This is similar to WellModelMatrixAdapter, with the difference that + here we assume a parallel ordering of rows, where ghost rows are + located after interior rows. + */ +template +class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + using field_type = typename V::field_type; + using PressureMatrix = Dune::BCRSMatrix>; + + //! 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 + virtual void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + wellOper_.apply(x, y); + comm_.project(y); + } + + // y += \alpha * A * x + virtual void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + wellOper_.applyscaleadd(alpha, x, y); + comm_.project(y); + } + + 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(); + } + + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::overlapping; + } + + +protected: + + const M& A_ ; + const C& comm_ ; + const LinearOperatorExtra& wellOper_; +}; + + + + + } // namespace Opm #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index bdadc38d6f7..2b4a0233a37 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -369,7 +369,21 @@ class MixedBiCGSTABSolver:public InverseOperator } else if constexpr (std::is_same_v, Opm::WellModelGhostLastMatrixAdapter>) { - OPM_THROW(std::invalid_argument, "Opm::WellModelGhostLastMatrixAdapter\n"); + //OPM_THROW(std::invalid_argument, "Opm::WellModelGhostLastMatrixAdapter\n"); + 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 OptimizedScalarProductType = GhostLastScalarProduct; + scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); + } } //initialize bicgstab solver from Dune From 8cc0aaea1be2482525e7800888f73b5866d2a239 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 9 Jul 2026 13:12:44 -0500 Subject: [PATCH 17/39] mixed: post-rebase fixes --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 5 +++-- opm/simulators/linalg/mixed/SolverAdapter.hpp | 16 +--------------- 2 files changed, 4 insertions(+), 17 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index ac69e33bc5a..cea2753b2ab 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -32,8 +32,6 @@ class MixedMatrixWrapper //! @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); @@ -57,6 +55,7 @@ void MixedMatrixWrapper:: mv(const Vector& x, Vector& y) const { // mixed-precision block spmv (y = M.x) + int const b = block_size; if constexpr(b==1){printf("MixedMatrixWrapper::mv does not support block size == 1!\n");getchar();} 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]); @@ -98,6 +97,7 @@ void MixedMatrixWrapper:: umv(const Vector& x, Vector& y) const { // mixed-precision block spmv with update (y += M.x) + int const b = block_size; if constexpr(b==1){printf("MixedMatrixWrapper::umv does not support block size == 1!\n");getchar();} 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); @@ -114,6 +114,7 @@ void MixedMatrixWrapper:: usmv(double alpha, const Vector& x, Vector& y) const { // scaled mixed-precision block spmv with update (y += alpha * M.x) + int const b = block_size; if constexpr(b==1){printf("MixedMatrixWrapper::usmv does not support block size == 1!\n");getchar();} 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); diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 2b4a0233a37..941d07a4784 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -231,7 +231,7 @@ class MixedBiCGSTABSolver:public InverseOperator using AbstractScalarProductType = Dune::ScalarProduct; static constexpr auto block_size = Vector::block_type::dimension; - using MixedMatrixType = Opm::MixedMatrixWrapper; + using MixedMatrixType = Opm::MixedMatrixWrapper; //using MixedOperatorType = MixedOperator::type; //! @brief constructor @@ -320,20 +320,6 @@ 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. - using MatrixType = std::remove_const_tgetmat())>>; - if constexpr (std::is_same_v, Dune::MatrixAdapter>) - { - } - else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) - { - } - else - { - OPM_THROW(std::invalid_argument, "MixedBiCGSTABSolver only supports Dune::MatrixAdapter and Opm::GhostLastMatrixAdapter\n"); - } - //initialize mixed operator and optimized scalar product double_operator_ = op; using MatrixType = std::remove_const_tgetmat())>>; From a0bdd0351b75f1b33918e2a935393690825d2ec5 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Wed, 15 Jul 2026 13:49:30 -0500 Subject: [PATCH 18/39] mixed: move scalar products to separate file --- opm/simulators/linalg/ScalarProducts.hpp | 170 +++++++++++++++ opm/simulators/linalg/mixed/SolverAdapter.hpp | 204 +----------------- 2 files changed, 172 insertions(+), 202 deletions(-) create mode 100644 opm/simulators/linalg/ScalarProducts.hpp diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp new file mode 100644 index 00000000000..07a61bce40c --- /dev/null +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -0,0 +1,170 @@ +#ifndef OPM_SCALAR_PRODUCTS_HEADER_INCLUDED +#define OPM_SCALAR_PRODUCTS_HEADER_INCLUDED + +namespace Dune +{ + +template +class GhostLastScalarProduct : public ScalarProduct +{ + public: + + static constexpr auto block_size = Vector::block_type::dimension; + + /*! + * \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(); + int verify = verifyLocalCount(); + if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n"); + } + + /*! + * \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. + It is assumed that the vectors are consistent on the interior+border + partition. + */ + virtual 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_; + + // 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;icommunicator(); + double result = cc.sum(agg[0]); + return result; + } + + /*! \brief Norm of a right-hand side vector. + The vector must be consistent on the interior+border partition + */ + virtual double norm (const Vector& x) const override + { + return sqrt(dot(x,x)); + } + + //! Category of the scalar product (see SolverCategory::Category) + virtual SolverCategory::Category category() const override + { + return _category; + } + + private: + std::shared_ptr _communication; + SolverCategory::Category _category; + int count_; + + 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; + } + + 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 + } + +}; + + + +//! @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 + virtual 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)); + } +}; + +} + +#endif //OPM_SCALAR_PRODUCTS_HEADER_INCLUDED + diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 941d07a4784..994010c6d65 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -9,215 +9,16 @@ #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)); - } -}; - - -template -class GhostLastScalarProduct : public ScalarProduct -{ - public: - - static constexpr auto block_size = Vector::block_type::dimension; - - //! \brief The type of the vector to compute the scalar product on. - //! - //! E.g. BlockVector or another type fulfilling the ISTL - //! vector interface. - - /*! - * \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(); - int verify = verifyLocalCount(); - if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n"); - } - - /*! - * \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. - It is assumed that the vectors are consistent on the interior+border - partition. - */ - void setLocalCount(int count){ count_ = count; } - - virtual 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_; - - //double result(0); - - // 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;icommunicator(); - double result = cc.sum(agg[0]); - return result; - } - - /*! \brief Norm of a right-hand side vector. - The vector must be consistent on the interior+border partition - */ - virtual double norm (const Vector& x) const override - { - return sqrt(dot(x,x)); - } - - //! Category of the scalar product (see SolverCategory::Category) - virtual SolverCategory::Category category() const override - { - return _category; - } - - private: - std::shared_ptr _communication; - SolverCategory::Category _category; - int count_; - - //int getLocalCount(const Comm& comm) const - 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; - } - - 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 - } - -}; - - - -/* -//! @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; - //using type = Dune::AssembledLinearOperator; - //using type = Opm::WellModelMatrixAdapter; -}; -*/ -//! @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 @@ -232,7 +33,6 @@ class MixedBiCGSTABSolver:public InverseOperator static constexpr auto block_size = Vector::block_type::dimension; using MixedMatrixType = Opm::MixedMatrixWrapper; - //using MixedOperatorType = MixedOperator::type; //! @brief constructor //! From f9a6cc260cbf92240d38f5125993921b234da9b8 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Wed, 15 Jul 2026 14:59:24 -0500 Subject: [PATCH 19/39] mixed: scalar products documentation --- opm/simulators/linalg/ScalarProducts.hpp | 52 +++++++++++++++--------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp index 07a61bce40c..4a7612db1cd 100644 --- a/opm/simulators/linalg/ScalarProducts.hpp +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -4,14 +4,20 @@ 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: + ///Exctract block size from vector type static constexpr auto block_size = Vector::block_type::dimension; - /*! + /*! \brief constructor * \param com The communication object for syncing overlap and copy * data points. * \param cat parallel solver category (nonoverlapping or overlapping) @@ -19,12 +25,12 @@ class GhostLastScalarProduct : public ScalarProduct GhostLastScalarProduct (std::shared_ptr com, SolverCategory::Category cat) : _communication(com), _category(cat) { - count_ = getLocalCount(); - int verify = verifyLocalCount(); + 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) @@ -35,8 +41,8 @@ class GhostLastScalarProduct : public ScalarProduct {} /*! \brief Dot product of two vectors. - It is assumed that the vectors are consistent on the interior+border - partition. + * \param vx first input vector + * \param vy second input vector */ virtual double dot (const Vector& vx, const Vector& vy) const override { @@ -58,7 +64,7 @@ class GhostLastScalarProduct : public ScalarProduct for(int j=0;j<2;j++) agg[j]+=agg[j+2]; for(int j=0;j<1;j++) agg[j]+=agg[j+1]; - // trailing end + // loop-peeling of trailing end for(int j=N;j return result; } - /*! \brief Norm of a right-hand side vector. - The vector must be consistent on the interior+border partition + /*! \brief Vector L2-norm. + * \param vx input vector */ - virtual double norm (const Vector& x) const override + virtual double norm (const Vector& vx) const override { - return sqrt(dot(x,x)); + return sqrt(dot(vx,vx)); } //! Category of the scalar product (see SolverCategory::Category) @@ -86,6 +92,8 @@ class GhostLastScalarProduct : public ScalarProduct SolverCategory::Category _category; int count_; + /*! \brief Count number of local cells. + */ int getLocalCount() const { int count = 0; @@ -97,6 +105,8 @@ class GhostLastScalarProduct : public ScalarProduct return count; } + /*! \brief Infer number of local cells from largest local index. + */ int verifyLocalCount() const { auto indexSet = _communication->indexSet(); @@ -121,9 +131,8 @@ class GhostLastScalarProduct : public ScalarProduct -//! @brief Optimized sequential scalar product. -//! -//! @tparam Vector block-vector class with data stored as contiguous double array +/// A sequential scalar product that ignores block structure of the vector +/// to facilitate well-known optimization techniques template class SeqOptmizedProduct : public Dune::SeqScalarProduct { @@ -132,7 +141,10 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct // extract block size static constexpr auto block_size = Vector::block_type::dimension; - // Compute the dot product + /*! \brief Dot product of two vectors. + * \param vx first input vector + * \param vy second input vector + */ virtual double dot(const Vector& vx, const Vector& vy) const override { // access underlying data @@ -152,15 +164,17 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct for(int j=0;j<2;j++) agg[j]+=agg[j+2]; for(int j=0;j<1;j++) agg[j]+=agg[j+1]; - // trailing end + // loop-peeling of trailing end for(int j=N;jdot(x, x)); + /*! \brief Vector L2-norm. + * \param vx input vector + */ + virtual double norm(const Vector& vx) const override { + return std::sqrt(this->dot(vx, vx)); } }; From 463c17558135ed7464155669e53e8c63f45b76fe Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 11:31:32 -0500 Subject: [PATCH 20/39] mixed: matrix wrapper documentation --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 120 +++++++++++++----- 1 file changed, 89 insertions(+), 31 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index cea2753b2ab..dcd4cdd9dd8 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,10 +21,6 @@ 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 @@ -40,9 +35,39 @@ 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 + virtual 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 + virtual 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 + virtual 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;} private: @@ -50,13 +75,21 @@ class MixedMatrixWrapper 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) int const b = block_size; - if constexpr(b==1){printf("MixedMatrixWrapper::mv does not support block size == 1!\n");getchar();} + 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]); @@ -92,30 +125,43 @@ mv(const Vector& x, Vector& y) const } } +//! @brief mixed-precision block-sparse matrix-vector multiplication +//! with update (y += M.x) +//! +//! @note hand-optimized versions are provided for block-sizes +//! 2,3, and 4. A generic implementation for block-sizes > 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) int const b = block_size; - if constexpr(b==1){printf("MixedMatrixWrapper::umv does not support block size == 1!\n");getchar();} + 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 - { - printf("MixedMatrixWrapper::umv only supports block sizes < 4!\n"); - getchar(); - } + 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) int const b = block_size; - if constexpr(b==1){printf("MixedMatrixWrapper::usmv does not support block size == 1!\n");getchar();} + 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); @@ -152,6 +198,12 @@ usmv(double alpha, const Vector& x, Vector& y) const } } +//! @brief update matrix entries +//! +//! @note downcasts from double precision and transposes +//! each non-zero block entry +//! +//! @param data pointer to double precision data template void MixedMatrixWrapper:: update(double const *data) @@ -165,12 +217,17 @@ update(double const *data) for(int i=0;iflt[bb*k + i] = B[i]; } - - // downcast to single precision - //bsr_downcast(M_); } +//! @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 { @@ -193,8 +250,10 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator comm_.project(y); } + // accessor to matix object virtual const M& getmat() const override { return A_; } + // solver category Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::overlapping; @@ -205,15 +264,15 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator const C& comm_; }; - -/*! - \brief Adapter to combine a matrix and another linear operator into - a combined linear operator. - - This is similar to WellModelMatrixAdapter, with the difference that - here we assume a parallel ordering of rows, where ghost rows are - located after interior rows. - */ +//! @brief Adapter to combine a matrix with another linear operator while +//! taking 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 WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator { @@ -245,6 +304,7 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato comm_.project(y); } + // accessor to matix object const M& getmat() const override { return A_; } void addWellPressureEquations(PressureMatrix& jacobian, @@ -266,6 +326,7 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato return wellOper_.getNumberOfExtraEquations(); } + // solver category Dune::SolverCategory::Category category() const override { return Dune::SolverCategory::overlapping; @@ -280,8 +341,5 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato }; - - - } // namespace Opm #endif // OPM_MIXED_MATRIX_HEADER_INCLUDED From 9d08c72755ccd90efcb7078bc7cb0b9c92762730 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 11:56:38 -0500 Subject: [PATCH 21/39] mixed: move linear operators to separate file --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 121 ---------------- opm/simulators/linalg/mixed/Operators.hpp | 130 ++++++++++++++++++ opm/simulators/linalg/mixed/SolverAdapter.hpp | 1 + 3 files changed, 131 insertions(+), 121 deletions(-) create mode 100644 opm/simulators/linalg/mixed/Operators.hpp diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index dcd4cdd9dd8..173d08f06a6 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -220,126 +220,5 @@ update(double const *data) } -//! @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 - virtual void apply( const V& x, V& y ) const override - { - A_.mv(x,y); - comm_.project(y); - } - - // y += \alpha * A * x - virtual void applyscaleadd (double alpha, const V& x, V& y) const override - { - A_.usmv(alpha,x,y); - comm_.project(y); - } - - // accessor to matix object - virtual const M& getmat() const override { return A_; } - - // solver category - Dune::SolverCategory::Category category() const override - { - return Dune::SolverCategory::overlapping; - } - -private: - const M& A_; - const C& comm_; -}; - -//! @brief Adapter to combine a matrix with another linear operator while -//! taking 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 WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator -{ -public: - using field_type = typename V::field_type; - using PressureMatrix = Dune::BCRSMatrix>; - - //! 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 - virtual void apply( const V& x, V& y ) const override - { - A_.mv(x,y); - wellOper_.apply(x, y); - comm_.project(y); - } - - // y += \alpha * A * x - virtual void applyscaleadd (double alpha, const V& x, V& y) const override - { - A_.usmv(alpha,x,y); - wellOper_.applyscaleadd(alpha, x, y); - comm_.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 C& comm_ ; - const LinearOperatorExtra& wellOper_; -}; - - } // 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..ba65342ca05 --- /dev/null +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -0,0 +1,130 @@ +#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 + virtual void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + comm_.project(y); + } + + // y += \alpha * A * x + virtual void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + comm_.project(y); + } + + // accessor to matix object + virtual const M& getmat() const override { return A_; } + + // solver category + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::overlapping; + } + +private: + const M& A_; + const C& comm_; +}; + + +//! @brief Adapter to combine a matrix with another linear operator while +//! taking 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 WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator +{ +public: + using field_type = typename V::field_type; + using PressureMatrix = Dune::BCRSMatrix>; + + //! 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 + virtual void apply( const V& x, V& y ) const override + { + A_.mv(x,y); + wellOper_.apply(x, y); + comm_.project(y); + } + + // y += \alpha * A * x + virtual void applyscaleadd (double alpha, const V& x, V& y) const override + { + A_.usmv(alpha,x,y); + wellOper_.applyscaleadd(alpha, x, y); + comm_.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 C& comm_ ; + const LinearOperatorExtra& wellOper_; +}; + +} // namespace Opm +#endif // OPM_MIXED_OPERATORS_HEADER_INCLUDED diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 994010c6d65..0780bfaa4e1 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -12,6 +12,7 @@ #include #include +#include #include From 3cbcb894816c7a58c9ceba23c7ddf8bf15655ff0 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 13:18:51 -0500 Subject: [PATCH 22/39] mixed: preconditioner documentation --- .../linalg/mixed/PreconditionerWrapper.hpp | 113 +++++++++++++++++- 1 file changed, 109 insertions(+), 4 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 73d89baf311..6bcb1ab2840 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -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_; @@ -85,15 +98,63 @@ class MixedPreconditioner : public Dune::PreconditionerWithUpdate 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 () @@ -109,7 +170,7 @@ update () for(int i=0;idbl[NN*k + i] = B[i]; } - if constexpr(N==1){printf("MixedPreconditioner::update does not support block size == 1!\n");getchar();} + if constexpr(N==1){OPM_THROW(std::invalid_argument, "MixedMatrixPreconditioner::update does not support block size == 1!\n");} else if constexpr(N==2) use_dilu_ ? prec_dilu_factorize2(prec_, mixed_matrix_) : prec_ilu0_factorize2(prec_, mixed_matrix_); else if constexpr(N==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); else if constexpr(N==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); @@ -195,6 +256,14 @@ update () 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) @@ -202,7 +271,7 @@ apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) x=y; int const b = block_size; - if constexpr(b==1){printf("MixedPreconditioner::apply does not support block size == 1!\n");getchar();} + 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]); @@ -248,6 +317,12 @@ apply ([[maybe_unused]] X& x, [[maybe_unused]] const Y& y) } } +//! @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) @@ -263,6 +338,12 @@ matvec_mul(double *y, float const *A, double const * x) for(int i=0;i void MixedPreconditioner:: matvec_mulsub(double *y, float const *A, double const * x) @@ -278,6 +359,10 @@ matvec_mulsub(double *y, float const *A, double const * x) for(int i=0;i void MixedPreconditioner:: mat_copy(double *C, double const * A) @@ -287,6 +372,13 @@ mat_copy(double *C, double const * A) for(int i=0;i void MixedPreconditioner:: mat_inv(double *invA, const double *A) @@ -312,6 +404,11 @@ mat_inv(double *invA, const double *A) mat_copy(invA,T); } +//! @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 template void MixedPreconditioner:: mat_mulsub(double *C, double const *A, double const * B) @@ -330,6 +427,10 @@ 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 template void MixedPreconditioner:: mat_rmul(double *C, double const *A) @@ -349,6 +450,10 @@ mat_rmul(double *C, double const *A) mat_copy(C,T); } +//! @brief In-place matrix-matrix multiplication (C = A.C) +//! +//! @param C column-major input/output matrix +//! @param A right column-major input matrix template void MixedPreconditioner:: mat_lmul(double const *A, double *C) From 61e9b8657e8a6bcc9a43a6a7e94a4de9bbc3b921 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 14:25:59 -0500 Subject: [PATCH 23/39] mixed: solver adapter documentation --- opm/simulators/linalg/mixed/SolverAdapter.hpp | 83 +++++++++---------- 1 file changed, 41 insertions(+), 42 deletions(-) diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 0780bfaa4e1..c2cb751405f 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -52,18 +52,19 @@ class MixedBiCGSTABSolver:public InverseOperator const int& verbosity, const Comm &comm) { -#if 1 int halo; int nrows; int nnz=0; auto &A = op->getmat(); + // 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()]; @@ -76,15 +77,12 @@ class MixedBiCGSTABSolver:public InverseOperator // number of nonzeros for local cells int irow=0; - //for(auto row=A.begin(); row!=A.end(); row++) for(auto row=A.begin(); row.index() < nrows; row++) { if(local_[irow++]==1) for(auto col = row->begin(); col != row->end(); col++) nnz++; - //else break; // This line is used to verify that all local nodes have indices lower than all ghost nodes. This appears to be true! - //irow++; } } - +/* printf("nnz = %d\n",nnz); printf("A.nonzeroes() = %ld\n",A.nonzeroes()); @@ -92,13 +90,8 @@ class MixedBiCGSTABSolver:public InverseOperator printf("halo = %d\n",halo); printf("total = %ld\n",A.N()); //getchar(); -#endif +*/ // Access matrix data from double precision operator -#if 0 - auto &A = op->getmat(); - int nrows = A.N(); - int nnz = A.nonzeroes(); -#endif double_data_ = &A[0][0][0][0]; //allocate mixed matrix @@ -121,42 +114,42 @@ class MixedBiCGSTABSolver:public InverseOperator irow++; } - //initialize mixed operator and optimized scalar product + // 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>) { - //OPM_THROW(std::invalid_argument, "Dune::MatrixAdapter\n"); using MixedOperatorType = Dune::MatrixAdapter; mixed_operator_ = std::make_shared(*mixed_matrix_); - using OptimizedProductType = SeqOptmizedProduct; - scalar_product_ = std::make_shared(); - } - else if constexpr (std::is_same_v, Opm::GhostLastMatrixAdapter>) - { - //OPM_THROW(std::invalid_argument, "Opm::GhostLastMatrixAdapter\n"); - //using MixedOperatorType = Dune::OverlappingSchwarzOperator; - using MixedOperatorType = Opm::MixedGhostLastMatrixAdapter; - mixed_operator_ = std::make_shared(*mixed_matrix_,comm); - using OptimizedScalarProductType = GhostLastScalarProduct; - scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); - //scalar_product_ = sp; + using ScalarProductType = SeqOptmizedProduct; + scalar_product_ = std::make_shared(); } + // serial runs with separate linear operator for wells, i.e. Opm::WellModelMatrixAdapter else if constexpr (std::is_same_v, Opm::WellModelMatrixAdapter>) { - //OPM_THROW(std::invalid_argument, "Opm::WellModelMatrixAdapter\n"); using MixedOperatorType = Opm::WellModelMatrixAdapter; using WellOperatorType = Opm::LinearOperatorExtra; const WellOperatorType &wellOper = op->getwellOper(); mixed_operator_ = std::make_shared(*mixed_matrix_, wellOper); - using OptimizedProductType = SeqOptmizedProduct; - scalar_product_ = std::make_shared(); - //scalar_product_ = sp; + + using ScalarProductType = SeqOptmizedProduct; + scalar_product_ = std::make_shared(); } + // 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>) + { + using MixedOperatorType = Opm::MixedGhostLastMatrixAdapter; + mixed_operator_ = std::make_shared(*mixed_matrix_,comm); + + using ScalarProductType = GhostLastScalarProduct; + scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); + } + // 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>) { - //OPM_THROW(std::invalid_argument, "Opm::WellModelGhostLastMatrixAdapter\n"); using MixedOperatorType = Opm::WellModelMixedGhostLastMatrixAdapter; using WellOperatorType = Opm::LinearOperatorExtra; const WellOperatorType &wellOper = op->getwellOper(); @@ -168,10 +161,12 @@ class MixedBiCGSTABSolver:public InverseOperator } else { - using OptimizedScalarProductType = GhostLastScalarProduct; - scalar_product_ = std::make_shared(comm,Dune::SolverCategory::overlapping); + 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>( @@ -184,6 +179,14 @@ class MixedBiCGSTABSolver:public InverseOperator } + //! @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 @@ -193,6 +196,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; @@ -201,10 +205,14 @@ 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: + //! @brief Count number of ghost cells + //! + //! @param comm communicator object int getHaloCount(const Comm& comm) const { int count = 0; @@ -213,16 +221,7 @@ class MixedBiCGSTABSolver:public InverseOperator for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) { if (idx->local().attribute()!=1) count++; // count ghost indices -/* - // This code snippet is used to check wheter or not all ghost indices occur last in the index set - // In general, that is NOT the case - if (idx->local().attribute()==1) count++; // count local indices - else - { - printf("debug: %d: %ld %d\n",count, idx->local().local(), idx->local().attribute()); - break; - } -*/ + int i=idx->local().local(); // tag local indices local_[i] = (idx->local().attribute()==1) ? 1 : 0; } From 3b11ef9d4751461d10f95fcfb51fc44f6571a8c0 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 15:53:24 -0500 Subject: [PATCH 24/39] mixed: improved naming convention --- opm/simulators/linalg/FlexibleSolver_impl.hpp | 14 +++++++------- opm/simulators/linalg/setupPropertyTree.cpp | 10 +++++----- 2 files changed, 12 insertions(+), 12 deletions(-) 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/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index 76ea086a6d3..ea540561536 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -368,7 +368,7 @@ getSolverString(const FlowLinearSolverParameters& p) else { return {"bicgstab"}; - //return {"mixed-precision"}; + //return {"mixed-bicgstab"}; } } @@ -534,7 +534,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; } @@ -547,7 +547,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; } @@ -561,7 +561,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; } @@ -574,7 +574,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; } From 2e0f0475cf5b13e10befa3232ca8836d36c36c21 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 16 Jul 2026 16:32:12 -0500 Subject: [PATCH 25/39] mixed: update README file --- opm/simulators/linalg/mixed/README.md | 98 ++++++++++++++++++++++----- 1 file changed, 80 insertions(+), 18 deletions(-) diff --git a/opm/simulators/linalg/mixed/README.md b/opm/simulators/linalg/mixed/README.md index 44ae8572bb1..5ecd8e01d01 100644 --- a/opm/simulators/linalg/mixed/README.md +++ b/opm/simulators/linalg/mixed/README.md @@ -1,22 +1,17 @@ # Mixed-precision linear solvers This folder contains mixed-precision building blocks for Krylov subspace methods -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 implementations are specialized for 3x3 block-sparse matrices due to their -importance in reservoir simulation. The original implementation only works in serial. -The parallel implementation wraps the optimized mixed-precision matrix-vector multiplication -and ILU0/DILU preconditioners into building blocks for the ISTL-based implementation of -bicgstab. This sacrifices some performance for improved modularity. Extending the work to -block-sparse matrices of arbitrary block size is work in progress. - -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. +for block-sparse linear systems of equations with highly optimized implementations +for select block-sizes. The mixed-precision sparse matrix-vector multiplications +(SPMV) and preconditioners (ILU0/DILU) are combined to provide ILU0/DILU and +CPR+AMG preconditioned bicgstab algorithms. In the latter case, only the second +stage of the CPR algorithm is performed in mixed-precision. Moreover, the algorithms +leverage an improved scalar product implementation that takes advantage of the fact +that for parallel runs all ghost cells are sorted after local cells. +The current implementations work for both serial and parallel runs and for any block +size > 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 \ --matrix-add-well-contributions=true \ @@ -25,6 +20,73 @@ OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ --linear-solver-max-iter=1024 \ $@ ``` +Currently, a JSON specification file is required to activate mixed-precision CPR+AMG, i.e. +use the wrapper script below +``` +OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ + --linear-solver=../mixed-cprw.json \ + $@ +``` +and modify 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. +Moreover, only block sizes 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 numa --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! From 423871db15282f522a74dae0cbf976dcd4220998 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Sat, 25 Jul 2026 11:00:35 -0500 Subject: [PATCH 26/39] mixed: replace project calls in linear operators --- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 2 ++ opm/simulators/linalg/mixed/Operators.hpp | 35 ++++++++++++++++--- 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 173d08f06a6..8076d666d0b 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -70,6 +70,8 @@ class MixedMatrixWrapper //! @brief access column index pointer int *colidx(){return M_->colidx;} + int nrows() const {return M_->nrows;} + private: int nnz_; bsr_matrix *M_; diff --git a/opm/simulators/linalg/mixed/Operators.hpp b/opm/simulators/linalg/mixed/Operators.hpp index ba65342ca05..b9172b98bbc 100644 --- a/opm/simulators/linalg/mixed/Operators.hpp +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -17,6 +17,8 @@ template class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator { public: + // extract block size + static constexpr auto block_size = V::block_type::dimension; //! constructor: just store a reference to matrix and communicator MixedGhostLastMatrixAdapter (const M& A, const C& comm) : A_( A ), comm_(comm) {} @@ -25,14 +27,14 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator virtual void apply( const V& x, V& y ) const override { A_.mv(x,y); - comm_.project(y); + ghostLast_project(y); } // y += \alpha * A * x virtual void applyscaleadd (double alpha, const V& x, V& y) const override { A_.usmv(alpha,x,y); - comm_.project(y); + ghostLast_project(y); } // accessor to matix object @@ -45,6 +47,15 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator } 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>; + // 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, @@ -79,7 +93,8 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato { A_.mv(x,y); wellOper_.apply(x, y); - comm_.project(y); + //comm_.project(y); + ghostLast_project(y); } // y += \alpha * A * x @@ -87,7 +102,8 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato { A_.usmv(alpha,x,y); wellOper_.applyscaleadd(alpha, x, y); - comm_.project(y); + //comm_.project(y); + ghostLast_project(y); } // accessor to matix object @@ -124,6 +140,17 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato const M& A_ ; const C& comm_ ; const LinearOperatorExtra& wellOper_; + +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 Date: Sat, 25 Jul 2026 11:02:43 -0500 Subject: [PATCH 27/39] mixed: move scalar product to c for avx2 support --- CMakeLists_files.cmake | 1 + opm/simulators/linalg/ScalarProducts.hpp | 14 ++++++-- opm/simulators/linalg/mixed/dot.c | 46 ++++++++++++++++++++++++ opm/simulators/linalg/mixed/dot.h | 12 +++++++ 4 files changed, 71 insertions(+), 2 deletions(-) create mode 100644 opm/simulators/linalg/mixed/dot.c create mode 100644 opm/simulators/linalg/mixed/dot.h 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/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp index 4a7612db1cd..d0b20ad4b57 100644 --- a/opm/simulators/linalg/ScalarProducts.hpp +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -1,6 +1,8 @@ #ifndef OPM_SCALAR_PRODUCTS_HEADER_INCLUDED #define OPM_SCALAR_PRODUCTS_HEADER_INCLUDED +#include + namespace Dune { @@ -53,7 +55,7 @@ class GhostLastScalarProduct : public ScalarProduct // total array length int NN = block_size*count_; - +#if 0 // unroll loop in multiples of 8 int n=NN/8; int N=8*n; @@ -71,6 +73,11 @@ class GhostLastScalarProduct : public ScalarProduct auto cc = _communication->communicator(); double result = cc.sum(agg[0]); return result; +#else + auto cc = _communication->communicator(); + return cc.sum(vec_dot(x,y,NN)); +#endif + //return cc.sum(vec_dot(x,y,NN)); } /*! \brief Vector L2-norm. @@ -153,7 +160,7 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct // total array length int NN = block_size*vx.N(); - +#if 1 // unroll loop in multiples of 8 int n=NN/8; int N=8*n; @@ -168,6 +175,9 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct for(int j=N;j + +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 Date: Mon, 27 Jul 2026 20:40:43 -0500 Subject: [PATCH 28/39] mixed: enable mixed-cprw solver option --- opm/simulators/linalg/setupPropertyTree.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index ea540561536..cb9d23ca394 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; @@ -380,14 +380,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); From 81c7e379d71a45846d9ae3cfe4e427d1c993d29b Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 28 Jul 2026 00:01:38 -0500 Subject: [PATCH 29/39] mixed: update README file --- opm/simulators/linalg/mixed/README.md | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/opm/simulators/linalg/mixed/README.md b/opm/simulators/linalg/mixed/README.md index 5ecd8e01d01..4d814a22786 100644 --- a/opm/simulators/linalg/mixed/README.md +++ b/opm/simulators/linalg/mixed/README.md @@ -20,14 +20,24 @@ OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ --linear-solver-max-iter=1024 \ $@ ``` -Currently, a JSON specification file is required to activate mixed-precision CPR+AMG, i.e. -use the wrapper script below +Similarly, a sample wrapper for running the simulator with mixed-precision CPR+AMG++BiCGSTAB +is given by +``` bash +OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --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 1 --map-by numa --bind-to core build/bin/flow \ --linear-solver=../mixed-cprw.json \ $@ ``` -and modify the following `mixed-cprw.json` file to your liking +and modifying the following `mixed-cprw.json` file to your liking ``` { "maxiter": "1024", From 0899deb1f244aeaebfcb9dcde828288af877f953 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Fri, 21 Aug 2026 09:47:59 -0500 Subject: [PATCH 30/39] mixed: documentation and clean-up --- opm/simulators/linalg/mixed/SolverAdapter.hpp | 10 +- opm/simulators/linalg/mixed/bslv.c | 52 +-- opm/simulators/linalg/mixed/bsr.c | 1 - opm/simulators/linalg/mixed/bsr.h | 1 - opm/simulators/linalg/mixed/dot.c | 12 +- opm/simulators/linalg/mixed/matvec.h | 428 +++++++++--------- opm/simulators/linalg/mixed/prec.c | 113 ----- opm/simulators/linalg/mixed/prec.h | 2 - 8 files changed, 219 insertions(+), 400 deletions(-) diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index c2cb751405f..e51ed746fc0 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -82,15 +82,7 @@ class MixedBiCGSTABSolver:public InverseOperator if(local_[irow++]==1) for(auto col = row->begin(); col != row->end(); col++) nnz++; } } -/* - printf("nnz = %d\n",nnz); - printf("A.nonzeroes() = %ld\n",A.nonzeroes()); - - printf("local = %d\n",nrows); - printf("halo = %d\n",halo); - printf("total = %ld\n",A.N()); - //getchar(); -*/ + // Access matrix data from double precision operator double_data_ = &A[0][0][0][0]; diff --git a/opm/simulators/linalg/mixed/bslv.c b/opm/simulators/linalg/mixed/bslv.c index 2fa2fef4254..ed3c69f57bf 100644 --- a/opm/simulators/linalg/mixed/bslv.c +++ b/opm/simulators/linalg/mixed/bslv.c @@ -103,7 +103,6 @@ double __attribute__((noinline)) vec_inner2(const double *a, const double *b, in int bslv_pbicgstab3m(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; @@ -128,53 +127,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;jtol; int max_iter = mem->max_iter; int n = mem->n; @@ -199,46 +190,39 @@ int bslv_pbicgstab4m(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;jtol; int max_iter = mem->max_iter; int n = mem->n; @@ -271,46 +254,39 @@ 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;jnnz; diff --git a/opm/simulators/linalg/mixed/bsr.h b/opm/simulators/linalg/mixed/bsr.h index 5384dba8c20..33787061eb7 100644 --- a/opm/simulators/linalg/mixed/bsr.h +++ b/opm/simulators/linalg/mixed/bsr.h @@ -79,7 +79,6 @@ void bsr_vmspumv4(bsr_matrix *A, const double *x, double *y, double alpha); * @param x Pointer to input vector. * @param y Pointer to output vector. */ - void bsr_vmspmv2(bsr_matrix *A, const double *x, double *y); void bsr_vmspmv3(bsr_matrix *A, const double *x, double *y); void bsr_vmspmv4(bsr_matrix *A, const double *x, double *y); diff --git a/opm/simulators/linalg/mixed/dot.c b/opm/simulators/linalg/mixed/dot.c index ba952936403..49eeb06c181 100644 --- a/opm/simulators/linalg/mixed/dot.c +++ b/opm/simulators/linalg/mixed/dot.c @@ -4,7 +4,6 @@ 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; @@ -23,9 +22,9 @@ double vec_dot (double const *x, double const *y, int NN) double vec_bdot (double const *x, double const *y, int NN) { - // unroll loop in multiples of 8 - //int n=NN/8; + // assumes vectors are padded by zeros + // to nearest multiple of 8 doubles int N=8*((NN+7)/8); double agg[8]; for(int i=0;i<8;i++) agg[i]=0.0; @@ -34,13 +33,6 @@ double vec_bdot (double const *x, double const *y, int NN) for(int j=0;j<4;j++) agg[j]+=agg[j+4]; for(int j=0;j<2;j++) agg[j]+=agg[j+2]; for(int j=0;j<1;j++) agg[j]+=agg[j+1]; -/* - // loop-peeling of trailing end - for(int j=N;j [0 1 4 5] - row1 = _mm_loadh_pi(_mm_loadl_pi(row1, (__m64*)(A+8)), (__m64*)(A +12)); // [8 9 x x] -> [8 9 C D] - - row0 = _mm_shuffle_ps(tmp1, row1, 0x88); // [0 4 8 C] 0b10 00 10 00 - row1 = _mm_shuffle_ps(row1, tmp1, 0xDD); // [9 D 1 5] 0b11 01 11 01 -*/ - tmp1 = _mm_unpacklo_ps(col0,col1); // [0 4 1 5] - row1 = _mm_unpacklo_ps(col2,col3); // [8 C 9 D] - - row0 = _mm_shuffle_ps(tmp1, row1, 0x44); // [0 4 8 C] 0b 01 00 01 00 - row1 = _mm_shuffle_ps(row1, tmp1, 0xEE); // [9 D 1 5] 0b 11 10 11 10 // notice flipped order -/* - tmp1 = _mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(A+ 2)), (__m64*)(A+ 6)); // [2 3 x x] -> [2 3 6 7] - row3 = _mm_loadh_pi(_mm_loadl_pi(row3, (__m64*)(A+10)), (__m64*)(A+14)); // [A B x x] -> [A B E F] - - row2 = _mm_shuffle_ps(tmp1, row3, 0x88); // [2 6 A E] 0b 10 00 10 00 - row3 = _mm_shuffle_ps(row3, tmp1, 0xDD); // [B F 3 7] 0b 11 01 11 01 -*/ - tmp1 = _mm_unpackhi_ps(col0,col1); // [2 6 3 7] 0b 01 00 01 00 - row3 = _mm_unpackhi_ps(col2,col3); // [A E B F] 0b 11 10 11 10 - - row2 = _mm_shuffle_ps(tmp1, row3, 0x44); // [2 6 A E] 0b 01 00 01 00 - row3 = _mm_shuffle_ps(row3, tmp1, 0xEE); // [B F 3 7] 0b 11 10 11 10 // notice flipped order - // ----------------------------------------------- - tmp1 = _mm_mul_ps(row2, row3); // [2B 6F A3 E7] - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); // 1011 0001 [2B 6F A3 E7] -> [6F 2B E7 A3] - - col0 = _mm_mul_ps(row1, tmp1); // [96F D2B 1E7 5A3] - col1 = _mm_mul_ps(row0, tmp1); // [06F 42B 8E7 CA3] - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); // 0100 1110 [2B 6F A3 E7] -> [A3 E7 2B 6F] - - col0 = _mm_sub_ps(_mm_mul_ps(row1, tmp1), col0);// [9A3 DE7 12B 56F] - [96F D2B 1E7 5A3] = [9(A3-6F) D(E7-2B) 1(2B-E7) 5(6F-A3)] - col1 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col1);// [0A3 4E7 82B C6F] - [06F 42B 8E7 CA3] = [0(A3-6F) 4(E7-2B) 8(2B-E7) C(6F-A3)] - col1 = _mm_shuffle_ps(col1, col1, 0x4E); // [8(2B-E7) C(6F-A3) 0(A3-6F) 4(E7-2B)] - // ----------------------------------------------- - tmp1 = _mm_mul_ps(row1, row2); - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); - - col0 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col0); - col3 = _mm_mul_ps(row0, tmp1); - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); - - col0 = _mm_sub_ps(col0, _mm_mul_ps(row3, tmp1)); - col3 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col3); - col3 = _mm_shuffle_ps(col3, col3, 0x4E); - // ----------------------------------------------- - tmp1 = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3); - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); - row2 = _mm_shuffle_ps(row2, row2, 0x4E); - - col0 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col0); - col2 = _mm_mul_ps(row0, tmp1); - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); - - col0 = _mm_sub_ps(col0, _mm_mul_ps(row2, tmp1)); - col2 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col2); - col2 = _mm_shuffle_ps(col2, col2, 0x4E); - // ----------------------------------------------- - tmp1 = _mm_mul_ps(row0, row1); - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); - - col2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col2); - col3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), col3); - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); - - col2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), col2); - col3 = _mm_sub_ps(col3, _mm_mul_ps(row2, tmp1)); - // ----------------------------------------------- - tmp1 = _mm_mul_ps(row0, row3); - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); - - col1 = _mm_sub_ps(col1, _mm_mul_ps(row2, tmp1)); - col2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col2); - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); - - col1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col1); - col2 = _mm_sub_ps(col2, _mm_mul_ps(row1, tmp1)); - // ----------------------------------------------- - tmp1 = _mm_mul_ps(row0, row2); - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); - - col1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col1); - col3 = _mm_sub_ps(col3, _mm_mul_ps(row1, tmp1)); - - tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); - - col1 = _mm_sub_ps(col1, _mm_mul_ps(row3, tmp1)); - col3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col3); -// ----------------------------------------------- - det = _mm_mul_ps(row0, col0); - det = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det); - det = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det); - tmp1 = _mm_rcp_ss(det); - - det = _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1))); - det = _mm_shuffle_ps(det, det, 0x00); - -/* - _mm_store_ps(A+ 0, col0); - _mm_store_ps(A+ 4, col1); - _mm_store_ps(A+ 8, col2); - _mm_store_ps(A+12, col3); - return; -*/ -/* - col0 = _mm_mul_ps(det, col0); - _mm_storel_pi((__m64*)(A+0), col0); - _mm_storeh_pi((__m64*)(A+2), col0); - - col1 = _mm_mul_ps(det, col1); - _mm_storel_pi((__m64*)(A+4), col1); - _mm_storeh_pi((__m64*)(A+6), col1); - - col2 = _mm_mul_ps(det, col2); - _mm_storel_pi((__m64*)(A+ 8), col2); - _mm_storeh_pi((__m64*)(A+10), col2); - - col3 = _mm_mul_ps(det, col3); - _mm_storel_pi((__m64*)(A+12), col3); - _mm_storeh_pi((__m64*)(A+14), col3); -*/ + double M[16]; + for(int k=0;k<16;k++) M[k] = A[k]; - _mm_store_ps(A+ 0,_mm_mul_ps(det, col0)); - _mm_store_ps(A+ 4,_mm_mul_ps(det, col1)); - _mm_store_ps(A+ 8,_mm_mul_ps(det, col2)); - _mm_store_ps(A+12,_mm_mul_ps(det, col3)); + for(int k=0;k<4;k++) + { + double scale=-1.0/M[5*k]; + for(int i=0;i<4;i++) M[i+4*k] *= i==k?0:scale; // scale column k + for(int j=0;j<4;j++) + { + if (j==k) continue; + for(int i=0;i<4;i++) M[i+4*j] += i==k?0:M[i+4*k]*M[k+4*j]; //sweep + } + scale=-scale; + for(int j=0;j<4;j++) M[k+4*j] *= scale; // scale row k + M[5*k] = scale; + } + for(int k=0;k<16;k++) invA[k] = M[k]; } -// AVX2 double-precision translation of Intel's SSE single-precision 4x4 matrix inverse using Cramer's Rule + +/** + * @brief Matrix inverse for 4x4 matrix. + * + * @note This is an avx2 double-precision translation of Intel's + * SSE single-precision 4x4 matrix inverse using Cramer's + * rule. + * + * @param invA Pointer to inverse matrix. + * @param A Pointer to input matrix. + */ void mat4_vinv(double *invA, double const *A) { __m256d col0, col1, col2, col3; __m256d row0, row1, row2, row3; __m256d tmp1, det; -/* - 0 4 8 C - 1 5 9 D - 2 6 A E - 3 7 B F -*/ - // ----------------------------------------------- - // extract rows from columns - // ----------------------------------------------- + + // load columns col0 = _mm256_load_pd(A+ 0); col1 = _mm256_load_pd(A+ 4); col2 = _mm256_load_pd(A+ 8); col3 = _mm256_load_pd(A+12); - tmp1 = _mm256_unpacklo_pd(col0,col1); // [0 4 1 5] //all maps have to be updated due to per-lane unpack - row2 = _mm256_unpacklo_pd(col2,col3); // [8 C 9 D] + tmp1 = _mm256_unpacklo_pd(col0,col1); + row2 = _mm256_unpacklo_pd(col2,col3); - row0 = _mm256_permute2f128_pd(tmp1, row2, 0x20); // [0 4 8 C] 0b 00 10 00 00 extract lower pairs - row2 = _mm256_permute2f128_pd(row2, tmp1, 0x13); // [9 D 1 5] 0b 00 01 00 11 extract and swap upper pairs + row0 = _mm256_permute2f128_pd(tmp1, row2, 0x20); + row2 = _mm256_permute2f128_pd(row2, tmp1, 0x13); - tmp1 = _mm256_unpackhi_pd(col0,col1); // [2 6 3 7] 0b 01 00 01 00 - row1 = _mm256_unpackhi_pd(col2,col3); // [A E B F] 0b 11 10 11 10 + tmp1 = _mm256_unpackhi_pd(col0,col1); + row1 = _mm256_unpackhi_pd(col2,col3); - row3 = _mm256_permute2f128_pd(tmp1, row1, 0x13); // [2 6 A E] 0b 01 00 01 00 - row1 = _mm256_permute2f128_pd(row1, tmp1, 0x20); // [B F 3 7] 0b 11 10 11 10 // notice flipped order - // ----------------------------------------------- - tmp1 = _mm256_mul_pd(row2, row3); // [2B 6F A3 E7] - tmp1 = _mm256_permute_pd(tmp1, 0x05); // 1011 0001 [2B 6F A3 E7] -> [6F 2B E7 A3] reverse upper and lower pairs + row3 = _mm256_permute2f128_pd(tmp1, row1, 0x13); + row1 = _mm256_permute2f128_pd(row1, tmp1, 0x20); + // --------------------------------------------- + tmp1 = _mm256_mul_pd(row2, row3); + tmp1 = _mm256_permute_pd(tmp1, 0x05); - col0 = _mm256_mul_pd(row1, tmp1); // [96F D2B 1E7 5A3] - col1 = _mm256_mul_pd(row0, tmp1); // [06F 42B 8E7 CA3] + col0 = _mm256_mul_pd(row1, tmp1); + col1 = _mm256_mul_pd(row0, tmp1); tmp1 = _mm256_permute4x64_pd(tmp1,0x4E); - col0 = _mm256_fmsub_pd(row1,tmp1,col0); // [9A3 DE7 12B 56F] - [96F D2B 1E7 5A3] = [9(A3-6F) D(E7-2B) 1(2B-E7) 5(6F-A3)] - col1 = _mm256_fmsub_pd(row0,tmp1,col1); // [0A3 4E7 82B C6F] - [06F 42B 8E7 CA3] = [0(A3-6F) 4(E7-2B) 8(2B-E7) C(6F-A3)] + col0 = _mm256_fmsub_pd(row1,tmp1,col0); + col1 = _mm256_fmsub_pd(row0,tmp1,col1); col1 = _mm256_permute4x64_pd(col1,0x4E); // ----------------------------------------------- tmp1 = _mm256_mul_pd(row1, row2); @@ -498,6 +355,125 @@ void mat4_vinv(double *invA, double const *A) _mm256_store_pd(invA+12,_mm256_mul_pd(det, col3)); } + +/** + * @brief Matrix inverse for 4x4 matrix. + * + * @note This is a slightly modifed version of Intel's SSE + * single-precision 4x4 matrix inverse using Cramer's + * rule. + * + * @param invA Pointer to inverse matrix. + * @param A Pointer to input matrix. + */ +void mat4_inv2(float *A) +{ + __m128 col0, col1, col2, col3; + __m128 row0, row1, row2, row3; + __m128 tmp1, det; + + // load columns + col0 = _mm_load_ps(A+ 0); + col1 = _mm_load_ps(A+ 4); + col2 = _mm_load_ps(A+ 8); + col3 = _mm_load_ps(A+12); + + tmp1 = _mm_unpacklo_ps(col0,col1); + row1 = _mm_unpacklo_ps(col2,col3); + + row0 = _mm_shuffle_ps(tmp1, row1, 0x44); + row1 = _mm_shuffle_ps(row1, tmp1, 0xEE); + + tmp1 = _mm_unpackhi_ps(col0,col1); + row3 = _mm_unpackhi_ps(col2,col3); + + row2 = _mm_shuffle_ps(tmp1, row3, 0x44); + row3 = _mm_shuffle_ps(row3, tmp1, 0xEE); + // --------------------------------------------- + tmp1 = _mm_mul_ps(row2, row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col0 = _mm_mul_ps(row1, tmp1); + col1 = _mm_mul_ps(row0, tmp1); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col0 = _mm_sub_ps(_mm_mul_ps(row1, tmp1), col0); + col1 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col1); + col1 = _mm_shuffle_ps(col1, col1, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row1, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col0 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col0); + col3 = _mm_mul_ps(row0, tmp1); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col0 = _mm_sub_ps(col0, _mm_mul_ps(row3, tmp1)); + col3 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col3); + col3 = _mm_shuffle_ps(col3, col3, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + row2 = _mm_shuffle_ps(row2, row2, 0x4E); + + col0 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col0); + col2 = _mm_mul_ps(row0, tmp1); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col0 = _mm_sub_ps(col0, _mm_mul_ps(row2, tmp1)); + col2 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), col2); + col2 = _mm_shuffle_ps(col2, col2, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col2); + col3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), col3); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), col2); + col3 = _mm_sub_ps(col3, _mm_mul_ps(row2, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col1 = _mm_sub_ps(col1, _mm_mul_ps(row2, tmp1)); + col2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col2); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), col1); + col2 = _mm_sub_ps(col2, _mm_mul_ps(row1, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + + col1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), col1); + col3 = _mm_sub_ps(col3, _mm_mul_ps(row1, tmp1)); + + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + + col1 = _mm_sub_ps(col1, _mm_mul_ps(row3, tmp1)); + col3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), col3); +// ----------------------------------------------- + det = _mm_mul_ps(row0, col0); + det = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det); + det = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det); + tmp1 = _mm_rcp_ss(det); + + det = _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1))); + det = _mm_shuffle_ps(det, det, 0x00); + + _mm_store_ps(A+ 0,_mm_mul_ps(det, col0)); + _mm_store_ps(A+ 4,_mm_mul_ps(det, col1)); + _mm_store_ps(A+ 8,_mm_mul_ps(det, col2)); + _mm_store_ps(A+12,_mm_mul_ps(det, col3)); +} + #ifdef __cplusplus } #endif diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index d2ad01ce650..b1ca1d283d3 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -772,51 +772,7 @@ void prec_ilu0_factorize4(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) -{ - const int b=3; - double z[3]; - for(int k=0;k<3;k++) z[k]=0; - for(int c=0;cL; @@ -1107,72 +1063,3 @@ void prec_info(prec_t *P) bsr_info(P->U); } -void prec_test() -{ -//#if 0 - // verify 2x2 inverse and matrix-matrix multiplications - double A[4] = {1,0.2,0.3,4}; - double B[4] = {1,0.2,0.3,4}; - double C[4] = {1,0.2,0.3,4}; - double I[4] = {1,0,0,1}; - - mat_show(A,2,"A"); - - mat2_inv(A,A); - mat_show(A,2,"A"); - - mat2_vfms(I,A,B); - mat_show(I,2,"I"); - - mat2_rmul(B,A); - mat_show(B,2,"B"); - - mat2_lmul(A,C); - mat_show(C,2,"C"); -//#endif -#if 0 - // verify 4x4 inverse and matrix-matrix multiplications - double AA[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double BB[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double CC[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double II[16] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; - mat_show(AA,4,"AA"); - - mat4_inv2(AA,AA); - mat_show(AA,4,"AA"); - - mat4_vfms(II,AA,BB); - mat_show(II,4,"II"); - - mat4_rmul(BB,AA); - mat_show(BB,4,"BB"); - - mat4_lmul(AA,CC); - mat_show(CC,4,"CC"); -#endif -#if 0 - // verify 4x4 inverse and matrix-matrix multiplications - double A[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double B[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double C[16] = {1,0.2,0.3,0.4, 0.5,6,0.7,0.8, 0.9,1.0,11,1.2, 1.3,1.4,1.5,16}; - double I[16] = {1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1}; -/* - mat_fshow(A,4,"A"); - mat4_inv2(A); - mat_fshow(A,4,"A"); -*/ - mat_show(A,4,"A"); - mat4_vinv(A,A); - mat_show(A,4,"A"); - - mat4_vfms(I,A,B); - mat_show(I,4,"I"); - - mat4_rmul(B,A); - mat_show(B,4,"B"); - - mat4_lmul(A,C); - mat_show(C,4,"C"); -#endif -} - diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index cca6c114804..9c1f9752538 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -113,8 +113,6 @@ void prec_downcast(prec_t *P); */ void prec_info(prec_t *P); - -void prec_test(); #ifdef __cplusplus } #endif From 56752452ded8d83ffbf5125f3dcb4f3b13ab6edc Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Fri, 21 Aug 2026 12:04:36 -0500 Subject: [PATCH 31/39] mixed: delete unnecessary virtual keywords --- opm/simulators/linalg/ScalarProducts.hpp | 2 +- opm/simulators/linalg/mixed/MatrixWrapper.hpp | 6 +++--- opm/simulators/linalg/mixed/Operators.hpp | 10 +++++----- opm/simulators/linalg/mixed/SolverAdapter.hpp | 3 ++- opm/simulators/linalg/mixed/matvec.h | 3 +++ 5 files changed, 14 insertions(+), 10 deletions(-) diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp index d0b20ad4b57..37a647ea25f 100644 --- a/opm/simulators/linalg/ScalarProducts.hpp +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -160,7 +160,7 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct // total array length int NN = block_size*vx.N(); -#if 1 +#if 0 // unroll loop in multiples of 8 int n=NN/8; int N=8*n; diff --git a/opm/simulators/linalg/mixed/MatrixWrapper.hpp b/opm/simulators/linalg/mixed/MatrixWrapper.hpp index 8076d666d0b..8ca328332f1 100644 --- a/opm/simulators/linalg/mixed/MatrixWrapper.hpp +++ b/opm/simulators/linalg/mixed/MatrixWrapper.hpp @@ -47,14 +47,14 @@ class MixedMatrixWrapper //! //! @param x input vector //! @param y output vector - virtual void mv(const Vector& x, Vector& y) const; + 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 - virtual void umv(const Vector& x, Vector& y) const; + void umv(const Vector& x, Vector& y) const; //! @brief block-sparse matrix-vector multiplication with //! scaled update (y += alpha * M.x) @@ -62,7 +62,7 @@ class MixedMatrixWrapper //! @param alpha scaling factor //! @param x input vector //! @param y output vector - virtual void usmv(double alpha, const Vector& x, Vector& y) const; + void usmv(double alpha, const Vector& x, Vector& y) const; //! @brief access row offset pointer int *rowptr(){return M_->rowptr;} diff --git a/opm/simulators/linalg/mixed/Operators.hpp b/opm/simulators/linalg/mixed/Operators.hpp index b9172b98bbc..19861c71108 100644 --- a/opm/simulators/linalg/mixed/Operators.hpp +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -24,21 +24,21 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator MixedGhostLastMatrixAdapter (const M& A, const C& comm) : A_( A ), comm_(comm) {} // y = A * x - virtual void apply( const V& x, V& y ) const override + void apply( const V& x, V& y ) const override { A_.mv(x,y); ghostLast_project(y); } // y += \alpha * A * x - virtual void applyscaleadd (double alpha, const V& x, V& y) const override + void applyscaleadd (double alpha, const V& x, V& y) const override { A_.usmv(alpha,x,y); ghostLast_project(y); } // accessor to matix object - virtual const M& getmat() const override { return A_; } + const M& getmat() const override { return A_; } // solver category Dune::SolverCategory::Category category() const override @@ -89,7 +89,7 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato {} // y = A * x - virtual void apply( const V& x, V& y ) const override + void apply( const V& x, V& y ) const override { A_.mv(x,y); wellOper_.apply(x, y); @@ -98,7 +98,7 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato } // y += \alpha * A * x - virtual void applyscaleadd (double alpha, const V& x, V& y) const override + void applyscaleadd (double alpha, const V& x, V& y) const override { A_.usmv(alpha,x,y); wellOper_.applyscaleadd(alpha, x, y); diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index e51ed746fc0..1f72c6c5d82 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -53,7 +53,8 @@ class MixedBiCGSTABSolver:public InverseOperator const Comm &comm) { int halo; - int nrows; + //int nrows; + size_t nrows; int nnz=0; auto &A = op->getmat(); diff --git a/opm/simulators/linalg/mixed/matvec.h b/opm/simulators/linalg/mixed/matvec.h index 11cefa64a73..4bd692e99cc 100644 --- a/opm/simulators/linalg/mixed/matvec.h +++ b/opm/simulators/linalg/mixed/matvec.h @@ -4,6 +4,8 @@ extern "C" { #endif +// Developer helpers to visualize matrices +#if 0 void mat_fshow(float const *A, int n, char const *name) { printf("%s = [\n",name); @@ -25,6 +27,7 @@ void mat_show(double const *A, int n, char const *name) } printf("]\n\n"); } +#endif /** * @brief In-place right matrix-matrix multiplication for 2x2 matrices. From fa38182198b38d40501e2d3c488a3cda715bd93e Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Fri, 21 Aug 2026 13:38:13 -0500 Subject: [PATCH 32/39] mixed: minor improvements --- opm/simulators/linalg/ScalarProducts.hpp | 25 +++++++++++-------- opm/simulators/linalg/mixed/Operators.hpp | 5 ++-- opm/simulators/linalg/mixed/SolverAdapter.hpp | 2 +- opm/simulators/linalg/mixed/wrapper.hpp | 7 +++--- 4 files changed, 22 insertions(+), 17 deletions(-) diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp index 37a647ea25f..6de6a4c2637 100644 --- a/opm/simulators/linalg/ScalarProducts.hpp +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -16,9 +16,6 @@ class GhostLastScalarProduct : public ScalarProduct { public: - ///Exctract block size from vector type - static constexpr auto block_size = Vector::block_type::dimension; - /*! \brief constructor * \param com The communication object for syncing overlap and copy * data points. @@ -46,7 +43,7 @@ class GhostLastScalarProduct : public ScalarProduct * \param vx first input vector * \param vy second input vector */ - virtual double dot (const Vector& vx, const Vector& vy) const override + double dot (const Vector& vx, const Vector& vy) const override { // access underlying data @@ -83,7 +80,7 @@ class GhostLastScalarProduct : public ScalarProduct /*! \brief Vector L2-norm. * \param vx input vector */ - virtual double norm (const Vector& vx) const override + double norm (const Vector& vx) const override { return sqrt(dot(vx,vx)); } @@ -95,6 +92,10 @@ class GhostLastScalarProduct : public ScalarProduct } 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_; @@ -145,14 +146,11 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct { public: - // extract block size - static constexpr auto block_size = Vector::block_type::dimension; - /*! \brief Dot product of two vectors. * \param vx first input vector * \param vy second input vector */ - virtual double dot(const Vector& vx, const Vector& vy) const override + double dot(const Vector& vx, const Vector& vy) const override { // access underlying data double const *x = &vx[0][0]; @@ -183,9 +181,16 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct /*! \brief Vector L2-norm. * \param vx input vector */ - virtual double norm(const Vector& vx) const override { + 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; + + }; } diff --git a/opm/simulators/linalg/mixed/Operators.hpp b/opm/simulators/linalg/mixed/Operators.hpp index 19861c71108..8e0b32d989c 100644 --- a/opm/simulators/linalg/mixed/Operators.hpp +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -17,9 +17,6 @@ template class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator { public: - // extract block size - static constexpr auto block_size = V::block_type::dimension; - //! constructor: just store a reference to matrix and communicator MixedGhostLastMatrixAdapter (const M& A, const C& comm) : A_( A ), comm_(comm) {} @@ -47,6 +44,8 @@ class MixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperator } private: + // extract block size + static constexpr auto block_size = V::block_type::dimension; void ghostLast_project( V& y ) const { diff --git a/opm/simulators/linalg/mixed/SolverAdapter.hpp b/opm/simulators/linalg/mixed/SolverAdapter.hpp index 1f72c6c5d82..335b331f97e 100644 --- a/opm/simulators/linalg/mixed/SolverAdapter.hpp +++ b/opm/simulators/linalg/mixed/SolverAdapter.hpp @@ -80,7 +80,7 @@ class MixedBiCGSTABSolver:public InverseOperator int irow=0; for(auto row=A.begin(); row.index() < nrows; row++) { - if(local_[irow++]==1) for(auto col = row->begin(); col != row->end(); col++) nnz++; + nnz += local_[irow++] ? std::distance(row->begin(), row->end()) : 0; } } diff --git a/opm/simulators/linalg/mixed/wrapper.hpp b/opm/simulators/linalg/mixed/wrapper.hpp index c8e0ee37077..23906480edd 100644 --- a/opm/simulators/linalg/mixed/wrapper.hpp +++ b/opm/simulators/linalg/mixed/wrapper.hpp @@ -14,9 +14,6 @@ class MixedSolver : public InverseOperator { public: - // extract block size - static constexpr auto block_size = X::block_type::dimension; - MixedSolver(const M &A, double tol, int maxiter, bool use_dilu) { // verify that well contributions are added to the matrix @@ -104,6 +101,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_; From 023cf73dcc476dee4afdec59aa5461206e1702b2 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Fri, 21 Aug 2026 21:10:38 -0500 Subject: [PATCH 33/39] mixed: more clean-up --- opm/simulators/linalg/ScalarProducts.hpp | 40 +------------------ opm/simulators/linalg/mixed/SolverAdapter.hpp | 1 - 2 files changed, 2 insertions(+), 39 deletions(-) diff --git a/opm/simulators/linalg/ScalarProducts.hpp b/opm/simulators/linalg/ScalarProducts.hpp index 6de6a4c2637..4d26a866717 100644 --- a/opm/simulators/linalg/ScalarProducts.hpp +++ b/opm/simulators/linalg/ScalarProducts.hpp @@ -52,29 +52,9 @@ class GhostLastScalarProduct : public ScalarProduct // total array length int NN = block_size*count_; -#if 0 - // 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;icommunicator(); - double result = cc.sum(agg[0]); - return result; -#else + auto cc = _communication->communicator(); return cc.sum(vec_dot(x,y,NN)); -#endif - //return cc.sum(vec_dot(x,y,NN)); } /*! \brief Vector L2-norm. @@ -158,24 +138,8 @@ class SeqOptmizedProduct : public Dune::SeqScalarProduct // total array length int NN = block_size*vx.N(); -#if 0 - // 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 const Comm &comm) { int halo; - //int nrows; size_t nrows; int nnz=0; From 802d6c6ceb58c798017814de138285ad98cd78fc Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Fri, 21 Aug 2026 23:11:45 -0500 Subject: [PATCH 34/39] mixed: aligned and buffered allocations --- opm/simulators/linalg/mixed/bsr.c | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/opm/simulators/linalg/mixed/bsr.c b/opm/simulators/linalg/mixed/bsr.c index f88da4fe56d..49234209cdd 100644 --- a/opm/simulators/linalg/mixed/bsr.c +++ b/opm/simulators/linalg/mixed/bsr.c @@ -38,6 +38,14 @@ void bsr_free(bsr_matrix *A) A=NULL; } +inline void * buffered_alloc(size_t alignment, size_t size) +{ + // round up to nearest integer multiple of alignment + size_t real_size = alignment*((size + alignment - 1)/alignment); + + return aligned_alloc(alignment, real_size); +} + void bsr_init(bsr_matrix *A, int nrows, int nnz, int b) { A->nrows=nrows; @@ -48,8 +56,8 @@ 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 = aligned_alloc(64,b*b*nnz*sizeof(double)); - A->flt = aligned_alloc(64,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); From 46bebb4b5415ec8eebe1a900c2dcfd67add8f0c4 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Mon, 24 Aug 2026 13:32:38 -0500 Subject: [PATCH 35/39] mixed: consolidate ilu0 and dilu factorization --- .../linalg/mixed/PreconditionerWrapper.hpp | 6 +- opm/simulators/linalg/mixed/bslv.c | 6 +- opm/simulators/linalg/mixed/prec.c | 223 +----------------- opm/simulators/linalg/mixed/prec.h | 7 +- 4 files changed, 16 insertions(+), 226 deletions(-) diff --git a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp index 6bcb1ab2840..55ed976f7fb 100644 --- a/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp +++ b/opm/simulators/linalg/mixed/PreconditionerWrapper.hpp @@ -171,9 +171,9 @@ update () } if constexpr(N==1){OPM_THROW(std::invalid_argument, "MixedMatrixPreconditioner::update does not support block size == 1!\n");} - else if constexpr(N==2) use_dilu_ ? prec_dilu_factorize2(prec_, mixed_matrix_) : prec_ilu0_factorize2(prec_, mixed_matrix_); - else if constexpr(N==3) use_dilu_ ? prec_dilu_factorize(prec_, mixed_matrix_) : prec_ilu0_factorize(prec_, mixed_matrix_); - else if constexpr(N==4) use_dilu_ ? prec_dilu_factorize4(prec_, mixed_matrix_) : prec_ilu0_factorize4(prec_, mixed_matrix_); + 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_; diff --git a/opm/simulators/linalg/mixed/bslv.c b/opm/simulators/linalg/mixed/bslv.c index ed3c69f57bf..633e48b50f6 100644 --- a/opm/simulators/linalg/mixed/bslv.c +++ b/opm/simulators/linalg/mixed/bslv.c @@ -119,7 +119,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 + prec_ilu0_factorize3(P,A,mem->use_dilu); // choose dilu or ilu0 prec_downcast(P); vec_fill(x_j,0.0,n); @@ -182,7 +182,7 @@ int bslv_pbicgstab4m(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_factorize4(P,A) : prec_ilu0_factorize4(P,A); // choose dilu or ilu0 + prec_ilu0_factorize4(P,A,mem->use_dilu); // choose dilu or ilu0 prec_downcast(P); vec_fill(x_j,0.0,n); @@ -246,7 +246,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); diff --git a/opm/simulators/linalg/mixed/prec.c b/opm/simulators/linalg/mixed/prec.c index b1ca1d283d3..d61916b3d5d 100644 --- a/opm/simulators/linalg/mixed/prec.c +++ b/opm/simulators/linalg/mixed/prec.c @@ -325,218 +325,7 @@ void mat3_vfms(double *C, double const *A, double const *B) } } -void prec_dilu_factorize2(prec_t *P, bsr_matrix *A) -{ - int nrows = A->nrows; - int b = A->b; - int bb = b*b; - - 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++) - { - int j=A->colidx[k]; - if(jrowptr[j]; - vec_copy4(L->dbl + bb*kL, A->dbl + bb*k); - L->rowptr[j]++; - } - else if(j==i) // struct-copy of D - { - vec_copy4(D->dbl + bb*i, A->dbl + bb*k); - } - else if(j>i) // struct-copy of U - { - vec_copy4(U->dbl + bb*kU, A->dbl + bb*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 - double scale[4]; //hard-coded to 2x2 blocks - for(int i=0;inrows;i++) - { - 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 - mat2_rmul(L->dbl+k*bb,scale); - - //update diagonal of U - int j=L->colidx[k]; - mat2_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); - - //scale row i of U - mat2_lmul(scale,U->dbl+k*bb); - - //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(); - } - } - } - } -} - - -void prec_dilu_factorize(prec_t *P, bsr_matrix *A) -{ - int nrows = A->nrows; - int b = A->b; - int bb = b*b; - - 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++) - { - int j=A->colidx[k]; - if(jrowptr[j]; - vec_copy9(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); - } - else if(j>i) // struct-copy of U - { - vec_copy9(U->dbl + bb*kU, A->dbl + bb*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 - double scale[9]; //hard-coded to 3x3 blocks for now - 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 - for(int k=L->rowptr[i];krowptr[i+1];k++) - { - //scale column i of L - mat3_rmul(L->dbl+k*bb,scale); - - //update diagonal of U - int j=L->colidx[k]; - mat3_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); - - //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(); - } - } - } - } -} - -void prec_dilu_factorize4(prec_t *P, bsr_matrix *A) -{ - int nrows = A->nrows; - int b = A->b; - int bb = b*b; - - 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++) - { - 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++; - } - } - } - // reset rowptr of L - for(int i=nrows;i>0;i--) L->rowptr[i]=L->rowptr[i-1]; - L->rowptr[0]=0; - - // Factorizing - double scale[16]; //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); - - //update diagonal of U - int j=L->colidx[k]; - mat4_vfms(D->dbl+j*bb,L->dbl+k*bb,U->dbl+k*bb); - - //scale row i of U - mat4_lmul(scale,U->dbl+k*bb); - - //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(); - } - } - } - } -} - -void prec_ilu0_factorize2(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; @@ -576,7 +365,7 @@ void prec_ilu0_factorize2(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[4]; //hard-coded to 2x2 blocks for(int i=0;inrows;i++) { @@ -615,7 +404,7 @@ void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A) } } -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; @@ -655,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++) { @@ -694,7 +483,7 @@ void prec_ilu0_factorize(prec_t *P, bsr_matrix *A) } } -void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A) +void prec_ilu0_factorize4(prec_t *P, bsr_matrix *A, bool use_dilu) { int nrows = A->nrows; @@ -735,7 +524,7 @@ void prec_ilu0_factorize4(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[16] __attribute__((aligned(64))); //hard-coded to 4x4 blocks for(int i=0;inrows;i++) { diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index 9c1f9752538..cd4e2941f99 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. @@ -73,9 +74,9 @@ void prec_dilu_factorize4(prec_t *P, bsr_matrix *A); * @param P Pointer preconditioner object. * @param A Pointer to bsr matrix. */ -void prec_ilu0_factorize2(prec_t *P, bsr_matrix *A); -void prec_ilu0_factorize(prec_t *P, bsr_matrix *A); -void prec_ilu0_factorize4(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. From 72a166d41931eec2e5a72049a8a56ca0a67ae0fa Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 25 Aug 2026 09:29:33 -0500 Subject: [PATCH 36/39] mixed: generic legacy bicgstab implementation --- opm/simulators/linalg/mixed/bslv.c | 106 +++++++----------------- opm/simulators/linalg/mixed/bslv.h | 12 ++- opm/simulators/linalg/mixed/prec.h | 12 +-- opm/simulators/linalg/mixed/wrapper.hpp | 6 +- 4 files changed, 44 insertions(+), 92 deletions(-) diff --git a/opm/simulators/linalg/mixed/bslv.c b/opm/simulators/linalg/mixed/bslv.c index 633e48b50f6..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,7 +125,7 @@ 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; @@ -119,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; - prec_ilu0_factorize3(P,A,mem->use_dilu); // 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); @@ -134,15 +158,15 @@ int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x 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_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_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_pbicgstab4m(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; - - double * restrict e = mem->e; - const double *r0 = b; - //const double * restrict r0 = mem->dtmp[0]; //access randomly initialized one-dimensional shadow space - double * restrict p_j = mem->dtmp[1]; - double * restrict q_j = mem->dtmp[2]; - double * restrict r_j = mem->dtmp[3]; - double * restrict s_j = mem->dtmp[4]; - double * restrict t_j = mem->dtmp[5]; - double * restrict v_j = mem->dtmp[6]; - double * restrict x_j = x; - - prec_t * restrict P = mem->P; - prec_ilu0_factorize4(P,A,mem->use_dilu); // choose dilu or ilu0 - prec_downcast(P); - - vec_fill(x_j,0.0,n); - vec_copy(r_j,b,n); - vec_copy(p_j,b,n); - - vec_copy(q_j,p_j,n); - double norm_0 = sqrt(vec_inner2(r_j,r_j,n)); - - double rho_j = vec_inner2(r0,r_j,n); - int j; - for(j=0;jtol; @@ -292,7 +252,3 @@ int bslv_pbicgstab3d(bslv_memory *mem, bsr_matrix *A, const double *b, double *x } - - - - diff --git a/opm/simulators/linalg/mixed/bslv.h b/opm/simulators/linalg/mixed/bslv.h index ca4e5af4a83..550ec26fc9b 100644 --- a/opm/simulators/linalg/mixed/bslv.h +++ b/opm/simulators/linalg/mixed/bslv.h @@ -33,6 +33,15 @@ struct bslv_memory // pointer to preconditioner prec_t *P; + + // spmv function pointer + void (*bsr_spmv)(bsr_matrix *A, const double *x, double *y); + + // ILU0/DILU factorization function pointer + void (*prec_factorize)(prec_t *P, bsr_matrix *A, bool use_dilu); + + // ILU0/DILU apply function pointer + void (*prec_apply)(prec_t *P, double *x); } bslv_memory; @@ -81,8 +90,7 @@ void bslv_init(bslv_memory *mem, double tol, int max_iter, bsr_matrix const *A, * * @return Number of linear iterations. */ -int bslv_pbicgstab3m(bslv_memory *mem, bsr_matrix *A, const double *b, double *x); -int bslv_pbicgstab4m(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); /** * @brief Preconditioned bicgstab in double-precision. diff --git a/opm/simulators/linalg/mixed/prec.h b/opm/simulators/linalg/mixed/prec.h index cd4e2941f99..4c58a528fe5 100644 --- a/opm/simulators/linalg/mixed/prec.h +++ b/opm/simulators/linalg/mixed/prec.h @@ -59,17 +59,7 @@ void prec_init(prec_t *P, bsr_matrix const *A); int prec_analyze(bsr_matrix *M, int (*offsets)[3]); /** - * @brief DILU factorization. - * - * @param P Pointer preconditioner object. - * @param A Pointer to bsr matrix. - */ -void prec_dilu_factorize2(prec_t *P, bsr_matrix *A); -void prec_dilu_factorize(prec_t *P, bsr_matrix *A); -void prec_dilu_factorize4(prec_t *P, bsr_matrix *A); - -/** - * @brief ILU0 factorization. + * @brief ILU0/DILU factorization. * * @param P Pointer preconditioner object. * @param A Pointer to bsr matrix. diff --git a/opm/simulators/linalg/mixed/wrapper.hpp b/opm/simulators/linalg/mixed/wrapper.hpp index 23906480edd..14d60c9ade9 100644 --- a/opm/simulators/linalg/mixed/wrapper.hpp +++ b/opm/simulators/linalg/mixed/wrapper.hpp @@ -26,7 +26,7 @@ class MixedSolver : public InverseOperator int b = A[0][0].N(); // verify that block size is 3x3 or 4x4 - if (b<3 || b>4) {OPM_THROW(std::logic_error, "Legacy mixed precision only supports 3x3 and 4x4 blocks.");} + 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(); @@ -79,9 +79,7 @@ class MixedSolver : public InverseOperator bsr_downcast(jacobian_); // solve linear system - int count = 0; - if constexpr(N==3) count = bslv_pbicgstab3m(mem_, jacobian_, &b[0][0], &x[0][0]); - else if constexpr(N==4) count = bslv_pbicgstab4m(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); From 42fd03a2cf4ed2011da47b68328b8ff04ca035cf Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Tue, 25 Aug 2026 09:30:13 -0500 Subject: [PATCH 37/39] mixed: removing vec_bdot, i.e. buffered dot product --- opm/simulators/linalg/mixed/dot.c | 16 ---------------- opm/simulators/linalg/mixed/dot.h | 1 - 2 files changed, 17 deletions(-) diff --git a/opm/simulators/linalg/mixed/dot.c b/opm/simulators/linalg/mixed/dot.c index 49eeb06c181..2762d50139e 100644 --- a/opm/simulators/linalg/mixed/dot.c +++ b/opm/simulators/linalg/mixed/dot.c @@ -20,19 +20,3 @@ double vec_dot (double const *x, double const *y, int NN) return agg[0]; } -double vec_bdot (double const *x, double const *y, int NN) -{ - // unroll loop in multiples of 8 - // assumes vectors are padded by zeros - // to nearest multiple of 8 doubles - int N=8*((NN+7)/8); - double agg[8]; - for(int i=0;i<8;i++) agg[i]=0.0; - for(int i=0;i Date: Mon, 31 Aug 2026 09:32:07 -0500 Subject: [PATCH 38/39] mixed: match declaration and init order --- opm/simulators/linalg/mixed/Operators.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/simulators/linalg/mixed/Operators.hpp b/opm/simulators/linalg/mixed/Operators.hpp index 8e0b32d989c..62a5d2de804 100644 --- a/opm/simulators/linalg/mixed/Operators.hpp +++ b/opm/simulators/linalg/mixed/Operators.hpp @@ -137,8 +137,8 @@ class WellModelMixedGhostLastMatrixAdapter : public Dune::AssembledLinearOperato protected: const M& A_ ; - const C& comm_ ; const LinearOperatorExtra& wellOper_; + const C& comm_ ; private: From cb60954ac53dc1b2e18039541df7688b5b99a562 Mon Sep 17 00:00:00 2001 From: Kjetil B Haugen Date: Thu, 3 Sep 2026 09:23:48 -0500 Subject: [PATCH 39/39] mixed: update README file --- opm/simulators/linalg/mixed/README.md | 18 ++++++++++-------- opm/simulators/linalg/setupPropertyTree.cpp | 1 - 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/opm/simulators/linalg/mixed/README.md b/opm/simulators/linalg/mixed/README.md index 4d814a22786..9f208ab4a01 100644 --- a/opm/simulators/linalg/mixed/README.md +++ b/opm/simulators/linalg/mixed/README.md @@ -1,9 +1,9 @@ # Mixed-precision linear solvers This folder contains mixed-precision building blocks for Krylov subspace methods for block-sparse linear systems of equations with highly optimized implementations -for select block-sizes. The mixed-precision sparse matrix-vector multiplications +for select block-sizes. The mixed-precision sparse matrix-vector multiplications (SPMV) and preconditioners (ILU0/DILU) are combined to provide ILU0/DILU and -CPR+AMG preconditioned bicgstab algorithms. In the latter case, only the second +CPR+AMG preconditioned BiCGSTAB algorithms. In the latter case, only the second stage of the CPR algorithm is performed in mixed-precision. Moreover, the algorithms leverage an improved scalar product implementation that takes advantage of the fact that for parallel runs all ghost cells are sorted after local cells. @@ -13,7 +13,7 @@ size > 1. However, only block-sizes 2,3, and 4 benefit from hand-optimized imple 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 \ @@ -23,7 +23,7 @@ OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ Similarly, a sample wrapper for running the simulator with mixed-precision CPR+AMG++BiCGSTAB is given by ``` 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=false \ --linear-solver=mixed-cprw \ --linear-solver-reduction=1e-3 \ @@ -33,7 +33,7 @@ OMP_NUM_THREADS=1 mpirun -np 1 --map-by numa --bind-to core build/bin/flow \ 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 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 \ --linear-solver=../mixed-cprw.json \ $@ ``` @@ -88,10 +88,12 @@ and modifying the following `mixed-cprw.json` file to your liking 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. -Moreover, only block sizes 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 +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 numa --bind-to core build/bin/flow \ +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 \ diff --git a/opm/simulators/linalg/setupPropertyTree.cpp b/opm/simulators/linalg/setupPropertyTree.cpp index cb9d23ca394..41fd56c2433 100644 --- a/opm/simulators/linalg/setupPropertyTree.cpp +++ b/opm/simulators/linalg/setupPropertyTree.cpp @@ -368,7 +368,6 @@ getSolverString(const FlowLinearSolverParameters& p) else { return {"bicgstab"}; - //return {"mixed-bicgstab"}; } }