Adds LSQR and randomised solvers

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

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. 

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 ################

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)

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)