SSL UpperBidiagonalization.h Interaktion und
PortierbarkeitC

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2013-2014 Gael Guennebaud <gael.guennebaud@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_BIDIAGONALIZATION_H
#define EIGEN_BIDIAGONALIZATION_H

namespace Eigen {

namespace internal {
// UpperBidiagonalization will probably be replaced by a Bidiagonalization class, don't want to make it stable API.
// At the same time, it's useful to keep for now as it's about the only thing that is testing the BandMatrix class.

template<typename _MatrixType> class UpperBidiagonalization
{
  public:

    typedef _MatrixType MatrixType;
    enum {
      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
      ColsAtCompileTimeMinusOne = internal::decrement_size<ColsAtCompileTime>::ret
    };
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
    typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
    typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
    typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0, RowMajor> BidiagonalType;
    typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType;
    typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType;
    typedef HouseholderSequence<
              const MatrixType,
              const typename internal::remove_all<typename Diagonal<const MatrixType,0>::ConjugateReturnType>::type
            > HouseholderUSequenceType;
    typedef HouseholderSequence<
              const typename internal::remove_all<typename MatrixType::ConjugateReturnType>::type,
              Diagonal<const MatrixType,1>,
              OnTheRight
            > HouseholderVSequenceType;

    /**
    * \brief Default Constructor.
    *
    * The default constructor is useful in cases in which the user intends to
    * perform decompositions via Bidiagonalization::compute(const MatrixType&).
    */
    UpperBidiagonalization() : m_householder(), m_bidiagonal(), m_isInitialized(false) {}

    explicit UpperBidiagonalization(const MatrixType& matrix)
      : m_householder(matrix.rows(), matrix.cols()),
        m_bidiagonal(matrix.cols(), matrix.cols()),
        m_isInitialized(false)
    {
      compute(matrix);
    }

    UpperBidiagonalization& compute(const MatrixType& matrix);
    UpperBidiagonalization& computeUnblocked(const MatrixType& matrix);

    const MatrixType& householder() const { return m_householder; }
    const BidiagonalType& bidiagonal() const { return m_bidiagonal; }

    const HouseholderUSequenceType householderU() const
    {
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
      return HouseholderUSequenceType(m_householder, m_householder.diagonal().conjugate());
    }

    const HouseholderVSequenceType householderV() // const here gives nasty errors and i'm lazy
    {
      eigen_assert(m_isInitialized && "UpperBidiagonalization is not initialized.");
      return HouseholderVSequenceType(m_householder.conjugate(), m_householder.const_derived().template diagonal<1>())
             .setLength(m_householder.cols()-1)
             .setShift(1);
    }

  protected:
    MatrixType m_householder;
    BidiagonalType m_bidiagonal;
    bool m_isInitialized;
};

// Standard upper bidiagonalization without fancy optimizations
// This version should be faster for small matrix size
template<typename MatrixType>
void upperbidiagonalization_inplace_unblocked(MatrixType& mat,
                                              typename MatrixType::RealScalar *diagonal,
                                              typename MatrixType::RealScalar *upper_diagonal,
                                              typename MatrixType::Scalar* tempData = 0)
{
  typedef typename MatrixType::Scalar Scalar;

  Index rows = mat.rows();
  Index cols = mat.cols();

  typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxRowsAtCompileTime,1> TempType;
  TempType tempVector;
  if(tempData==0)
  {
    tempVector.resize(rows);
    tempData = tempVector.data();
  }

  for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k)
  {
    Index remainingRows = rows - k;
    Index remainingCols = cols - k - 1;

    // construct left householder transform in-place in A
    mat.col(k).tail(remainingRows)
       .makeHouseholderInPlace(mat.coeffRef(k,k), diagonal[k]);
    // apply householder transform to remaining part of A on the left
    mat.bottomRightCorner(remainingRows, remainingCols)
       .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), mat.coeff(k,k), tempData);

    if(k == cols-1) break;

    // construct right householder transform in-place in mat
    mat.row(k).tail(remainingCols)
       .makeHouseholderInPlace(mat.coeffRef(k,k+1), upper_diagonal[k]);
    // apply householder transform to remaining part of mat on the left
    mat.bottomRightCorner(remainingRows-1, remainingCols)
       .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).adjoint(), mat.coeff(k,k+1), tempData);
  }
}

/** \internal
  * Helper routine for the block reduction to upper bidiagonal form.
  *
  * Let's partition the matrix A:
  *
  *      | A00 A01 |
  *  A = |         |
  *      | A10 A11 |
  *
  * This function reduces to bidiagonal form the left \c rows x \a blockSize vertical panel [A00/A10]
  * and the \a blockSize x \c cols horizontal panel [A00 A01] of the matrix \a A. The bottom-right block A11
  * is updated using matrix-matrix products:
  *   A22 -= V * Y^T - X * U^T
  * where V and U contains the left and right Householder vectors. U and V are stored in A10, and A01
  * respectively, and the update matrices X and Y are computed during the reduction.
  *
  */
