Xpetra_Operator.hpp

Go to the documentation of this file.

// @HEADER
//
// ***********************************************************************
//
//             Xpetra: A linear algebra interface package
//                  Copyright 2012 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// 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
//                    Jeremie Gaidamour (jngaida@sandia.gov)
//                    Jonathan Hu       (jhu@sandia.gov)
//                    Ray Tuminaro      (rstumin@sandia.gov)
//
// ***********************************************************************
//
// @HEADER

// WARNING: This code is experimental. Backwards compatibility should not be expected.

#ifndef XPETRA_OPERATOR_HPP
#define XPETRA_OPERATOR_HPP

#include <Kokkos_DefaultNode.hpp>
#include <Kokkos_DefaultKernels.hpp>

#include "Xpetra_ConfigDefs.hpp"
#include "Xpetra_Exceptions.hpp"

#include "Xpetra_MultiVector.hpp"
#include "Xpetra_CrsGraph.hpp"
#include "Xpetra_CrsMatrix.hpp"
#include "Xpetra_CrsMatrixFactory.hpp"
#include "Xpetra_OperatorView.hpp"
#include "Xpetra_StridedMap.hpp"
#include "Xpetra_StridedMapFactory.hpp"

#include <Teuchos_SerialDenseMatrix.hpp>
#include <Teuchos_Hashtable.hpp>

namespace Xpetra {

  typedef std::string viewLabel_t;

  template <class Scalar, 
            class LocalOrdinal  = int, 
            class GlobalOrdinal = LocalOrdinal, 
            class Node          = Kokkos::DefaultNode::DefaultNodeType, 
            class LocalMatOps   = typename Kokkos::DefaultKernels<Scalar,LocalOrdinal,Node>::SparseOps > //TODO: or BlockSparseOp ?
  class Operator : virtual public Teuchos::Describable {
    
    typedef Xpetra::Map<LocalOrdinal, GlobalOrdinal, Node> Map;
    typedef Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> CrsMatrix;
    typedef Xpetra::CrsGraph<LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> CrsGraph;
#ifdef HAVE_XPETRA_TPETRA
    typedef Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> TpetraCrsMatrix;
#endif
    typedef Xpetra::CrsMatrixFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node, LocalMatOps> CrsMatrixFactory;
    typedef Xpetra::OperatorView<LocalOrdinal, GlobalOrdinal, Node> OperatorView;

  public:
  


    Operator() { }

    virtual ~Operator() { }



    void CreateView(viewLabel_t viewLabel, const RCP<const Map> & rowMap, const RCP<const Map> & colMap) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == true, Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.CreateView(): a view labeled '" + viewLabel + "' already exist.");
      RCP<OperatorView> view = rcp(new OperatorView(rowMap, colMap));
      operatorViewTable_.put(viewLabel, view);
    }
    
    void CreateView(const viewLabel_t viewLabel, const RCP<Operator> & A, bool transposeA = false, const RCP<Operator> & B = Teuchos::null, bool transposeB = false) {
            
      RCP<const Map> domainMap = Teuchos::null;
      RCP<const Map> rangeMap  = Teuchos::null;

      // A has strided Maps
      if(A->IsView(viewLabel)) {
        Xpetra::viewLabel_t oldView = A->SwitchToView(viewLabel); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
        rangeMap = (transposeA) ? A->getColMap() : A->getRowMap();
        domainMap = (transposeA) ? A->getRowMap() : A->getColMap(); // overwrite if B != Teuchos::null
        oldView = A->SwitchToView(oldView);
      } else rangeMap = (transposeA) ? A->getDomainMap() : A->getRangeMap();
      
      // B has strided Maps
      if(B != Teuchos::null ) {
        if(B->IsView(viewLabel)) {
          Xpetra::viewLabel_t oldView = B->SwitchToView(viewLabel); // note: "stridedMaps are always non-overlapping (correspond to range and domain maps!)
          domainMap = (transposeB) ? B->getRowMap() : B->getColMap();
          oldView = B->SwitchToView(oldView);
        } else domainMap = (transposeB) ? B->getRangeMap() : B->getDomainMap();
      }
      
      
      if(IsView(viewLabel)) RemoveView(viewLabel);
      
      CreateView(viewLabel, rangeMap, domainMap);   
    }

