TbbTsqr_TbbRecursiveTsqr.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_TbbRecursiveTsqr_hpp
#define __TSQR_TbbRecursiveTsqr_hpp

#include <Tsqr_ApplyType.hpp>
#include <Tsqr_CacheBlocker.hpp>
#include <Tsqr_SequentialTsqr.hpp>
#include <TbbTsqr_Partitioner.hpp>

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


namespace TSQR {
  namespace TBB {

    template< class LocalOrdinal, class Scalar >
    class TbbRecursiveTsqr {
    public:
      TbbRecursiveTsqr (const size_t num_cores = 1,
      const size_t cache_block_size = 0);

      size_t ncores() const { return ncores_; }

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

      typedef typename SequentialTsqr< LocalOrdinal, Scalar >::FactorOutput SeqOutput;
      typedef std::vector< std::vector< Scalar > > ParOutput;
      typedef typename std::pair< std::vector< SeqOutput >, ParOutput > FactorOutput;

      void
      cache_block (const LocalOrdinal nrows,
       const LocalOrdinal ncols, 
       Scalar A_out[],
       const Scalar A_in[],
       const LocalOrdinal lda_in) const;

      void
      un_cache_block (const LocalOrdinal nrows,
          const LocalOrdinal ncols,
          Scalar A_out[],
          const LocalOrdinal lda_out,       
          const Scalar A_in[]) const;

      FactorOutput 
      factor (const LocalOrdinal nrows,
        const LocalOrdinal ncols, 
        Scalar A[],
        const LocalOrdinal lda,
        Scalar R[],
        const LocalOrdinal ldr,
        const bool contiguous_cache_blocks = false);

      void
      apply (const std::string& op,
       const LocalOrdinal nrows,
       const LocalOrdinal ncols_C,
       Scalar C[],
       const LocalOrdinal ldc,
       const LocalOrdinal ncols_Q,
       const Scalar Q[],
       const LocalOrdinal ldq,
       const FactorOutput& factor_output,
       const bool contiguous_cache_blocks = false);

      void 
      explicit_Q (const LocalOrdinal nrows,
      const LocalOrdinal ncols_Q_in,
      const Scalar Q_in[],
      const LocalOrdinal ldq_in,
      const LocalOrdinal ncols_Q_out,
      Scalar Q_out[],
      const LocalOrdinal ldq_out,
      const FactorOutput& factor_output,
      const bool contiguous_cache_blocks = false);

    private:
      size_t ncores_;
      TSQR::SequentialTsqr< LocalOrdinal, Scalar > seq_;
      Partitioner< LocalOrdinal, Scalar > partitioner_;

      typedef MatView< LocalOrdinal, Scalar > mat_view;
      typedef ConstMatView< LocalOrdinal, Scalar > const_mat_view;
      typedef std::pair< const_mat_view, const_mat_view > const_split_t;
      typedef std::pair< mat_view, mat_view > split_t;
      typedef std::pair< const_mat_view, mat_view > top_blocks_t;
      typedef std::vector< top_blocks_t > array_top_blocks_t;

      void
      explicit_Q_helper (const size_t P_first, 
       const size_t P_last,
       MatView< LocalOrdinal, Scalar >& Q_out,
       const bool contiguous_cache_blocks);

      MatView< LocalOrdinal, Scalar >
      factor_helper (const size_t P_first, 
         const size_t P_last,
         const size_t depth,
         MatView< LocalOrdinal, Scalar > A,
         std::vector< SeqOutput >& seq_outputs,
         ParOutput& par_outputs,
         Scalar R[],
         const LocalOrdinal ldr,
         const bool contiguous_cache_blocks);

      bool
      apply_helper_empty (const size_t P_first,
        const size_t P_last,
        const_mat_view &Q,
        mat_view& C) const;

      void
      build_partition_array (const size_t P_first,
           const size_t P_last,
           array_top_blocks_t& top_blocks,
           const_mat_view& Q,
           mat_view& C,
           const bool contiguous_cache_blocks) const;

      void
      apply_helper (const size_t P_first, 
        const size_t P_last,
        const_mat_view Q,
        mat_view C,
        array_top_blocks_t& top_blocks, 
        const FactorOutput& factor_output,
        const bool contiguous_cache_blocks);

      std::pair< ConstMatView< LocalOrdinal, Scalar >, MatView< LocalOrdinal, Scalar > >
      apply_transpose_helper (const std::string& op,
            const size_t P_first, 
            const size_t P_last,
            const_mat_view Q,
            mat_view C,
            const FactorOutput& factor_output,
            const bool contiguous_cache_blocks);

      void 
      factor_pair (const size_t P_top,
       const size_t P_bot,
       mat_view& A_top,
       mat_view& A_bot,
       std::vector< std::vector< Scalar > >& par_outputs,
       const bool contiguous_cache_blocks);

      void
      apply_pair (const std::string& trans,
      const size_t P_top,
      const size_t P_bot,
      const_mat_view& Q_bot,
      const std::vector< std::vector< Scalar > >& tau_arrays,
      mat_view& C_top,
      mat_view& C_bot,
      const bool contiguous_cache_blocks);

      void 
      cache_block_helper (MatView< LocalOrdinal, Scalar >& A_out,
        ConstMatView< LocalOrdinal, Scalar >& A_in,
        const size_t P_first,
        const size_t P_last) const;

      void 
      un_cache_block_helper (MatView< LocalOrdinal, Scalar >& A_out,
           const ConstMatView< LocalOrdinal, Scalar >& A_in,
           const size_t P_first,
           const size_t P_last) const;

    }; // class TbbRecursiveTsqr
  } // namespace TBB
} // namespace TSQR

#include <TSQR/TBB/TbbRecursiveTsqr_Def.hpp>

#endif // __TSQR_TbbRecursiveTsqr_hpp