Scaling.h
Go to the documentation of this file.
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 
00044 namespace Eigen {
00045 using std::abs; 
00046 template<typename _MatrixType>
00047 class IterScaling
00048 {
00049   public:
00050     typedef _MatrixType MatrixType; 
00051     typedef typename MatrixType::Scalar Scalar;
00052     typedef typename MatrixType::Index Index;
00053     
00054   public:
00055     IterScaling() { init(); }
00056     
00057     IterScaling(const MatrixType& matrix)
00058     {
00059       init();
00060       compute(matrix);
00061     }
00062     
00063     ~IterScaling() { }
00064     
00072     void compute (const MatrixType& mat)
00073     {
00074       int m = mat.rows(); 
00075       int n = mat.cols();
00076       eigen_assert((m>0 && m == n) && "Please give a non - empty matrix");
00077       m_left.resize(m); 
00078       m_right.resize(n);
00079       m_left.setOnes();
00080       m_right.setOnes();
00081       m_matrix = mat;
00082       VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors
00083       Dr.resize(m); Dc.resize(n);
00084       DrRes.resize(m); DcRes.resize(n);
00085       double EpsRow = 1.0, EpsCol = 1.0;
00086       int its = 0; 
00087       do
00088       { // Iterate until the infinite norm of each row and column is approximately 1
00089         // Get the maximum value in each row and column
00090         Dr.setZero(); Dc.setZero();
00091         for (int k=0; k<m_matrix.outerSize(); ++k)
00092         {
00093           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00094           {
00095             if ( Dr(it.row()) < abs(it.value()) )
00096               Dr(it.row()) = abs(it.value());
00097             
00098             if ( Dc(it.col()) < abs(it.value()) )
00099               Dc(it.col()) = abs(it.value());
00100           }
00101         }
00102         for (int i = 0; i < m; ++i) 
00103         {
00104           Dr(i) = std::sqrt(Dr(i));
00105           Dc(i) = std::sqrt(Dc(i));
00106         }
00107         // Save the scaling factors 
00108         for (int i = 0; i < m; ++i) 
00109         {
00110           m_left(i) /= Dr(i);
00111           m_right(i) /= Dc(i);
00112         }
00113         // Scale the rows and the columns of the matrix
00114         DrRes.setZero(); DcRes.setZero(); 
00115         for (int k=0; k<m_matrix.outerSize(); ++k)
00116         {
00117           for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
00118           {
00119             it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
00120             // Accumulate the norms of the row and column vectors   
00121             if ( DrRes(it.row()) < abs(it.value()) )
00122               DrRes(it.row()) = abs(it.value());
00123             
00124             if ( DcRes(it.col()) < abs(it.value()) )
00125               DcRes(it.col()) = abs(it.value());
00126           }
00127         }  
00128         DrRes.array() = (1-DrRes.array()).abs();
00129         EpsRow = DrRes.maxCoeff();
00130         DcRes.array() = (1-DcRes.array()).abs();
00131         EpsCol = DcRes.maxCoeff();
00132         its++;
00133       }while ( (EpsRow >m_tol || EpsCol > m_tol) && (its < m_maxits) );
00134       m_isInitialized = true;
00135     }
00141     void computeRef (MatrixType& mat)
00142     {
00143       compute (mat);
00144       mat = m_matrix;
00145     }
00148     VectorXd& LeftScaling()
00149     {
00150       return m_left;
00151     }
00152     
00155     VectorXd& RightScaling()
00156     {
00157       return m_right;
00158     }
00159     
00162     void setTolerance(double tol)
00163     {
00164       m_tol = tol; 
00165     }
00166       
00167   protected:
00168     
00169     void init()
00170     {
00171       m_tol = 1e-10;
00172       m_maxits = 5;
00173       m_isInitialized = false;
00174     }
00175     
00176     MatrixType m_matrix;
00177     mutable ComputationInfo m_info; 
00178     bool m_isInitialized; 
00179     VectorXd m_left; // Left scaling vector
00180     VectorXd m_right; // m_right scaling vector
00181     double m_tol; 
00182     int m_maxits; // Maximum number of iterations allowed
00183 };
00184 }
00185 #endif


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:38:51