Tsqr_CacheBlocker.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_CacheBlocker_hpp
#define __TSQR_CacheBlocker_hpp

#include <Tsqr_CacheBlockingStrategy.hpp>
#include <Tsqr_MatView.hpp>
#include <Tsqr_Util.hpp>

// #include <iostream>
#include <sstream>
#include <stdexcept>

// #define CACHE_BLOCKER_DEBUG 1
// #ifdef CACHE_BLOCKER_DEBUG
// #  undef CACHE_BLOCKER_DEBUG
// #endif // CACHE_BLOCKER_DEBUG
#ifdef CACHE_BLOCKER_DEBUG
#  include <iostream>
using std::cerr;
using std::endl;
#endif // CACHE_BLOCKER_DEBUG


namespace TSQR {

  template< class Ordinal, class Scalar > // class CBS=CacheBlockingStrategy< Ordinal, Scalar > 
  class CacheBlocker {
  private:
    typedef MatView< Ordinal, Scalar > mat_view;
    typedef ConstMatView< Ordinal, Scalar > const_mat_view;

    void
    validate () 
    {
      if (nrows_cache_block_ < ncols_)
  {
    std::ostringstream os;
    os << "Cache block is too small: only " << nrows_cache_block_
       << " rows fit, but need at least as many rows as columns (= "
       << ncols_ << ")";
    throw std::logic_error (os.str());
  }
    }

  public:
    CacheBlocker (const Ordinal num_rows,
      const Ordinal num_cols,
      const CacheBlockingStrategy< Ordinal, Scalar >& strategy) :
      nrows_ (num_rows), 
      ncols_ (num_cols), 
      strategy_ (strategy)
    {
#ifdef CACHE_BLOCKER_DEBUG
      cerr << "CacheBlocker:" << endl
     << "# rows = " << nrows() << endl
     << "# cols = " << ncols() << endl
     << "Cache block size (bytes) = " << strategy_.cache_block_size()
     << endl << endl;
#endif
      const Ordinal temp_nrows_cache_block = 
  strategy_.cache_block_num_rows (ncols());
#ifdef CACHE_BLOCKER_DEBUG
      cerr << "Strategy says: # rows per cache block = "
     << temp_nrows_cache_block << endl;
#endif
      nrows_cache_block_ = temp_nrows_cache_block;
      validate ();
    }

    CacheBlocker (const CacheBlocker& rhs) :
      nrows_ (rhs.nrows()), 
      ncols_ (rhs.ncols()), 
      strategy_ (rhs.strategy_)
    {}

    CacheBlocker& operator= (const CacheBlocker& rhs) {
      nrows_ = rhs.nrows();
      ncols_ = rhs.ncols();
      strategy_ = rhs.strategy_;
      return *this;
    }

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

    Ordinal nrows () const { return nrows_; }

    Ordinal ncols () const { return ncols_; }

    template< class MatrixViewType >
    MatrixViewType
    split_top_block (MatrixViewType& A, const bool contiguous_cache_blocks) const
    {
      typedef typename MatrixViewType::ordinal_type ordinal_type;
      const ordinal_type nrows_top = 
  strategy_.top_block_split_nrows (A.nrows(), ncols(), 
           nrows_cache_block());
      // split_top() modifies A
      return A.split_top (nrows_top, contiguous_cache_blocks);
    }

    template< class MatrixViewType >
    MatrixViewType
    top_block (const MatrixViewType& A, const bool contiguous_cache_blocks) const
    {
      typedef typename MatrixViewType::ordinal_type ordinal_type;
      const ordinal_type nrows_top = 
  strategy_.top_block_split_nrows (A.nrows(), ncols(), 
           nrows_cache_block());
      MatrixViewType A_copy (A);
      return A_copy.split_top (nrows_top, contiguous_cache_blocks);
    }

    template< class MatrixViewType >
    MatrixViewType
    split_bottom_block (MatrixViewType& A, const bool contiguous_cache_blocks) const
    {
      typedef typename MatrixViewType::ordinal_type ordinal_type;
      const ordinal_type nrows_bottom = 
  strategy_.bottom_block_split_nrows (A.nrows(), ncols(), 
              nrows_cache_block());
      // split_bottom() modifies A
      return A.split_bottom (nrows_bottom, contiguous_cache_blocks);
    }

