Algorithm to equilibrate rows and columns of a square matrix

1 files changed, 200 insertions, 0 deletions

diff --git a/unsupported/Eigen/src/IterativeSolvers/Scaling.h b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
new file mode 100644
index 000000000..4aad69d0e
--- /dev/null
+++ b/unsupported/Eigen/src/IterativeSolvers/Scaling.h
@@ -0,0 +1,200 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.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_SCALING_H
+#define EIGEN_SCALING_H
+/**
+  * \ingroup IterativeSolvers_Module
+  * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
+  * 
+  * This class can be used as a preprocessing tool to accelerate the convergence of iterative methods 
+  * 
+  * This feature is  useful to limit the pivoting amount during LU/ILU factorization
+  * The  scaling strategy as presented here preserves the symmetry of the problem
+  * NOTE It is assumed that the matrix does not have empty row or column, 
+  * 
+  * Example with key steps 
+  * \code
+  * VectorXd x(n), b(n);
+  * SparseMatrix<double> A;
+  * // fill A and b;
+  * Scaling<SparseMatrix<double> > scal; 
+  * // Compute the left and right scaling vectors. The matrix is equilibrated at output
+  * scal.computeRef(A); 
+  * // Scale the right hand side
+  * b = scal.LeftScaling().cwiseProduct(b); 
+  * // Now, solve the equilibrated linear system with any available solver
+  * 
+  * // Scale back the computed solution
+  * x = scal.RightScaling().cwiseProduct(x); 
+  * \endcode
+  * 
+  * \tparam _MatrixType the type of the matrix. It should be a real square sparsematrix
+  * 
+  * References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552
+  * 
+  * \sa \ref IncompleteLUT 
+  */
+using std::abs; 
+using namespace Eigen;
+template<typename _MatrixType>
+class Scaling
+{
+  public:
+    typedef _MatrixType MatrixType; 
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::Index Index;
+    
+  public:
+    Scaling() { init(); }
+    
+    Scaling(const MatrixType& matrix)
+    {
+      init();
+      compute(matrix);
+    }
+    
+    ~Scaling() { }
+    
+    /** 
+     * Compute the left and right diagonal matrices to scale the input matrix @p mat
+     * 
+     * FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal. 
+     * 
+     * \sa LeftScaling() RightScaling()
+     */
+    void compute (const MatrixType& mat)
+    {
+      int m = mat.rows(); 
+      int n = mat.cols();
+      assert((m>0 && m == n) && "Please give a non - empty matrix");
+      m_left.resize(m); 
+      m_right.resize(n);
+      m_left.setOnes();
+      m_right.setOnes();
+      m_matrix = mat;
+      VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
+      Dr.resize(m); Dc.resize(n);
+      DrRes.resize(m); DcRes.resize(n);
+      double EpsRow = 1.0, EpsCol = 1.0;
+      int its = 0; 
+      do
+      { // Iterate until the infinite norm of each row and column is approximately 1
+        // Get the maximum value in each row and column
+        Dr.setZero(); Dc.setZero();
+        for (int k=0; k<m_matrix.outerSize(); ++k)
+        {
+          for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
+          {
+            if ( Dr(it.row()) < abs(it.value()) )
+              Dr(it.row()) = abs(it.value());
+            
+            if ( Dc(it.col()) < abs(it.value()) )
+              Dc(it.col()) = abs(it.value());
+          }
+        }
+        for (int i = 0; i < m; ++i) 
+        {
+          Dr(i) = std::sqrt(Dr(i));
+          Dc(i) = std::sqrt(Dc(i));
+        }
+        // Save the scaling factors 
+        for (int i = 0; i < m; ++i) 
+        {
+          m_left(i) /= Dr(i);
+          m_right(i) /= Dc(i);
+        }
+        // Scale the rows and the columns of the matrix
+        DrRes.setZero(); DcRes.setZero(); 
+        for (int k=0; k<m_matrix.outerSize(); ++k)
+        {
+          for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
+          {
+            it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
+            // Accumulate the norms of the row and column vectors   
+            if ( DrRes(it.row()) < abs(it.value()) )
+              DrRes(it.row()) = abs(it.value());
+            
+            if ( DcRes(it.col()) < abs(it.value()) )
+              DcRes(it.col()) = abs(it.value());
+          }
+        }  
+        DrRes.array() = (1-DrRes.array()).abs();
+        EpsRow = DrRes.maxCoeff();
+        DcRes.array() = (1-DcRes.array()).abs();
+        EpsCol = DcRes.maxCoeff();
+        its++;
+      }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
+      m_isInitialized = true;
+    }
+    /** Compute the left and right vectors to scale the vectors
+     * the input matrix is scaled with the computed vectors at output
+     * 
+     * \sa compute()
+     */
+    void computeRef (MatrixType& mat)
+    {
+      compute (mat);
+      mat = m_matrix;
+    }
+    /** Get the vector to scale the rows of the matrix 
+     */
+    VectorXd& LeftScaling()
+    {
+      return m_left;
+    }
+    
+    /** Get the vector to scale the columns of the matrix 
+     */
+    VectorXd& RightScaling()
+    {
+      return m_right;
+    }
+    
+    /** Set the tolerance for the convergence of the iterative scaling algorithm
+     */
+    void setTolerance(double tol)
+    {
+      m_tol = tol; 
+    }
+      
+  protected:
+    
+    void init()
+    {
+      m_tol = 1e-10;
+      m_maxits = 5;
+      m_isInitialized = false;
+    }
+    
+    MatrixType m_matrix;
+    mutable ComputationInfo m_info; 
+    bool m_isInitialized; 
+    VectorXd m_left; // Left scaling vector
+    VectorXd m_right; // m_right scaling vector
+    double m_tol; 
+    int m_maxits; // Maximum number of iterations allowed
+};
+
+#endif
\ No newline at end of file