TsqrAdaptor.hpp

Go to the documentation of this file.

// @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_Trilinos_TsqrAdaptor_hpp
#define __TSQR_Trilinos_TsqrAdaptor_hpp


#include "AnasaziConfigDefs.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"

#include "TsqrTypeAdaptor.hpp"
#include "TsqrCommFactory.hpp"

#include "Tsqr_GlobalVerify.hpp"
#include "Tsqr_ScalarTraits.hpp"

#include <stdexcept>
#include <sstream>


namespace TSQR {
  namespace Trilinos {

    template< class S, class LO, class GO, class MV >
    class TsqrAdaptor {
    public:
      typedef S   scalar_type;
      typedef LO  local_ordinal_type;
      typedef GO  global_ordinal_type;
      typedef MV  multivector_type;

      typedef typename TSQR::ScalarTraits< scalar_type >::magnitude_type magnitude_type;

      typedef TsqrTypeAdaptor< S, LO, GO, MV >      type_adaptor;
      typedef typename type_adaptor::factory_type   factory_type;

      typedef typename type_adaptor::node_tsqr_type node_tsqr_type;
      typedef typename type_adaptor::node_tsqr_ptr  node_tsqr_ptr;

      typedef typename type_adaptor::comm_type      comm_type;
      typedef typename type_adaptor::comm_ptr       comm_ptr;

      typedef typename type_adaptor::dist_tsqr_type dist_tsqr_type;
      typedef typename type_adaptor::dist_tsqr_ptr  dist_tsqr_ptr;

      typedef typename type_adaptor::tsqr_type      tsqr_type;
      typedef typename type_adaptor::tsqr_ptr       tsqr_ptr;

      typedef typename tsqr_type::FactorOutput      factor_output_type;
      typedef Teuchos::SerialDenseMatrix< LO, S >   dense_matrix_type;
      typedef Teuchos::RCP< MessengerBase< S > >    scalar_messenger_ptr;
      typedef Teuchos::RCP< MessengerBase< LO > >   ordinal_messenger_ptr;

      virtual ~TsqrAdaptor() {}

      virtual factor_output_type
      factor (multivector_type& A, 
        dense_matrix_type& R,
        const bool contiguousCacheBlocks = false)
      {
  local_ordinal_type nrowsLocal, ncols, LDA;
  fetchDims (A, nrowsLocal, ncols, LDA);
  // This is guaranteed to be _correct_ for any Node type, but
  // won't necessary be efficient.  The desired model is that
  // A_local requires no copying.
  Teuchos::ArrayRCP< scalar_type > A_local = fetchNonConstView (A);

  // Reshape R if necessary.  This operation zeros out all the
  // entries of R, which is what we want anyway.
  if (R.numRows() != ncols || R.numCols() != ncols)
    {
      if (0 != R.shape (ncols, ncols))
        throw std::runtime_error ("Failed to reshape matrix R");
    }
  return pTsqr_->factor (nrowsLocal, ncols, A_local.get(), LDA, 
             R.values(), R.stride(), contiguousCacheBlocks);
      }

      virtual void 
      explicitQ (const multivector_type& Q_in, 
     const factor_output_type& factorOutput,
     multivector_type& Q_out, 
     const bool contiguousCacheBlocks = false)
      {
  using Teuchos::ArrayRCP;

  local_ordinal_type nrowsLocal, ncols_in, LDQ_in;
  fetchDims (Q_in, nrowsLocal, ncols_in, LDQ_in);
  local_ordinal_type nrowsLocal_out, ncols_out, LDQ_out;
  fetchDims (Q_out, nrowsLocal_out, ncols_out, LDQ_out);

  if (nrowsLocal_out != nrowsLocal)
    {
      std::ostringstream os;
      os << "TSQR explicit Q: input Q factor\'s node-local part has a di"
        "fferent number of rows (" << nrowsLocal << ") than output Q fac"
        "tor\'s node-local part (" << nrowsLocal_out << ").";
      throw std::runtime_error (os.str());
    }
  ArrayRCP< const scalar_type > pQin = fetchConstView (Q_in);
  ArrayRCP< scalar_type > pQout = fetchNonConstView (Q_out);
  pTsqr_->explicit_Q (nrowsLocal, 
          ncols_in, pQin.get(), LDQ_in, 
          factorOutput,
          ncols_out, pQout.get(), LDQ_out,
          contiguousCacheBlocks);
      }

      local_ordinal_type
      revealRank (multivector_type& Q,
      dense_matrix_type& R,
      const magnitude_type relativeTolerance,
      const bool contiguousCacheBlocks = false) const
      {
  using Teuchos::ArrayRCP;

  local_ordinal_type nrowsLocal, ncols, ldqLocal;
  fetchDims (Q, nrowsLocal, ncols, ldqLocal);

  ArrayRCP< scalar_type > Q_ptr = fetchNonConstView (Q);
  return pTsqr_->reveal_rank (nrowsLocal, ncols, 
            Q_ptr.get(), ldqLocal,
            R.values(), R.stride(), 
            relativeTolerance, 
            contiguousCacheBlocks);
      }

