diff --git a/MathLib/LinAlg/Eigen/EigenTools.cpp b/MathLib/LinAlg/Eigen/EigenTools.cpp index 3321d1b67fd13c47353617585514095acd57f9f2..3d1e43c215734047b6a1121b98565bb67a61b4aa 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