From 3115b7868fb094e3adfa45d2d6f5c7403d97a25c Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Tue, 25 Mar 2014 09:58:07 +0100
Subject: [PATCH] Moved ISTL right and left hand side construction for ISTL to
 solveSystem

cherry-picked from add-scalarproduct and ported.
---
 opm/core/linalg/LinearSolverIstl.cpp | 47 ++++++++++++++--------------
 opm/core/linalg/LinearSolverIstl.hpp | 10 +++---
 2 files changed, 28 insertions(+), 29 deletions(-)

diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp
index 822782904..5355d05b4 100644
--- a/opm/core/linalg/LinearSolverIstl.cpp
+++ b/opm/core/linalg/LinearSolverIstl.cpp
@@ -154,30 +154,12 @@ namespace Opm
                 A[ri][ja[i]] = sa[i];
             }
         }
-        // System RHS
-        Vector b(size);
-        std::copy(rhs, rhs + size, b.begin());
-        // System solution
-        Vector x(size);
-        x = 0.0;
-
-        if (linsolver_save_system_)
-        {
-            // Save system to files.
-            writeMatrixToMatlab(A, linsolver_save_filename_ + "-mat");
-            std::string rhsfile(linsolver_save_filename_ + "-rhs");
-            std::ofstream rhsf(rhsfile.c_str());
-            rhsf.precision(15);
-            rhsf.setf(std::ios::scientific | std::ios::showpos);
-            std::copy(b.begin(), b.end(),
-                      std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
-        }
 
         int maxit = linsolver_max_iterations_;
         if (maxit == 0) {
             maxit = 5000;
         }
-#if HAVE_MPI
+        #if HAVE_MPI
         if(comm.type()==typeid(ParallelISTLInformation))
         {
             typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
@@ -187,7 +169,7 @@ namespace Opm
             Dune::OverlappingSchwarzOperator<Mat,Vector,Vector, Comm>
                 opA(A, istlComm);
             Dune::OverlappingSchwarzScalarProduct<Vector,Comm> sp(istlComm);
-            return solveSystem(opA, x, solution, b, sp, istlComm, maxit);
+            return solveSystem(opA, solution, rhs, sp, istlComm, maxit);
         }
         else
 #endif
@@ -195,15 +177,34 @@ namespace Opm
             Dune::SeqScalarProduct<Vector> sp;
             Dune::Amg::SequentialInformation comm;
             Operator opA(A);
-            return solveSystem(opA, x, solution, b, sp, comm, maxit);
+            return solveSystem(opA, solution, rhs, sp, comm, maxit);
         }
     }
 
-template<class O, class V, class S, class C>
+    template<class O, class S, class C>
     LinearSolverInterface::LinearSolverReport
-    LinearSolverIstl::solveSystem (O& opA, V& x, double* solution, V& b,
+    LinearSolverIstl::solveSystem (O& opA, double* solution, const double* rhs,
                                    S& sp, const C& comm, int maxit) const
     {
+                // System RHS
+        Vector b(opA.getmat().N());
+        std::copy(rhs, rhs + size, b.begin());
+        // System solution
+        Vector x(opA.getmat().M());
+        x = 0.0;
+
+        if (linsolver_save_system_)
+        {
+            // Save system to files.
+            writeMatrixToMatlab(opA.getmat(), linsolver_save_filename_ + "-mat");
+            std::string rhsfile(linsolver_save_filename_ + "-rhs");
+            std::ofstream rhsf(rhsfile.c_str());
+            rhsf.precision(15);
+            rhsf.setf(std::ios::scientific | std::ios::showpos);
+            std::copy(b.begin(), b.end(),
+                      std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
+        }
+
         LinearSolverReport res;
         switch (linsolver_type_) {
         case CG_ILU0:
diff --git a/opm/core/linalg/LinearSolverIstl.hpp b/opm/core/linalg/LinearSolverIstl.hpp
index de32783f3..4a948c087 100644
--- a/opm/core/linalg/LinearSolverIstl.hpp
+++ b/opm/core/linalg/LinearSolverIstl.hpp
@@ -88,16 +88,14 @@ namespace Opm
     private:
         /// \brief Solve the linear system using ISTL
         /// \param[in] opA The linear operator of the system to solve.
-        /// \param[in,out] x The vector with the initial guess that will store
-        ///                 the computed solution
         /// \param[out]    solution C array for storing the solution vector.
-        /// \param[in]     b The vector containing the right hand side of the system.
+        /// \param[in]     rhs C array containing the right hand side.
         /// \param[in]     sp The scalar product to use.
         /// \param[in]     comm The information about the parallel domain decomposition.
         /// \param[in]     maxit The maximum number of iterations allowed.
-        template<class O, class V, class S, class C>
-        LinearSolverReport solveSystem(O& opA, V& x, double* solution,
-                                       V& b, S& sp, const C& comm, int maxit) const;
+        template<class O, class S, class C>
+        LinearSolverReport solveSystem(O& opA, double* solution, const double *rhs,
+                                       S& sp, const C& comm, int maxit) const;
 
         double linsolver_residual_tolerance_;
         int linsolver_verbosity_;