unsupported/Eigen/src/IterativeSolvers/Scaling.h
Go to the documentation of this file.
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 
44 namespace Eigen {
45 using std::abs;
46 template<typename _MatrixType>
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 
64 
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  }
141  void computeRef (MatrixType& mat)
142  {
143  compute (mat);
144  mat = m_matrix;
145  }
148  VectorXd& LeftScaling()
149  {
150  return m_left;
151  }
152 
155  VectorXd& RightScaling()
156  {
157  return m_right;
158  }
159 
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;
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
IntermediateState sqrt(const Expression &arg)
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
#define eigen_assert(x)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:04