Commit 67336069 authored by Jack Poulson's avatar Jack Poulson

Adding a control structure for refined solves.

parent a26e92a1
*.swp
build*/
docs/_build/*
subprojects/*
!subprojects/*.wrap
......@@ -78,11 +78,12 @@ right_hand_sides.Resize(num_rows, num_rhs);
factorization.Solve(&right_hand_sides.view);
// Alternatively, one can solve using iterative-refinement, e.g., using:
const Int max_refine_iters = 3;
const double relative_tol = 1e-15;
const bool verbose = true;
catamari::RefinedSolveControl<double> refined_solve_control;
refined_solve_control.relative_tol = 1e-15;
refined_solve_control.max_iters = 3;
refined_solve_control.verbose = true;
factorization.RefinedSolve(
matrix, relative_tol, max_refine_iters, verbose, &right_hand_sides.view);
matrix, refined_solve_control, &right_hand_sides.view);
```
One can also browse the [example/](https://gitlab.com/hodge_star/catamari/tree/master/example) folder for complete examples (e.g., for [solving 3D Helmholtz equations](https://gitlab.com/hodge_star/catamari/blob/master/example/helmholtz_3d_pml.cc) with PML boundary conditions discretized using trilinear hexahedral elements using a complex LDL^T factorization).
......
......@@ -60,11 +60,12 @@ Usage of the sparse-direct solver through the
factorization.Solve(&right_hand_sides.view);
// Alternatively, one can solve using iterative-refinement, e.g., using:
const Int max_refine_iters = 3;
const double relative_tol = 1e-15;
const bool verbose = true;
catamari::RefinedSolveControl<double> refined_solve_control;
refined_solve_control.relative_tol = 1e-15;
refined_solve_control.max_iters = 3;
refined_solve_control.verbose = true;
factorization.RefinedSolve(
matrix, relative_tol, max_refine_iters, verbose, &right_hand_sides.view);
matrix, refined_solve_control, &right_hand_sides.view);
One can also browse the `example/ <https://gitlab.com/hodge_star/catamari/tree/master/example>`_ folder for complete examples (e.g., for `solving 3D Helmholtz equations <https://gitlab.com/hodge_star/catamari/blob/master/example/helmholtz_3d_pml.cc>`_ with PML boundary conditions discretized using trilinear hexahedral elements using a complex LDL^T factorization).
......
......@@ -821,17 +821,16 @@ Experiment RunTest(SpeedProfile profile, const double& omega,
// Solve the linear system with iterative refinement.
{
// TODO(Jack Poulson): Make these parameters configurable.
const Real relative_tol = 1e-15;
const Int max_refine_iters = 3;
const bool verbose = true;
catamari::RefinedSolveControl<Real> refined_solve_control;
refined_solve_control.verbose = true;
if (print_progress) {
std::cout << " Running iteratively-refined solve..." << std::endl;
}
BlasMatrix<Field> solution = right_hand_sides;
timer.Start();
ldl_factorization.RefinedSolve(matrix, relative_tol, max_refine_iters,
verbose, &solution.view);
ldl_factorization.RefinedSolve(
matrix, refined_solve_control, &solution.view);
experiment.refined_solve_seconds = timer.Stop();
if (print_progress) {
......
......@@ -988,17 +988,16 @@ Experiment RunTest(SpeedProfile profile, const double& omega,
// Solve the linear systems using iterative refinement.
{
// TODO(Jack Poulson): Make these parameters configurable.
const Real relative_tol = 1e-15;
const Int max_refine_iters = 3;
const bool verbose = true;
catamari::RefinedSolveControl<Real> refined_solve_control;
refined_solve_control.verbose = true;
if (print_progress) {
std::cout << " Running iteratively-refined solve..." << std::endl;
}
BlasMatrix<Field> solution = right_hand_sides;
timer.Start();
ldl_factorization.RefinedSolve(matrix, relative_tol, max_refine_iters,
verbose, &solution.view);
ldl_factorization.RefinedSolve(
matrix, refined_solve_control, &solution.view);
experiment.refined_solve_seconds = timer.Stop();
if (print_progress) {
......
......@@ -116,13 +116,12 @@ void LDLFactorization<Field>::Solve(
template <class Field>
Int LDLFactorization<Field>::RefinedSolve(
const CoordinateMatrix<Field>& matrix, ComplexBase<Field> relative_tol,
Int max_refine_iters, bool verbose,
const CoordinateMatrix<Field>& matrix,
const RefinedSolveControl<Real>& control,
BlasMatrixView<Field>* right_hand_sides) const {
typedef ComplexBase<Field> Real;
const Int num_rows = matrix.NumRows();
const Int num_rhs = right_hand_sides->width;
if (max_refine_iters <= 0) {
if (control.max_iters <= 0) {
Solve(right_hand_sides);
return 0;
}
......@@ -155,7 +154,7 @@ Int LDLFactorization<Field>::RefinedSolve(
}
error_norms[j] = MaxNorm(column.ToConst());
if (verbose) {
if (control.verbose) {
std::cout << "Original relative error " << j << ": "
<< error_norms[j] / rhs_orig_norms[j] << std::endl;
}
......@@ -180,11 +179,11 @@ Int LDLFactorization<Field>::RefinedSolve(
const Real error_norm = error_norms[j];
const Real rhs_orig_norm = rhs_orig_norms[j];
const Real relative_error = error_norm / rhs_orig_norm;
if (relative_error <= relative_tol) {
if (verbose) {
if (relative_error <= control.relative_tol) {
if (control.verbose) {
std::cout << "Relative error " << j << " (" << j_active
<< "): " << relative_error << " <= " << relative_tol
<< std::endl;
<< "): " << relative_error << " <= "
<< control.relative_tol << std::endl;
}
} else {
active_indices[num_remaining++] = j;
......@@ -233,7 +232,7 @@ Int LDLFactorization<Field>::RefinedSolve(
column(i, 0) = rhs_orig(i, j) - image(i, j_active);
}
const Real new_error_norm = MaxNorm(column.ToConst());
if (verbose) {
if (control.verbose) {
std::cout << "Refined relative error " << j << ": "
<< new_error_norm / rhs_orig_norms[j] << std::endl;
}
......@@ -243,7 +242,7 @@ Int LDLFactorization<Field>::RefinedSolve(
}
error_norms[j] = new_error_norm;
active_indices[num_remaining++] = j;
} else if (verbose) {
} else if (control.verbose) {
std::cout << "Right-hand side " << j << "(" << j_active << ") diverged."
<< std::endl;
}
......@@ -251,7 +250,7 @@ Int LDLFactorization<Field>::RefinedSolve(
active_indices.Resize(num_remaining);
++refine_iter;
if (refine_iter >= max_refine_iters || active_indices.Empty()) {
if (refine_iter >= control.max_iters || active_indices.Empty()) {
break;
}
}
......
......@@ -8,6 +8,7 @@
#ifndef CATAMARI_LDL_H_
#define CATAMARI_LDL_H_
#include <limits>
#include <memory>
#include "catamari/ldl/scalar_ldl.hpp"
......@@ -48,10 +49,27 @@ struct LDLControl {
}
};
// Configuration options for solving linear systems with iterative refinement.
template <typename Real>
struct RefinedSolveControl {
// The desired relative error (in the max norm) of the solution. Iteration
// stops early if this tolerance is achieved.
Real relative_tol = 10 * std::numeric_limits<Real>::epsilon();
// The maximum number of iterations of iterative refinement to perform.
Int max_iters = 3;
// Whether
bool verbose = false;
};
// A wrapper for the scalar and supernodal factorization data structures.
template <class Field>
class LDLFactorization {
public:
// The underlying real datatype of the scalar type.
typedef ComplexBase<Field> Real;
// Whether or not a supernodal factorization was used. If it is true, only
// 'supernodal_factorization' should be non-null, and vice versa.
bool is_supernodal;
......@@ -77,8 +95,8 @@ class LDLFactorization {
// Solves a set of linear systems using iterative refinement.
Int RefinedSolve(const CoordinateMatrix<Field>& matrix,
ComplexBase<Field> relative_tol, Int max_refine_iters,
bool verbose, BlasMatrixView<Field>* right_hand_sides) const;
const RefinedSolveControl<Real>& control,
BlasMatrixView<Field>* right_hand_sides) const;
// Solves a set of linear systems using the lower-triangular factor.
void LowerTriangularSolve(BlasMatrixView<Field>* right_hand_sides) const;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment