$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 00006 /* NOTE The functions of this file have been adapted from the GMM++ library */ 00007 00008 //======================================================================== 00009 // 00010 // Copyright (C) 2002-2007 Yves Renard 00011 // 00012 // This file is a part of GETFEM++ 00013 // 00014 // Getfem++ is free software; you can redistribute it and/or modify 00015 // it under the terms of the GNU Lesser General Public License as 00016 // published by the Free Software Foundation; version 2.1 of the License. 00017 // 00018 // This program is distributed in the hope that it will be useful, 00019 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00021 // GNU Lesser General Public License for more details. 00022 // You should have received a copy of the GNU Lesser General Public 00023 // License along with this program; if not, write to the Free Software 00024 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, 00025 // USA. 00026 // 00027 //======================================================================== 00028 00029 #include "../../../../Eigen/src/Core/util/NonMPL2.h" 00030 00031 #ifndef EIGEN_CONSTRAINEDCG_H 00032 #define EIGEN_CONSTRAINEDCG_H 00033 00034 #include <Eigen/Core> 00035 00036 namespace Eigen { 00037 00038 namespace internal { 00039 00046 template <typename CMatrix, typename CINVMatrix> 00047 void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV) 00048 { 00049 // optimisable : copie de la ligne, precalcul de C * trans(C). 00050 typedef typename CMatrix::Scalar Scalar; 00051 typedef typename CMatrix::Index Index; 00052 // FIXME use sparse vectors ? 00053 typedef Matrix<Scalar,Dynamic,1> TmpVec; 00054 00055 Index rows = C.rows(), cols = C.cols(); 00056 00057 TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows); 00058 Scalar rho, rho_1, alpha; 00059 d.setZero(); 00060 00061 typedef Triplet<double> T; 00062 std::vector<T> tripletList; 00063 00064 for (Index i = 0; i < rows; ++i) 00065 { 00066 d[i] = 1.0; 00067 rho = 1.0; 00068 e.setZero(); 00069 r = d; 00070 p = d; 00071 00072 while (rho >= 1e-38) 00073 { /* conjugate gradient to compute e */ 00074 /* which is the i-th row of inv(C * trans(C)) */ 00075 l = C.transpose() * p; 00076 q = C * l; 00077 alpha = rho / p.dot(q); 00078 e += alpha * p; 00079 r += -alpha * q; 00080 rho_1 = rho; 00081 rho = r.dot(r); 00082 p = (rho/rho_1) * p + r; 00083 } 00084 00085 l = C.transpose() * e; // l is the i-th row of CINV 00086 // FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse 00087 for (Index j=0; j<l.size(); ++j) 00088 if (l[j]<1e-15) 00089 tripletList.push_back(T(i,j,l(j))); 00090 00091 00092 d[i] = 0.0; 00093 } 00094 CINV.setFromTriplets(tripletList.begin(), tripletList.end()); 00095 } 00096 00097 00098 00104 template<typename TMatrix, typename CMatrix, 00105 typename VectorX, typename VectorB, typename VectorF> 00106 void constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x, 00107 const VectorB& b, const VectorF& f, IterationController &iter) 00108 { 00109 using std::sqrt; 00110 typedef typename TMatrix::Scalar Scalar; 00111 typedef typename TMatrix::Index Index; 00112 typedef Matrix<Scalar,Dynamic,1> TmpVec; 00113 00114 Scalar rho = 1.0, rho_1, lambda, gamma; 00115 Index xSize = x.size(); 00116 TmpVec p(xSize), q(xSize), q2(xSize), 00117 r(xSize), old_z(xSize), z(xSize), 00118 memox(xSize); 00119 std::vector<bool> satured(C.rows()); 00120 p.setZero(); 00121 iter.setRhsNorm(sqrt(b.dot(b))); // gael vect_sp(PS, b, b) 00122 if (iter.rhsNorm() == 0.0) iter.setRhsNorm(1.0); 00123 00124 SparseMatrix<Scalar,RowMajor> CINV(C.rows(), C.cols()); 00125 pseudo_inverse(C, CINV); 00126 00127 while(true) 00128 { 00129 // computation of residual 00130 old_z = z; 00131 memox = x; 00132 r = b; 00133 r += A * -x; 00134 z = r; 00135 bool transition = false; 00136 for (Index i = 0; i < C.rows(); ++i) 00137 { 00138 Scalar al = C.row(i).dot(x) - f.coeff(i); 00139 if (al >= -1.0E-15) 00140 { 00141 if (!satured[i]) 00142 { 00143 satured[i] = true; 00144 transition = true; 00145 } 00146 Scalar bb = CINV.row(i).dot(z); 00147 if (bb > 0.0) 00148 // FIXME: we should allow that: z += -bb * C.row(i); 00149 for (typename CMatrix::InnerIterator it(C,i); it; ++it) 00150 z.coeffRef(it.index()) -= bb*it.value(); 00151 } 00152 else 00153 satured[i] = false; 00154 } 00155 00156 // descent direction 00157 rho_1 = rho; 00158 rho = r.dot(z); 00159 00160 if (iter.finished(rho)) break; 00161 00162 if (iter.noiseLevel() > 0 && transition) std::cerr << "CCG: transition\n"; 00163 if (transition || iter.first()) gamma = 0.0; 00164 else gamma = (std::max)(0.0, (rho - old_z.dot(z)) / rho_1); 00165 p = z + gamma*p; 00166 00167 ++iter; 00168 // one dimensionnal optimization 00169 q = A * p; 00170 lambda = rho / q.dot(p); 00171 for (Index i = 0; i < C.rows(); ++i) 00172 { 00173 if (!satured[i]) 00174 { 00175 Scalar bb = C.row(i).dot(p) - f[i]; 00176 if (bb > 0.0) 00177 lambda = (std::min)(lambda, (f.coeff(i)-C.row(i).dot(x)) / bb); 00178 } 00179 } 00180 x += lambda * p; 00181 memox -= x; 00182 } 00183 } 00184 00185 } // end namespace internal 00186 00187 } // end namespace Eigen 00188 00189 #endif // EIGEN_CONSTRAINEDCG_H