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  */
44namespace Eigen {
45using std::abs;
46template<typename _MatrixType>
47class 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