Tsqr.hpp

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright (2010) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef __TSQR_Tsqr_hpp
#define __TSQR_Tsqr_hpp

#include <Tsqr_ApplyType.hpp>
#include <Tsqr_Matrix.hpp>
#include <Tsqr_MessengerBase.hpp>
#include <Tsqr_DistTsqr.hpp>
#include <Tsqr_SequentialTsqr.hpp>
#include <Tsqr_ScalarTraits.hpp>
#include <Tsqr_Util.hpp>


namespace TSQR {

  template< class LocalOrdinal, 
      class Scalar, 
      class NodeTsqrType = SequentialTsqr< LocalOrdinal, Scalar >,
      class DistTsqrType = DistTsqr< LocalOrdinal, Scalar > >
  class Tsqr {
  public:
    typedef MatView< LocalOrdinal, Scalar > matview_type;
    typedef ConstMatView< LocalOrdinal, Scalar > const_matview_type;
    typedef Matrix< LocalOrdinal, Scalar > matrix_type;

    typedef Scalar scalar_type;
    typedef LocalOrdinal ordinal_type;
    typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;

    typedef NodeTsqrType node_tsqr_type;
    typedef DistTsqrType dist_tsqr_type;
    typedef typename Teuchos::RCP< node_tsqr_type > node_tsqr_ptr;
    typedef typename Teuchos::RCP< dist_tsqr_type > dist_tsqr_ptr;
    typedef typename dist_tsqr_type::rank_type rank_type;

    typedef typename node_tsqr_type::FactorOutput NodeOutput;
    typedef typename dist_tsqr_type::FactorOutput DistOutput;
    typedef std::pair< NodeOutput, DistOutput > FactorOutput;

    Tsqr (const node_tsqr_ptr& nodeTsqr, const dist_tsqr_ptr& distTsqr) :
      nodeTsqr_ (nodeTsqr), distTsqr_ (distTsqr)
    {}

    size_t cache_block_size() const { return nodeTsqr_->cache_block_size(); }

    bool QR_produces_R_factor_with_nonnegative_diagonal () const {
      return nodeTsqr_->QR_produces_R_factor_with_nonnegative_diagonal() &&
  distTsqr_->QR_produces_R_factor_with_nonnegative_diagonal();
    }


    void
    factorExplicit (matview_type& A, 
        matview_type& Q, 
        matview_type& R, 
        const bool contiguousCacheBlocks = false)
    {
      // Sanity checks for matrix dimensions
      if (A.nrows() < A.ncols())
  throw std::logic_error ("A must have at least as many rows as columns");
      else if (A.nrows() != Q.nrows())
  throw std::logic_error ("A and Q must have same number of rows");
      else if (R.nrows() < R.ncols())
  throw std::logic_error ("R must have at least as many rows as columns");
      else if (A.ncols() != R.ncols())
  throw std::logic_error ("A and R must have the same number of columns");

      // Checks for quick exit, based on matrix dimensions
      if (Q.ncols() == ordinal_type(0))
  return;

      R.fill (scalar_type (0));
      NodeOutput nodeResults = 
  nodeTsqr_->factor (A.nrows(), A.ncols(), A.get(), A.lda(), 
         R.get(), R.lda(), contiguousCacheBlocks);
      nodeTsqr_->fill_with_zeros (Q.nrows(), Q.ncols(), Q.get(), Q.lda(), 
          contiguousCacheBlocks);
      matview_type Q_top_block = 
  nodeTsqr_->top_block (Q, contiguousCacheBlocks);
      if (Q_top_block.nrows() < R.ncols())
  throw std::logic_error("Q has too few rows");
      matview_type Q_top (R.ncols(), Q.ncols(), 
        Q_top_block.get(), Q_top_block.lda());
      distTsqr_->factorExplicit (R, Q_top);
      nodeTsqr_->apply (ApplyType::NoTranspose, 
      A.nrows(), A.ncols(), A.get(), A.lda(), nodeResults,
      Q.ncols(), Q.get(), Q.lda(), contiguousCacheBlocks);
    }


