48 #ifndef EIGEN_CONSTRAINEDCG_H
49 #define EIGEN_CONSTRAINEDCG_H
63 template <
typename CMatrix,
typename CINVMatrix>
67 typedef typename CMatrix::Scalar Scalar;
68 typedef typename CMatrix::Index Index;
70 typedef Matrix<Scalar,Dynamic,1> TmpVec;
72 Index rows = C.rows(), cols = C.cols();
74 TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows);
75 Scalar rho, rho_1, alpha;
79 for (Index i = 0; i < rows; ++i)
90 l = C.transpose() * p;
92 alpha = rho / p.dot(q);
97 p = (rho/rho_1) * p + r;
100 l = C.transpose() * e;
102 for (Index j=0; j<l.size(); ++j)
104 CINV.fill(i,j) = l[j];
118 template<
typename TMatrix,
typename CMatrix,
119 typename VectorX,
typename VectorB,
typename VectorF>
123 typedef typename TMatrix::Scalar Scalar;
124 typedef typename TMatrix::Index Index;
125 typedef Matrix<Scalar,Dynamic,1> TmpVec;
127 Scalar rho = 1.0, rho_1, lambda, gamma;
128 Index xSize = x.size();
129 TmpVec p(xSize), q(xSize), q2(xSize),
130 r(xSize), old_z(xSize), z(xSize),
132 std::vector<bool> satured(C.rows());
134 iter.setRhsNorm(
sqrt(b.dot(b)));
135 if (iter.rhsNorm() == 0.0) iter.setRhsNorm(1.0);
137 SparseMatrix<Scalar,RowMajor> CINV(C.rows(), C.cols());
148 bool transition =
false;
149 for (Index i = 0; i < C.rows(); ++i)
151 Scalar al = C.row(i).dot(x) - f.coeff(i);
159 Scalar bb = CINV.row(i).dot(z);
162 for (
typename CMatrix::InnerIterator it(C,i); it; ++it)
163 z.coeffRef(it.index()) -= bb*it.value();
173 if (iter.finished(rho))
break;
175 if (iter.noiseLevel() > 0 && transition) std::cerr <<
"CCG: transition\n";
176 if (transition || iter.first()) gamma = 0.0;
177 else gamma = (std::max)(0.0, (rho - old_z.dot(z)) / rho_1);
183 lambda = rho / q.dot(p);
184 for (Index i = 0; i < C.rows(); ++i)
188 Scalar bb = C.row(i).dot(p) - f[i];
190 lambda = (std::min)(lambda, (f.coeff(i)-C.row(i).dot(x)) / bb);
202 #endif // EIGEN_CONSTRAINEDCG_H