![]() |
Eigen-unsupported
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_ITERSCALING_H 00011 #define EIGEN_ITERSCALING_H 00012 00013 namespace Eigen { 00014 00047 template<typename _MatrixType> 00048 class IterScaling 00049 { 00050 public: 00051 typedef _MatrixType MatrixType; 00052 typedef typename MatrixType::Scalar Scalar; 00053 typedef typename MatrixType::Index Index; 00054 00055 public: 00056 IterScaling() { init(); } 00057 00058 IterScaling(const MatrixType& matrix) 00059 { 00060 init(); 00061 compute(matrix); 00062 } 00063 00064 ~IterScaling() { } 00065 00073 void compute (const MatrixType& mat) 00074 { 00075 using std::abs; 00076 int m = mat.rows(); 00077 int n = mat.cols(); 00078 eigen_assert((m>0 && m == n) && "Please give a non - empty matrix"); 00079 m_left.resize(m); 00080 m_right.resize(n); 00081 m_left.setOnes(); 00082 m_right.setOnes(); 00083 m_matrix = mat; 00084 VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors 00085 Dr.resize(m); Dc.resize(n); 00086 DrRes.resize(m); DcRes.resize(n); 00087 double EpsRow = 1.0, EpsCol = 1.0; 00088 int its = 0; 00089 do 00090 { // Iterate until the infinite norm of each row and column is approximately 1 00091 // Get the maximum value in each row and column 00092 Dr.setZero(); Dc.setZero(); 00093 for (int k=0; k<m_matrix.outerSize(); ++k) 00094 { 00095 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) 00096 { 00097 if ( Dr(it.row()) < abs(it.value()) ) 00098 Dr(it.row()) = abs(it.value()); 00099 00100 if ( Dc(it.col()) < abs(it.value()) ) 00101 Dc(it.col()) = abs(it.value()); 00102 } 00103 } 00104 for (int i = 0; i < m; ++i) 00105 { 00106 Dr(i) = std::sqrt(Dr(i)); 00107 Dc(i) = std::sqrt(Dc(i)); 00108 } 00109 // Save the scaling factors 00110 for (int i = 0; i < m; ++i) 00111 { 00112 m_left(i) /= Dr(i); 00113 m_right(i) /= Dc(i); 00114 } 00115 // Scale the rows and the columns of the matrix 00116 DrRes.setZero(); DcRes.setZero(); 00117 for (int k=0; k<m_matrix.outerSize(); ++k) 00118 { 00119 for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it) 00120 { 00121 it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) ); 00122 // Accumulate the norms of the row and column vectors 00123 if ( DrRes(it.row()) < abs(it.value()) ) 00124 DrRes(it.row()) = abs(it.value()); 00125 00126 if ( DcRes(it.col()) < abs(it.value()) ) 00127 DcRes(it.col()) = abs(it.value()); 00128 } 00129 } 00130 DrRes.array() = (1-DrRes.array()).abs(); 00131 EpsRow = DrRes.maxCoeff(); 00132 DcRes.array() = (1-DcRes.array()).abs(); 00133 EpsCol = DcRes.maxCoeff(); 00134 its++; 00135 }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) ); 00136 m_isInitialized = true; 00137 } 00143 void computeRef (MatrixType& mat) 00144 { 00145 compute (mat); 00146 mat = m_matrix; 00147 } 00150 VectorXd& LeftScaling() 00151 { 00152 return m_left; 00153 } 00154 00157 VectorXd& RightScaling() 00158 { 00159 return m_right; 00160 } 00161 00164 void setTolerance(double tol) 00165 { 00166 m_tol = tol; 00167 } 00168 00169 protected: 00170 00171 void init() 00172 { 00173 m_tol = 1e-10; 00174 m_maxits = 5; 00175 m_isInitialized = false; 00176 } 00177 00178 MatrixType m_matrix; 00179 mutable ComputationInfo m_info; 00180 bool m_isInitialized; 00181 VectorXd m_left; // Left scaling vector 00182 VectorXd m_right; // m_right scaling vector 00183 double m_tol; 00184 int m_maxits; // Maximum number of iterations allowed 00185 }; 00186 } 00187 #endif