Stokhos_TotalOrderBasis.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                           Stokhos Package
//
// Copyright 2009 NTESS and the Stokhos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef STOKHOS_TOTAL_ORDER_BASIS_HPP
#define STOKHOS_TOTAL_ORDER_BASIS_HPP

#include "Teuchos_RCP.hpp"

#include "Stokhos_ProductBasis.hpp"
#include "Stokhos_OneDOrthogPolyBasis.hpp"
#include "Stokhos_ProductBasisUtils.hpp"

namespace Stokhos {

  template <typename ordinal_type, typename value_type,
      typename coeff_compare_type = 
      TotalOrderLess<MultiIndex<ordinal_type> > >
  class TotalOrderBasis : 
    public ProductBasis<ordinal_type,value_type> {
  public:


    TotalOrderBasis(
      const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,
 value_type> > >& bases,
      const value_type& sparse_tol = 1.0e-12,
      const coeff_compare_type& coeff_compare = coeff_compare_type());

    virtual ~TotalOrderBasis();



    ordinal_type order() const;

    ordinal_type dimension() const;

    virtual ordinal_type size() const;


    virtual const Teuchos::Array<value_type>& norm_squared() const;

    virtual const value_type& norm_squared(ordinal_type i) const;


    virtual 
    Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > 
    computeTripleProductTensor() const;

    virtual 
    Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > 
    computeLinearTripleProductTensor() const;

    virtual value_type evaluateZero(ordinal_type i) const;


    virtual void evaluateBases(
      const Teuchos::ArrayView<const value_type>& point,
      Teuchos::Array<value_type>& basis_vals) const;

    virtual void print(std::ostream& os) const;

    virtual const std::string& getName() const;





    virtual const MultiIndex<ordinal_type>& term(ordinal_type i) const;


    virtual ordinal_type index(const MultiIndex<ordinal_type>& term) const;


    Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, 
                 value_type> > > 
    getCoordinateBases() const;

    virtual MultiIndex<ordinal_type> getMaxOrders() const;


  private:

    // Prohibit copying
    TotalOrderBasis(const TotalOrderBasis&);

    // Prohibit Assignment
    TotalOrderBasis& operator=(const TotalOrderBasis& b);
    
  protected:

    typedef MultiIndex<ordinal_type> coeff_type;
    typedef std::map<coeff_type,ordinal_type,coeff_compare_type> coeff_set_type;
    typedef Teuchos::Array<coeff_type> coeff_map_type;

    std::string name;

    ordinal_type p;

    ordinal_type d;

    ordinal_type sz;

    Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > > bases;

    value_type sparse_tol;

    coeff_type max_orders;

    coeff_set_type basis_set;

    coeff_map_type basis_map;

    Teuchos::Array<value_type> norms;

    mutable Teuchos::Array< Teuchos::Array<value_type> > basis_eval_tmp;

  }; // class TotalOrderBasis

