27 #include "Komplex_LinearProblem.h"
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)
41 template<
typename Scalar>
44 this->free_permutation_data();
47 template<
typename Scalar>
53 template<
typename Scalar>
70 template<
typename Scalar>
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;
88 aztec.SetAztecOption(AZ_solver, az_solver);
91 template<
typename Scalar>
94 this->tolerance = tol;
97 template<
typename Scalar>
100 this->max_iters = iters;
103 template<
typename Scalar>
108 this->precond_yes =
true;
112 template<
typename Scalar>
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)
127 az_precond = AZ_none;
131 az_precond = AZ_none;
133 this->precond_yes = (az_precond != AZ_none);
134 aztec.SetAztecOption(AZ_precond, az_precond);
137 template<
typename Scalar>
140 aztec.SetAztecOption(option, value);
143 template<
typename Scalar>
146 aztec.SetAztecParam(param, value);
149 template<
typename Scalar>
163 template<
typename Scalar>
169 template<
typename Scalar>
173 this->free_permutation_data();
176 template<
typename Scalar>
180 this->free_permutation_data();
183 template<
typename Scalar>
187 int ndof_per_pde = ndof / this->n_eq;
191 this->row_perm = malloc_with_check<AztecOOSolver<Scalar>,
int>(ndof,
this);
192 this->col_perm = malloc_with_check<AztecOOSolver<Scalar>,
int>(ndof,
this);
194 for (
int i = 0; i < ndof; i++)
196 this->row_perm[i] = j%ndof + jc;
197 this->col_perm[i] = k%ndof + kc;
218 assert(m !=
nullptr);
219 assert(rhs !=
nullptr);
220 assert(m->size == rhs->size);
224 aztec.SetAztecOption(AZ_output, AZ_summary);
230 if (node_wise_ordering)
233 free_permutation_data();
235 create_permutation_vectors();
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);
258 aztec.SetUserMatrix(final_matrix->mat);
261 Epetra_Vector *final_rhs;
264 EpetraExt::Permutation<Epetra_MultiVector> Pv(Copy, rhs->vec->Map(), row_perm);
265 final_rhs = Pv(*rhs->vec)(0);
268 final_rhs = rhs->vec;
270 aztec.SetRHS(final_rhs);
272 Epetra_Vector x(final_matrix->mat->DomainMap());
280 pc->create(final_matrix);
282 aztec.SetPrecOperator(pc->get_obj());
291 assert(pc->get_obj());
293 ML_Epetra::MultiLevelPreconditioner* op_ml =
dynamic_cast<ML_Epetra::MultiLevelPreconditioner*
>(pc->get_obj());
294 assert(op_ml->IsPreconditionerComputed());
300 aztec.Iterate(this->max_iters, this->tolerance);
303 this->time = this->accumulated();
305 free_with_check(this->sln);
306 this->sln = malloc_with_check<AztecOOSolver<double>,
double>(final_matrix->size,
this);
308 memset(this->sln, 0, final_matrix->size *
sizeof(
double));
312 for (
unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[col_perm[i]];
314 for (
unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[i];
320 assert(m !=
nullptr);
321 assert(rhs !=
nullptr);
322 assert(m->size == rhs->size);
326 aztec.SetAztecOption(AZ_output, AZ_summary);
328 if (this->get_verbose_output())
329 aztec.SetAztecOption(AZ_output, AZ_all);
335 if (node_wise_ordering)
338 free_permutation_data();
340 create_permutation_vectors();
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);
363 aztec.SetUserMatrix(final_matrix->mat);
366 Epetra_Vector *final_rhs;
369 EpetraExt::Permutation<Epetra_MultiVector> Pv(Copy, rhs->vec->Map(), row_perm);
370 final_rhs = Pv(*rhs->vec)(0);
373 final_rhs = rhs->vec;
375 aztec.SetRHS(final_rhs);
377 Epetra_Vector x(final_matrix->mat->DomainMap());
382 for (
unsigned int i = 0; i < m->size; i++)
383 x[i] = initial_guess[row_perm[i]];
385 for (
unsigned int i = 0; i < m->size; i++)
386 x[i] = initial_guess[i];
396 pc->create(final_matrix);
398 aztec.SetPrecOperator(pc->get_obj());
407 assert(pc->get_obj());
409 ML_Epetra::MultiLevelPreconditioner* op_ml =
dynamic_cast<ML_Epetra::MultiLevelPreconditioner*
>(pc->get_obj());
410 assert(op_ml->IsPreconditionerComputed());
416 aztec.Iterate(this->max_iters, this->tolerance);
419 this->time = this->accumulated();
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));
427 for (
unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[col_perm[i]];
429 for (
unsigned int i = 0; i < final_matrix->size; i++) this->sln[i] = x[i];
435 assert(m !=
nullptr);
436 assert(rhs !=
nullptr);
437 assert(m->size == rhs->size);
441 aztec.SetAztecOption(AZ_output, AZ_none);
443 double c0r = 1.0, c0i = 0.0;
444 double c1r = 0.0, c1i = 1.0;
446 Epetra_Vector xr(*rhs->std_map);
447 Epetra_Vector xi(*rhs->std_map);
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);
454 aztec.Iterate(this->max_iters, this->tolerance);
456 kp.ExtractSolution(xr, xi);
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>));
463 for (
unsigned int i = 0; i < m->size; i++)
464 this->sln[i] = std::complex<double>(xr[i], xi[i]);
472 assert(m !=
nullptr);
473 assert(rhs !=
nullptr);
474 assert(m->size == rhs->size);
478 aztec.SetAztecOption(AZ_output, AZ_none);
480 double c0r = 1.0, c0i = 0.0;
481 double c1r = 0.0, c1i = 1.0;
483 Epetra_Vector xr(*rhs->std_map);
484 Epetra_Vector xi(*rhs->std_map);
488 for (
unsigned int i = 0; i < m->size; i++)
490 xr[i] = initial_guess[i].real();
491 xi[i] = initial_guess[i].imag();
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);
500 aztec.Iterate(this->max_iters, this->tolerance);
502 kp.ExtractSolution(xr, xi);
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>));
509 for (
unsigned int i = 0; i < m->size; i++)
510 this->sln[i] = std::complex<double>(xr[i], xi[i]);
514 template<
typename Scalar>
517 return aztec.NumIters();
520 template<
typename Scalar>
523 return aztec.TrueResidual();
General namespace for the Hermes library.
Encapsulation of AztecOO linear solver.
Abstract class to define interface for preconditioners.
Abstract class for Epetra preconditioners.
AztecOOSolver class as an interface to AztecOO.
Abstract middle-class for solvers that work in a loop of a kind (iterative, multigrid, ...)
existing factorization data.
File containing functionality for investigating call stack.
Abstract class for defining solver interface.
MatrixStructureReuseScheme
Abstract class for defining interface for iterative solvers. Internal, though utilizable for defining...
pattern as the one already factorized.