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:
|
---|
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
|
---|