43 #include "Teuchos_TimeMonitor.hpp" 47 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
49 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
50 const Teuchos::RCP<const Epetra_Map>& base_map_,
51 const Teuchos::RCP<const Epetra_Map>& sg_map_,
52 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& prec_factory_,
53 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
54 label(
"Stokhos Approximate Gauss-Seidel Preconditioner"),
57 epetraCijk(epetraCijk_),
60 prec_factory(prec_factory_),
65 Cijk(epetraCijk->getParallelCijk()),
68 only_use_linear(
true),
72 scale_op = params_->get(
"Scale Operator by Inverse Basis Norms",
true);
73 symmetric = params_->get(
"Symmetric Gauss-Seidel",
false);
88 sg_poly = sg_op->getSGPolynomial();
89 mean_prec = prec_factory->compute(sg_poly->getCoeffPtr(0));
90 label = std::string(
"Stokhos Approximate Gauss-Seidel Preconditioner:\n") +
91 std::string(
" ***** ") +
92 std::string(mean_prec->Label());
99 useTranspose = UseTheTranspose;
100 sg_op->SetUseTranspose(UseTheTranspose);
101 mean_prec->SetUseTranspose(UseTheTranspose);
110 return sg_op->Apply(Input, Result);
117 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 118 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Approximate Gauss-Seidel Time");
124 bool made_copy =
false;
131 if (mat_vec_tmp == Teuchos::null || mat_vec_tmp->NumVectors() != m)
133 if (rhs_block == Teuchos::null || rhs_block->NumVectors() != m)
135 Teuchos::rcp(
new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
138 EpetraExt::BlockMultiVector input_block(
View, *base_map, *input);
139 EpetraExt::BlockMultiVector result_block(
View, *base_map, Result);
141 result_block.PutScalar(0.0);
143 int k_limit = sg_poly->size();
145 k_limit = sg_poly->basis()->dimension() + 1;
146 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
148 rhs_block->Update(1.0, input_block, 0.0);
151 i_it!=Cijk->i_end(); ++i_it) {
154 Teuchos::RCP<Epetra_MultiVector> res_i = result_block.GetBlock(i);
157 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 158 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total AGS Deterministic Preconditioner Time");
160 mean_prec->ApplyInverse(*(rhs_block->GetBlock(i)), *res_i);
163 int i_gid = epetraCijk->GRID(i);
165 k_it != Cijk->k_end(i_it); ++k_it) {
167 if (k!=0 && k<k_limit) {
168 bool do_mat_vec =
false;
170 j_it != Cijk->j_end(k_it); ++j_it) {
172 int j_gid = epetraCijk->GCID(
j);
174 bool on_proc = epetraCijk->myGRID(j_gid);
182 (*sg_poly)[k].Apply(*res_i, *mat_vec_tmp);
184 j_it != Cijk->j_end(k_it); ++j_it) {
186 int j_gid = epetraCijk->GCID(
j);
187 double c = value(j_it);
195 bool on_proc = epetraCijk->myGRID(j_gid);
197 rhs_block->GetBlock(
j)->Update(-c, *mat_vec_tmp, 1.0);
210 i_it!=Cijk->i_rend(); ++i_it) {
213 Teuchos::RCP<Epetra_MultiVector> res_i = result_block.GetBlock(i);
216 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 217 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total AGS Deterministic Preconditioner Time");
219 mean_prec->ApplyInverse(*(rhs_block->GetBlock(i)), *res_i);
222 int i_gid = epetraCijk->GRID(i);
224 k_it != Cijk->k_end(i_it); ++k_it) {
226 if (k!=0 && k<k_limit) {
227 bool do_mat_vec =
false;
229 j_it != Cijk->j_end(k_it); ++j_it) {
231 int j_gid = epetraCijk->GCID(
j);
233 bool on_proc = epetraCijk->myGRID(j_gid);
241 (*sg_poly)[k].Apply(*res_i, *mat_vec_tmp);
243 j_it != Cijk->j_end(k_it); ++j_it) {
245 int j_gid = epetraCijk->GCID(
j);
246 double c = value(j_it);
250 bool on_proc = epetraCijk->myGRID(j_gid);
252 rhs_block->GetBlock(
j)->Update(-c, *mat_vec_tmp, 1.0);
272 return sg_op->NormInf();
280 return const_cast<char*
>(label.c_str());
294 return sg_op->HasNormInf();
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
ApproxGaussSeidelPreconditioner(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 ~ApproxGaussSeidelPreconditioner()
Destructor.
Bi-directional iterator for traversing a sparse array.
Bi-directional reverse iterator for traversing a sparse array.
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 ...
virtual const char * Label() const
Returns a character string describing the operator.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
bool symmetric
Use symmetric Gauss-Seidel.
bool scale_op
Flag indicating whether operator be scaled with <^2>
bool only_use_linear
Limit Gauss-Seidel loop to linear terms.
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.
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.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix 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 ...