Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
NOX_Epetra_LinearSystem_SGGS.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 "EpetraExt_BlockVector.h"
51#include "EpetraExt_BlockUtility.h"
52#include "Teuchos_ParameterList.hpp"
53#include "Teuchos_TimeMonitor.hpp"
54
55NOX::Epetra::LinearSystemSGGS::
56LinearSystemSGGS(
57 Teuchos::ParameterList& printingParams,
58 Teuchos::ParameterList& linearSolverParams,
59 const Teuchos::RCP<NOX::Epetra::LinearSystem>& det_solver_,
60 const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq,
61 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac,
62 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
63 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data,
64 const Teuchos::RCP<Epetra_Operator>& J,
65 const Teuchos::RCP<const Epetra_Map>& base_map_,
66 const Teuchos::RCP<const Epetra_Map>& sg_map_,
67 const Teuchos::RCP<NOX::Epetra::Scaling> s):
68 det_solver(det_solver_),
69 sg_basis(sg_basis_),
70 epetraCijk(sg_parallel_data->getEpetraCijk()),
71 is_stoch_parallel(epetraCijk->isStochasticParallel()),
72 stoch_row_map(epetraCijk->getStochasticRowMap()),
73 Cijk(epetraCijk->getParallelCijk()),
74 jacInterfacePtr(iJac),
75 base_map(base_map_),
76 sg_map(sg_map_),
77 scaling(s),
78 utils(printingParams),
79 is_parallel(epetraCijk->isStochasticParallel())
80{
81 sg_op = Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(J, true);
82 sg_poly = sg_op->getSGPolynomial();
83
84 sg_df_block = Teuchos::rcp(new EpetraExt::BlockVector(*base_map, *sg_map));
85 sg_y_block = Teuchos::rcp(new EpetraExt::BlockVector(*base_map, *sg_map));
86 kx = Teuchos::rcp(new Epetra_Vector(*base_map));
87
88 only_use_linear = linearSolverParams.get("Only Use Linear Terms", false);
89 k_limit = sg_poly->size();
90 int dim = sg_basis->dimension();
91 if (only_use_linear && sg_poly->size() > dim+1)
92 k_limit = dim + 1;
93
94 if (is_parallel) {
95 Teuchos::RCP<const Epetra_BlockMap> stoch_col_map =
96 epetraCijk->getStochasticColMap();
97 sg_col_map =
98 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(*base_map,
99 *stoch_col_map,
100 sg_map->Comm()));
101 col_exporter = Teuchos::rcp(new Epetra_Export(*sg_col_map, *sg_map));
102 sg_df_col = Teuchos::rcp(new EpetraExt::BlockVector(*base_map,
103 *sg_col_map));
104 sg_df_tmp = Teuchos::rcp(new EpetraExt::BlockVector(*base_map, *sg_map));
105 }
106}
107
108NOX::Epetra::LinearSystemSGGS::
109~LinearSystemSGGS()
110{
111}
112
113bool NOX::Epetra::LinearSystemSGGS::
114applyJacobian(const NOX::Epetra::Vector& input,
115 NOX::Epetra::Vector& result) const
116{
117 sg_op->SetUseTranspose(false);
118 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
119
120 return (status == 0);
121}
122
123bool NOX::Epetra::LinearSystemSGGS::
124applyJacobianTranspose(const NOX::Epetra::Vector& input,
125 NOX::Epetra::Vector& result) const
126{
127 sg_op->SetUseTranspose(true);
128 int status = sg_op->Apply(input.getEpetraVector(), result.getEpetraVector());
129 sg_op->SetUseTranspose(false);
130
131 return (status == 0);
132}
133
134bool NOX::Epetra::LinearSystemSGGS::
135applyJacobianInverse(Teuchos::ParameterList &params,
136 const NOX::Epetra::Vector &input,
137 NOX::Epetra::Vector &result)
138{
139 int max_iter = params.get("Max Iterations",100 );
140 double sg_tol = params.get("Tolerance", 1e-12);
141 bool scale_op = params.get("Scale Operator by Inverse Basis Norms", true);
142
143 // Extract blocks
144 EpetraExt::BlockVector sg_dx_block(View, *base_map, result.getEpetraVector());
145 EpetraExt::BlockVector sg_f_block(View, *base_map, input.getEpetraVector());
146 sg_dx_block.PutScalar(0.0);
147
148 // Compute initial residual norm
149 double norm_f,norm_df;
150 sg_f_block.Norm2(&norm_f);
151 sg_op->Apply(sg_dx_block, *sg_df_block);
152 sg_df_block->Update(-1.0, sg_f_block, 1.0);
153 sg_df_block->Norm2(&norm_df);
154
155 Teuchos::RCP<Epetra_Vector> f, df, dx;
156 const Teuchos::Array<double>& norms = sg_basis->norm_squared();
157
158 // sg_df_block stores the right-hand-sides for the Gauss-Seidel iteration
159 // sg_y_block stores the residual, i.e., A*x-b
160 // sg_df_col stores off-processor contributions to the RHS and residual
161 // In parallel, this is essentially domain decomposition using Gauss-Seidel
162 // in each domain and Jacobi across domains.
163
164 int iter = 0;
165 while (((norm_df/norm_f)>sg_tol) && (iter<max_iter)) {
166 TEUCHOS_FUNC_TIME_MONITOR("Total global solve Time");
167 iter++;
168
169 sg_y_block->Update(1.0, sg_f_block, 0.0);
170 if (is_parallel)
171 sg_df_col->PutScalar(0.0);
172
173 for (Cijk_type::i_iterator i_it=Cijk->i_begin();
174 i_it!=Cijk->i_end(); ++i_it) {
175 int i = Stokhos::index(i_it);
176 f = sg_f_block.GetBlock(i);
177 df = sg_df_block->GetBlock(i);
178 dx = sg_dx_block.GetBlock(i);
179
180 dx->PutScalar(0.0);
181 Teuchos::ParameterList& det_solver_params =
182 params.sublist("Deterministic Solver Parameters");
183 NOX::Epetra::Vector nox_df(df, NOX::Epetra::Vector::CreateView);
184 NOX::Epetra::Vector nox_dx(dx, NOX::Epetra::Vector::CreateView);
185 // Solve linear system
186 {
187 TEUCHOS_FUNC_TIME_MONITOR("Total deterministic solve Time");
188 det_solver->applyJacobianInverse(det_solver_params, nox_df, nox_dx);
189 }
190
191 df->Update(1.0, *f, 0.0);
192
193 for (Cijk_type::ik_iterator k_it = Cijk->k_begin(i_it);
194 k_it != Cijk->k_end(i_it); ++k_it) {
195 int k = Stokhos::index(k_it);
196 if (k!=0 && k<k_limit) {
197 (*sg_poly)[k].Apply(*dx, *kx);
198 for (Cijk_type::ikj_iterator j_it = Cijk->j_begin(k_it);
199 j_it != Cijk->j_end(k_it); ++j_it) {
200 int j = Stokhos::index(j_it);
201 int j_gid = epetraCijk->GCID(j);
202 double c = Stokhos::value(j_it);
203 if (scale_op)
204 c /= norms[j_gid];
205 bool owned = epetraCijk->myGRID(j_gid);
206 if (!is_parallel || owned) {
207 sg_df_block->GetBlock(j)->Update(-c, *kx, 1.0);
208 sg_y_block->GetBlock(j)->Update(-c, *kx, 1.0);
209 }
210 else
211 sg_df_col->GetBlock(j)->Update(-c, *kx, 1.0);
212 }
213 }
214 }
215
216 (*sg_poly)[0].Apply(*dx, *kx);
217 sg_y_block->GetBlock(i)->Update(-1.0, *kx, 1.0);
218
219 } //End of k loop
220
221 // Add in contributions from off-process
222 if (is_parallel) {
223 sg_df_tmp->Export(*sg_df_col, *col_exporter, InsertAdd);
224 sg_df_block->Update(1.0, *sg_df_tmp, 1.0);
225 sg_y_block->Update(1.0, *sg_df_tmp, 1.0);
226 }
227
228 sg_y_block->Norm2(&norm_df);
229 utils.out() << "\t Gauss-Seidel relative residual norm at iteration "
230 << iter <<" is " << norm_df/norm_f << std::endl;
231 } //End of iter loop
232
233 //return status;
234 return true;
235}
236
237bool NOX::Epetra::LinearSystemSGGS::
238applyRightPreconditioning(bool useTranspose,
239 Teuchos::ParameterList& params,
240 const NOX::Epetra::Vector& input,
241 NOX::Epetra::Vector& result) const
242{
243 return false;
244}
245
246Teuchos::RCP<NOX::Epetra::Scaling> NOX::Epetra::LinearSystemSGGS::
247getScaling()
248{
249 return scaling;
250}
251
252void NOX::Epetra::LinearSystemSGGS::
253resetScaling(const Teuchos::RCP<NOX::Epetra::Scaling>& s)
254{
255 scaling = s;
256 return;
257}
258
259bool NOX::Epetra::LinearSystemSGGS::
260computeJacobian(const NOX::Epetra::Vector& x)
261{
262 bool success = jacInterfacePtr->computeJacobian(x.getEpetraVector(),
263 *sg_op);
264 sg_poly = sg_op->getSGPolynomial();
265 det_solver->setJacobianOperatorForSolve(sg_poly->getCoeffPtr(0));
266 return success;
267}
268
269bool NOX::Epetra::LinearSystemSGGS::
270createPreconditioner(const NOX::Epetra::Vector& x,
271 Teuchos::ParameterList& p,
272 bool recomputeGraph) const
273{
274 EpetraExt::BlockVector sg_x_block(View, *base_map, x.getEpetraVector());
275 bool success =
276 det_solver->createPreconditioner(*(sg_x_block.GetBlock(0)),
277 p.sublist("Deterministic Solver Parameters"),
278 recomputeGraph);
279
280 return success;
281}
282
283bool NOX::Epetra::LinearSystemSGGS::
284destroyPreconditioner() const
285{
286 return det_solver->destroyPreconditioner();
287}
288
289bool NOX::Epetra::LinearSystemSGGS::
290recomputePreconditioner(const NOX::Epetra::Vector& x,
291 Teuchos::ParameterList& linearSolverParams) const
292{
293 EpetraExt::BlockVector sg_x_block(View, *base_map, x.getEpetraVector());
294 bool success =
295 det_solver->recomputePreconditioner(
296 *(sg_x_block.GetBlock(0)),
297 linearSolverParams.sublist("Deterministic Solver Parameters"));
298
299 return success;
300
301}
302
303NOX::Epetra::LinearSystem::PreconditionerReusePolicyType
304NOX::Epetra::LinearSystemSGGS::
305getPreconditionerPolicy(bool advanceReuseCounter)
306{
307 return det_solver->getPreconditionerPolicy(advanceReuseCounter);
308}
309
310bool NOX::Epetra::LinearSystemSGGS::
311isPreconditionerConstructed() const
312{
313 return det_solver->isPreconditionerConstructed();
314}
315
316bool NOX::Epetra::LinearSystemSGGS::
317hasPreconditioner() const
318{
319 return det_solver->hasPreconditioner();
320}
321
322Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
323getJacobianOperator() const
324{
325 return sg_op;
326}
327
328Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
329getJacobianOperator()
330{
331 return sg_op;
332}
333
334Teuchos::RCP<const Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
335getGeneratedPrecOperator() const
336{
337 return Teuchos::null;
338}
339
340Teuchos::RCP<Epetra_Operator> NOX::Epetra::LinearSystemSGGS::
341getGeneratedPrecOperator()
342{
343 return Teuchos::null;
344}
345
346void NOX::Epetra::LinearSystemSGGS::
347setJacobianOperatorForSolve(const
348 Teuchos::RCP<const Epetra_Operator>& solveJacOp)
349{
350 Teuchos::RCP<const Stokhos::SGOperator> const_sg_op =
351 Teuchos::rcp_dynamic_cast<const Stokhos::SGOperator>(solveJacOp, true);
352 sg_op = Teuchos::rcp_const_cast<Stokhos::SGOperator>(const_sg_op);
353 sg_poly = sg_op->getSGPolynomial();
354}
355
356void NOX::Epetra::LinearSystemSGGS::
357setPrecOperatorForSolve(const Teuchos::RCP<const Epetra_Operator>& solvePrecOp)
358{
359}
expr expr expr dx(i, j)
Abstract base class for multivariate orthogonal polynomials.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265