diff --git a/libRALFit/doc/common/errors.rst b/libRALFit/doc/common/errors.rst index b0b35b16..a03319f2 100644 --- a/libRALFit/doc/common/errors.rst +++ b/libRALFit/doc/common/errors.rst @@ -40,6 +40,12 @@ 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. + * - -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 43cfb631..facae896 100644 --- a/libRALFit/doc/common/options.rst +++ b/libRALFit/doc/common/options.rst @@ -84,6 +84,12 @@ .. |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. + +.. |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/CMakeLists.txt b/libRALFit/src/CMakeLists.txt index e3b93242..77318df1 100644 --- a/libRALFit/src/CMakeLists.txt +++ b/libRALFit/src/CMakeLists.txt @@ -14,10 +14,16 @@ 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) +# 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") @@ -92,6 +98,11 @@ 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}" "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 target_link_libraries("RAL_NLLS_DTRS_${precision}" "RAL_NLLS_TYPES_${precision}" @@ -113,19 +124,29 @@ 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}" "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() # Post ################ 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..0a2ece1a --- /dev/null +++ b/libRALFit/src/external/CMakeLists.txt @@ -0,0 +1,11 @@ +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) + +# 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 new file mode 100644 index 00000000..09d87af1 --- /dev/null +++ b/libRALFit/src/external/dct.f90 @@ -0,0 +1,95 @@ +! Interface to DCT implementations. +#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 :: FFTW_FORWARD + parameter (FFTW_FORWARD = -1) + integer :: FFTW_ESTIMATE + parameter (FFTW_ESTIMATE = 64) + + interface dct + 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 + ! 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, work) + ! Compute the 1D DCT across the columns of A using the FFTZ library. + ! `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 + + plan = 0 + + do i = 1, n + call dct_fftz_1d(A(:, i), m, work) + end do + + 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 new file mode 100644 index 00000000..ee1bae11 --- /dev/null +++ b/libRALFit/src/external/lsqr.f90 @@ -0,0 +1,1500 @@ +! 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. +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +#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, 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 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) + + ! 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 PREC(scal)(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 PREC(scal)(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 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. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call PREC(scal)(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 PREC(scal)(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 PREC(scal)(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 PREC(scal)(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 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. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call PREC(scal)(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 PREC(scal)(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 PREC(scal)(n, (one/keep%alpha), v, 1) + call PREC(scal)(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 PREC(scal)(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 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. & + keep%itn_count .lt. keep%localVecs) & + call localVEnqueue(n, v, keep) + end if + + call PREC(scal)(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 PREC(scal)(n, (one/keep%alpha), v, 1) + call PREC(scal)(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 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 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 PREC(axpy)(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 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_internal.F90 b/libRALFit/src/ral_nlls_internal.F90 index 396895c4..c64bc512 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 implicit none @@ -22,13 +23,12 @@ 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 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 @@ -1803,8 +1803,11 @@ SUBROUTINE dogleg(J,f,hf,g,n,m,Delta,d,normd,options,inform,w) 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) + 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 + w%d_gn = -w%solve_LLS_ws%temp(1:n) case default inform%status = NLLS_ERROR_DOGLEG_MODEL goto 100 @@ -1867,6 +1870,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 @@ -1892,15 +1898,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(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) @@ -1957,19 +1959,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(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) @@ -1982,25 +1974,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(:,:) Do i = 1, n w%B(i,i) = w%B(i,i) + lam + w%p1(i) = -v(i) 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(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 @@ -2087,6 +2069,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 @@ -2105,11 +2094,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(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) End If if (inform%status == 0) then ! A is symmetric positive definite.... @@ -2119,6 +2109,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 @@ -2126,10 +2120,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) @@ -2255,8 +2245,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(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) end if if (inform%status /= 0) goto 100 @@ -2279,7 +2270,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') @@ -2387,8 +2378,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(w%LtL,d,n,n,inform,w%lls_ws,options,.true.) else factorization_done = .false. end if @@ -2435,10 +2428,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 @@ -2447,6 +2436,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 @@ -2491,7 +2484,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 @@ -2518,6 +2512,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 @@ -2558,8 +2557,10 @@ 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 + d(1:n) = -v(1:n) + call solve_LLS(w%AplusSigma,d,n,n,inform,w%lls_ws,options,.true.) + if (inform%status == 0) then region = 2 sigma = sigma + sigma_shift @@ -2591,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 @@ -2722,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(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 @@ -2737,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 @@ -2887,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 @@ -2908,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(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 @@ -2924,54 +2937,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) ! ----------------------------------------------------------------- @@ -3579,76 +3544,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(:) @@ -4338,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 -! ... End If End Subroutine check_options diff --git a/libRALFit/src/ral_nlls_linear.F90 b/libRALFit/src/ral_nlls_linear.F90 new file mode 100644 index 00000000..45483529 --- /dev/null +++ b/libRALFit/src/ral_nlls_linear.F90 @@ -0,0 +1,467 @@ +! 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" + +module MODULE_PREC(ral_nlls_linear) + 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, dct1d + + implicit none + + private + public :: solve_LLS + +contains + ! todo: Jacobian argument + subroutine solve_LLS(A, b, 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 +! 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 +! w - workspace structure +! options - NLLS_options structure +! pd - logical, true if A is known to be positive definite +! ----------------------------------------------------------------- + implicit none + real(wp), intent(inout), contiguous :: A(:,:), b(:) + integer, intent(in) :: n, m + type(NLLS_inform), intent(inout) :: inform + type(solve_LLS_work), intent(inout) :: w + 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 (n == m) then + call solve_gesv(A,b,n,inform,w) + else + call solve_gels(A,b,n,m,inform,w,options) + 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) + 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 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 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) + end if + + end subroutine solve_LLS + + subroutine solve_gels(A,b,n,m,inform,w,options) +! Wrapper around LAPACK's ?gels + implicit none + real(wp), intent(inout), contiguous :: A(:,:), b(:) + INTEGER, INTENT(IN) :: n, m + type(NLLS_inform), INTENT(INOUT) :: inform + type(NLLS_options), Intent(In) :: 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 + lda = m + ldb = max(m,n) + 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, 1, A, lda, b, 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_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 + + 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' + 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(INOUT), contiguous :: A + REAL(wp), DIMENSION(:), INTENT(INOUT), contiguous :: b + INTEGER, INTENT(IN) :: n + type(NLLS_inform), INTENT(INOUT) :: inform + + 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' + end if + + end subroutine solve_posv + + 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 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), intent(inout) :: inform + type(LLS_lsqr_work), intent(inout) :: w + type(NLLS_options), intent(in) :: options + real(wp), intent(in), optional, contiguous :: precon(:,:) + logical, intent(in), optional :: init_guess + + integer :: action, use_precon, rows, cols + character(len=1) :: trans, no_trans + real(wp) :: norm_a + type(lsqr_options) :: lsqr_opt + + 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, 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 + w%z = 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)(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 + 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(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)(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)(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 + 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)(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 + inform%external_name = 'lapack_?trtrs' + return + end if + w%v = w%v + w%ATu + case (2) ! compute u = u + Az + 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) + 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%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 + 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 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 + type(solve_LLS_work), Intent(inout) :: w + type(NLLS_options), Intent(In) :: options + type(lsqr_options) :: lsqr_opt + + integer :: i, lwork, action, tsize, s + real(wp) :: norm_a + + lwork = size(w%work) + lsqr_opt%stop_test = 2 + + 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%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, 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)(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' + return + end if + + ! 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 + + ! calculate initial guess from sketch-and-solve: x_0 = R^{-1} Q^T Sb + 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' + 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) + + ! 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 solve_rand + + 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 + + ! 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 + 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)) + + end subroutine estimate_norm + + 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(:,:), Sb(:) + type(NLLS_options), intent(in) :: options + type(LLS_rand_work), intent(inout) :: w + type(NLLS_inform), intent(inout) :: inform + + select case (options%sketch_method) + case (1) ! random sampling of rows + call select_random_rows(A, b, n, sketch_size, SA, Sb, w) + case (2) ! random projection using DCT + 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, 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(:,:), Sb(:) + type(LLS_rand_work), intent(inout) :: w + + integer :: i, idx + + 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_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, 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(:,:) + 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 + 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, :) + w%temp_1(i) = -b(i) + else + w%DCT_A(i, :) = A(i, :) + w%temp_1(i) = b(i) + end if + end do + + ! apply discrete cosine transform to each column + 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)) + 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 + ! 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 + +end module MODULE_PREC(ral_nlls_linear) diff --git a/libRALFit/src/ral_nlls_workspaces.F90 b/libRALFit/src/ral_nlls_workspaces.F90 index d433f79d..3b090ba8 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 @@ -48,6 +49,9 @@ 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 + 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 @@ -143,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) @@ -266,7 +265,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. @@ -285,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 (sketch-and-precondition) + + INTEGER :: lls_solver = 1 + +! If using randomised linear solver, what sketch size? +! this must be manually set by the user as it depends on `m` + + INTEGER :: sketch_size = -1 + +! 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 !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -593,20 +615,34 @@ 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(:) end type evaluate_model_work + type, public :: LLS_lsqr_work + logical :: allocated = .false. + ! variables for LSQR + real(wp), allocatable :: u(:), v(:), z(:), x(:), ATu(:) + + ! 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 :: temp_1(:), temp_2(:), Sb(:) + 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 logical :: allocated = .false. - real(wp), allocatable :: temp(:), work(:), Jlls(:) + 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 @@ -649,13 +685,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 +701,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(:,:) @@ -734,6 +771,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 @@ -1258,7 +1296,13 @@ 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), 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 @@ -1277,10 +1321,87 @@ 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 ) + 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), 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 + 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) + if(allocated(w%ATu)) deallocate(w%ATu, 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 + + 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(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" + 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%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%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 + subroutine setup_workspace_evaluate_model(n,m,w,options,inform) implicit none integer, intent(in) :: n, m @@ -1337,7 +1458,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 +1497,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 +1681,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..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) @@ -1424,7 +1455,22 @@ 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_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 + + 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 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) @@ -1448,12 +1494,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 cdaace7d..576522de 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 @@ -1411,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 @@ -1434,6 +1435,149 @@ 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(4), F_results(4) + 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, 4 + ! we use 4 for solver 3 with sketch method 2 + call generate_data_example(params) + 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 + 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) @@ -1767,7 +1911,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 @@ -1781,7 +1925,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 @@ -1813,6 +1957,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 @@ -1934,10 +2079,11 @@ 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 - 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 @@ -2048,7 +2194,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 @@ -2149,12 +2295,37 @@ subroutine all_eig_symm_tests(options,fails) end subroutine all_eig_symm_tests - subroutine solve_LLS_tests(options,fails) + 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 integer, intent(out) :: fails - real(wp), allocatable :: J(:), f(:), d(:), Jd(:) + real(wp), allocatable :: J(:,:), J_copy(:,:), f(:), d(:), Jd(:), JT(:,:), JT_copy(:,:) real(wp) :: normerror integer :: n,m type( nlls_workspace ) :: work @@ -2165,20 +2336,37 @@ subroutine solve_LLS_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 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 ] - 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.) + 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 + 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 + + 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' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2186,13 +2374,14 @@ subroutine solve_LLS_tests(options,fails) end if ! check answer call mult_J(J,n,m,d,Jd,.True.) - normerror = norm2(Jd + f) - if ( normerror > 1.0e-12_wp ) then + normerror = norm2(Jd - f) + 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 fails = fails + 1 end if + status%status = 0 ! J = [ 1.0_wp, 6.0_wp, ! 2.0_wp, 7.0_wp, @@ -2200,41 +2389,59 @@ 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 ] - 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 + + JT_copy = JT + d = f + + options%fortran_jacobian = .false. + call solve_LLS(JT_copy,d,n,m,status, & + work%calculate_step_ws%dogleg_ws%solve_lls_ws,options,.false.) - call solve_LLS(J,f,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.False.) 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) - if ( normerror > 1.0e-12_wp ) then + normerror = norm2(Jd - f) + 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 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,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.True.) + 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 write (*, *) '[id:3] solve_LLS test failed: wrong error message returned' write(*,*) 'status = ', status%status, " (expected ",NLLS_ERROR_FROM_EXTERNAL,")" @@ -2242,20 +2449,324 @@ 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,n,m,d,status, & - work%calculate_step_ws%dogleg_ws%solve_LLS_ws,.True.) + 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' 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 + options%lls_solver = 1 ! LAPACK + + 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, 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 (.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 + end if + + call reset_default_options(options) + + 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%lls_solver = 1 ! LAPACK + 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(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 (.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 + end if + + A = 0.0_wp + b(1) = 6.0_wp + b(2) = 3.0_wp + 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 + fails = fails + 1 + end if + + call nlls_finalize(work,options,status) + call reset_default_options(options) + + 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 ), 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 + + 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,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,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-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) + end subroutine solve_LLS_lsqr_tests - end subroutine solve_LLS_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 ), 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 + + 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) + ! 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,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-12_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 + end do + + deallocate(J, J_copy, f, d, Jd) + call nlls_finalize(work,options,inform) + call reset_default_options(options) + + end subroutine solve_LLS_randomised_tests subroutine findbeta_tests(options,fails) @@ -2279,7 +2790,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 @@ -2323,7 +2834,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 @@ -2332,7 +2843,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 @@ -2342,7 +2853,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 @@ -2382,7 +2893,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 @@ -2392,7 +2903,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 @@ -2403,7 +2914,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 @@ -2414,7 +2925,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 @@ -2439,7 +2950,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 @@ -2448,7 +2959,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 @@ -2457,7 +2968,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 @@ -2528,7 +3039,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 @@ -2556,7 +3067,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 @@ -2598,126 +3109,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 @@ -2744,7 +3135,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 @@ -2783,7 +3174,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 @@ -2821,7 +3212,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 @@ -2893,7 +3284,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 @@ -2908,7 +3299,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 @@ -2981,7 +3372,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 @@ -3025,7 +3416,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 @@ -3099,12 +3490,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