    FactorOutput
    factor (const LocalOrdinal nrows_local,
      const LocalOrdinal ncols, 
      Scalar A_local[],
      const LocalOrdinal lda_local,
      Scalar R[],
      const LocalOrdinal ldr,
      const bool contiguousCacheBlocks = false)
    {
      MatView< LocalOrdinal, Scalar > R_view (ncols, ncols, R, ldr);
      R_view.fill (Scalar(0));
      NodeOutput nodeResults = 
  nodeTsqr_->factor (nrows_local, ncols, A_local, lda_local, 
        R_view.get(), R_view.lda(), 
        contiguousCacheBlocks);
      DistOutput distResults = distTsqr_->factor (R_view);
      return std::make_pair (nodeResults, distResults);
    }

    void
    apply (const std::string& op,
     const LocalOrdinal nrows_local,
     const LocalOrdinal ncols_Q,
     const Scalar Q_local[],
     const LocalOrdinal ldq_local,
     const FactorOutput& factor_output,
     const LocalOrdinal ncols_C,
     Scalar C_local[],
     const LocalOrdinal ldc_local,
     const bool contiguousCacheBlocks = false)
    {
      ApplyType applyType (op);

      // This determines the order in which we apply the intranode
      // part of the Q factor vs. the internode part of the Q factor.
      const bool transposed = applyType.transposed();

      // View of this node's local part of the matrix C.
      matview_type C_view (nrows_local, ncols_C, C_local, ldc_local);

      // Identify top ncols_C by ncols_C block of C_local.
      // top_block() will set C_top_view to have the correct leading
      // dimension, whether or not cache blocks are stored
      // contiguously.
      //
      // C_top_view is the topmost cache block of C_local.  It has at
      // least as many rows as columns, but it likely has more rows
      // than columns.
      matview_type C_view_top_block = 
  nodeTsqr_->top_block (C_view, contiguousCacheBlocks);

      // View of the topmost ncols_C by ncols_C block of C.
      matview_type C_top_view (ncols_C, ncols_C, C_view_top_block.get(), 
             C_view_top_block.lda());

      if (! transposed) 
  {
    // C_top (small compact storage) gets a deep copy of the top
    // ncols_C by ncols_C block of C_local.
          matrix_type C_top (C_top_view);

    // Compute in place on all processors' C_top blocks.
    distTsqr_->apply (applyType, C_top.ncols(), ncols_Q, C_top.get(), 
          C_top.lda(), factor_output.second);

    // Copy the result from C_top back into the top ncols_C by
    // ncols_C block of C_local.
    C_top_view.copy (C_top);

    // Apply the local Q factor (in Q_local and
    // factor_output.first) to C_local.
    nodeTsqr_->apply (applyType, nrows_local, ncols_Q, 
          Q_local, ldq_local, factor_output.first, 
          ncols_C, C_local, ldc_local,
          contiguousCacheBlocks);
  }
      else
  {
    // Apply the (transpose of the) local Q factor (in Q_local
    // and factor_output.first) to C_local.
    nodeTsqr_->apply (applyType, nrows_local, ncols_Q, 
          Q_local, ldq_local, factor_output.first, 
          ncols_C, C_local, ldc_local,
          contiguousCacheBlocks);

    // C_top (small compact storage) gets a deep copy of the top
    // ncols_C by ncols_C block of C_local.
          matrix_type C_top (C_top_view);

    // Compute in place on all processors' C_top blocks.
    distTsqr_->apply (applyType, ncols_C, ncols_Q, C_top.get(), 
          C_top.lda(), factor_output.second);

    // Copy the result from C_top back into the top ncols_C by
    // ncols_C block of C_local.
    C_top_view.copy (C_top);
  }
    }

    void
    explicit_Q (const LocalOrdinal nrows_local,
    const LocalOrdinal ncols_Q_in,
    const Scalar Q_local_in[],
    const LocalOrdinal ldq_local_in,
    const FactorOutput& factorOutput,
    const LocalOrdinal ncols_Q_out,
    Scalar Q_local_out[],
    const LocalOrdinal ldq_local_out,
    const bool contiguousCacheBlocks = false)
    {
      nodeTsqr_->fill_with_zeros (nrows_local, ncols_Q_out, Q_local_out,
         ldq_local_out, contiguousCacheBlocks);
      const rank_type myRank = distTsqr_->rank();
      if (myRank == 0)
  {
    // View of this node's local part of the matrix Q_out.
    matview_type Q_out_view (nrows_local, ncols_Q_out, 
           Q_local_out, ldq_local_out);

    // View of the topmost cache block of Q_out.  It is
    // guaranteed to have at least as many rows as columns.
    matview_type Q_out_top = 
      nodeTsqr_->top_block (Q_out_view, contiguousCacheBlocks);

    // Fill (topmost cache block of) Q_out with the first
    // ncols_Q_out columns of the identity matrix.
    for (ordinal_type j = 0; j < ncols_Q_out; ++j)
      Q_out_top(j, j) = Scalar(1);
  }
      apply ("N", nrows_local, 
       ncols_Q_in, Q_local_in, ldq_local_in, factorOutput,
       ncols_Q_out, Q_local_out, ldq_local_out,
       contiguousCacheBlocks);
    }

