ShyLU  Version of the Day
Ifpack_ShyLU.cpp
Go to the documentation of this file.
1 
2 //@HEADER
3 // ************************************************************************
4 //
5 // ShyLU: Hybrid preconditioner package
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ************************************************************************
41 //@HEADER
42 
51 #include "Ifpack_ShyLU.h"
52 #include "Teuchos_Time.hpp"
53 #include "Ifpack_config.h"
54 
55 #include <IQRSolver.h>
56 
57 //#define DUMP_MATRICES
58 
59 Ifpack_ShyLU::Ifpack_ShyLU(Epetra_CrsMatrix* A):
60  A_(A),
61  IsParallel_(true),
62  IsInitialized_(false),
63  IsComputed_(false),
64  Label_("Ifpack_ShyLU"),
65  NumCompute_(0),
66  NumApplyInverse_(0),
67  Time_(A_->Comm())
68 {
69 }
70 
71 void Ifpack_ShyLU::Destroy()
72 {
73  if (IsInitialized_)
74  {
75  //delete A_;
76  //delete partitioner_;
77  //delete rd_;
78 
79  // I would rather explicitly delete
80  /*delete slu_sym_.LP;
81  delete slu_sym_.Solver;
82  delete slu_sym_.C;
83  delete slu_sym_.R;
84  delete slu_sym_.D;
85  delete slu_sym_.G;
86  delete slu_sym_.Sg;
87  delete slu_sym_.prober;*/
88 
89  if (slu_config_.schurSolver == "Amesos")
90  {
91  delete slu_data_.dsolver;
92  }
93  else if (slu_config_.libName == "Belos")
94  delete slu_data_.innersolver;
95 
96  delete[] slu_data_.DRowElems;
97  delete[] slu_data_.SRowElems;
98  delete[] slu_data_.DColElems;
99  delete[] slu_data_.gvals;
100  }
101  if (IsComputed_)
102  {
103  // I would rather explicitly delete
104  // delete slu_data_.Sbar;
105  slu_data_.localSbargraph = Teuchos::null;
106  slu_data_.guided_prober = Teuchos::null;
107  }
108 }
109 
111 {
112  if(Comm().NumProc() != 1)
113  IsParallel_ = true;
114  else
115  IsParallel_ = false;
116 
117  // TODO:
118  // Need to enable partitioning here in Initialize once moving to Belos
119 
120  // Cannot call this method , need the partitioner around TODO : Can we
121  // avoid this
122  //A_ = balanceAndRedistribute(A_, List_);
123 
124  // ==================== Symbolic factorization =========================
125  // 1. Partition and redistribute [
126  //partitioner_ = new Isorropia::Epetra::Partitioner(A_, List_, false);
127  //partitioner_->partition();
128 
129  //rd_ = new Isorropia::Epetra::Redistributor(partitioner_);
130  //Epetra_CrsMatrix *newA;
131  //rd_->redistribute(*A_, newA);
132  //A_ = newA;
133  // ]
134 
135  /*double Sdiagfactor = 0.05; // hard code the diagonals
136  shylu_factor(A_, 1, LP_, Solver_, C_, Dnr_, DRowElems_, Snr_, SRowElems_,
137  Sbar_, Sdiagfactor);*/
138  slu_config_.sym = Teuchos::getParameter<int>(List_,
139  "Symmetry");
140  slu_config_.libName = Teuchos::getParameter<string>(List_,
141  "Outer Solver Library");
142  string schurApproxMethod = Teuchos::getParameter<string>(List_,
143  "Schur Approximation Method");
144  slu_config_.schurSolver = Teuchos::getParameter<string>(List_,
145  "Schur Complement Solver");
146  slu_config_.schurAmesosSolver = List_.get<string>("Schur Amesos Solver",
147  "Amesos_Klu");
148  //slu_config_.diagonalBlockSolver =
149  //Teuchos::getParameter<string>(List_, "Diagonal Block Solver");
150  slu_config_.diagonalBlockSolver = List_.get<string>("Diagonal Block Solver",
151  "Amesos_Klu");
152 
153 
154  slu_config_.relative_threshold = 0.0;
155  slu_config_.Sdiagfactor = 0.05;
156  if (schurApproxMethod == "A22AndBlockDiagonals")
157  {
158  slu_config_.schurApproxMethod = 1;
159  slu_config_.Sdiagfactor = Teuchos::getParameter<double>(List_,
160  "Diagonal Factor");
161  }
162  else if (schurApproxMethod == "Threshold")
163  {
164  slu_config_.schurApproxMethod = 2;
165  slu_config_.relative_threshold = Teuchos::getParameter<double>(List_,
166  "Relative Threshold");
167  }
168  else if (schurApproxMethod == "Guided Probing")
169  {
170  slu_config_.schurApproxMethod = 3;
171  slu_config_.reset_iter = List_.get<int>("Schur Recompute Iteration", 10);
172  slu_config_.relative_threshold = Teuchos::getParameter<double>(List_,
173  "Relative Threshold");
174  }
175  if (schurApproxMethod == "IQR")
176  {
177  slu_config_.schurSolver = "IQR";
178  slu_config_.schurApproxMethod = 4;
179 
180  List_.set<bool>("Use full IQR", true, "");
181  slu_data_.iqrSolver.reset(new IQR::IQRSolver(List_));
182  }
183  if (schurApproxMethod == "G")
184  {
185  slu_config_.schurSolver = "G";
186  slu_config_.schurApproxMethod = 5;
187 
188  List_.set<bool>("Use full IQR", false, "");
189  slu_data_.iqrSolver.reset(new IQR::IQRSolver(List_));
190  }
191 
192  slu_config_.schurPreconditioner = List_.get<string>("Schur Preconditioner",
193  "ILU stand-alone");
194  slu_config_.silent_subiter = List_.get<bool>("Silent subiterations",
195  true);
196 
197  slu_config_.inner_tolerance = Teuchos::getParameter<double>(List_,
198  "Inner Solver Tolerance");
199  slu_config_.inner_maxiters = Teuchos::getParameter<int>(List_,
200  "Inner Solver MaxIters");
201  slu_config_.overlap = List_.get<int>("Schur Preconditioner Overlap", 0);
202  string sep_type = Teuchos::getParameter<string>(List_,
203  "Separator Type");
204  // mfh 25 May 2015: dl is unused (don't declare it, or you'll get
205  // build warnings). However, two-argument get() has side effects,
206  // so I don't want to comment it out.
207  //
208  //int dl = List_.get<int>("Debug Level", 0);
209  (void) List_.get<int>("Debug Level", 0);
210  //slu_config_.dm.setDebugLevel(dl);
211 
212  if (sep_type == "Wide")
213  slu_config_.sep_type = 1;
214  else
215  slu_config_.sep_type = 2;
216 
217  shylu_symbolic_factor(A_, &slu_sym_, &slu_data_, &slu_config_);
218 
219  if (slu_config_.schurSolver == "Amesos")
220  {
221  slu_data_.LP2 = Teuchos::rcp(new Epetra_LinearProblem());
222  }
223  else
224  {
225  if (slu_config_.libName == "Belos")
226  slu_data_.innersolver = new AztecOO() ;
227  else
228  slu_data_.innersolver = NULL;
229  }
230 
231  IsInitialized_ = true;
232  return 0;
233 }
234 
235 int Ifpack_ShyLU::SetParameters(Teuchos::ParameterList& parameterlist)
236 {
237 #ifdef HAVE_IFPACK_DYNAMIC_FACTORY
238  List_ = parameterlist.sublist("ShyLU list", true);
239 #else
240  List_ = parameterlist;
241 #endif
242  return 0;
243 }
244 
246 {
247  Teuchos::Time ftime("setup time");
248  ftime.start();
249 
250  slu_data_.num_compute = NumCompute_;
251 
252  shylu_factor(A_, &slu_sym_, &slu_data_, &slu_config_);
253 
254  ftime.stop();
255  IsComputed_ = true;
256  NumCompute_++;
257  //cout << " Done with the compute" << endl ;
258  return 0;
259 }
260 
262 {
263  // Dummy function, To show the error in AztecOO, This works
264  //cout << "Entering JustTryIt" << endl;
265  AztecOO *solver;
266  solver = slu_data_.innersolver;
267  //cout << solver_ << endl;
268  Epetra_Map BsMap(-1, slu_data_.Snr, slu_data_.SRowElems, 0, A_->Comm());
269  Epetra_MultiVector Xs(BsMap, 1);
270  Epetra_MultiVector Bs(BsMap, 1);
271  Xs.PutScalar(0.0);
272  solver->SetLHS(&Xs);
273  solver->SetRHS(&Bs);
274  solver->Iterate(30, 1e-10);
275  return 0;
276 }
277 
278 int Ifpack_ShyLU::ApplyInverse(const Epetra_MultiVector& X,
279  Epetra_MultiVector& Y) const
280 {
281 #ifdef DUMP_MATRICES
282  if (NumApplyInverse_ == 0)
283  {
284  EpetraExt::MultiVectorToMatlabFile("X.mat", X);
285  }
286 #endif
287  //cout << "Entering ApplyInvers" << endl;
288 
289  shylu_solve(&slu_sym_, &slu_data_, &slu_config_, X, Y);
290 #ifdef DUMP_MATRICES
291  if (NumApplyInverse_ == 0)
292  {
293  EpetraExt::MultiVectorToMatlabFile("Y.mat", Y);
294  }
295 #endif
296  NumApplyInverse_++;
297  //cout << "Leaving ApplyInvers" << endl;
298  return 0;
299 }
300 
302 double Ifpack_ShyLU::Condest(const Ifpack_CondestType CT,
303  const int MaxIters, const double Tol, Epetra_RowMatrix* Matrix_in)
304 {
305  return -1.0;
306 }
307 
309 ostream& Ifpack_ShyLU::Print(ostream& os) const
310 {
311  os << " !!!!!!!!! " << endl;
312  return os;
313 }
int shylu_solve(shylu_symbolic *ssym, shylu_data *data, shylu_config *config, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
Call solve on multiple RHS.
int Compute()
Compute ILU factors L and U using the specified parameters.
int ApplyInverse(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const
Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y...
Encapsulates the IQR inexact solver functionality.
Use ShyLU as a preconditioner within IFPACK.
double Condest(const Ifpack_CondestType CT=Ifpack_Cheap, const int MaxIters=1550, const double Tol=1e-9, Epetra_RowMatrix *Matrix_in=0)
Computes the estimated condition number and returns the value.
int SetParameters(Teuchos::ParameterList &parameterlist)
Set parameters using a Teuchos::ParameterList object.
const Epetra_Comm & Comm() const
Returns the Epetra_BlockMap object associated with the range of this matrix operator.
Definition: Ifpack_ShyLU.h:233
int shylu_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Main function call into ShylU.
int JustTryIt()
Returns the computed estimated condition number, or -1.0 if not computed.
int shylu_symbolic_factor(Epetra_CrsMatrix *A, shylu_symbolic *ssym, shylu_data *data, shylu_config *config)
Call symbolic factorization on matrix.
int Initialize()
Initialize the preconditioner, does not touch matrix values.
Ifpack_ShyLU(Epetra_CrsMatrix *A)
Constructor.
virtual ostream & Print(ostream &os) const
Prints on stream basic information about this object.