Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
NOX_Epetra_LinearSystem_SGJacobi.cpp
Go to the documentation of this file.
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stokhos Package
6// Copyright (2009) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44
46
47#include "Epetra_Vector.h"
48#include "NOX_Epetra_Interface_Jacobian.H"
49#include "EpetraExt_BlockVector.h"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_TimeMonitor.hpp"
53
54NOX::Epetra::LinearSystemSGJacobi::
55LinearSystemSGJacobi(
56 Teuchos::ParameterList& printingParams,
57 Teuchos::ParameterList& linearSolverParams_,
58 const Teuchos::RCP<NOX::Epetra::LinearSystem>& det_solver_,
59 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
60 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
61 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis,
62 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data,
63 const Teuchos::RCP<Epetra_Operator>& J,
64 const Teuchos::RCP<const Epetra_Map>& base_map_,
65 const Teuchos::RCP<const Epetra_Map>& sg_map_,
66 const Teuchos::RCP<NOX::Epetra::Scaling> s):
67 det_solver(det_solver_),
68 epetraCijk(sg_parallel_data->getEpetraCijk()),
69 jacInterfacePtr(iJac),
70 base_map(base_map_),
71 sg_map(sg_map_),
72 scaling(s),
73 utils(printingParams)
74{
75 sg_op = Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(J, true);
76 sg_poly = sg_op->getSGPolynomial();
77
78 sg_df_block =
79 Teuchos::rcp(new EpetraExt::BlockVector(*base_map, *sg_map));
80 kx = Teuchos::rcp(new Epetra_Vector(*base_map));
81
82 Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
83 Teuchos::rcp(&(linearSolverParams_.sublist("Jacobi SG Operator")), false);
84 sgOpParams->set("Include Mean", false);
85 if (!sgOpParams->isParameter("Only Use Linear Terms"))
86 sgOpParams->set("Only Use Linear Terms", false);
87
88 // Build new parallel Cijk if we are only using the linear terms, Cijk
89 // is distributed over proc's, and Cijk includes more than just the linear
90 // terms (so we have the right column map; otherwise we will be importing
91 // much more than necessary)
92 if (sgOpParams->get<bool>("Only Use Linear Terms") &&
93 epetraCijk->isStochasticParallel()) {
94 int dim = sg_basis->dimension();
95 if (epetraCijk->getKEnd() > dim+1)
96 epetraCijk =
97 Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(
98 *epetraCijk, 1, dim+1));
99
100 }
101
102 Stokhos::SGOperatorFactory sg_op_factory(sgOpParams);
103 Teuchos::RCP<const EpetraExt::MultiComm> sg_comm =
104 sg_parallel_data->getMultiComm();
105 mat_free_op =
106 sg_op_factory.build(sg_comm, sg_basis, epetraCijk,
107 base_map, base_map, sg_map, sg_map);
108}
109
110NOX::Epetra::LinearSystemSGJacobi::
111~LinearSystemSGJacobi()
112{
113}
114
115bool NOX::Epetra::LinearSystemSGJacobi::
116applyJacobian(const NOX::Epetra::Vector& input,
117 NOX::Epetra::Vector& result) const
118{
119 sg_op->SetUseTranspose(false);
120 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
121
122 return (status == 0);
123}
124
125bool NOX::Epetra::LinearSystemSGJacobi::
126applyJacobianTranspose(const NOX::Epetra::Vector& input,
127 NOX::Epetra::Vector& result) const
128{
129 sg_op->SetUseTranspose(true);
130 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
131 sg_op->SetUseTranspose(false);
132
133 return (status == 0);
134}
135
136bool NOX::Epetra::LinearSystemSGJacobi::
137applyJacobianInverse(Teuchos::ParameterList &params,
138 const NOX::Epetra::Vector &input,
139 NOX::Epetra::Vector &result)
140{
141 int max_iter = params.get("Max Iterations",100 );
142 double sg_tol = params.get("Tolerance", 1e-12);
143
144 // Extract blocks
145 EpetraExt::BlockVector sg_dx_block(View, *base_map, result.getEpetraVector());
146 EpetraExt::BlockVector sg_f_block(View, *base_map, input.getEpetraVector());
147 sg_dx_block.PutScalar(0.0);
148
149 // Compute initial residual norm
150 double norm_f,norm_df;
151 sg_f_block.Norm2(&norm_f);
152 sg_op->Apply(sg_dx_block, *sg_df_block);
153 sg_df_block->Update(-1.0, sg_f_block, 1.0);
154 sg_df_block->Norm2(&norm_df);
155
156 Teuchos::RCP<Epetra_Vector> df, dx;
157 Teuchos::ParameterList& det_solver_params =
158 params.sublist("Deterministic Solver Parameters");
159
160 int myBlockRows = epetraCijk->numMyRows();
161 int iter = 0;
162 while (((norm_df/norm_f)>sg_tol) && (iter<max_iter)) {
163 TEUCHOS_FUNC_TIME_MONITOR("Total global solve Time");
164 iter++;
165
166 // Compute RHS
167 if (iter == 0)
168 sg_df_block->Update(1.0, sg_f_block, 0.0);
169 else {
170 mat_free_op->Apply(sg_dx_block, *sg_df_block);
171 sg_df_block->Update(1.0, sg_f_block, -1.0);
172 }
173
174 for (int i=0; i<myBlockRows; i++) {
175 df = sg_df_block->GetBlock(i);
176 dx = sg_dx_block.GetBlock(i);
177 NOX::Epetra::Vector nox_df(df, NOX::Epetra::Vector::CreateView);
178 NOX::Epetra::Vector nox_dx(dx, NOX::Epetra::Vector::CreateView);
179
180 (*sg_poly)[0].Apply(*dx, *kx);
181
182 dx->PutScalar(0.0);
183 // Solve linear system
184 {
185 TEUCHOS_FUNC_TIME_MONITOR("Total deterministic solve Time");
186 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
187 }
188
189 df->Update(-1.0, *kx, 1.0);
190 }
191
192 sg_df_block->Norm2(&norm_df);
193 utils.out() << "\t Jacobi relative residual norm at iteration "
194 << iter <<" is " << norm_df/norm_f << std::endl;
195 } //End of iter loop
196
197 //return status;
198 return true;
199}
200
201bool NOX::Epetra::LinearSystemSGJacobi::
202applyRightPreconditioning(bool useTranspose,
203 Teuchos::ParameterList& params,
204 const NOX::Epetra::Vector& input,
205 NOX::Epetra::Vector& result) const
206{
207 return false;
208}
209
210Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemSGJacobi::
211getScaling()
212{
213 return scaling;
214}
215
216void NOX::Epetra::LinearSystemSGJacobi::
217resetScaling(const Teuchos::RCP<NOX::Epetra::Scaling>& s)
218{
219 scaling = s;
220 return;
221}
222
223bool NOX::Epetra::LinearSystemSGJacobi::
224computeJacobian(const NOX::Epetra::Vector& x)
225{
226 bool success = jacInterfacePtr->computeJacobian(x.getEpetraVector(),
227 *sg_op);
228 sg_poly = sg_op->getSGPolynomial();
229 det_solver->setJacobianOperatorForSolve(sg_poly->getCoeffPtr(0));
230 mat_free_op->setupOperator(sg_poly);
231 return success;
232}
233
234bool NOX::Epetra::LinearSystemSGJacobi::
235createPreconditioner(const NOX::Epetra::Vector& x,
236 Teuchos::ParameterList& p,
237 bool recomputeGraph) const
238{
239 EpetraExt::BlockVector sg_x_block(View, *base_map, x.getEpetraVector());
240 bool success =
241 det_solver->createPreconditioner(*(sg_x_block.GetBlock(0)),
242 p.sublist("Deterministic Solver Parameters"),
243 recomputeGraph);
244
245 return success;
246}
247
248bool NOX::Epetra::LinearSystemSGJacobi::
249destroyPreconditioner() const
250{
251 return det_solver->destroyPreconditioner();
252}
253
254bool NOX::Epetra::LinearSystemSGJacobi::
255recomputePreconditioner(const NOX::Epetra::Vector& x,
256 Teuchos::ParameterList& linearSolverParams) const
257{
258 EpetraExt::BlockVector sg_x_block(View, *base_map, x.getEpetraVector());
259 bool success =
260 det_solver->recomputePreconditioner(
261 *(sg_x_block.GetBlock(0)),
262 linearSolverParams.sublist("Deterministic Solver Parameters"));
263
264 return success;
265}
266
267NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
268NOX::Epetra::LinearSystemSGJacobi::
269getPreconditionerPolicy(bool advanceReuseCounter)
270{
271 return det_solver->getPreconditionerPolicy(advanceReuseCounter);
272}
273
274bool NOX::Epetra::LinearSystemSGJacobi::
275isPreconditionerConstructed() const
276{
277 return det_solver->isPreconditionerConstructed();
278}
279
280bool NOX::Epetra::LinearSystemSGJacobi::
281hasPreconditioner() const
282{
283 return det_solver->hasPreconditioner();
284}
285
286Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
287getJacobianOperator() const
288{
289 return sg_op;
290}
291
292Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
293getJacobianOperator()
294{
295 return sg_op;
296}
297
298Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
299getGeneratedPrecOperator() const
300{
301 return Teuchos::null;
302}
303
304Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGJacobi::
305getGeneratedPrecOperator()
306{
307 return Teuchos::null;
308}
309
310void NOX::Epetra::LinearSystemSGJacobi::
311setJacobianOperatorForSolve(const
312 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
313{
314 Teuchos::RCP<const Stokhos::SGOperator> const_sg_op =
315 Teuchos::rcp_dynamic_cast<const Stokhos::SGOperator>(solveJacOp, true);
316 sg_op = Teuchos::rcp_const_cast<Stokhos::SGOperator>(const_sg_op);
317 sg_poly = sg_op->getSGPolynomial();
318}
319
320void NOX::Epetra::LinearSystemSGJacobi::
321setPrecOperatorForSolve(const
322 Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
323{
324}
expr expr expr dx(i, j)
Abstract base class for multivariate orthogonal polynomials.
Factory for generating stochastic Galerkin preconditioners.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition csr_vector.h:265