43#include "Teuchos_TimeMonitor.hpp"
45#include "EpetraExt_BlockMultiVector.h"
46#include "Teuchos_Assert.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>& mean_prec_factory_,
56 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& G_prec_factory_,
57 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
58 label(
"Stokhos Kronecker Product Preconditioner"),
61 epetraCijk(epetraCijk_),
64 mean_prec_factory(mean_prec_factory_),
65 G_prec_factory(G_prec_factory_),
71 Cijk(epetraCijk->getParallelCijk()),
73 only_use_linear(false)
75 scale_op =
params->get(
"Scale Operator by Inverse Basis Norms",
true);
101 sg_poly = sg_op->getSGPolynomial();
102 mean_prec = mean_prec_factory->compute(sg_poly->getCoeffPtr(0));
103 label = std::string(
"Stokhos Kronecker Product Preconditioner:\n") +
104 std::string(
" ***** ") +
105 std::string(mean_prec->Label());
108 Teuchos::RCP<const Epetra_CrsGraph> graph = epetraCijk->getStochasticGraph();
111 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
113 Teuchos::RCP<Epetra_CrsMatrix> A0 =
114 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_poly->getCoeffPtr(0),
true);
115 double traceAB0 = MatrixTrace(*A0, *A0);
118 if (only_use_linear) {
119 int dim = sg_basis->dimension();
120 k_end = Cijk->find_k(dim+1);
124 Teuchos::RCP<Epetra_CrsMatrix> A_k =
125 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(sg_poly->getCoeffPtr(k),
127 double traceAB = MatrixTrace(*A_k, *A0);
129 j_it != Cijk->j_end(k_it); ++j_it) {
130 int j = epetraCijk->GCID(index(j_it));
132 i_it != Cijk->i_end(j_it); ++i_it) {
133 int i = epetraCijk->GRID(index(i_it));
134 double c = value(i_it)*traceAB/traceAB0;
137 G->SumIntoGlobalValues(i, 1, &c, &
j);
144 G_prec = G_prec_factory->compute(G);
146 label = std::string(
"Stokhos Kronecker Product Preconditioner:\n") +
147 std::string(
" ***** ") +
148 std::string(mean_prec->Label()) + std::string(
"\n") +
149 std::string(
" ***** ") +
150 std::string(G_prec->Label());
157 useTranspose = UseTranspose;
158 mean_prec->SetUseTranspose(useTranspose);
159 TEUCHOS_TEST_FOR_EXCEPTION(
160 UseTranspose ==
true, std::logic_error,
161 "Stokhos::KroneckerProductPreconditioner::SetUseTranspose(): " <<
162 "Preconditioner does not support transpose!" << std::endl);
171 EpetraExt::BlockMultiVector sg_input(
View, *base_map, Input);
172 EpetraExt::BlockMultiVector sg_result(
View, *base_map, Result);
176 int vecLen = sg_input.GetBlock(0)->MyLength();
177 int m = sg_input.NumVectors();
179 if (result_MVT == Teuchos::null || result_MVT->NumVectors() != vecLen*m) {
186 for (
int i=0; i<NumMyElements; i++) {
187 mean_prec->Apply(*(sg_input.GetBlock(i)), *(sg_result.GetBlock(i)));
190 Teuchos::RCP<Epetra_MultiVector> x;
191 for (
int irow=0 ; irow<NumMyElements; irow++) {
192 x = sg_result.GetBlock(irow);
193 for (
int vcol=0; vcol<m; vcol++) {
194 for (
int icol=0; icol<vecLen; icol++) {
195 (*result_MVT)[m*vcol+icol][irow] = (*x)[vcol][icol];
201 G_prec->Apply(*result_MVT, *result_MVT);
203 for (
int irow=0; irow<NumMyElements; irow++) {
204 x = sg_result.GetBlock(irow);
205 for (
int vcol=0; vcol<m; vcol++) {
206 for (
int icol=0; icol<vecLen; icol++) {
207 (*x)[vcol][icol] = (*result_MVT)[m*vcol+icol][irow];
219#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
220 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Kronecker Product Prec Time");
223 EpetraExt::BlockMultiVector sg_input(
View, *base_map, Input);
224 EpetraExt::BlockMultiVector sg_result(
View, *base_map, Result);
228 int vecLen = sg_input.GetBlock(0)->MyLength();
229 int m = sg_input.NumVectors();
231 if (result_MVT == Teuchos::null || result_MVT->NumVectors() != vecLen*m) {
238 Teuchos::RCP<Epetra_MultiVector> x;
239 for (
int irow=0 ; irow<NumMyElements; irow++) {
240 x = sg_input.GetBlock(irow);
241 for (
int vcol=0; vcol<m; vcol++) {
242 for (
int icol=0; icol<vecLen; icol++) {
243 (*result_MVT)[m*vcol+icol][irow] = (*x)[vcol][icol];
250#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
251 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: G Preconditioner Apply Inverse");
253 G_prec->ApplyInverse(*result_MVT, *result_MVT);
256 for (
int irow=0; irow<NumMyElements; irow++) {
257 x = sg_result.GetBlock(irow);
258 for (
int vcol=0; vcol<m; vcol++) {
259 for (
int icol=0; icol<vecLen; icol++) {
260 (*x)[vcol][icol] = (*result_MVT)[m*vcol+icol][irow];
267#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
268 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Mean Preconditioner Apply Inverse");
270 for (
int i=0; i<NumMyElements; i++) {
271 mean_prec->ApplyInverse(*(sg_result.GetBlock(i)),
272 *(sg_result.GetBlock(i)));
284 return mean_prec->NormInf() * G_prec->NormInf();
292 return const_cast<char*
>(label.c_str());
306 return mean_prec->NormInf() && G_prec->NormInf();
333 double traceAB = 0.0;
334 for (
int i=0; i<n; i++) {
336 for (
int j=0;
j<m;
j++) {
337 traceAB += A[i][
j]*B[i][
j];
int NumMyElements() const
int NumMyEntries(int Row) const
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 double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
virtual bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
bool only_use_linear
Limit construction of G to linear terms.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
virtual ~KroneckerProductPreconditioner()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > params
Preconditioner parameters.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
bool scale_op
Flag indicating whether operator be scaled with <\psi_i^2>
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
double MatrixTrace(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B) const
Compute trace of matrix A'B.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
KroneckerProductPreconditioner(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 > &mean_prec_factory, const Teuchos::RCP< Stokhos::AbstractPreconditionerFactory > &G_prec_factory, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Constructor.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
virtual int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
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.
Abstract base class for multivariate orthogonal polynomials.
Bi-directional iterator for traversing a sparse array.