43#include "Teuchos_TimeMonitor.hpp"
50 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
52 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
53 const Teuchos::RCP<const Epetra_Map>& base_map_,
54 const Teuchos::RCP<const Epetra_Map>& sg_map_,
55 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& prec_factory_,
56 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
57 label(
"Stokhos Approximate Schur Complement Preconditioner"),
60 epetraCijk(epetraCijk_),
63 prec_factory(prec_factory_),
68 Cijk(epetraCijk->getParallelCijk()),
71 upper_block_Cijk(P+1),
72 lower_block_Cijk(P+1),
75 only_use_linear(
true),
79 TEUCHOS_TEST_FOR_EXCEPTION(
80 epetraCijk->isStochasticParallel(), std::logic_error,
81 "Stokhos::ApproxSchurComplementPreconditioner does not support " <<
82 "a parallel stochastic distribution.");
84 scale_op = params_->get(
"Scale Operator by Inverse Basis Norms",
true);
85 symmetric = params_->get(
"Symmetric Gauss-Seidel",
false);
95 int nj =
Cijk->num_j(k);
101 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > prod_basis =
102 Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<int,double> >(
sg_basis,
true);
103 int d = prod_basis->dimension();
104 MultiIndex<int> term(d);
105 for (
int p=0; p<=
P; p++) {
125 Teuchos::Array<int>::iterator col_it =
132 double c = value(i_it);
133 Teuchos::Array<int>::iterator row_it =
141 else if (p_col < p_row) {
148 for (
int p=0; p<=
P; p++) {
165 sg_poly = sg_op->getSGPolynomial();
166 mean_prec = prec_factory->compute(sg_poly->getCoeffPtr(0));
167 label = std::string(
"Stokhos Approximate Schur Complement Preconditioner:\n")
168 + std::string(
" ***** ") + std::string(mean_prec->Label());
175 useTranspose = UseTranspose;
176 sg_op->SetUseTranspose(UseTranspose);
177 mean_prec->SetUseTranspose(UseTranspose);
186 return sg_op->Apply(Input, Result);
193#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
194 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Approximate Schur Complement Time");
200 bool made_copy =
false;
201 if (Input.Values() == Result.Values()) {
208 if (rhs_block == Teuchos::null || rhs_block->NumVectors() != m)
210 Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
211 if (tmp == Teuchos::null || tmp->NumVectors() != m*max_num_mat_vec)
214 j_ptr.resize(m*max_num_mat_vec);
215 mj_indices.resize(m*max_num_mat_vec);
218 EpetraExt::BlockMultiVector input_block(
View, *base_map, *input);
219 EpetraExt::BlockMultiVector result_block(
View, *base_map, Result);
221 result_block.PutScalar(0.0);
224 rhs_block->Update(1.0, input_block, 0.0);
230 for (
int l=P; l>=1; l--) {
232 divide_diagonal_block(block_indices[l], block_indices[l+1],
233 *rhs_block, result_block);
236 multiply_block(upper_block_Cijk[l], -1.0, result_block, *rhs_block);
240 divide_diagonal_block(0, 1, *rhs_block, result_block);
242 for (
int l=1; l<=P; l++) {
244 multiply_block(lower_block_Cijk[l], -1.0, result_block, *rhs_block);
247 divide_diagonal_block(block_indices[l], block_indices[l+1],
248 *rhs_block, result_block);
261 return sg_op->NormInf();
269 return const_cast<char*
>(label.c_str());
283 return sg_op->HasNormInf();
311 const EpetraExt::BlockMultiVector& Input,
312 EpetraExt::BlockMultiVector& Result)
const
316 int m = Input.NumVectors();
317 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
321 k_end = cijk->find_k(sg_basis()->dimension() + 1);
326 int nj = cijk->num_j(k_it);
331 for (
int mm=0; mm<m; mm++) {
332 j_ptr[l*m+mm] = (*(Input.GetBlock(
j)))[mm];
333 mj_indices[l*m+mm] = l*m+mm;
339 (*sg_poly)[k].Apply(input_tmp, result_tmp);
346 double c = value(i_it);
349 for (
int mm=0; mm<m; mm++)
350 (*Result.GetBlock(i))(mm)->
Update(alpha*c, *result_tmp(l*m+mm), 1.0);
361 const EpetraExt::BlockMultiVector& Input,
362 EpetraExt::BlockMultiVector& Result)
const
364 for (
int i=row_begin; i<row_end; i++) {
365#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
366 TEUCHOS_FUNC_TIME_MONITOR(
367 "Stokhos: ASC Deterministic Preconditioner Time");
369 mean_prec->ApplyInverse(*(Input.GetBlock(i)), *(Result.GetBlock(i)));
virtual int Apply(const Epetra_MultiVector &Input, Epetra_MultiVector &Result) const
Returns the result of a Epetra_Operator applied to a Epetra_MultiVector Input in Result as described ...
bool only_use_linear
Limit Gauss-Seidel loop to linear terms.
Teuchos::Array< Teuchos::RCP< Cijk_type > > lower_block_Cijk
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
void multiply_block(const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &cijk, double alpha, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
int P
Total polynomial order.
virtual const char * Label() const
Returns a character string describing the operator.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual ~ApproxSchurComplementPreconditioner()
Destructor.
Teuchos::Array< int > block_indices
Starting block indices.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
Teuchos::RCP< const Cijk_type > Cijk
Pointer to triple product.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
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 ...
bool symmetric
Use symmetric Gauss-Seidel.
int max_num_mat_vec
Maximum number of matvecs in Apply.
void divide_diagonal_block(int row_begin, int row_end, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
ApproxSchurComplementPreconditioner(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_Map > &base_map, const Teuchos::RCP< const Epetra_Map > &sg_map, const Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > &prec_factory, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Constructor.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
Teuchos::Array< Teuchos::RCP< Cijk_type > > upper_block_Cijk
Triple product tensor for each sub-block.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
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.