Teuchos_SerialQRDenseSolver_UQ_PCE.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                           Stokhos Package
//                 Copyright (2009) 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.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
//
// ***********************************************************************
// @HEADER

#ifndef _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_
#define _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_

#include <new>
#include "Teuchos_SerialQRDenseSolver.hpp"

#include "Sacado_UQ_PCE.hpp"

namespace Teuchos {

  template<typename OrdinalType, typename Storage>
  class SerialQRDenseSolver<OrdinalType, Sacado::UQ::PCE<Storage> >
    : public CompObject,
      public Object
  {

  public:

    typedef Sacado::UQ::PCE<Storage> ScalarType;
    typedef typename ScalarTraits<ScalarType>::magnitudeType MagnitudeType;


    SerialQRDenseSolver();

    virtual ~SerialQRDenseSolver() {}




    int setMatrix(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& A);


    int setVectors(const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& X,
                   const RCP<SerialDenseMatrix<OrdinalType, ScalarType> >& B);




    void factorWithEquilibration(bool flag) {
      base_QR_.factorWithEquilibration(flag);
    }

    void solveWithTranspose(bool flag) {
      base_QR_.solveWithTranspose(flag);
    }

    void solveWithTransposeFlag(Teuchos::ETransp trans) {
      base_QR_.solveWithTranspose(trans);
    }





    int factor() { return base_QR_.factor(); }


    int solve() { return base_QR_.solve(); }


    int computeEquilibrateScaling() {
      return base_QR_.computeEquilibrateScaling();
    }


    int equilibrateMatrix() { return base_QR_.equilibrateMatrix(); }


    int equilibrateRHS() { return base_QR_.equilibrateRHS(); }


    int unequilibrateLHS() { return base_QR_.unequilibrateLHS(); }


    int formQ();


    int formR();


    int multiplyQ (ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C);


    int solveR (ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C);



    bool transpose() { return base_QR_.transpose(); }

    bool factored() { return base_QR_.factored(); }

    bool equilibratedA() { return base_QR_.equilibratedA(); }

    bool equilibratedB() { return base_QR_.equilibratedB(); }

    bool shouldEquilibrate() { return base_QR_.shouldEquilibrate(); }

    bool solved() { return base_QR_.solved(); }

    bool formedQ() { return base_QR_.formedQ(); }

    bool formedR() { return base_QR_.formedR();}




    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getMatrix() const {return(Matrix_);}

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getFactoredMatrix()  const {return(Factor_);}

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getQ()  const {return(FactorQ_);}

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getR()  const {return(FactorR_);}

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getLHS() const {return(LHS_);}

    RCP<SerialDenseMatrix<OrdinalType, ScalarType> > getRHS() const {return(RHS_);}

    OrdinalType numRows()  const {return(M_);}

    OrdinalType numCols()  const {return(N_);}

    std::vector<ScalarType> tau()  const { return base_QR_.tau(); }

    MagnitudeType ANORM()  const { return base_QR_.ANORM(); }



    void Print(std::ostream& os) const;
  protected:

    typedef typename ScalarType::value_type BaseScalarType;
    typedef SerialQRDenseSolver<OrdinalType,BaseScalarType> BaseQRType;
    typedef SerialDenseMatrix<OrdinalType, BaseScalarType> BaseMatrixType;
    typedef SerialDenseMatrix<OrdinalType, ScalarType> MatrixType;

    void resetMatrix();
    void resetVectors();
    RCP<BaseMatrixType> createBaseMatrix(const RCP<MatrixType>& mat) const;
    RCP<MatrixType> createMatrix(const RCP<BaseMatrixType>& base_mat) const;
    BaseQRType base_QR_;

    OrdinalType M_;
    OrdinalType N_;
    OrdinalType SacadoSize_;

    RCP<MatrixType> Matrix_;
    RCP<MatrixType> LHS_;
    RCP<MatrixType> RHS_;
    RCP<MatrixType> Factor_;
    RCP<MatrixType> FactorQ_;
    RCP<MatrixType> FactorR_;

    RCP<BaseMatrixType> Base_Matrix_;
    RCP<BaseMatrixType> Base_LHS_;
    RCP<BaseMatrixType> Base_RHS_;
    RCP<BaseMatrixType> Base_Factor_;
    RCP<BaseMatrixType> Base_FactorQ_;
    RCP<BaseMatrixType> Base_FactorR_;

  private:

    // SerialQRDenseSolver copy constructor (put here because we don't want user access)
    SerialQRDenseSolver(const SerialQRDenseSolver& Source);
    SerialQRDenseSolver & operator=(const SerialQRDenseSolver& Source);

  };

