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:
56     IterScaling() { init(); }
57 
58     IterScaling(const MatrixType& matrix)
59     {
60       init();
61       compute(matrix);
62     }
63 
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      */
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           Dc(i) = std::sqrt(Dc(i));
108         }
109         // Save the scaling factors
110         for (int i = 0; i < m; ++i)
111         {
112           m_left(i) /= Dr(i);
113           m_right(i) /= Dc(i);
114         }
115         // Scale the rows and the columns of the matrix
116         DrRes.setZero(); DcRes.setZero();
117         for (int k=0; k<m_matrix.outerSize(); ++k)
118         {
119           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
120           {
121             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
122             // Accumulate the norms of the row and column vectors
123             if ( DrRes(it.row()) < abs(it.value()) )
124               DrRes(it.row()) = abs(it.value());
125 
126             if ( DcRes(it.col()) < abs(it.value()) )
127               DcRes(it.col()) = abs(it.value());
128           }
129         }
130         DrRes.array() = (1-DrRes.array()).abs();
131         EpsRow = DrRes.maxCoeff();
132         DcRes.array() = (1-DcRes.array()).abs();
133         EpsCol = DcRes.maxCoeff();
134         its++;
135       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
136       m_isInitialized = true;
137     }
138     /** Compute the left and right vectors to scale the vectors
139      * the input matrix is scaled with the computed vectors at output
140      *
141      * \sa compute()
142      */
143     void computeRef (MatrixType& mat)
144     {
145       compute (mat);
146       mat = m_matrix;
147     }
148     /** Get the vector to scale the rows of the matrix
149      */
150     VectorXd& LeftScaling()
151     {
152       return m_left;
153     }
154 
155     /** Get the vector to scale the columns of the matrix
156      */
157     VectorXd& RightScaling()
158     {
159       return m_right;
160     }
161 
162     /** Set the tolerance for the convergence of the iterative scaling algorithm
163      */
164     void setTolerance(double tol)
165     {
166       m_tol = tol;
167     }
168 
169   protected:
170 
171     void init()
172     {
173       m_tol = 1e-10;
174       m_maxits = 5;
175       m_isInitialized = false;
176     }
177 
178     MatrixType m_matrix;
179     mutable ComputationInfo m_info;
180     bool m_isInitialized;
181     VectorXd m_left; // Left scaling vector
182     VectorXd m_right; // m_right scaling vector
183     double m_tol;
184     int m_maxits; // Maximum number of iterations allowed
185 };
186 }
187 #endif
188