    void
    Q_times_B (const LocalOrdinal nrows,
         const LocalOrdinal ncols,
         Scalar Q[],
         const LocalOrdinal ldq,
         const Scalar B[],
         const LocalOrdinal ldb,
         const bool contiguousCacheBlocks = false) const
    {
      // This requires no internode communication.  However, the work
      // is not redundant, since each MPI process has a different Q.
      nodeTsqr_->Q_times_B (nrows, ncols, Q, ldq, B, ldb, 
         contiguousCacheBlocks);
      // We don't need a barrier after this method, unless users are
      // doing mean MPI_Get() things.
    }

    LocalOrdinal
    reveal_R_rank (const LocalOrdinal ncols,
       Scalar R[],
       const LocalOrdinal ldr,
       Scalar U[],
       const LocalOrdinal ldu,
       const magnitude_type tol) const 
    {
      // Forward the request to the intranode TSQR implementation.
      // This work is performed redundantly on all MPI processes.
      //
      // FIXME (mfh 26 Aug 2010) Could be a problem if your cluster is
      // heterogeneous, because then you might obtain different
      // answers (due to rounding error) on different processors.
      return nodeTsqr_->reveal_R_rank (ncols, R, ldr, U, ldu, tol);
    }

    LocalOrdinal
    reveal_rank (const LocalOrdinal nrows,
     const LocalOrdinal ncols,
     Scalar Q[],
     const LocalOrdinal ldq,
     Scalar R[],
     const LocalOrdinal ldr,
     const magnitude_type tol,
     const bool contiguousCacheBlocks = false) const
    {
      // Take the easy exit if available.
      if (ncols == 0)
  return 0;

      //
      // FIXME (mfh 16 Jul 2010) We _should_ compute the SVD of R (as
      // the copy B) on Proc 0 only.  This would ensure that all
      // processors get the same SVD and rank (esp. in a heterogeneous
      // computing environment).  For now, we just do this computation
      // redundantly, and hope that all the returned rank values are
      // the same.
      //
      matrix_type U (ncols, ncols, Scalar(0));
      const ordinal_type rank = 
  reveal_R_rank (ncols, R, ldr, U.get(), U.lda(), tol);
      
      if (rank < ncols)
  {
    // cerr << ">>> Rank of R: " << rank << " < ncols=" << ncols << endl;
    // cerr << ">>> Resulting U:" << endl;
    // print_local_matrix (cerr, ncols, ncols, R, ldr);
    // cerr << endl;

    // If R is not full rank: reveal_R_rank() already computed
    // the SVD \f$R = U \Sigma V^*\f$ of (the input) R, and
    // overwrote R with \f$\Sigma V^*\f$.  Now, we compute \f$Q
    // := Q \cdot U\f$, respecting cache blocks of Q.
    Q_times_B (nrows, ncols, Q, ldq, U.get(), U.lda(), 
         contiguousCacheBlocks);
  }
      return rank;
    }

    void
    cache_block (const LocalOrdinal nrows_local,
     const LocalOrdinal ncols,
     Scalar A_local_out[],
     const Scalar A_local_in[],
     const LocalOrdinal lda_local_in) const
    {
      nodeTsqr_->cache_block (nrows_local, ncols, 
            A_local_out, 
            A_local_in, lda_local_in);
    }

    void
    un_cache_block (const LocalOrdinal nrows_local,
        const LocalOrdinal ncols,
        Scalar A_local_out[],
        const LocalOrdinal lda_local_out,       
        const Scalar A_local_in[]) const
    {
      nodeTsqr_->un_cache_block (nrows_local, ncols, 
         A_local_out, lda_local_out, 
         A_local_in);
    }

  private:
    node_tsqr_ptr nodeTsqr_;
    dist_tsqr_ptr distTsqr_;
  }; // class Tsqr

} // namespace TSQR

#endif // __TSQR_Tsqr_hpp