Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ApproxJacobiPreconditioner.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
43#include "Teuchos_TimeMonitor.hpp"
45
48 const Teuchos::RCP<const EpetraExt::MultiComm>& sg_comm_,
49 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
50 const Teuchos::RCP<const Stokhos::EpetraSparse3Tensor>& epetraCijk_,
51 const Teuchos::RCP<const Epetra_Map>& base_map_,
52 const Teuchos::RCP<const Epetra_Map>& sg_map_,
53 const Teuchos::RCP<Stokhos::AbstractPreconditionerFactory>& prec_factory_,
54 const Teuchos::RCP<Teuchos::ParameterList>& params_) :
55 label("Stokhos Approximate Jacobi Preconditioner"),
56 sg_comm(sg_comm_),
57 sg_basis(sg_basis_),
58 epetraCijk(epetraCijk_),
59 base_map(base_map_),
60 sg_map(sg_map_),
61 prec_factory(prec_factory_),
62 mean_prec(),
63 useTranspose(false),
64 num_iter(2),
65 sg_op(),
66 sg_poly(),
67 rhs_block()
68{
69 num_iter = params_->get("Number of Jacobi Iterations", 2);
70
71 Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
72 Teuchos::rcp(&(params_->sublist("Jacobi SG Operator")), false);
73 sgOpParams->set("Include Mean", false);
74 if (!sgOpParams->isParameter("Only Use Linear Terms"))
75 sgOpParams->set("Only Use Linear Terms", true);
76
77 // Build new parallel Cijk if we are only using the linear terms, Cijk
78 // is distributed over proc's, and Cijk includes more than just the linear
79 // terms (so we have the right column map; otherwise we will be importing
80 // much more than necessary)
81 if (sgOpParams->get<bool>("Only Use Linear Terms") &&
82 epetraCijk->isStochasticParallel()) {
83 int dim = sg_basis->dimension();
84 if (epetraCijk->getKEnd() > dim+1)
85 epetraCijk =
86 Teuchos::rcp(new EpetraSparse3Tensor(*epetraCijk, 1, dim+1));
87
88 }
89 Stokhos::SGOperatorFactory sg_op_factory(sgOpParams);
90 mat_free_op = sg_op_factory.build(sg_comm, sg_basis, epetraCijk,
92}
93
96{
97}
98
99void
101setupPreconditioner(const Teuchos::RCP<Stokhos::SGOperator>& sg_op_,
102 const Epetra_Vector& x)
103{
104 sg_op = sg_op_;
105 sg_poly = sg_op->getSGPolynomial();
106 mean_prec = prec_factory->compute(sg_poly->getCoeffPtr(0));
107 label = std::string("Stokhos Approximate Jacobi Preconditioner:\n") +
108 std::string(" ***** ") +
109 std::string(mean_prec->Label());
110 mat_free_op->setupOperator(sg_poly);
111}
112
113int
115SetUseTranspose(bool UseTranspose)
116{
117 useTranspose = UseTranspose;
118 sg_op->SetUseTranspose(UseTranspose);
119 mat_free_op->SetUseTranspose(UseTranspose);
120 mean_prec->SetUseTranspose(UseTranspose);
121
122 return 0;
123}
124
125int
127Apply(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
128{
129 return sg_op->Apply(Input, Result);
130}
131
132int
134ApplyInverse(const Epetra_MultiVector& Input, Epetra_MultiVector& Result) const
135{
136#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
137 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Approximate Jacobi Time");
138#endif
139
140 // We have to be careful if Input and Result are the same vector.
141 // If this is the case, the only possible solution is to make a copy
142 const Epetra_MultiVector *input = &Input;
143 bool made_copy = false;
144 if (Input.Values() == Result.Values()) {
145 input = new Epetra_MultiVector(Input);
146 made_copy = true;
147 }
148
149 int m = input->NumVectors();
150 if (rhs_block == Teuchos::null || rhs_block->NumVectors() != m)
151 rhs_block =
152 Teuchos::rcp(new EpetraExt::BlockMultiVector(*base_map, *sg_map, m));
153
154 // Extract blocks
155 EpetraExt::BlockMultiVector input_block(View, *base_map, *input);
156 EpetraExt::BlockMultiVector result_block(View, *base_map, Result);
157
158 int myBlockRows = epetraCijk->numMyRows();
159 result_block.PutScalar(0.0);
160 for (int iter=0; iter<num_iter; iter++) {
161
162 // Compute RHS
163 if (iter == 0)
164 rhs_block->Update(1.0, input_block, 0.0);
165 else {
166 mat_free_op->Apply(result_block, *rhs_block);
167 rhs_block->Update(1.0, input_block, -1.0);
168 }
169
170 // Apply deterministic preconditioner
171 for(int i=0; i<myBlockRows; i++) {
172#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
173 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total AJ Deterministic Preconditioner Time");
174#endif
175 mean_prec->ApplyInverse(*(rhs_block->GetBlock(i)),
176 *(result_block.GetBlock(i)));
177 }
178
179 }
180
181 if (made_copy)
182 delete input;
183
184 return 0;
185}
186
187double
189NormInf() const
190{
191 return sg_op->NormInf();
192}
193
194
195const char*
197Label () const
198{
199 return const_cast<char*>(label.c_str());
200}
201
202bool
204UseTranspose() const
205{
206 return useTranspose;
207}
208
209bool
211HasNormInf() const
212{
213 return sg_op->HasNormInf();
214}
215
216const Epetra_Comm &
218Comm() const
219{
220 return *sg_comm;
221}
222const Epetra_Map&
224OperatorDomainMap() const
225{
226 return *sg_map;
227}
228
229const Epetra_Map&
231OperatorRangeMap() const
232{
233 return *sg_map;
234}
int NumVectors() const
double * Values() 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 const Epetra_Map & OperatorDomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Stores SG parallel communicator.
Teuchos::RCP< Stokhos::SGOperator > mat_free_op
SG operator to implement SG mat-vec.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)
Setup preconditioner.
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 bool HasNormInf() const
Returns true if the this object can provide an approximate Inf-norm, false otherwise.
virtual const char * Label() const
Returns a character string describing the operator.
virtual bool UseTranspose() const
Returns the current UseTranspose setting.
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.
Teuchos::RCP< const Epetra_Map > base_map
Stores base map.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Stores Epetra Cijk tensor.
virtual double NormInf() const
Returns an approximate infinity norm of the operator matrix.
virtual const Epetra_Map & OperatorRangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
virtual const Epetra_Comm & Comm() const
Returns a reference to the Epetra_Comm communicator associated with this operator.
ApproxJacobiPreconditioner(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 > &params)
Constructor.
Teuchos::RCP< const Epetra_Map > sg_map
Stores SG map.
Abstract base class for multivariate orthogonal polynomials.
Factory for generating stochastic Galerkin preconditioners.
virtual Teuchos::RCP< Stokhos::SGOperator > build(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 > &domain_base_map, const Teuchos::RCP< const Epetra_Map > &range_base_map, const Teuchos::RCP< const Epetra_Map > &domain_sg_map, const Teuchos::RCP< const Epetra_Map > &range_sg_map)
Build preconditioner operator.