blob: 9400e81750575f2379c76ecfd37e08f09a05a5d5 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012 Desire NUENTSA WAKAM A; * // fill A and b; * IterScaling > scal; * // Compute the left and right scaling vectors. The matrix is equilibrated at output * scal.computeRef(A); * // Scale the right hand side * b = scal.LeftScaling().cwiseProduct(b); * // Now, solve the equilibrated linear system with any available solver * * // Scale back the computed solution * x = scal.RightScaling().cwiseProduct(x); * \endcode * * \tparam MatrixType_ the type of the matrix. It should be a real square sparsematrix * * References : D. Ruiz and B. Ucar, A Symmetry Preserving Algorithm for Matrix Scaling, INRIA Research report RR-7552 * * \sa \ref IncompleteLUT */ template class IterScaling { public: typedef MatrixType_ MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; public: IterScaling() { init(); } IterScaling(const MatrixType& matrix) { init(); compute(matrix); } ~IterScaling() { } /** * Compute the left and right diagonal matrices to scale the input matrix @p mat * * FIXME This algorithm will be modified such that the diagonal elements are permuted on the diagonal. * * \sa LeftScaling() RightScaling() */ void compute (const MatrixType& mat) { using std::abs; int m = mat.rows(); int n = mat.cols(); eigen_assert((m>0 && m == n) && "Please give a non - empty matrix"); m_left.resize(m); m_right.resize(n); m_left.setOnes(); m_right.setOnes(); m_matrix = mat; VectorXd Dr, Dc, DrRes, DcRes; // Temporary Left and right scaling vectors Dr.resize(m); Dc.resize(n); DrRes.resize(m); DcRes.resize(n); double EpsRow = 1.0, EpsCol = 1.0; int its = 0; do { // Iterate until the infinite norm of each row and column is approximately 1 // Get the maximum value in each row and column Dr.setZero(); Dc.setZero(); for (int k=0; km_tol || EpsCol > m_tol) && (its < m_maxits) ); m_isInitialized = true; } /** Compute the left and right vectors to scale the vectors * the input matrix is scaled with the computed vectors at output * * \sa compute() */ void computeRef (MatrixType& mat) { compute (mat); mat = m_matrix; } /** Get the vector to scale the rows of the matrix */ VectorXd& LeftScaling() { return m_left; } /** Get the vector to scale the columns of the matrix */ VectorXd& RightScaling() { return m_right; } /** Set the tolerance for the convergence of the iterative scaling algorithm */ void setTolerance(double tol) { m_tol = tol; } protected: void init() { m_tol = 1e-10; m_maxits = 5; m_isInitialized = false; } MatrixType m_matrix; mutable ComputationInfo m_info; bool m_isInitialized; VectorXd m_left; // Left scaling vector VectorXd m_right; // m_right scaling vector double m_tol; int m_maxits; // Maximum number of iterations allowed }; } #endif