From f9efd72ecc5b4da5653402c0622e8374bd150178 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Wed, 17 Oct 2012 12:40:43 +0200
Subject: [PATCH] Fix output in case of LAPACK error.

Make copy of matrix before calling dgesv, since it will overwrite it.
---
 opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp | 4 +++-
 opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp | 1 +
 opm/core/utility/VelocityInterpolation.cpp                    | 4 +++-
 3 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
index 718f2cf6..404452e0 100644
--- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
+++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
@@ -465,6 +465,7 @@ namespace Opm
         tof_coeff_ = &tof_coeff[0];
         rhs_.resize(num_basis);
         jac_.resize(num_basis*num_basis);
+        orig_jac_.resize(num_basis*num_basis);
         basis_.resize(num_basis);
         basis_nb_.resize(num_basis);
         grad_basis_.resize(num_basis*grid_.dimensions);
@@ -617,6 +618,7 @@ namespace Opm
         std::vector<MAT_SIZE_T> piv(num_basis);
         MAT_SIZE_T ldb = num_basis;
         MAT_SIZE_T info = 0;
+        orig_jac_ = jac_;
         dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
         if (info != 0) {
             // Print the local matrix and rhs.
@@ -624,7 +626,7 @@ namespace Opm
                       << " with A = \n";
             for (int row = 0; row < n; ++row) {
                 for (int col = 0; col < n; ++col) {
-                    std::cerr << "    " << jac_[row + n*col];
+                    std::cerr << "    " << orig_jac_[row + n*col];
                 }
                 std::cerr << '\n';
             }
diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
index 68307e1c..9a2e3b9c 100644
--- a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
+++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
@@ -91,6 +91,7 @@ namespace Opm
         double* tof_coeff_;
         std::vector<double> rhs_;   // single-cell right-hand-side
         std::vector<double> jac_;   // single-cell jacobian
+        std::vector<double> orig_jac_;   // single-cell jacobian (copy)
         // Below: storage for quantities needed by solveSingleCell().
         std::vector<double> coord_;
         std::vector<double> basis_;
diff --git a/opm/core/utility/VelocityInterpolation.cpp b/opm/core/utility/VelocityInterpolation.cpp
index 76b9aba7..b235fc77 100644
--- a/opm/core/utility/VelocityInterpolation.cpp
+++ b/opm/core/utility/VelocityInterpolation.cpp
@@ -178,6 +178,7 @@ namespace Opm
         // for each corner.
         const int dim = grid_.dimensions;
         std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering!
+        std::vector<double> orig_N(dim*dim); // Normals matrix. Fortran ordering!
         std::vector<double> f(dim);     // Flux vector.
         std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve
         const int num_cells = grid_.number_of_cells;
@@ -203,6 +204,7 @@ namespace Opm
                 MAT_SIZE_T lda = n;
                 MAT_SIZE_T ldb = n;
                 MAT_SIZE_T info = 0;
+                orig_N = N;
                 dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info);
                 if (info != 0) {
                     // Print the local matrix and rhs.
@@ -210,7 +212,7 @@ namespace Opm
                               << " with N = \n";
                     for (int row = 0; row < n; ++row) {
                         for (int col = 0; col < n; ++col) {
-                            std::cerr << "    " << N[row + n*col];
+                            std::cerr << "    " << orig_N[row + n*col];
                         }
                         std::cerr << '\n';
                     }