Port unsupported constrained CG to Eigen3
(grafted from 4cd4be97a7165e6e45ee60aee23b9342af03c491
)
diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
index b83bf7a..d3c4c65 100644
--- a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
+++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
@@ -62,7 +62,9 @@
Scalar rho, rho_1, alpha;
d.setZero();
- CINV.startFill(); // FIXME estimate the number of non-zeros
+ typedef Triplet<double> T;
+ std::vector<T> tripletList;
+
for (Index i = 0; i < rows; ++i)
{
d[i] = 1.0;
@@ -88,11 +90,12 @@
// FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse
for (Index j=0; j<l.size(); ++j)
if (l[j]<1e-15)
- CINV.fill(i,j) = l[j];
+ tripletList.push_back(T(i,j,l(j)));
+
d[i] = 0.0;
}
- CINV.endFill();
+ CINV.setFromTriplets(tripletList.begin(), tripletList.end());
}
@@ -107,6 +110,7 @@
void constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x,
const VectorB& b, const VectorF& f, IterationController &iter)
{
+ using std::sqrt;
typedef typename TMatrix::Scalar Scalar;
typedef typename TMatrix::Index Index;
typedef Matrix<Scalar,Dynamic,1> TmpVec;
