Quelle Scaling.h Sprache: C

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Desire NUENTSA WAKAM <desire.nuentsa_wakam@inria.fr
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_ITERSCALING_H
#define EIGEN_ITERSCALING_H

namespace Eigen {

/**
  * \ingroup IterativeSolvers_Module
  * \brief iterative scaling algorithm to equilibrate rows and column norms in matrices
  *
  * This class can be used as a preprocessing tool to accelerate the convergence of iterative methods
  *
  * This feature is  useful to limit the pivoting amount during LU/ILU factorization
  * The  scaling strategy as presented here preserves the symmetry of the problem
  * NOTE It is assumed that the matrix does not have empty row or column,
  *
  * Example with key steps
  * \code
  * VectorXd x(n), b(n);
  * SparseMatrix<double> A;
  * // fill A and b;
  * IterScaling<SparseMatrix<double> > 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<typename _MatrixType>
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; k<m_matrix.outerSize(); ++k)
        {
          for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
          {
            if ( Dr(it.row()) < abs(it.value()) )
              Dr(it.row()) = abs(it.value());

            if ( Dc(it.col()) < abs(it.value()) )
              Dc(it.col()) = abs(it.value());
          }
        }
        for (int i = 0; i < m; ++i)
        {
          Dr(i) = std::sqrt(Dr(i));
        }
        for (int i = 0; i < n; ++i)
        {
          Dc(i) = std::sqrt(Dc(i));
        }
        // Save the scaling factors
        for (int i = 0; i < m; ++i)
        {
          m_left(i) /= Dr(i);
        }
        for (int i = 0; i < n; ++i)
        {
          m_right(i) /= Dc(i);
        }
        // Scale the rows and the columns of the matrix
        DrRes.setZero(); DcRes.setZero();
        for (int k=0; k<m_matrix.outerSize(); ++k)
        {
          for (typename MatrixType::InnerIterator it(m_matrix, k); it; ++it)
          {
            it.valueRef() = it.value()/( Dr(it.row()) * Dc(it.col()) );
            // Accumulate the norms of the row and column vectors
            if ( DrRes(it.row()) < abs(it.value()) )
              DrRes(it.row()) = abs(it.value());

            if ( DcRes(it.col()) < abs(it.value()) )
              DcRes(it.col()) = abs(it.value());
          }
        }
        DrRes.array() = (1-DrRes.array()).abs();
        EpsRow = DrRes.maxCoeff();
        DcRes.array() = (1-DcRes.array()).abs();
        EpsCol = DcRes.maxCoeff();
        its++;
      }while ( (EpsRow >m_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

quality63%

¤ Dauer der Verarbeitung: 0.14 Sekunden (vorverarbeitet) ¤

Wurzel

Suchen

Beweissystem der NASA

Beweissystem Isabelle

NIST Cobol Testsuite

Cephes Mathematical Library

Wiener Entwicklungsmethode

Haftungshinweis

Die Informationen auf dieser Webseite wurden nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit, noch Qualität der bereit gestellten Informationen zugesichert.

Bemerkung:

Die farbliche Syntaxdarstellung ist noch experimentell.