From 6ca5f167b2ee38fce97b979c4ecef950b5c7c24e Mon Sep 17 00:00:00 2001
From: Tong Dong Qiu <tongdong.qiu@bigdataccelerate.com>
Date: Fri, 25 Feb 2022 11:06:50 +0100
Subject: [PATCH] Allow BISAI to use jacMatrix

---
 opm/simulators/linalg/bda/opencl/BISAI.cpp | 44 +++++++++++++---------
 opm/simulators/linalg/bda/opencl/BISAI.hpp |  5 +--
 2 files changed, 29 insertions(+), 20 deletions(-)

diff --git a/opm/simulators/linalg/bda/opencl/BISAI.cpp b/opm/simulators/linalg/bda/opencl/BISAI.cpp
index fcf04d2395..695045ca59 100644
--- a/opm/simulators/linalg/bda/opencl/BISAI.cpp
+++ b/opm/simulators/linalg/bda/opencl/BISAI.cpp
@@ -79,23 +79,29 @@ std::vector<int> buildCsrToCscOffsetMap(std::vector<int> colPointers, std::vecto
 template <unsigned int block_size>
 bool BISAI<block_size>::analyze_matrix(BlockedMatrix *mat)
 {
-    const unsigned int bs = block_size;
-
-    this->N = mat->Nb * bs;
-    this->Nb = mat->Nb;
-    this->nnz = mat->nnzbs * bs * bs;
-    this->nnzb = mat->nnzbs;
-
-    bilu0->analyze_matrix(mat);
-
-    return true;
+    return analyze_matrix(mat, nullptr);
 }
 
 template <unsigned int block_size>
 bool BISAI<block_size>::analyze_matrix(BlockedMatrix *mat, BlockedMatrix *jacMat)
 {
-    (void) jacMat;
-    return analyze_matrix(mat);
+    const unsigned int bs = block_size;
+    auto *m = mat;
+
+    if (jacMat) {
+        m = jacMat;
+    }
+
+    this->N = m->Nb * bs;
+    this->Nb = m->Nb;
+    this->nnz = m->nnzbs * bs * bs;
+    this->nnzb = m->nnzbs;
+
+    if (jacMat) {
+        return bilu0->analyze_matrix(mat, jacMat);
+    } else {
+        return bilu0->analyze_matrix(mat);
+    }
 }
 
 template <unsigned int block_size>
@@ -169,7 +175,7 @@ void BISAI<block_size>::buildUpperSubsystemsStructures(){
 }
 
 template <unsigned int block_size>
-bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat)
+bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat, BlockedMatrix *jacMat)
 {
     const unsigned int bs = block_size;
 
@@ -179,7 +185,11 @@ bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat)
 
     Dune::Timer t_preconditioner;
 
-    bilu0->create_preconditioner(mat);
+    if (jacMat) {
+        bilu0->create_preconditioner(mat, jacMat);
+    } else {
+        bilu0->create_preconditioner(mat);
+    }
 
     std::call_once(initialize, [&]() {
         std::tie(colPointers, rowIndices, diagIndex) = bilu0->get_preconditioner_structure();
@@ -261,11 +271,11 @@ bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat)
 
     return true;
 }
+
 template <unsigned int block_size>
-bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat, BlockedMatrix *jacMat)
+bool BISAI<block_size>::create_preconditioner(BlockedMatrix *mat)
 {
-    (void) jacMat;
-    return create_preconditioner(mat);
+    return create_preconditioner(mat, nullptr);
 }
 
 template <unsigned int block_size>
diff --git a/opm/simulators/linalg/bda/opencl/BISAI.hpp b/opm/simulators/linalg/bda/opencl/BISAI.hpp
index bcaace1b97..f0fe852422 100644
--- a/opm/simulators/linalg/bda/opencl/BISAI.hpp
+++ b/opm/simulators/linalg/bda/opencl/BISAI.hpp
@@ -117,13 +117,12 @@ public:
 
     // analysis, find reordering if specified
     bool analyze_matrix(BlockedMatrix *mat) override;
-    
     bool analyze_matrix(BlockedMatrix *mat, BlockedMatrix *jacMat) override;
+
     // ilu_decomposition
     bool create_preconditioner(BlockedMatrix *mat) override;
-
     bool create_preconditioner(BlockedMatrix *mat, BlockedMatrix *jacMat) override;
-    
+
     // apply preconditioner, x = prec(y)
     void apply(const cl::Buffer& y, cl::Buffer& x) override;