Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
NOX_Epetra_LinearSystem_MPBD.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_Map.h"
48#include "Epetra_Vector.h"
49#include "NOX_Epetra_Interface_Jacobian.H"
50#include "Teuchos_ParameterList.hpp"
51#include "Teuchos_TimeMonitor.hpp"
52#include "Teuchos_Assert.hpp"
53
54NOX::Epetra::LinearSystemMPBD::
55LinearSystemMPBD(
56 Teuchos::ParameterList& printingParams,
57 Teuchos::ParameterList& linearSolverParams,
58 const Teuchos::RCP<NOX::Epetra::LinearSystem>& block_solver_,
59 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
60 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
61 const Teuchos::RCP<Epetra_Operator>& J,
62 const Teuchos::RCP<const Epetra_Map>& base_map_,
63 const Teuchos::RCP<NOX::Epetra::Scaling> s):
64 block_solver(block_solver_),
65 jacInterfacePtr(iJac),
66 base_map(base_map_),
67 scaling(s),
68 utils(printingParams)
69{
70 mp_op = Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(J, true);
71 block_ops = mp_op->getMPOps();
72 num_mp_blocks = block_ops->size();
73
74 std::string prec_strategy = linearSolverParams.get("Preconditioner Strategy",
75 "Standard");
76 if (prec_strategy == "Standard")
77 precStrategy = STANDARD;
78 else if (prec_strategy == "Mean")
79 precStrategy = MEAN;
80 else if (prec_strategy == "On the fly")
81 precStrategy = ON_THE_FLY;
82 else
83 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
84 "Invalid preconditioner strategy " << prec_strategy);
85
86 if (precStrategy == STANDARD)
87 precs.resize(num_mp_blocks);
88}
89
90NOX::Epetra::LinearSystemMPBD::
91~LinearSystemMPBD()
92{
93}
94
95bool NOX::Epetra::LinearSystemMPBD::
96applyJacobian(const NOX::Epetra::Vector& input,
97 NOX::Epetra::Vector& result) const
98{
99 mp_op->SetUseTranspose(false);
100 int status = mp_op->Apply(input.getEpetraVector(), result.getEpetraVector());
101
102 return (status == 0);
103}
104
105bool NOX::Epetra::LinearSystemMPBD::
106applyJacobianTranspose(const NOX::Epetra::Vector& input,
107 NOX::Epetra::Vector& result) const
108{
109 mp_op->SetUseTranspose(true);
110 int status = mp_op->Apply(input.getEpetraVector(), result.getEpetraVector());
111 mp_op->SetUseTranspose(false);
112
113 return (status == 0);
114}
115
116bool NOX::Epetra::LinearSystemMPBD::
117applyJacobianInverse(Teuchos::ParameterList &params,
118 const NOX::Epetra::Vector &input,
119 NOX::Epetra::Vector &result)
120{
121 TEUCHOS_FUNC_TIME_MONITOR("Total deterministic solve Time");
122
123 // Extract blocks
124 EpetraExt::BlockVector input_block(View, *base_map,
125 input.getEpetraVector());
126 EpetraExt::BlockVector result_block(View, *base_map,
127 result.getEpetraVector());
128 result_block.PutScalar(0.0);
129
130
131 Teuchos::ParameterList& block_solver_params =
132 params.sublist("Deterministic Solver Parameters");
133
134 // Solve block linear systems
135 bool final_status = true;
136 bool status;
137 for (int i=0; i<num_mp_blocks; i++) {
138 NOX::Epetra::Vector nox_input(input_block.GetBlock(i),
139 NOX::Epetra::Vector::CreateView);
140 NOX::Epetra::Vector nox_result(result_block.GetBlock(i),
141 NOX::Epetra::Vector::CreateView);
142
143 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
144
145 if (precStrategy == STANDARD)
146 block_solver->setPrecOperatorForSolve(precs[i]);
147 else if (precStrategy == ON_THE_FLY) {
148 block_solver->createPreconditioner(*(prec_x->GetBlock(i)),
149 block_solver_params, false);
150 }
151
152 status = block_solver->applyJacobianInverse(block_solver_params, nox_input,
153 nox_result);
154 final_status = final_status && status;
155 }
156
157 return final_status;
158}
159
160bool NOX::Epetra::LinearSystemMPBD::
161applyRightPreconditioning(bool useTranspose,
162 Teuchos::ParameterList& params,
163 const NOX::Epetra::Vector& input,
164 NOX::Epetra::Vector& result) const
165{
166 return false;
167}
168
169Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemMPBD::
170getScaling()
171{
172 return scaling;
173}
174
175void NOX::Epetra::LinearSystemMPBD::
176resetScaling(const Teuchos::RCP<NOX::Epetra::Scaling>& s)
177{
178 scaling = s;
179 return;
180}
181
182bool NOX::Epetra::LinearSystemMPBD::
183computeJacobian(const NOX::Epetra::Vector& x)
184{
185 bool success = jacInterfacePtr->computeJacobian(x.getEpetraVector(),
186 *mp_op);
187 block_ops = mp_op->getMPOps();
188 //block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
189 return success;
190}
191
192bool NOX::Epetra::LinearSystemMPBD::
193createPreconditioner(const NOX::Epetra::Vector& x,
194 Teuchos::ParameterList& p,
195 bool recomputeGraph) const
196{
197 EpetraExt::BlockVector mp_x_block(View, *base_map, x.getEpetraVector());
198 Teuchos::ParameterList& solverParams =
199 p.sublist("Deterministic Solver Parameters");
200 bool total_success = true;
201 if (precStrategy == STANDARD) {
202 for (int i=0; i<num_mp_blocks; i++) {
203 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
204 if (precs[i] != Teuchos::null)
205 block_solver->setPrecOperatorForSolve(precs[i]);
206 bool success =
207 block_solver->createPreconditioner(*(mp_x_block.GetBlock(i)),
208 solverParams, recomputeGraph);
209 precs[i] = block_solver->getGeneratedPrecOperator();
210 total_success = total_success && success;
211 }
212 }
213
214 else if (precStrategy == MEAN) {
215 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
216 bool success =
217 block_solver->createPreconditioner(*(mp_x_block.GetBlock(0)),
218 solverParams, recomputeGraph);
219 total_success = total_success && success;
220 }
221
222 else if (precStrategy == ON_THE_FLY) {
223 if (prec_x == Teuchos::null)
224 prec_x = Teuchos::rcp(new EpetraExt::BlockVector(mp_x_block));
225 else
226 *prec_x = mp_x_block;
227 }
228
229 return total_success;
230}
231
232bool NOX::Epetra::LinearSystemMPBD::
233destroyPreconditioner() const
234{
235 return block_solver->destroyPreconditioner();
236}
237
238bool NOX::Epetra::LinearSystemMPBD::
239recomputePreconditioner(const NOX::Epetra::Vector& x,
240 Teuchos::ParameterList& p) const
241{
242 EpetraExt::BlockVector mp_x_block(View, *base_map, x.getEpetraVector());
243 Teuchos::ParameterList& solverParams =
244 p.sublist("Deterministic Solver Parameters");
245 bool total_success = true;
246
247 if (precStrategy == STANDARD) {
248 for (int i=0; i<num_mp_blocks; i++) {
249 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
250 if (precs[i] != Teuchos::null)
251 block_solver->setPrecOperatorForSolve(precs[i]);
252 bool success =
253 block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(i)),
254 solverParams);
255 precs[i] = block_solver->getGeneratedPrecOperator();
256 total_success = total_success && success;
257 }
258 }
259
260 else if (precStrategy == MEAN) {
261 block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
262 bool success =
263 block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(0)),
264 solverParams);
265 total_success = total_success && success;
266 }
267
268 else if (precStrategy == ON_THE_FLY) {
269 if (prec_x == Teuchos::null)
270 prec_x = Teuchos::rcp(new EpetraExt::BlockVector(mp_x_block));
271 else
272 *prec_x = mp_x_block;
273 }
274
275 return total_success;
276}
277
278NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
279NOX::Epetra::LinearSystemMPBD::
280getPreconditionerPolicy(bool advanceReuseCounter)
281{
282 return block_solver->getPreconditionerPolicy(advanceReuseCounter);
283}
284
285bool NOX::Epetra::LinearSystemMPBD::
286isPreconditionerConstructed() const
287{
288 return block_solver->isPreconditionerConstructed();
289}
290
291bool NOX::Epetra::LinearSystemMPBD::
292hasPreconditioner() const
293{
294 return block_solver->hasPreconditioner();
295}
296
297Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
298getJacobianOperator() const
299{
300 return mp_op;
301}
302
303Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
304getJacobianOperator()
305{
306 return mp_op;
307}
308
309Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
310getGeneratedPrecOperator() const
311{
312 if (precStrategy == MEAN)
313 return block_solver->getGeneratedPrecOperator();
314 else
315 TEUCHOS_TEST_FOR_EXCEPTION(
316 true, std::logic_error,
317 "Cannot call getGeneratedPrecOperator() unless prec strategy is Mean");
318 return Teuchos::null;
319}
320
321Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemMPBD::
322getGeneratedPrecOperator()
323{
324 if (precStrategy == MEAN)
325 return block_solver->getGeneratedPrecOperator();
326 else
327 TEUCHOS_TEST_FOR_EXCEPTION(
328 true, std::logic_error,
329 "Cannot call getGeneratedPrecOperator() unless prec strategy is Mean");
330 return Teuchos::null;
331}
332
333void NOX::Epetra::LinearSystemMPBD::
334setJacobianOperatorForSolve(const
335 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
336{
337 Teuchos::RCP<const Stokhos::BlockDiagonalOperator> const_mp_op =
338 Teuchos::rcp_dynamic_cast<const Stokhos::BlockDiagonalOperator>(solveJacOp,
339 true);
340 mp_op = Teuchos::rcp_const_cast<Stokhos::BlockDiagonalOperator>(const_mp_op);
341 block_ops = mp_op->getMPOps();
342}
343
344void NOX::Epetra::LinearSystemMPBD::
345setPrecOperatorForSolve(const Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
346{
347 if (precStrategy == MEAN)
348 block_solver->setPrecOperatorForSolve(solvePrecOp);
349 else
350 TEUCHOS_TEST_FOR_EXCEPTION(
351 true, std::logic_error,
352 "Cannot call setPrecOperatorForSolve() unless prec strategy is Mean");
353}
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265