blob: 9b3eb53e06cc675f157b9a18991a383b947b21c9 [file] [log] [blame]
Brian Silverman72890c22015-09-19 14:37:37 -04001// 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
Austin Schuh189376f2018-12-20 22:11:15 +110012
13namespace Eigen {
14
Brian Silverman72890c22015-09-19 14:37:37 -040015/**
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 */
Brian Silverman72890c22015-09-19 14:37:37 -040047template<typename _MatrixType>
48class 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 {
Austin Schuh189376f2018-12-20 22:11:15 +110075 using std::abs;
Brian Silverman72890c22015-09-19 14:37:37 -040076 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));
Austin Schuhc55b0172022-02-20 17:52:35 -0800107 }
108 for (int i = 0; i < n; ++i)
109 {
Brian Silverman72890c22015-09-19 14:37:37 -0400110 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);
Austin Schuhc55b0172022-02-20 17:52:35 -0800116 }
117 for (int i = 0; i < n; ++i)
118 {
Brian Silverman72890c22015-09-19 14:37:37 -0400119 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 */
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 */
156 VectorXd& LeftScaling()
157 {
158 return m_left;
159 }
160
161 /** Get the vector to scale the columns of the matrix
162 */
163 VectorXd& RightScaling()
164 {
165 return m_right;
166 }
167
168 /** Set the tolerance for the convergence of the iterative scaling algorithm
169 */
170 void setTolerance(double tol)
171 {
172 m_tol = tol;
173 }
174
175 protected:
176
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