template<typename MatrixType>
void upperbidiagonalization_blocked_helper(MatrixType& A,
                                           typename MatrixType::RealScalar *diagonal,
                                           typename MatrixType::RealScalar *upper_diagonal,
                                           Index bs,
                                           Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
                                                      traits<MatrixType>::Flags & RowMajorBit> > X,
                                           Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic,
                                                      traits<MatrixType>::Flags & RowMajorBit> > Y)
{
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  typedef typename NumTraits<RealScalar>::Literal Literal;
  enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
  typedef InnerStride<int(StorageOrder) == int(ColMajor) ? 1 : Dynamic> ColInnerStride;
  typedef InnerStride<int(StorageOrder) == int(ColMajor) ? Dynamic : 1> RowInnerStride;
  typedef Ref<Matrix<Scalar, Dynamic, 1>, 0, ColInnerStride>    SubColumnType;
  typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, RowInnerStride>    SubRowType;
  typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder > > SubMatType;

  Index brows = A.rows();
  Index bcols = A.cols();

  Scalar tau_u, tau_u_prev(0), tau_v;

  for(Index k = 0; k < bs; ++k)
  {
    Index remainingRows = brows - k;
    Index remainingCols = bcols - k - 1;

    SubMatType X_k1( X.block(k,0, remainingRows,k) );
    SubMatType V_k1( A.block(k,0, remainingRows,k) );

    // 1 - update the k-th column of A
    SubColumnType v_k = A.col(k).tail(remainingRows);
          v_k -= V_k1 * Y.row(k).head(k).adjoint();
    if(k) v_k -= X_k1 * A.col(k).head(k);

    // 2 - construct left Householder transform in-place
    v_k.makeHouseholderInPlace(tau_v, diagonal[k]);

    if(k+1<bcols)
    {
      SubMatType Y_k  ( Y.block(k+1,0, remainingCols, k+1) );
      SubMatType U_k1 ( A.block(0,k+1, k,remainingCols) );

      // this eases the application of Householder transforAions
      // A(k,k) will store tau_v later
      A(k,k) = Scalar(1);

      // 3 - Compute y_k^T = tau_v * ( A^T*v_k - Y_k-1*V_k-1^T*v_k - U_k-1*X_k-1^T*v_k )
      {
        SubColumnType y_k( Y.col(k).tail(remainingCols) );

        // let's use the beginning of column k of Y as a temporary vector
        SubColumnType tmp( Y.col(k).head(k) );
        y_k.noalias()  = A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k; // bottleneck
        tmp.noalias()  = V_k1.adjoint()  * v_k;
        y_k.noalias() -= Y_k.leftCols(k) * tmp;
        tmp.noalias()  = X_k1.adjoint()  * v_k;
        y_k.noalias() -= U_k1.adjoint()  * tmp;
        y_k *= numext::conj(tau_v);
      }

      // 4 - update k-th row of A (it will become u_k)
      SubRowType u_k( A.row(k).tail(remainingCols) );
      u_k = u_k.conjugate();
      {
        u_k -= Y_k * A.row(k).head(k+1).adjoint();
        if(k) u_k -= U_k1.adjoint() * X.row(k).head(k).adjoint();
      }

      // 5 - construct right Householder transform in-place
      u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]);

      // this eases the application of Householder transformations
      // A(k,k+1) will store tau_u later
      A(k,k+1) = Scalar(1);

      // 6 - Compute x_k = tau_u * ( A*u_k - X_k-1*U_k-1^T*u_k - V_k*Y_k^T*u_k )
      {
        SubColumnType x_k ( X.col(k).tail(remainingRows-1) );

        // let's use the beginning of column k of X as a temporary vectors
        // note that tmp0 and tmp1 overlaps
        SubColumnType tmp0 ( X.col(k).head(k) ),
                      tmp1 ( X.col(k).head(k+1) );

        x_k.noalias()   = A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose(); // bottleneck
        tmp0.noalias()  = U_k1 * u_k.transpose();
        x_k.noalias()  -= X_k1.bottomRows(remainingRows-1) * tmp0;
        tmp1.noalias()  = Y_k.adjoint() * u_k.transpose();
        x_k.noalias()  -= A.block(k+1,0, remainingRows-1,k+1) * tmp1;
        x_k *= numext::conj(tau_u);
        tau_u = numext::conj(tau_u);
        u_k = u_k.conjugate();
      }

      if(k>0) A.coeffRef(k-1,k) = tau_u_prev;
      tau_u_prev = tau_u;
    }
    else
      A.coeffRef(k-1,k) = tau_u_prev;

    A.coeffRef(k,k) = tau_v;
  }

  if(bs<bcols)
    A.coeffRef(bs-1,bs) = tau_u_prev;

  // update A22
  if(bcols>bs && brows>bs)
  {
    SubMatType A11( A.bottomRightCorner(brows-bs,bcols-bs) );
    SubMatType A10( A.block(bs,0, brows-bs,bs) );
    SubMatType A01( A.block(0,bs, bs,bcols-bs) );
    Scalar tmp = A01(bs-1,0);
    A01(bs-1,0) = Literal(1);
    A11.noalias() -= A10 * Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint();
    A11.noalias() -= X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A01;
    A01(bs-1,0) = tmp;
  }
}

