TbbTsqr.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_TbbTsqr_hpp
#define __TSQR_TbbTsqr_hpp

#include <TbbTsqr_TbbParallelTsqr.hpp>
#include <Tsqr_TimeStats.hpp>
#include <Teuchos_Time.hpp>
// #include <TbbRecursiveTsqr.hpp>

#include <stdexcept>
#include <string>
#include <utility> // std::pair
#include <vector>


namespace TSQR {
  namespace TBB {

    template< class LocalOrdinal, class Scalar, class TimerType = Teuchos::Time >
    class TbbTsqr {
    private:
      // Note: this is NOT a use of the pImpl idiom.  TbbRecursiveTsqr
      // is a nonparallel implementation that emulates the control
      // flow of the parallel implementation TbbParallelTsqr.  The
      // latter depends on the Intel Threading Building Blocks
      // library.
      //
      //TbbRecursiveTsqr< LocalOrdinal, Scalar > impl_;
      TbbParallelTsqr< LocalOrdinal, Scalar, TimerType > impl_;

      // Collected running statistcs on various computations
      mutable TimeStats factorStats_, applyStats_, explicitQStats_, cacheBlockStats_, unCacheBlockStats_;

      // Timers for various computations
      mutable TimerType factorTimer_, applyTimer_, explicitQTimer_, cacheBlockTimer_, unCacheBlockTimer_;

    public:
      typedef Scalar scalar_type;
      typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
      typedef LocalOrdinal ordinal_type;
      // typedef typename TbbRecursiveTsqr< LocalOrdinal, Scalar >::FactorOutput FactorOutput;
      typedef typename TbbParallelTsqr< LocalOrdinal, Scalar, TimerType >::FactorOutput FactorOutput;

      size_t ncores() const { return impl_.ncores(); }

      size_t cache_block_size() const { return impl_.cache_block_size(); }

      TbbTsqr (const size_t numCores,
         const size_t cacheBlockSize = 0) :
  impl_ (numCores, cacheBlockSize),
  factorTimer_ ("TbbTsqr::factor"),
  applyTimer_ ("TbbTsqr::apply"),
  explicitQTimer_ ("TbbTsqr::explicit_Q"),
  cacheBlockTimer_ ("TbbTsqr::cache_block"),
  unCacheBlockTimer_ ("TbbTsqr::un_cache_block")
      {}

      static bool QR_produces_R_factor_with_nonnegative_diagonal() {
  typedef TbbParallelTsqr< LocalOrdinal, Scalar, TimerType > impl_type;
  return impl_type::QR_produces_R_factor_with_nonnegative_diagonal();
      }

      void
      cache_block (const LocalOrdinal nrows,
       const LocalOrdinal ncols, 
       Scalar A_out[],
       const Scalar A_in[],
       const LocalOrdinal lda_in) const
      {
  cacheBlockTimer_.start(true);
  impl_.cache_block (nrows, ncols, A_out, A_in, lda_in);
  cacheBlockStats_.update (cacheBlockTimer_.stop());
      }

      void
      un_cache_block (const LocalOrdinal nrows,
          const LocalOrdinal ncols,
          Scalar A_out[],
          const LocalOrdinal lda_out,       
          const Scalar A_in[]) const
      {
  unCacheBlockTimer_.start(true);
  impl_.un_cache_block (nrows, ncols, A_out, lda_out, A_in);
  unCacheBlockStats_.update (unCacheBlockTimer_.stop());
      }

      void
      fill_with_zeros (const LocalOrdinal nrows,
           const LocalOrdinal ncols,
           Scalar C[],
           const LocalOrdinal ldc, 
           const bool contiguous_cache_blocks = false) const
      {
  impl_.fill_with_zeros (nrows, ncols, C, ldc, contiguous_cache_blocks);
      }

      template< class MatrixViewType >
      MatrixViewType
      top_block (const MatrixViewType& C, 
     const bool contiguous_cache_blocks = false) const
      {
  return impl_.top_block (C, contiguous_cache_blocks);
      }

