Skip to content
Snippets Groups Projects
Commit 183dc717 authored by Christoph Lehmann's avatar Christoph Lehmann
Browse files

[MaL] faster apply known solution

parent 8806fc19
No related branches found
No related tags found
1 merge request!888[MaL] faster applyKnownSolution
......@@ -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
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment