added an experimental IterativeSolvers module (currently in unsupported)

with a constrained conjugate gradient algorithm adapted from GMM++/ITL.
This algorithm is needed for Step.

4 files changed, 408 insertions, 1 deletions

diff --git a/unsupported/Eigen/CMakeLists.txt b/unsupported/Eigen/CMakeLists.txt
index 08b2a38a6..309ed3ee3 100644
--- a/unsupported/Eigen/CMakeLists.txt
+++ b/unsupported/Eigen/CMakeLists.txt
@@ -1,4 +1,4 @@
-set(Eigen_HEADERS AdolcForward)
+set(Eigen_HEADERS AdolcForward IterativeSolvers)
 
 install(FILES
   ${Eigen_HEADERS}
diff --git a/unsupported/Eigen/IterativeSolvers b/unsupported/Eigen/IterativeSolvers
new file mode 100644
index 000000000..e60ce4d1e
--- /dev/null
+++ b/unsupported/Eigen/IterativeSolvers
@@ -0,0 +1,49 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_ITERATIVE_SOLVERS_MODULE_H
+#define EIGEN_ITERATIVE_SOLVERS_MODULE_H
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+  * \defgroup IterativeSolvers_Module Iterative solvers module
+  * This module aims to provide various iterative linear and non linear solver algorithms.
+  * It currently provides:
+  *  - a constrained conjugate gradient
+  *
+  * \code
+  * #include <unsupported/Eigen/IterativeSolvers>
+  * \endcode
+  */
+//@{
+  
+#include "src/IterativeSolvers/IterationController.h"
+#include "src/IterativeSolvers/ConstrainedConjGrad.h"
+
+//@}
+
+}
+
+#endif // EIGEN_ITERATIVE_SOLVERS_MODULE_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
new file mode 100644
index 000000000..c00dcfddf
--- /dev/null
+++ b/unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h
@@ -0,0 +1,192 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/* NOTE The functions of this file have been adapted from the GMM++ library */
+
+//========================================================================
+//
+// Copyright (C) 2002-2007 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++ is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation; version 2.1 of the License.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+// You should have received a copy of the GNU Lesser General Public
+// License along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301,
+// USA.
+//
+//========================================================================
+
+#ifndef EIGEN_CONSTRAINEDCG_H
+#define EIGEN_CONSTRAINEDCG_H
+
+#include <Eigen/Core>
+
+/** \ingroup IterativeSolvers_Module
+  * Compute the pseudo inverse of the non-square matrix C such that
+  * \f$ CINV = (C * C^T)^{-1} * C \f$ based on a conjugate gradient method.
+  *
+  * This function is internally used by ei_constrained_cg.
+  */
+template <typename CMatrix, typename CINVMatrix>
+void ei_pseudo_inverse(const CMatrix &C, CINVMatrix &CINV)
+{
+  // optimisable : copie de la ligne, precalcul de C * trans(C).
+  typedef typename CMatrix::Scalar Scalar;
+  // FIXME use sparse vectors ?
+  typedef Matrix<Scalar,Dynamic,1> TmpVec;
+
+  int rows = C.rows(), cols = C.cols();
+
+  TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows);
+  Scalar rho, rho_1, alpha;
+  d.setZero();
+
+  CINV.startFill(); // FIXME estimate the number of non-zeros
+  for (int i = 0; i < rows; ++i)
+  {
+    d[i] = 1.0;
+    rho = 1.0;
+    e.setZero();
+    r = d;
+    p = d;
+
+    while (rho >= 1e-38)
+    { /* conjugate gradient to compute e             */
+      /* which is the i-th row of inv(C * trans(C))  */
+      l = C.transpose() * p;
+      q = C * l;
+      alpha = rho / p.dot(q);
+      e +=  alpha * p;
+      r += -alpha * q;
+      rho_1 = rho;
+      rho = r.dot(r);
+      p = (rho/rho_1) * p + r;
+    }
+
+    l = C.transpose() * e; // l is the i-th row of CINV
+    // FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse
+    for (int j=0; j<l.size(); ++j)
+      if (l[j]<1e-15)
+        CINV.fill(i,j) = l[j];
+
+    d[i] = 0.0;
+  }
+  CINV.endFill();
+}
+
+
+
+/** \ingroup IterativeSolvers_Module
+  * Constrained conjugate gradient
+  *
+  * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx <= f @\$
+  */
+template<typename TMatrix, typename CMatrix,
+         typename VectorX, typename VectorB, typename VectorF>
+void ei_constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x,
+                       const VectorB& b, const VectorF& f, IterationController &iter)
+{
+  typedef typename TMatrix::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1>  TmpVec;
+
+  Scalar rho = 1.0, rho_1, lambda, gamma;
+  int xSize = x.size();
+  TmpVec  p(xSize), q(xSize), q2(xSize),
+          r(xSize), old_z(xSize), z(xSize),
+          memox(xSize);
+  std::vector<bool> satured(C.rows());
+  p.setZero();
+  iter.setRhsNorm(ei_sqrt(b.dot(b))); // gael vect_sp(PS, b, b)
+  if (iter.rhsNorm() == 0.0) iter.setRhsNorm(1.0);
+
+  SparseMatrix<Scalar,RowMajor> CINV(C.rows(), C.cols());
+  ei_pseudo_inverse(C, CINV);
+
+  while(true)
+  {
+    // computation of residual
+    old_z = z;
+    memox = x;
+    r = b;
+    r += A * -x;
+    z = r;
+    bool transition = false;
+    for (int i = 0; i < C.rows(); ++i)
+    {
+      Scalar al = C.row(i).dot(x) - f.coeff(i);
+      if (al >= -1.0E-15)
+      {
+        if (!satured[i])
+        {
+          satured[i] = true;
+          transition = true;
+        }
+        Scalar bb = CINV.row(i).dot(z);
+        if (bb > 0.0)
+          // FIXME: we should allow that: z += -bb * C.row(i);
+          for (typename CMatrix::InnerIterator it(C,i); it; ++it)
+            z.coeffRef(it.index()) -= bb*it.value();
+      }
+      else
+        satured[i] = false;
+    }
+
+    // descent direction
+    rho_1 = rho;
+    rho = r.dot(z);
+
+    if (iter.finished(rho)) break;
+
+    if (iter.noiseLevel() > 0 && transition) std::cerr << "CCG: transition\n";
+    if (transition || iter.first()) gamma = 0.0;
+    else gamma = std::max(0.0, (rho - old_z.dot(z)) / rho_1);
+    p = z + gamma*p;
+
+    ++iter;
+    // one dimensionnal optimization
+    q = A * p;
+    lambda = rho / q.dot(p);
+    for (int i = 0; i < C.rows(); ++i)
+    {
+      if (!satured[i])
+      {
+        Scalar bb = C.row(i).dot(p) - f[i];
+        if (bb > 0.0)
+          lambda = std::min(lambda, (f.coeff(i)-C.row(i).dot(x)) / bb);
+      }
+    }
+    x += lambda * p;
+    memox -= x;
+  }
+}
+
+#endif // EIGEN_CONSTRAINEDCG_H
diff --git a/unsupported/Eigen/src/IterativeSolvers/IterationController.h b/unsupported/Eigen/src/IterativeSolvers/IterationController.h
new file mode 100644
index 000000000..a0ed3efde
--- /dev/null
+++ b/unsupported/Eigen/src/IterativeSolvers/IterationController.h
@@ -0,0 +1,166 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+/* NOTE The class IterationController has been adapted from the iteration
+ *      class of the GMM++ and ITL libraries.
+ */
+
+//=======================================================================
+// Copyright (C) 1997-2001
+// Authors: Andrew Lumsdaine <lums@osl.iu.edu> 
+//          Lie-Quan Lee     <llee@osl.iu.edu>
+//
+// This file is part of the Iterative Template Library
+//
+// You should have received a copy of the License Agreement for the
+// Iterative Template Library along with the software;  see the
+// file LICENSE.  
+//
+// Permission to modify the code and to distribute modified code is
+// granted, provided the text of this NOTICE is retained, a notice that
+// the code was modified is included with the above COPYRIGHT NOTICE and
+// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
+// file is distributed with the modified code.
+//
+// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
+// By way of example, but not limitation, Licensor MAKES NO
+// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
+// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
+// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
+// OR OTHER RIGHTS.
+//=======================================================================
+
+//========================================================================
+//
+// Copyright (C) 2002-2007 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++ is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Lesser General Public License as
+// published by the Free Software Foundation; version 2.1 of the License.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Lesser General Public License for more details.
+// You should have received a copy of the GNU Lesser General Public
+// License along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301,
+// USA.
+//
+//========================================================================
+
+#ifndef EIGEN_ITERATION_CONTROLLER_H
+#define EIGEN_ITERATION_CONTROLLER_H
+
+/** \ingroup IterativeSolvers_Module
+  * \class IterationController
+  *
+  * \brief Controls the iterations of the iterative solvers
+  *
+  * This class has been adapted from the iteration class of GMM++ and ITL libraries.
+  *
+  */
+class IterationController
+{
+  protected :
+    double m_rhsn;        ///< Right hand side norm
+    size_t m_maxiter;     ///< Max. number of iterations
+    int m_noise;          ///< if noise > 0 iterations are printed
+    double m_resmax;      ///< maximum residual
+    double m_resminreach, m_resadd;
+    size_t m_nit;         ///< iteration number
+    double m_res;         ///< last computed residual
+    bool m_written;
+    void (*m_callback)(const IterationController&);
+  public :
+
+    void init()
+    {
+      m_nit = 0; m_res = 0.0; m_written = false;
+      m_resminreach = 1E50; m_resadd = 0.0;
+      m_callback = 0;
+    }
+
+    IterationController(double r = 1.0E-8, int noi = 0, size_t mit = size_t(-1))
+      : m_rhsn(1.0), m_maxiter(mit), m_noise(noi), m_resmax(r) { init(); }
+
+    void operator ++(int) { m_nit++; m_written = false; m_resadd += m_res; }
+    void operator ++() { (*this)++; }
+
+    bool first() { return m_nit == 0; }
+
+    /* get/set the "noisyness" (verbosity) of the solvers */
+    int noiseLevel() const { return m_noise; }
+    void setNoiseLevel(int n) { m_noise = n; }
+    void reduceNoiseLevel() { if (m_noise > 0) m_noise--; }
+
+    double maxResidual() const { return m_resmax; }
+    void setMaxResidual(double r) { m_resmax = r; }
+
+    double residual() const { return m_res; }
+
+    /* change the user-definable callback, called after each iteration */
+    void setCallback(void (*t)(const IterationController&))
+    {
+      m_callback = t;
+    }
+
+    size_t iteration() const { return m_nit; }
+    void setIteration(size_t i) { m_nit = i; }
+
+    size_t maxIterarions() const { return m_maxiter; }
+    void setMaxIterations(size_t i) { m_maxiter = i; }
+
+    double rhsNorm() const { return m_rhsn; }
+    void setRhsNorm(double r) { m_rhsn = r; }
+
+    bool converged() const { return m_res <= m_rhsn * m_resmax; }
+    bool converged(double nr)
+    {
+      m_res = ei_abs(nr); 
+      m_resminreach = std::min(m_resminreach, m_res);
+      return converged();
+    }
+    template<typename VectorType> bool converged(const VectorType &v)
+    { return converged(v.squaredNorm()); }
+
+    bool finished(double nr)
+    {
+      if (m_callback) m_callback(*this);
+      if (m_noise > 0 && !m_written)
+      {
+        converged(nr);
+        m_written = true;
+      }
+      return (m_nit >= m_maxiter || converged(nr));
+    }
+    template <typename VectorType>
+    bool finished(const MatrixBase<VectorType> &v)
+    { return finished(double(v.squaredNorm())); }
+
+};
+
+#endif // EIGEN_ITERATION_CONTROLLER_H