    void PrintViews(Teuchos::FancyOStream &out) const {
      int last = out.getOutputToRootOnly();
      Teuchos::OSTab tab(out);
      out.setOutputToRootOnly(0);
      Teuchos::Array<viewLabel_t> viewLabels;
      Teuchos::Array<RCP<OperatorView> > viewList;
      operatorViewTable_.arrayify(viewLabels,viewList);
      out << "views associated with this operator" << std::endl;
      for (int i=0; i<viewLabels.size(); ++i)
        out << viewLabels[i] << std::endl;
      out.setOutputToRootOnly(last);
    }

    
    void RemoveView(const viewLabel_t viewLabel) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.RemoveView(): view '" + viewLabel + "' does not exist.");
      TEUCHOS_TEST_FOR_EXCEPTION(viewLabel == GetDefaultViewLabel(), Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.RemoveView(): view '" + viewLabel + "' is the default view and cannot be removed.");
      operatorViewTable_.remove(viewLabel);
    }
    
    const viewLabel_t SwitchToView(const viewLabel_t viewLabel) {
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.SwitchToView(): view '" + viewLabel + "' does not exist.");
      viewLabel_t oldViewLabel = GetCurrentViewLabel();
      currentViewLabel_ = viewLabel;
      return oldViewLabel;
    }

    bool IsView(const viewLabel_t viewLabel) const {
      return operatorViewTable_.containsKey(viewLabel);
    }

    const viewLabel_t SwitchToDefaultView() { return SwitchToView(GetDefaultViewLabel()); }

    const viewLabel_t & GetDefaultViewLabel() const { return defaultViewLabel_; }

    const viewLabel_t & GetCurrentViewLabel() const { return currentViewLabel_; }
    




    virtual void insertGlobalValues(GlobalOrdinal globalRow, const ArrayView<const GlobalOrdinal> &cols, const ArrayView<const Scalar> &vals) = 0;


    virtual void insertLocalValues(LocalOrdinal localRow, const ArrayView<const LocalOrdinal> &cols, const ArrayView<const Scalar> &vals) = 0;




    virtual void fillComplete(const RCP<const Map> &domainMap, const RCP<const Map> &rangeMap, const RCP<ParameterList> &params = null) =0;

    //TODO : Get ride of "Tpetra"::OptimizeOption
    virtual void fillComplete(const RCP<ParameterList> &params=null) =0;




    virtual const RCP<const Map> & getRowMap() const { return getRowMap(GetCurrentViewLabel()); }

    virtual const RCP<const Map> & getRowMap(viewLabel_t viewLabel) const { 
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.GetRowMap(): view '" + viewLabel + "' does not exist.");
      return operatorViewTable_.get(viewLabel)->GetRowMap(); 
    }

    virtual const RCP<const Map> & getColMap() const { return getColMap(GetCurrentViewLabel()); }

    virtual const RCP<const Map> & getColMap(viewLabel_t viewLabel) const { 
      TEUCHOS_TEST_FOR_EXCEPTION(operatorViewTable_.containsKey(viewLabel) == false, Xpetra::Exceptions::RuntimeError, "Xpetra::Operator.GetColMap(): view '" + viewLabel + "' does not exist.");
      return operatorViewTable_.get(viewLabel)->GetColMap(); 
    }


    virtual global_size_t getGlobalNumRows() const =0;


    virtual global_size_t getGlobalNumCols() const =0;

    virtual size_t getNodeNumRows() const =0;

    virtual global_size_t getGlobalNumEntries() const =0;

    virtual size_t getNodeNumEntries() const =0;


    virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const =0;


    virtual global_size_t getGlobalNumDiags() const =0;


    virtual size_t getNodeNumDiags() const =0;


    virtual size_t getGlobalMaxNumRowEntries() const =0;


    virtual size_t getNodeMaxNumRowEntries() const =0;

    virtual bool isLocallyIndexed() const =0;

    virtual bool isGloballyIndexed() const =0;

    virtual bool isFillComplete() const =0;


    virtual void getLocalRowCopy(LocalOrdinal LocalRow, 
                                 const ArrayView<LocalOrdinal> &Indices, 
                                 const ArrayView<Scalar> &Values,
                                 size_t &NumEntries
                                 ) const =0;


    virtual void getGlobalRowView(GlobalOrdinal GlobalRow, ArrayView<const GlobalOrdinal> &indices, ArrayView<const Scalar> &values) const =0;


    virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView<const LocalOrdinal> &indices, ArrayView<const Scalar> &values) const =0;


    virtual void getLocalDiagCopy(Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &diag) const =0;

    virtual typename ScalarTraits<Scalar>::magnitudeType getFrobeniusNorm() const = 0;





    //TODO virtual=0 // TODO: Add default parameters ?