      virtual void 
      cacheBlock (const multivector_type& A_in, 
      multivector_type& A_out)
      {
  using Teuchos::ArrayRCP;

  local_ordinal_type nrowsLocal, ncols, LDA_in;
  fetchDims (A_in, nrowsLocal, ncols, LDA_in);
  local_ordinal_type nrowsLocal_out, ncols_out, LDA_out;
  fetchDims (A_out, nrowsLocal_out, ncols_out, LDA_out);

  if (nrowsLocal_out != nrowsLocal)
    {
      std::ostringstream os;
      os << "TSQR cache block: the input matrix\'s node-local part has a"
        " different number of rows (" << nrowsLocal << ") than the outpu"
        "t matrix\'s node-local part (" << nrowsLocal_out << ").";
      throw std::runtime_error (os.str());
    }
  else if (ncols_out != ncols)
    {
      std::ostringstream os;
      os << "TSQR cache block: the input matrix\'s node-local part has a"
        " different number of columns (" << ncols << ") than the output "
        "matrix\'s node-local part (" << ncols_out << ").";
      throw std::runtime_error (os.str());
    }
  ArrayRCP< const scalar_type > pA_in = fetchConstView (A_in);
  ArrayRCP< scalar_type > pA_out = fetchNonConstView (A_out);
  pTsqr_->cache_block (nrowsLocal, ncols, pA_out.get(), 
           pA_in.get(), LDA_in);
      }

      virtual void 
      unCacheBlock (const multivector_type& A_in, 
        multivector_type& A_out)
      {
  using Teuchos::ArrayRCP;

  local_ordinal_type nrowsLocal, ncols, LDA_in;
  fetchDims (A_in, nrowsLocal, ncols, LDA_in);
  local_ordinal_type nrowsLocal_out, ncols_out, LDA_out;
  fetchDims (A_out, nrowsLocal_out, ncols_out, LDA_out);

  if (nrowsLocal_out != nrowsLocal)
    {
      std::ostringstream os;
      os << "TSQR un-cache-block: the input matrix\'s node-local part ha"
        "s a different number of rows (" << nrowsLocal << ") than the ou"
        "tput matrix\'s node-local part (" << nrowsLocal_out << ").";
      throw std::runtime_error (os.str());
    }
  else if (ncols_out != ncols)
    {
      std::ostringstream os;
      os << "TSQR cache block: the input matrix\'s node-local part has a"
        " different number of columns (" << ncols << ") than the output "
        "matrix\'s node-local part (" << ncols_out << ").";
      throw std::runtime_error (os.str());
    }
  ArrayRCP< const scalar_type > pA_in = fetchConstView (A_in);
  ArrayRCP< scalar_type > pA_out = fetchNonConstView (A_out);
  pTsqr_->un_cache_block (nrowsLocal, ncols, pA_out.get(), 
        LDA_out, pA_in.get());
      }

      virtual std::pair< magnitude_type, magnitude_type >
      verify (const multivector_type& A,
        const multivector_type& Q,
        const Teuchos::SerialDenseMatrix< local_ordinal_type, scalar_type >& R)
      {
  using Teuchos::ArrayRCP;

  local_ordinal_type nrowsLocal_A, ncols_A, LDA;
  local_ordinal_type nrowsLocal_Q, ncols_Q, LDQ;
  fetchDims (A, nrowsLocal_A, ncols_A, LDA);
  fetchDims (Q, nrowsLocal_Q, ncols_Q, LDQ);
  if (nrowsLocal_A != nrowsLocal_Q)
    throw std::runtime_error ("A and Q must have same number of rows");
  else if (ncols_A != ncols_Q)
    throw std::runtime_error ("A and Q must have same number of columns");
  else if (ncols_A != R.numCols())
    throw std::runtime_error ("A and R must have same number of columns");
  else if (R.numRows() < R.numCols())
    throw std::runtime_error ("R must have no fewer rows than columns");

  // Const views suffice for verification
  ArrayRCP< const scalar_type > A_ptr = fetchConstView (A);
  ArrayRCP< const scalar_type > Q_ptr = fetchConstView (Q);
  return global_verify (nrowsLocal_A, ncols_A, A_ptr.get(), LDA,
            Q_ptr.get(), LDQ, R.values(), R.stride(), 
            pScalarMessenger_.get());
      }

    protected:
      void 
      init (const multivector_type& mv,
      const Teuchos::ParameterList& plist)
      {
  // This is done in a multivector type - dependent way.
  fetchMessengers (mv, pScalarMessenger_, pOrdinalMessenger_);

  factory_type factory;
  // plist and pScalarMessenger_ are inputs.  Construct *pTsqr_.
  factory.makeTsqr (plist, pScalarMessenger_, pTsqr_);
      }

    private:
      virtual void 
      fetchDims (const multivector_type& A, 
     local_ordinal_type& nrowsLocal, 
     local_ordinal_type& ncols, 
     local_ordinal_type& LDA) const = 0;

      virtual Teuchos::ArrayRCP< scalar_type > 
      fetchNonConstView (multivector_type& A) const = 0;

      virtual Teuchos::ArrayRCP< const scalar_type > 
      fetchConstView (const multivector_type& A) const = 0;

      virtual void
      fetchMessengers (const multivector_type& mv,
           scalar_messenger_ptr& pScalarMessenger,
           ordinal_messenger_ptr& pOrdinalMessenger) const = 0;

      scalar_messenger_ptr pScalarMessenger_;

      ordinal_messenger_ptr pOrdinalMessenger_;

      tsqr_ptr pTsqr_;
    };

  } // namespace Trilinos
} // namespace TSQR

#endif // __TSQR_Trilinos_TsqrAdaptor_hpp