HermesCommon  3.0
aztecoo_solver.cpp
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
4 // Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #include "config.h"
23 #ifdef HAVE_AZTECOO
24 #include "aztecoo_solver.h"
25 #include "callstack.h"
26 #ifdef HAVE_KOMPLEX
27 #include "Komplex_LinearProblem.h"
28 #endif
29 
30 namespace Hermes
31 {
32  namespace Solvers
33  {
34  template<typename Scalar>
36  : IterSolver<Scalar>(m, rhs), LoopSolver<Scalar>(m, rhs), m(m), rhs(rhs), final_matrix(nullptr), P(nullptr), Q(nullptr), row_perm(nullptr), col_perm(nullptr)
37  {
38  pc = nullptr;
39  }
40 
41  template<typename Scalar>
43  {
44  this->free_permutation_data();
45  }
46 
47  template<typename Scalar>
49  {
50  this->free();
51  }
52 
53  template<typename Scalar>
55  {
56  if (row_perm)
57  {
58  assert(col_perm);
59  assert(Q);
60  assert(final_matrix);
61 
62  delete row_perm;
63  delete col_perm;
64  delete P;
65  delete Q;
66  delete final_matrix;
67  }
68  }
69 
70  template<typename Scalar>
71  void AztecOOSolver<Scalar>::set_solver(const char *name)
72  {
73  int az_solver;
74  if (name)
75  {
76  if (strcasecmp(name, "gmres") == 0) az_solver = AZ_gmres;
77  else if (strcasecmp(name, "cg") == 0) az_solver = AZ_cg;
78  else if (strcasecmp(name, "cgs") == 0) az_solver = AZ_cgs;
79  else if (strcasecmp(name, "tfqmr") == 0) az_solver = AZ_tfqmr;
80  else if (strcasecmp(name, "bicgstab") == 0) az_solver = AZ_bicgstab;
81  else az_solver = AZ_gmres;
82  }
83  else
84  {
85  az_solver = AZ_gmres;
86  }
87 
88  aztec.SetAztecOption(AZ_solver, az_solver);
89  }
90 
91  template<typename Scalar>
93  {
94  this->tolerance = tol;
95  }
96 
97  template<typename Scalar>
99  {
100  this->max_iters = iters;
101  }
102 
103  template<typename Scalar>
105  {
106  this->pc = dynamic_cast<EpetraPrecond<Scalar>*>(pc);
107  if (this->pc)
108  this->precond_yes = true;
109  // TODO: else warn that a wrong type of preconditioner has been used.
110  }
111 
112  template<typename Scalar>
113  void AztecOOSolver<Scalar>::set_precond(const char *name)
114  {
115  int az_precond;
116  if (name)
117  {
118  if (strcasecmp(name, "none") == 0)
119  az_precond = AZ_none;
120  else if (strcasecmp(name, "jacobi") == 0)
121  az_precond = AZ_Jacobi;
122  else if (strcasecmp(name, "neumann") == 0)
123  az_precond = AZ_Neumann;
124  else if (strcasecmp(name, "least-squares") == 0)
125  az_precond = AZ_ls;
126  else
127  az_precond = AZ_none;
128  }
129  else
130  //asi by to melo byt nastaveno
131  az_precond = AZ_none;
132 
133  this->precond_yes = (az_precond != AZ_none);
134  aztec.SetAztecOption(AZ_precond, az_precond);
135  }
136 
137  template<typename Scalar>
138  void AztecOOSolver<Scalar>::set_option(int option, int value)
139  {
140  aztec.SetAztecOption(option, value);
141  }
142 
143  template<typename Scalar>
144  void AztecOOSolver<Scalar>::set_param(int param, double value)
145  {
146  aztec.SetAztecParam(param, value);
147  }
148 
149  template<typename Scalar>
151  {
153  /*
154  if (reuse_scheme == HERMES_REUSE_MATRIX_STRUCTURE_COMPLETELY || reuse_scheme == HERMES_REUSE_MATRIX_REORDERING_AND_SCALING)
155  aztec.SetAztecOption(AZ_pre_calc, AZ_reuse);
156  else if (reuse_scheme == HERMES_REUSE_MATRIX_REORDERING)
157  aztec.SetAztecOption(AZ_pre_calc, AZ_recalc);
158  else
159  aztec.SetAztecOption(AZ_pre_calc, AZ_calc);
160  */
161  }
162 
163  template<typename Scalar>
165  {
166  return m->size;
167  }
168 
169  template<typename Scalar>
171  {
173  this->free_permutation_data();
174  }
175 
176  template<typename Scalar>
178  {
180  this->free_permutation_data();
181  }
182 
183  template<typename Scalar>
185  {
186  int ndof = m->size;
187  int ndof_per_pde = ndof / this->n_eq;
188  int j = 0, jc = 0;
189  int k = 0, kc = 0;
190 
191  this->row_perm = malloc_with_check<AztecOOSolver<Scalar>, int>(ndof, this);
192  this->col_perm = malloc_with_check<AztecOOSolver<Scalar>, int>(ndof, this);
193 
194  for (int i = 0; i < ndof; i++)
195  {
196  this->row_perm[i] = j%ndof + jc;
197  this->col_perm[i] = k%ndof + kc;
198 
199  j += ndof_per_pde;
200  if (j == ndof)
201  {
202  jc++;
203  j = 0;
204  }
205 
206  k += this->n_eq;
207  if (k == ndof)
208  {
209  kc++;
210  k = 0;
211  }
212  }
213  }
214 
215  template<>
217  {
218  assert(m != nullptr);
219  assert(rhs != nullptr);
220  assert(m->size == rhs->size);
221 
222  // no output
223  // AZ_all | AZ_warnings | AZ_last | AZ_summary
224  aztec.SetAztecOption(AZ_output, AZ_summary);
225 
226  // setup the problem
227  if (reuse_scheme == HERMES_CREATE_STRUCTURE_FROM_SCRATCH)
228  //if (aztec.GetAztecOption(AZ_pre_calc) == AZ_calc)
229  {
230  if (node_wise_ordering)
231  {
232  if (row_perm)
233  free_permutation_data();
234 
235  create_permutation_vectors();
236 
237  // NOTE: RowMap() == RangeMap() == ColMap() == DomainMap()
238  P = new EpetraExt::Permutation<Epetra_CrsMatrix>(Copy, m->mat->RowMap(), row_perm);
239  Q = new EpetraExt::Permutation<Epetra_CrsMatrix>(Copy, m->mat->ColMap(), col_perm);
240 
241  final_matrix = new EpetraMatrix<double>((*P)((*Q)(*m->mat, true)));
242 
243  // NOTE: According to Trilinos docs, final_matrix should be fill_completed by now. However, when row permutation is performed,
244  // the input matrix is not finalized (unlike the case of column permutation) - this may possibly be a bug in Trilinos.
245  // Also, when doing a row permutation, all maps of the source matrix are set to a temporary permutation
246  // map in order to export the values from the source matrix to the target matrix in a permuted way. Then,
247  // according to code documentation in EpetraExt_Permutation_impl.h, the row indexing for the new_ permuted
248  // matrix (specified by the target matrix' RowMap_) is set to the original indexing (the original matrix'
249  // RowMap_). However, the other maps of the new_ matrix are not reset (and still point to the temporary
250  // permutation map), which causes another set of problems. Hence, replaceMap() implementation in EpetraExt_Permutation_impl.h
251  // has to be changed to call FillComplete with correct maps as input (calling FillComplete here doesn't suffice, since it
252  // does not prevent creating a non-trivial Exporter, which is the root cause of the problems (and which cannot be deleted
253  // from outside.)
254  }
255  else
256  final_matrix = m;
257 
258  aztec.SetUserMatrix(final_matrix->mat);
259  }
260 
261  Epetra_Vector *final_rhs;
262  if (row_perm)
263  {
264  EpetraExt::Permutation<Epetra_MultiVector> Pv(Copy, rhs->vec->Map(), row_perm);
265  final_rhs = Pv(*rhs->vec)(0);
266  }
267  else
268  final_rhs = rhs->vec;
269 
270  aztec.SetRHS(final_rhs);
271 
272  Epetra_Vector x(final_matrix->mat->DomainMap());
273  aztec.SetLHS(&x);
274 
275  if (pc != nullptr)
276  {
277  if (reuse_scheme == HERMES_CREATE_STRUCTURE_FROM_SCRATCH)
278  //if(aztec.GetAztecOption(AZ_pre_calc) == AZ_calc)
279  {
280  pc->create(final_matrix);
281  pc->compute();
282  aztec.SetPrecOperator(pc->get_obj());
283  }
284  else if (reuse_scheme == HERMES_REUSE_MATRIX_REORDERING || reuse_scheme == HERMES_REUSE_MATRIX_REORDERING_AND_SCALING)
285  //else if (aztec.GetAztecOption(AZ_pre_calc) == AZ_recalc)
286  {
287  pc->recompute();
288  }
289  else
290  {
291  assert(pc->get_obj());
292 #ifdef HAVE_ML
293  ML_Epetra::MultiLevelPreconditioner* op_ml = dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(pc->get_obj());
294  assert(op_ml->IsPreconditionerComputed());
295 #endif
296  }
297  }
298 
299  // solve it
300  aztec.Iterate(this->max_iters, this->tolerance);
301 
302  this->tick();
303  this->time = this->accumulated();
304 
305  free_with_check(this->sln);
306  this->sln = malloc_with_check<AztecOOSolver<double>, double>(final_matrix->size, this);
307 
308  memset(this->sln, 0, final_matrix->size * sizeof(double));
309 
310  // copy the solution into sln vector
311  if (col_perm)
312  for (unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[col_perm[i]];
313  else
314  for (unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[i];
315  }
316 
317  template<>
318  void AztecOOSolver<double>::solve(double *initial_guess)
319  {
320  assert(m != nullptr);
321  assert(rhs != nullptr);
322  assert(m->size == rhs->size);
323 
324  // no output
325  // AZ_all | AZ_warnings | AZ_last | AZ_summary
326  aztec.SetAztecOption(AZ_output, AZ_summary);
327 
328  if (this->get_verbose_output())
329  aztec.SetAztecOption(AZ_output, AZ_all);
330 
331  // setup the problem
332  if (reuse_scheme == HERMES_CREATE_STRUCTURE_FROM_SCRATCH)
333  //if (aztec.GetAztecOption(AZ_pre_calc) == AZ_calc)
334  {
335  if (node_wise_ordering)
336  {
337  if (row_perm)
338  free_permutation_data();
339 
340  create_permutation_vectors();
341 
342  // NOTE: RowMap() == RangeMap() == ColMap() == DomainMap()
343  P = new EpetraExt::Permutation<Epetra_CrsMatrix>(Copy, m->mat->RowMap(), row_perm);
344  Q = new EpetraExt::Permutation<Epetra_CrsMatrix>(Copy, m->mat->ColMap(), col_perm);
345 
346  final_matrix = new EpetraMatrix<double>((*P)((*Q)(*m->mat, true)));
347 
348  // NOTE: According to Trilinos docs, final_matrix should be fill_completed by now. However, when row permutation is performed,
349  // the input matrix is not finalized (unlike the case of column permutation) - this may possibly be a bug in Trilinos.
350  // Also, when doing a row permutation, all maps of the source matrix are set to a temporary permutation
351  // map in order to export the values from the source matrix to the target matrix in a permuted way. Then,
352  // according to code documentation in EpetraExt_Permutation_impl.h, the row indexing for the new_ permuted
353  // matrix (specified by the target matrix' RowMap_) is set to the original indexing (the original matrix'
354  // RowMap_). However, the other maps of the new_ matrix are not reset (and still point to the temporary
355  // permutation map), which causes another set of problems. Hence, replaceMap() implementation in EpetraExt_Permutation_impl.h
356  // has to be changed to call FillComplete with correct maps as input (calling FillComplete here doesn't suffice, since it
357  // does not prevent creating a non-trivial Exporter, which is the root cause of the problems (and which cannot be deleted
358  // from outside.)
359  }
360  else
361  final_matrix = m;
362 
363  aztec.SetUserMatrix(final_matrix->mat);
364  }
365 
366  Epetra_Vector *final_rhs;
367  if (row_perm)
368  {
369  EpetraExt::Permutation<Epetra_MultiVector> Pv(Copy, rhs->vec->Map(), row_perm);
370  final_rhs = Pv(*rhs->vec)(0);
371  }
372  else
373  final_rhs = rhs->vec;
374 
375  aztec.SetRHS(final_rhs);
376 
377  Epetra_Vector x(final_matrix->mat->DomainMap());
378 
379  if (initial_guess)
380  {
381  if (row_perm)
382  for (unsigned int i = 0; i < m->size; i++)
383  x[i] = initial_guess[row_perm[i]];
384  else
385  for (unsigned int i = 0; i < m->size; i++)
386  x[i] = initial_guess[i];
387  }
388 
389  aztec.SetLHS(&x);
390 
391  if (pc != nullptr)
392  {
393  if (reuse_scheme == HERMES_CREATE_STRUCTURE_FROM_SCRATCH)
394  //if(aztec.GetAztecOption(AZ_pre_calc) == AZ_calc)
395  {
396  pc->create(final_matrix);
397  pc->compute();
398  aztec.SetPrecOperator(pc->get_obj());
399  }
400  else if (reuse_scheme == HERMES_REUSE_MATRIX_REORDERING || reuse_scheme == HERMES_REUSE_MATRIX_REORDERING_AND_SCALING)
401  //else if (aztec.GetAztecOption(AZ_pre_calc) == AZ_recalc)
402  {
403  pc->recompute();
404  }
405  else
406  {
407  assert(pc->get_obj());
408 #ifdef HAVE_ML
409  ML_Epetra::MultiLevelPreconditioner* op_ml = dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(pc->get_obj());
410  assert(op_ml->IsPreconditionerComputed());
411 #endif
412  }
413  }
414 
415  // solve it
416  aztec.Iterate(this->max_iters, this->tolerance);
417 
418  this->tick();
419  this->time = this->accumulated();
420 
421  free_with_check(this->sln);
422  this->sln = malloc_with_check<AztecOOSolver<double>, double>(final_matrix->size, this);
423  memset(this->sln, 0, final_matrix->size * sizeof(double));
424 
425  // copy the solution into sln vector
426  if (col_perm)
427  for (unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[col_perm[i]];
428  else
429  for (unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[i];
430  }
431  template<>
433  {
434 #ifdef HAVE_KOMPLEX
435  assert(m != nullptr);
436  assert(rhs != nullptr);
437  assert(m->size == rhs->size);
438 
439  // no output
440  // AZ_all | AZ_warnings | AZ_last | AZ_summary
441  aztec.SetAztecOption(AZ_output, AZ_none);
442 
443  double c0r = 1.0, c0i = 0.0;
444  double c1r = 0.0, c1i = 1.0;
445 
446  Epetra_Vector xr(*rhs->std_map);
447  Epetra_Vector xi(*rhs->std_map);
448 
449  Komplex_LinearProblem kp(c0r, c0i, *m->mat, c1r, c1i, *m->mat_im, xr, xi, *rhs->vec, *rhs->vec_im);
450  Epetra_LinearProblem *lp = kp.KomplexProblem();
451  aztec.SetProblem(*lp, true);
452 
453  // solve it
454  aztec.Iterate(this->max_iters, this->tolerance);
455 
456  kp.ExtractSolution(xr, xi);
457 
458  free_with_check(this->sln);
459  this->sln = malloc_with_check<AztecOOSolver<std::complex<double> >, std::complex<double> >(m->size, this);
460  memset(this->sln, 0, m->size * sizeof(std::complex<double>));
461 
462  // copy the solution into sln vector
463  for (unsigned int i = 0; i < m->size; i++)
464  this->sln[i] = std::complex<double>(xr[i], xi[i]);
465 #endif
466  }
467 
468  template<>
469  void AztecOOSolver<std::complex<double> >::solve(std::complex<double>* initial_guess)
470  {
471 #ifdef HAVE_KOMPLEX
472  assert(m != nullptr);
473  assert(rhs != nullptr);
474  assert(m->size == rhs->size);
475 
476  // no output
477  // AZ_all | AZ_warnings | AZ_last | AZ_summary
478  aztec.SetAztecOption(AZ_output, AZ_none);
479 
480  double c0r = 1.0, c0i = 0.0;
481  double c1r = 0.0, c1i = 1.0;
482 
483  Epetra_Vector xr(*rhs->std_map);
484  Epetra_Vector xi(*rhs->std_map);
485 
486  if (initial_guess)
487  {
488  for (unsigned int i = 0; i < m->size; i++)
489  {
490  xr[i] = initial_guess[i].real();
491  xi[i] = initial_guess[i].imag();
492  }
493  }
494 
495  Komplex_LinearProblem kp(c0r, c0i, *m->mat, c1r, c1i, *m->mat_im, xr, xi, *rhs->vec, *rhs->vec_im);
496  Epetra_LinearProblem *lp = kp.KomplexProblem();
497  aztec.SetProblem(*lp, true);
498 
499  // solve it
500  aztec.Iterate(this->max_iters, this->tolerance);
501 
502  kp.ExtractSolution(xr, xi);
503 
504  free_with_check(this->sln);
505  this->sln = malloc_with_check<AztecOOSolver<std::complex<double> >, std::complex<double> >(m->size, this);
506  memset(this->sln, 0, m->size * sizeof(std::complex<double>));
507 
508  // copy the solution into sln vector
509  for (unsigned int i = 0; i < m->size; i++)
510  this->sln[i] = std::complex<double>(xr[i], xi[i]);
511 #endif
512  }
513 
514  template<typename Scalar>
516  {
517  return aztec.NumIters();
518  }
519 
520  template<typename Scalar>
522  {
523  return aztec.TrueResidual();
524  }
525 
526  template class HERMES_API AztecOOSolver < double > ;
527  template class HERMES_API AztecOOSolver < std::complex<double> > ;
528  }
529 }
530 #endif
General namespace for the Hermes library.
Encapsulation of AztecOO linear solver.
Abstract class to define interface for preconditioners.
Definition: precond.h:57
Abstract class for Epetra preconditioners.
Definition: precond.h:66
AztecOOSolver class as an interface to AztecOO.
Abstract middle-class for solvers that work in a loop of a kind (iterative, multigrid, ...)
File containing functionality for investigating call stack.
Abstract class for defining solver interface.
Abstract class for defining interface for iterative solvers. Internal, though utilizable for defining...
pattern as the one already factorized.