/** \internal
  *
  * Implementation of a block-bidiagonal reduction.
  * It is based on the following paper:
  *   The Design of a Parallel Dense Linear Algebra Software Library: Reduction to Hessenberg, Tridiagonal, and Bidiagonal Form.
  *   by Jaeyoung Choi, Jack J. Dongarra, David W. Walker. (1995)
  *   section 3.3
  */
template<typename MatrixType, typename BidiagType>
void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagonal,
                                            Index maxBlockSize=32,
                                            typename MatrixType::Scalar* /*tempData*/ = 0)
{
  typedef typename MatrixType::Scalar Scalar;
  typedef Block<MatrixType,Dynamic,Dynamic> BlockType;

  Index rows = A.rows();
  Index cols = A.cols();
  Index size = (std::min)(rows, cols);

  // X and Y are work space
  enum { StorageOrder = traits<MatrixType>::Flags & RowMajorBit };
  Matrix<Scalar,
         MatrixType::RowsAtCompileTime,
         Dynamic,
         StorageOrder,
         MatrixType::MaxRowsAtCompileTime> X(rows,maxBlockSize);
  Matrix<Scalar,
         MatrixType::ColsAtCompileTime,
         Dynamic,
         StorageOrder,
         MatrixType::MaxColsAtCompileTime> Y(cols,maxBlockSize);
  Index blockSize = (std::min)(maxBlockSize,size);

  Index k = 0;
  for(k = 0; k < size; k += blockSize)
  {
    Index bs = (std::min)(size-k,blockSize);  // actual size of the block
    Index brows = rows - k;                   // rows of the block
    Index bcols = cols - k;                   // columns of the block

    // partition the matrix A:
    //
    //      | A00 A01 A02 |
    //      |             |
    // A  = | A10 A11 A12 |
    //      |             |
    //      | A20 A21 A22 |
    //
    // where A11 is a bs x bs diagonal block,
    // and let:
    //      | A11 A12 |
    //  B = |         |
    //      | A21 A22 |

    BlockType B = A.block(k,k,brows,bcols);

    // This stage performs the bidiagonalization of A11, A21, A12, and updating of A22.
    // Finally, the algorithm continue on the updated A22.
    //
    // However, if B is too small, or A22 empty, then let's use an unblocked strategy
    if(k+bs==cols || bcols<48) // somewhat arbitrary threshold
    {
      upperbidiagonalization_inplace_unblocked(B,
                                               &(bidiagonal.template diagonal<0>().coeffRef(k)),
                                               &(bidiagonal.template diagonal<1>().coeffRef(k)),
                                               X.data()
                                              );
      break; // We're done
    }
    else
    {
      upperbidiagonalization_blocked_helper<BlockType>( B,
                                                        &(bidiagonal.template diagonal<0>().coeffRef(k)),
                                                        &(bidiagonal.template diagonal<1>().coeffRef(k)),
                                                        bs,
                                                        X.topLeftCorner(brows,bs),
                                                        Y.topLeftCorner(bcols,bs)
                                                      );
    }
  }
}

template<typename _MatrixType>
UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::computeUnblocked(const _MatrixType& matrix)
{
  Index rows = matrix.rows();
  Index cols = matrix.cols();
  EIGEN_ONLY_USED_FOR_DEBUG(cols);

  eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols.");

  m_householder = matrix;

  ColVectorType temp(rows);

  upperbidiagonalization_inplace_unblocked(m_householder,
                                           &(m_bidiagonal.template diagonal<0>().coeffRef(0)),
                                           &(m_bidiagonal.template diagonal<1>().coeffRef(0)),
                                           temp.data());

  m_isInitialized = true;
  return *this;
}

template<typename _MatrixType>
UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix)
{
  Index rows = matrix.rows();
  Index cols = matrix.cols();
  EIGEN_ONLY_USED_FOR_DEBUG(rows);
  EIGEN_ONLY_USED_FOR_DEBUG(cols);

  eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols.");

  m_householder = matrix;
  upperbidiagonalization_inplace_blocked(m_householder, m_bidiagonal);

  m_isInitialized = true;
  return *this;
}

#if 0
/** \return the Householder QR decomposition of \c *this.
  *
  * \sa class Bidiagonalization
  */
template<typename Derived>
const UpperBidiagonalization<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::bidiagonalization() const
{
  return UpperBidiagonalization<PlainObject>(eval());
}
#endif

} // end namespace internal

} // end namespace Eigen

#endif // EIGEN_BIDIAGONALIZATION_H

quality97%

¤ Diese beiden folgenden Angebotsgruppen bietet das Unternehmen0.6Angebot ¤

*Eine klare Vorstellung vom Zielzustand

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.

SSL UpperBidiagonalization.h Interaktion und PortierbarkeitC

¤ Diese beiden folgenden Angebotsgruppen bietet das Unternehmen0.6Angebot ¤

Haftungshinweis

Bemerkung:

SSL UpperBidiagonalization.h Interaktion und
PortierbarkeitC