ShyLU  Version of the Day
shylu_solve.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 "shylu_util.h"
52 #include "shylu.h"
53 
54 #include <Ifpack.h>
55 
56 static int shylu_dist_solve(
57  shylu_symbolic *ssym,
58  shylu_data *data,
59  shylu_config *config,
60  const Epetra_MultiVector& X,
61  Epetra_MultiVector& Y
62 )
63 {
64  int err;
65  AztecOO *solver = 0;
66  assert(X.Map().SameAs(Y.Map()));
67  //assert(X.Map().SameAs(A_->RowMap()));
68  const Epetra_MultiVector *newX;
69  newX = &X;
70  //rd_->redistribute(X, newX);
71 
72  int nvectors = newX->NumVectors();
73 
74  // May have to use importer/exporter
75  Epetra_Map BsMap(-1, data->Snr, data->SRowElems, 0, X.Comm());
76  Epetra_Map BdMap(-1, data->Dnr, data->DRowElems, 0, X.Comm());
77 
78  Epetra_MultiVector Bs(BsMap, nvectors);
79  Epetra_Import BsImporter(BsMap, newX->Map());
80 
81  assert(BsImporter.SourceMap().SameAs(newX->Map()));
82  assert((newX->Map()).SameAs(BsImporter.SourceMap()));
83 
84  Bs.Import(*newX, BsImporter, Insert);
85  Epetra_MultiVector Xs(BsMap, nvectors);
86 
87  Epetra_SerialComm LComm; // Use Serial Comm for the local vectors.
88  Epetra_Map LocalBdMap(-1, data->Dnr, data->DRowElems, 0, LComm);
89  Epetra_MultiVector localrhs(LocalBdMap, nvectors);
90  Epetra_MultiVector locallhs(LocalBdMap, nvectors);
91 
92  Epetra_MultiVector Z(BdMap, nvectors);
93 
94  Epetra_MultiVector Bd(BdMap, nvectors);
95  Epetra_Import BdImporter(BdMap, newX->Map());
96  assert(BdImporter.SourceMap().SameAs(newX->Map()));
97  assert((newX->Map()).SameAs(BdImporter.SourceMap()));
98  Bd.Import(*newX, BdImporter, Insert);
99 
100  int lda;
101  double *values;
102  err = Bd.ExtractView(&values, &lda);
103  assert (err == 0);
104  int nrows = ssym->C->RowMap().NumMyElements();
105 
106  // copy to local vector //TODO: OMP ?
107  assert(lda == nrows);
108  for (int v = 0; v < nvectors; v++)
109  {
110  for (int i = 0; i < nrows; i++)
111  {
112  err = localrhs.ReplaceMyValue(i, v, values[i+v*lda]);
113  assert (err == 0);
114  }
115  }
116 
117  // TODO : Do we need to reset the lhs and rhs here ?
118  if (config->amesosForDiagonal)
119  {
120  ssym->LP->SetRHS(&localrhs);
121  ssym->LP->SetLHS(&locallhs);
122  ssym->Solver->Solve();
123  }
124  else
125  {
126  ssym->ifSolver->ApplyInverse(localrhs, locallhs);
127  }
128 
129  err = locallhs.ExtractView(&values, &lda);
130  assert (err == 0);
131 
132  // copy to distributed vector //TODO: OMP ?
133  assert(lda == nrows);
134  for (int v = 0; v < nvectors; v++)
135  {
136  for (int i = 0; i < nrows; i++)
137  {
138  err = Z.ReplaceMyValue(i, v, values[i+v*lda]);
139  assert (err == 0);
140  }
141  }
142 
143  Epetra_MultiVector temp1(BsMap, nvectors);
144  ssym->R->Multiply(false, Z, temp1);
145  Bs.Update(-1.0, temp1, 1.0);
146 
147  Xs.PutScalar(0.0);
148 
149  Epetra_LinearProblem Problem(data->Sbar.get(), &Xs, &Bs);
150  if (config->schurSolver == "Amesos")
151  {
152  Amesos_BaseSolver *solver2 = data->dsolver;
153  data->LP2->SetLHS(&Xs);
154  data->LP2->SetRHS(&Bs);
155  //cout << "Calling solve *****************************" << endl;
156  solver2->Solve();
157  //cout << "Out of solve *****************************" << endl;
158  }
159  else
160  {
161  if (config->libName == "Belos")
162  {
163  solver = data->innersolver;
164  solver->SetLHS(&Xs);
165  solver->SetRHS(&Bs);
166  }
167  else
168  {
169  // See the comment above on why we are not able to reuse the solver
170  // when outer solve is AztecOO as well.
171  solver = new AztecOO();
172  //solver.SetPrecOperator(precop_);
173  solver->SetAztecOption(AZ_solver, AZ_gmres);
174  // Do not use AZ_none
175  solver->SetAztecOption(AZ_precond, AZ_dom_decomp);
176  //solver->SetAztecOption(AZ_precond, AZ_none);
177  //solver->SetAztecOption(AZ_precond, AZ_Jacobi);
179  //solver->SetAztecOption(AZ_overlap, 3);
180  //solver->SetAztecOption(AZ_subdomain_solve, AZ_ilu);
181  //solver->SetAztecOption(AZ_output, AZ_all);
182  //solver->SetAztecOption(AZ_diagnostics, AZ_all);
183  solver->SetProblem(Problem);
184  }
185 
186  // What should be a good inner_tolerance :-) ?
187  solver->Iterate(config->inner_maxiters, config->inner_tolerance);
188  }
189 
190  Epetra_MultiVector temp(BdMap, nvectors);
191  ssym->C->Multiply(false, Xs, temp);
192  temp.Update(1.0, Bd, -1.0);
193 
194  //Epetra_SerialComm LComm; // Use Serial Comm for the local vectors.
195  //Epetra_Map LocalBdMap(-1, data->Dnr, data->DRowElems, 0, LComm);
196  //Epetra_MultiVector localrhs(LocalBdMap, nvectors);
197  //Epetra_MultiVector locallhs(LocalBdMap, nvectors);
198 
199  //int lda;
200  //double *values;
201  err = temp.ExtractView(&values, &lda);
202  assert (err == 0);
203  //int nrows = data->Cptr->RowMap().NumMyElements();
204 
205  // copy to local vector //TODO: OMP ?
206  assert(lda == nrows);
207  for (int v = 0; v < nvectors; v++)
208  {
209  for (int i = 0; i < nrows; i++)
210  {
211  err = localrhs.ReplaceMyValue(i, v, values[i+v*lda]);
212  assert (err == 0);
213  }
214  }
215 
216  if (config->amesosForDiagonal)
217  {
218  ssym->LP->SetRHS(&localrhs);
219  ssym->LP->SetLHS(&locallhs);
220  ssym->Solver->Solve();
221  }
222  else
223  {
224  ssym->ifSolver->ApplyInverse(localrhs, locallhs);
225  }
226 
227  err = locallhs.ExtractView(&values, &lda);
228  assert (err == 0);
229 
230  // copy to distributed vector //TODO: OMP ?
231  assert(lda == nrows);
232  for (int v = 0; v < nvectors; v++)
233  {
234  for (int i = 0; i < nrows; i++)
235  {
236  err = temp.ReplaceMyValue(i, v, values[i+v*lda]);
237  assert (err == 0);
238  }
239  }
240 
241  // For checking faults
242  //if (NumApplyInverse_ == 5) temp.ReplaceMyValue(0, 0, 0.0);
243 
244  Epetra_Export XdExporter(BdMap, Y.Map());
245  Y.Export(temp, XdExporter, Insert);
246 
247  Epetra_Export XsExporter(BsMap, Y.Map());
248  Y.Export(Xs, XsExporter, Insert);
249 
250  if (config->libName == "Belos" || config->schurSolver == "Amesos")
251  {
252  // clean up
253  }
254  else
255  {
256  delete solver;
257  }
258  return 0;
259 }
260 
261 static int shylu_local_solve(
262  shylu_symbolic *ssym,
263  shylu_data *data,
264  shylu_config *config,
265  const Epetra_MultiVector& X,
266  Epetra_MultiVector& Y
267 )
268 {
269  int err;
270 #ifndef NDEBUG
271  int nvectors = X.NumVectors();
272  assert (nvectors == data->localrhs->NumVectors());
273 #endif // NDEBUG
274 
275  // Initialize the X vector for iterative solver
276  data->Xs->PutScalar(0.0);
277 
278  // Get local portion of X
279  data->localrhs->Import(X, *(data->BdImporter), Insert);
280 
281  // locallhs is z in paper
282  if (config->amesosForDiagonal) {
283  ssym->OrigLP->SetRHS((data->localrhs).getRawPtr());
284  ssym->OrigLP->SetLHS((data->locallhs).getRawPtr());
285  ssym->ReIdx_LP->fwd();
286  ssym->Solver->Solve();
287  }
288  else {
289  ssym->ifSolver->ApplyInverse(*(data->localrhs), *(data->locallhs));
290  }
291 
292  err = ssym->R->Multiply(false, *(data->locallhs), *(data->temp1));
293  assert (err == 0);
294 
295  // Export temp1 to a dist vector - temp2
296  data->temp2->Import(*(data->temp1), *(data->DistImporter), Insert);
297 
298  //Epetra_MultiVector Bs(SMap, nvectors); // b_2 - R * z in ShyLU paper
299  data->Bs->Import(X, *(data->BsImporter), Insert);
300  data->Bs->Update(-1.0, *(data->temp2), 1.0);
301 
302  AztecOO *solver = 0;
303  Epetra_LinearProblem Problem(data->Sbar.get(),
304  (data->Xs).getRawPtr(), (data->Bs).getRawPtr());
305  if ((config->schurSolver == "G") || (config->schurSolver == "IQR"))
306  {
307  IFPACK_CHK_ERR(data->iqrSolver->Solve(*(data->schur_op),
308  *(data->Bs), *(data->Xs)));
309  }
310  else if (config->schurSolver == "Amesos")
311  {
312  Amesos_BaseSolver *solver2 = data->dsolver;
313  data->OrigLP2->SetLHS((data->Xs).getRawPtr());
314  data->OrigLP2->SetRHS((data->Bs).getRawPtr());
315  data->ReIdx_LP2->fwd();
316  //cout << "Calling solve *****************************" << endl;
317  solver2->Solve();
318  //cout << "Out of solve *****************************" << endl;
319  }
320  else
321  {
322  if (config->libName == "Belos")
323  {
324  solver = data->innersolver;
325  solver->SetLHS((data->Xs).getRawPtr());
326  solver->SetRHS((data->Bs).getRawPtr());
327  }
328  else
329  {
330  // See the comment above on why we are not able to reuse the solver
331  // when outer solve is AztecOO as well.
332  solver = new AztecOO();
333  //solver.SetPrecOperator(precop_);
334  solver->SetAztecOption(AZ_solver, AZ_gmres);
335  // Do not use AZ_none
336  solver->SetAztecOption(AZ_precond, AZ_dom_decomp);
337  //solver->SetAztecOption(AZ_precond, AZ_none);
338  //solver->SetAztecOption(AZ_precond, AZ_Jacobi);
340  //solver->SetAztecOption(AZ_overlap, 3);
341  //solver->SetAztecOption(AZ_subdomain_solve, AZ_ilu);
342  //solver->SetAztecOption(AZ_output, AZ_all);
343  //solver->SetAztecOption(AZ_diagnostics, AZ_all);
344  solver->SetProblem(Problem);
345  }
346 
347  // What should be a good inner_tolerance :-) ?
348  solver->Iterate(config->inner_maxiters, config->inner_tolerance);
349  }
350 
351  // Import Xs locally
352  data->LocalXs->Import(*(data->Xs), *(data->XsImporter), Insert);
353 
354  err = ssym->C->Multiply(false, *(data->LocalXs), *(data->temp3));
355  assert (err == 0);
356  data->temp3->Update(1.0, *(data->localrhs), -1.0);
357 
358  if (config->amesosForDiagonal) {
359  ssym->OrigLP->SetRHS((data->temp3).getRawPtr());
360  ssym->OrigLP->SetLHS((data->locallhs).getRawPtr());
361  ssym->ReIdx_LP->fwd();
362  ssym->Solver->Solve();
363  }
364  else {
365  ssym->ifSolver->ApplyInverse(*(data->temp3), *(data->locallhs));
366  }
367 
368  Y.Export(*(data->locallhs), *(data->XdExporter), Insert);
369  Y.Export(*(data->LocalXs), *(data->XsExporter), Insert);
370 
371  if (config->libName == "Belos" || config->schurSolver == "Amesos")
372  {
373  // clean up
374  }
375  else
376  {
377  delete solver;
378  }
379  return 0;
380 }
381 
383  shylu_symbolic *ssym,
384  shylu_data *data,
385  shylu_config *config,
386  const Epetra_MultiVector& X,
387  Epetra_MultiVector& Y
388 )
389 {
390  if (config->sep_type != 1)
391  shylu_dist_solve(ssym, data, config, X, Y);
392  else
393  shylu_local_solve(ssym, data, config, X, Y);
394 
395  return 0;
396 }
int shylu_solve(shylu_symbolic *ssym, shylu_data *data, shylu_config *config, const Epetra_MultiVector &X, Epetra_MultiVector &Y)
Call solve on multiple RHS.
T * get() const
Utilities for ShyLU.
Main data structure holding needed offset and temp variables.
Definition: shylu.h:108
Main header file of ShyLU (Include main user calls)