ConstrainedConjGrad.h
Go to the documentation of this file.
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 // This Source Code Form is subject to the terms of the Mozilla
00007 // Public License v. 2.0. If a copy of the MPL was not distributed
00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00009 
00010 /* NOTE The functions of this file have been adapted from the GMM++ library */
00011 
00012 //========================================================================
00013 //
00014 // Copyright (C) 2002-2007 Yves Renard
00015 //
00016 // This file is a part of GETFEM++
00017 //
00018 // Getfem++ is free software; you can redistribute it and/or modify
00019 // it under the terms of the GNU Lesser General Public License as
00020 // published by the Free Software Foundation; version 2.1 of the License.
00021 //
00022 // This program is distributed in the hope that it will be useful,
00023 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00024 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00025 // GNU Lesser General Public License for more details.
00026 // You should have received a copy of the GNU Lesser General Public
00027 // License along with this program; if not, write to the Free Software
00028 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301,
00029 // USA.
00030 //
00031 //========================================================================
00032 
00033 #include "../../../../Eigen/src/Core/util/NonMPL2.h"
00034 
00035 #ifndef EIGEN_CONSTRAINEDCG_H
00036 #define EIGEN_CONSTRAINEDCG_H
00037 
00038 #include <Eigen/Core>
00039 
00040 namespace Eigen { 
00041 
00042 namespace internal {
00043 
00050 template <typename CMatrix, typename CINVMatrix>
00051 void pseudo_inverse(const CMatrix &C, CINVMatrix &CINV)
00052 {
00053   // optimisable : copie de la ligne, precalcul de C * trans(C).
00054   typedef typename CMatrix::Scalar Scalar;
00055   typedef typename CMatrix::Index Index;
00056   // FIXME use sparse vectors ?
00057   typedef Matrix<Scalar,Dynamic,1> TmpVec;
00058 
00059   Index rows = C.rows(), cols = C.cols();
00060 
00061   TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows);
00062   Scalar rho, rho_1, alpha;
00063   d.setZero();
00064 
00065   CINV.startFill(); // FIXME estimate the number of non-zeros
00066   for (Index i = 0; i < rows; ++i)
00067   {
00068     d[i] = 1.0;
00069     rho = 1.0;
00070     e.setZero();
00071     r = d;
00072     p = d;
00073 
00074     while (rho >= 1e-38)
00075     { /* conjugate gradient to compute e             */
00076       /* which is the i-th row of inv(C * trans(C))  */
00077       l = C.transpose() * p;
00078       q = C * l;
00079       alpha = rho / p.dot(q);
00080       e +=  alpha * p;
00081       r += -alpha * q;
00082       rho_1 = rho;
00083       rho = r.dot(r);
00084       p = (rho/rho_1) * p + r;
00085     }
00086 
00087     l = C.transpose() * e; // l is the i-th row of CINV
00088     // FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse
00089     for (Index j=0; j<l.size(); ++j)
00090       if (l[j]<1e-15)
00091         CINV.fill(i,j) = l[j];
00092 
00093     d[i] = 0.0;
00094   }
00095   CINV.endFill();
00096 }
00097 
00098 
00099 
00105 template<typename TMatrix, typename CMatrix,
00106          typename VectorX, typename VectorB, typename VectorF>
00107 void constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x,
00108                        const VectorB& b, const VectorF& f, IterationController &iter)
00109 {
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


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:21