33 #include "../../../../Eigen/src/Core/util/NonMPL2.h"
35 #ifndef EIGEN_CONSTRAINEDCG_H
36 #define EIGEN_CONSTRAINEDCG_H
50 template <
typename CMatrix,
typename CINVMatrix>
54 typedef typename CMatrix::Scalar Scalar;
55 typedef typename CMatrix::Index Index;
59 Index rows = C.rows(), cols = C.cols();
61 TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows);
62 Scalar rho, rho_1, alpha;
66 for (Index i = 0; i < rows; ++i)
77 l = C.transpose() * p;
79 alpha = rho / p.dot(q);
84 p = (rho/rho_1) * p + r;
87 l = C.transpose() * e;
89 for (Index j=0; j<l.size(); ++j)
91 CINV.fill(i,j) = l[j];
105 template<
typename TMatrix,
typename CMatrix,
106 typename VectorX,
typename VectorB,
typename VectorF>
110 typedef typename TMatrix::Scalar Scalar;
111 typedef typename TMatrix::Index Index;
114 Scalar rho = 1.0, rho_1, lambda, gamma;
115 Index xSize = x.size();
116 TmpVec p(xSize), q(xSize), q2(xSize),
117 r(xSize), old_z(xSize), z(xSize),
119 std::vector<bool> satured(C.rows());
121 iter.setRhsNorm(sqrt(b.dot(b)));
122 if (iter.rhsNorm() == 0.0) iter.setRhsNorm(1.0);
135 bool transition =
false;
136 for (Index i = 0; i < C.rows(); ++i)
138 Scalar al = C.row(i).dot(x) - f.coeff(i);
146 Scalar bb = CINV.row(i).dot(z);
149 for (
typename CMatrix::InnerIterator it(C,i); it; ++it)
150 z.coeffRef(it.index()) -= bb*it.value();
160 if (iter.finished(rho))
break;
162 if (iter.noiseLevel() > 0 && transition) std::cerr <<
"CCG: transition\n";
163 if (transition || iter.first()) gamma = 0.0;
164 else gamma = (std::max)(0.0, (rho - old_z.dot(z)) / rho_1);
170 lambda = rho / q.dot(p);
171 for (Index i = 0; i < C.rows(); ++i)
175 Scalar bb = C.row(i).dot(p) - f[i];
177 lambda = (std::min)(lambda, (f.coeff(i)-C.row(i).dot(x)) / bb);
189 #endif // EIGEN_CONSTRAINEDCG_H