44#include "Epetra_config.h"
45#include "EpetraExt_BlockMultiVector.h"
46#include "EpetraExt_MatrixMatrix.h"
50 const Teuchos::RCP<const Epetra_Map>& domain_base_map_,
51 const Teuchos::RCP<const Epetra_Map>& range_base_map_,
52 const Teuchos::RCP<const Epetra_Map>& domain_sg_map_,
53 const Teuchos::RCP<const Epetra_Map>& range_sg_map_,
56 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops_)
57 : label(
"Stokhos Diagonal Operator"),
58 domain_base_map(domain_base_map_),
59 range_base_map(range_base_map_),
60 domain_sg_map(domain_sg_map_),
61 range_sg_map(range_sg_map_),
66 expansion_size(sg_basis->size()),
67 num_blocks(block_ops->size()),
68 input_block(expansion_size),
69 result_block(expansion_size),
81 const Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >& ops)
86Teuchos::RCP<const Stokhos::EpetraOperatorOrthogPoly >
92Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly >
101 useTranspose = UseTranspose;
102 for (
int i=0; i<num_blocks; i++)
103 (*block_ops)[i].SetUseTranspose(useTranspose);
109Stokhos::DiagEpetraOp::Apply(std::vector< Teuchos::RCP< const Epetra_CrsMatrix> >& sg_J_all, std::vector< Teuchos::RCP< Epetra_CrsMatrix> >& sg_Kkk_all)
const
115 for (
int i=0;i<num_blocks+1;i++) {
116 Teuchos::RCP<Epetra_CrsMatrix> sg_J_poly_Crs =
117 Teuchos::rcp_dynamic_cast< Epetra_CrsMatrix>((*block_ops).getCoeffPtr(i),
true);
118 sg_J_all.push_back(sg_J_poly_Crs);
128 for(
int k=0;k<expansion_size;k++) {
129 Teuchos::RCP<Epetra_CrsMatrix> Kkk =
131 sg_Kkk_all.push_back(Kkk);
132 sg_Kkk_all[k]->PutScalar(0.0);
136 for (
int k=0; k<expansion_size; k++) {
137 int nj = Cijk->num_j(k);
138 const Teuchos::Array<int>& j_indices = Cijk->Jindices(k);
139 for (
int jj=0; jj<nj; jj++) {
140 int j = j_indices[jj];
142 const Teuchos::Array<double>& cijk_values = Cijk->values(k,jj);
143 const Teuchos::Array<int>& i_indices = Cijk->Iindices(k,jj);
144 int ni = i_indices.size();
145 for (
int ii=0; ii<ni; ii++) {
146 int i = i_indices[ii];
147 if (i<num_blocks+1) {
148 double cikk = cijk_values[ii];
149 EpetraExt::MatrixMatrix::Add((*sg_J_all[i]),
false, cikk, *(sg_Kkk_all[k]), 1.0);
204 throw "DiagEpetraOp::ApplyInverse not defined!";
218 return const_cast<char*
>(label.c_str());
236 return domain_base_map->Comm();
241 return *domain_sg_map;
247 return *range_sg_map;
virtual void reset(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &ops)
Reset operator blocks.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual Teuchos::RCP< const Stokhos::EpetraOperatorOrthogPoly > getOperatorBlocks() const
Get operator blocks.
DiagEpetraOp(const Teuchos::RCP< const Epetra_Map > &domain_base_map_, const Teuchos::RCP< const Epetra_Map > &range_base_map_, const Teuchos::RCP< const Epetra_Map > &domain_sg_map_, const Teuchos::RCP< const Epetra_Map > &range_sg_map_, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &Cijk, const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &ops)
Constructor.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual const char * Label() const
Returns a character string describing the operator.
virtual int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of the inverse of the operator applied to a Epetra_MultiVector Input in Result as ...
virtual ~DiagEpetraOp()
Destructor.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual int Apply(std::vector< Teuchos::RCP< const Epetra_CrsMatrix > > &sg_J_all, std::vector< Teuchos::RCP< Epetra_CrsMatrix > > &sg_Kkk_all) const
Returns Diagonal blocks of SG matrix when PC coefficients of the SG matrix are given.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
Abstract base class for multivariate orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.