cart-elc

Scaling.h (5853B)

---

 // 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
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #ifndef EIGEN_ITERSCALING_H
 #define EIGEN_ITERSCALING_H
 
 namespace Eigen {
 
 /**
   * \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;
   * IterScaling<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 
   */
 template<typename _MatrixType>
 class IterScaling
 {
   public:
     typedef _MatrixType MatrixType; 
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::Index Index;
     
   public:
     IterScaling() { init(); }
     
     IterScaling(const MatrixType& matrix)
     {
       init();
       compute(matrix);
     }
     
     ~IterScaling() { }
     
     /** 
      * 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)
     {
       using std::abs;
       int m = mat.rows(); 
       int n = mat.cols();
       eigen_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));
         }
         for (int i = 0; i < n; ++i) 
         {
           Dc(i) = std::sqrt(Dc(i));
         }
         // Save the scaling factors 
         for (int i = 0; i < m; ++i) 
         {
           m_left(i) /= Dr(i);
         }
         for (int i = 0; i < n; ++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