      FactorOutput
      factor (const LocalOrdinal nrows,
        const LocalOrdinal ncols, 
        Scalar A[],
        const LocalOrdinal lda,
        Scalar R[],
        const LocalOrdinal ldr,
        const bool contiguous_cache_blocks = false)
      {
  factorTimer_.start(true);
  return impl_.factor (nrows, ncols, A, lda, R, ldr, contiguous_cache_blocks);
  factorStats_.update (factorTimer_.stop());
      }

      void
      apply (const ApplyType& apply_type,
       const LocalOrdinal nrows,
       const LocalOrdinal ncols_Q,
       const Scalar Q[],
       const LocalOrdinal ldq,
       const FactorOutput& factor_output,
       const LocalOrdinal ncols_C,
       Scalar C[],
       const LocalOrdinal ldc,
       const bool contiguous_cache_blocks = false)
      {
  applyTimer_.start(true);
  impl_.apply (apply_type, nrows, ncols_Q, Q, ldq, factor_output, 
         ncols_C, C, ldc, contiguous_cache_blocks);
  applyStats_.update (applyTimer_.stop());
      }

      void
      explicit_Q (const LocalOrdinal nrows,
      const LocalOrdinal ncols_Q_in,
      const Scalar Q_in[],
      const LocalOrdinal ldq_in,
      const FactorOutput& factor_output,
      const LocalOrdinal ncols_Q_out,
      Scalar Q_out[],
      const LocalOrdinal ldq_out,
      const bool contiguous_cache_blocks = false)
      {
  explicitQTimer_.start(true);
  impl_.explicit_Q (nrows, ncols_Q_in, Q_in, ldq_in, factor_output,
        ncols_Q_out, Q_out, ldq_out, contiguous_cache_blocks);
  explicitQStats_.update (explicitQTimer_.stop());
      }

      void
      Q_times_B (const LocalOrdinal nrows,
     const LocalOrdinal ncols,
     Scalar Q[],
     const LocalOrdinal ldq,
     const Scalar B[],
     const LocalOrdinal ldb,
     const bool contiguous_cache_blocks = false) const
      {
  impl_.Q_times_B (nrows, ncols, Q, ldq, B, ldb, contiguous_cache_blocks);
      }

      LocalOrdinal
      reveal_R_rank (const LocalOrdinal ncols,
         Scalar R[],
         const LocalOrdinal ldr,
         Scalar U[],
         const LocalOrdinal ldu,
         const magnitude_type tol) const 
      {
  return impl_.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 contiguous_cache_blocks = false)
      {
  return impl_.reveal_rank (nrows, ncols, Q, ldq, R, ldr, tol, 
          contiguous_cache_blocks);
      }

      double
      min_seq_factor_timing () const { return impl_.min_seq_factor_timing(); }
      double
      max_seq_factor_timing () const { return impl_.max_seq_factor_timing(); }
      double
      min_seq_apply_timing () const { return impl_.min_seq_apply_timing(); }
      double
      max_seq_apply_timing () const { return impl_.max_seq_apply_timing(); }

      void getStats (std::vector< TimeStats >& stats) {
  const int numStats = 5;
  stats.resize (numStats);
  stats[0] = factorStats_;
  stats[1] = applyStats_;
  stats[2] = explicitQStats_;
  stats[3] = cacheBlockStats_;
  stats[4] = unCacheBlockStats_;
      }

      void getStatsLabels (std::vector< std::string >& labels) {
  const int numStats = 5;
  labels.resize (numStats);
  labels[0] = factorTimer_.name();
  labels[1] = applyTimer_.name();
  labels[2] = explicitQTimer_.name();
  labels[3] = cacheBlockTimer_.name();
  labels[4] = unCacheBlockTimer_.name();
      }

    }; // class TbbTsqr

  } // namespace TBB
} // namespace TSQR

#endif // __TSQR_TbbTsqr_hpp