Skip to content

Commit

Permalink
Update PGMRES solver
Browse files Browse the repository at this point in the history
  • Loading branch information
@jchristopherson
jchristopherson committed Feb 13, 2024
1 parent 3766ef1 commit 3d03dfc
Show file tree
Hide file tree
Showing 4 changed files with 148 additions and 12 deletions.
63 changes: 63 additions & 0 deletions src/linalg.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,7 @@ module linalg
public :: assignment(=)
public :: transpose
public :: sparse_direct_solve
public :: pgmres_solver
public :: LA_NO_OPERATION
public :: LA_TRANSPOSE
public :: LA_HERMITIAN_TRANSPOSE
Expand Down Expand Up @@ -5482,6 +5483,56 @@ module subroutine sort_eigen_dbl(vals, vecs, ascend, err)
module procedure :: csr_solve_sparse_direct
end interface

!> @brief A preconditioned GMRES solver.
!!
!! @par Syntax
!! @code{.f90}
!! subroutine pgmres_solver( &
!! class(csr_matrix) a, &
!! class(msr_matrix) lu, &
!! integer(int32) ju(:), &
!! real(real64) b(:), &
!! real(real64) x(:), &
!! optional integer(int32) im, &
!! optional real(real64) tol, &
!! optional integer(int32) maxits, &
!! optional integer(int32) iout, &
!! optional class(errors) err)
!! @endcode
!!
!! @param[in] a The original N-by-N matrix.
!! @param[in] lu The N-by-N LU-factored matrix of the approximation to the
!! system as output by @ref lu_factor.
!! @param[in] ju The N-element U row tracking array output by
!! @ref lu_factor.
!! @param[in,out] b On input, the N-element right-hand-side array. On
!! output, this array is overwritten as it is used as in-place storage
!! by the PGMRES algorithm.
!! @param[out] x The N-element solution array.
!! @param[in] im An optional parameter specifying the size of the Krylov
!! subspace. This value should not exceed 50.
!! @param[in] tol An optional parameter specifying the convergence tolerance
!! against which the Euclidean norm of the residual is checked. The
!! default value is the square root of machine precision.
!! @param[in] maxits An optional parameter specifying the maximum number
!! of iterations allowed. The default is 100.
!! @param[in] iout An optional parameter used to specify the device to
!! which status updates will be written. If no updates are requested,
!! a value less than or equal to zero should be supplied. The default
!! is zero such that no updates will be provided.
!! @param[in,out] err An optional errors-based object that if provided can
!! be used to retrieve information relating to any errors encountered
!! during execution. If not provided, a default implementation of the
!! errors class is used internally to provide error handling. Possible
!! errors and warning messages that may be encountered are as follows.
!! - LA_ARRAY_SIZE_ERROR: Occurs if any of the arrays and/or matrices are
!! not sized correctly.
!! - LA_OUT_OF_MEMORY_ERROR: Occurs if there is an issue with internal
!! memory allocations.
interface pgmres_solver
module procedure :: csr_pgmres_solver
end interface

interface
module function csr_get_element(this, i, j) result(rst)
class(csr_matrix), intent(in) :: this
Expand Down Expand Up @@ -5675,6 +5726,18 @@ module subroutine csr_lu_solve(lu, ju, b, x, err)
real(real64), intent(out), dimension(:) :: x
class(errors), intent(inout), optional, target :: err
end subroutine

module subroutine csr_pgmres_solver(a, lu, ju, b, x, im, tol, maxits, &
iout, err)
class(csr_matrix), intent(in) :: a
class(msr_matrix), intent(in) :: lu
integer(int32), intent(in), dimension(:) :: ju
real(real64), intent(inout), dimension(:) :: b
real(real64), intent(out), dimension(:) :: x
integer(int32), intent(in), optional :: im, maxits, iout
real(real64), intent(in), optional :: tol
class(errors), intent(inout), optional, target :: err
end subroutine
end interface

! ------------------------------------------------------------------------------
Expand Down
49 changes: 37 additions & 12 deletions src/linalg_sparse.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1329,12 +1329,12 @@ module subroutine csr_lu_solve(lu, ju, b, x, err)
! ******************************************************************************
! ITERATIVE SOLVERS
! ------------------------------------------------------------------------------
subroutine pgmres_solver(a, lu, ju, b, x, im, tol, maxits, iout, err)
module subroutine csr_pgmres_solver(a, lu, ju, b, x, im, tol, maxits, iout, err)
! Arguments
class(csr_matrix), intent(in) :: a
class(msr_matrix), intent(in) :: lu
integer(int32), intent(in), dimension(:) :: ju
real(real64), intent(in), dimension(:) :: b
real(real64), intent(inout), dimension(:) :: b
real(real64), intent(out), dimension(:) :: x
integer(int32), intent(in), optional :: im, maxits, iout
real(real64), intent(in), optional :: tol
Expand All @@ -1348,6 +1348,7 @@ subroutine pgmres_solver(a, lu, ju, b, x, im, tol, maxits, iout, err)
type(errors), target :: deferr

