43 #include "Teuchos_TimeMonitor.hpp" 44 #include "Epetra_LocalMap.h" 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);
98 const Epetra_Vector&
x)
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();
112 G = Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *graph));
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);
169 Apply(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 171 EpetraExt::BlockMultiVector sg_input(View, *base_map, Input);
172 EpetraExt::BlockMultiVector sg_result(View, *base_map, Result);
174 const Epetra_Map& G_map = G->RowMap();
175 int NumMyElements = G_map.NumMyElements();
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) {
180 result_MVT = Teuchos::rcp(
new Epetra_MultiVector(G_map, 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];
217 ApplyInverse(
const Epetra_MultiVector& Input, Epetra_MultiVector& Result)
const 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);
226 const Epetra_Map& G_map = G->RowMap();
227 int NumMyElements = G_map.NumMyElements();
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) {
232 result_MVT = Teuchos::rcp(
new Epetra_MultiVector(G_map, 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();
331 MatrixTrace(
const Epetra_CrsMatrix& A,
const Epetra_CrsMatrix& B)
const {
332 int n = A.NumMyRows();
333 double traceAB = 0.0;
334 for (
int i=0; i<n; i++) {
335 int m = A.NumMyEntries(i);
336 for (
int j=0;
j<m;
j++) {
337 traceAB += A[i][
j]*B[i][
j];
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
bool only_use_linear
Limit construction of G to linear terms.
virtual ~KroneckerProductPreconditioner()
Destructor.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
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.
Bi-directional iterator for traversing a sparse array.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
bool scale_op
Flag indicating whether operator be scaled with <^2>
Teuchos::RCP< Teuchos::ParameterList > params
Preconditioner parameters.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerking basis.
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.
double MatrixTrace(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B) const
Compute trace of matrix A'B.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
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 int SetUseTranspose(bool UseTranspose)
Set to true if the transpose of the operator is requested.
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.