    template< class MatrixViewType >
    void
    fill_with_zeros (MatrixViewType A,
         const bool contiguous_cache_blocks) const
    {
      // Note: if the cache blocks are stored contiguously, A.lda()
      // won't be the correct leading dimension of A, but it won't
      // matter: we only ever operate on A_cur here, and A_cur's
      // leading dimension is set correctly by split_top_block().
      while (! A.empty())
  {
    // This call modifies A, but that's OK since we passed the
    // input view by copy, not by reference.
    MatrixViewType A_cur = split_top_block (A, contiguous_cache_blocks);
    A_cur.fill (Scalar(0));
  }
    }

    void
    cache_block (const Ordinal num_rows,
     const Ordinal num_cols,
     Scalar A_out[],
     const Scalar A_in[],
     const Ordinal lda_in) const
    {
      // We say "*_rest" because it points to the remaining part of
      // the matrix left to cache block; at the beginning, the
      // "remaining" part is the whole matrix, but that will change as
      // the algorithm progresses.
      const_mat_view A_in_rest (num_rows, num_cols, A_in, lda_in);
      // Leading dimension doesn't matter since A_out will be cache blocked.
      mat_view A_out_rest (num_rows, num_cols, A_out, lda_in);

      while (! A_in_rest.empty())
  {
    if (A_out_rest.empty())
      throw std::logic_error("A_out_rest is empty, but A_in_rest is not");

    // This call modifies A_in_rest.
    const_mat_view A_in_cur = split_top_block (A_in_rest, false);

    // This call modifies A_out_rest.
    mat_view A_out_cur = split_top_block (A_out_rest, true);

    copy_matrix (A_in_cur.nrows(), num_cols, A_out_cur.get(), 
           A_out_cur.lda(), A_in_cur.get(), A_in_cur.lda());
  }
    }

    void
    un_cache_block (const Ordinal num_rows,
        const Ordinal num_cols,
        Scalar A_out[],
        const Ordinal lda_out,        
        const Scalar A_in[]) const
    {
      // We say "*_rest" because it points to the remaining part of
      // the matrix left to cache block; at the beginning, the
      // "remaining" part is the whole matrix, but that will change as
      // the algorithm progresses.
      //
      // Leading dimension doesn't matter since A_in is cache blocked.
      const_mat_view A_in_rest (num_rows, num_cols, A_in, lda_out);
      mat_view A_out_rest (num_rows, num_cols, A_out, lda_out);

      while (! A_in_rest.empty())
  {
    if (A_out_rest.empty())
      throw std::logic_error("A_out_rest is empty, but A_in_rest is not");

    // This call modifies A_in_rest.
    const_mat_view A_in_cur = split_top_block (A_in_rest, true);

    // This call modifies A_out_rest.
    mat_view A_out_cur = split_top_block (A_out_rest, false);

    copy_matrix (A_in_cur.nrows(), num_cols, A_out_cur.get(), 
           A_out_cur.lda(), A_in_cur.get(), A_in_cur.lda());
  }
    }


    void
    fill_with_zeros (const Ordinal num_rows,
         const Ordinal num_cols,
         Scalar A[],
         const Ordinal lda, 
         const bool contiguous_cache_blocks)
    {
      // We say "A_rest" because it points to the remaining part of
      // the matrix left to process; at the beginning, the "remaining"
      // part is the whole matrix, but that will change as the
      // algorithm progresses.
      //
      // Note: if the cache blocks are stored contiguously, lda won't
      // be the correct leading dimension of A, but it won't matter:
      // we only ever operate on A_cur here, and A_cur's leading
      // dimension is set correctly by A_rest.split_top().
      mat_view A_rest (num_rows, num_cols, A, lda);

      while (! A_rest.empty())
  {
    // This call modifies A_rest.
    mat_view A_cur = split_top_block (A_rest, contiguous_cache_blocks);
    A_cur.fill (Scalar(0));
  }
    }


  private:
    Ordinal nrows_;
    Ordinal ncols_;
    CacheBlockingStrategy< Ordinal, Scalar > strategy_;
    Ordinal nrows_cache_block_;

    size_t nrows_cache_block () const { return nrows_cache_block_; }
  };

} // namespace TSQR


#endif // __TSQR_CacheBlocker_hpp