From b48969f7fe69b380a68da2e3e490fff4b7565e65 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 26 Jan 2026 09:34:16 +0000 Subject: [PATCH 01/16] added linear solver module --- libRALFit/src/CMakeLists.txt | 5 ++ libRALFit/src/ral_nlls_linear.F90 | 135 ++++++++++++++++++++++++++++++ 2 files changed, 140 insertions(+) create mode 100644 libRALFit/src/ral_nlls_linear.F90 diff --git a/libRALFit/src/CMakeLists.txt b/libRALFit/src/CMakeLists.txt index e3b93242..23e48031 100644 --- a/libRALFit/src/CMakeLists.txt +++ b/libRALFit/src/CMakeLists.txt @@ -14,6 +14,7 @@ set(source_files nag_export_mod.F90 ral_nlls_fd.F90 ral_nlls_internal.F90 + ral_nlls_linear.F90 ral_nlls.F90 ral_nlls_ciface.F90 CACHE INTERNAL "source_files" FORCE) @@ -92,6 +93,8 @@ function(generate_ral_nlls_modules precision obj_list) "NAG_EXPORT_MOD_${precision}" "RAL_NLLS_WORKSPACES_${precision}") add_dependencies("RAL_NLLS_${precision}" "RAL_NLLS_INTERNAL_${precision}") add_dependencies("RAL_NLLS_CIFACE_${precision}" "RAL_NLLS_${precision}") + add_dependencies("RAL_NLLS_INTERNAL_${precision}" "RAL_NLLS_LINEAR_${precision}") + add_dependencies("RAL_NLLS_LINEAR_${precision}" "RAL_NLLS_WORKSPACES_${precision}") # Replicate tree to be compatible with Ninja target_link_libraries("RAL_NLLS_DTRS_${precision}" "RAL_NLLS_TYPES_${precision}" @@ -113,6 +116,8 @@ function(generate_ral_nlls_modules precision obj_list) "NAG_EXPORT_MOD_${precision}" "RAL_NLLS_WORKSPACES_${precision}") target_link_libraries("RAL_NLLS_${precision}" "RAL_NLLS_INTERNAL_${precision}") target_link_libraries("RAL_NLLS_CIFACE_${precision}" "RAL_NLLS_${precision}") + target_link_libraries("RAL_NLLS_INTERNAL_${precision}" "RAL_NLLS_LINEAR_${precision}") + target_link_libraries("RAL_NLLS_LINEAR_${precision}" "RAL_NLLS_WORKSPACES_${precision}") endfunction() # Post ################ diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 new file mode 100644 index 00000000..3fb2e108 --- /dev/null +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -0,0 +1,135 @@ +! ral_nlls_linear :: linear solvers for internal use in RALFit + +#include "preprocessor.FPP" + +module MODULE_PREC(ral_nlls_linear) + use MODULE_PREC(ral_nlls_workspaces), only: wp, solve_LLS_work, NLLS_inform, & + NLLS_options, NLLS_ERROR_WORKSPACE_ERROR, & + NLLS_ERROR_FROM_EXTERNAL + implicit none + + private + public :: solve_LLS, solve_LLS_nocopy + +contains + subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) +! ----------------------------------------------------------------- +! solve_LLS, a subroutine to solve a linear least squares problem +! ----------------------------------------------------------------- +! Solves the linear least squares problem Ax = b +! using the method specified in options%lls_solver +! Input: +! A - LHS matrix of the least squares problem (m x n) +! b - The RHS, overwritten with x on output (m x 1) +! n - Number of columns in A +! m - Number of rows in A +! inform - NLLS_inform structure +! options - NLLS_options structure +! w - workspace structure +! pd - logical, true if A is known to be positive definite +! Output: +! A_out - output copy of the original A matrix (m x n) +! x - solution vector (n x 1) +! inform - NLLS_inform structure +! work - updated workspace structure +! ----------------------------------------------------------------- + implicit none + real(wp), dimension(:), intent(in) :: A + real(wp), dimension(:), intent(in) :: b + real(wp), dimension(:), intent(out) :: A_out + real(wp), dimension(:), intent(out) :: x + integer, intent(in) :: n, m + type(NLLS_inform), intent(inout) :: inform + type(NLLS_options), Intent(In) :: options + type(solve_LLS_work), intent(inout) :: w + logical, intent(in) :: pd + + A_out = A + x = b + + call solve_LLS_nocopy(A_out, x, n, m, inform, w, options, pd) + + end subroutine solve_LLS + + subroutine solve_LLS_nocopy(A, b, n, m, inform, w, options, pd) +! linear solver core which overwrites A and b. + implicit none + real(wp), dimension(:), intent(inout) :: A + real(wp), dimension(:), intent(inout) :: b + integer, intent(in) :: n, m + type(NLLS_inform), intent(inout) :: inform + type(NLLS_options), intent(in) :: options + type(solve_LLS_work), intent(inout) :: w + logical, intent(in) :: pd + + select case (options%lls_solver) + case (1) ! LAPACK + if (pd) then + call solve_posv(A,b,n,inform,options) + else + if (.not. w%allocated) then + inform%status = NLLS_ERROR_WORKSPACE_ERROR + goto 100 + end if + call solve_gels(A,b,n,m,inform,w,options) + end if + end select + + if (inform%external_return /= 0 ) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + end if + +100 continue + end subroutine solve_LLS_nocopy + + subroutine solve_gels(A,b,n,m,inform,w,options) +! Wrapper around LAPACK's ?gels + implicit none + REAL(wp), DIMENSION(:), INTENT(IN) :: A + REAL(wp), DIMENSION(:), INTENT(IN) :: b + INTEGER, INTENT(IN) :: n, m + type(NLLS_inform), INTENT(INOUT) :: inform + type(NLLS_options), Intent(In) :: options + + integer, Parameter :: nrhs = 1 + integer :: lwork, lda, ldb + type( solve_LLS_work ), Intent(inout) :: w + + lwork = size(w%work) + + if (options%Fortran_Jacobian) then + lda = m + ldb = max(m,n) + call PREC(gels)('N', m, n, nrhs, A, lda, w%temp, ldb, w%work, lwork, & + inform%external_return) + else + lda = n + ldb = max(m,n) + call PREC(gels)('T', n, m, nrhs, A, lda, w%temp, ldb, w%work, lwork, & + inform%external_return) + end if + if (inform%external_return /= 0 ) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?gels' + end if + + end subroutine solve_gels + + subroutine solve_posv(A,b,n,inform,options) +! Wrapper around LAPACK's ?posv for positive definite systems + implicit none + REAL(wp), DIMENSION(:), INTENT(IN) :: A + REAL(wp), DIMENSION(:), INTENT(IN) :: b + INTEGER, INTENT(IN) :: n + type(NLLS_inform), INTENT(INOUT) :: inform + type(NLLS_options), Intent(In) :: options + + call PREC(posv)('N', n, 1, A, n, b, n, inform%external_return) + if (inform%external_return /= 0 ) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?posv' + end if + + end subroutine solve_posv + +end module MODULE_PREC(ral_nlls_linear) From eb7bbfa1a24477db341aa5925503b779a7d34726 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 26 Jan 2026 09:36:20 +0000 Subject: [PATCH 02/16] changed dogleg to use solver module --- libRALFit/src/ral_nlls_internal.F90 | 56 +++-------------------------- libRALFit/test/unit_test_mod.F90 | 26 +++++++++----- 2 files changed, 23 insertions(+), 59 deletions(-) diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index 396895c4..e3fb0789 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -14,6 +14,7 @@ module MODULE_PREC(ral_nlls_internal) Use MODULE_PREC(ral_nlls_printing) Use MODULE_PREC(ral_nlls_bounds) Use MODULE_PREC(ral_nlls_types), Only: PREC(nrm2), PREC(dot) + use MODULE_PREC(ral_nlls_linear), only: solve_LLS, solve_LLS_nocopy implicit none @@ -28,7 +29,7 @@ module MODULE_PREC(ral_nlls_internal) public :: calculate_step, evaluate_model public :: update_trust_region_radius, apply_second_order_info, rank_one_update public :: test_convergence, calculate_rho - public :: solve_LLS, shift_matrix + public :: shift_matrix public :: dogleg, more_sorensen, generate_scaling, solve_newton_tensor, aint_tr public :: switch_to_quasi_newton, evaltensor_J public :: setup_iparams_type, free_iparams_type @@ -1802,9 +1803,10 @@ SUBROUTINE dogleg(J,f,hf,g,n,m,Delta,d,normd,options,inform,w) ! Solve the linear problem... select case (options%model) case (1) - ! linear model... - call solve_LLS(J,f,n,m,w%d_gn,inform,w%solve_LLS_ws,options%Fortran_Jacobian) + ! linear model... + call solve_LLS(J,f,w%solve_LLS_ws%Jlls,w%d_gn,n,m,inform,w%solve_LLS_ws,options,.false.) if ( inform%status /= 0 ) goto 100 + w%d_gn(:) = -w%d_gn(:) case default inform%status = NLLS_ERROR_DOGLEG_MODEL goto 100 @@ -2924,54 +2926,6 @@ subroutine regularization_solver(A,v,n,m,Delta,num_successful_steps,& End If end subroutine regularization_solver - SUBROUTINE solve_LLS(J,f,n,m,d_gn,inform,w,Fortran_Jacobian) -! ----------------------------------------------------------------- -! solve_LLS, a subroutine to solve a linear least squares problem -! ----------------------------------------------------------------- - - Implicit None - REAL(wp), DIMENSION(:), INTENT(IN) :: J - REAL(wp), DIMENSION(:), INTENT(IN) :: f - INTEGER, INTENT(IN) :: n, m - REAL(wp), DIMENSION(:), INTENT(OUT) :: d_gn - type(NLLS_inform), INTENT(INOUT) :: inform - logical, Intent(In) :: Fortran_Jacobian - - integer, Parameter :: nrhs = 1 - integer :: lwork, lda, ldb - type( solve_LLS_work ), Intent(inout) :: w - - If (.not. w%allocated) Then - inform%status = NLLS_ERROR_WORKSPACE_ERROR - goto 100 - End If - - w%temp(1:m) = f(1:m) - lwork = size(w%work) - - w%Jlls(:) = J(:) - If (Fortran_Jacobian) Then - lda = m - ldb = max(m,n) - call PREC(gels)('N', m, n, nrhs, w%Jlls, lda, w%temp, ldb, w%work, lwork, & - inform%external_return) - else - lda = n - ldb = max(m,n) - call PREC(gels)('T', n, m, nrhs, w%Jlls, lda, w%temp, ldb, w%work, lwork, & - inform%external_return) - end If - if (inform%external_return .ne. 0 ) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?gels' - Go To 100 - end if - - d_gn = -w%temp(1:n) - -100 continue - END SUBROUTINE solve_LLS - SUBROUTINE findbeta(a, b, Delta, beta, inform) ! ----------------------------------------------------------------- diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index cdaace7d..bcf043ce 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -10,6 +10,7 @@ module MODULE_PREC(unit_test_mod) use MODULE_PREC(ral_nlls) use MODULE_PREC(ral_nlls_internal) + use MODULE_PREC(ral_nlls_linear) use MODULE_PREC(ral_nlls_workspaces) implicit none @@ -2177,8 +2178,10 @@ subroutine solve_LLS_tests(options,fails) 6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp ] f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] - call solve_LLS(J,f,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.True.) + call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + d,n,m,status, & + work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + d(:) = -d(:) if ( status%status .ne. 0 ) then write (*, *) '[id:1] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2204,8 +2207,12 @@ subroutine solve_LLS_tests(options,fails) 8.0_wp, 4.0_wp, 9.0_wp, 5.0_wp, 10.0_wp ] f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] - call solve_LLS(J,f,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.False.) + + call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + d,n,m,status, & + work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + d(:) = -d(:) + if ( status%status .ne. 0 ) then write (*, *) '[id:2] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2233,8 +2240,10 @@ subroutine solve_LLS_tests(options,fails) J = 1.0_wp J(1:m) = 0.0_wp f = 1.0_wp - call solve_LLS(J,f,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.True.) + call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + d,n,m,status, & + work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + if ( status%status .ne. NLLS_ERROR_FROM_EXTERNAL ) then write (*, *) '[id:3] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2244,8 +2253,9 @@ subroutine solve_LLS_tests(options,fails) call nlls_finalize(work, options, status) - call solve_LLS(J,f,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.True.) + call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + d,n,m,status, & + work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if (status%status .ne. NLLS_ERROR_WORKSPACE_ERROR) then write(*,*) 'Error: workspace error not flagged when workspaces not setup' fails = fails + 1 From 9bf2fd8bd9200a9ba38a23f8c79c4d2d5d88850b Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 26 Jan 2026 13:32:28 +0000 Subject: [PATCH 03/16] integrated minus_solve functionalities into solve_LLS --- libRALFit/log | 0 libRALFit/src/ral_nlls_internal.F90 | 234 ++++++++------------ libRALFit/src/ral_nlls_linear.F90 | 91 +++++--- libRALFit/src/ral_nlls_workspaces.F90 | 54 +---- libRALFit/test/nlls_test.F90 | 14 +- libRALFit/test/unit_test_mod.F90 | 293 +++++++++++++------------- 6 files changed, 307 insertions(+), 379 deletions(-) create mode 100644 libRALFit/log diff --git a/libRALFit/log b/libRALFit/log new file mode 100644 index 00000000..e69de29b diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index e3fb0789..ac4422fe 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -23,7 +23,6 @@ module MODULE_PREC(ral_nlls_internal) public :: nlls_solve, nlls_iterate, nlls_finalize, nlls_strerror public :: solve_galahad, findbeta, mult_j public :: mult_jt, matmult_inner - public :: minus_solve_spd, minus_solve_spd_nocopy, minus_solve_general public :: matmult_outer, outer_product, min_eig_symm, max_eig, all_eig_symm public :: remove_workspaces, setup_workspaces public :: calculate_step, evaluate_model @@ -1461,7 +1460,6 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ options, inform) end if subproblem_method = 2 ! try aint next - call dogleg(J,f,hf,g,n,m,Delta,d,norm_S_d, & options,inform,w%dogleg_ws) if (inform%status == 0) subproblem_success = .true. @@ -1803,10 +1801,12 @@ SUBROUTINE dogleg(J,f,hf,g,n,m,Delta,d,normd,options,inform,w) ! Solve the linear problem... select case (options%model) case (1) - ! linear model... - call solve_LLS(J,f,w%solve_LLS_ws%Jlls,w%d_gn,n,m,inform,w%solve_LLS_ws,options,.false.) + ! linear model... + w%solve_LLS_ws%Jlls = reshape(J, [n,m]) + w%solve_LLS_ws%temp(1:m) = f(1:m) + call solve_LLS_nocopy(w%solve_LLS_ws%Jlls,w%solve_LLS_ws%temp,n,m,inform,w%solve_LLS_ws,options,.false.) if ( inform%status /= 0 ) goto 100 - w%d_gn(:) = -w%d_gn(:) + w%d_gn = -w%solve_LLS_ws%temp(1:n) case default inform%status = NLLS_ERROR_DOGLEG_MODEL goto 100 @@ -1869,6 +1869,9 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) real(wp) :: obj_p0, obj_p1, obj_p0_gn, obj_p1_gn REAL(wp) :: norm_p0, tau, lam, eta Character(Len=80) :: rec(1) + logical :: pd ! whether A is positive-definite + + pd = (options%model == 1) If ( .not. w%allocated ) Then inform%status = NLLS_ERROR_WORKSPACE_ERROR @@ -1894,15 +1897,11 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) w%B(i,i) = 1.0_wp end do - select case (options%model) - case (1) - call minus_solve_spd(A,v,w%LtL,w%p0,n,inform) - if (inform%status /= 0) goto 100 - case default - ! note: v is mult by -1 in minus_solve_general - call minus_solve_general(A,v,w%p0,n,inform,w%minus_solve_general_ws) - if (inform%status /= 0) goto 100 - end select + w%LtL(:, :) = A(:, :) + w%p0(:) = -v(:) + call solve_LLS_nocopy(w%LtL,w%p0,n,n,inform,w%solve_LLS_ws,options,pd) + if (inform%status /= 0) goto 100 + call matrix_norm(w%p0,w%B,norm_p0) @@ -1959,19 +1958,9 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) call matmult_outer( w%By_hardcase, size_hard2, n, w%M1_small) w%M0_small(:,:) = A(:,:) + lam*w%B(:,:) + w%M1_small ! solve Hq + g = 0 for q - select case (options%model) - case (1) - ! note: v is mult by -1 in minus_solve_spd - call minus_solve_spd(w%M0_small,v,w%LtL,w%q,n,inform) - if (inform%status /= 0) goto 100 - case default - ! note: v is mult by -1 in minus_solve_general - call minus_solve_general(w%M0_small,v,w%q,n,inform,w%minus_solve_general_ws) - if (inform%status /= 0) goto 100 - end select - ! note -- a copy of the matrix is taken on entry to the solve routines - ! (I think..) and inside...fix - + w%q(:) = -v(:) + call solve_LLS_nocopy(w%M0_small, w%q, n, n, inform, w%solve_LLS_ws, options, pd) + if (inform%status /= 0) goto 100 ! find max eta st ||q + eta v(:,1)||_B = Delta call findbeta(w%q,w%y_hardcase(:,1),Delta,eta,inform) @@ -1984,25 +1973,15 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) w%p1(:) = w%q(:) + eta * w%y_hardcase(:,1) else - ! Solve (A+lam*w%B)x=-v: - ! 1. we can use w%B to actually store A + lam * w%B - ! 2. b=v is then inverted inside of `minus_solve_spd` or `minus_solve_general` + ! Solve (A+lam*B)x=-v: + ! 1. we can use w%B to actually store A + lam * B w%B(:,:) = A(:,:) + w%p1(:) = -v(:) Do i = 1, n w%B(i,i) = w%B(i,i) + lam End Do - select case (options%model) - case (1) - ! note: v is mult by -1 in minus_solve_spd - call minus_solve_spd(w%B,v,w%LtL,w%p1,n,inform) - if (inform%status /= 0) goto 100 - case default - ! note: v is mult by -1 in minus_solve_general - call minus_solve_general(w%B,v,w%p1,n,inform,w%minus_solve_general_ws) - if (inform%status /= 0) goto 100 - end select - ! note -- a copy of the matrix is taken on entry to the solve routines - ! and inside...fix + call solve_LLS_nocopy(w%B, w%p1, n, n, inform, w%solve_LLS_ws, options, pd) + if (inform%status /= 0) goto 100 end if ! get obj_p1: the value of the model at p1 @@ -2089,6 +2068,13 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) integer :: i, no_restarts Character(Len=80) :: rec(2) + ! reset error calls from previous method + ! else a previous method failure will force + ! min_eig_symm no matter what + inform%status = 0 + inform%external_return = 0 + inform%external_name = '' + ! The code finds ! d = arg min_p v^T p + 0.5 * p^T A p ! s.t. ||p|| \leq Delta @@ -2107,11 +2093,12 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) w%AplusSigma(1:n,1:n) = A(1:n,1:n) If (options%force_min_eig_symm) Then - ! Skip minus_solve_spd_nocopy and jump directly to min_eig_symm + ! Skip solve and jump directly to min_eig_symm inform%status = 1 Else - ! note: v is mult by -1 in minus_solve_spd_nocopy - call minus_solve_spd_nocopy(w%AplusSigma,v,d,n,inform) + ! solve (A+sigma)d = -v + d(1:n) = -v(1:n) + call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) End If if (inform%status == 0) then ! A is symmetric positive definite.... @@ -2121,6 +2108,10 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) Call Printmsg(5,.False.,options,1,rec) End If else + If (buildmsg(5,.False.,options)) Then + Write(rec(1), Fmt=6010) inform%external_return + Call Printmsg(5,.False.,options,1,rec) + End If ! reset the error calls -- handled in the code.... inform%status = 0 inform%external_return = 0 @@ -2128,10 +2119,6 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) call min_eig_symm(A,n,sigma,w%y1,options,inform,w%min_eig_symm_ws) if (inform%status /= 0) goto 100 sigma = -(sigma - local_ms_shift) - If (buildmsg(5,.False.,options)) Then - Write(rec(1), Fmt=6010) - Call Printmsg(5,.False.,options,1,rec) - End If ! find a shift that makes (A + sigma I) positive definite, ! and solve (A + sigma I) (-v) = d call check_shift_and_solve(n,A,w%AplusSigma,v,sigma,d,options,inform,w) @@ -2225,7 +2212,7 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) Write(rec(1), Fmt=5040) Call Printmsg(5,.False.,options,1,rec) End If - goto 100 + goto 100 end if w%q(:) = d ! w%q = R'\d @@ -2257,8 +2244,9 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) End If sigma = sigma + sigma_shift call shift_matrix(A,sigma,w%AplusSigma,n) - ! note: v is mult by -1 in minus_solve_spd_nocopy - call minus_solve_spd_nocopy(w%AplusSigma,v,d,n,inform) + ! solve (A + sigma)d = -v + d(1:n) = -v(1:n) + call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) end if if (inform%status /= 0) goto 100 @@ -2281,7 +2269,7 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) 5040 FORMAT('Leaving More-Sorensen') ! print_level >= 3 6000 FORMAT('A is symmetric positive definite') -6010 FORMAT('Trying a shift of sigma = ',ES12.4) +6010 FORMAT('A is not symmetric PD: non-PD principal minor ',i4) 6020 FORMAT('A + sigma I is symmetric positive definite') 6030 FORMAT('We''re within the trust region radius initially') 6035 FORMAT('We''re within the trust region radius') @@ -2389,8 +2377,10 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) ! d = -A\v if (.not. factorization_done) then call shift_matrix(A,sigma,w%AplusSigma,n) - ! note: v is mult by -1 in minus_solve_spd - call minus_solve_spd(w%AplusSigma,v,w%LtL,d,n,inform) + ! solve (A + sigma)d = -v + w%LtL(1:n,1:n) = w%AplusSigma(1:n,1:n) + d(1:n) = -v(1:n) + call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) else factorization_done = .false. end if @@ -2437,10 +2427,6 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) Call printmsg(5,.False.,options,1,rec) End If location_of_breakdown = inform%external_return - ! reset the error calls -- handled in the code.... - inform%status = 0 - inform%external_return = 0 - inform%external_name = REPEAT( ' ', 80 ) region = 1 ! N region_char = 'N' sigma_l = sigma @@ -2449,6 +2435,10 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) Call printmsg(5,.False.,options,1,rec) End If end if + ! reset the error calls -- handled in the code.... + inform%status = 0 + inform%external_return = 0 + inform%external_name = REPEAT( ' ', 80 ) if (region .ne. 1 ) then ! in F w%q(:) = d ! w%q = R'\d @@ -2493,7 +2483,8 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) Call printmsg(5,.False.,options,1,rec) End If end if - else ! not in F + else if (location_of_breakdown > 0) then + ! not in F ! find an alpha and y1 such that ! (A + sigmaI + alpha e_k e_k^T)v = 0 !!$ write(*,*) 'location of breakdown = ', location_of_breakdown @@ -2520,6 +2511,11 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) end do nrmy1 = PREC(nrm2)(n,w%y1,1) ! w%y1(1:n) sigma_l = max(sigma_l, sigma + indef_delta/nrmy1) + else + ! `dposv` found internal error (this is a developer error!) + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_return = location_of_breakdown + goto 100 end if ! check for termination @@ -2560,8 +2556,11 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) End If elseif (region == 3) then call shift_matrix(A,sigma + sigma_shift,w%AplusSigma,n) - ! note: v is mult by -1 in minus_solve_spd - call minus_solve_spd(w%AplusSigma,v,w%LtL,d,n,inform) + ! solve (A + sigma)d = -v + w%LtL(1:n,1:n) = w%AplusSigma(1:n,1:n) + d(1:n) = -v(1:n) + call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) + if (inform%status == 0) then region = 2 sigma = sigma + sigma_shift @@ -2593,6 +2592,10 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) Write(rec(1), Fmt=5010) i, region_char, nd, sigma_l, sigma, sigma_u Call printmsg(5,.False.,options,1,rec) End If + ! reset the error calls -- handled in the code.... + inform%status = 0 + inform%external_return = 0 + inform%external_name = REPEAT( ' ', 80 ) end do @@ -2724,14 +2727,17 @@ subroutine check_shift_and_solve(n,A,AplusSigma,v,sigma,d,options,inform,w) no_shifts = 0 successful_shift = .false. do while( .not. successful_shift ) + If (buildmsg(5,.False.,options)) Then + Write(rec(1), Fmt=99999) sigma + Call Printmsg(5,.False.,options,1,rec) + end if call shift_matrix(A,sigma,w%AplusSigma,n) - ! note: v is mult by -1 in minus_solve_spd_nocopy - call minus_solve_spd_nocopy(w%AplusSigma,v,d,n,inform) - if ( inform%status /= 0 ) then - ! reset the error calls -- handled in the code.... - inform%status = 0 - inform%external_return = 0 - inform%external_name = REPEAT( ' ', 80 ) + ! solve (A + sigma)d = -v + d(:) = -v(:) + call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + if ( inform%status == 0 ) then + successful_shift = .true. + else no_shifts = no_shifts + 1 If ( no_shifts >= 10 ) Then inform%status = NLLS_ERROR_MS_TOO_MANY_SHIFTS @@ -2739,15 +2745,18 @@ subroutine check_shift_and_solve(n,A,AplusSigma,v,sigma,d,options,inform,w) End If sigma = sigma + (10.0_wp**no_shifts) * options%more_sorensen_shift If (buildmsg(5,.False.,options)) Then - Write(rec(1), Fmt=99999) sigma + Write(rec(1), Fmt=99998) inform%external_return Call Printmsg(5,.False.,options,1,rec) - End If - else - successful_shift = .true. + end if + ! reset the error calls -- handled in the code.... + inform%status = 0 + inform%external_return = 0 + inform%external_name = REPEAT( ' ', 80 ) end if end do 100 continue +99998 FORMAT('A + sigma I has non-positive principal minor', i4) 99999 FORMAT('Trying a shift of sigma = ',ES12.4) end subroutine check_shift_and_solve @@ -2889,7 +2898,7 @@ subroutine regularization_solver(A,v,n,m,Delta,num_successful_steps,& ! regularization_solver ! Solve the regularized subproblem using ! a home-rolled algorithm - ! Note: inform%status is set by minus_solve_spd + ! Note: inform%status is set by solve_LLS !-------------------------------------------- Implicit None @@ -2910,8 +2919,10 @@ subroutine regularization_solver(A,v,n,m,Delta,num_successful_steps,& if ( reg_order == 2.0_wp ) then call shift_matrix(A,reg_param,w%AplusSigma,n) - ! note: v is mult by -1 in minus_solve_spd - call minus_solve_spd(w%AplusSigma,v,w%LtL,d,n,inform) + ! solve (A+sigma)d = -v + w%LtL(1:n, 1:n) = w%AplusSigma(1:n,1:n) + d(1:n) = -v(1:n) + call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) ! informa%status is passed along, this routine exits here else ! Feature not yet implemented, this should have been caught in @@ -3435,6 +3446,7 @@ subroutine mult_Jt(J,n,m,x,Jtx,Fortran_Jacobian) End If end subroutine mult_Jt + subroutine scale_J_by_weights(J,n,m,weights,options) Implicit None real(wp), intent(inout) :: J(*) @@ -3533,76 +3545,6 @@ subroutine get_element_of_matrix(J,m,ii,jj,Jij) Jij = J(ii + (jj-1)*m) end subroutine get_element_of_matrix - subroutine minus_solve_spd(A,b,LtL,x,n,inform) - Implicit None - REAL(wp), intent(in), contiguous :: A(:,:) - REAL(wp), intent(in), contiguous :: b(:) - REAL(wp), intent(out), contiguous :: LtL(:,:) - REAL(wp), intent(out), contiguous :: x(:) - integer, intent(in) :: n - type( nlls_inform), intent(inout) :: inform - inform%status = 0 - ! Wrapper for the lapack subroutine ?posv - ! Given A and b, solves - ! A x = -b - LtL(1:n,1:n) = A(1:n,1:n) - x(1:n) = -b(1:n) - call PREC(posv)('L', n, 1, LtL, n, x, n, inform%external_return) - if (inform%external_return .ne. 0) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?posv' - end if - end subroutine minus_solve_spd - - subroutine minus_solve_spd_nocopy(A,b,x,n,inform) - Implicit None - REAL(wp), intent(inout), contiguous :: A(:,:) - REAL(wp), intent(in), contiguous :: b(:) - REAL(wp), intent(out), contiguous :: x(:) - integer, intent(in) :: n - type( nlls_inform), intent(inout) :: inform - inform%status = 0 - ! Wrapper for the lapack subroutine ?posv - ! NOTE: A will be destroyed - ! Given A and b, solves - ! A x = -b - x(1:n) = -b(1:n) - call PREC(posv)('L', n, 1, A, n, x, n, inform%external_return) - if (inform%external_return .ne. 0) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?posv' - end if - end subroutine minus_solve_spd_nocopy - - subroutine minus_solve_general(A,b,x,n,inform,w) - Implicit None - REAL(wp), intent(in), contiguous :: A(:,:) - REAL(wp), intent(in), contiguous :: b(:) - REAL(wp), intent(out), contiguous :: x(:) - integer, intent(in) :: n - type( nlls_inform ), intent(inout) :: inform - type( minus_solve_general_work ),Intent(inout) :: w - ! Wrapper for the lapack subroutine ?gesv - ! NOTE: A would be destroyed - ! Given A and b, solves - ! A x = -b - - If (.not. w%allocated) Then - inform%status = NLLS_ERROR_WORKSPACE_ERROR - goto 100 - End If - - w%A(1:n,1:n) = A(1:n,1:n) - x(1:n) = -b(1:n) - call PREC(gesv)( n, 1, w%A, n, w%ipiv, x, n, inform%external_return) - if (inform%external_return .ne. 0 ) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?gesv' - end if - -100 continue - end subroutine minus_solve_general - subroutine matrix_norm(x,A,norm_A_x) Implicit None REAL(wp), intent(in) :: A(:,:), x(:) diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 3fb2e108..68fd1863 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -5,13 +5,14 @@ module MODULE_PREC(ral_nlls_linear) use MODULE_PREC(ral_nlls_workspaces), only: wp, solve_LLS_work, NLLS_inform, & NLLS_options, NLLS_ERROR_WORKSPACE_ERROR, & - NLLS_ERROR_FROM_EXTERNAL + NLLS_ERROR_FROM_EXTERNAL, NLLS_ERROR_BAD_LLS_SOLVER implicit none private public :: solve_LLS, solve_LLS_nocopy contains + ! todo: Jacobian argument subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) ! ----------------------------------------------------------------- ! solve_LLS, a subroutine to solve a linear least squares problem @@ -19,8 +20,8 @@ subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) ! Solves the linear least squares problem Ax = b ! using the method specified in options%lls_solver ! Input: -! A - LHS matrix of the least squares problem (m x n) -! b - The RHS, overwritten with x on output (m x 1) +! A - LHS matrix of the least squares problem +! b - The RHS, overwritten with x on output ! n - Number of columns in A ! m - Number of rows in A ! inform - NLLS_inform structure @@ -28,35 +29,40 @@ subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) ! w - workspace structure ! pd - logical, true if A is known to be positive definite ! Output: -! A_out - output copy of the original A matrix (m x n) -! x - solution vector (n x 1) +! A_out - output copy of the original A matrix +! x - solution vector ! inform - NLLS_inform structure ! work - updated workspace structure ! ----------------------------------------------------------------- implicit none - real(wp), dimension(:), intent(in) :: A - real(wp), dimension(:), intent(in) :: b - real(wp), dimension(:), intent(out) :: A_out - real(wp), dimension(:), intent(out) :: x integer, intent(in) :: n, m + real(wp), intent(in) :: A(n,m) + real(wp), intent(in) :: b(m) + real(wp), intent(out) :: A_out(n,m) + real(wp), intent(out) :: x(n) type(NLLS_inform), intent(inout) :: inform type(NLLS_options), Intent(In) :: options type(solve_LLS_work), intent(inout) :: w logical, intent(in) :: pd + real(wp) :: b_out(m) + A_out = A - x = b + b_out = b - call solve_LLS_nocopy(A_out, x, n, m, inform, w, options, pd) + call solve_LLS_nocopy(A_out, b_out, n, m, inform, w, options, pd) + if (inform%status == 0) then + x(1:n) = b_out(1:n) + end if end subroutine solve_LLS subroutine solve_LLS_nocopy(A, b, n, m, inform, w, options, pd) ! linear solver core which overwrites A and b. implicit none - real(wp), dimension(:), intent(inout) :: A - real(wp), dimension(:), intent(inout) :: b integer, intent(in) :: n, m + real(wp), dimension(m,n), intent(inout) :: A + real(wp), dimension(m), intent(inout) :: b type(NLLS_inform), intent(inout) :: inform type(NLLS_options), intent(in) :: options type(solve_LLS_work), intent(inout) :: w @@ -65,33 +71,34 @@ subroutine solve_LLS_nocopy(A, b, n, m, inform, w, options, pd) select case (options%lls_solver) case (1) ! LAPACK if (pd) then - call solve_posv(A,b,n,inform,options) + call solve_posv(A,b,n,inform) else - if (.not. w%allocated) then - inform%status = NLLS_ERROR_WORKSPACE_ERROR - goto 100 + if (n == m) then + call solve_gesv(A,b,n,inform) + else + if (.not. w%allocated) then + inform%status = NLLS_ERROR_WORKSPACE_ERROR + goto 100 + end if + call solve_gels(A,b,n,m,inform,w,options) end if - call solve_gels(A,b,n,m,inform,w,options) end if + case default + inform%status = NLLS_ERROR_BAD_LLS_SOLVER end select - if (inform%external_return /= 0 ) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - end if - -100 continue +100 continue end subroutine solve_LLS_nocopy subroutine solve_gels(A,b,n,m,inform,w,options) ! Wrapper around LAPACK's ?gels implicit none - REAL(wp), DIMENSION(:), INTENT(IN) :: A - REAL(wp), DIMENSION(:), INTENT(IN) :: b + REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A + REAL(wp), DIMENSION(:), INTENT(INOUT) :: b INTEGER, INTENT(IN) :: n, m type(NLLS_inform), INTENT(INOUT) :: inform type(NLLS_options), Intent(In) :: options - integer, Parameter :: nrhs = 1 integer :: lwork, lda, ldb type( solve_LLS_work ), Intent(inout) :: w @@ -100,12 +107,12 @@ subroutine solve_gels(A,b,n,m,inform,w,options) if (options%Fortran_Jacobian) then lda = m ldb = max(m,n) - call PREC(gels)('N', m, n, nrhs, A, lda, w%temp, ldb, w%work, lwork, & + call PREC(gels)('N', m, n, 1, A, lda, b, ldb, w%work, lwork, & inform%external_return) else lda = n ldb = max(m,n) - call PREC(gels)('T', n, m, nrhs, A, lda, w%temp, ldb, w%work, lwork, & + call PREC(gels)('T', n, m, 1, A, lda, b, ldb, w%work, lwork, & inform%external_return) end if if (inform%external_return /= 0 ) then @@ -115,16 +122,34 @@ subroutine solve_gels(A,b,n,m,inform,w,options) end subroutine solve_gels - subroutine solve_posv(A,b,n,inform,options) + subroutine solve_gesv(A,b,n,inform) +! Wrapper around LAPACK's ?gesv + implicit none + REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A + REAL(wp), DIMENSION(:), INTENT(INOUT) :: b + INTEGER, INTENT(IN) :: n + type(NLLS_inform), INTENT(INOUT) :: inform + + ! NB: we never actually use ipiv in the library + integer, dimension(n) :: ipiv + + call PREC(gesv)(n, 1, A, n, ipiv, b, n, inform%external_return) + if (inform%external_return /= 0 ) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?gesv' + end if + + end subroutine solve_gesv + + subroutine solve_posv(A,b,n,inform) ! Wrapper around LAPACK's ?posv for positive definite systems implicit none - REAL(wp), DIMENSION(:), INTENT(IN) :: A - REAL(wp), DIMENSION(:), INTENT(IN) :: b + REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A + REAL(wp), DIMENSION(:), INTENT(INOUT) :: b INTEGER, INTENT(IN) :: n type(NLLS_inform), INTENT(INOUT) :: inform - type(NLLS_options), Intent(In) :: options - call PREC(posv)('N', n, 1, A, n, b, n, inform%external_return) + call PREC(posv)('L', n, 1, A, n, b, n, inform%external_return) if (inform%external_return /= 0 ) then inform%status = NLLS_ERROR_FROM_EXTERNAL inform%external_name = 'lapack_?posv' diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index d433f79d..f6736bd5 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -48,6 +48,7 @@ module MODULE_PREC(ral_nlls_workspaces) Integer, Parameter, Public :: NLLS_ERROR_BAD_BOX_BOUNDS = -18 Integer, Parameter, Public :: NLLS_ERROR_BAD_JACOBIAN = -19 Integer, Parameter, Public :: NLLS_ERROR_BAD_WEIGHTS = -20 + Integer, Parameter, Public :: NLLS_ERROR_BAD_LLS_SOLVER = -21 ! dogleg errors Integer, Parameter, Public :: NLLS_ERROR_DOGLEG_MODEL = -101 @@ -266,7 +267,7 @@ module MODULE_PREC(ral_nlls_workspaces) ! use eigendecomposition in subproblem solve? LOGICAL :: use_ews_subproblem = .TRUE. - ! This forces to call min_eig_symm without previously calling minus_solve_spd_nocopy + ! This forces to call min_eig_symm without previously calling solve_LLS ! This option is used for code coverage and can be hidden from user. Logical :: force_min_eig_symm = .FALSE. @@ -593,12 +594,6 @@ module MODULE_PREC(ral_nlls_workspaces) real(wp), allocatable :: nullevs(:,:) end type max_eig_work - type, public :: minus_solve_general_work ! workspace for subroutine minus_solve_general - logical :: allocated = .false. - real(wp), allocatable :: A(:,:) - integer, allocatable :: ipiv(:) - end type minus_solve_general_work - type, public :: evaluate_model_work ! workspace for subroutine evaluate_model logical :: allocated = .false. real(wp), allocatable :: Jd(:), dH(:), Hd(:), dHd(:) @@ -606,7 +601,7 @@ module MODULE_PREC(ral_nlls_workspaces) type, public :: solve_LLS_work ! workspace for subroutine solve_LLS logical :: allocated = .false. - real(wp), allocatable :: temp(:), work(:), Jlls(:) + real(wp), allocatable :: temp(:), work(:), Jlls(:,:) end type solve_LLS_work type, public :: min_eig_symm_work ! workspace for subroutine min_eig_work @@ -649,13 +644,15 @@ module MODULE_PREC(ral_nlls_workspaces) type, public :: regularization_solver_work ! workspace for subroutine regularization_solver logical :: allocated = .false. real(wp), allocatable :: AplusSigma(:,:),LtL(:,:) + type( solve_LLS_work ) :: lls_ws end type regularization_solver_work type, public :: more_sorensen_work ! workspace for subroutine more_sorensen logical :: allocated = .false. real(wp), allocatable :: LtL(:,:), AplusSigma(:,:) real(wp), allocatable :: q(:), y1(:) - type( min_eig_symm_work ) :: min_eig_symm_ws + type(min_eig_symm_work) :: min_eig_symm_ws + type(solve_LLS_work) :: lls_ws real(wp), allocatable :: norm_work(:) end type more_sorensen_work @@ -663,8 +660,7 @@ module MODULE_PREC(ral_nlls_workspaces) logical :: allocated = .false. type( max_eig_work ) :: max_eig_ws type( evaluate_model_work ) :: evaluate_model_ws - type( minus_solve_general_work ) :: minus_solve_general_ws - ! type( minus_solve_spd_work ) :: minus_solve_spd_ws + type( solve_LLS_work ) :: solve_LLS_ws REAL(wp), allocatable :: LtL(:,:), B(:,:), p0(:), p1(:) REAL(wp), allocatable :: M0(:,:), M1(:,:), y(:), gtg(:,:), q(:) REAL(wp), allocatable :: M0_small(:,:), M1_small(:,:) @@ -1258,7 +1254,7 @@ subroutine setup_workspace_solve_LLS(n,m,w,options,inform) inform%status = 0 lwork = max(1, min(m,n) + max(min(m,n), 1)*4) - allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(n*m), stat = inform%alloc_status) + allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(n,m), stat = inform%alloc_status) If (inform%alloc_status /= 0) Then Call remove_workspace_solve_LLS(w,options) inform%status = NLLS_ERROR_ALLOCATION @@ -1337,7 +1333,7 @@ subroutine setup_workspace_AINT_tr(n,m,w,options,inform) if (inform%status /= 0) goto 100 ! setup space for the solve routine if ((options%model .ne. 1)) then - call setup_workspace_minus_solve_general(n,m,w%minus_solve_general_ws,options,inform) + call setup_workspace_solve_LLS(n,m,w%solve_LLS_ws,options,inform) if (inform%status /= 0 ) goto 100 end if @@ -1376,7 +1372,7 @@ subroutine remove_workspace_AINT_tr(w,options) call remove_workspace_evaluate_model(w%evaluate_model_ws,options) ! setup space for the solve routine if (options%model .ne. 1) then - call remove_workspace_minus_solve_general(w%minus_solve_general_ws,options) + call remove_workspace_solve_LLS(w%solve_LLS_ws,options) end if w%allocated = .false. @@ -1560,36 +1556,6 @@ subroutine remove_workspace_max_eig(w,options) w%allocated = .false. end subroutine remove_workspace_max_eig - subroutine setup_workspace_minus_solve_general(n, m, w, options, inform) - implicit none - integer, intent(in) :: n, m - type( minus_solve_general_work ), INTENT( INOUT) :: w - type( nlls_options ), intent(in) :: options - type( nlls_inform), intent(inout) :: inform - - inform%status = 0 - allocate( w%A(n,n),w%ipiv(n),stat=inform%alloc_status) - If (inform%alloc_status /= 0) Then - Call remove_workspace_minus_solve_general(w,options) - inform%status = NLLS_ERROR_ALLOCATION - inform%bad_alloc = "setup_workspace_minus_solve_general" - Else - w%allocated = .true. - End If - end subroutine setup_workspace_minus_solve_general - - subroutine remove_workspace_minus_solve_general(w, options) - implicit none - type( minus_solve_general_work ), INTENT( INOUT) :: w - type( nlls_options ), intent(in) :: options - Integer :: ierr_dummy - - if(allocated( w%A )) deallocate( w%A, stat=ierr_dummy ) - if(allocated( w%ipiv )) deallocate( w%ipiv, stat=ierr_dummy ) - - w%allocated = .false. - end subroutine remove_workspace_minus_solve_general - subroutine setup_workspace_solve_galahad(n,m,w,options,inform) implicit none integer, intent(in) :: n,m diff --git a/libRALFit/test/nlls_test.F90 b/libRALFit/test/nlls_test.F90 index 30658c9a..677f15eb 100644 --- a/libRALFit/test/nlls_test.F90 +++ b/libRALFit/test/nlls_test.F90 @@ -1424,7 +1424,13 @@ program nlls_test call all_eig_symm_tests(options,fails) no_errors_helpers = no_errors_helpers + fails - call solve_LLS_tests(options,fails) + call solve_LLS_dgels_tests(options,fails) + no_errors_helpers = no_errors_helpers + fails + + call solve_LLS_dposv_tests(options,fails) + no_errors_helpers = no_errors_helpers + fails + + call solve_LLS_dgesv_tests(options,fails) no_errors_helpers = no_errors_helpers + fails call findbeta_tests(options,fails) @@ -1448,12 +1454,6 @@ program nlls_test call switch_to_quasi_newton_tests(options,fails) no_errors_helpers = no_errors_helpers + fails - call minus_solve_spd_tests(options,fails) - no_errors_helpers = no_errors_helpers + fails - - call minus_solve_general_tests(options,fails) - no_errors_helpers = no_errors_helpers + fails - call matmult_inner_tests(options,fails) no_errors_helpers = no_errors_helpers + fails diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index bcf043ce..b3da03fa 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -1768,7 +1768,7 @@ subroutine more_sorensen_tests(options,fails) call more_sorensen(A,g,n,m,Delta,d,normd,options,status,& work%calculate_step_ws%more_sorensen_ws) if (status%status .ne. 0) then - write(*,*) 'Error: unexpected error in more-sorensen test with non-zero shift' + write(*,*) 'Error: unexpected error in more-sorensen test with non-zero shift and nd /= Delta' write(*,*) 'status = ', status%status, ' returned' fails = fails + 1 status%status = 0 @@ -1782,7 +1782,7 @@ subroutine more_sorensen_tests(options,fails) call more_sorensen(A,g,n,m,Delta,d,normd,options,status,& work%calculate_step_ws%more_sorensen_ws) if (status%status .ne. 0) then - write(*,*) 'Error: unexpected error in more-sorensen test with non-zero shift' + write(*,*) 'Error: unexpected error in more-sorensen test with non-zero shift and nd = Delta' write(*,*) 'status = ', status%status, ' returned' fails = fails + 1 status%status = 0 @@ -1814,6 +1814,7 @@ subroutine more_sorensen_tests(options,fails) work%calculate_step_ws%more_sorensen_ws) if (status%status .ne. NLLS_ERROR_MS_MAXITS) then write(*,*) 'Error: Expected maximum iterations error in more_sorensen' + write(*,*) 'status = ', status%status, ' returned' fails = fails + 1 end if status%status = 0 @@ -1935,6 +1936,7 @@ subroutine trust_region_subproblem_tests(options,fails) if ( (status%status .ne. 0)) then write(*,*) 'Error: unexpected error in ', method_name write(*,*) '(status = ',status%status,')' + write(*,*) '(external return = ',status%external_return,')' write(*,*) 'TR Book Example ',problem_name fails = fails + 1 status%status = 0 @@ -2150,12 +2152,12 @@ subroutine all_eig_symm_tests(options,fails) end subroutine all_eig_symm_tests - subroutine solve_LLS_tests(options,fails) + subroutine solve_LLS_dgels_tests(options,fails) type( nlls_options ), intent(inout) :: options integer, intent(out) :: fails - real(wp), allocatable :: J(:), f(:), d(:), Jd(:) + real(wp), allocatable :: J(:,:), f(:), d(:), Jd(:), JT(:,:) real(wp) :: normerror integer :: n,m type( nlls_workspace ) :: work @@ -2167,21 +2169,19 @@ subroutine solve_LLS_tests(options,fails) iw%iw_ptr => iw options%nlls_method = 1 ! dogleg - n = 2 m = 5 call setup_workspaces(work,n,m,options,status) - allocate(J(n*m), f(m), d(n), Jd(m)) - J = [ 1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp, & - 6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp ] + allocate(J(m,n), JT(m,n), f(m), d(n), Jd(m)) + J = reshape([ 1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp, & + 6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp ], [m, n]) f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & d,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) - d(:) = -d(:) if ( status%status .ne. 0 ) then write (*, *) '[id:1] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2189,13 +2189,14 @@ subroutine solve_LLS_tests(options,fails) end if ! check answer call mult_J(J,n,m,d,Jd,.True.) - normerror = norm2(Jd + f) + normerror = norm2(Jd - f) if ( normerror > 1.0e-12_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS test failed: wrong solution returned' write(*,*) '||Jd - f|| = ', normerror fails = fails + 1 end if + status%status = 0 ! J = [ 1.0_wp, 6.0_wp, ! 2.0_wp, 7.0_wp, @@ -2203,42 +2204,41 @@ subroutine solve_LLS_tests(options,fails) ! 4.0_wp, 9.0_wp, ! 5.0_wp, 10.0_wp ] ! Jacobian Transpose: - J = [ 1.0_wp, 6.0_wp, 2.0_wp, 7.0_wp, 3.0_wp, & - 8.0_wp, 4.0_wp, 9.0_wp, 5.0_wp, 10.0_wp ] + JT = reshape([ 1.0_wp, 6.0_wp, 2.0_wp, 7.0_wp, 3.0_wp, & + 8.0_wp, 4.0_wp, 9.0_wp, 5.0_wp, 10.0_wp ], [n,m]) f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] - - call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + options%fortran_jacobian = .false. + call solve_LLS(jt,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & d,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) - d(:) = -d(:) if ( status%status .ne. 0 ) then - write (*, *) '[id:2] solve_LLS test failed: wrong error message returned' + write (*, *) '[id:2] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" fails = fails + 1 - end if + else ! check answer using J and not JT!!! - J = [ 1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp, & - 6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp ] call mult_J(J,n,m,d,Jd,.True.) - normerror = norm2(Jd + f) + normerror = norm2(Jd - f) if ( normerror > 1.0e-12_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS transpose test failed: wrong solution returned' write(*,*) '||J^Td - f|| = ', normerror fails = fails + 1 end if + end if + options%fortran_jacobian = .true. - n = 65 - m = 60 - deallocate(J,f,Jd,d) - allocate(J(n*m), f(m), d(n)) + n = 65 + m = 60 + deallocate(J,JT,f,Jd,d) + allocate(J(m,n), f(m), d(n)) call setup_workspaces(work,n,m,options,status) - J = 1.0_wp - J(1:m) = 0.0_wp + J(:,:) = 1.0_wp + J(:,1) = 0.0_wp f = 1.0_wp call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & d,n,m,status, & @@ -2251,21 +2251,136 @@ subroutine solve_LLS_tests(options,fails) end if status%status = 0 - call nlls_finalize(work, options, status) + call remove_workspaces(work, options) - call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & + call solve_LLS(j,f,j, & d,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if (status%status .ne. NLLS_ERROR_WORKSPACE_ERROR) then write(*,*) 'Error: workspace error not flagged when workspaces not setup' fails = fails + 1 end if + deallocate(J,f,d) + + + call reset_default_options(options) + + end subroutine solve_LLS_dgels_tests + + subroutine solve_LLS_dposv_tests(options,fails) + + type( nlls_options ), intent(inout) :: options + integer, intent(out) :: fails + + real(wp), allocatable :: A(:,:), b(:), LtL(:,:), x_calc(:), x_true(:), tol + integer :: n + type( nlls_inform ) :: status + type(solve_LLS_work) :: w + + continue + if (wp == np) then + tol = 1.0e-12_wp + else + tol = 5.0e-7_wp + end if + + fails = 0 + + n = 2 + allocate(A(n,n), b(n), LtL(n,n), x_calc(n), x_true(n)) + A = 0.0_wp + A(1,1) = 4.0_wp + A(2,1) = 1.0_wp + A(1,2) = 1.0_wp + A(2,2) = 2.0_wp + x_true(1) = 1.0_wp + x_true(2) = 1.0_wp + b(1) = 5.0_wp + b(2) = 3.0_wp + call solve_LLS(A, b, LtL, x_calc, n, n, status, w, options, .true.) + if (status%status .ne. 0) then + write(*,*) 'Error: solve_LLS info = ', status%status, ' returned from solve_LLS' + fails = fails + 1 + else if (norm2(x_calc-x_true) > tol) then + write(*,*) 'Error: incorrect value returned from solve_LLS for `dposv` case' + write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol + fails = fails + 1 + end if + ! test also _nocopy variant ! + call solve_LLS_nocopy(A, b, n, n, status, w, options, .true.) + if (status%status .ne. 0) then + write(*,*) 'Error: solve_LLS info = ', status%status, ' returned from solve_LLS' + fails = fails + 1 + else if (norm2(x_calc-x_true) > tol) then + write(*,*) 'Error: incorrect value returned from solve_LLS for `dposv` case (nocopy)' + write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol + fails = fails + 1 + end if call reset_default_options(options) - end subroutine solve_LLS_tests + end subroutine solve_LLS_dposv_tests + + subroutine solve_LLS_dgesv_tests(options,fails) + + type( nlls_options ), intent(inout) :: options + integer, intent(out) :: fails + + real(wp), allocatable :: A(:,:), b(:), x_calc(:), x_true(:) + integer :: n,m + type( nlls_inform ) :: status + type( nlls_workspace ) :: work + type( nlls_workspace ), Target :: iw + + fails = 0 + work%iw_ptr => iw + iw%iw_ptr => iw + + n = 2 + m = 3 + + allocate(A(n,n), b(n), x_calc(n), x_true(n)) + + options%nlls_method = 2 + options%model = 2 + call setup_workspaces(work,n,m,options,status) + A = 0.0_wp + A(1,1) = 4.0_wp + A(2,1) = 1.0_wp + A(1,2) = 2.0_wp + A(2,2) = 2.0_wp + x_true(1) = 1.0_wp + x_true(2) = 1.0_wp + x_calc(1) = 6.0_wp + x_calc(2) = 3.0_wp + call solve_LLS_nocopy(A, x_calc, n, n, status, & + work%calculate_step_ws%AINT_tr_ws%solve_LLS_ws, options, .false.) + if (status%status .ne. 0) then + write(*,*) 'Error: info = ', status%status, ' returned from solve_LLS' + fails = fails + 1 + status%status = 0 + else if (norm2(x_true-x_calc) > 1e-12) then + write(*,*) 'Error: incorrect value returned from solve_LLS for `dgesv` case' + write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', 1e-12 + fails = fails + 1 + end if + + A = 0.0_wp + b(1) = 6.0_wp + b(2) = 3.0_wp + call solve_LLS_nocopy(A, b, n, n, status, & + work%calculate_step_ws%AINT_tr_ws%solve_LLS_ws, options, .false.) + if (status%status .ne. NLLS_ERROR_FROM_EXTERNAL) then + write(*,*) 'Error: expected error return from solve_LLS, got info = ', status%status + fails = fails + 1 + end if + + call nlls_finalize(work,options,status) + call reset_default_options(options) + + end subroutine solve_LLS_dgesv_tests subroutine findbeta_tests(options,fails) @@ -2608,126 +2723,6 @@ subroutine switch_to_quasi_newton_tests(options,fails) end subroutine switch_to_quasi_newton_tests - - subroutine minus_solve_spd_tests(options,fails) - - type( nlls_options ), intent(inout) :: options - integer, intent(out) :: fails - - real(wp), allocatable :: A(:,:), b(:), LtL(:,:), x_calc(:), x_true(:), tol - integer :: n - type( nlls_inform ) :: status - - continue - if (wp == np) then - tol = 1.0e-12_wp - else - tol = 5.0e-7_wp - end if - - fails = 0 - - n = 2 - allocate(A(n,n), b(n), LtL(n,n), x_calc(n), x_true(n)) - A = 0.0_wp - A(1,1) = 4.0_wp - A(2,1) = 1.0_wp - A(1,2) = 1.0_wp - A(2,2) = 2.0_wp - x_true(1) = 1.0_wp - x_true(2) = 1.0_wp - b(1) = 5.0_wp - b(2) = 3.0_wp - ! note: b is inverted in minus_solve_spd, so we need to invert it (previously) - b(:) = -b(:) - call minus_solve_spd(A,b,LtL,x_calc,n,status) - if (status%status .ne. 0) then - write(*,*) 'Error: minus_solve_spd info = ', status%status, ' returned from minus_solve_spd' - fails = fails + 1 - else if (norm2(x_calc-x_true) > tol) then - write(*,*) 'Error: incorrect value returned from minus_solve_spd' - write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol - fails = fails + 1 - end if - - ! test also _nocopy variant ! - call minus_solve_spd_nocopy(A,b,x_calc,n,status) - if (status%status .ne. 0) then - write(*,*) 'Error: minus_solve_spd_nocopy info = ', status%status, ' returned from minus_solve_spd' - fails = fails + 1 - else if (norm2(x_calc-x_true) > tol) then - write(*,*) 'Error: incorrect value returned from minus_solve_spd_nocopy' - write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol - fails = fails + 1 - end if - - call reset_default_options(options) - - end subroutine minus_solve_spd_tests - - subroutine minus_solve_general_tests(options,fails) - - type( nlls_options ), intent(inout) :: options - integer, intent(out) :: fails - - real(wp), allocatable :: A(:,:), b(:), x_calc(:), x_true(:) - integer :: n,m - type( nlls_inform ) :: status - type( nlls_workspace ) :: work - type( nlls_workspace ), Target :: iw - - fails = 0 - work%iw_ptr => iw - iw%iw_ptr => iw - - n = 2 - m = 3 - - allocate(A(n,n), b(n), x_calc(n), x_true(n)) - - options%nlls_method = 2 - options%model = 2 - call setup_workspaces(work,n,m,options,status) - - A = 0.0_wp - A(1,1) = 4.0_wp - A(2,1) = 1.0_wp - A(1,2) = 2.0_wp - A(2,2) = 2.0_wp - x_true(1) = 1.0_wp - x_true(2) = 1.0_wp - b(1) = 6.0_wp - b(2) = 3.0_wp - ! note: b is inverted in minus_solve_general, so we need to invert it (previously) - b(:) = -b(:) - call minus_solve_general(A,b,x_calc,n,status,& - work%calculate_step_ws%AINT_tr_ws%minus_solve_general_ws) - if (status%status .ne. 0) then - write(*,*) 'Error: info = ', status%status, ' returned from minus_solve_general' - fails = fails + 1 - status%status = 0 - else if (norm2(x_true-x_calc) > 1e-12) then - write(*,*) 'Error: incorrect value returned from minus_solve_general' - fails = fails + 1 - end if - - A = 0.0_wp - b(1) = 6.0_wp - b(2) = 3.0_wp - ! note: b is inverted in minus_solve_general, so we need to invert it (previously) - b(:) = -b(:) - call minus_solve_general(A,b,x_calc,n,status,& - work%calculate_step_ws%AINT_tr_ws%minus_solve_general_ws) - if (status%status .ne. NLLS_ERROR_FROM_EXTERNAL) then - write(*,*) 'Error: expected error return from minus_solve_general, got info = ', status%status - fails = fails + 1 - end if - - call nlls_finalize(work,options,status) - call reset_default_options(options) - - end subroutine minus_solve_general_tests - subroutine matmult_inner_tests(options,fails) type( nlls_options ), intent(inout) :: options From c66658a6b23c09fbecd74fdd241b59250b7682e1 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Tue, 10 Feb 2026 07:43:23 +0000 Subject: [PATCH 04/16] some review fixes --- libRALFit/doc/common/errors.rst | 2 ++ libRALFit/log | 0 libRALFit/src/ral_nlls_internal.F90 | 9 +++---- libRALFit/src/ral_nlls_linear.F90 | 4 +++ libRALFit/test/unit_test_mod.F90 | 38 +++++++++++++++++++++++------ 5 files changed, 41 insertions(+), 12 deletions(-) delete mode 100644 libRALFit/log diff --git a/libRALFit/doc/common/errors.rst b/libRALFit/doc/common/errors.rst index b0b35b16..079da402 100644 --- a/libRALFit/doc/common/errors.rst +++ b/libRALFit/doc/common/errors.rst @@ -40,6 +40,8 @@ Possible values are: - One or more elements in the Jacobian appear to be wrong * - -20 - Weights vector must be sufficiently positive + * - -21 + - Unsupported LLS solver (``lls_solver``) specified in options. * - -101 - Unsupported model in dogleg (``nlls_method=1``). * - -201 diff --git a/libRALFit/log b/libRALFit/log deleted file mode 100644 index e69de29b..00000000 diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index ac4422fe..52573709 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -1460,6 +1460,7 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ options, inform) end if subproblem_method = 2 ! try aint next + call dogleg(J,f,hf,g,n,m,Delta,d,norm_S_d, & options,inform,w%dogleg_ws) if (inform%status == 0) subproblem_success = .true. @@ -1976,9 +1977,9 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) ! Solve (A+lam*B)x=-v: ! 1. we can use w%B to actually store A + lam * B w%B(:,:) = A(:,:) - w%p1(:) = -v(:) Do i = 1, n w%B(i,i) = w%B(i,i) + lam + w%p1(i) = -v(i) End Do call solve_LLS_nocopy(w%B, w%p1, n, n, inform, w%solve_LLS_ws, options, pd) if (inform%status /= 0) goto 100 @@ -2212,7 +2213,7 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) Write(rec(1), Fmt=5040) Call Printmsg(5,.False.,options,1,rec) End If - goto 100 + goto 100 end if w%q(:) = d ! w%q = R'\d @@ -2557,9 +2558,8 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) elseif (region == 3) then call shift_matrix(A,sigma + sigma_shift,w%AplusSigma,n) ! solve (A + sigma)d = -v - w%LtL(1:n,1:n) = w%AplusSigma(1:n,1:n) d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) if (inform%status == 0) then region = 2 @@ -3446,7 +3446,6 @@ subroutine mult_Jt(J,n,m,x,Jtx,Fortran_Jacobian) End If end subroutine mult_Jt - subroutine scale_J_by_weights(J,n,m,weights,options) Implicit None real(wp), intent(inout) :: J(*) diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 68fd1863..45e67ea3 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -1,3 +1,7 @@ +! Copyright (c) 2020, The Science and Technology Facilities Council (STFC) +! All rights reserved. +! Copyright (C) 2020 Numerical Algorithms Group (NAG). All rights reserved. +! Copyright (C) 2024 Advanced Micro Devices, Inc. All rights reserved. ! ral_nlls_linear :: linear solvers for internal use in RALFit #include "preprocessor.FPP" diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index b3da03fa..094f0cfa 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -2174,10 +2174,22 @@ subroutine solve_LLS_dgels_tests(options,fails) call setup_workspaces(work,n,m,options,status) - allocate(J(m,n), JT(m,n), f(m), d(n), Jd(m)) - J = reshape([ 1.0_wp, 2.0_wp, 3.0_wp, 4.0_wp, 5.0_wp, & - 6.0_wp, 7.0_wp, 8.0_wp, 9.0_wp, 10.0_wp ], [m, n]) - f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] + allocate(J(m,n), JT(n,m), f(m), d(n), Jd(m)) + J(1, 1) = 1.0_wp + J(2, 1) = 2.0_wp + J(3, 1) = 3.0_wp + J(4, 1) = 4.0_wp + J(5, 1) = 5.0_wp + J(1, 2) = 6.0_wp + J(2, 2) = 7.0_wp + J(3, 2) = 8.0_wp + J(4, 2) = 9.0_wp + J(5, 2) = 10.0_wp + f(1) = 7.0_wp + f(2) = 9.0_wp + f(3) = 11.0_wp + f(4) = 13.0_wp + f(5) = 15.0_wp call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & d,n,m,status, & @@ -2204,9 +2216,21 @@ subroutine solve_LLS_dgels_tests(options,fails) ! 4.0_wp, 9.0_wp, ! 5.0_wp, 10.0_wp ] ! Jacobian Transpose: - JT = reshape([ 1.0_wp, 6.0_wp, 2.0_wp, 7.0_wp, 3.0_wp, & - 8.0_wp, 4.0_wp, 9.0_wp, 5.0_wp, 10.0_wp ], [n,m]) - f = [ 7.0_wp, 9.0_wp, 11.0_wp, 13.0_wp, 15.0_wp ] + JT(1, 1) = 1.0_wp + JT(2, 1) = 6.0_wp + JT(1, 2) = 2.0_wp + JT(2, 2) = 7.0_wp + JT(1, 3) = 3.0_wp + JT(2, 3) = 8.0_wp + JT(1, 4) = 4.0_wp + JT(2, 4) = 9.0_wp + JT(1, 5) = 5.0_wp + JT(2, 5) = 10.0_wp + f(1) = 7.0_wp + f(2) = 9.0_wp + f(3) = 11.0_wp + f(4) = 13.0_wp + f(5) = 15.0_wp options%fortran_jacobian = .false. call solve_LLS(jt,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & From d397f8e12508e83515ee31f9b4f4daeed9a05cfb Mon Sep 17 00:00:00 2001 From: alexhroom Date: Fri, 13 Feb 2026 08:58:31 +0000 Subject: [PATCH 05/16] removed solve_LLS with copy --- libRALFit/src/ral_nlls_internal.F90 | 22 ++++++------ libRALFit/src/ral_nlls_linear.F90 | 53 ++++++----------------------- libRALFit/test/unit_test_mod.F90 | 43 ++++++++++------------- 3 files changed, 39 insertions(+), 79 deletions(-) diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index 52573709..a31a0ed6 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -14,7 +14,7 @@ module MODULE_PREC(ral_nlls_internal) Use MODULE_PREC(ral_nlls_printing) Use MODULE_PREC(ral_nlls_bounds) Use MODULE_PREC(ral_nlls_types), Only: PREC(nrm2), PREC(dot) - use MODULE_PREC(ral_nlls_linear), only: solve_LLS, solve_LLS_nocopy + use MODULE_PREC(ral_nlls_linear), only: solve_LLS implicit none @@ -1805,7 +1805,7 @@ SUBROUTINE dogleg(J,f,hf,g,n,m,Delta,d,normd,options,inform,w) ! linear model... w%solve_LLS_ws%Jlls = reshape(J, [n,m]) w%solve_LLS_ws%temp(1:m) = f(1:m) - call solve_LLS_nocopy(w%solve_LLS_ws%Jlls,w%solve_LLS_ws%temp,n,m,inform,w%solve_LLS_ws,options,.false.) + call solve_LLS(w%solve_LLS_ws%Jlls,w%solve_LLS_ws%temp,n,m,inform,w%solve_LLS_ws,options,.false.) if ( inform%status /= 0 ) goto 100 w%d_gn = -w%solve_LLS_ws%temp(1:n) case default @@ -1900,7 +1900,7 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) w%LtL(:, :) = A(:, :) w%p0(:) = -v(:) - call solve_LLS_nocopy(w%LtL,w%p0,n,n,inform,w%solve_LLS_ws,options,pd) + call solve_LLS(w%LtL,w%p0,n,n,inform,w%solve_LLS_ws,options,pd) if (inform%status /= 0) goto 100 @@ -1960,7 +1960,7 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) w%M0_small(:,:) = A(:,:) + lam*w%B(:,:) + w%M1_small ! solve Hq + g = 0 for q w%q(:) = -v(:) - call solve_LLS_nocopy(w%M0_small, w%q, n, n, inform, w%solve_LLS_ws, options, pd) + call solve_LLS(w%M0_small, w%q, n, n, inform, w%solve_LLS_ws, options, pd) if (inform%status /= 0) goto 100 ! find max eta st ||q + eta v(:,1)||_B = Delta @@ -1981,7 +1981,7 @@ SUBROUTINE AINT_tr(J,A,f,X,v,hf,n,m,Delta,d,normd,options,inform,w) w%B(i,i) = w%B(i,i) + lam w%p1(i) = -v(i) End Do - call solve_LLS_nocopy(w%B, w%p1, n, n, inform, w%solve_LLS_ws, options, pd) + call solve_LLS(w%B, w%p1, n, n, inform, w%solve_LLS_ws, options, pd) if (inform%status /= 0) goto 100 end if @@ -2099,7 +2099,7 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) Else ! solve (A+sigma)d = -v d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) End If if (inform%status == 0) then ! A is symmetric positive definite.... @@ -2247,7 +2247,7 @@ subroutine more_sorensen_ew(A,v,n,m,Delta,d,nd,options,inform,w) call shift_matrix(A,sigma,w%AplusSigma,n) ! solve (A + sigma)d = -v d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) end if if (inform%status /= 0) goto 100 @@ -2381,7 +2381,7 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) ! solve (A + sigma)d = -v w%LtL(1:n,1:n) = w%AplusSigma(1:n,1:n) d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) else factorization_done = .false. end if @@ -2559,7 +2559,7 @@ subroutine more_sorensen_noew(A,v,n,m,Delta,d,nd,options,inform,w) call shift_matrix(A,sigma + sigma_shift,w%AplusSigma,n) ! solve (A + sigma)d = -v d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) if (inform%status == 0) then region = 2 @@ -2734,7 +2734,7 @@ subroutine check_shift_and_solve(n,A,AplusSigma,v,sigma,d,options,inform,w) call shift_matrix(A,sigma,w%AplusSigma,n) ! solve (A + sigma)d = -v d(:) = -v(:) - call solve_LLS_nocopy(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) if ( inform%status == 0 ) then successful_shift = .true. else @@ -2922,7 +2922,7 @@ subroutine regularization_solver(A,v,n,m,Delta,num_successful_steps,& ! solve (A+sigma)d = -v w%LtL(1:n, 1:n) = w%AplusSigma(1:n,1:n) d(1:n) = -v(1:n) - call solve_LLS_nocopy(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) + call solve_LLS(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) ! informa%status is passed along, this routine exits here else ! Feature not yet implemented, this should have been caught in diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 45e67ea3..d5a52635 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -13,11 +13,11 @@ module MODULE_PREC(ral_nlls_linear) implicit none private - public :: solve_LLS, solve_LLS_nocopy + public :: solve_LLS contains ! todo: Jacobian argument - subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) + subroutine solve_LLS(A, b, n, m, inform, w, options, pd) ! ----------------------------------------------------------------- ! solve_LLS, a subroutine to solve a linear least squares problem ! ----------------------------------------------------------------- @@ -25,51 +25,20 @@ subroutine solve_LLS(A,b,A_out,x,n,m,inform,w,options,pd) ! using the method specified in options%lls_solver ! Input: ! A - LHS matrix of the least squares problem -! b - The RHS, overwritten with x on output +! b - The RHS, overwritten with result x on output ! n - Number of columns in A ! m - Number of rows in A ! inform - NLLS_inform structure -! options - NLLS_options structure ! w - workspace structure +! options - NLLS_options structure ! pd - logical, true if A is known to be positive definite -! Output: -! A_out - output copy of the original A matrix -! x - solution vector -! inform - NLLS_inform structure -! work - updated workspace structure ! ----------------------------------------------------------------- implicit none + real(wp), intent(inout), contiguous :: A(:,:), b(:) integer, intent(in) :: n, m - real(wp), intent(in) :: A(n,m) - real(wp), intent(in) :: b(m) - real(wp), intent(out) :: A_out(n,m) - real(wp), intent(out) :: x(n) type(NLLS_inform), intent(inout) :: inform - type(NLLS_options), Intent(In) :: options type(solve_LLS_work), intent(inout) :: w - logical, intent(in) :: pd - - real(wp) :: b_out(m) - - A_out = A - b_out = b - - call solve_LLS_nocopy(A_out, b_out, n, m, inform, w, options, pd) - if (inform%status == 0) then - x(1:n) = b_out(1:n) - end if - - end subroutine solve_LLS - - subroutine solve_LLS_nocopy(A, b, n, m, inform, w, options, pd) -! linear solver core which overwrites A and b. - implicit none - integer, intent(in) :: n, m - real(wp), dimension(m,n), intent(inout) :: A - real(wp), dimension(m), intent(inout) :: b - type(NLLS_inform), intent(inout) :: inform type(NLLS_options), intent(in) :: options - type(solve_LLS_work), intent(inout) :: w logical, intent(in) :: pd select case (options%lls_solver) @@ -92,13 +61,12 @@ subroutine solve_LLS_nocopy(A, b, n, m, inform, w, options, pd) end select 100 continue - end subroutine solve_LLS_nocopy + end subroutine solve_LLS subroutine solve_gels(A,b,n,m,inform,w,options) ! Wrapper around LAPACK's ?gels implicit none - REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A - REAL(wp), DIMENSION(:), INTENT(INOUT) :: b + real(wp), intent(inout), contiguous :: A(:,:), b(:) INTEGER, INTENT(IN) :: n, m type(NLLS_inform), INTENT(INOUT) :: inform type(NLLS_options), Intent(In) :: options @@ -129,8 +97,7 @@ end subroutine solve_gels subroutine solve_gesv(A,b,n,inform) ! Wrapper around LAPACK's ?gesv implicit none - REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A - REAL(wp), DIMENSION(:), INTENT(INOUT) :: b + real(wp), intent(inout), contiguous :: A(:,:), b(:) INTEGER, INTENT(IN) :: n type(NLLS_inform), INTENT(INOUT) :: inform @@ -148,8 +115,8 @@ end subroutine solve_gesv subroutine solve_posv(A,b,n,inform) ! Wrapper around LAPACK's ?posv for positive definite systems implicit none - REAL(wp), DIMENSION(:,:), INTENT(INOUT) :: A - REAL(wp), DIMENSION(:), INTENT(INOUT) :: b + REAL(wp), DIMENSION(:,:), INTENT(INOUT), contiguous :: A + REAL(wp), DIMENSION(:), INTENT(INOUT), contiguous :: b INTEGER, INTENT(IN) :: n type(NLLS_inform), INTENT(INOUT) :: inform diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index 094f0cfa..d8805bcc 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -2157,7 +2157,7 @@ subroutine solve_LLS_dgels_tests(options,fails) type( nlls_options ), intent(inout) :: options integer, intent(out) :: fails - real(wp), allocatable :: J(:,:), f(:), d(:), Jd(:), JT(:,:) + real(wp), allocatable :: J(:,:), J_copy(:,:), f(:), d(:), Jd(:), JT(:,:), JT_copy(:,:) real(wp) :: normerror integer :: n,m type( nlls_workspace ) :: work @@ -2174,7 +2174,9 @@ subroutine solve_LLS_dgels_tests(options,fails) call setup_workspaces(work,n,m,options,status) - allocate(J(m,n), JT(n,m), f(m), d(n), Jd(m)) + allocate(J(m,n), J_copy(m,n), & + JT(n,m), JT_copy(n,m), & + f(m), d(n), Jd(m)) J(1, 1) = 1.0_wp J(2, 1) = 2.0_wp J(3, 1) = 3.0_wp @@ -2191,8 +2193,10 @@ subroutine solve_LLS_dgels_tests(options,fails) f(4) = 13.0_wp f(5) = 15.0_wp - call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & - d,n,m,status, & + J_copy = J + d = f + + call solve_LLS(J_copy,d,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if ( status%status .ne. 0 ) then write (*, *) '[id:1] solve_LLS test failed: wrong error message returned' @@ -2232,9 +2236,11 @@ subroutine solve_LLS_dgels_tests(options,fails) f(4) = 13.0_wp f(5) = 15.0_wp + JT_copy = JT + d = f + options%fortran_jacobian = .false. - call solve_LLS(jt,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & - d,n,m,status, & + call solve_LLS(JT_copy,d,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if ( status%status .ne. 0 ) then @@ -2264,8 +2270,7 @@ subroutine solve_LLS_dgels_tests(options,fails) J(:,:) = 1.0_wp J(:,1) = 0.0_wp f = 1.0_wp - call solve_LLS(j,f,work%calculate_step_ws%dogleg_ws%solve_lls_ws%jlls, & - d,n,m,status, & + call solve_LLS(j,f,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if ( status%status .ne. NLLS_ERROR_FROM_EXTERNAL ) then @@ -2277,8 +2282,7 @@ subroutine solve_LLS_dgels_tests(options,fails) call remove_workspaces(work, options) - call solve_LLS(j,f,j, & - d,n,m,status, & + call solve_LLS(j,f,n,m,status, & work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) if (status%status .ne. NLLS_ERROR_WORKSPACE_ERROR) then write(*,*) 'Error: workspace error not flagged when workspaces not setup' @@ -2321,27 +2325,16 @@ subroutine solve_LLS_dposv_tests(options,fails) x_true(2) = 1.0_wp b(1) = 5.0_wp b(2) = 3.0_wp - call solve_LLS(A, b, LtL, x_calc, n, n, status, w, options, .true.) + call solve_LLS(A, b, n, n, status, w, options, .true.) if (status%status .ne. 0) then write(*,*) 'Error: solve_LLS info = ', status%status, ' returned from solve_LLS' fails = fails + 1 - else if (norm2(x_calc-x_true) > tol) then + else if (norm2(b-x_true) > tol) then write(*,*) 'Error: incorrect value returned from solve_LLS for `dposv` case' write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol fails = fails + 1 end if - ! test also _nocopy variant ! - call solve_LLS_nocopy(A, b, n, n, status, w, options, .true.) - if (status%status .ne. 0) then - write(*,*) 'Error: solve_LLS info = ', status%status, ' returned from solve_LLS' - fails = fails + 1 - else if (norm2(x_calc-x_true) > tol) then - write(*,*) 'Error: incorrect value returned from solve_LLS for `dposv` case (nocopy)' - write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol - fails = fails + 1 - end if - call reset_default_options(options) end subroutine solve_LLS_dposv_tests @@ -2379,7 +2372,7 @@ subroutine solve_LLS_dgesv_tests(options,fails) x_true(2) = 1.0_wp x_calc(1) = 6.0_wp x_calc(2) = 3.0_wp - call solve_LLS_nocopy(A, x_calc, n, n, status, & + call solve_LLS(A, x_calc, n, n, status, & work%calculate_step_ws%AINT_tr_ws%solve_LLS_ws, options, .false.) if (status%status .ne. 0) then write(*,*) 'Error: info = ', status%status, ' returned from solve_LLS' @@ -2394,7 +2387,7 @@ subroutine solve_LLS_dgesv_tests(options,fails) A = 0.0_wp b(1) = 6.0_wp b(2) = 3.0_wp - call solve_LLS_nocopy(A, b, n, n, status, & + call solve_LLS(A, b, n, n, status, & work%calculate_step_ws%AINT_tr_ws%solve_LLS_ws, options, .false.) if (status%status .ne. NLLS_ERROR_FROM_EXTERNAL) then write(*,*) 'Error: expected error return from solve_LLS, got info = ', status%status From e34cebc7b9d4ff64cbf9e877769aeb456c40240f Mon Sep 17 00:00:00 2001 From: alexhroom Date: Fri, 13 Feb 2026 09:27:23 +0000 Subject: [PATCH 06/16] added test for BAD_LLS_SOLVER error --- libRALFit/test/nlls_test.F90 | 3 +++ libRALFit/test/unit_test_mod.F90 | 28 ++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+) diff --git a/libRALFit/test/nlls_test.F90 b/libRALFit/test/nlls_test.F90 index 677f15eb..e6c0658d 100644 --- a/libRALFit/test/nlls_test.F90 +++ b/libRALFit/test/nlls_test.F90 @@ -1424,6 +1424,9 @@ program nlls_test call all_eig_symm_tests(options,fails) no_errors_helpers = no_errors_helpers + fails + call solve_LLS_general_tests(options,fails) + no_errors_helpers = no_errors_helpers + fails + call solve_LLS_dgels_tests(options,fails) no_errors_helpers = no_errors_helpers + fails diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index d8805bcc..3f41567c 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -2152,6 +2152,31 @@ subroutine all_eig_symm_tests(options,fails) end subroutine all_eig_symm_tests + subroutine solve_LLS_general_tests(options,fails) + ! Tests for solve_LLS that do not depend on any actual solvers, just the logic of the routine itself + type( nlls_options ), intent(inout) :: options + integer, intent(out) :: fails + + ! test that a BAD_LLS_SOLVER error is returned if an invalid solver is passed in + real(wp), allocatable :: J(:,:), d(:) + integer :: n,m + type( solve_LLS_work ) :: work + type( nlls_inform ) :: inform + + fails = 0 + options%lls_solver = -1 ! invalid solver + + call solve_LLS(J,d,n,m,inform,work,options,.false.) + + if (inform%status .ne. NLLS_ERROR_BAD_LLS_SOLVER) then + write(*,*) 'Error: expected BAD_LLS_SOLVER error when passing uninitialized work to solve_LLS' + write(*,*) 'status = ', inform%status, ' returned' + fails = fails + 1 + end if + + end subroutine solve_LLS_general_tests + + subroutine solve_LLS_dgels_tests(options,fails) type( nlls_options ), intent(inout) :: options @@ -2168,6 +2193,7 @@ subroutine solve_LLS_dgels_tests(options,fails) work%iw_ptr => iw iw%iw_ptr => iw + options%lls_solver = 1 ! LAPACK options%nlls_method = 1 ! dogleg n = 2 m = 5 @@ -2313,6 +2339,7 @@ subroutine solve_LLS_dposv_tests(options,fails) end if fails = 0 + options%lls_solver = 1 ! LAPACK n = 2 allocate(A(n,n), b(n), LtL(n,n), x_calc(n), x_true(n)) @@ -2359,6 +2386,7 @@ subroutine solve_LLS_dgesv_tests(options,fails) allocate(A(n,n), b(n), x_calc(n), x_true(n)) + options%lls_solver = 1 ! LAPACK options%nlls_method = 2 options%model = 2 call setup_workspaces(work,n,m,options,status) From 4e782a75e23102f2d20017f9069371d0ac064726 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 2 Mar 2026 09:48:49 +0000 Subject: [PATCH 07/16] changed inequality logic --- libRALFit/test/unit_test_mod.F90 | 62 ++++++++++++++++---------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index 3f41567c..f7414e8b 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -1412,7 +1412,7 @@ subroutine c_fortran_tests(options, fails, tol_type) if ( c_status%iter == status%iter ) then resvec_error = norm2(c_status%resvec(1:c_status%iter+1) - & status%resvec(1:status%iter+1)) - if (resvec_error > abstol) then + if (.not. resvec_error < abstol) then write(*,*) 'Warning: fortran and c resvec differ!' write(*,*) 'Different resvecs in 2-norm for' write(*,*) 'NLLS_METHOD = ', options%nlls_method @@ -1940,7 +1940,7 @@ subroutine trust_region_subproblem_tests(options,fails) write(*,*) 'TR Book Example ',problem_name fails = fails + 1 status%status = 0 - elseif ( normd - Delta > 1e-3 ) then + elseif ( .not. normd - Delta < 1e-3 ) then write(*,*) 'Error: answer returned outside the TR Radius' write(*,*) 'TR Book Example ', trim(problem_name),' using method ',method_name write(*,*) 'Delta = ', Delta, '||d|| = ', normd @@ -2051,7 +2051,7 @@ subroutine solve_galahad_tests(options,fails) if (i == 1) then ! check result lies within the trust region - if ( abs(dot_product(d,d) - Delta**2) > tol ) then + if ( .not. abs(dot_product(d,d) - Delta**2) < tol ) then write(*,*) testname,'failed' write(*,*) 'Delta = ', Delta, '||d|| = ', sqrt(dot_product(d,d)) write(*,*) 'sq diff = ', abs(dot_product(d,d) - Delta**2), 'tol = ', tol @@ -2232,7 +2232,7 @@ subroutine solve_LLS_dgels_tests(options,fails) ! check answer call mult_J(J,n,m,d,Jd,.True.) normerror = norm2(Jd - f) - if ( normerror > 1.0e-12_wp ) then + if ( .not. normerror < 1.0e-12_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS test failed: wrong solution returned' write(*,*) '||Jd - f|| = ', normerror @@ -2277,7 +2277,7 @@ subroutine solve_LLS_dgels_tests(options,fails) ! check answer using J and not JT!!! call mult_J(J,n,m,d,Jd,.True.) normerror = norm2(Jd - f) - if ( normerror > 1.0e-12_wp ) then + if ( .not. normerror < 1.0e-12_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS transpose test failed: wrong solution returned' write(*,*) '||J^Td - f|| = ', normerror @@ -2356,7 +2356,7 @@ subroutine solve_LLS_dposv_tests(options,fails) if (status%status .ne. 0) then write(*,*) 'Error: solve_LLS info = ', status%status, ' returned from solve_LLS' fails = fails + 1 - else if (norm2(b-x_true) > tol) then + else if (.not. norm2(b-x_true) < tol) then write(*,*) 'Error: incorrect value returned from solve_LLS for `dposv` case' write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', tol fails = fails + 1 @@ -2406,7 +2406,7 @@ subroutine solve_LLS_dgesv_tests(options,fails) write(*,*) 'Error: info = ', status%status, ' returned from solve_LLS' fails = fails + 1 status%status = 0 - else if (norm2(x_true-x_calc) > 1e-12) then + else if (.not. norm2(x_true-x_calc) < 1e-12) then write(*,*) 'Error: incorrect value returned from solve_LLS for `dgesv` case' write(*,*) 'diff norm 2 = ', norm2(x_calc-x_true), 'tol = ', 1e-12 fails = fails + 1 @@ -2449,7 +2449,7 @@ subroutine findbeta_tests(options,fails) if (status%status .ne. 0) then write(*,*) 'error -- findbeta did not work: info /= 0' fails = fails + 1 - else if ( ( norm2( a + beta * b ) - 10.0_wp ) > 1e-12 ) then + else if ( .not. ( norm2( a + beta * b ) - 10.0_wp ) < 1e-12 ) then write(*,*) 'error -- findbeta did not work' write(*,*) '|| x + beta y|| = ', norm2( (a + beta * b)-10.0_wp) fails = fails + 1 @@ -2493,7 +2493,7 @@ subroutine calculate_rho_tests(options,fails) normfnew = 1.0_wp md = 1.5_wp call calculate_rho(normf, normfnew, md, rho,options) - if ( abs(rho - 3.0_wp) > 1e-10) then + if ( .not. abs(rho - 3.0_wp) < 1e-10) then write(*,*) 'Unexpected answer from calculate_rho' write(*,*) 'Expected 3.0, got ', rho fails = fails + 1 @@ -2502,7 +2502,7 @@ subroutine calculate_rho_tests(options,fails) ! now, let's check one is returned if alpha = beta normfnew = 2.0_wp call calculate_rho(normf, normfnew, md, rho,options) - if (abs(rho - 1.0_wp) > 1e-10) then + if (.not. abs(rho - 1.0_wp) < 1e-10) then write(*,*) 'Unexpected answer from calculate_rho' write(*,*) 'Expected 1.0, got ', rho fails = fails + 1 @@ -2512,7 +2512,7 @@ subroutine calculate_rho_tests(options,fails) ! finally, check that 1 is returned if denominator = 0 md = 2.0_wp call calculate_rho(normf, normfnew, md, rho,options) - if (abs(rho - 1.0_wp) > 1e-10) then + if (.not. abs(rho - 1.0_wp) < 1e-10) then write(*,*) 'Unexpected answer from calculate_rho' write(*,*) 'Expected 1.0, got ', rho fails = fails + 1 @@ -2552,7 +2552,7 @@ subroutine update_trust_region_radius_tests(options,fails) ! check if rho reduced... rho = options%eta_success_but_reduce - 0.5_wp call update_trust_region_radius(rho,options,status,work) - if ( work%Delta >= 100.0_wp ) then + if ( .not. work%Delta <= 100.0_wp ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not decrease as expected: delta = ', work%Delta fails = fails + 1 @@ -2562,7 +2562,7 @@ subroutine update_trust_region_radius_tests(options,fails) ! check if rho stays the same... rho = (options%eta_success_but_reduce + options%eta_very_successful) / 2 call update_trust_region_radius(rho,options,status,work) - if ( abs(work%Delta - 100.0_wp) > 1e-12 ) then + if ( .not. abs(work%Delta - 100.0_wp) < 1e-12 ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not stay the same: Delta = ', work%Delta fails = fails + 1 @@ -2573,7 +2573,7 @@ subroutine update_trust_region_radius_tests(options,fails) rho = (options%eta_very_successful + options%eta_too_successful) / 2 work%norm_S_d = 100.0_wp call update_trust_region_radius(rho,options,status,work) - if ( work%Delta <= 100.0_wp ) then + if ( .not. work%Delta >= 100.0_wp ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not incease: delta = ', work%Delta fails = fails + 1 @@ -2584,7 +2584,7 @@ subroutine update_trust_region_radius_tests(options,fails) ! check if rho stays the same because too successful... rho = options%eta_too_successful + 1.0_wp call update_trust_region_radius(rho,options,status,work) - if ( abs(work%Delta - 100.0_wp) > 1e-12 ) then + if ( .not. abs(work%Delta - 100.0_wp) < 1e-12 ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not stay the same: delta = ', work%Delta fails = fails + 1 @@ -2609,7 +2609,7 @@ subroutine update_trust_region_radius_tests(options,fails) ! check if rho stays the same because too successful... rho = options%eta_too_successful + 1.0_wp call update_trust_region_radius(rho,options,status,work) - if ( abs(work%Delta - 100.0_wp) > 1e-12 ) then + if ( .not. abs(work%Delta - 100.0_wp) < 1e-12 ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not stay the same: delta = ', work%Delta fails = fails + 1 @@ -2618,7 +2618,7 @@ subroutine update_trust_region_radius_tests(options,fails) rho = options%eta_success_but_reduce - 0.5_wp call update_trust_region_radius(rho,options,status,work) - if ( work%Delta >= 100.0_wp ) then + if ( .not. work%Delta <= 100.0_wp ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not decrease as expected: delta = ', work%Delta fails = fails + 1 @@ -2627,7 +2627,7 @@ subroutine update_trust_region_radius_tests(options,fails) rho = options%eta_successful - 10.0_wp call update_trust_region_radius(rho,options,status,work) - if ( work%Delta >= 100.0_wp ) then + if ( .not. work%Delta <= 100.0_wp ) then write(*,*) 'Unexpected answer from update_trust_region_radius' write(*,*) 'Delta did not decrease as expected: delta = ', work%Delta fails = fails + 1 @@ -2698,7 +2698,7 @@ subroutine mult_J_tests(options,fails) x = 1.0_wp J = [ 1.0 , 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 ] call mult_J(J,m,n,x,Jx,.True.) - if ( norm2( Jx - [16.0, 20.0 ] ) > 1e-12) then + if ( .not. norm2( Jx - [16.0, 20.0 ] ) < 1e-12) then write(*,*) 'error :: mult_J test failed' fails = fails + 1 end if @@ -2726,7 +2726,7 @@ subroutine mult_Jt_tests(options,fails) x = 1.0_wp J = [ 1.0 , 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0 ] call mult_Jt(J,n,m,x,Jtx,.True.) - if ( norm2( Jtx - [10.0, 26.0 ] ) > 1e-12) then + if ( .not. norm2( Jtx - [10.0, 26.0 ] ) < 1e-12) then write(*,*) 'error :: mult_Jt test failed' fails = fails + 1 end if @@ -2794,7 +2794,7 @@ subroutine matmult_inner_tests(options,fails) do i = 1,n diff(i) = norm2(AtA(:,i) - AtA_expected(:,i)) end do - if (norm2(diff) > 1e-10) then + if (.not. norm2(diff) < 1e-10) then write(*,*) 'error :: matmult_inner test failed' fails = fails + 1 end if @@ -2833,7 +2833,7 @@ subroutine matmult_outer_tests(options,fails) do i = 1,m diff(i) = norm2(AAt(:,i) - AAt_expected(:,i)) end do - if (norm2(diff) > 1e-10) then + if (.not. norm2(diff) < 1e-10) then write(*,*) 'error :: matmult_outer test failed' fails = fails + 1 end if @@ -2871,7 +2871,7 @@ subroutine outer_product_tests(options,fails) do i = 1, n diff(i) = norm2(xxt(i,:) - xxt_exact(i,:)) end do - if (norm2(diff) > 1e-12) then + if (.not. norm2(diff) < 1e-12) then write(*,*) 'error :: outer_product test failed' fails = fails +1 end if @@ -2943,7 +2943,7 @@ subroutine min_eig_symm_tests(options,fails) call min_eig_symm(A,n,ew,ev,options,status, & work%calculate_step_ws%more_sorensen_ws%min_eig_symm_ws) - if ( (abs( ew + 6.0_wp ) > tol ).or.(status%status .ne. 0) ) then + if ( .not. (abs( ew + 6.0_wp ) < tol ).or.(status%status .ne. 0) ) then write(*,*) 'error :: min_eig_symm test failed -- wrong eig found' write(*,*) 'diff = ', abs( ew + 6.0_wp ), 'tol = ', tol fails = fails +1 @@ -2958,7 +2958,7 @@ subroutine min_eig_symm_tests(options,fails) errnorm = errnorm + (avec(row) - ew * ev(row))**2 end do errnorm = sqrt(errnorm) - if (errnorm > tol) then + if (.not. errnorm < tol) then write(*,*) 'error :: min_eig_symm test failed -- not an eigenvector' write(*,*) 'diff = ', errnorm, 'tol = ', tol fails = fails +1 @@ -3031,7 +3031,7 @@ subroutine max_eig_tests(options,fails) if ( status%status .ne. 0 ) then write(*,*) 'error :: max_eig test failed, status = ', status%status fails = fails + 1 - elseif ( (abs( ew - 10.0_wp) > tol) ) then + elseif ( .not. (abs( ew - 10.0_wp) < tol) ) then write(*,*) 'error :: max_eig test failed, incorrect answer' write(*,*) 'expected 10.0, got ', ew fails = fails + 1 @@ -3075,7 +3075,7 @@ subroutine max_eig_tests(options,fails) diff(i) = diff(i) + (tmpA(2) - ew * tmpB(2))**2 diff(i) = sqrt(diff(i)) end do - if (norm2(diff) > 1e-10) then + if (.not. norm2(diff) < 1e-10) then write(*,*) 'error :: hard case of max_eig test failed - wrong vectors returned' write(*,*) 'diff = ', diff fails = fails + 1 @@ -3149,12 +3149,12 @@ subroutine shift_matrix_tests(options,fails) AplusSigma = 0.0_wp sigma = 5.0_wp call shift_matrix(A,sigma,AplusSigma,n) - if ( ( (AplusSigma(1,1)-6.0_wp) > 1e-12) .or. & - ((AplusSigma(2,2) - 6.0_wp) > 1e-12) ) then + if ( .not. ( (AplusSigma(1,1)-6.0_wp) < 1e-12) .and. & + (.not. (AplusSigma(2,2) - 6.0_wp) < 1e-12) ) then write(*,*) 'Error: incorrect return from shift_matrix' fails = fails + 1 - elseif ( ( (AplusSigma(1,2)-1.0_wp) > 1e-12) .or. & - ((AplusSigma(2,1) - 1.0_wp) > 1e-12) ) then + elseif ( .not. ( (AplusSigma(1,2)-1.0_wp) < 1e-12) .and. & + (.not.(AplusSigma(2,1) - 1.0_wp) < 1e-12) ) then write(*,*) 'Error: incorrect return from shift_matrix' fails = fails + 1 end if From 3f1fb412bf9d5c144b83b986971e75651f7b6037 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Wed, 11 Mar 2026 12:14:48 +0000 Subject: [PATCH 08/16] documented lls_solver in options --- libRALFit/doc/common/options.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/libRALFit/doc/common/options.rst b/libRALFit/doc/common/options.rst index 43cfb631..4306798f 100644 --- a/libRALFit/doc/common/options.rst +++ b/libRALFit/doc/common/options.rst @@ -84,6 +84,8 @@ .. |scale_require_increase| replace:: specifies whether or not to require :math:`{\tt D}_{i,i}` to increase before updating it. +.. |lls_solver| replace:: specifies which linear least squares solver to use internally. + .. |more_sorensen_maxits| replace:: if ``nlls_method = 3``, specifies the maximum number of iterations allowed in the More-Sorensen method. .. |more_sorensen_shift| replace:: if ``nlls_method = 3``, specifies the shift to be used in the More-Sorensen method. From 8c14bb19438944865605dd4ffa3b5cf32f13e32d Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 23 Mar 2026 09:19:50 +0000 Subject: [PATCH 09/16] replaced temporary alloc with ?lacpy --- libRALFit/src/ral_nlls_internal.F90 | 2 +- libRALFit/src/ral_nlls_workspaces.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index a31a0ed6..fd004599 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -1803,7 +1803,7 @@ SUBROUTINE dogleg(J,f,hf,g,n,m,Delta,d,normd,options,inform,w) select case (options%model) case (1) ! linear model... - w%solve_LLS_ws%Jlls = reshape(J, [n,m]) + call PREC(lacpy)('A', m, n, J, m, w%solve_LLS_ws%Jlls, m) w%solve_LLS_ws%temp(1:m) = f(1:m) call solve_LLS(w%solve_LLS_ws%Jlls,w%solve_LLS_ws%temp,n,m,inform,w%solve_LLS_ws,options,.false.) if ( inform%status /= 0 ) goto 100 diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index f6736bd5..21d5d665 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -1254,7 +1254,7 @@ subroutine setup_workspace_solve_LLS(n,m,w,options,inform) inform%status = 0 lwork = max(1, min(m,n) + max(min(m,n), 1)*4) - allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(n,m), stat = inform%alloc_status) + allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(m,n), stat = inform%alloc_status) If (inform%alloc_status /= 0) Then Call remove_workspace_solve_LLS(w,options) inform%status = NLLS_ERROR_ALLOCATION From a72e24c78ef126af400825b17b3f2c418c120cf9 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 23 Mar 2026 11:55:44 +0000 Subject: [PATCH 10/16] removed ipiv stack-alloc --- libRALFit/src/ral_nlls_linear.F90 | 10 ++++------ libRALFit/src/ral_nlls_workspaces.F90 | 4 +++- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index d5a52635..8f8470cd 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -47,7 +47,7 @@ subroutine solve_LLS(A, b, n, m, inform, w, options, pd) call solve_posv(A,b,n,inform) else if (n == m) then - call solve_gesv(A,b,n,inform) + call solve_gesv(A,b,n,inform,w) else if (.not. w%allocated) then inform%status = NLLS_ERROR_WORKSPACE_ERROR @@ -94,17 +94,15 @@ subroutine solve_gels(A,b,n,m,inform,w,options) end subroutine solve_gels - subroutine solve_gesv(A,b,n,inform) + subroutine solve_gesv(A,b,n,inform,w) ! Wrapper around LAPACK's ?gesv implicit none real(wp), intent(inout), contiguous :: A(:,:), b(:) INTEGER, INTENT(IN) :: n type(NLLS_inform), INTENT(INOUT) :: inform + type(solve_LLS_work), intent(inout) :: w - ! NB: we never actually use ipiv in the library - integer, dimension(n) :: ipiv - - call PREC(gesv)(n, 1, A, n, ipiv, b, n, inform%external_return) + call PREC(gesv)(n, 1, A, n, w%ipiv, b, n, inform%external_return) if (inform%external_return /= 0 ) then inform%status = NLLS_ERROR_FROM_EXTERNAL inform%external_name = 'lapack_?gesv' diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index 21d5d665..f5c1e949 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -602,6 +602,7 @@ module MODULE_PREC(ral_nlls_workspaces) type, public :: solve_LLS_work ! workspace for subroutine solve_LLS logical :: allocated = .false. real(wp), allocatable :: temp(:), work(:), Jlls(:,:) + integer, allocatable :: ipiv(:) end type solve_LLS_work type, public :: min_eig_symm_work ! workspace for subroutine min_eig_work @@ -1254,7 +1255,7 @@ subroutine setup_workspace_solve_LLS(n,m,w,options,inform) inform%status = 0 lwork = max(1, min(m,n) + max(min(m,n), 1)*4) - allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(m,n), stat = inform%alloc_status) + allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(m,n), w%ipiv(n), stat = inform%alloc_status) If (inform%alloc_status /= 0) Then Call remove_workspace_solve_LLS(w,options) inform%status = NLLS_ERROR_ALLOCATION @@ -1273,6 +1274,7 @@ subroutine remove_workspace_solve_LLS(w,options) if(allocated( w%temp )) deallocate( w%temp, stat=ierr_dummy ) if(allocated( w%work )) deallocate( w%work, stat=ierr_dummy ) if(allocated( w%Jlls )) deallocate( w%Jlls, stat=ierr_dummy ) + if(allocated( w%ipiv )) deallocate( w%ipiv, stat=ierr_dummy ) w%allocated = .false. end subroutine remove_workspace_solve_LLS From 537f2ced2cbf6b2531c12981b7d94b28b5bc675a Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 2 Feb 2026 14:33:35 +0000 Subject: [PATCH 11/16] added LSQR and randomised solvers --- libRALFit/src/CMakeLists.txt | 20 +- libRALFit/src/external/CMakeLists.txt | 9 + libRALFit/src/external/dct.f90 | 35 + libRALFit/src/external/lsqr.f90 | 1494 +++++++++++++++++++++++ libRALFit/src/external/preprocessor.FPP | 1 + libRALFit/src/ral_nlls_internal.F90 | 4 +- libRALFit/src/ral_nlls_linear.F90 | 306 ++++- libRALFit/src/ral_nlls_workspaces.F90 | 122 +- libRALFit/test/nlls_test.F90 | 6 + libRALFit/test/unit_test_mod.F90 | 169 +++ 10 files changed, 2147 insertions(+), 19 deletions(-) create mode 100644 libRALFit/src/external/CMakeLists.txt create mode 100644 libRALFit/src/external/dct.f90 create mode 100644 libRALFit/src/external/lsqr.f90 create mode 120000 libRALFit/src/external/preprocessor.FPP diff --git a/libRALFit/src/CMakeLists.txt b/libRALFit/src/CMakeLists.txt index 23e48031..77318df1 100644 --- a/libRALFit/src/CMakeLists.txt +++ b/libRALFit/src/CMakeLists.txt @@ -19,6 +19,11 @@ set(source_files ral_nlls_ciface.F90 CACHE INTERNAL "source_files" FORCE) +# subdirectory for external routines e.g. LSQR +add_subdirectory(external) +set_property(SOURCE ral_nlls_linear.F90 APPEND PROPERTY COMPILE_OPTIONS "-I./external") +set_property(SOURCE ral_nlls_workspaces.F90 APPEND PROPERTY COMPILE_OPTIONS "-I./external") + # Suppress warnings ############ if (SUPPRESS_INTERNAL_WONTFIX) if(CMAKE_Fortran_COMPILER_ID STREQUAL "GNU") @@ -94,6 +99,9 @@ function(generate_ral_nlls_modules precision obj_list) add_dependencies("RAL_NLLS_${precision}" "RAL_NLLS_INTERNAL_${precision}") add_dependencies("RAL_NLLS_CIFACE_${precision}" "RAL_NLLS_${precision}") add_dependencies("RAL_NLLS_INTERNAL_${precision}" "RAL_NLLS_LINEAR_${precision}") + add_dependencies("RAL_NLLS_LINEAR_${precision}" "dct_module_${precision}") + add_dependencies("RAL_NLLS_LINEAR_${precision}" "lsqr_reverse_${precision}") + add_dependencies("RAL_NLLS_WORKSPACES_${precision}" "lsqr_reverse_${precision}") add_dependencies("RAL_NLLS_LINEAR_${precision}" "RAL_NLLS_WORKSPACES_${precision}") # Replicate tree to be compatible with Ninja @@ -117,6 +125,9 @@ function(generate_ral_nlls_modules precision obj_list) target_link_libraries("RAL_NLLS_${precision}" "RAL_NLLS_INTERNAL_${precision}") target_link_libraries("RAL_NLLS_CIFACE_${precision}" "RAL_NLLS_${precision}") target_link_libraries("RAL_NLLS_INTERNAL_${precision}" "RAL_NLLS_LINEAR_${precision}") + target_link_libraries("RAL_NLLS_LINEAR_${precision}" "dct_module_${precision}") + target_link_libraries("RAL_NLLS_LINEAR_${precision}" "lsqr_reverse_${precision}") + target_link_libraries("RAL_NLLS_WORKSPACES_${precision}" "lsqr_reverse_${precision}") target_link_libraries("RAL_NLLS_LINEAR_${precision}" "RAL_NLLS_WORKSPACES_${precision}") endfunction() @@ -124,13 +135,18 @@ endfunction() set(nlls_obj_list "") if(SINGLE_PRECISION) # Generate single-precision library - generate_ral_nlls_modules("single" obj_list) + set(precision "single") else() # Generate double-precision library - generate_ral_nlls_modules("double" obj_list) + set(precision "double") endif() +generate_ral_nlls_modules(${precision} obj_list) + add_library (ral_nlls-static STATIC ${obj_list}) add_library (ral_nlls SHARED ${obj_list}) +target_link_libraries(ral_nlls) +target_link_libraries(ral_nlls lsqr_reverse_${precision}) +target_link_libraries(ral_nlls dct_module_${precision}) target_link_libraries(ral_nlls ${LIBS} m ${CMAKE_DL_LIBS}) target_include_directories(ral_nlls PUBLIC ${CMAKE_SOURCE_DIR}/include ${CMAKE_BINARY_DIR}/src) # End Post ################ diff --git a/libRALFit/src/external/CMakeLists.txt b/libRALFit/src/external/CMakeLists.txt new file mode 100644 index 00000000..348c3b77 --- /dev/null +++ b/libRALFit/src/external/CMakeLists.txt @@ -0,0 +1,9 @@ +if(SINGLE_PRECISION) + set(precision "single") +else() + set(precision "double") +endif() +add_library(lsqr_reverse_${precision} OBJECT lsqr.f90) +add_library(dct_module_${precision} OBJECT dct.f90) + +target_link_libraries(dct_module_${precision} fftw3) diff --git a/libRALFit/src/external/dct.f90 b/libRALFit/src/external/dct.f90 new file mode 100644 index 00000000..35412a03 --- /dev/null +++ b/libRALFit/src/external/dct.f90 @@ -0,0 +1,35 @@ +! Interface to DCT implementations. +#include "../preprocessor.FPP" + +module MODULE_PREC(dct_module) + implicit none + + ! these parameters are defined in FFTW, but we define them here to avoid including the FFTW header + integer(8) :: FFTW_REDFT10 + parameter (FFTW_REDFT10 = 4) + integer(8) :: FFTW_ESTIMATE + parameter (FFTW_ESTIMATE = 64) + integer(4), parameter :: wp = kind( 0.0d+0 ) + + interface dct + module procedure dct_fftz + end interface dct + + contains + + subroutine dct_fftz(A, n, m) + ! Compute the 1D DCT across the columns of A using the FFTZ library. + real(8), intent(inout), contiguous :: A(:,:) + integer, intent(in) :: n, m + + integer :: i, plan + + call dfftw_plan_r2r_1d(plan, m, A(:, 1), A(:, 1), FFTW_REDFT10, FFTW_ESTIMATE) + do i = 1, n + call dfftw_execute_r2r(plan, A(:, i), A(:, i)) + end do + call dfftw_destroy_plan(plan) + + end subroutine dct_fftz + +end module MODULE_PREC(dct_module) diff --git a/libRALFit/src/external/lsqr.f90 b/libRALFit/src/external/lsqr.f90 new file mode 100644 index 00000000..0406b50c --- /dev/null +++ b/libRALFit/src/external/lsqr.f90 @@ -0,0 +1,1494 @@ +! July 2025 +! Reverse communication variant of LSQR. +! Optionally allows split or left preconditioning of the normal equations. +! The default stopping is to use the stopping criteria of Papez and Tichy. +! The user can test for convergence +! (or terminate after a chosen number of iterations). + +! Also, allows for (partial or selective) one-sided reorthogonalization. + +! Here, we base notation on the paper +! Preconditioning of LSQR and CGLS: Variants, Properties and Relations. +! Havelková E, Hnětynková I. +! In: European Conference on Numerical Mathematics and Advanced Applications +! (2023), pp. 415-424). Springer Nature Switzerland. +! We implement RP-LSQR, Modified RP-LSQR and LP-LSQR. + +! Stopping criteria of Papez and Tichy is described in +! Estimating error norms in CG-like algorithms for least-squares and +! least-norm problems. Jan Papez · Petr Tichy. Numerical Algorithms +! (2024) 97:1--28. + +! Interface designed to be as for CGLS code, but with inclusion of +! an optional damping parameter (regularization). The action parameter +! has the same meaning so easy for user to call LSQR or CGLS +! using same driver. + +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +! LSQR solves +! min_x || b - Ax || +! +! min_x ||( A )*x - ( b ) || (regularized LS) +! ||( damp*I ) ( 0 ) || + +! The parameter damp is intended to help regularize +! ill-conditioned systems, by preventing the true solution from +! being very large. + +! If M = PP^T then right preconditioned LSQR solves +! min_y ||b - AP^{-1} xhat|| +! + +! Right preconditioning corresponds to symmetric (split) preconditioning +! of the normal equations +! +! P^{-1} A^T A P^{-T} xhat = P^{-1} A^T b, P^T x_{LS} = xhat + +! One choice is P = L, an incomplete factor of the normal matrix A^T A. +! This can be computed using HSL_MI35. +! +! The code also offers factorization-free preconditioned LSQR. This solves +! the LS problem by solving the left preconditioned normal equations +! +! M^{-1} A^T A x_{LS} = M^{-1} A^T b + +! where M \approx A^T A. Again, M = LL^T can be used but the factored form of +! M is NOT necessary. + + +! Note that x is not an input parameter. + ! If some initial estimate x0 is known and if damp = 0, + ! one could proceed as follows: + ! + ! 1. Compute a residual vector r0 = b - A*x0. + ! 2. Use LSQR to solve the system A*dx = r0. + ! 3. Add the correction dx to obtain a final solution x = x0 + dx. + ! + ! This requires that x0 be available before the first call + ! to LSQR and after the final call. +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +module lsqr_reverse_double + + implicit none + + private + public :: lsqr_keep, lsqr_options, lsqr_inform + public :: lsqr, lsqr_free + + integer(4), parameter :: ip = kind( 0 ) + integer(4), parameter :: wp = kind( 0.0d+0 ) + integer, parameter :: long = selected_int_kind(18) + real(wp), parameter :: zero = 0.0_wp, one = 1.0_wp + + ! error flags(explained below) + integer(ip), parameter :: lsqr_stop_m_oor = -1 + integer(ip), parameter :: lsqr_stop_allocation = -2 + integer(ip), parameter :: lsqr_stop_deallocation = -3 + integer(ip), parameter :: lsqr_stop_itnlim = -4 + + ! warning flags(explained below) + integer(ip), parameter :: lsqr_stop_anorm = 1 + integer(ip), parameter :: lsqr_stop_x0 = 2 + integer(ip), parameter :: lsqr_stop_stagnate = 4 + ! positive values are added if more than one warning issued + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! Make interfaces generic. + interface lsqr + module procedure lsqr_double + end interface + + interface lsqr_free + module procedure lsqr_free_double ! deallocates components of keep + end interface + +! ------------------------------------------------------------------ + + ! The matrices A, P and M are treated as linear operators that are + ! accessed by reverse communication, that is, requests for matrix-vector + ! products with A, P or M are passed back to the user. + ! + ! LSQR uses an iterative method to approximate the solution. + ! The number of iterations required to reach a certain accuracy + ! depends strongly on the scaling of the problem and on the preconditioner. + + ! Unless better information is known, we recommend that + ! the columns of A are prescaled so that they all have + ! the same Euclidean norm (e.g., 1.0). + ! + +! ------------------------------------------------------------------ + ! control data type + ! These may be set before first call and must not be altered between calls. + type :: lsqr_options + + integer(ip) :: itnlim = -1 ! max number of iterations. + ! This must be set because Papez and Tichy stopping requires arrays of + ! size the number of iterations. + ! If it is negative then we will use itnlim = n. + + integer(ip) :: test_frequency = 1 + ! Return to user for convergence testing every test_frequency iterations. + ! Values less than 1 are treated as 1 + + integer(ip) :: localSize = 0 + ! No. of vectors for local reorthogonalization. Not used + ! if options%orthog_choice = 0 + ! 0 No reorthogonalization is performed. + ! >0 This many n-vectors "v" are saved for reorthogonalizing. + ! If localSize = min(m,n) then full one-sided reorth'n is done. + ! min(options%localSize,m,n) vectors will be allocated or, if + ! itnlim > 0, min(options%localSize,options%itnlim,m,n). + ! See also options%orthog_tol + + integer(ip) :: stop_test = 1 ! set up data required for + ! using the Papez and Tichy stopping test. + ! The user performs the test for convergence + + ! stop_test = 2: this is as for stop_test = 1 except that + ! the code performs the Papez and Tichy stopping test. + ! For this, the user is required to supply an estimate of ||A||_2 using + ! the optional input argument anorm. If this is not present, stop_test=1 + ! is used. Note: stop_test = 2 should NOT be used if precon = 0 + ! and split preconditioning is used (so that A is replaced by A*L^{-1}). + ! In this case, the user must test for convergence (otherwise, the + ! stopping test is for the preconditioned problem). + + ! if stop_test = 0, then this data is not set up (saves memory and work but + ! user then has to decide how to terminate the computation). + + ! other values treated as stop_test = 1 + + real(wp) :: delta = sqrt(epsilon(1.0_wp)) ! convergence tolerance + ! if stop_test = 2 + + integer(ip) :: orthog_choice = 0 ! controls reorthog strategy + ! 0 = no reorthogonalisation + ! 1 = reorthogonalise using MGS against previous options%localSize vectors + ! 2 = reorthogonalise using MGS against first options%localSize vectors + !-1 = reorthogonalise using MGS2 against previous options%localSize vectors + !-2 = reorthogonalise using MGS2 against first options%localSize vectors + ! All other values are treated as zero + ! MGS = modified Gram Schmidt + ! MGS2 = two applications of MGS (can help but can add significant overhead) + + real(wp) :: orthog_tol = zero ! if orthog_tol > zero then + ! selective reorthogonalization (ie only reorthog v against a + ! previous vector if inner product between them is greater + ! than orthog_tol). MGS2 is not used with this option (so abs(orthog_choice) + ! is used). + + real(wp) :: stagnation = zero ! stagnation occurs + ! if the update to x on some iteration is less than options%stagnation + + ! The following are used by Papez and Tichy stopping (these values were used + ! in their paper). Only accessed on first call. Values .le. zero are reset + ! to the default. + real(wp) :: tau = 0.25_wp ! + real(wp) :: tol = 0.0001_wp ! + + end type lsqr_options + +! ------------------------------------------------------------------ + + ! information data type + type :: lsqr_inform + + integer(ip) :: flag ! Gives reason for termination + ! negative indicates failure + + ! lsqr_stop_m_oor : n < 1 or m < 1 + ! lsqr_stop_allocation : An array allocation failed. + ! lsqr_stop_deallocation : An array deallocation failed. + ! lsqr_stop_itnlim : iteration limit reached + ! x holds the current solution. + + ! warning flags(explained below) + ! lsqr_stop_anorm : options%stop_test = 2 only. + ! Either the user did not supply anorm or anorm.le.0 + ! lsqr_stop_x0 : x = 0 is the exact solution. + ! No iterations were performed. + ! lsqr_stop_stagnate : computation of LS solution has stagnated. + ! User should test for convergence (options%stop_test.ne.2) + ! and then terminate (even if convergence test not satisfied) + + integer(ip) :: itn ! The number of iterations performed + + integer :: stat ! Fortran stat parameter + + real(wp) :: estim ! set if options%stop_test = 1 or 2. In this case, holds + ! backward error estimate that is required by Papez-Tichy stopping test. + + real :: time_RO ! accumulate the time for performing reorthogonalization. + + end type lsqr_inform + +! ------------------------------------------------------------------ + +! Derived type to ensure local variables are saved safely in the +! reverse communication + type :: lsqr_keep + private + + logical :: damped, localVQueueFull + + integer(ip) :: branch = 0 + integer(ip) :: flag = 0 + integer(ip) :: itn_count ! iteration count + integer(ip) :: itnlim + integer(ip) :: localPointer, localVecs, orthog_choice + integer(ip) :: precon + integer(ip) :: stop_test ! copy of options%stop_test + integer(ip) :: test_frequency ! testing frequency + + real(wp) :: alpha, beta, rho, rhobar + real(wp) :: backward + real(wp) :: anorm, bnorm, xnorm + real(wp) :: orthog_tol ! copy of options%orthog_tol + real(wp) :: damp + real(wp) :: phibar + real(wp) :: stagnation ! copy of options%stagnation + real(wp) :: t1, t2 + + real(wp), allocatable :: w(:) + real(wp), allocatable :: localV(:,:) + + ! These are only needed if stop_test = 1 or 2 + ! (Papez and Tichy stopping condition) + integer(ip) :: d, ell, estimend + real(wp) :: tau, tol ! copies of options%tau, options%tol + + real(wp), allocatable, dimension(:) :: delay + real(wp), allocatable, dimension(:) :: delta + real(wp), allocatable, dimension(:) :: curve + real(wp), allocatable, dimension(:) :: estim + + end type lsqr_keep + +contains + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & + options, inform, damp, anorm) + + integer(ip), intent(in) :: precon ! This parameter controls + ! the use of preconditioning. It must be set on the first call + ! to one of the following values: + + ! 0: No preconditioning. Also can be used for right preconditioning by + ! replacing products with A by products with A* P^{-T} + ! where the preconditioner is M = PP^T. In this case, + ! the least squares solution must be recovered by the user + ! as x_LS = P^T x, where x is returned by lsqr (if testing is + ! to be done for the unpreconditioned problem, this will need to be + ! performed each time action = 4 is returned). + ! The advantage compared to precon = 1 is that the array z is not + ! needed (and so can be of size 1). + ! 1: Right (split) preconditioning. Requires the preconditioner + ! is in factored form M = PP^T eg if A^T*A = LL^T (Cholesky fact) + ! can use P = L + ! 2: Left preconditioning of normal equations. M does not need to + ! be in factored form. + + ! Note: precon is only accessed on the first call (the call with + ! action = 0). If the user wishes to try different options, then + ! the computation must be restarted. + + integer(ip), intent(inout) :: action ! This parameter controls + ! the action. On initial entry, action must be set to 0 + ! and the initial guess for the solution must be in the array x + ! (if no guess is available, the user should set x(1:n) = 0). + ! On subsequent calls, action should be unchanged by the user. + ! On return, action determines what the user is required to do. + ! Possible values and their consequences are: + ! + ! 0. Computation has terminated. + ! action = 0 is returned if an error has been encountered. + ! Also returned if options%stop_test = 2 and convergence + ! has been achieved + ! + ! 1. precon = 0 or 2 : the user must compute v = v + A^T *u + ! precon = 1 : the user must compute v = v + P^{-1} A^T *u + ! The user must then re-call the subroutine without any other arguments. + ! The vectors u and v are available in the arrays u and v (see below). + ! This action is returned once per iteration. + + ! 2. precon = 0 : the user must compute u = u + A *v + ! precon = 1 or 2 : the user must compute u = u + A *z + ! The user must then re-call the subroutine without any other arguments. + ! The vectors u, v and z are available in the arrays u, v and z (see below). + ! This action is returned once per iteration. + + ! 3. precon = 1 : the user must compute z = P^{-T}*v + ! precon = 2 : the user must compute z = M^{-1}*v + ! The user must then re-call the subroutine without any other arguments. + ! The vectors z and v are available in the arrays z and v (see below). + ! This action is returned once per iteration. + + ! 4: The user may test for convergence. + ! Re-call the subroutine WITHOUT alterning any of the arguments. + ! This return does not happen if options%stop_test = 2 + + integer(ip), intent(in) :: m ! the number of rows in A. + + integer(ip), intent(in) :: n ! the number of columns in A. + + real(wp), intent(inout) :: u(m) ! Prior to the first call, u + ! must hold the rhs vector b. u is then the vector used to communicate u + ! when further action is required (see action, above). + + real(wp), intent(inout) :: v(n) ! The vector used to communicate v + ! when further action is required (see action, above). + + real(wp), intent(inout) :: z(:) ! The vector used to communicate z + ! when further action is required (see action, above). + ! z is not accessed if precon = 0 (so can be size 1 in this case; + ! otherwise it should be length n) + + real(wp), intent(inout) :: x(n) ! Does not need to be set by the user. + ! Returns the computed solution x + + type ( lsqr_keep ), intent( inout) :: keep + ! A variable of type lsqr_keep that is used to + ! preserve internal data between reverse-communication + ! calls. The components are private. + + type ( lsqr_options ), intent( inout) :: options + ! A variable of type lsqr_options that is used to control the options. + ! It should not be altered by the user between calls. + + type ( lsqr_inform ), intent( inout) :: inform + ! A variable of type lsqr_info_type that is used to hold information. + ! It should not be altered by the user (so as to preserve the + ! information between calls). + + real(wp),intent(in), optional :: damp ! The damping parameter. + ! Only accessed on the first call. damp should be > 0 + ! Note: The work per iteration and the storage needed + ! by LSQR are the same for all values of damp. + + real(wp),intent(in), optional :: anorm ! The norm of A. + ! must be present and hold an estimate of the 2 norm of A + ! if options%stop_test = 2 + + + real(wp) :: ddot, dnrm2 + ! Local arrays and variables + integer :: dst, st + real(wp) :: backward + real(wp) :: cs, cs1, phi, psi, rhbar1, rho, sn, sn1 + real(wp) :: theta + real(wp) :: xnorm + + !------------------------------------------------------------------- + + ! on first call, initialize keep%branch and keep%flag + ! And initial components of inform so they are not undefined. + if (action.eq.0) then + keep%branch = 0 + keep%flag = 0 + + inform%flag = 0 + inform%itn = 0 + inform%stat = 0 + inform%estim = -1.0_wp + inform%time_RO = 0.0 + end if + + ! Immediate return if we have already had an error + if (keep%flag.lt.0) then + inform%flag = keep%flag + action = 0 + return + end if + + ! also terminate if we things have stagnated. + if (keep%flag.ge.lsqr_stop_stagnate) then + action = 0 + return + end if + + ! on other calls, jump to the appropriate place after reverse-communication + + if (precon.eq.0) then + ! LSQR without preconditioning (or split preconditioning with + ! A replaced by AP^{-T}) + select case ( keep%branch ) + case ( 1 ) + go to 20 + case ( 2 ) + go to 40 + case ( 3 ) + go to 50 + case ( 4 ) + go to 60 + end select + + else if (precon.eq.1) then + ! modified RP-LSQR + select case ( keep%branch ) + case ( 1 ) + go to 120 + case ( 2 ) + go to 140 + case ( 3 ) + go to 150 + case ( 4 ) + go to 160 + case ( 5 ) + go to 125 + case ( 6 ) + go to 170 + end select + + else if (precon.eq.2) then + ! LP-LSQR + select case ( keep%branch ) + case ( 1 ) + go to 220 + case ( 2 ) + go to 240 + case ( 3 ) + go to 250 + case ( 4 ) + go to 280 + case ( 5 ) + go to 225 + case ( 6 ) + go to 270 + end select + end if + + ! Initialize. + keep%damped = .false. + if (present(damp)) then + keep%damp = damp + keep%damped = keep%damp > zero + end if + + keep%itnlim = options%itnlim + if (options%itnlim < 1) keep%itnlim = n + + keep%stop_test = options%stop_test + if (keep%stop_test.eq.2) then + if (.not. present(anorm)) then + keep%stop_test = 1 + keep%flag = lsqr_stop_anorm + else if (anorm .le. zero) then + keep%stop_test = 1 + keep%flag = lsqr_stop_anorm + else + keep%anorm = anorm + end if + else if (keep%stop_test .ne. 0) then + keep%stop_test = 1 + end if + inform%flag = keep%flag + + keep%itn_count = 0 + keep%test_frequency = options%test_frequency + if (keep%test_frequency .lt. 1) keep%test_frequency = 1 + + !!! orthogonalisation controls + keep%orthog_choice = options%orthog_choice + if (abs(keep%orthog_choice).gt.2) keep%orthog_choice = 0 + + + keep%localVecs = min(options%localSize,m,n) + if (options%itnlim > 0) & + keep%localVecs = min(keep%localVecs, options%itnlim) + if (keep%orthog_choice.eq.0) keep%localVecs = 0 + + keep%orthog_tol = max(zero, options%orthog_tol) + ! Do not use MGS2 if we have orthog_tol non zero + if (keep%orthog_tol.ne.zero) keep%orthog_choice = abs(keep%orthog_choice) + !!!! + + keep%stagnation = max(zero, options%stagnation) + + keep%tau = options%tau + if (keep%tau .le. zero) keep%tau = 0.25_wp + + keep%tol = options%tol + if (keep%tol .le. zero) keep%tol = 0.0001_wp + + st = 0 + dst = 0 + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! quick check of m and n + if (n.lt.1 .or. m.lt.1) then + keep%flag = lsqr_stop_m_oor + inform%flag = lsqr_stop_m_oor + return + end if + + ! if precon is out-of-range, run without preconditioning + keep%precon = precon + if (precon.lt.1 .or. precon.gt.2) keep%precon = 0 + + ! allocate workarrays + + if (.not. allocated(keep%w)) then + allocate( keep%w(n), stat = st) + if (st /= 0) go to 10 + else if (size(keep%w).lt.n) then + deallocate (keep%w, stat = dst) + if (dst /= 0) then + go to 10 + else + allocate( keep%w(n), stat = st) + if (st /= 0) go to 10 + end if + end if + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! set up buffer for optional reorthogonalization + if (keep%localVecs > 0) then + if (.not. allocated(keep%localV)) then + allocate( keep%localV(n,keep%localVecs), stat = st) + if (st /= 0) go to 10 + else if (size(keep%localV,1).lt.n .or. & + size(keep%localV,2).lt.keep%localVecs) then + deallocate (keep%localV, stat = dst) + if (dst /= 0) go to 10 + allocate( keep%localV(n,keep%localVecs), stat = st) + if (st /= 0) go to 10 + end if + end if + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! set up for Papez and Tichy stopping + if (keep%stop_test.ne.0) then + keep%ell = 1 + keep%d = 0 + keep%estimend = -1 + + if (allocated(keep%curve)) deallocate(keep%curve, stat = dst) + if (dst /= 0) go to 10 + if (allocated(keep%delay)) deallocate(keep%delay, stat = dst) + if (dst /= 0) go to 10 + if (allocated(keep%delta)) deallocate(keep%delta, stat = dst) + if (dst /= 0) go to 10 + if (allocated(keep%estim)) deallocate(keep%estim, stat = dst) + if (dst /= 0) go to 10 + + allocate (keep%curve(keep%itnlim), keep%delay(keep%itnlim), & + keep%delta(keep%itnlim), keep%estim(keep%itnlim), stat = st) + if (st /= 0) go to 10 + end if + +10 continue + if (st.ne.0) then + inform%stat = st + inform%flag = lsqr_stop_allocation + return + else if (dst.ne.0) then + inform%stat = dst + inform%flag = lsqr_stop_deallocation + return + end if + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + keep%bnorm = dnrm2 (m, u, 1) + if (keep%bnorm .eq. zero) then + ! Exit if b = 0. Solution is x_{LS} = 0. + x(1:n) = zero + keep%flag = keep%flag + lsqr_stop_x0 + inform%flag = keep%flag + return + end if + + !------------------------------------------------------------------- + ! Set up the first vectors u_1 and v_1 for the bidiagonalization. + ! These satisfy beta*u = r_0 = b, alpha*v = A^T *u. + ! u holds r_0 = b - A*x_0. + !------------------------------------------------------------------- + + keep%beta = dnrm2 (m, u, 1) + + call dscal (m, (one/keep%beta), u, 1) ! u_1 = u / beta + + ! User must compute v = A^T * u (precon = 0 or 2), + ! or v = P^{-1} A^T *u (precon = 1) + + ! initialise v = 0 (as user may have routine for computing v = v + A^T * u) + v(1:n) = zero + + action = 1 + keep%branch = 1 + return + + !=================================================================== + !=================================================================== + ! precon = 0 : This is the no preconditioning variant + ! or split conditioning with A replaced by AP^{-T} + + 20 continue + + ! v = A^T *u_1 has been performed by the user + + keep%alpha = dnrm2 (n, v, 1) + + if (keep%alpha .eq. zero) then + ! terminate A'b = 0 + action = 0 + return + end if + + call dscal (n, (one/keep%alpha), v, 1) + + ! w_1 = v_1 + keep%w = v + + ! Initialization for local reorthogonalization. + if (keep%localVecs > 0) then + keep%localPointer = 1 + keep%localVQueueFull = .false. + keep%localV(1:n,1) = v(1:n) + end if + + ! Initialize variables for 1st iteration. + keep%rhobar = keep%alpha + keep%phibar = keep%beta + + !=================================================================== + ! Main iteration loop. + !=================================================================== + + 30 continue + + keep%itn_count = keep%itn_count + 1 + inform%itn = keep%itn_count + + !---------------------------------------------------------------- + ! Perform the next step of the bidiagonalization to obtain the + ! next beta, u, alpha, v. These satisfy + ! beta*u = A*v - alpha*u, + ! alpha*v = (A)'*u - beta*v. + !---------------------------------------------------------------- + call dscal (m,(- keep%alpha), u, 1) ! - alpha*u_j + + ! Compute u = u + A*v + action = 2 + keep%branch = 2 + return + + 40 continue + + ! user has computed u_{j+1} = A*v_j - alpha *u_j. Normalise. + keep%beta = dnrm2 (m, u, 1) + + if (keep%beta > zero) then + call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + if (keep%localVecs > 0) then + ! Check whether to store v_j for local reorthog'n of v_{j+1} + if (abs(keep%orthog_choice).eq.1 .or. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + + ! Compute v = v + A'*u + action = 1 + keep%branch = 3 + return + end if + + 50 continue + + ! user has computed v_{j+1} = A^T*u_{j+1} - beta *v_j. + if (keep%beta > zero) then + if (keep%localVecs > 0) & ! Local reorthogonalization of new v. + call localVOrtho(n, v, keep, inform%time_RO) + + keep%alpha = dnrm2 (n, v, 1) + if (keep%alpha > zero) call dscal (n, (one/keep%alpha), v, 1) + + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the damping parameter. + ! This alters the diagonal (rhobar) of the lower-bidiagonal matrix. + !---------------------------------------------------------------- + rhbar1 = keep%rhobar + psi = zero + if (keep%damped) then + rhbar1 = d2norm(keep%rhobar, keep%damp) + cs1 = keep%rhobar/rhbar1 + sn1 = keep%damp /rhbar1 + psi = sn1 * keep%phibar + keep%phibar = cs1 * keep%phibar + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the subdiagonal element (beta) + ! of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. + !---------------------------------------------------------------- + rho = d2norm(rhbar1, keep%beta) + cs = rhbar1 /rho + sn = keep%beta /rho + theta = sn*keep%alpha + keep%rhobar = - cs*keep%alpha + phi = cs*keep%phibar + keep%phibar = sn*keep%phibar + + !---------------------------------------------------------------- + ! Update x + ! --------------------------------------------------------------- + keep%t1 = phi/rho + keep%t2 = - theta/rho + x(1:n) = keep%t1* keep%w(1:n) + x(1:n) + + ! the following enables Papez and Tichy stopping test + if (keep%stop_test.ne.0) then + call adaptive(keep%itn_count,phi*phi,keep) + if (keep%estimend > 0) then + inform%estim = keep%estim(keep%estimend) + end if + end if + + !---------------------------------------------------------------- + ! Test for convergence. And also check for stagnation. + !---------------------------------------------------------------- + + if (keep%itn_count >= keep%itnlim) then ! iteration count exceeded + inform%flag = lsqr_stop_itnlim + keep%flag = lsqr_stop_itnlim + action = 0 + return + end if + + if (abs(keep%t1).le.keep%stagnation) then + ! Stagnated. + if (keep%stop_test.ne.2) then + ! The user should check for convergence. + if (keep%flag .le. lsqr_stop_stagnate) then + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + end if + action = 4 + keep%branch = 4 + return + else + ! Perform Papez and Tichy testing + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + else + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + action = 0 + return + end if + end if + end if + + ! see if it is time to return to test convergence + if (mod(keep%itn_count,keep%test_frequency) == 0) then + if (keep%stop_test.ne.2) then + action = 4 + keep%branch = 4 + return + end if + end if + if (keep%stop_test.eq.2) then + ! Perform Papez and Tichy testing (do this every iteration + ! as involves little extra work) + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + end if + end if + + 60 continue + + !---------------------------------------------------------------- + ! Update w + ! --------------------------------------------------------------- + + keep%w(1:n) = keep%t2* keep%w(1:n) + v(1:n) + + ! cycle for next iteration + go to 30 + + !=================================================================== + !=================================================================== + ! This is the modified RP-LSQR preconditioning variant (precon = 1) + +120 continue + + ! v = P^{-1} A^T *u has been performed by the user + + keep%alpha = dnrm2 (n, v, 1) + + if (keep%alpha .eq. zero) then + ! terminate P^{-1} A^T b = 0 + action = 0 + return + end if + + call dscal (n, (one/keep%alpha), v, 1) + + !---------------------------------------------------------------- + ! user to compute z_1 = P^{-T} v_1 + !---------------------------------------------------------------- + + action = 3 + keep%branch = 5 + return + + 125 continue + + ! d_1 = z_1 + keep%w = z + + ! Initialization for local reorthogonalization. + if (keep%localVecs > 0) then + keep%localPointer = 1 + keep%localVQueueFull = .false. + keep%localV(1:n,1) = v(1:n) + end if + + ! Initialize variables for 1st iteration. + keep%rhobar = keep%alpha + keep%phibar = keep%beta + + !=================================================================== + ! Main iteration loop. + !=================================================================== + +130 continue + + keep%itn_count = keep%itn_count + 1 + inform%itn = keep%itn_count + + !---------------------------------------------------------------- + ! Perform the next step of the bidiagonalization to obtain the + ! next beta, u, alpha, v. These satisfy + ! beta*u = A*z - alpha*u, + ! alpha*v = P^{-1} A^T*u - beta*v. + !---------------------------------------------------------------- + call dscal (m,(-keep%alpha), u, 1) ! -alpha*u_j + + ! Compute u = u + A*z + action = 2 + keep%branch = 2 + return + +140 continue + + ! user has computed u_{j+1} = A*z_j - alpha *u_j. Normalise. + keep%beta = dnrm2 (m, u, 1) + + if (keep%beta > zero) then + call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + if (keep%localVecs > 0) then + ! Check whether to store v_j for local reorthog'n of v_{j+1} + if (abs(keep%orthog_choice).eq.1 .or. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + + ! Compute v = v + P^{-1} A^T*u + action = 1 + keep%branch = 3 + return + end if + +150 continue + + ! user has computed v_{j+1} = A^T*u_{j+1} - beta *v_j. + if (keep%beta > zero) then + + if (keep%localVecs > 0) & ! Local reorthogonalization of new v. + call localVOrtho(n, v, keep, inform%time_RO) + + keep%alpha = dnrm2 (n, v, 1) + if (keep%alpha > zero) call dscal (n, (one/keep%alpha), v, 1) + + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the damping parameter. + ! This alters the diagonal (rhobar) of the lower-bidiagonal matrix. + !---------------------------------------------------------------- + rhbar1 = keep%rhobar + psi = zero + if (keep%damped) then + rhbar1 = d2norm(keep%rhobar, keep%damp) + cs1 = keep%rhobar/rhbar1 + sn1 = keep%damp /rhbar1 + psi = sn1 * keep%phibar + keep%phibar = cs1 * keep%phibar + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the subdiagonal element (beta) + ! of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. + !---------------------------------------------------------------- + rho = d2norm(rhbar1, keep%beta) + cs = rhbar1 /rho + sn = keep%beta /rho + theta = sn*keep%alpha + keep%rhobar = - cs*keep%alpha + phi = cs*keep%phibar + keep%phibar = sn*keep%phibar + + !---------------------------------------------------------------- + ! Update x + ! --------------------------------------------------------------- + keep%t1 = phi/rho + keep%t2 = - theta/rho + + x(1:n) = keep%t1* keep%w(1:n) + x(1:n) + + ! the following enables Papez and Tichy stopping test + if (keep%stop_test.ne.0) then + call adaptive(keep%itn_count,phi*phi,keep) + if (keep%estimend > 0) & + inform%estim = keep%estim(keep%estimend) + end if + + !---------------------------------------------------------------- + ! Test for convergence. And also check for stagnation. + !---------------------------------------------------------------- + + if (keep%itn_count >= keep%itnlim) then ! iteration count exceeded + inform%flag = lsqr_stop_itnlim + keep%flag = lsqr_stop_itnlim + action = 0 + return + end if + + if (abs(keep%t1).le.keep%stagnation) then + ! Stagnated. + if (keep%stop_test.ne.2) then + ! The user should check for convergence. + if (keep%flag .le. lsqr_stop_stagnate) then + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + end if + action = 4 + keep%branch = 4 + return + else + ! Perform Papez and Tichy testing + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + else + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + action = 0 + return + end if + end if + end if + + ! see if it is time to return to test convergence + if (mod(keep%itn_count,keep%test_frequency) == 0) then + if (keep%stop_test.ne.2) then + action = 4 + keep%branch = 4 + return + end if + end if + if (keep%stop_test.eq.2) then + ! Perform Papez and Tichy testing (do this every iteration + ! as involves little extra work) + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + end if + end if + + 160 continue + !---------------------------------------------------------------- + ! User to compute z = P^{-T} *v + !---------------------------------------------------------------- + action = 3 + keep%branch = 6 + return + + 170 continue + + !---------------------------------------------------------------- + ! Update d + ! --------------------------------------------------------------- + + keep%w(1:n) = keep%t2* keep%w(1:n) + z(1:n) + + ! cycle for next iteration + go to 130 + + + !=================================================================== + !=================================================================== + ! This is the LP-LSQR (left) preconditioning variant (precon = 2). + ! Uses M^{-1} (factored form not needed + +220 continue + + ! v = A^T *u has been performed by the user + + !---------------------------------------------------------------- + ! user to compute z_1 = M^{-1} v_1 + !---------------------------------------------------------------- + + action = 3 + keep%branch = 5 + return + + 225 continue + + keep%alpha = sqrt( ddot (n, v, 1, v, 1)) + + if (keep%alpha .eq. zero) then + ! terminate A^T b = 0 + action = 0 + return + end if + + call dscal (n, (one/keep%alpha), v, 1) + call dscal (n, (one/keep%alpha), z, 1) + + ! d_1 = z_1 + keep%w = z + + ! Initialization for local reorthogonalization. + if (keep%localVecs > 0) then + keep%localPointer = 1 + keep%localVQueueFull = .false. + keep%localV(1:n,1) = v(1:n) + end if + + ! Initialize variables for 1st iteration. + keep%rhobar = keep%alpha + keep%phibar = keep%beta + + !=================================================================== + ! Main iteration loop. + !=================================================================== + +230 continue + + keep%itn_count = keep%itn_count + 1 + inform%itn = keep%itn_count + + !---------------------------------------------------------------- + ! Perform the next step of the bidiagonalization to obtain the + ! next beta, u, alpha, v. These satisfy + ! beta*u = A*z - alpha*u, + ! alpha*v = (A)'*u - beta*v. + !---------------------------------------------------------------- + call dscal (m,(- keep%alpha), u, 1) ! - alpha*u_j + ! Compute u = u + A*z + action = 2 + keep%branch = 2 + return + +240 continue + + ! user has computed u_{j+1} = A*z_j - alpha *u_j. Normalise. + keep%beta = dnrm2 (m, u, 1) + + if (keep%beta > zero) then + call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + if (keep%localVecs > 0) then + ! Check whether to store v_j for local reorthog'n of v_{j+1} + if (abs(keep%orthog_choice).eq.1 .or. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + + ! Compute v = v + (A)'*u + action = 1 + keep%branch = 3 + return + end if + +250 continue + + !---------------------------------------------------------------- + ! User to compute z = M^{-1} *v + !---------------------------------------------------------------- + if (keep%localVecs > 0) & ! Local-reorthogonalization of new v. + call localVOrtho(n, v, keep, inform%time_RO) + + action = 3 + keep%branch = 6 + return + + 270 continue + + ! user has computed v_{j+1} = A^T*u_{j+1} - beta *v_j and + ! z_{j+1} = M^{-1} *v_{j+1} + if (keep%beta > zero) then + + keep%alpha = sqrt( ddot(n, v, 1, z, 1)) + if (keep%alpha > zero) then + call dscal (n, (one/keep%alpha), v, 1) + call dscal (n, (one/keep%alpha), z, 1) + end if + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the damping parameter. + ! This alters the diagonal (rhobar) of the lower-bidiagonal matrix. + !---------------------------------------------------------------- + rhbar1 = keep%rhobar + psi = zero + if (keep%damped) then + rhbar1 = d2norm(keep%rhobar, keep%damp) + cs1 = keep%rhobar/rhbar1 + sn1 = keep%damp /rhbar1 + psi = sn1 * keep%phibar + keep%phibar = cs1 * keep%phibar + end if + + !---------------------------------------------------------------- + ! Use a plane rotation to eliminate the subdiagonal element (beta) + ! of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix. + !---------------------------------------------------------------- + rho = d2norm(rhbar1, keep%beta) + cs = rhbar1 /rho + sn = keep%beta /rho + theta = sn*keep%alpha + keep%rhobar = - cs*keep%alpha + phi = cs*keep%phibar + keep%phibar = sn*keep%phibar + + !---------------------------------------------------------------- + ! Update x + ! --------------------------------------------------------------- + keep%t1 = phi/rho + keep%t2 = - theta/rho + + x(1:n) = keep%t1* keep%w(1:n) + x(1:n) + + ! the following enables Papez and Tichy stopping test + if (keep%stop_test.ne.0) then + call adaptive(keep%itn_count,phi*phi,keep) + if (keep%estimend > 0) & + inform%estim = keep%estim(keep%estimend) + end if + + !---------------------------------------------------------------- + ! Test for convergence. And also check for stagnation. + !---------------------------------------------------------------- + + if (keep%itn_count >= keep%itnlim) then ! iteration count exceeded + inform%flag = lsqr_stop_itnlim + keep%flag = lsqr_stop_itnlim + action = 0 + return + end if + + if (abs(keep%t1).le.keep%stagnation) then + ! Stagnated. + if (keep%stop_test.ne.2) then + ! The user should check for convergence. + if (keep%flag .le. lsqr_stop_stagnate) then + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + end if + action = 4 + keep%branch = 4 + return + else + ! Perform Papez and Tichy testing + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + else + keep%flag = keep%flag + lsqr_stop_stagnate + inform%flag = keep%flag + action = 0 + return + end if + end if + end if + + ! see if it is time to return to test convergence + if (mod(keep%itn_count,keep%test_frequency) == 0) then + if (keep%stop_test.ne.2) then + action = 4 + keep%branch = 4 + return + end if + end if + if (keep%stop_test.eq.2) then + ! Perform Papez and Tichy testing (do this every iteration + ! as involves little extra work) + xnorm = dnrm2 (n, x, 1) + backward = inform%estim /(keep%anorm*xnorm + keep%bnorm) + ! write (*,*) xnorm, backward + if (backward > 0 .and. backward < options%delta) then + ! convergence achieved + action = 0 + return + end if + end if + + 280 continue + !---------------------------------------------------------------- + ! Update d + ! --------------------------------------------------------------- + + keep%w(1:n) = keep%t2* keep%w(1:n) + z(1:n) + + ! cycle for next iteration + go to 230 + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function d2norm( a, b ) + + real(wp) :: d2norm + real(wp), intent(in) :: a, b + + !------------------------------------------------------------------- + ! d2norm returns sqrt( a**2 + b**2 ) + ! with precautions to avoid overflow. + !------------------------------------------------------------------- + + intrinsic :: abs, sqrt + real(wp) :: scale + + d2norm = zero + scale = abs(a) + abs(b) + if (scale > zero) d2norm = scale*sqrt((a/scale)**2 + (b/scale)**2) + + end function d2norm + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine localVEnqueue(n, v, keep) + + integer(ip), intent(in) :: n + real(wp), intent(in) :: v(n) + type ( lsqr_keep ), intent(inout) :: keep + + ! Store v into the circular buffer keep%localV. + + if (keep%localPointer < keep%localVecs) then + keep%localPointer = keep%localPointer + 1 + else + keep%localPointer = 1 + keep%localVQueueFull = .true. + end if + keep%localV(1:n,keep%localPointer) = v + + end subroutine localVEnqueue + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine localVOrtho(n, v, keep, time) + + ! Perform local reorthogonalization of current v + ! against previously computed v's (which are orthonormal + ! so that ddot(n,keep%localV(1:n,j),1,keep%localV(1:n,j),1) = 1 + ! and does not need computing). + ! Modified Gram Schmidt reorthog. is used (MGS2 if keep%orthog_choice < 0) + ! Notes: Found it is more efficient to use blas ddot (rather than dot_product). + ! Using precompiled blas can lead to different performance (LSQR iteration counts) + + integer(ip), intent(in) :: n + real(wp), intent(inout) :: v(n) + type ( lsqr_keep ), intent(in) :: keep + real, intent(inout) :: time + + real(wp) :: d, ddot + integer(ip) :: j, limit + integer(long) :: time1, time2, rate_t + + ! set limit to be number of vectors we are orthogonalising against + if (keep%localVQueueFull) then + limit = keep%localVecs + else + limit = keep%localPointer + end if + + ! time the cost of RO + call system_clock(time1, rate_t) + + if (keep%orthog_tol > zero) then + ! selective reorthogonalisation + do j = 1, limit + d = ddot(n,v,1,keep%localV(1:n,j),1) + if (d .gt. keep%orthog_tol) call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + end do + else + do j = 1, limit + d = ddot(n,v,1,keep%localV(1:n,j),1) + call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + if (keep%orthog_choice .lt. 0) then + ! second application of MGS + d = ddot(n,v,1,keep%localV(1:n,j),1) + call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + end if + end do + end if + + call system_clock(time2) + + time = time + (time2 - time1)/real(rate_t) + + end subroutine localVOrtho + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + end subroutine lsqr_double + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine lsqr_free_double (keep, stat) + + ! Routine to deallocate components of keep. Problem is indicated + ! by nonzero stat value. No printing. + + type ( lsqr_keep ), intent( inout) :: keep + integer(ip), intent(out) :: stat ! set to Fortran stat parameter + ! Non zero only if a deallocation fails + integer(ip) :: st + + stat = 0 + + deallocate (keep%w, stat=st) + if (st.ne.0) stat = st + + if (allocated(keep%localV)) then + deallocate (keep%localV, stat=st) + if (st.ne.0) stat = st + end if + + ! If Papez and Tichy stopping was used, additional arrays to be dellocated + if (keep%stop_test.ne.0) then + deallocate (keep%delay, stat=st) + if (st.ne.0) stat = st + deallocate (keep%delta, stat=st) + if (st.ne.0) stat = st + deallocate (keep%curve, stat=st) + if (st.ne.0) stat = st + deallocate (keep%estim, stat=st) + if (st.ne.0) stat = st + end if + + keep%flag = 0 ! make sure error flag is reset as it is tested on the first + ! call to lsqr and we don't want to pick up an old value + + end subroutine lsqr_free_double + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine adaptive(k,add,keep) + + ! This is to get the information for the Papez and Tichy stopping + ! criteria + + integer, intent(in) :: k ! iteration number + real(wp), intent(in) :: add + + type(lsqr_keep), intent(inout) :: keep + ! components that are changed: + ! curve(1:k) + ! d + ! delay(k) + ! ell + ! estim + ! estimend + + integer :: ell + real(wp) :: den, num, s + + ell = keep%ell + + keep%curve(1:k-1) = keep%curve(1:k-1) + add + keep%curve(k) = add + + if (k > 1) then + + call find_s(k, keep%curve, keep%delta, ell, & + keep%tol, s) + + num = s * add + den = sum(keep%delta(ell:k-1)) + + do while (keep%d >= 0 .and. num <= keep%tau*den) + keep%delay(ell) = keep%d + keep%estim(ell) = sqrt(den + add) + keep%estimend = ell + ell = ell + 1 + keep%d = keep%d - 1 + den = sum(keep%delta(ell:k-1)) + end do + ! update for next call + keep%d = keep%d + 1 + keep%ell = ell + end if + keep%delta(k) = add + + end subroutine adaptive + + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + subroutine find_s(k,curve,delta,ell,tol,s) + +! This is used by subroutine adaptive + + integer, intent(in) :: k,ell + real(wp), intent(in) :: curve(k),delta(k) + real(wp), intent(in) :: tol + real(wp), intent(out) :: s +! + integer :: i,ind + real(wp) :: temp +! + ind = 1 + do i = k,1,-1 + temp = curve(ell)/curve(i) + if (temp <= tol) ind = i + end do + + s = zero + do i = ind, k-1 + s = max(s, curve(i)/delta(i)) + end do +! + end subroutine find_s + +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +end module lsqr_reverse_double diff --git a/libRALFit/src/external/preprocessor.FPP b/libRALFit/src/external/preprocessor.FPP new file mode 120000 index 00000000..852dbbdc --- /dev/null +++ b/libRALFit/src/external/preprocessor.FPP @@ -0,0 +1 @@ +libRALFit/src/preprocessor.FPP \ No newline at end of file diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index fd004599..0333dc48 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -4233,8 +4233,8 @@ Subroutine check_options(opt, inform) inform%status = NLLS_ERROR_PRINT_LEVEL ElseIf (opt%box_linesearch_type<1 .Or. opt%box_linesearch_type>2) Then inform%status = NLLS_ERROR_UNSUPPORTED_LINESEARCH -! ElseIf -! ... + ElseIf (opt%lls_solver == 2 .and. opt%sketch_size <= 0) then + inform%status = NLLS_ERROR_BAD_SKETCH_SIZE End If End Subroutine check_options diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 8f8470cd..2f79fbd4 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -7,9 +7,13 @@ #include "preprocessor.FPP" module MODULE_PREC(ral_nlls_linear) - use MODULE_PREC(ral_nlls_workspaces), only: wp, solve_LLS_work, NLLS_inform, & - NLLS_options, NLLS_ERROR_WORKSPACE_ERROR, & - NLLS_ERROR_FROM_EXTERNAL, NLLS_ERROR_BAD_LLS_SOLVER + use MODULE_PREC(ral_nlls_workspaces), only: wp, solve_LLS_work, LLS_lsqr_work, LLS_rand_work, & + NLLS_inform, NLLS_options, NLLS_ERROR_WORKSPACE_ERROR, & + NLLS_ERROR_FROM_EXTERNAL, NLLS_ERROR_BAD_LLS_SOLVER, & + NLLS_ERROR_BAD_SKETCH_METHOD, NLLS_ERROR_BAD_SKETCH_SIZE + use MODULE_PREC(lsqr_reverse), only: lsqr, lsqr_options, lsqr_inform, lsqr_keep, lsqr_free + use MODULE_PREC(dct_module), only: dct + implicit none private @@ -41,7 +45,7 @@ subroutine solve_LLS(A, b, n, m, inform, w, options, pd) type(NLLS_options), intent(in) :: options logical, intent(in) :: pd - select case (options%lls_solver) + select case (options%lls_solver) case (1) ! LAPACK if (pd) then call solve_posv(A,b,n,inform) @@ -49,18 +53,40 @@ subroutine solve_LLS(A, b, n, m, inform, w, options, pd) if (n == m) then call solve_gesv(A,b,n,inform,w) else - if (.not. w%allocated) then - inform%status = NLLS_ERROR_WORKSPACE_ERROR - goto 100 - end if call solve_gels(A,b,n,m,inform,w,options) end if end if + case (2) ! use LSQR + if (.not. w%lsqr_ws%allocated) then + inform%status = NLLS_ERROR_WORKSPACE_ERROR + return + end if + call solve_lsqr(A, b, n, m, inform, w%lsqr_ws, options) + ! with LAPACK fallback... + if (inform%status /= 0) then + if (options%allow_fallback_method) then + inform%status = 0 + call solve_gels(A,b,n,m,inform,w,options) + else + return + end if + end if + case (3) ! Randomised method (Blendenpik) + if (.not. w%lsqr_ws%allocated .or. .not. w%rand_ws%allocated) then + inform%status = NLLS_ERROR_WORKSPACE_ERROR + return + end if + call blendenpik(A, b, n, m, options%sketch_size, inform, w, options) + ! if Blendenpik fails, fall back to LAPACK case default inform%status = NLLS_ERROR_BAD_LLS_SOLVER end select + ! if LSQR or Blendenpik fails, fallback to LAPACK GELS + if (inform%status /= 0 .and. options%lls_solver > 1 .and. options%allow_fallback_method) then + inform%status = 0 + call solve_gels(A,b,n,m,inform,w,options) + end if -100 continue end subroutine solve_LLS subroutine solve_gels(A,b,n,m,inform,w,options) @@ -73,7 +99,10 @@ subroutine solve_gels(A,b,n,m,inform,w,options) integer :: lwork, lda, ldb type( solve_LLS_work ), Intent(inout) :: w - + if (.not. w%allocated) then + inform%status = NLLS_ERROR_WORKSPACE_ERROR + return + end if lwork = size(w%work) if (options%Fortran_Jacobian) then @@ -126,4 +155,261 @@ subroutine solve_posv(A,b,n,inform) end subroutine solve_posv + subroutine solve_lsqr(A, b, n, m, inform, w, options) + ! non-preconditioned LSQR iterative solver for large sparse systems + ! Implementation by Jennifer Scott, using the Papez-Tichy stopping criterion + ! This is a reverse communication implementation, so we need + ! a wrapper to handle the arithmetic ourselves + integer, intent(in) :: m, n + real(wp), intent(inout), contiguous :: A(:, :), b(:) + type(NLLS_inform) :: inform + type(LLS_lsqr_work), intent(inout) :: w + type(NLLS_options), intent(in) :: options + type(lsqr_options) :: lsqr_opt + + integer :: action + real(wp) :: norm_a + + lsqr_opt%stop_test = 2 + lsqr_opt%itnlim = 5 * n + + ! stopping criterion requires an estimate of the spectral norm of A, + ! which we can get from a few iterations of the power method + call estimate_norm(A, m, n, norm_a, 15, w) + + w%u = b + w%v = 0.0_wp + w%z = 0.0_wp + w%x = 0.0_wp + + action = 0 + + ! first call to lsqr initializes the algorithm and returns the first action + call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) + + do while (action /= 0) + select case (action) + case (1) ! compute v = v + A^T u + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 1.0_wp, w%v, 1) + case (2) ! compute u = u + A v + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 1.0_wp, w%u, 1) + end select + call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) + end do + + if (w%inform%flag /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_return = w%inform%flag + inform%external_name = 'lsqr' + return + end if + + + b(1:n) = w%x(1:n) + + end subroutine solve_lsqr + + subroutine blendenpik(A, b, n, m, sketch_size, inform, w, options) + ! Blendenpik (sketch-and-precondition) randomised solver + ! Source: + ! P. Avron, P. Maymounkov, S. Toledo, + ! "Blendenpik: Supercharging LAPACK's Least-Squares Solver" + ! URL: https://pdos.csail.mit.edu/~petar/papers/blendenpik-v1.pdf + real(wp), intent(inout), contiguous :: A(:,:), b(:) + integer, intent(in) :: sketch_size, m, n + type(NLLS_inform), INTENT(INOUT) :: inform + type(solve_LLS_work), Intent(inout) :: w + type(NLLS_options), Intent(In) :: options + type(lsqr_options) :: lsqr_opt + + integer :: i, lwork, action + real(wp) :: norm_a + + lwork = size(w%work) + lsqr_opt%stop_test = 2 + + w%rand_ws%tau = 0.0_wp + w%rand_ws%R = 0.0_wp + w%rand_ws%SM = 0.0_wp + w%rand_ws%ATu = 0.0_wp + w%rand_ws%DCT_A = 0.0_wp + + ! get sketch matrix SM + call sketch(A, m, n, sketch_size, w%rand_ws%SM, options, w%rand_ws, inform) + if (inform%status /= 0) then + return + end if + + ! now take QR decomposition of SM + call PREC(geqrf)(sketch_size, n, w%rand_ws%SM, sketch_size, w%rand_ws%tau, w%work, lwork, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?geqrf' + return + end if + + ! todo: initial guess! + ! extract R from QR decomposition + do i = 1, n + w%rand_ws%R(i, i:n) = w%rand_ws%SM(i, i:n) + end do + + call estimate_norm(A, m, n, norm_a, 15, w%lsqr_ws) + + ! first, we run LSQR without preconditioning to + w%lsqr_ws%x = 0.0_wp + + ! solve using LSQR + ! first call to lsqr initializes the algorithm and returns the first action + w%lsqr_ws%u = b + w%lsqr_ws%v = 0.0_wp + w%lsqr_ws%z = 0.0_wp + action = 0 + call lsqr(1, action, m, n, w%lsqr_ws%u, w%lsqr_ws%v, w%lsqr_ws%z, & + w%lsqr_ws%x, w%lsqr_ws%keep, lsqr_opt, w%lsqr_ws%inform, & + anorm=norm_a) + do while (action /= 0) + select case (action) + case (1) ! compute v = v + P^{-1} A^T u + ! we don't want to explicitly invert R (for numerical stability), so we do this in two steps: + ! first compute ATu = A^T u, then let t = R^{-1} ATu (i.e. solve R t = ATu) + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%lsqr_ws%u, 1, 0.0_wp, w%rand_ws%ATu, 1) + call PREC(trtrs)('U', 'N', 'N', n, 1, w%rand_ws%R, n, w%rand_ws%ATu, n, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?trtrs' + return + end if + w%lsqr_ws%v = w%lsqr_ws%v + w%rand_ws%ATu + case (2) ! compute u = u + Az + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%lsqr_ws%z, 1, 1.0_wp, w%lsqr_ws%u, 1) + case (3) ! compute z = P^{-T} v (or equivalently, solve R^T z = v) + w%lsqr_ws%z = w%lsqr_ws%v + call PREC(trtrs)('U', 'T', 'N', n, 1, w%rand_ws%R, n, w%lsqr_ws%z, n, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?trtrs' + return + end if + end select + call lsqr(1, action, m, n, w%lsqr_ws%u, w%lsqr_ws%v, w%lsqr_ws%z, & + w%lsqr_ws%x, w%lsqr_ws%keep, lsqr_opt, w%lsqr_ws%inform, & + anorm=norm_a) + end do + if (w%lsqr_ws%inform%flag /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_return = w%lsqr_ws%inform%flag + inform%external_name = 'lsqr' + return + end if + + b(1:n) = w%lsqr_ws%x(1:n) + + end subroutine blendenpik + + subroutine estimate_norm(A, m, n, norm_a, num_iters, w) + ! power method to estimate spectral norm of A + real(wp), intent(in), contiguous :: A(:,:) + integer, intent(in) :: m, n, num_iters + real(wp), intent(out) :: norm_a + type(LLS_lsqr_work), intent(inout) :: w + + integer :: i + + ! start with a random vector + call random_number(w%v) + + ! each iteration of the power method, we multiply by A^T A and renormalize + ! this is done in two steps to be better-conditioned + do i = 1, num_iters + w%v = w%v / norm2(w%v) + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 0.0_wp, w%u, 1) + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 0.0_wp, w%v, 1) + end do + + norm_a = sqrt(norm2(w%v)) + + end subroutine estimate_norm + + subroutine sketch(A, m, n, sketch_size, SA, options, w, inform) + ! Given a matrix A, return a "sketch" SA of size `sketch_size x n`. + real(wp), intent(in), contiguous :: A(:,:) + integer, intent(in) :: m, n, sketch_size + real(wp), intent(out), contiguous :: SA(:,:) + type(NLLS_options), intent(in) :: options + type(LLS_rand_work), intent(inout) :: w + type(NLLS_inform), intent(inout) :: inform + + if (sketch_size <= 0 .or. sketch_size > m) then + inform%status = NLLS_ERROR_BAD_SKETCH_SIZE + return + end if + + select case (options%sketch_method) + case (1) ! random sampling of rows + call select_random_rows(A, n, sketch_size, SA, w) + case (2) ! random projection using DCT + call dct_sketch(A, m, n, sketch_size, SA, w) + case default + inform%status = NLLS_ERROR_BAD_SKETCH_METHOD + end select + + end subroutine sketch + + subroutine select_random_rows(A, n, s, SA, w) + ! Select `s` random rows of A to form SA + real(wp), intent(in), contiguous :: A(:,:) + integer, intent(in) :: n, s + real(wp), intent(out) :: SA(:,:) + type(LLS_rand_work), intent(inout) :: w + + integer :: i, idx + + call random_number(w%temp) + + ! to select `s` random rows, we generate `m` random numbers in [0,1), + ! and take the indices of the smallest `s` numbers + do i = 1, s + idx = minloc(w%temp, dim=1) + w%temp(idx) = 1.0_wp + SA(i, :) = A(idx, :) + end do + end subroutine select_random_rows + + subroutine dct_sketch(A, m, n, s, SA, w) + ! Create a sketch of A by taking a random subsample of rows + ! applying random sign flips and a DCT to reduce coherence + ! (i.e. the probability that the sketch is rank deficient) + real(wp), intent(in), contiguous :: A(:,:) + integer, intent(in) :: m, n, s + real(wp), intent(inout) :: SA(:,:) + type(LLS_rand_work), intent(inout) :: w + + integer :: i + real(wp) :: rand_no + + w%DCT_A = 0.0_wp + + ! create DA, which is A with random sign flips + ! i.e. D is a random diagonal matrix with +/- 1 on the diagonal + do i = 1, m + call random_number(rand_no) + if (rand_no < 0.5_wp) then + w%DCT_A(i, :) = -A(i, :) + else + w%DCT_A(i, :) = A(i, :) + end if + end do + + ! apply discrete cosine transform to each column + call dct(w%DCT_A, n, m) + + ! and normalise DCT + w%DCT_A = w%DCT_A / sqrt(2.0_wp * real(m, wp)) + + ! and finaly select `s` random rows of DCT_A to form SA + call select_random_rows(w%DCT_A, n, s, SA, w) + + end subroutine dct_sketch + end module MODULE_PREC(ral_nlls_linear) diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index f5c1e949..190ded6f 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -12,6 +12,7 @@ module MODULE_PREC(ral_nlls_workspaces) Use MODULE_PREC(ral_nlls_types), Only: wp, np, lp, params_base_type, & eval_f_type, eval_j_type, & eval_hf_type, eval_hp_type + use MODULE_PREC(lsqr_reverse), Only: lsqr_options, lsqr_inform, lsqr_keep, lsqr_free implicit none @@ -49,6 +50,8 @@ module MODULE_PREC(ral_nlls_workspaces) Integer, Parameter, Public :: NLLS_ERROR_BAD_JACOBIAN = -19 Integer, Parameter, Public :: NLLS_ERROR_BAD_WEIGHTS = -20 Integer, Parameter, Public :: NLLS_ERROR_BAD_LLS_SOLVER = -21 + Integer, Parameter, Public :: NLLS_ERROR_BAD_SKETCH_SIZE = -22 + Integer, Parameter, Public :: NLLS_ERROR_BAD_SKETCH_METHOD = -23 ! dogleg errors Integer, Parameter, Public :: NLLS_ERROR_DOGLEG_MODEL = -101 @@ -144,11 +147,6 @@ module MODULE_PREC(ral_nlls_workspaces) LOGICAL :: allow_fallback_method = .true. - -! which linear least squares solver should we use? - - INTEGER :: lls_solver = 1 - ! overall convergence tolerances. The iteration will terminate when the ! norm of the gradient of the objective function is smaller than ! MAX( %stop_g_absolute, %stop_g_relative * norm of the initial gradient) @@ -286,6 +284,29 @@ module MODULE_PREC(ral_nlls_workspaces) logical :: setup_workspaces = .true. logical :: remove_workspaces = .true. +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +!!! L I N E A R S O L V E R C O N T R O L S !!! +!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + +! which linear least squares solver should we use? +! 1 LAPACK solver +! 2 LSQR (iterative) +! 3 Randomised solver (Blendenpik) + + INTEGER :: lls_solver = 1 + +! If using randomised linear solver (Blendenpik), what sketch size? +! this must be manually set by the user as it depends on `m` + + INTEGER :: sketch_size = 0 + +! If using randomised solver, which sketching method should we use? +! 1 Uniform subsample +! 2 Subsample with DCT row mixing +! + INTEGER :: sketch_method = 1 + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! M O R E - S O R E N S E N C O N T R O L S !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -599,10 +620,28 @@ module MODULE_PREC(ral_nlls_workspaces) real(wp), allocatable :: Jd(:), dH(:), Hd(:), dHd(:) end type evaluate_model_work + type, public :: LLS_lsqr_work + logical :: allocated = .false. + ! variables for LSQR + real(wp), allocatable :: u(:), v(:), z(:), x(:) + + ! LSQR reverse communication objects + type(lsqr_keep) :: keep + type(lsqr_inform) :: inform + end type LLS_lsqr_work + + type, public :: LLS_rand_work + logical :: allocated = .false. + real(wp), allocatable :: tau(:), temp(:), ATu(:) + real(wp), allocatable :: R(:,:), SM(:,:), DCT_A(:,:) + end type LLS_rand_work + type, public :: solve_LLS_work ! workspace for subroutine solve_LLS logical :: allocated = .false. real(wp), allocatable :: temp(:), work(:), Jlls(:,:) integer, allocatable :: ipiv(:) + type(LLS_lsqr_work) :: lsqr_ws + type(LLS_rand_work) :: rand_ws end type solve_LLS_work type, public :: min_eig_symm_work ! workspace for subroutine min_eig_work @@ -731,6 +770,7 @@ module MODULE_PREC(ral_nlls_workspaces) Public :: params_base_type Public :: eval_f_type, eval_j_type, eval_hf_type, eval_hp_type public :: setup_workspaces, remove_workspaces + public :: setup_workspace_solve_LLS, remove_workspace_solve_LLS public :: setup_workspace_dogleg, setup_workspace_AINT_tr public :: setup_workspace_more_sorensen, setup_workspace_solve_galahad public :: setup_workspace_regularization_solver @@ -1256,6 +1296,12 @@ subroutine setup_workspace_solve_LLS(n,m,w,options,inform) inform%status = 0 lwork = max(1, min(m,n) + max(min(m,n), 1)*4) allocate( w%temp(max(m,n)),w%work(lwork),w%Jlls(m,n), w%ipiv(n), stat = inform%alloc_status) + if (options%LLS_solver /= 1) then + call setup_workspace_LLS_lsqr(w%lsqr_ws, n, m, inform) + end if + if (options%LLS_solver == 3) then + call setup_workspace_LLS_rand(w%rand_ws, n, m, options, inform) + end if If (inform%alloc_status /= 0) Then Call remove_workspace_solve_LLS(w,options) inform%status = NLLS_ERROR_ALLOCATION @@ -1275,10 +1321,76 @@ subroutine remove_workspace_solve_LLS(w,options) if(allocated( w%work )) deallocate( w%work, stat=ierr_dummy ) if(allocated( w%Jlls )) deallocate( w%Jlls, stat=ierr_dummy ) if(allocated( w%ipiv )) deallocate( w%ipiv, stat=ierr_dummy ) + call remove_workspace_LLS_lsqr(w%lsqr_ws) + call remove_workspace_LLS_rand(w%rand_ws) w%allocated = .false. end subroutine remove_workspace_solve_LLS + subroutine setup_workspace_LLS_lsqr(w, n, m, inform) + implicit none + type( LLS_lsqr_work ), intent(out) :: w + integer, intent(in) :: n, m + type( nlls_inform ), intent(inout) :: inform + + inform%status = 0 + allocate(w%u(m), w%v(n), w%z(n), w%x(n), stat=inform%alloc_status) + If (inform%alloc_status /= 0) Then + inform%bad_alloc = "setup_workspace_LLS_lsqr" + inform%status = NLLS_ERROR_ALLOCATION + return + End If + + w%allocated = .true. + end subroutine setup_workspace_LLS_lsqr + + subroutine remove_workspace_LLS_lsqr(w) + implicit none + type( LLS_lsqr_work ), intent(out) :: w + Integer :: ierr_dummy + + if(allocated(w%u)) deallocate(w%u, stat=ierr_dummy) + if(allocated(w%v)) deallocate(w%v, stat=ierr_dummy) + if(allocated(w%z)) deallocate(w%z, stat=ierr_dummy) + if(allocated(w%x)) deallocate(w%x, stat=ierr_dummy) + + w%allocated = .false. + end subroutine remove_workspace_LLS_lsqr + + subroutine setup_workspace_LLS_rand(w, n, m, options, inform) + implicit none + type( LLS_rand_work ), intent(out) :: w + integer, intent(in) :: n, m + type( nlls_options ), intent(in) :: options + type( nlls_inform ), intent(inout) :: inform + + inform%status = 0 + allocate(w%tau(min(options%sketch_size, n)), w%temp(m), w%R(n,n), & + w%SM(options%sketch_size, n), w%ATu(n), w%DCT_A(m,n), stat=inform%alloc_status) + If (inform%alloc_status /= 0) Then + inform%bad_alloc = "setup_workspace_LLS_rand" + inform%status = NLLS_ERROR_ALLOCATION + return + End If + + w%allocated = .true. + end subroutine setup_workspace_LLS_rand + + subroutine remove_workspace_LLS_rand(w) + implicit none + type( LLS_rand_work ), intent(out) :: w + Integer :: ierr_dummy + + if(allocated(w%tau)) deallocate(w%tau, stat=ierr_dummy) + if(allocated(w%temp)) deallocate(w%temp, stat=ierr_dummy) + if(allocated(w%R)) deallocate(w%R, stat=ierr_dummy) + if(allocated(w%SM)) deallocate(w%SM, stat=ierr_dummy) + if(allocated(w%ATu)) deallocate(w%ATu, stat=ierr_dummy) + if(allocated(w%DCT_A)) deallocate(w%DCT_A, stat=ierr_dummy) + + w%allocated = .false. + end subroutine remove_workspace_LLS_rand + subroutine setup_workspace_evaluate_model(n,m,w,options,inform) implicit none integer, intent(in) :: n, m diff --git a/libRALFit/test/nlls_test.F90 b/libRALFit/test/nlls_test.F90 index e6c0658d..ffbf9f77 100644 --- a/libRALFit/test/nlls_test.F90 +++ b/libRALFit/test/nlls_test.F90 @@ -1436,6 +1436,12 @@ program nlls_test call solve_LLS_dgesv_tests(options,fails) no_errors_helpers = no_errors_helpers + fails + call solve_LLS_lsqr_tests(options,fails) + no_errors_helpers = no_errors_helpers + fails + + call solve_LLS_randomised_tests(options,fails) + no_errors_helpers = no_errors_helpers + fails + call findbeta_tests(options,fails) no_errors_helpers = no_errors_helpers + fails diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index f7414e8b..38e7118f 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -2427,6 +2427,175 @@ subroutine solve_LLS_dgesv_tests(options,fails) end subroutine solve_LLS_dgesv_tests + subroutine solve_LLS_lsqr_tests(options, fails) + type( nlls_options ), intent(inout) :: options + integer, intent(out) :: fails + + real(wp), allocatable :: J(:,:), J_copy(:,:), Jd(:), f(:), d(:), x_calc(:), x_true(:) + real(wp) :: normerror + integer :: n, m, i, k, l + type( nlls_workspace ) :: work + type( nlls_workspace ), Target :: iw + type( nlls_inform ) :: inform + + fails = 0 + work%iw_ptr => iw + iw%iw_ptr => iw + + n = 5 + m = 30 + + fails = 0 + + ! setup without LSQR option, so LSQR workspace isn't initialised + options%lls_solver = 1 + allocate(J(m,n), J_copy(m,n), f(m), d(m), Jd(m), x_calc(n), x_true(n)) + call setup_workspaces(work,n,m,options,inform) + + ! this is mostly a generic version of the matrices in the other tests; + ! i.e. J = [ 1, 2, 3, 4, 5...] reshaped to (n, m) + ! but in this case we need a bigger matrix to ensure we get a decent sketch + ! and we add a multiplicative factor to make the matrix full-rank + do k = 1, m*n + i = ((k-1) / n) + 1 + l = mod(k-1, n) + 1 + J(i, l) = real(k + i*l, wp) + end do + + do k = 1, m + f(k) = 7.0_wp + real(k, wp) * 2.0_wp + end do + + J_copy = J + d = f + options%nlls_method = 1 ! dogleg (just so we know what workspace we're in) + options%lls_solver = 2 ! LSQR + options%allow_fallback_method = .false. + + call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + if ( inform%status .ne. NLLS_ERROR_WORKSPACE_ERROR ) then + write (*, *) 'solve_LLS LSQR test failed: wrong error message returned when workspace unallocated' + write(*,*) 'status = ', inform%status, " (expected ",NLLS_ERROR_WORKSPACE_ERROR,")" + fails = fails + 1 + end if + + call remove_workspaces(work,options) + call setup_workspaces(work,n,m,options,inform) + + + ! test correct residual + call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + if ( inform%status .ne. 0 ) then + write (*, *) 'solve_LLS LSQR test failed: error message returned' + write(*,*) 'status = ', inform%status, " (expected 0)" + if (inform%status == NLLS_ERROR_FROM_EXTERNAL) then + write(*,*) 'external return = ', inform%external_return + write(*,*) 'external name = ', inform%external_name + end if + fails = fails + 1 + end if + + ! check answer + call mult_J(J,n,m,d,Jd,.True.) + normerror = norm2(Jd - f) + if ( .not. normerror < 1.0e-12_wp ) then + ! wrong answer, as data chosen to fit + write(*,*) 'solve_LLS LSQR test failed: wrong solution returned' + write(*,*) '||Jd - f|| = ', normerror + fails = fails + 1 + end if + + call nlls_finalize(work,options,inform) + call reset_default_options(options) + end subroutine solve_LLS_lsqr_tests + + subroutine solve_LLS_randomised_tests(options,fails) + type( nlls_options ), intent(inout) :: options + integer, intent(out) :: fails + + real(wp), allocatable :: J(:,:), J_copy(:,:), Jd(:), d(:), f(:) + real(wp) :: normerror + integer :: n,m,k,i,l,sketch_method + type( nlls_workspace ) :: work + type( nlls_workspace ), Target :: iw + type( nlls_inform ) :: inform + + fails = 0 + work%iw_ptr => iw + iw%iw_ptr => iw + + n = 5 + m = 60 + + options%nlls_method = 1 ! dogleg (just so we know what workspace we're in) + options%lls_solver = 3 ! randomised solver + options%allow_fallback_method = .false. + + allocate(J(m,n), J_copy(m,n), Jd(m), d(m), f(m)) + + call setup_workspaces(work,n,m,options,inform) + + ! ensure solver gives error if sketch size is not a positive integer + options%sketch_size = 0 + call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + if ( inform%status .ne. NLLS_ERROR_BAD_SKETCH_SIZE ) then + write (*, *) 'solve_LLS randomised test failed: wrong error message returned' + write(*,*) 'status = ', inform%status, " (expected ",NLLS_ERROR_BAD_SKETCH_SIZE,")" + fails = fails + 1 + end if + inform%status = 0 + + call remove_workspaces(work,options) + + options%sketch_size = 10 ! reset to valid sketch size + + ! this is mostly a generic version of the matrices in the other tests; + ! i.e. J = [ 1, 2, 3, 4, 5...] reshaped to (n, m) + ! but in this case we need a bigger matrix to ensure we get a decent sketch + ! and we add a multiplicative factor to make the matrix full-rank + do k = 1, m*n + i = ((k-1) / n) + 1 + l = mod(k-1, n) + 1 + J(i, l) = real(k + i*l, wp) + end do + + do k = 1, m + f(k) = 7.0_wp + real(k, wp) * 2.0_wp + end do + + do sketch_method = 1, 2 + options%sketch_method = sketch_method + call setup_workspaces(work,n,m,options,inform) + + J_copy = J + d = f + call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + if ( inform%status .ne. 0 ) then + write (*, *) 'solve_LLS randomised test failed: error message returned' + write(*,*) 'status = ', inform%status, " (expected 0)" + write(*,*) 'sketch method = ', sketch_method + if (inform%status == NLLS_ERROR_FROM_EXTERNAL) then + write(*,*) 'external return = ', inform%external_return + write(*,*) 'external name = ', inform%external_name + end if + fails = fails + 1 + end if + call mult_J(J,n,m,d,Jd,.True.) + normerror = norm2(Jd - f) + if ( .not. normerror < 1.0e-11_wp ) then + ! wrong answer, as data chosen to fit + write(*,*) 'solve_LLS randomised test failed: wrong solution returned' + write(*,*) '||Jd - f|| = ', normerror + write(*,*) 'sketch method = ', sketch_method + fails = fails + 1 + end if + call remove_workspaces(work,options) + end do + call nlls_finalize(work,options,inform) + call reset_default_options(options) + + end subroutine solve_LLS_randomised_tests + subroutine findbeta_tests(options,fails) type( nlls_options ), intent(inout) :: options From 0fced069b181c3ce91ad978e4c79d1b0c0000204 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 9 Mar 2026 07:34:46 +0000 Subject: [PATCH 12/16] improvements to external routines --- libRALFit/src/external/CMakeLists.txt | 4 +- libRALFit/src/external/dct.f90 | 74 ++++++++++++++++++++++--- libRALFit/src/external/lsqr.f90 | 50 +++++++++-------- libRALFit/src/external/preprocessor.FPP | 25 ++++++++- libRALFit/src/ral_nlls_linear.F90 | 29 +++++++--- libRALFit/src/ral_nlls_workspaces.F90 | 5 +- 6 files changed, 144 insertions(+), 43 deletions(-) mode change 120000 => 100644 libRALFit/src/external/preprocessor.FPP diff --git a/libRALFit/src/external/CMakeLists.txt b/libRALFit/src/external/CMakeLists.txt index 348c3b77..0a2ece1a 100644 --- a/libRALFit/src/external/CMakeLists.txt +++ b/libRALFit/src/external/CMakeLists.txt @@ -6,4 +6,6 @@ endif() add_library(lsqr_reverse_${precision} OBJECT lsqr.f90) add_library(dct_module_${precision} OBJECT dct.f90) -target_link_libraries(dct_module_${precision} fftw3) +# NOTE this currently requires the new dev build of FFTZ +# i.e. version 5.3 or above +target_link_libraries(dct_module_${precision} fftw3xc_wrapper_fftz) diff --git a/libRALFit/src/external/dct.f90 b/libRALFit/src/external/dct.f90 index 35412a03..c24a9eca 100644 --- a/libRALFit/src/external/dct.f90 +++ b/libRALFit/src/external/dct.f90 @@ -2,34 +2,90 @@ #include "../preprocessor.FPP" module MODULE_PREC(dct_module) + implicit none +#ifdef SINGLE_PRECISION + integer, parameter :: wp = selected_real_kind(6) ! c_float +#else + integer, parameter :: wp = selected_real_kind(15) ! c_double +#endif + ! these parameters are defined in FFTW, but we define them here to avoid including the FFTW header - integer(8) :: FFTW_REDFT10 - parameter (FFTW_REDFT10 = 4) - integer(8) :: FFTW_ESTIMATE + integer :: FFTW_FORWARD + parameter (FFTW_FORWARD = -1) + integer :: FFTW_ESTIMATE parameter (FFTW_ESTIMATE = 64) - integer(4), parameter :: wp = kind( 0.0d+0 ) interface dct module procedure dct_fftz end interface dct + + !interface fftw_plan_dft_1d + ! module subroutine fftw_plan_dft_1d(plan, n, in, out, sign, flags) bind(C) + ! use iso_c_binding, only: c_int, c_float, c_double + ! integer(c_int) :: plan + ! integer(c_int), value :: n + ! complex(wp) :: in(*) + ! complex(wp) :: out(*) + ! integer(c_int), value :: sign + ! integer(c_int), value :: flags + ! end subroutine fftw_plan_dft_1d + !end interface fftw_plan_dft_1d contains - subroutine dct_fftz(A, n, m) + subroutine dct_fftz(A, n, m, work) ! Compute the 1D DCT across the columns of A using the FFTZ library. - real(8), intent(inout), contiguous :: A(:,:) + ! `work` is expected to be a complex array of size m. + real(wp), intent(inout), contiguous :: A(:,:) + complex(wp), intent(inout) :: work(:) integer, intent(in) :: n, m integer :: i, plan - call dfftw_plan_r2r_1d(plan, m, A(:, 1), A(:, 1), FFTW_REDFT10, FFTW_ESTIMATE) + plan = 0 + do i = 1, n - call dfftw_execute_r2r(plan, A(:, i), A(:, i)) + call dct_fftz_1d(A(:, i), m, work) end do - call dfftw_destroy_plan(plan) end subroutine dct_fftz + subroutine dct_fftz_1d(x, m, work) + ! Compute the 1D DCT of x using the FFTZ library. + ! This uses an FFT and a phase factor to compute the DCT coefficients, + ! because it's easier to get FFT codes than DCT ones. + ! `work` is expected to be a complex array of size m, + ! and `plan` is an integer that can be reused across calls + ! to avoid the overhead of creating a new plan each time. + real(wp), intent(inout), contiguous :: x(:) + complex(wp), intent(inout) :: work(:) + complex(wp) :: phase_factor + integer, intent(in) :: m + + integer :: i, plan + + ! phase factor for converting FFT output into DCT coefficients + phase_factor = complex(0.0_wp, -3.14159_wp/(2.0_wp*m)) + + ! The DCT can be computed using an FFT by creating a new array where the first m/2 elements are the even-indexed + ! elements of x and the next m/2 elements are the odd-indexed elements of x (in reverse). + ! Then we can compute the DCT using an FFT on this new array. + do i = 1, m/2 + work(i) = complex(x(2*i - 1), 0.0_wp) + work(m/2 + i) = complex(x(2*(m/2 - i + 1)), 0.0_wp) + end do + + call dfftw_plan_dft_1d(plan, m, work, work, FFTW_FORWARD, FFTW_ESTIMATE) + call dfftw_execute_dft(plan, work, work) + + do i = 1, m + x(i) = real(work(i) * exp(phase_factor * (i-1)), wp) + end do + + call dfftw_destroy_plan(plan) + + end subroutine dct_fftz_1d + end module MODULE_PREC(dct_module) diff --git a/libRALFit/src/external/lsqr.f90 b/libRALFit/src/external/lsqr.f90 index 0406b50c..ee1bae11 100644 --- a/libRALFit/src/external/lsqr.f90 +++ b/libRALFit/src/external/lsqr.f90 @@ -68,17 +68,23 @@ ! This requires that x0 be available before the first call ! to LSQR and after the final call. !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +#include "preprocessor.FPP" module lsqr_reverse_double implicit none +#ifdef SINGLE_PRECISION + integer, parameter :: wp = selected_real_kind(6) ! c_float +#else + integer, parameter :: wp = selected_real_kind(15) ! c_double +#endif + private public :: lsqr_keep, lsqr_options, lsqr_inform public :: lsqr, lsqr_free integer(4), parameter :: ip = kind( 0 ) - integer(4), parameter :: wp = kind( 0.0d+0 ) integer, parameter :: long = selected_int_kind(18) real(wp), parameter :: zero = 0.0_wp, one = 1.0_wp @@ -615,7 +621,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & keep%beta = dnrm2 (m, u, 1) - call dscal (m, (one/keep%beta), u, 1) ! u_1 = u / beta + call PREC(scal)(m, (one/keep%beta), u, 1) ! u_1 = u / beta ! User must compute v = A^T * u (precon = 0 or 2), ! or v = P^{-1} A^T *u (precon = 1) @@ -644,7 +650,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & return end if - call dscal (n, (one/keep%alpha), v, 1) + call PREC(scal)(n, (one/keep%alpha), v, 1) ! w_1 = v_1 keep%w = v @@ -675,7 +681,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & ! beta*u = A*v - alpha*u, ! alpha*v = (A)'*u - beta*v. !---------------------------------------------------------------- - call dscal (m,(- keep%alpha), u, 1) ! - alpha*u_j + call PREC(scal)(m,(- keep%alpha), u, 1) ! - alpha*u_j ! Compute u = u + A*v action = 2 @@ -688,7 +694,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & keep%beta = dnrm2 (m, u, 1) if (keep%beta > zero) then - call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + call PREC(scal)(m, (one/keep%beta), u, 1) ! u = u / beta if (keep%localVecs > 0) then ! Check whether to store v_j for local reorthog'n of v_{j+1} if (abs(keep%orthog_choice).eq.1 .or. & @@ -696,7 +702,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & call localVEnqueue(n, v, keep) end if - call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + call PREC(scal)(n, (-keep%beta), v, 1) ! v = -beta *v_j ! Compute v = v + A'*u action = 1 @@ -712,7 +718,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & call localVOrtho(n, v, keep, inform%time_RO) keep%alpha = dnrm2 (n, v, 1) - if (keep%alpha > zero) call dscal (n, (one/keep%alpha), v, 1) + if (keep%alpha > zero) call PREC(scal)(n, (one/keep%alpha), v, 1) end if @@ -844,7 +850,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & return end if - call dscal (n, (one/keep%alpha), v, 1) + call PREC(scal)(n, (one/keep%alpha), v, 1) !---------------------------------------------------------------- ! user to compute z_1 = P^{-T} v_1 @@ -885,7 +891,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & ! beta*u = A*z - alpha*u, ! alpha*v = P^{-1} A^T*u - beta*v. !---------------------------------------------------------------- - call dscal (m,(-keep%alpha), u, 1) ! -alpha*u_j + call PREC(scal)(m,(-keep%alpha), u, 1) ! -alpha*u_j ! Compute u = u + A*z action = 2 @@ -898,7 +904,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & keep%beta = dnrm2 (m, u, 1) if (keep%beta > zero) then - call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + call PREC(scal)(m, (one/keep%beta), u, 1) ! u = u / beta if (keep%localVecs > 0) then ! Check whether to store v_j for local reorthog'n of v_{j+1} if (abs(keep%orthog_choice).eq.1 .or. & @@ -906,7 +912,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & call localVEnqueue(n, v, keep) end if - call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + call PREC(scal)(n, (-keep%beta), v, 1) ! v = -beta *v_j ! Compute v = v + P^{-1} A^T*u action = 1 @@ -923,7 +929,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & call localVOrtho(n, v, keep, inform%time_RO) keep%alpha = dnrm2 (n, v, 1) - if (keep%alpha > zero) call dscal (n, (one/keep%alpha), v, 1) + if (keep%alpha > zero) call PREC(scal)(n, (one/keep%alpha), v, 1) end if @@ -1074,8 +1080,8 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & return end if - call dscal (n, (one/keep%alpha), v, 1) - call dscal (n, (one/keep%alpha), z, 1) + call PREC(scal)(n, (one/keep%alpha), v, 1) + call PREC(scal)(n, (one/keep%alpha), z, 1) ! d_1 = z_1 keep%w = z @@ -1106,7 +1112,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & ! beta*u = A*z - alpha*u, ! alpha*v = (A)'*u - beta*v. !---------------------------------------------------------------- - call dscal (m,(- keep%alpha), u, 1) ! - alpha*u_j + call PREC(scal)(m,(- keep%alpha), u, 1) ! - alpha*u_j ! Compute u = u + A*z action = 2 keep%branch = 2 @@ -1118,7 +1124,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & keep%beta = dnrm2 (m, u, 1) if (keep%beta > zero) then - call dscal (m, (one/keep%beta), u, 1) ! u = u / beta + call PREC(scal)(m, (one/keep%beta), u, 1) ! u = u / beta if (keep%localVecs > 0) then ! Check whether to store v_j for local reorthog'n of v_{j+1} if (abs(keep%orthog_choice).eq.1 .or. & @@ -1126,7 +1132,7 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & call localVEnqueue(n, v, keep) end if - call dscal (n, (-keep%beta), v, 1) ! v = -beta *v_j + call PREC(scal)(n, (-keep%beta), v, 1) ! v = -beta *v_j ! Compute v = v + (A)'*u action = 1 @@ -1154,8 +1160,8 @@ subroutine lsqr_double (precon, action, m, n, u, v, z, x, keep, & keep%alpha = sqrt( ddot(n, v, 1, z, 1)) if (keep%alpha > zero) then - call dscal (n, (one/keep%alpha), v, 1) - call dscal (n, (one/keep%alpha), z, 1) + call PREC(scal)(n, (one/keep%alpha), v, 1) + call PREC(scal)(n, (one/keep%alpha), z, 1) end if end if @@ -1350,16 +1356,16 @@ subroutine localVOrtho(n, v, keep, time) ! selective reorthogonalisation do j = 1, limit d = ddot(n,v,1,keep%localV(1:n,j),1) - if (d .gt. keep%orthog_tol) call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + if (d .gt. keep%orthog_tol) call PREC(axpy)(n,-d,keep%localV(1:n,j),1,v,1) end do else do j = 1, limit d = ddot(n,v,1,keep%localV(1:n,j),1) - call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + call PREC(axpy)(n,-d,keep%localV(1:n,j),1,v,1) if (keep%orthog_choice .lt. 0) then ! second application of MGS d = ddot(n,v,1,keep%localV(1:n,j),1) - call daxpy(n,-d,keep%localV(1:n,j),1,v,1) + call PREC(axpy)(n,-d,keep%localV(1:n,j),1,v,1) end if end do end if diff --git a/libRALFit/src/external/preprocessor.FPP b/libRALFit/src/external/preprocessor.FPP deleted file mode 120000 index 852dbbdc..00000000 --- a/libRALFit/src/external/preprocessor.FPP +++ /dev/null @@ -1 +0,0 @@ -libRALFit/src/preprocessor.FPP \ No newline at end of file diff --git a/libRALFit/src/external/preprocessor.FPP b/libRALFit/src/external/preprocessor.FPP new file mode 100644 index 00000000..f1199554 --- /dev/null +++ b/libRALFit/src/external/preprocessor.FPP @@ -0,0 +1,24 @@ +! Copyright (C) 2024-2026 Advanced Micro Devices, Inc. All rights reserved. + +#define CONCAT_(a) a +#define CONCAT(a, b) CONCAT_(a)CONCAT_(b) + +#if SINGLE_PRECISION + +! BLAS nomenclature +#define PREC(fn) CONCAT(s,fn) +! Module precision nomenclature +#define MODULE_PREC(fn) CONCAT(fn,_single) +! CIFACE nomenclature +#define IFACE_PREC(fn) CONCAT(fn,_s) + +#else + + ! BLAS nomenclature +#define PREC(fn) CONCAT(d,fn) +! Module precision nomenclature +#define MODULE_PREC(fn) CONCAT(fn,_double) +! CIFACE nomenclature +#define IFACE_PREC(fn) CONCAT(fn,_d) + +#endif \ No newline at end of file diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 2f79fbd4..0b94f3d8 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -45,16 +45,19 @@ subroutine solve_LLS(A, b, n, m, inform, w, options, pd) type(NLLS_options), intent(in) :: options logical, intent(in) :: pd + ! if A is positive definite, we need to use Cholesky solver + ! as routines rely on its test for positive-definiteness + if (pd) then + call solve_posv(A, b, n, inform) + return + end if + select case (options%lls_solver) case (1) ! LAPACK - if (pd) then - call solve_posv(A,b,n,inform) + if (n == m) then + call solve_gesv(A,b,n,inform,w) else - if (n == m) then - call solve_gesv(A,b,n,inform,w) - else - call solve_gels(A,b,n,m,inform,w,options) - end if + call solve_gels(A,b,n,m,inform,w,options) end if case (2) ! use LSQR if (.not. w%lsqr_ws%allocated) then @@ -181,6 +184,14 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options) w%v = 0.0_wp w%z = 0.0_wp w%x = 0.0_wp + w%ATu = 0.0_wp + + if (present(init_guess)) then + if (init_guess) then + ! update `u` for LSQR to be b - A x_0, so that we can start LSQR with this initial guess + call PREC(gemv)('N', m, n, -1.0_wp, A, m, w%x, 1, 1.0_wp, w%u, 1) + end if + end if action = 0 @@ -256,7 +267,6 @@ subroutine blendenpik(A, b, n, m, sketch_size, inform, w, options) call estimate_norm(A, m, n, norm_a, 15, w%lsqr_ws) - ! first, we run LSQR without preconditioning to w%lsqr_ws%x = 0.0_wp ! solve using LSQR @@ -402,7 +412,8 @@ subroutine dct_sketch(A, m, n, s, SA, w) end do ! apply discrete cosine transform to each column - call dct(w%DCT_A, n, m) + call dct(w%DCT_A, n, m, w%dct_work) + call dct1d(w%temp_1, m, w%dct_work) ! and normalise DCT w%DCT_A = w%DCT_A / sqrt(2.0_wp * real(m, wp)) diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index 190ded6f..f3202bbe 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -634,6 +634,7 @@ module MODULE_PREC(ral_nlls_workspaces) logical :: allocated = .false. real(wp), allocatable :: tau(:), temp(:), ATu(:) real(wp), allocatable :: R(:,:), SM(:,:), DCT_A(:,:) + complex(wp), allocatable :: dct_work(:) end type LLS_rand_work type, public :: solve_LLS_work ! workspace for subroutine solve_LLS @@ -1366,7 +1367,8 @@ subroutine setup_workspace_LLS_rand(w, n, m, options, inform) inform%status = 0 allocate(w%tau(min(options%sketch_size, n)), w%temp(m), w%R(n,n), & - w%SM(options%sketch_size, n), w%ATu(n), w%DCT_A(m,n), stat=inform%alloc_status) + w%SM(options%sketch_size, n), w%ATu(n), w%DCT_A(m,n), & + w%dct_work(m), stat=inform%alloc_status) If (inform%alloc_status /= 0) Then inform%bad_alloc = "setup_workspace_LLS_rand" inform%status = NLLS_ERROR_ALLOCATION @@ -1387,6 +1389,7 @@ subroutine remove_workspace_LLS_rand(w) if(allocated(w%SM)) deallocate(w%SM, stat=ierr_dummy) if(allocated(w%ATu)) deallocate(w%ATu, stat=ierr_dummy) if(allocated(w%DCT_A)) deallocate(w%DCT_A, stat=ierr_dummy) + if(allocated(w%dct_work)) deallocate(w%dct_work, stat=ierr_dummy) w%allocated = .false. end subroutine remove_workspace_LLS_rand From 756154dd8b384c664512b9da9ee537eb2c06bec0 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Mon, 9 Mar 2026 14:38:21 +0000 Subject: [PATCH 13/16] use sketch-and-solve initial guess --- libRALFit/src/external/dct.f90 | 4 + libRALFit/src/ral_nlls_linear.F90 | 192 ++++++++++++++------------ libRALFit/src/ral_nlls_workspaces.F90 | 16 +-- 3 files changed, 119 insertions(+), 93 deletions(-) diff --git a/libRALFit/src/external/dct.f90 b/libRALFit/src/external/dct.f90 index c24a9eca..09d87af1 100644 --- a/libRALFit/src/external/dct.f90 +++ b/libRALFit/src/external/dct.f90 @@ -21,6 +21,10 @@ module MODULE_PREC(dct_module) module procedure dct_fftz end interface dct + interface dct1d + module procedure dct_fftz_1d + end interface dct1d + !interface fftw_plan_dft_1d ! module subroutine fftw_plan_dft_1d(plan, n, in, out, sign, flags) bind(C) ! use iso_c_binding, only: c_int, c_float, c_double diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 0b94f3d8..6e14e2d3 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -12,7 +12,7 @@ module MODULE_PREC(ral_nlls_linear) NLLS_ERROR_FROM_EXTERNAL, NLLS_ERROR_BAD_LLS_SOLVER, & NLLS_ERROR_BAD_SKETCH_METHOD, NLLS_ERROR_BAD_SKETCH_SIZE use MODULE_PREC(lsqr_reverse), only: lsqr, lsqr_options, lsqr_inform, lsqr_keep, lsqr_free - use MODULE_PREC(dct_module), only: dct + use MODULE_PREC(dct_module), only: dct, dct1d implicit none @@ -158,20 +158,26 @@ subroutine solve_posv(A,b,n,inform) end subroutine solve_posv - subroutine solve_lsqr(A, b, n, m, inform, w, options) + subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) ! non-preconditioned LSQR iterative solver for large sparse systems ! Implementation by Jennifer Scott, using the Papez-Tichy stopping criterion ! This is a reverse communication implementation, so we need - ! a wrapper to handle the arithmetic ourselves + ! a wrapper to handle the matrix products ourselves + ! Optional arguments: + ! precon: if present, a preconditioner to apply in the matrix products + ! init_guess: if true, assume w%x has been initialised with an initial guess + ! currently, we assume the preconditioner is upper-triangular of size n x n integer, intent(in) :: m, n real(wp), intent(inout), contiguous :: A(:, :), b(:) - type(NLLS_inform) :: inform + type(NLLS_inform), intent(inout) :: inform type(LLS_lsqr_work), intent(inout) :: w type(NLLS_options), intent(in) :: options - type(lsqr_options) :: lsqr_opt + real(wp), intent(in), optional, contiguous :: precon(:,:) + logical, intent(in), optional :: init_guess - integer :: action + integer :: action, use_precon real(wp) :: norm_a + type(lsqr_options) :: lsqr_opt lsqr_opt%stop_test = 2 lsqr_opt%itnlim = 5 * n @@ -183,30 +189,67 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options) w%u = b w%v = 0.0_wp w%z = 0.0_wp - w%x = 0.0_wp w%ATu = 0.0_wp if (present(init_guess)) then if (init_guess) then ! update `u` for LSQR to be b - A x_0, so that we can start LSQR with this initial guess call PREC(gemv)('N', m, n, -1.0_wp, A, m, w%x, 1, 1.0_wp, w%u, 1) + else + w%x = 0.0_wp end if + else + w%x = 0.0_wp end if action = 0 + if (present(precon)) then + use_precon = 1 + else + use_precon = 0 + end if ! first call to lsqr initializes the algorithm and returns the first action - call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) - - do while (action /= 0) - select case (action) - case (1) ! compute v = v + A^T u - call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 1.0_wp, w%v, 1) - case (2) ! compute u = u + A v - call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 1.0_wp, w%u, 1) + call lsqr(use_precon, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) + + if (use_precon == 0) then + do while (action /= 0) + select case (action) + case (1) ! compute v = v + A^T u + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 1.0_wp, w%v, 1) + case (2) ! compute u = u + A v + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 1.0_wp, w%u, 1) + end select + call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) + end do + else + do while (action /= 0) + select case (action) + case (1) ! compute v = v + P^{-1} A^T u + ! we don't want to explicitly invert R (for numerical stability), so we do this in two steps: + ! first compute ATu = A^T u, then let t = R^{-1} ATu (i.e. solve R t = ATu) + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 0.0_wp, w%ATu, 1) + call PREC(trtrs)('U', 'N', 'N', n, 1, precon, n, w%ATu, n, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?trtrs' + return + end if + w%v = w%v + w%ATu + case (2) ! compute u = u + Az + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%z, 1, 1.0_wp, w%u, 1) + case (3) ! compute z = P^{-T} v (or equivalently, solve R^T z = v) + w%z = w%v + call PREC(trtrs)('U', 'T', 'N', n, 1, precon, n, w%z, n, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?trtrs' + return + end if end select - call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) - end do + call lsqr(1, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) + end do + end if if (w%inform%flag /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL @@ -233,87 +276,56 @@ subroutine blendenpik(A, b, n, m, sketch_size, inform, w, options) type(NLLS_options), Intent(In) :: options type(lsqr_options) :: lsqr_opt - integer :: i, lwork, action + integer :: i, lwork, action, tsize real(wp) :: norm_a lwork = size(w%work) lsqr_opt%stop_test = 2 - w%rand_ws%tau = 0.0_wp + w%rand_ws%temp_1 = 0.0_wp w%rand_ws%R = 0.0_wp w%rand_ws%SM = 0.0_wp - w%rand_ws%ATu = 0.0_wp + w%rand_ws%Sb = 0.0_wp w%rand_ws%DCT_A = 0.0_wp - ! get sketch matrix SM - call sketch(A, m, n, sketch_size, w%rand_ws%SM, options, w%rand_ws, inform) + ! get sketch matrix SM and Sb + call sketch(A, b, m, n, sketch_size, w%rand_ws%SM, w%rand_ws%Sb, & + options, w%rand_ws, inform) if (inform%status /= 0) then return end if ! now take QR decomposition of SM - call PREC(geqrf)(sketch_size, n, w%rand_ws%SM, sketch_size, w%rand_ws%tau, w%work, lwork, inform%external_return) + call PREC(geqrf)(sketch_size, n, w%rand_ws%SM, sketch_size, w%rand_ws%temp_1, w%work, lwork, inform%external_return) if (inform%external_return /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL inform%external_name = 'lapack_?geqrf' return end if - ! todo: initial guess! - ! extract R from QR decomposition + ! extract R from the output of geqrf (the upper triangle of the output matrix) do i = 1, n w%rand_ws%R(i, i:n) = w%rand_ws%SM(i, i:n) end do - call estimate_norm(A, m, n, norm_a, 15, w%lsqr_ws) - - w%lsqr_ws%x = 0.0_wp - - ! solve using LSQR - ! first call to lsqr initializes the algorithm and returns the first action - w%lsqr_ws%u = b - w%lsqr_ws%v = 0.0_wp - w%lsqr_ws%z = 0.0_wp - action = 0 - call lsqr(1, action, m, n, w%lsqr_ws%u, w%lsqr_ws%v, w%lsqr_ws%z, & - w%lsqr_ws%x, w%lsqr_ws%keep, lsqr_opt, w%lsqr_ws%inform, & - anorm=norm_a) - do while (action /= 0) - select case (action) - case (1) ! compute v = v + P^{-1} A^T u - ! we don't want to explicitly invert R (for numerical stability), so we do this in two steps: - ! first compute ATu = A^T u, then let t = R^{-1} ATu (i.e. solve R t = ATu) - call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%lsqr_ws%u, 1, 0.0_wp, w%rand_ws%ATu, 1) - call PREC(trtrs)('U', 'N', 'N', n, 1, w%rand_ws%R, n, w%rand_ws%ATu, n, inform%external_return) - if (inform%external_return /= 0) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?trtrs' - return - end if - w%lsqr_ws%v = w%lsqr_ws%v + w%rand_ws%ATu - case (2) ! compute u = u + Az - call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%lsqr_ws%z, 1, 1.0_wp, w%lsqr_ws%u, 1) - case (3) ! compute z = P^{-T} v (or equivalently, solve R^T z = v) - w%lsqr_ws%z = w%lsqr_ws%v - call PREC(trtrs)('U', 'T', 'N', n, 1, w%rand_ws%R, n, w%lsqr_ws%z, n, inform%external_return) - if (inform%external_return /= 0) then - inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_name = 'lapack_?trtrs' - return - end if - end select - call lsqr(1, action, m, n, w%lsqr_ws%u, w%lsqr_ws%v, w%lsqr_ws%z, & - w%lsqr_ws%x, w%lsqr_ws%keep, lsqr_opt, w%lsqr_ws%inform, & - anorm=norm_a) - end do - if (w%lsqr_ws%inform%flag /= 0) then + ! calculate initial guess from sketch-and-solve: x_0 = R^{-1} Q^T Sb + call PREC(ormqr)('L', 'T', sketch_size, 1, n, w%rand_ws%SM, sketch_size, w%rand_ws%temp_1, & + w%rand_ws%Sb, sketch_size, w%work, lwork, inform%external_return) + if (inform%external_return /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL - inform%external_return = w%lsqr_ws%inform%flag - inform%external_name = 'lsqr' + inform%external_name = 'lapack_?ormqr' return end if + call PREC(trtrs)('U', 'N', 'N', n, 1, w%rand_ws%R, n, w%rand_ws%Sb, n, inform%external_return) + if (inform%external_return /= 0) then + inform%status = NLLS_ERROR_FROM_EXTERNAL + inform%external_name = 'lapack_?trtrs' + return + end if + w%lsqr_ws%x(1:n) = w%rand_ws%Sb(1:n) - b(1:n) = w%lsqr_ws%x(1:n) + ! now use this initial guess in LSQR to solve the original system, using R as a preconditioner + call solve_lsqr(A, b, n, m, inform, w%lsqr_ws, options, precon=w%rand_ws%R, init_guess=.true.) end subroutine blendenpik @@ -341,11 +353,13 @@ subroutine estimate_norm(A, m, n, norm_a, num_iters, w) end subroutine estimate_norm - subroutine sketch(A, m, n, sketch_size, SA, options, w, inform) - ! Given a matrix A, return a "sketch" SA of size `sketch_size x n`. + subroutine sketch(A, b, m, n, sketch_size, SA, Sb, options, w, inform) + ! Given a matrix A, return a "sketch" SA of size `sketch_size x n`, + ! as well as the corresponding sketch Sb of b. real(wp), intent(in), contiguous :: A(:,:) + real(wp), intent(inout), contiguous :: b(:) integer, intent(in) :: m, n, sketch_size - real(wp), intent(out), contiguous :: SA(:,:) + real(wp), intent(out), contiguous :: SA(:,:), Sb(:) type(NLLS_options), intent(in) :: options type(LLS_rand_work), intent(inout) :: w type(NLLS_inform), intent(inout) :: inform @@ -357,45 +371,49 @@ subroutine sketch(A, m, n, sketch_size, SA, options, w, inform) select case (options%sketch_method) case (1) ! random sampling of rows - call select_random_rows(A, n, sketch_size, SA, w) + call select_random_rows(A, b, n, sketch_size, SA, Sb, w) case (2) ! random projection using DCT - call dct_sketch(A, m, n, sketch_size, SA, w) + call dct_sketch(A, b, m, n, sketch_size, SA, Sb, w) case default inform%status = NLLS_ERROR_BAD_SKETCH_METHOD end select end subroutine sketch - subroutine select_random_rows(A, n, s, SA, w) + subroutine select_random_rows(A, b, n, s, SA, Sb, w) ! Select `s` random rows of A to form SA real(wp), intent(in), contiguous :: A(:,:) + real(wp), intent(inout), contiguous :: b(:) integer, intent(in) :: n, s - real(wp), intent(out) :: SA(:,:) + real(wp), intent(out) :: SA(:,:), Sb(:) type(LLS_rand_work), intent(inout) :: w integer :: i, idx - call random_number(w%temp) + call random_number(w%temp_2) ! to select `s` random rows, we generate `m` random numbers in [0,1), ! and take the indices of the smallest `s` numbers do i = 1, s - idx = minloc(w%temp, dim=1) - w%temp(idx) = 1.0_wp + idx = minloc(w%temp_2, dim=1) + w%temp_2(idx) = 1.0_wp SA(i, :) = A(idx, :) + Sb(i) = b(idx) end do + end subroutine select_random_rows - subroutine dct_sketch(A, m, n, s, SA, w) + subroutine dct_sketch(A, b, m, n, s, SA, Sb, w) ! Create a sketch of A by taking a random subsample of rows ! applying random sign flips and a DCT to reduce coherence ! (i.e. the probability that the sketch is rank deficient) real(wp), intent(in), contiguous :: A(:,:) - integer, intent(in) :: m, n, s - real(wp), intent(inout) :: SA(:,:) + real(wp), intent(inout), contiguous :: b(:) + integer, intent(in) :: m, n, s + real(wp), intent(inout) :: SA(:,:), Sb(:) type(LLS_rand_work), intent(inout) :: w - integer :: i + integer :: i real(wp) :: rand_no w%DCT_A = 0.0_wp @@ -406,8 +424,10 @@ subroutine dct_sketch(A, m, n, s, SA, w) call random_number(rand_no) if (rand_no < 0.5_wp) then w%DCT_A(i, :) = -A(i, :) + w%temp_1(i) = -b(i) else w%DCT_A(i, :) = A(i, :) + w%temp_1(i) = b(i) end if end do @@ -417,9 +437,11 @@ subroutine dct_sketch(A, m, n, s, SA, w) ! and normalise DCT w%DCT_A = w%DCT_A / sqrt(2.0_wp * real(m, wp)) + w%temp_1 = w%temp_1 / sqrt(2.0_wp * real(m, wp)) ! and finaly select `s` random rows of DCT_A to form SA - call select_random_rows(w%DCT_A, n, s, SA, w) + ! and the corresponding entries of b to form Sb + call select_random_rows(w%DCT_A, w%temp_1, n, s, SA, Sb, w) end subroutine dct_sketch diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index f3202bbe..ab86ac64 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -623,7 +623,7 @@ module MODULE_PREC(ral_nlls_workspaces) type, public :: LLS_lsqr_work logical :: allocated = .false. ! variables for LSQR - real(wp), allocatable :: u(:), v(:), z(:), x(:) + real(wp), allocatable :: u(:), v(:), z(:), x(:), ATu(:) ! LSQR reverse communication objects type(lsqr_keep) :: keep @@ -632,7 +632,7 @@ module MODULE_PREC(ral_nlls_workspaces) type, public :: LLS_rand_work logical :: allocated = .false. - real(wp), allocatable :: tau(:), temp(:), ATu(:) + real(wp), allocatable :: temp_1(:), temp_2(:), Sb(:) real(wp), allocatable :: R(:,:), SM(:,:), DCT_A(:,:) complex(wp), allocatable :: dct_work(:) end type LLS_rand_work @@ -1335,7 +1335,7 @@ subroutine setup_workspace_LLS_lsqr(w, n, m, inform) type( nlls_inform ), intent(inout) :: inform inform%status = 0 - allocate(w%u(m), w%v(n), w%z(n), w%x(n), stat=inform%alloc_status) + allocate(w%u(m), w%v(n), w%z(n), w%x(n), w%ATu(n), stat=inform%alloc_status) If (inform%alloc_status /= 0) Then inform%bad_alloc = "setup_workspace_LLS_lsqr" inform%status = NLLS_ERROR_ALLOCATION @@ -1354,6 +1354,7 @@ subroutine remove_workspace_LLS_lsqr(w) if(allocated(w%v)) deallocate(w%v, stat=ierr_dummy) if(allocated(w%z)) deallocate(w%z, stat=ierr_dummy) if(allocated(w%x)) deallocate(w%x, stat=ierr_dummy) + if(allocated(w%ATu)) deallocate(w%ATu, stat=ierr_dummy) w%allocated = .false. end subroutine remove_workspace_LLS_lsqr @@ -1366,8 +1367,8 @@ subroutine setup_workspace_LLS_rand(w, n, m, options, inform) type( nlls_inform ), intent(inout) :: inform inform%status = 0 - allocate(w%tau(min(options%sketch_size, n)), w%temp(m), w%R(n,n), & - w%SM(options%sketch_size, n), w%ATu(n), w%DCT_A(m,n), & + allocate(w%temp_1(m), w%temp_2(m), w%Sb(options%sketch_size), w%R(n,n), & + w%SM(options%sketch_size, n), w%DCT_A(m,n), & w%dct_work(m), stat=inform%alloc_status) If (inform%alloc_status /= 0) Then inform%bad_alloc = "setup_workspace_LLS_rand" @@ -1383,11 +1384,10 @@ subroutine remove_workspace_LLS_rand(w) type( LLS_rand_work ), intent(out) :: w Integer :: ierr_dummy - if(allocated(w%tau)) deallocate(w%tau, stat=ierr_dummy) - if(allocated(w%temp)) deallocate(w%temp, stat=ierr_dummy) + if(allocated(w%temp_1)) deallocate(w%temp_1, stat=ierr_dummy) + if(allocated(w%temp_2)) deallocate(w%temp_2, stat=ierr_dummy) if(allocated(w%R)) deallocate(w%R, stat=ierr_dummy) if(allocated(w%SM)) deallocate(w%SM, stat=ierr_dummy) - if(allocated(w%ATu)) deallocate(w%ATu, stat=ierr_dummy) if(allocated(w%DCT_A)) deallocate(w%DCT_A, stat=ierr_dummy) if(allocated(w%dct_work)) deallocate(w%dct_work, stat=ierr_dummy) From 5e85c74f2c5903d259969f7141ef071bd0d6929d Mon Sep 17 00:00:00 2001 From: alexhroom Date: Wed, 11 Mar 2026 12:19:05 +0000 Subject: [PATCH 14/16] improved documentation --- libRALFit/doc/common/errors.rst | 4 ++++ libRALFit/doc/common/options.rst | 4 ++++ libRALFit/src/ral_nlls_linear.F90 | 27 ++++++++++----------------- libRALFit/src/ral_nlls_workspaces.F90 | 4 ++-- libRALFit/test/unit_test_mod.F90 | 2 +- 5 files changed, 21 insertions(+), 20 deletions(-) diff --git a/libRALFit/doc/common/errors.rst b/libRALFit/doc/common/errors.rst index 079da402..a03319f2 100644 --- a/libRALFit/doc/common/errors.rst +++ b/libRALFit/doc/common/errors.rst @@ -42,6 +42,10 @@ Possible values are: - Weights vector must be sufficiently positive * - -21 - Unsupported LLS solver (``lls_solver``) specified in options. + * - -22 + - Bad sketch size (``sketch_size``) specified in options for randomised method. + * - -23 + - Unsupported value of ``sketch_method`` specified in options for randomised method. * - -101 - Unsupported model in dogleg (``nlls_method=1``). * - -201 diff --git a/libRALFit/doc/common/options.rst b/libRALFit/doc/common/options.rst index 4306798f..facae896 100644 --- a/libRALFit/doc/common/options.rst +++ b/libRALFit/doc/common/options.rst @@ -86,6 +86,10 @@ .. |lls_solver| replace:: specifies which linear least squares solver to use internally. +.. |sketch_size| replace:: if using randomised solver (``lls_solver = 3``), specifies the size of the sketch to use. + +.. |sketch_method| replace:: if using randomised solver (``lls_solver = 3``), specifies the method to use to perform the sketch. + .. |more_sorensen_maxits| replace:: if ``nlls_method = 3``, specifies the maximum number of iterations allowed in the More-Sorensen method. .. |more_sorensen_shift| replace:: if ``nlls_method = 3``, specifies the shift to be used in the More-Sorensen method. diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 6e14e2d3..4a44d649 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -65,26 +65,16 @@ subroutine solve_LLS(A, b, n, m, inform, w, options, pd) return end if call solve_lsqr(A, b, n, m, inform, w%lsqr_ws, options) - ! with LAPACK fallback... - if (inform%status /= 0) then - if (options%allow_fallback_method) then - inform%status = 0 - call solve_gels(A,b,n,m,inform,w,options) - else - return - end if - end if - case (3) ! Randomised method (Blendenpik) + case (3) ! Randomised method (sketch-and-precondition) if (.not. w%lsqr_ws%allocated .or. .not. w%rand_ws%allocated) then inform%status = NLLS_ERROR_WORKSPACE_ERROR return end if - call blendenpik(A, b, n, m, options%sketch_size, inform, w, options) - ! if Blendenpik fails, fall back to LAPACK + call solve_rand(A, b, n, m, options%sketch_size, inform, w, options) case default inform%status = NLLS_ERROR_BAD_LLS_SOLVER end select - ! if LSQR or Blendenpik fails, fallback to LAPACK GELS + ! if LSQR or randomised fails, fallback to LAPACK GELS if (inform%status /= 0 .and. options%lls_solver > 1 .and. options%allow_fallback_method) then inform%status = 0 call solve_gels(A,b,n,m,inform,w,options) @@ -263,12 +253,15 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) end subroutine solve_lsqr - subroutine blendenpik(A, b, n, m, sketch_size, inform, w, options) - ! Blendenpik (sketch-and-precondition) randomised solver - ! Source: + subroutine solve_rand(A, b, n, m, sketch_size, inform, w, options) + ! Sketch-and-precondition randomised solver with sketch-and-solve initialisation + ! Sources: ! P. Avron, P. Maymounkov, S. Toledo, ! "Blendenpik: Supercharging LAPACK's Least-Squares Solver" ! URL: https://pdos.csail.mit.edu/~petar/papers/blendenpik-v1.pdf + ! Meier, M., Nakatsukasa, Y., Townsend, A., Webb, M. + ! "Are sketch-and-precondition linear solvers numerically stable?" + ! URL: https://arxiv.org/abs/2302.07202 real(wp), intent(inout), contiguous :: A(:,:), b(:) integer, intent(in) :: sketch_size, m, n type(NLLS_inform), INTENT(INOUT) :: inform @@ -327,7 +320,7 @@ subroutine blendenpik(A, b, n, m, sketch_size, inform, w, options) ! now use this initial guess in LSQR to solve the original system, using R as a preconditioner call solve_lsqr(A, b, n, m, inform, w%lsqr_ws, options, precon=w%rand_ws%R, init_guess=.true.) - end subroutine blendenpik + end subroutine solve_rand subroutine estimate_norm(A, m, n, norm_a, num_iters, w) ! power method to estimate spectral norm of A diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index ab86ac64..7e3a348a 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -291,11 +291,11 @@ module MODULE_PREC(ral_nlls_workspaces) ! which linear least squares solver should we use? ! 1 LAPACK solver ! 2 LSQR (iterative) -! 3 Randomised solver (Blendenpik) +! 3 Randomised solver (sketch-and-precondition) INTEGER :: lls_solver = 1 -! If using randomised linear solver (Blendenpik), what sketch size? +! If using randomised linear solver, what sketch size? ! this must be manually set by the user as it depends on `m` INTEGER :: sketch_size = 0 diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index 38e7118f..4b2be866 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -2582,7 +2582,7 @@ subroutine solve_LLS_randomised_tests(options,fails) end if call mult_J(J,n,m,d,Jd,.True.) normerror = norm2(Jd - f) - if ( .not. normerror < 1.0e-11_wp ) then + if ( .not. normerror < 1.0e-12_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS randomised test failed: wrong solution returned' write(*,*) '||Jd - f|| = ', normerror From 2a88286e7546febfada0e0efc950e41a7e3540a1 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Fri, 20 Mar 2026 16:00:41 +0000 Subject: [PATCH 15/16] added main test --- libRALFit/test/nlls_test.F90 | 31 +++++++ libRALFit/test/unit_test_mod.F90 | 137 +++++++++++++++++++++++++++++++ 2 files changed, 168 insertions(+) diff --git a/libRALFit/test/nlls_test.F90 b/libRALFit/test/nlls_test.F90 index ffbf9f77..97eeee55 100644 --- a/libRALFit/test/nlls_test.F90 +++ b/libRALFit/test/nlls_test.F90 @@ -172,6 +172,37 @@ program nlls_test end do end do + ! Test linear solvers are consistent + do tr_update = 1,2 + do nlls_method = 1,4 + call reset_default_options(options) + options%print_options = .True. + options%allow_fallback_method = .false. + options%nlls_method = nlls_method + options%tr_update_strategy = tr_update + options%model = 1 + options%exact_second_derivatives = .true. + options%output_progress_vectors = .true. + call print_line(options%out) + write(options%out,*) "Testing linear solver consistency" + write(options%out,*) "tr_update_strategy = ", options%tr_update_strategy + write(options%out,*) "nlls_method = ", options%nlls_method + write(options%out,*) "model = ", options%model + write(options%out,*) "Tolerance type for resvec = ", tol_type + call print_line(options%out) + if (nlls_method == 4) then + do inner_method = 1,3 + ! check the tests with c and fortran jacobians + ! pass individually, and give consistent results. + options%inner_method = inner_method + call linear_solve_tests(options,no_errors_main, tol_type) + end do + else + call linear_solve_tests(options,no_errors_main, tol_type) + end if + end do + end do + ! dogleg, no fallback call reset_default_options(options) diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index 4b2be866..7c0468ef 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -1435,6 +1435,143 @@ subroutine c_fortran_tests(options, fails, tol_type) fails = fails + 1 end if end subroutine c_fortran_tests + + subroutine linear_solve_tests(options, fails, tol_type) + Implicit None + type( NLLS_options ), intent(inout) :: options + integer, intent(inout) :: fails + Character(len=*), intent(in) :: tol_type + real(wp), allocatable :: x(:) + type( user_type ), target :: params + type( NLLS_inform ) :: status, c_status + real(wp) :: resvec_error, solver_type_error, abstol, reltol + real(wp) :: C_results(3), F_results(3) + integer :: n, lls_solver + Logical :: Ok, any_fails + + Continue + + abstol = 2.0e-8_wp + reltol = -1.0_wp + + any_fails = .false. + + if (trim(tol_type) == 'both') then + reltol = 1.0e-8_wp + end if + + n = 2 + + allocate( x(n) ) + options%print_level = 5 + + do lls_solver = 1, 3 + call generate_data_example(params) + options%lls_solver = lls_solver + options%fortran_jacobian = .true. + call solve_basic(X,params,options,status) + if ( status%status .ne. 0 ) then + write(*,*) 'nlls_solve failed to converge:' + write(*,*) status%error_message + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'LLS_SOLVER = ', options%lls_solver + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + write(*,*) 'info%status = ', status%status + write(*,*) 'scale? ', options%scale + fails = fails + 1 + any_fails = .true. + else + F_results(lls_solver) = norm2(status%resvec(1:status%iter+1)) + end if + + ! run the same test with a c jacobian + options%fortran_jacobian = .false. + call solve_basic_c(X,params,options,c_status) + if ( c_status%status .ne. 0 ) then + write(*,*) 'nlls_solve (c jac) failed to converge:' + write(*,*) c_status%error_message + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'LLS_SOLVER = ', options%lls_solver + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + write(*,*) 'info%status = ', c_status%status + write(*,*) 'scale? ', options%scale + fails = fails + 1 + any_fails = .true. + else + C_results(lls_solver) = norm2(c_status%resvec(1:c_status%iter+1)) + end if + + ! check that the results are consistent with + ! both c and fortran jacobians + if ( c_status%iter == status%iter ) then + resvec_error = norm2(c_status%resvec(1:c_status%iter+1) - & + status%resvec(1:status%iter+1)) + if (.not. resvec_error < abstol) then + write(*,*) 'Warning: fortran and c resvec differ!' + write(*,*) 'Different resvecs in 2-norm for' + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'LLS_SOLVER = ', options%lls_solver + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + ! Report all entries and show how much they differ in L2 and relative L2 + ! Ok will inform is at least reltol is met + Call check_resvec(c_status%resvec(1:c_status%iter+1), & + status%resvec(1:status%iter+1), atol=abstol, reltol=reltol, prn=.true., Ok=Ok) + fails = fails + merge(0,1,Ok) + end if + else + write(*,*) 'error: fortran and c jacobians' + write(*,*) 'took different numbers of iterations: F iters=', status%iter, "C iters=", c_status%iter + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'LLS_SOLVER = ', options%lls_solver + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + fails = fails + 1 + end if + end do + if (.not. any_fails) then + do lls_solver = 1, 2 + do n = lls_solver + 1, 3 + solver_type_error = abs(C_results(lls_solver) - C_results(n)) + if (.not. solver_type_error < abstol) then + write(*,*) 'Warning: different linear solvers give different results' + write(*,*) 'Different resvecs in 2-norm for C Jacobian in' + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + write(*,*) 'LLS_SOLVER 1 = ', lls_solver, " residual ", C_results(lls_solver) + write(*,*) 'LLS_SOLVER 2 = ', n, " residual ", C_results(n) + fails = fails + 1 + end if + end do + end do + + do lls_solver = 1, 2 + do n = lls_solver + 1, 3 + solver_type_error = abs(F_results(lls_solver) - F_results(n)) + if (.not. solver_type_error < abstol) then + write(*,*) 'Warning: different linear solvers give different results' + write(*,*) 'Different resvecs in 2-norm for F Jacobian in' + write(*,*) 'NLLS_METHOD = ', options%nlls_method + write(*,*) 'MODEL = ', options%model + write(*,*) 'TR_UPDATE = ', options%tr_update_strategy + write(*,*) 'inner_method = ', options%inner_method + write(*,*) 'LLS_SOLVER 1 = ', lls_solver, " residual ", F_results(lls_solver) + write(*,*) 'LLS_SOLVER 2 = ', n, " residual ", F_results(n) + fails = fails + 1 + end if + end do + end do + end if + + end subroutine linear_solve_tests subroutine dogleg_tests(options,fails) From 00f35a2482f525cf0c9bbead845b253933ccd880 Mon Sep 17 00:00:00 2001 From: alexhroom Date: Fri, 20 Mar 2026 16:59:03 +0000 Subject: [PATCH 16/16] made new solvers work with C jacobian --- libRALFit/src/ral_nlls_internal.F90 | 2 - libRALFit/src/ral_nlls_linear.F90 | 72 +++++++++++++++-------- libRALFit/src/ral_nlls_workspaces.F90 | 14 ++++- libRALFit/test/unit_test_mod.F90 | 85 +++++++++++++++++++-------- 4 files changed, 120 insertions(+), 53 deletions(-) diff --git a/libRALFit/src/ral_nlls_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index 0333dc48..c64bc512 100644 --- a/libRALFit/src/ral_nlls_internal.F90 +++ b/libRALFit/src/ral_nlls_internal.F90 @@ -4233,8 +4233,6 @@ Subroutine check_options(opt, inform) inform%status = NLLS_ERROR_PRINT_LEVEL ElseIf (opt%box_linesearch_type<1 .Or. opt%box_linesearch_type>2) Then inform%status = NLLS_ERROR_UNSUPPORTED_LINESEARCH - ElseIf (opt%lls_solver == 2 .and. opt%sketch_size <= 0) then - inform%status = NLLS_ERROR_BAD_SKETCH_SIZE End If End Subroutine check_options diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 index 4a44d649..45483529 100644 --- a/libRALFit/src/ral_nlls_linear.F90 +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -165,7 +165,8 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) real(wp), intent(in), optional, contiguous :: precon(:,:) logical, intent(in), optional :: init_guess - integer :: action, use_precon + integer :: action, use_precon, rows, cols + character(len=1) :: trans, no_trans real(wp) :: norm_a type(lsqr_options) :: lsqr_opt @@ -174,7 +175,19 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) ! stopping criterion requires an estimate of the spectral norm of A, ! which we can get from a few iterations of the power method - call estimate_norm(A, m, n, norm_a, 15, w) + call estimate_norm(A, m, n, norm_a, 15, w, options%fortran_jacobian) + + if (options%fortran_jacobian) then + rows = m + cols = n + trans = 'T' + no_trans = 'N' + else + rows = n + cols = m + trans = 'N' + no_trans = 'T' + end if w%u = b w%v = 0.0_wp @@ -184,7 +197,7 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) if (present(init_guess)) then if (init_guess) then ! update `u` for LSQR to be b - A x_0, so that we can start LSQR with this initial guess - call PREC(gemv)('N', m, n, -1.0_wp, A, m, w%x, 1, 1.0_wp, w%u, 1) + call PREC(gemv)(no_trans, rows, cols, -1.0_wp, A, rows, w%x, 1, 1.0_wp, w%u, 1) else w%x = 0.0_wp end if @@ -206,9 +219,9 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) do while (action /= 0) select case (action) case (1) ! compute v = v + A^T u - call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 1.0_wp, w%v, 1) + call PREC(gemv)(trans, rows, cols, 1.0_wp, A, rows, w%u, 1, 1.0_wp, w%v, 1) case (2) ! compute u = u + A v - call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 1.0_wp, w%u, 1) + call PREC(gemv)(no_trans, rows, cols, 1.0_wp, A, rows, w%v, 1, 1.0_wp, w%u, 1) end select call lsqr(0, action, m, n, w%u, w%v, w%z, w%x, w%keep, lsqr_opt, w%inform, anorm=norm_a) end do @@ -218,7 +231,7 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) case (1) ! compute v = v + P^{-1} A^T u ! we don't want to explicitly invert R (for numerical stability), so we do this in two steps: ! first compute ATu = A^T u, then let t = R^{-1} ATu (i.e. solve R t = ATu) - call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 0.0_wp, w%ATu, 1) + call PREC(gemv)(trans, rows, cols, 1.0_wp, A, rows, w%u, 1, 0.0_wp, w%ATu, 1) call PREC(trtrs)('U', 'N', 'N', n, 1, precon, n, w%ATu, n, inform%external_return) if (inform%external_return /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL @@ -227,7 +240,7 @@ subroutine solve_lsqr(A, b, n, m, inform, w, options, precon, init_guess) end if w%v = w%v + w%ATu case (2) ! compute u = u + Az - call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%z, 1, 1.0_wp, w%u, 1) + call PREC(gemv)(no_trans, rows, cols, 1.0_wp, A, rows, w%z, 1, 1.0_wp, w%u, 1) case (3) ! compute z = P^{-T} v (or equivalently, solve R^T z = v) w%z = w%v call PREC(trtrs)('U', 'T', 'N', n, 1, precon, n, w%z, n, inform%external_return) @@ -269,7 +282,7 @@ subroutine solve_rand(A, b, n, m, sketch_size, inform, w, options) type(NLLS_options), Intent(In) :: options type(lsqr_options) :: lsqr_opt - integer :: i, lwork, action, tsize + integer :: i, lwork, action, tsize, s real(wp) :: norm_a lwork = size(w%work) @@ -281,15 +294,24 @@ subroutine solve_rand(A, b, n, m, sketch_size, inform, w, options) w%rand_ws%Sb = 0.0_wp w%rand_ws%DCT_A = 0.0_wp + if (sketch_size == -1) then ! default + s = 4 * n + else if (sketch_size <= n .or. sketch_size > m) then + write(*,*) "Sketch size must be between n and m; setting default of 4n" + s = 4 * n + else + s = sketch_size + end if + ! get sketch matrix SM and Sb - call sketch(A, b, m, n, sketch_size, w%rand_ws%SM, w%rand_ws%Sb, & + call sketch(A, b, m, n, s, w%rand_ws%SM, w%rand_ws%Sb, & options, w%rand_ws, inform) if (inform%status /= 0) then return end if ! now take QR decomposition of SM - call PREC(geqrf)(sketch_size, n, w%rand_ws%SM, sketch_size, w%rand_ws%temp_1, w%work, lwork, inform%external_return) + call PREC(geqrf)(s, n, w%rand_ws%SM, s, w%rand_ws%temp_1, w%work, lwork, inform%external_return) if (inform%external_return /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL inform%external_name = 'lapack_?geqrf' @@ -302,8 +324,8 @@ subroutine solve_rand(A, b, n, m, sketch_size, inform, w, options) end do ! calculate initial guess from sketch-and-solve: x_0 = R^{-1} Q^T Sb - call PREC(ormqr)('L', 'T', sketch_size, 1, n, w%rand_ws%SM, sketch_size, w%rand_ws%temp_1, & - w%rand_ws%Sb, sketch_size, w%work, lwork, inform%external_return) + call PREC(ormqr)('L', 'T', s, 1, n, w%rand_ws%SM, s, w%rand_ws%temp_1, & + w%rand_ws%Sb, s, w%work, lwork, inform%external_return) if (inform%external_return /= 0) then inform%status = NLLS_ERROR_FROM_EXTERNAL inform%external_name = 'lapack_?ormqr' @@ -322,12 +344,13 @@ subroutine solve_rand(A, b, n, m, sketch_size, inform, w, options) end subroutine solve_rand - subroutine estimate_norm(A, m, n, norm_a, num_iters, w) + subroutine estimate_norm(A, m, n, norm_a, num_iters, w, fortran_jacobian) ! power method to estimate spectral norm of A real(wp), intent(in), contiguous :: A(:,:) integer, intent(in) :: m, n, num_iters real(wp), intent(out) :: norm_a type(LLS_lsqr_work), intent(inout) :: w + logical, intent(in) :: fortran_jacobian integer :: i @@ -336,11 +359,19 @@ subroutine estimate_norm(A, m, n, norm_a, num_iters, w) ! each iteration of the power method, we multiply by A^T A and renormalize ! this is done in two steps to be better-conditioned - do i = 1, num_iters - w%v = w%v / norm2(w%v) - call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 0.0_wp, w%u, 1) - call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 0.0_wp, w%v, 1) - end do + if (fortran_jacobian) then + do i = 1, num_iters + w%v = w%v / norm2(w%v) + call PREC(gemv)('N', m, n, 1.0_wp, A, m, w%v, 1, 0.0_wp, w%u, 1) + call PREC(gemv)('T', m, n, 1.0_wp, A, m, w%u, 1, 0.0_wp, w%v, 1) + end do + else + do i = 1, num_iters + w%v = w%v / norm2(w%v) + call PREC(gemv)('T', n, m, 1.0_wp, A, n, w%v, 1, 0.0_wp, w%u, 1) + call PREC(gemv)('N', n, m, 1.0_wp, A, n, w%u, 1, 0.0_wp, w%v, 1) + end do + end if norm_a = sqrt(norm2(w%v)) @@ -357,11 +388,6 @@ subroutine sketch(A, b, m, n, sketch_size, SA, Sb, options, w, inform) type(LLS_rand_work), intent(inout) :: w type(NLLS_inform), intent(inout) :: inform - if (sketch_size <= 0 .or. sketch_size > m) then - inform%status = NLLS_ERROR_BAD_SKETCH_SIZE - return - end if - select case (options%sketch_method) case (1) ! random sampling of rows call select_random_rows(A, b, n, sketch_size, SA, Sb, w) diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index 7e3a348a..3b090ba8 100644 --- a/libRALFit/src/ral_nlls_workspaces.F90 +++ b/libRALFit/src/ral_nlls_workspaces.F90 @@ -298,7 +298,7 @@ module MODULE_PREC(ral_nlls_workspaces) ! If using randomised linear solver, what sketch size? ! this must be manually set by the user as it depends on `m` - INTEGER :: sketch_size = 0 + INTEGER :: sketch_size = -1 ! If using randomised solver, which sketching method should we use? ! 1 Uniform subsample @@ -1366,9 +1366,17 @@ subroutine setup_workspace_LLS_rand(w, n, m, options, inform) type( nlls_options ), intent(in) :: options type( nlls_inform ), intent(inout) :: inform + integer :: s + + if (options%sketch_size <= n .or. options%sketch_size > m) then + s = 4 * n + else + s = options%sketch_size + end if + inform%status = 0 - allocate(w%temp_1(m), w%temp_2(m), w%Sb(options%sketch_size), w%R(n,n), & - w%SM(options%sketch_size, n), w%DCT_A(m,n), & + allocate(w%temp_1(m), w%temp_2(m), w%Sb(s), w%R(n,n), & + w%SM(s, n), w%DCT_A(m,n), & w%dct_work(m), stat=inform%alloc_status) If (inform%alloc_status /= 0) Then inform%bad_alloc = "setup_workspace_LLS_rand" diff --git a/libRALFit/test/unit_test_mod.F90 b/libRALFit/test/unit_test_mod.F90 index 7c0468ef..576522de 100644 --- a/libRALFit/test/unit_test_mod.F90 +++ b/libRALFit/test/unit_test_mod.F90 @@ -1445,7 +1445,7 @@ subroutine linear_solve_tests(options, fails, tol_type) type( user_type ), target :: params type( NLLS_inform ) :: status, c_status real(wp) :: resvec_error, solver_type_error, abstol, reltol - real(wp) :: C_results(3), F_results(3) + real(wp) :: C_results(4), F_results(4) integer :: n, lls_solver Logical :: Ok, any_fails @@ -1465,9 +1465,15 @@ subroutine linear_solve_tests(options, fails, tol_type) allocate( x(n) ) options%print_level = 5 - do lls_solver = 1, 3 + do lls_solver = 1, 4 + ! we use 4 for solver 3 with sketch method 2 call generate_data_example(params) - options%lls_solver = lls_solver + if (lls_solver == 4) then + options%lls_solver = 3 + options%sketch_method = 2 + else + options%lls_solver = lls_solver + end if options%fortran_jacobian = .true. call solve_basic(X,params,options,status) if ( status%status .ne. 0 ) then @@ -2571,13 +2577,15 @@ subroutine solve_LLS_lsqr_tests(options, fails) real(wp), allocatable :: J(:,:), J_copy(:,:), Jd(:), f(:), d(:), x_calc(:), x_true(:) real(wp) :: normerror integer :: n, m, i, k, l - type( nlls_workspace ) :: work + type( nlls_workspace ), target :: work type( nlls_workspace ), Target :: iw + type( solve_lls_work ), pointer :: ws type( nlls_inform ) :: inform fails = 0 work%iw_ptr => iw iw%iw_ptr => iw + ws => work%calculate_step_ws%dogleg_ws%solve_lls_ws n = 5 m = 30 @@ -2609,7 +2617,7 @@ subroutine solve_LLS_lsqr_tests(options, fails) options%lls_solver = 2 ! LSQR options%allow_fallback_method = .false. - call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + call solve_LLS(J_copy,d,n,m,inform,ws,options,.false.) if ( inform%status .ne. NLLS_ERROR_WORKSPACE_ERROR ) then write (*, *) 'solve_LLS LSQR test failed: wrong error message returned when workspace unallocated' write(*,*) 'status = ', inform%status, " (expected ",NLLS_ERROR_WORKSPACE_ERROR,")" @@ -2621,7 +2629,7 @@ subroutine solve_LLS_lsqr_tests(options, fails) ! test correct residual - call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + call solve_LLS(J_copy,d,n,m,inform,ws,options,.false.) if ( inform%status .ne. 0 ) then write (*, *) 'solve_LLS LSQR test failed: error message returned' write(*,*) 'status = ', inform%status, " (expected 0)" @@ -2635,12 +2643,50 @@ subroutine solve_LLS_lsqr_tests(options, fails) ! check answer call mult_J(J,n,m,d,Jd,.True.) normerror = norm2(Jd - f) - if ( .not. normerror < 1.0e-12_wp ) then + if ( .not. normerror < 1.0e-11_wp ) then ! wrong answer, as data chosen to fit write(*,*) 'solve_LLS LSQR test failed: wrong solution returned' write(*,*) '||Jd - f|| = ', normerror fails = fails + 1 end if + + deallocate(J, J_copy) + allocate(J(n,m), J_copy(n,m)) + + ! now check for C jacobian + do k = 1, m*n + i = ((k-1) / n) + 1 + l = mod(k-1, n) + 1 + J(l, i) = real(k + i*l, wp) + end do + + J_copy = J + d = f + + options%fortran_jacobian = .false. + + call solve_LLS(J_copy,d,n,m,inform,ws,options,.false.) + if ( inform%status .ne. 0 ) then + write (*, *) 'solve_LLS LSQR test failed (C Jacobian): error message returned' + write(*,*) 'status = ', inform%status, " (expected 0)" + if (inform%status == NLLS_ERROR_FROM_EXTERNAL) then + write(*,*) 'external return = ', inform%external_return + write(*,*) 'external name = ', inform%external_name + end if + fails = fails + 1 + end if + + ! check answer + call mult_J(J,n,m,d,Jd,.false.) + normerror = norm2(Jd - f) + if ( .not. normerror < 1.0e-11_wp ) then + ! wrong answer, as data chosen to fit + write(*,*) 'solve_LLS LSQR test failed (C Jacobian): wrong solution returned' + write(*,*) '||Jd - f|| = ', normerror + fails = fails + 1 + end if + + deallocate(J, J_copy, f, d, Jd) call nlls_finalize(work,options,inform) call reset_default_options(options) @@ -2653,13 +2699,15 @@ subroutine solve_LLS_randomised_tests(options,fails) real(wp), allocatable :: J(:,:), J_copy(:,:), Jd(:), d(:), f(:) real(wp) :: normerror integer :: n,m,k,i,l,sketch_method - type( nlls_workspace ) :: work + type( nlls_workspace ), target :: work type( nlls_workspace ), Target :: iw + type( solve_lls_work ), pointer :: ws type( nlls_inform ) :: inform fails = 0 work%iw_ptr => iw iw%iw_ptr => iw + ws => work%calculate_step_ws%dogleg_ws%solve_lls_ws n = 5 m = 60 @@ -2667,25 +2715,11 @@ subroutine solve_LLS_randomised_tests(options,fails) options%nlls_method = 1 ! dogleg (just so we know what workspace we're in) options%lls_solver = 3 ! randomised solver options%allow_fallback_method = .false. + options%sketch_size = 10 allocate(J(m,n), J_copy(m,n), Jd(m), d(m), f(m)) call setup_workspaces(work,n,m,options,inform) - - ! ensure solver gives error if sketch size is not a positive integer - options%sketch_size = 0 - call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) - if ( inform%status .ne. NLLS_ERROR_BAD_SKETCH_SIZE ) then - write (*, *) 'solve_LLS randomised test failed: wrong error message returned' - write(*,*) 'status = ', inform%status, " (expected ",NLLS_ERROR_BAD_SKETCH_SIZE,")" - fails = fails + 1 - end if - inform%status = 0 - - call remove_workspaces(work,options) - - options%sketch_size = 10 ! reset to valid sketch size - ! this is mostly a generic version of the matrices in the other tests; ! i.e. J = [ 1, 2, 3, 4, 5...] reshaped to (n, m) ! but in this case we need a bigger matrix to ensure we get a decent sketch @@ -2706,7 +2740,7 @@ subroutine solve_LLS_randomised_tests(options,fails) J_copy = J d = f - call solve_LLS(J_copy,d,n,m,inform,work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) + call solve_LLS(J_copy,d,n,m,inform,ws,options,.false.) if ( inform%status .ne. 0 ) then write (*, *) 'solve_LLS randomised test failed: error message returned' write(*,*) 'status = ', inform%status, " (expected 0)" @@ -2726,8 +2760,9 @@ subroutine solve_LLS_randomised_tests(options,fails) write(*,*) 'sketch method = ', sketch_method fails = fails + 1 end if - call remove_workspaces(work,options) end do + + deallocate(J, J_copy, f, d, Jd) call nlls_finalize(work,options,inform) call reset_default_options(options)