  namespace details {

    // Helper traits class for converting between arrays of
    // Sacado::UQ::PCE and its corresponding value type
    template <typename Storage> struct PCEArrayHelper;

    template <typename Ordinal, typename Value, typename Device>
    struct PCEArrayHelper< Stokhos::DynamicStorage<Ordinal,Value,Device> >
    {
      typedef Stokhos::DynamicStorage<Ordinal,Value,Device> Storage;
      typedef Sacado::UQ::PCE<Storage> pce_type;
      typedef typename pce_type::cijk_type cijk_type;

      static Value* getValueArray(pce_type* v, Ordinal len) {
        if (len == 0)
          return 0;
        return v[0].coeff(); // Assume data is contiguous
      }

      static pce_type* getPCEArray(Value *v, Ordinal len, Ordinal vector_size){
        if (len == 0)
          return 0;
        pce_type* v_pce =
          static_cast<pce_type*>(operator new(len*sizeof(pce_type)));
        cijk_type cijk = Kokkos::getGlobalCijkTensor<cijk_type>();
        for (Ordinal i=0; i<len; ++i)
          new (v_pce+i) pce_type(cijk, vector_size, v+i*vector_size, false);
        return v_pce;
      }

    };

  }

  // Helper traits to distinguish work arrays for real and complex-valued datatypes.
  using namespace details;

//=============================================================================

template<typename OrdinalType, typename Storage>
SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
SerialQRDenseSolver()
  : CompObject(),
    Object("Teuchos::SerialQRDenseSolver"),
    M_(0),
    N_(0)
{
  resetMatrix();
}

//=============================================================================

template<typename OrdinalType, typename Storage>
void SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
resetVectors()
{
  LHS_ = Teuchos::null;
  RHS_ = Teuchos::null;
}
//=============================================================================

template<typename OrdinalType, typename Storage>
void SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
resetMatrix()
{
  resetVectors();
  M_ = 0;
  N_ = 0;
}
//=============================================================================

template<typename OrdinalType, typename Storage>
RCP<SerialDenseMatrix<OrdinalType, typename Sacado::UQ::PCE<Storage>::value_type> >
SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
createBaseMatrix(
  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& mat) const
{
  BaseScalarType* base_ptr =
    PCEArrayHelper<Storage>::getValueArray(
      mat->values(), mat->stride()*mat->numCols());
  RCP<BaseMatrixType> base_mat =
    rcp(new BaseMatrixType(Teuchos::View,
                           base_ptr,
                           mat->stride()*SacadoSize_,
                           mat->numRows()*SacadoSize_,
                           mat->numCols()));
  return base_mat;
}
//=============================================================================

template<typename OrdinalType, typename Storage>
RCP<SerialDenseMatrix<OrdinalType, Sacado::UQ::PCE<Storage> > >
SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
createMatrix(
  const RCP<SerialDenseMatrix<OrdinalType,BaseScalarType> >& base_mat) const
{
  ScalarType* ptr =
    PCEArrayHelper<Storage>::getPCEArray(
      base_mat->values(), base_mat->stride()*base_mat->numCols(), SacadoSize_);
  RCP<MatrixType> mat =
    rcp(new MatrixType(Teuchos::View,
                       ptr,
                       base_mat->stride()/SacadoSize_,
                       base_mat->numRows()/SacadoSize_,
                       base_mat->numCols()));
  return mat;
}
//=============================================================================

template<typename OrdinalType, typename Storage>
int SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
setMatrix(const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& A)
{
  TEUCHOS_TEST_FOR_EXCEPTION(
    A->numRows() < A->numCols(), std::invalid_argument,
    "SerialQRDenseSolver<T>::setMatrix: the matrix A must have A.numRows() >= A.numCols()!");

  resetMatrix();
  Matrix_ = A;
  Factor_ = A;
  FactorQ_ = A;
  FactorR_ = A;
  M_ = A->numRows();
  N_ = A->numCols();
  if (Storage::is_static)
    SacadoSize_ = Storage::static_size;
  else if (M_ > 0 && N_ > 0)
    SacadoSize_ = (*A)(0,0).size();
  else
    SacadoSize_ = 0;

  return base_QR_.setMatrix( createBaseMatrix(A) );
}
//=============================================================================

template<typename OrdinalType, typename Storage>
int SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
setVectors(
  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& X,
  const RCP<SerialDenseMatrix<OrdinalType,ScalarType> >& B)
{
  // Check that these new vectors are consistent
  TEUCHOS_TEST_FOR_EXCEPTION(
    B->numCols() != X->numCols(), std::invalid_argument,
    "SerialQRDenseSolver<T>::setVectors: X and B have different numbers of columns!");
  TEUCHOS_TEST_FOR_EXCEPTION(
    B->values()==0, std::invalid_argument,
    "SerialQRDenseSolver<T>::setVectors: B is an empty SerialDenseMatrix<T>!");
  TEUCHOS_TEST_FOR_EXCEPTION(
    X->values()==0, std::invalid_argument,
    "SerialQRDenseSolver<T>::setVectors: X is an empty SerialDenseMatrix<T>!");
  TEUCHOS_TEST_FOR_EXCEPTION(
    B->stride()<1, std::invalid_argument,
    "SerialQRDenseSolver<T>::setVectors: B has an invalid stride!");
  TEUCHOS_TEST_FOR_EXCEPTION(
    X->stride()<1, std::invalid_argument,
    "SerialQRDenseSolver<T>::setVectors: X has an invalid stride!");

  resetVectors();
  LHS_ = X;
  RHS_ = B;

  return base_QR_.setVectors( createBaseMatrix(X), createBaseMatrix(B) );
}

//=============================================================================

template<typename OrdinalType, typename Storage>
int SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
formQ() {
  int ret = base_QR_.formQ();
  FactorQ_ = createMatrix( base_QR_.getQ() );
  return(ret);
}

//=============================================================================

template<typename OrdinalType, typename Storage>
int SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
formR() {
  int ret = base_QR_.formR();
  FactorR_ = createMatrix( base_QR_.getR() );
  Factor_ = createMatrix( base_QR_.getFactoredMatrix() );
  return(ret);
}

//=============================================================================

template<typename OrdinalType, typename Storage>
int  SerialQRDenseSolver<OrdinalType, Sacado::UQ::PCE<Storage> >::
multiplyQ(ETransp transq, SerialDenseMatrix<OrdinalType, ScalarType> &C)
{
  return base_QR_.multiplyQ(transq, createBaseMatrix(C));
}

//=============================================================================

template<typename OrdinalType, typename Storage>
int  SerialQRDenseSolver<OrdinalType, Sacado::UQ::PCE<Storage> >::
solveR(ETransp transr, SerialDenseMatrix<OrdinalType, ScalarType> &C)
{
  return base_QR_.solveR(transr, createBaseMatrix(C));
}

//=============================================================================

template<typename OrdinalType, typename Storage>
void SerialQRDenseSolver<OrdinalType,Sacado::UQ::PCE<Storage> >::
Print(std::ostream& os) const {

  if (Matrix_!=Teuchos::null) os << "Solver Matrix"          << std::endl << *Matrix_ << std::endl;
  if (Factor_!=Teuchos::null) os << "Solver Factored Matrix" << std::endl << *Factor_ << std::endl;
  if (LHS_   !=Teuchos::null) os << "Solver LHS"             << std::endl << *LHS_    << std::endl;
  if (RHS_   !=Teuchos::null) os << "Solver RHS"             << std::endl << *RHS_    << std::endl;

}

} // namespace Teuchos

#endif /* _TEUCHOS_SERIALQRDENSESOLVER_UQ_PCE_HPP_ */