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