! Initialization
n = size(a, 1)
if (present(err)) then
errmgr => err
else
Expand All @@ -1356,7 +1357,7 @@ subroutine pgmres_solver(a, lu, ju, b, x, im, tol, maxits, iout, err)
if (present(im)) then
krylov = im
else
krylov = 10
krylov = min(n, 50)
end if
if (present(tol)) then
eps = tol
Expand All @@ -1373,40 +1374,64 @@ subroutine pgmres_solver(a, lu, ju, b, x, im, tol, maxits, iout, err)
else
io = 0
end if
n = size(a, 1)

! Input Checking
if (size(a, 2) /= n) then
! ERROR - input not square
call errmgr%report_error("csr_pgmres_solver", &
"The input matrix is not square.", LA_ARRAY_SIZE_ERROR)
return
end if
if (size(lu, 1) /= n .or. size(lu, 2) /= n) then
! ERROR: LU factored matrix incompatible size
call errmgr%report_error("csr_pgmres_solver", &
"The LU factored matrix size is not compatible with the " // &
"input matrix.", LA_ARRAY_SIZE_ERROR)
return
end if
if (size(b) /= n) then
! ERROR: RHS not sized correctly.
call errmgr%report_error("csr_pgmres_solver", &
"The output array dimension does not match the rest of the problem.", &
LA_ARRAY_SIZE_ERROR)
return
end if
if (size(x) /= n) then
! ERROR: Solution not sized correctly.
call errmgr%report_error("csr_pgmres_solver", &
"Inner matrix dimension mismatch.", LA_ARRAY_SIZE_ERROR)
return
end if
if (eps < epsilon(eps)) then
! ERROR: Tolerance too small
call errmgr%report_error("csr_pgmres_solver", &
"The convergence tolerance is too small.", LA_INVALID_INPUT_ERROR)
return
end if
if (mit < 1) then
! ERROR: Too few iterations allowed
call errmgr%report_error("csr_pgmres_solver", &
"Too few iterations allowed.", LA_INVALID_INPUT_ERROR)
return
end if
if (krylov < 1) then
! ERROR: Krylov subspace size too small
call errmgr%report_error("csr_pgmres_solver", &
"The requested Krylov subspace size is too small.", &
LA_INVALID_INPUT_ERROR)
return
end if

! Memory Allocation
allocate(vv(n,krylov+1), stat = flag)
if (flag /= 0) then
! ERROR: Memory allocation issue
call errmgr%report_error("csr_pgmres_solver", &
"Memory allocation error.", LA_OUT_OF_MEMORY_ERROR)
return
end if

! Process
call pgmres(n, krylov, b, x, vv, eps, mit, io, a%values, a%column_indices, &
a%row_indices, lu%values, lu%indices, ju, ierr)
if (ierr == 1) then
call errmgr%report_error("csr_pgmres_solver", &
"Convergence could not be achieved to the requested tolerance " // &
"in the allowed number of iterations.", LA_CONVERGENCE_ERROR)
return
end if
end subroutine

! ------------------------------------------------------------------------------
Expand Down
3 changes: 3 additions & 0 deletions tests/linalg_test.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -351,6 +351,9 @@ program main
rst = test_csr_lu_factor_1()
if (.not.rst) flag = 105

rst = test_pgmres_1()
if (.not.rst) flag = 106

! End
if (flag /= 0) stop flag
end program
45 changes: 45 additions & 0 deletions tests/test_sparse.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -631,5 +631,50 @@ function test_csr_lu_factor_1() result(rst)
end if
end function

! ------------------------------------------------------------------------------
function test_pgmres_1() result(rst)
! Arguments
logical :: rst

! Local Variables
integer(int32) :: ju(4), ipiv(4)
real(real64) :: dense(4, 4), m(4, 4), b(4), bc(4), x(4)
type(csr_matrix) :: a, am
type(msr_matrix) :: lu

! Initialization
rst = .true.
dense = reshape([ &
5.0d0, 0.0d0, 0.0d0, 0.0d0, &
0.0d0, 8.0d0, 0.0d0, 6.0d0, &
0.0d0, 0.0d0, 3.0d0, 0.0d0, &
0.0d0, 0.0d0, 0.0d0, 5.0d0], [4, 4])
m = reshape([ &
5.0d0, 0.0d0, 0.0d0, 0.0d0, &
0.0d0, 8.0d0, 0.0d0, 0.0d0, &
0.0d0, 0.0d0, 3.0d0, 0.0d0, &
0.0d0, 0.0d0, 0.0d0, 5.0d0], [4, 4])
a = dense
am = m
call random_number(b)
bc = b

! Compute the preconditioner
call lu_factor(am, lu, ju)

! Solve the sparse system
call pgmres_solver(a, lu, ju, b, x)

! Solve the dense system directly
call lu_factor(dense, ipiv)
call solve_lu(dense, ipiv, bc)

! Test
if (.not.assert(x, bc)) then
rst = .false.
print "(A)", "Test Failed: test_pgmres_1 -1"
end if
end function

! ------------------------------------------------------------------------------
end module

0 comments on commit 3d03dfc

Please sign in to comment.