From 183dc717612d7f68449373917d63c14f8d3f8a58 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Wed, 11 Nov 2015 12:54:46 +0100
Subject: [PATCH] [MaL] faster apply known solution
---
MathLib/LinAlg/Eigen/EigenTools.cpp | 43 +++++++++++++++--------------
1 file changed, 22 insertions(+), 21 deletions(-)
diff --git a/MathLib/LinAlg/Eigen/EigenTools.cpp b/MathLib/LinAlg/Eigen/EigenTools.cpp
index 3321d1b67fd..3d1e43c2157 100644
--- a/MathLib/LinAlg/Eigen/EigenTools.cpp
+++ b/MathLib/LinAlg/Eigen/EigenTools.cpp
@@ -21,35 +21,36 @@ void applyKnownSolution(EigenMatrix &A_, EigenVector &b_, const std::vector<std:
const std::vector<double> &vec_knownX_x, double /*penalty_scaling*/)
{
using SpMat = EigenMatrix::RawMatrixType;
+ static_assert(SpMat::IsRowMajor, "matrix is assumed to be row major!");
+
auto &A = A_.getRawMatrix();
auto &b = b_.getRawVector();
- const std::size_t n_rows = A.rows();
+
+ // A(k, j) = 0.
+ // set row to zero
+ for (auto row_id : vec_knownX_id)
+ for (SpMat::InnerIterator it(A,row_id); it; ++it) it.valueRef() = 0.0;
+
+ SpMat AT = A.transpose();
+
for (std::size_t ix=0; ix<vec_knownX_id.size(); ix++)
{
SpMat::Index const row_id = vec_knownX_id[ix];
auto const x = vec_knownX_x[ix];
- //A(k, j) = 0.
- for (SpMat::InnerIterator it(A,row_id); it; ++it)
- it.valueRef() = .0;
- //b_i -= A(i,k)*val, i!=k
- for (std::size_t i=0; i<n_rows; i++)
- for (SpMat::InnerIterator it(A,i); it; ++it)
- {
- if (it.col()!=row_id) continue;
- b[i] -= it.value()*x;
- }
- //b_k = val
+
+ // b_i -= A(i,k)*val, i!=k
+ // set column to zero, subtract from rhs
+ for (SpMat::InnerIterator it(AT, row_id); it; ++it)
+ {
+ b[it.col()] -= it.value()*x;
+ it.valueRef() = 0.0;
+ }
+
b[row_id] = x;
- //A(i, k) = 0., i!=k
- for (std::size_t i=0; i<n_rows; i++)
- for (SpMat::InnerIterator it(A,i); it; ++it)
- {
- if (it.col()!=row_id) continue;
- it.valueRef() = 0.0;
- }
- //A(k, k) = 1.0
- A.coeffRef(row_id, row_id) = 1.0; //=x
+ AT.coeffRef(row_id, row_id) = 1.0;
}
+
+ A = AT.transpose();
}
} // MathLib
--
GitLab