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 =
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++) {
162 const Epetra_Vector&
x)
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);
184 Apply(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 186 return sg_op->Apply(Input, Result);
191 ApplyInverse(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 193 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 194 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Approximate Schur Complement Time");
199 const Epetra_MultiVector *input = &Input;
200 bool made_copy =
false;
201 if (Input.Values() == Result.Values()) {
202 input =
new Epetra_MultiVector(Input);
207 int m = input->NumVectors();
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)
212 tmp = Teuchos::rcp(
new Epetra_MultiVector(*base_map,
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;
337 Epetra_MultiVector input_tmp(View, *base_map, &j_ptr[0], nj*m);
338 Epetra_MultiVector result_tmp(View, *tmp, &mj_indices[0], nj*m);
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 const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int max_num_mat_vec
Maximum number of matvecs in Apply.
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. ...
virtual ~ApproxSchurComplementPreconditioner()
Destructor.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
Teuchos::RCP< const Cijk_type > Cijk
Pointer to triple product.
virtual ordinal_type dimension() const =0
Return dimension of basis.
Bi-directional iterator for traversing a sparse array.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
A multidimensional index.
void multiply_block(const Teuchos::RCP< const Stokhos::Sparse3Tensor< int, double > > &cijk, double alpha, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
bool symmetric
Use symmetric Gauss-Seidel.
Teuchos::Array< Teuchos::RCP< Cijk_type > > upper_block_Cijk
Triple product tensor for each sub-block.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
bool only_use_linear
Limit Gauss-Seidel loop to linear terms.
Teuchos::Array< Teuchos::RCP< Cijk_type > > lower_block_Cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
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 int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this 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 void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
void divide_diagonal_block(int row_begin, int row_end, const EpetraExt::BlockMultiVector &Input, EpetraExt::BlockMultiVector &Result) const
bool scale_op
Flag indicating whether operator be scaled with <^2>
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual const char * Label() const
Returns a character string describing the operator.
int P
Total polynomial order.
Teuchos::Array< int > block_indices
Starting block indices.
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 ...