  // An approach for building a sparse 3-tensor only for lexicographically
  // ordered total order basis
  // To-do:
  //   * Remove the n_choose_k() calls
  //   * Remove the loops in the Cijk_1D_Iterator::increment() functions
  //   * Store the 1-D Cijk tensors in a compressed format and eliminate
  //     the implicit searches with getValue()
  //     * Instead of looping over (i,j,k) multi-indices we could just store
  //       the 1-D Cijk tensors as an array of (i,j,k,c) tuples.
  template <typename ordinal_type, 
      typename value_type>
  Teuchos::RCP< Sparse3Tensor<ordinal_type, value_type> >
  computeTripleProductTensorLTO(
    const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
    bool symmetric = false) {
#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
    TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
#endif
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::Array;

    typedef MultiIndex<ordinal_type> coeff_type;
    const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
    ordinal_type d = bases.size();
    //ordinal_type p = product_basis.order();
    Array<ordinal_type> basis_orders(d);
    for (int i=0; i<d; ++i)
      basis_orders[i] = bases[i]->order();

    // Create 1-D triple products
    Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
    for (ordinal_type i=0; i<d; i++) {
      Cijk_1d[i] = 
  bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
    }

    RCP< Sparse3Tensor<ordinal_type, value_type> > Cijk = 
      rcp(new Sparse3Tensor<ordinal_type, value_type>);
      
    // Create i, j, k iterators for each dimension
    typedef ProductBasisUtils::Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
    Array<Cijk_Iterator> Cijk_1d_iterators(d);
    coeff_type terms_i(d,0), terms_j(d,0), terms_k(d,0);
    Array<ordinal_type> sum_i(d,0), sum_j(d,0), sum_k(d,0);
    for (ordinal_type dim=0; dim<d; dim++) {
      Cijk_1d_iterators[dim] = Cijk_Iterator(bases[dim]->order(), symmetric);
    }

    ordinal_type I = 0;
    ordinal_type J = 0;
    ordinal_type K = 0;
    ordinal_type cnt = 0;
    bool stop = false;
    while (!stop) {

      // Fill out terms from 1-D iterators
      for (ordinal_type dim=0; dim<d; ++dim) {
  terms_i[dim] = Cijk_1d_iterators[dim].i;
  terms_j[dim] = Cijk_1d_iterators[dim].j;
  terms_k[dim] = Cijk_1d_iterators[dim].k;
      }

      // Compute global I,J,K
      /*
      ordinal_type II = lexicographicMapping(terms_i, p);
      ordinal_type JJ = lexicographicMapping(terms_j, p);
      ordinal_type KK = lexicographicMapping(terms_k, p);
      if (I != II || J != JJ || K != KK) {
  std::cout << "DIFF!!!" << std::endl;
  std::cout << terms_i << ":  I = " << I << ", II = " << II << std::endl;
  std::cout << terms_j << ":  J = " << J << ", JJ = " << JJ << std::endl;
  std::cout << terms_k << ":  K = " << K << ", KK = " << KK << std::endl;
      }
      */

      // Compute triple-product value
      value_type c = value_type(1.0);
      for (ordinal_type dim=0; dim<d; dim++) {
  c *= Cijk_1d[dim]->getValue(Cijk_1d_iterators[dim].i, 
            Cijk_1d_iterators[dim].j, 
            Cijk_1d_iterators[dim].k);
      }

      TEUCHOS_TEST_FOR_EXCEPTION(
  std::abs(c) <= 1.0e-12,
  std::logic_error,
  "Got 0 triple product value " << c 
  << ", I = " << I << " = " << terms_i
  << ", J = " << J << " = " << terms_j
  << ", K = " << K << " = " << terms_k
  << std::endl);
      
      // Add term to global Cijk
      Cijk->add_term(I,J,K,c);
      // Cijk->add_term(I,K,J,c);
      // Cijk->add_term(J,I,K,c);
      // Cijk->add_term(J,K,I,c);
      // Cijk->add_term(K,I,J,c);
      // Cijk->add_term(K,J,I,c);
      
      // Increment iterators to the next valid term
      ordinal_type cdim = d-1;
      bool inc = true;
      while (inc && cdim >= 0) {
  ordinal_type delta_i, delta_j, delta_k;
  bool more = 
    Cijk_1d_iterators[cdim].increment(delta_i, delta_j, delta_k);

  // Update number of terms used for computing global index
  if (cdim == d-1) {
    I += delta_i;
    J += delta_j;
    K += delta_k;
  }
  else {
    if (delta_i > 0) {
      for (ordinal_type ii=0; ii<delta_i; ++ii)
        I += 
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_i[cdim] -
      (Cijk_1d_iterators[cdim].i-ii)+d-cdim,
      d-cdim-1);
    }
    else {
      for (ordinal_type ii=0; ii<-delta_i; ++ii)
        I -= 
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_i[cdim] -
      (Cijk_1d_iterators[cdim].i+ii)+d-cdim-1,
      d-cdim-1);
    }

    if (delta_j > 0) {
      for (ordinal_type jj=0; jj<delta_j; ++jj)
        J +=  
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_j[cdim] -
      (Cijk_1d_iterators[cdim].j-jj)+d-cdim,
      d-cdim-1);
    }
    else {
      for (ordinal_type jj=0; jj<-delta_j; ++jj)
        J -=  
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_j[cdim] -
      (Cijk_1d_iterators[cdim].j+jj)+d-cdim-1,
      d-cdim-1);
    }
    
    if (delta_k > 0) {
      for (ordinal_type kk=0; kk<delta_k; ++kk)
        K +=  
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_k[cdim] -
      (Cijk_1d_iterators[cdim].k-kk)+d-cdim,
      d-cdim-1);
    }
    else {
      for (ordinal_type kk=0; kk<-delta_k; ++kk)
        K -=  
    Stokhos::n_choose_k(
      basis_orders[cdim+1]-sum_k[cdim] -
      (Cijk_1d_iterators[cdim].k+kk)+d-cdim-1,
      d-cdim-1);
    }
  }
  
  if (!more) {
    // If no more terms in this dimension, go to previous one
    --cdim;
  }
  else {
    // cdim has more terms, so reset iterators for all dimensions > cdim
    // adjusting max order based on sum of i,j,k for previous dims
    inc = false;

    for (ordinal_type dim=cdim+1; dim<d; ++dim) {

      // Update sums of orders for previous dimension
      sum_i[dim] = sum_i[dim-1] + Cijk_1d_iterators[dim-1].i;
      sum_j[dim] = sum_j[dim-1] + Cijk_1d_iterators[dim-1].j;
      sum_k[dim] = sum_k[dim-1] + Cijk_1d_iterators[dim-1].k;

      // Reset iterator for this dimension
      Cijk_1d_iterators[dim] = 
        Cijk_Iterator(basis_orders[dim]-sum_i[dim],
          basis_orders[dim]-sum_j[dim],
          basis_orders[dim]-sum_k[dim],
          symmetric);
    }
  }
      }
      
      if (cdim < 0)
  stop = true;
      
      cnt++;
    }
    
    Cijk->fillComplete();
    
    return Cijk;
  }

} // Namespace Stokhos

// Include template definitions
#include "Stokhos_TotalOrderBasisImp.hpp"

#endif