//     virtual void multiply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, Teuchos::ETransp trans, Scalar alpha, Scalar beta) const=0;





    virtual void apply(const MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & X, MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> &Y, 
                       Teuchos::ETransp mode = Teuchos::NO_TRANS,
                       Scalar alpha = ScalarTraits<Scalar>::one(),
                       Scalar beta = ScalarTraits<Scalar>::zero()) const =0;

    virtual const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getDomainMap() const =0;

    virtual const RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > getRangeMap() const =0;




    // TODO: describe of views can be done here

    //   /** \brief Return a simple one-line description of this object. */
    //   virtual std::string description() const =0;

    //   /** \brief Print the object with some verbosity level to an FancyOStream object. */
    //   virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const =0;





    virtual std::string description() const =0;

    virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const =0;

    // JG: Added:

    virtual RCP<const CrsGraph> getCrsGraph() const =0;

    // ----------------------------------------------------------------------------------
    // "TEMPORARY" VIEW MECHANISM
    // TODO: the view mechanism should be implemented as in MueMat.
    void SetFixedBlockSize(LocalOrdinal blksize) {

      TEUCHOS_TEST_FOR_EXCEPTION(isFillComplete() == false, Exceptions::RuntimeError, "Xpetra::Operator::SetFixedBlockSize(): operator is not filled and completed."); // TODO: do we need this? we just wanna "copy" the domain and range map

      std::vector<size_t> stridingInfo;
      stridingInfo.push_back(Teuchos::as<size_t>(blksize));
      LocalOrdinal stridedBlockId = -1;
            
      RCP<const Xpetra::StridedMap<LocalOrdinal, GlobalOrdinal, Node> > stridedRangeMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(getRangeMap()->lib(),
                                                    getRangeMap(),
                                                    stridingInfo,
                                                    stridedBlockId
                                                    );
      RCP<const Map> stridedDomainMap = Xpetra::StridedMapFactory<LocalOrdinal, GlobalOrdinal, Node>::Build(getDomainMap()->lib(),
                getDomainMap(),
                stridingInfo,
                stridedBlockId
                );

      if(IsView("stridedMaps") == true) RemoveView("stridedMaps");
      CreateView("stridedMaps", stridedRangeMap, stridedDomainMap); 
    }

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

    LocalOrdinal GetFixedBlockSize() const {
      if(IsView("stridedMaps")==true) {
        Teuchos::RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> > rangeMap = Teuchos::rcp_dynamic_cast<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >(getRowMap("stridedMaps"));
        Teuchos::RCP<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> > domainMap = Teuchos::rcp_dynamic_cast<const StridedMap<LocalOrdinal, GlobalOrdinal, Node> >(getColMap("stridedMaps"));
        TEUCHOS_TEST_FOR_EXCEPTION(rangeMap  == Teuchos::null, Exceptions::BadCast, "Xpetra::Operator::GetFixedBlockSize(): rangeMap is not of type StridedMap");
        TEUCHOS_TEST_FOR_EXCEPTION(domainMap == Teuchos::null, Exceptions::BadCast, "Xpetra::Operator::GetFixedBlockSize(): domainMap is not of type StridedMap");
        TEUCHOS_TEST_FOR_EXCEPTION(domainMap->getFixedBlockSize() != rangeMap->getFixedBlockSize(), Exceptions::RuntimeError, "Xpetra::Operator::GetFixedBlockSize(): block size of rangeMap and domainMap are different.");
        return Teuchos::as<LocalOrdinal>(domainMap->getFixedBlockSize()); // TODO: why LocalOrdinal?
      } else
        //TEUCHOS_TEST_FOR_EXCEPTION(false, Exceptions::RuntimeError, "Xpetra::Operator::GetFixedBlockSize(): no strided maps available."); // TODO remove this
        return 1;
    }; //TODO: why LocalOrdinal?

    // ----------------------------------------------------------------------------------

    protected:
      Teuchos::Hashtable<viewLabel_t, RCP<OperatorView> > operatorViewTable_; // hashtable storing the operator views (keys = view names, values = views).

      viewLabel_t defaultViewLabel_;  // label of the view associated with inital Operator construction
      viewLabel_t currentViewLabel_;  // label of the current view

  }; //class Operator

} //namespace Xpetra

#define XPETRA_OPERATOR_SHORT
#endif //XPETRA_OPERATOR_DECL_HPP