OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Ensure correct output in case of Lapack error.

Browse Source 
Make sure that we output the original values, since Lapack may overwrite those
used in the call.
This commit is contained in:
Atgeirr Flø Rasmussen 2012-12-17 11:07:40 +01:00
parent 389a943f34
commit b66920203d
1 changed files with 3 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
opm/core/utility/VelocityInterpolation.cpp
View File

@ -101,6 +101,7 @@ namespace Opm
std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering! std::vector<double> N(dim*dim); // Normals matrix. Fortran ordering!
std::vector<double> orig_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<double> f(dim); // Flux vector.
std::vector<double> orig_f(dim); // Flux vector.
std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve std::vector<MAT_SIZE_T> piv(dim); // For LAPACK solve
const SparseTable<WachspressCoord::CornerInfo>& all_ci = bcmethod_.cornerInfo(); const SparseTable<WachspressCoord::CornerInfo>& all_ci = bcmethod_.cornerInfo();
const std::vector<int>& adj_faces = bcmethod_.adjacentFaces(); const std::vector<int>& adj_faces = bcmethod_.adjacentFaces();
@ -129,6 +130,7 @@ namespace Opm
MAT_SIZE_T ldb = n; MAT_SIZE_T ldb = n;
MAT_SIZE_T info = 0; MAT_SIZE_T info = 0;
orig_N = N; orig_N = N;
orig_f = f;
dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info); dgesv_(&n, &nrhs, &N[0], &lda, &piv[0], &f[0], &ldb, &info);
if (info != 0) { if (info != 0) {
// Print the local matrix and rhs. // Print the local matrix and rhs.
@ -142,7 +144,7 @@ namespace Opm
} }
std::cerr << "and f = \n"; std::cerr << "and f = \n";
for (int row = 0; row < n; ++row) { for (int row = 0; row < n; ++row) {
std::cerr << " " << f[row] << '\n'; std::cerr << " " << orig_f[row] << '\n';
} }
THROW("Lapack error: " << info << " encountered in cell " << cell); THROW("Lapack error: " << info << " encountered in cell " << cell);
} }