Removal of Automatic Arrays

[talassio]

Remove automatic arrays from source

Code has TWO type of triggers for mallocs:

Type one

Calls of the form

Subroutine a(n)
      Integer, Intent(In) :: n
      Real(Kind=wp) :: vector(n)
...

Trigger mallocs and should be removed. A fix could be to propagate workspaces...

Subroutine a(n, vector)
      Integer, Intent(In) :: n
      Real(Kind=wp), Intent(Inout) :: vector(n)
...

Outstanding RALFit mallocs are

• GET_HI automatic array, EI [checkCoding-E-48] Ringfences GALAHAD code
• DRQS_SOLVE automatic array, C_SCALE [checkCoding-E-48]
• DRQS_SOLVE automatic array, H_SCALE [checkCoding-E-48]
• DTRS_SOLVE automatic array, C_SCALE [checkCoding-E-48]
• DTRS_SOLVE automatic array, H_SCALE [checkCoding-E-48]

Other mallocs are related to calling routines with expresion of the form -vector(1:n)

Call routine(n, -1.0_wp * vector(1:n))

mallocs for the result of -1.0_wp * vector(1:n).

Type two

Dummy arguments assumed shape/size translation

j(*) is addressed as a contiguous block of memory.

j(:) is a structure which varies by compiler

The j(:) memory metadata has to be traversed and flattened to contiguity before passing to j(*), which is why the temp is created.

Giving j(:) the Contiguous attribute should help, that would require it to always be passed 'sensible' Allocatable stuff without strides, and not Pointers.

Proposal for get_Hi

Issue partially fixed by #78, still need to address automatic array ei(m)

     Recursive subroutine get_Hi(n, m, X, params, i, Hi, eval_HF, inform, ei, weights)
       Implicit None
       integer, intent(in) :: n, m
       real(wp), intent(in) :: X(:)
       class( params_base_type ), intent(inout) :: params
       integer, intent(in) :: i
       real(wp), intent(out) :: Hi(:,:)
       procedure( eval_hf_type ) :: eval_HF
       type( nlls_inform ), intent( inout ) :: inform
       real( wp ), dimension( m ), intent( in ), optional :: weights
       real( wp ), Intent(Inout) :: ei( m ) ! <--- Avoid automatic arrays, trigger malloc on demand

       ei(:) = 0.0_wp
       If ( present(weights) ) Then
        ei(i) = weights(i)
       Else
        ei(i) = 1.0_wp
       End If

!      call eval_HF(inform%external_return, n, m, X, weights(1:m)*ei, Hi, params) !<---- avoid trigger to malloc 
!      for the storage of weights(:)*ei wich is double waste of space and work.
!      We should inform the user NOT do to the HP product blindly, first check non-zero entries of ei, which 
!      by construction here has only ONE non-zero.

       call eval_HF(inform%external_return, n, m, X, ei, Hi, params)

       inform%h_eval = inform%h_eval + 1
       If (inform%external_return/=0) Then
         inform%status = NLLS_ERROR_EVALUATION
         inform%external_name = "eval_HF"
         inform%external_return = 2513
       End iF
     end subroutine get_Hi

[talassio]

@tyronerees is already updating this

[talassio]

Updated list of automatic array allocation

Type one fixed by #78

• [x] line 293 in file nlls_iterate * OBVIOUS mult weight see #39
• [x] line 346 in file nlls_iterate * OBVIOUS mult weight see #39

Type one

