From 89ffc2f6b8fc127353151d34406e147b2f3085a2 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Wed, 24 Jan 2024 23:00:44 -0600 Subject: [PATCH] added the LSMR solver as an option Experimental See #16 --- fpm.toml | 1 + src/nlesolver_module.F90 | 32 ++++++++++++++++++++++++++++++++ test/sparse_test.f90 | 30 ++++++++++++++++++++++-------- 3 files changed, 55 insertions(+), 8 deletions(-) diff --git a/fpm.toml b/fpm.toml index d530b11..135d697 100644 --- a/fpm.toml +++ b/fpm.toml @@ -19,6 +19,7 @@ source-dir = "src" fmin = { git="https://github.com/jacobwilliams/fmin.git", tag="1.1.1" } LSQR = { git="https://github.com/jacobwilliams/LSQR", tag="1.1.0" } lusol = { git="https://github.com/jacobwilliams/lusol", tag="1.0.0" } +LSMR = { git="https://github.com/jacobwilliams/LSMR" } [install] library = true diff --git a/src/nlesolver_module.F90 b/src/nlesolver_module.F90 index 1de8926..f57b97f 100755 --- a/src/nlesolver_module.F90 +++ b/src/nlesolver_module.F90 @@ -42,6 +42,7 @@ module nlesolver_module use fmin_module, only: fmin use lsqr_module, only: lsqr_solver_ez use lusol_ez_module, only: solve, lusol_settings + use lsmrModule, only: lsmr_ez implicit none @@ -109,6 +110,7 @@ module nlesolver_module !! * 1 - assume dense (use dense solver) !! * 2 - assume sparse (use LSQR sparse solver). !! * 3 - assume sparse (use LUSOL sparse solver). + !! * 4 - assume sparse (use LSMR sparse solver). integer :: n_nonzeros = -1 !! number of nonzero Jacobian elements (used for `sparsity_mode > 1`) integer,dimension(:),allocatable :: irow !! sparsity pattern nonzero elements row indices. integer,dimension(:),allocatable :: icol !! sparsity pattern nonzero elements column indices @@ -366,6 +368,7 @@ subroutine initialize_nlesolver_variables(me,& !! Must specify `grad` for this mode. !! * 2 - assume sparse (use LSQR sparse solver). !! * 3 - assume sparse (use LUSOL sparse solver). + !! * 4 - assume sparse (use LSMR sparse solver). !! Must also specify `grad_sparse` and the `irow` and `icol` !! sparsity pattern for `mode>1`. integer,dimension(:),intent(in),optional :: irow !! sparsity pattern nonzero elements row indices. @@ -549,6 +552,8 @@ subroutine nlesolver_solver(me,x) integer :: i !! counter integer,dimension(:),allocatable :: idx, index_array !! for sparse indexing character(len=10) :: i_str !! string version of `i` for row string + real(wp) :: normA, condA, normr, normAr, normx !! for LSMR + integer :: localSize, itn !! for LSMR if (me%istat<0) return ! class was not initialized properly @@ -746,6 +751,33 @@ subroutine nlesolver_solver(me,x) lusol_options%method = me%lusol_method call solve(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info,& settings=lusol_options) + case (4) + + ! use LSMR solver + + ! TODO this should be an input + localSize = 0 + !localSize = min(me%m, me%n) + + me%conlim = 1.0_wp/(100*epsilon(1.0_wp)) + call lsmr_ez ( me%m, me%n, me%irow, me%icol, fjac_sparse, rhs, me%damp, & + me%atol, me%btol, me%conlim, me%itnlim, localSize, me%nout, & + p, info, itn, normA, condA, normr, normAr, normx ) + + ! check convergence: + select case (info) + case(4) + call me%set_status(istat = -1004, & + string = 'LSMR Error: The system appears to be ill-conditioned. istop =', i=info) + exit + case(5) + call me%set_status(istat = -1005, & + string = 'LSMR Error: The iteration limit was reached. istop =', i=info) + exit + case default + info = 0 + end select + end select ! check for errors: diff --git a/test/sparse_test.f90 b/test/sparse_test.f90 index baf4dbc..5fdc65d 100644 --- a/test/sparse_test.f90 +++ b/test/sparse_test.f90 @@ -13,6 +13,22 @@ program sparse_test integer,parameter :: max_iter = 100 real(wp),parameter :: tol = 1.0e-8_wp logical,parameter :: verbose = .false. + integer,dimension(3),parameter :: icol = [1,2,2] + integer,dimension(3),parameter :: irow = [1,1,2] + + integer :: f_evals + + call go(2, 'LSQR') + call go(3, 'LUSOL') + call go(4, 'LSMR') + + contains + + subroutine go(sparsity_mode, mode_name) + !! run the tests for the specified sparsity mode + + integer,intent(in) :: sparsity_mode + character(len=*),intent(in) :: mode_name !! name of the sparsity mode used type(nlesolver_type) :: solver real(wp) :: alpha @@ -22,21 +38,17 @@ program sparse_test integer :: istat !! Integer status code. character(len=:),allocatable :: message !! Text status message real(wp),dimension(n) :: x - integer :: f_evals integer :: i character(len=:),allocatable :: description real(wp) :: fmin_tol - integer,dimension(3),parameter :: icol = [1,2,2] - integer,dimension(3),parameter :: irow = [1,1,2] - fmin_tol = 1.0e-2_wp ! don't need a tight tol for this n_intervals = 2 alpha = 1.0_wp write(*,*) '' write(*,*) '***********************' - write(*,*) '* sparse_test *' + write(*,*) '* sparse_test : '//trim(mode_name) write(*,*) '***********************' write(*,*) '' do i = 1, 8 @@ -94,6 +106,8 @@ program sparse_test m = m, & max_iter = max_iter, & tol = tol, & + atol = tol, & + btol = tol, & func = func, & grad_sparse = grad_sparse, & step_mode = step_mode,& @@ -102,10 +116,10 @@ program sparse_test n_intervals = n_intervals, & fmin_tol = fmin_tol, & verbose = verbose,& - sparsity_mode = 3,& + sparsity_mode = sparsity_mode,& ! lsmr irow = irow,& icol = icol,& - damp = 1.0_wp) + damp = 0.0_wp) call solver%status(istat, message) write(*,'(I3,1X,A)') istat, message if (istat /= 0) error stop @@ -119,7 +133,7 @@ program sparse_test end do - contains + end subroutine go subroutine func(me,x,f) !! compute the function