Merge pull request #888 from chleh/eigen-known-solution

[MaL] faster applyKnownSolution

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