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   * \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:
IterScaling()55     IterScaling() { init(); }
56 
IterScaling(const MatrixType & matrix)57     IterScaling(const MatrixType& matrix)
58     {
59       init();
60       compute(matrix);
61     }
62 
~IterScaling()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      */
compute(const MatrixType & mat)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      */
computeRef(MatrixType & mat)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      */
LeftScaling()148     VectorXd& LeftScaling()
149     {
150       return m_left;
151     }
152 
153     /** Get the vector to scale the columns of the matrix
154      */
RightScaling()155     VectorXd& RightScaling()
156     {
157       return m_right;
158     }
159 
160     /** Set the tolerance for the convergence of the iterative scaling algorithm
161      */
setTolerance(double tol)162     void setTolerance(double tol)
163     {
164       m_tol = tol;
165     }
166 
167   protected:
168 
init()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