ShyLU  Version of the Day
IQRSolver.cpp
1 //@HEADER
2 // ************************************************************************
3 //
4 // ShyLU: Hybrid preconditioner package
5 // Copyright 2012 Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ************************************************************************
40 //@HEADER
41 
50 #include <IQRSolver.h>
51 
52 #include <cmath>
53 
54 #include <Ifpack_config.h>
55 
56 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
57 #include "Ifpack_DynamicFactory.h"
58 #else
59 #include "Ifpack.h"
60 #endif
61 
62 namespace IQR {
63 
64 IQRSolver::IQRSolver(const Teuchos::ParameterList& pList)
65  : gmresStateManager_(),
66  prec_(0),
67  isComputed_(false),
68  krylovDim_(0.0),
69  numIter_(0),
70  doScaling_(false),
71  useFullIQR_(true),
72  precType_("Amesos"),
73  precAmesosType_("Amesos_Klu")
74 {
75  // We shouldn't modify the external parameter list
76  Teuchos::ParameterList params = pList;
77 
78  // Default mode, G or IQR?
79  useFullIQR_ = params.get<bool>("Use full IQR", true);
80 
81  // Default parameters for IQR
82  krylovDim_ = params.get<double>("IQR Krylov Dim", 0.5);
83  numIter_ = params.get<int>("IQR Number Iterations", 0);
84  doScaling_ = params.get<bool>("IQR Scaling", false);
85  precType_ = params.get<string>("IQR Initial Prec Type", "Amesos");
86  precAmesosType_ = params.get<string>("IQR Initial Prec Amesos Type", "Amesos_Klu");
87 }
88 
89 IQRSolver::~IQRSolver()
90 {
91 }
92 
93 int IQRSolver::Compute(const ShyLU_Probing_Operator& S,
94  const Epetra_MultiVector& B,
95  Epetra_MultiVector & X)
96 {
97  // The preconditioner is used regardless of the value ofuseFullIQR_
98  Epetra_CrsMatrix* G = S.G_;
99 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
100  Ifpack_DynamicFactory Factory;
101 #else
102  Ifpack Factory;
103 #endif
104  prec_ = Teuchos::rcp(Factory.Create(precType_, G, 1));
105  Teuchos::ParameterList pList;
106  pList.set("amesos: solver type", precAmesosType_, "");
107  pList.set("Reindex", true);
108  prec_->SetParameters(pList);
109 
110  IFPACK_CHK_ERR(prec_->Initialize());
111  IFPACK_CHK_ERR(prec_->Compute());
112 
113  // Only for IQR
114  if (useFullIQR_) {
115  // Computation phase - only during the first outer GMRES iteration
116  int sSize = S.OperatorDomainMap().NumGlobalElements();
117  int kSize = std::floor(krylovDim_ * sSize);
118 
119  gmresStateManager_ = Teuchos::rcp(
120  new GMRESStateManager(S.OperatorDomainMap(),
121  kSize, false, doScaling_));
122  gmresStateManager_->isFirst = true;
123 
124  double tol = 1e-10;
126  IQR::GMRES<Epetra_Operator,
127  Epetra_MultiVector,
129  Ifpack_Preconditioner,
130  GMRESStateManager,
131  std::vector<double>, double>
132  (S, X, B, &L, &(*prec_), *(gmresStateManager_), kSize, tol);
133  gmresStateManager_->P2 = &(*prec_);
134 
135 // if (! Xs.Comm().MyPID()) {
136 // std::cout << "KSIZE: " << kSize
137 // << ", SSIZE: " << sSize
138 // << ", TOL: " << tol << std::endl;
139 // }
140  }
141  isComputed_ = true;
142 
143  return 0;
144 }
145 
146 int IQRSolver::Solve(const ShyLU_Probing_Operator& S,
147  const Epetra_MultiVector& B,
148  Epetra_MultiVector & X)
149 {
150  if (! isComputed_) {
151  IFPACK_CHK_ERR(Compute(S, B, X));
152  }
153  // Solve phase
154  if (! useFullIQR_) {
155  IFPACK_CHK_ERR(prec_->ApplyInverse(B, X));
156  } else {
157  if (numIter_ > 0) {
158  GMRESStateManager newGmresManager(S.OperatorDomainMap(),
159  numIter_+1, false);
160  double tol=1e-10;
162  IQR::GMRES<Epetra_Operator, Epetra_MultiVector,
163  IQR::IdPreconditioner, GMRESStateManager, GMRESStateManager,
164  std::vector<double>, double>
165  (S, X, B, &L, &*(gmresStateManager_), newGmresManager, numIter_, tol);
166  } else {
167  IFPACK_CHK_ERR(gmresStateManager_->ApplyInverse(B, X));
168  }
169  }
170 
171  return 0;
172 }
173 
174 } /* namespace IQR */
int Compute()
Compute ILU factors L and U using the specified parameters.
Encapsulates the IQR inexact solver functionality.