48#include "EpetraExt_BlockUtility.h"
52 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
54 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
55 const Teuchos::RCP<const Epetra_CrsGraph>& determ_graph,
56 const Teuchos::RCP<Teuchos::ParameterList>& params) :
62 epetraCijk(epetraCijk_),
63 Cijk(epetraCijk->getParallelCijk()),
67 only_use_linear(false),
68 determOffset_(EpetraExt::BlockUtility::CalculateOffset(determ_graph->RowMap()))
70 scale_op = params->get(
"Scale Operator by Inverse Basis Norms",
true);
83 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
93 if (!include_mean && index(k_begin) == 0)
95 if (only_use_linear) {
96 int dim = sg_basis->dimension();
97 k_end = Cijk->find_k(dim+1);
101 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
104 Teuchos::RCP<Epetra_RowMatrix> block =
105 Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(block_ops->getCoeffPtr(k),
108 j_it != Cijk->j_end(k_it); ++j_it) {
109 int j = epetraCijk->GCID(index(j_it));
111 i_it != Cijk->i_end(j_it); ++i_it) {
112 int i = epetraCijk->GRID(index(i_it));
113 double c = value(i_it);
116 this->SumIntoGlobalBlock_Deterministic(c, *block, i,
j);
122Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly >
129Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
139 const Epetra_Map & determMap = determBlock.RowMatrixRowMap();
140 const Epetra_Map & determColMap = determBlock.RowMatrixColMap();
146 int MaxIndices = determBlock.MaxNumEntries();
147 std::vector<int> Indices(MaxIndices);
148 std::vector<double> Values(MaxIndices);
152 for (
int i=0; i<determMap.NumMyElements(); i++) {
153 determBlock.ExtractMyRowCopy( i, MaxIndices, NumIndices, &Values[0], &Indices[0] );
156 for(
int l = 0; l < NumIndices; ++l ) {
157 Indices[l] = Col + COffset_*determColMap.GID(Indices[l]);
161 int BaseRow = determMap.GID(i);
162 ierr = this->SumIntoGlobalValues(ROffset_*BaseRow + Row, NumIndices, &Values[0], &Indices[0]);
163 if (ierr != 0) std::cout <<
"WARNING InterlacedOperator::SumIntoBlock_Deterministic SumIntoGlobalValues err = " << ierr <<
164 "\n\t Row " << ROffset_*BaseRow + Row <<
", Col start " << Indices[0] << std::endl;
CRS matrix of dense blocks.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
void SumIntoGlobalBlock_Deterministic(double alpha, const Epetra_RowMatrix &determBlock, int Row, int Col)
Sum into global matrix.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)
Setup operator.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
InterlacedOperator(const Teuchos::RCP< const EpetraExt::MultiComm > &sg_comm, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > &epetraCijk, const Teuchos::RCP< const Epetra_CrsGraph > &base_graph, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Constructor.
virtual ~InterlacedOperator()
Destructor.
bool only_use_linear
Flag indicating whether to only use linear terms.
bool include_mean
Flag indicating whether to include mean term.
Abstract base class for multivariate orthogonal polynomials.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
ji_sparse_array::const_iterator kj_iterator
Iterator for looping over j entries given k.