Merge pull request ralna#51 from ralna/add_fallbacks

Add fallbacks

libRALFit/include/ral_nlls.h

@@ -35,6 +35,7 @@ struct ral_nlls_options_d {
   int model; /* what model to use? */
   int type_of_method; /* what method to use? */
   int nlls_method; /* what nlls method to use? */
+  bool allow_fallback_method; /* switch nlls method if chosen one fails? */
   int lls_solver; /* which lls solver to use? */
   ral_nllspkgtype_d_ stop_g_absolute; /* absolute stopping tolerance for gradient*/
   ral_nllspkgtype_d_ stop_g_relative; /* relative stopping tolerance for gradient*/

libRALFit/src/ral_nlls_ciface.f90

@@ -23,6 +23,7 @@ module ral_nlls_ciface
      integer(C_INT) :: model
      integer(C_INT) :: type_of_method
      integer(C_INT) :: nlls_method
+     logical(c_bool) :: allow_fallback_method
      integer(C_INT) :: lls_solver
      real(wp) :: stop_g_absolute
      real(wp) :: stop_g_relative
@@ -151,6 +152,7 @@ subroutine copy_options_in(coptions, foptions, f_arrays)
     foptions%model = coptions%model
     foptions%type_of_method = coptions%type_of_method
     foptions%nlls_method = coptions%nlls_method
+    foptions%allow_fallback_method = coptions%allow_fallback_method
     foptions%lls_solver = coptions%lls_solver
     foptions%stop_g_absolute = coptions%stop_g_absolute
     foptions%stop_g_relative = coptions%stop_g_relative
@@ -299,6 +301,7 @@ subroutine ral_nlls_default_options_d(coptions) bind(C)
   coptions%model = foptions%model
   coptions%type_of_method = foptions%type_of_method
   coptions%nlls_method = foptions%nlls_method
+  coptions%allow_fallback_method = foptions%allow_fallback_method
   coptions%lls_solver = foptions%lls_solver
   coptions%stop_g_absolute = foptions%stop_g_absolute
   coptions%stop_g_relative = foptions%stop_g_relative

libRALFit/src/ral_nlls_internal.f90

@@ -811,11 +811,11 @@ recursive subroutine nlls_iterate(n, m, X,                   &
       End If
      Else If (buildmsg(3, .False., options).And.it_type=='R') Then
 !     Level 3: Same as level 2 but long banner version
-      If (mod(inform%iter, options%print_header)==0) Then
-        Write(rec(1), Fmt=8000)
-        Write(rec(2), Fmt=9000)
-        Write(rec(3), Fmt=8000)
-        Call printmsg(3, .False., options, 3, rec)
+        If (mod(w%iter, options%print_header)==0) Then
+           Write(rec(1), Fmt=8000)
+           Write(rec(2), Fmt=9000)
+           Write(rec(3), Fmt=8000)
+           Call printmsg(3, .False., options, 3, rec)
       End If
     End If
     If (buildmsg(2, .True., options).And.it_type == 'R') Then
@@ -907,7 +907,7 @@ end subroutine nlls_iterate
 
   subroutine nlls_finalize(w,options)
     implicit none
-    type( nlls_workspace ) :: w
+    type( nlls_workspace ), intent(inout) :: w
     type( nlls_options ) :: options
     ! reset all the scalars
     w%first_call = 1
@@ -1014,7 +1014,7 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ
     real(wp), intent(out) :: norm_2_d,norm_S_d
     TYPE( nlls_options ), INTENT( IN ) :: options
     TYPE( nlls_inform ), INTENT( INOUT ) :: inform
-    TYPE( calculate_step_work ) :: w
+    TYPE( calculate_step_work ), intent(inout) :: w
     Type( tenJ_type ), Intent(InOut) :: tenJ
     Type( NLLS_workspace ), Intent(InOut) :: inner_workspace
 
@@ -1023,7 +1023,9 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ
     logical :: scaling_used
     real(wp) :: normx
     Character(Len=85) :: rec(1)
-
+    integer :: num_methods_tried, subproblem_method
+    logical :: subproblem_success
+
 
     if (.not. w%allocated) then
       inform%status = NLLS_ERROR_WORKSPACE_ERROR
@@ -1136,69 +1138,148 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ
 
 
        ! (Gauss-)/(Quasi-)Newton method -- solve as appropriate...
+       subproblem_success = .false.
+       subproblem_method = options%nlls_method
+       num_methods_tried = 0
 
-       if ( options%type_of_method == 1) then
-          select case (options%nlls_method)
-          case (1) ! Powell's dogleg
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3000) 'dogleg'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call dogleg(J,f,hf,g,n,m,Delta,d,norm_S_d,options,inform,w%dogleg_ws)
-             if (inform%status /= 0) Go To 100
-          case (2) ! The AINT method
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3000) 'AINT_TR'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call AINT_TR(J,w%A,f,X,w%v,hf,n,m,Delta,d,norm_S_d,options,inform,w%AINT_tr_ws)
-             if (inform%status /= 0) Go To 100
-          case (3) ! More-Sorensen
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3000) 'More-Sorensen'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call more_sorensen(w%A,w%v,n,m,Delta,d,norm_S_d,options,inform,w%more_sorensen_ws)
-             if (inform%status /= 0) Go To 100
-          case (4) ! Galahad
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3000) 'Galahad DTRS'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call solve_galahad(w%A,w%v,n,m,Delta,num_successful_steps, &
-                  d,norm_S_d,w%reg_order,options,inform,w%solve_galahad_ws)
-             if (inform%status /= 0) Go To 100
-          case default
-             inform%status = NLLS_ERROR_UNSUPPORTED_METHOD
-             goto 100
-          end select ! nlls_method
-       elseif (options%type_of_method == 2) then
-          select case (options%nlls_method)
-          case (3) ! home-rolled regularization solver
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3020) 'RALFit Solver'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call regularization_solver(w%A,w%v,n,m,Delta,num_successful_steps, &
-                  d,norm_S_d,w%reg_order,options,inform,w%regularization_solver_ws)
-             if (inform%status /= 0) Go To 100
-          case(4) ! Galahad
-             If (buildmsg(5,.False.,options)) Then
-               Write(rec(1), Fmt=3020) 'Galahad DRQS'
-               Call printmsg(5,.False.,options,1,rec)
-             End If
-             call solve_galahad(w%A,w%v,n,m,Delta,num_successful_steps, &
-                  d,norm_S_d,w%reg_order,options,inform,w%solve_galahad_ws)
-             if (inform%status /= 0) Go To 100
-          case default
-             inform%status = NLLS_ERROR_UNSUPPORTED_METHOD
-             goto 100
-          end select ! nlls_method
-       else
-          inform%status = NLLS_ERROR_UNSUPPORTED_TYPE_METHOD
-          goto 100
-       end if ! type_of_method
+       do while (.not. subproblem_success)
+
+          if  (num_methods_tried>0) then
+             ! If we're using a method for the second time, or
+             ! we're not allowing a fallback, then bail...
+             if ( (subproblem_method == options%nlls_method) .or. & 
+                  (.not. options%allow_fallback_method)) Go To 100
+
+             if (buildmsg(5,.False.,options)) then
+                Write(rec(1), Fmt=3040)
+                Call printmsg(5,.False.,options,1,rec)
+             end if
+          end if
+
+          if ( options%type_of_method == 1) then
 
