1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam (at) inria.fr 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_ITERSCALING_H 11 #define EIGEN_ITERSCALING_H 12 /** 13 * \ingroup IterativeSolvers_Module 14 * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices 15 * 16 * This class can be used as a preprocessing tool to accelerate the convergence of iterative methods 17 * 18 * This feature is useful to limit the pivoting amount during LU/ILU factorization 19 * The scaling strategy as presented here preserves the symmetry of the problem 20 * NOTE It is assumed that the matrix does not have empty row or column, 21 * 22 * Example with key steps 23 * \code 24 * VectorXd x(n), b(n); 25 * SparseMatrix<double> A; 26 * // fill A and b; 27 * IterScaling<SparseMatrix<double> > scal; 28 * // Compute the left and right scaling vectors. The matrix is equilibrated at output 29 * scal.computeRef(A); 30 * // Scale the right hand side 31 * b = scal.LeftScaling().cwiseProduct(b); 32 * // Now, solve the equilibrated linear system with any available solver 33 * 34 * // Scale back the computed solution 35 * x = scal.RightScaling().cwiseProduct(x); 36 * \endcode 37 * 38 * \tparam _MatrixType the type of the matrix. It should be a real square sparsematrix 39 * 40 * References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552 41 * 42 * \sa \ref IncompleteLUT 43 */ 44 namespace Eigen { 45 using std::abs; 46 template<typename _MatrixType> 47 class IterScaling 48 { 49 public: 50 typedef _MatrixType MatrixType; 51 typedef typename MatrixType::Scalar Scalar; 52 typedef typename MatrixType::Index Index; 53 54 public: 55 IterScaling() { init(); } 56 57 IterScaling(const MatrixType& matrix) 58 { 59 init(); 60 compute(matrix); 61 } 62 63 ~IterScaling() { } 64 65 /** 66 * Compute the left and right diagonal matrices to scale the input matrix @p mat 67 * 68 * FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal. 69 * 70 * \sa LeftScaling() RightScaling() 71 */ 72 void compute (const MatrixType& mat) 73 { 74 int m = mat.rows(); 75 int n = mat.cols(); 76 eigen_assert((m>0 && m == n) && "Please give a non - empty matrix"); 77 m_left.resize(m); 78 m_right.resize(n); 79 m_left.setOnes(); 80 m_right.setOnes(); 81 m_matrix = mat; 82 VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors 83 Dr.resize(m); Dc.resize(n); 84 DrRes.resize(m); DcRes.resize(n); 85 double EpsRow = 1.0, EpsCol = 1.0; 86 int its = 0; 87 do 88 { // Iterate until the infinite norm of each row and column is approximately 1 89 // Get the maximum value in each row and column 90 Dr.setZero(); Dc.setZero(); 91 for (int k=0; k<m_matrix.outerSize(); ++k) 92 { 93 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) 94 { 95 if ( Dr(it.row()) < abs(it.value()) ) 96 Dr(it.row()) = abs(it.value()); 97 98 if ( Dc(it.col()) < abs(it.value()) ) 99 Dc(it.col()) = abs(it.value()); 100 } 101 } 102 for (int i = 0; i < m; ++i) 103 { 104 Dr(i) = std::sqrt(Dr(i)); 105 Dc(i) = std::sqrt(Dc(i)); 106 } 107 // Save the scaling factors 108 for (int i = 0; i < m; ++i) 109 { 110 m_left(i) /= Dr(i); 111 m_right(i) /= Dc(i); 112 } 113 // Scale the rows and the columns of the matrix 114 DrRes.setZero(); DcRes.setZero(); 115 for (int k=0; k<m_matrix.outerSize(); ++k) 116 { 117 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) 118 { 119 it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) ); 120 // Accumulate the norms of the row and column vectors 121 if ( DrRes(it.row()) < abs(it.value()) ) 122 DrRes(it.row()) = abs(it.value()); 123 124 if ( DcRes(it.col()) < abs(it.value()) ) 125 DcRes(it.col()) = abs(it.value()); 126 } 127 } 128 DrRes.array() = (1-DrRes.array()).abs(); 129 EpsRow = DrRes.maxCoeff(); 130 DcRes.array() = (1-DcRes.array()).abs(); 131 EpsCol = DcRes.maxCoeff(); 132 its++; 133 }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) ); 134 m_isInitialized = true; 135 } 136 /** Compute the left and right vectors to scale the vectors 137 * the input matrix is scaled with the computed vectors at output 138 * 139 * \sa compute() 140 */ 141 void computeRef (MatrixType& mat) 142 { 143 compute (mat); 144 mat = m_matrix; 145 } 146 /** Get the vector to scale the rows of the matrix 147 */ 148 VectorXd& LeftScaling() 149 { 150 return m_left; 151 } 152 153 /** Get the vector to scale the columns of the matrix 154 */ 155 VectorXd& RightScaling() 156 { 157 return m_right; 158 } 159 160 /** Set the tolerance for the convergence of the iterative scaling algorithm 161 */ 162 void setTolerance(double tol) 163 { 164 m_tol = tol; 165 } 166 167 protected: 168 169 void init() 170 { 171 m_tol = 1e-10; 172 m_maxits = 5; 173 m_isInitialized = false; 174 } 175 176 MatrixType m_matrix; 177 mutable ComputationInfo m_info; 178 bool m_isInitialized; 179 VectorXd m_left; // Left scaling vector 180 VectorXd m_right; // m_right scaling vector 181 double m_tol; 182 int m_maxits; // Maximum number of iterations allowed 183 }; 184 } 185 #endif 186