• [x] get_hi has an automatic array, EI * OBVIOUS see top comment in #76
• [x] line 131 in file aint_tr * OBVIOUS b=-v
• [x] line 137 in file aint_tr * OBVIOUS b=-v
• [x] line 158 in file aint_tr * OBVIOUS a=a+lam*w%b, b=-v
• [x] line 164 in file aint_tr * OBVIOUS a=a+lam*w%b, b=-v
• [x] line 74 in file aint_tr * OBVIOUS b=-v
• [x] line 79 in file aint_tr * OBVIOUS b=-v
• [x] line 125 in file more_sorensen_noew * OBVIOUS b=-v
• [x] line 232 in file more_sorensen_noew * OBVIOUS matmul(w%aplussigma,d)
• [x] line 313 in file more_sorensen_noew * OBVIOUS b=-v
• [x] line 18 in file matrix_norm * OBVIOUS matmul(a,x)
• [x] line 40 in file get_hi * OBVIOUS weights(1:m)*ei
• [x] line 43 in file apply_second_order_info * OBVIOUS weights(1:m)*w%f
• [x] line 57 in file evaluate_model * OBVIOUS x=f(1:m)+w%jd(1:m) NEXT one to FIX (easy)
• [x] line 40 in file regularization_solver * OBVIOUS b=-v
• [x] line 51 in file linpack_method * OBVIOUS dot_product(u,matmul(h,u))
• [x] line 49 in file check_shift_and_solve * OBVIOUS Call solve_spd_nocopy (a=w%aplussigma,b=-v
• [x] line 76 in file dogleg * OBVIOUStwo_norm(x=alpha*w%d_sd)
• [x] line 82 in file dogleg * OBVIOUS w%d_sd = alpha*w%d_sd (trivial, replace by loop)
• [x] line 72 in file more_sorensen_ew * OBVIOUS Callsolve_spd_nocopy(a=w%aplussigma,b=-v,
• [x] line 263 in file more_sorensen_ew * OBVIOUS Callsolve_spd_nocopy(a=w%aplussigma,b=-v
• [x] line 48 in file evaltensor_j # NOT OBVIOUS j((jj-1)*m+1:(jj-1)*m+params%m) = params%j((jj-1)*params%m+1:jj* replace by loop (stride=1)
• [x] line 55 in file evaltensor_f # NOT OBVIOUS f(1:params%m) = params%f(1:params%m) + params%tenj%js(1:params%m) replace by loop (stride=1)
• [x] line 59 in file evaltensor_f # NOT OBVIOUS f(params%m+1:params%m+n) = replace by loop (stride=1)

Ring-fenced GALAHAD code

• [x] line 146 in file drqs_solve_main * OBVIOUS maxval(abs(... fixable using a for loop
• [x] line 134 in file dtrs_solve_main * OBVIOUS maxval(abs(h(:n))) fixable using a for loop
• [x] line 69 in file dtrs_solve * OBVIOUS scale_h = maxval(abs(h)) replace by loop
• [x] line 83 in file dtrs_solve * OBVIOUS scale_c = maxval(abs(c)) replace by loop
• [x] line 70 in file drqs_solve * OBVIOUS scale_h = maxval(abs(h)) replace by loop
• [x] line 84 in file drqs_solve * OBVIOUS scale_c = maxval(abs(c)) replace by loop
• [x] line 459 in file dtrs_solve_main # NOT OBVIOUS maxloc?
• [x] line 765 in file drqs_solve_main # NOT OBVIOUS maxloc?

Type two

• [x] line 111 in file solve_newton_tensor (allocation: assumed-shape to assumed-size conversion)
• [x] line 143 in file solve_newton_tensor (allocation: assumed-shape to assumed-size conversion)
• [x] line 124 in file aint_tr (allocation: assumed-shape to assumed-size conversion)
• [x] line 125 in file aint_tr (allocation: assumed-shape to assumed-size conversion)
• [x] line 124 in file calculate_step (allocation: assumed-shape to assumed-size conversion)
• [x] line 130 in file calculate_step (allocation: assumed-shape to assumed-size conversion)
• [x] line 171 in file calculate_step (allocation: assumed-shape to assumed-size conversion)
• [x] line 27 in file two_norm (allocation: assumed-shape to assumed-size conversion)
• [x] line 77 in file more_sorensen_noew (allocation: assumed-shape to assumed-size conversion)
• [x] line 78 in file more_sorensen_noew (allocation: assumed-shape to assumed-size conversion)
• [x] line 82 in file more_sorensen_noew (allocation: assumed-shape to assumed-size conversion)
• [x] line 21 in file outer_product (allocation: assumed-shape to assumed-size conversion)
• [x] line 31 in file solve_general (allocation: assumed-shape to assumed-size conversion)
• [x] line 33 in file rank_one_update (allocation: assumed-shape to assumed-size conversion)
• [x] line 48 in file rank_one_update (allocation: assumed-shape to assumed-size conversion)
• [x] line 51 in file rank_one_update (allocation: assumed-shape to assumed-size conversion)
• [x] line 55 in file rank_one_update (allocation: assumed-shape to assumed-size conversion)
• [x] line 38 in file all_eig_symm (allocation: assumed-shape to assumed-size conversion)
• [x] line 48 in file apply_second_order_info (allocation: assumed-shape to assumed-size conversion)
• [x] line 44 in file get_hi (allocation: assumed-shape to assumed-size conversion)
• [x] line 43 in file reset_gradients (allocation: assumed-shape to assumed-size conversion)
• [x] line 55 in file reset_gradients (allocation: assumed-shape to assumed-size conversion)
• [x] line 64 in file reset_gradients (allocation: assumed-shape to assumed-size conversion)
• [x] line 56 in file evaluate_model (allocation: assumed-shape to assumed-size conversion)
• [x] line 78 in file evaluate_model (allocation: assumed-shape to assumed-size conversion)
• [x] line 47 in file linpack_method (allocation: assumed-shape to assumed-size conversion)
• [x] line 141 in file solve_galahad (allocation: assumed-shape to assumed-size conversion)
• [x] line 86 in file solve_galahad (allocation: assumed-shape to assumed-size conversion)
• [x] line 24 in file solve_spd_nocopy (allocation: assumed-shape to assumed-size conversion)
• [x] line 24 in file solve_spd (allocation: assumed-shape to assumed-size conversion)
• [x] line 56 in file max_eig (allocation: assumed-shape to assumed-size conversion)
• [x] line 66 in file min_eig_symm (allocation: assumed-shape to assumed-size conversion)
• [x] line 56 in file dogleg (allocation: assumed-shape to assumed-size conversion)

[talassio]

This will be closed by the box_logic pull request.