+             select case (subproblem_method)
+             case (1) ! Powell's dogleg
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3000) 'dogleg'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (.not. w%dogleg_ws%allocated) then
+                   call setup_workspace_dogleg(n,m,w%dogleg_ws, &
+                        options, inform)
+                end if
+                subproblem_method = 2 ! try aint next
+
+                call dogleg(J,f,hf,g,n,m,Delta,d,norm_S_d, &
+                     options,inform,w%dogleg_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case (2) ! The AINT method
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3000) 'AINT_TR'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (.not. w%aint_tr_ws%allocated) then
+                   call setup_workspace_aint_tr(n,m,w%AINT_tr_ws, &
+                        options, inform)
+                end if
+                subproblem_method = 3 ! try more-sorensen next
+
+                call AINT_TR(J,w%A,f,X,w%v,hf,n,m,Delta,d,norm_S_d, &
+                     options,inform,w%AINT_tr_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case (3) ! More-Sorensen
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3000) 'More-Sorensen'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (.not. w%more_sorensen_ws%allocated) then
+                   call setup_workspace_more_sorensen(n,m,w%more_sorensen_ws, &
+                        options, inform)
+                end if
+                subproblem_method = 4 ! try gltr next 
+
+                call more_sorensen(w%A,w%v,n,m,Delta,d,norm_S_d, &
+                     options,inform,w%more_sorensen_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case (4) ! Galahad
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3000) 'Galahad DTRS'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (use_second_derivatives) then
+                   subproblem_method = 2 ! try aint next
+                   ! dogleg doesn't work with Hessians...
+                else 
+                   subproblem_method = 1 ! try dogleg next
+                end if
+                if (.not. w%solve_galahad_ws%allocated) then
+                   call setup_workspace_solve_galahad(n,m,w%solve_galahad_ws, &
+                        options, inform)
+                end if
+
+                call solve_galahad(w%A,w%v,n,m,Delta,num_successful_steps, &
+                     d,norm_S_d,w%reg_order,options,inform,w%solve_galahad_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case default
+                ! this should never be reached, as it's checked in
+                ! setup_workspaces and the error should be flagged there
+                inform%status = NLLS_ERROR_UNSUPPORTED_METHOD
+                goto 100
+             end select ! nlls_method
+          elseif (options%type_of_method == 2) then
+
+             select case (options%nlls_method)
+             case (3) ! home-rolled regularization solver
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3020) 'RALFit Solver'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (.not. w%regularization_solver_ws%allocated) then
+                   call setup_workspace_regularization_solver(n,m, &
+                        w%regularization_solver_ws, options, inform)
+                end if               
+                subproblem_method = 4
+
+                call regularization_solver(w%A,w%v,n,m,Delta,      &
+                     num_successful_steps, d,norm_S_d,w%reg_order, &
+                     options,inform,w%regularization_solver_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case(4) ! Galahad
+                If (buildmsg(5,.False.,options)) Then
+                   Write(rec(1), Fmt=3020) 'Galahad DRQS'
+                   Call printmsg(5,.False.,options,1,rec)
+                End If
+                if (.not. w%solve_galahad_ws%allocated) then
+                   call setup_workspace_solve_galahad(n,m, &
+                        w%solve_galahad_ws, options, inform)
+                end if
+                subproblem_method = 3
+
+                call solve_galahad(w%A,w%v,n,m,Delta,num_successful_steps, &
+                     d,norm_S_d,w%reg_order,options,inform,w%solve_galahad_ws)
+                if (inform%status == 0) subproblem_success = .true.
+
+             case default
+                ! this should never be reached, as it's checked in
+                ! setup_workspaces and the error should be flagged there
+                inform%status = NLLS_ERROR_UNSUPPORTED_METHOD
+                goto 100
+             end select ! nlls_method
+
+          else
+             inform%status = NLLS_ERROR_UNSUPPORTED_TYPE_METHOD
+             goto 100
+          end if ! type_of_method
+          if (.not. subproblem_success) then 
+             num_methods_tried = num_methods_tried + 1
+          end if
+       end do
+
         ! reverse the scaling on the step
         if ( (scaling_used) ) then
            do i = 1, n
@@ -1241,6 +1322,7 @@ RECURSIVE SUBROUTINE calculate_step(J,f,hf,g,X,md,md_gn,n,m,use_second_derivativ
 3010 FORMAT('*** Subproblem solution found ***')
 3020 FORMAT('*** Solving the regularized subproblem using ',A,' ***')
 3030 FORMAT('*** Error or Subproblem solution NOT found ***')
+3040 FORMAT('*** Subproblem solution not found.  Trying next method ***')
 8000 Format(85('-'))
 
    END SUBROUTINE calculate_step

libRALFit/src/ral_nlls_pyiface.c

@@ -286,6 +286,22 @@ bool set_opts(struct ral_nlls_options *options, PyObject *pyoptions) {
          options->nlls_method = (int) v;
          continue;
       }
+      // bool: allow_fallback_method
+
+      if(strcmp(key_name, "allow_fallback_method")==0) {
+	int vint = PyObject_IsTrue(value); // 1 if true, 0 otherwise
+	printf("%d\n",vint);
+	if (vint == 1){
+	  options->allow_fallback_method=true;
+	}else if (vint == 0){
+	  options->allow_fallback_method=false;
+	}else{
+	  PyErr_SetString(PyExc_RuntimeError, "options['allow_fallback_method'] must be a bool.");
+	  return false;
+	}
+	continue;
+      }
+
       if(strcmp(key_name, "lls_solver")==0) {
 	long v = PyInt_AsLong(value);
 	if(v==-1 && PyErr_Occurred()) {

libRALFit/src/ral_nlls_workspaces.f90

@@ -118,6 +118,11 @@ module ral_nlls_workspaces
 
      INTEGER :: nlls_method = 4
 
+!   allow the algorithm to use a different subproblem solver if one fails
+
+     LOGICAL :: allow_fallback_method = .true.
+
+
 !  which linear least squares solver should we use?
 
      INTEGER :: lls_solver = 1
@@ -603,6 +608,9 @@ end subroutine eval_hp_type
   end type NLLS_workspace
 
   public :: setup_workspaces, remove_workspaces
+  public :: setup_workspace_dogleg, setup_workspace_AINT_tr
+  public :: setup_workspace_more_sorensen, setup_workspace_solve_galahad
+  public :: setup_workspace_regularization_solver
 
 contains
 
@@ -1010,7 +1018,7 @@ end subroutine setup_workspace_calculate_step
 
   recursive subroutine remove_workspace_calculate_step(w,options,tenJ, inner_workspace)
     implicit none
-    type( calculate_step_work ), intent(out) :: w
+    type( calculate_step_work ), intent(inout) :: w
     type( nlls_options ), intent(in) :: options
     type( tenJ_type), Intent(inout) :: tenJ
     Type( NLLS_workspace), Intent(InOut) :: inner_workspace
@@ -1030,28 +1038,23 @@ recursive subroutine remove_workspace_calculate_step(w,options,tenJ, inner_works
             w%solve_newton_tensor_ws, &
             options, tenJ, inner_workspace)
     else
-       if (options%type_of_method == 1) then
-          select case (options%nlls_method)
-          case (1) ! use the dogleg method
-             call remove_workspace_dogleg(w%dogleg_ws, options)
-          case(2) ! use the AINT method
-             call remove_workspace_AINT_tr(w%AINT_tr_ws, options)
-          case(3) ! More-Sorensen
-             call remove_workspace_more_sorensen(&
-                  w%more_sorensen_ws,options)
-          case (4) ! dtrs (Galahad)
-             call remove_workspace_solve_galahad(&
-                  w%solve_galahad_ws, options)
-          end select
-       elseif (options%type_of_method == 2) then
-          select case (options%nlls_method)
-          case (3) ! regularization_solver
-             call remove_workspace_regularization_solver(&
-                  w%regularization_solver_ws, options)
-          case (4) ! dtrs (Galahad)
-             call remove_workspace_solve_galahad(&
-                  w%solve_galahad_ws, options)
-          end select
+       if (w%dogleg_ws%allocated) then
+          call remove_workspace_dogleg(w%dogleg_ws, options)
+       end if
+       if (w%AINT_tr_ws%allocated) then
+          call remove_workspace_AINT_tr(w%AINT_tr_ws, options)
+       end if
+       if (w%more_sorensen_ws%allocated) then
+          call remove_workspace_more_sorensen(&
+               w%more_sorensen_ws,options)
+       end if
+       if (w%solve_galahad_ws%allocated) then
+          call remove_workspace_solve_galahad(&
+               w%solve_galahad_ws, options)
+       end if
+       if (w%regularization_solver_ws%allocated) then
+          call remove_workspace_regularization_solver(&
+               w%regularization_solver_ws, options)
        end if
     end if
     if (options%scale > 0) call remove_workspace_generate_scaling(w%generate_scaling_ws,options)