Main solver.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(nlesolver_type), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout), | dimension(:) | :: | x |
subroutine nlesolver_solver(me,x) implicit none class(nlesolver_type),intent(inout) :: me real(wp),dimension(:),intent(inout) :: x real(wp),dimension(:) ,allocatable :: fvec !! function vector real(wp),dimension(:,:),allocatable :: fjac !! jacobian matrix [dense] real(wp),dimension(:), allocatable :: fjac_sparse !! jacobian matrix [sparse] real(wp),dimension(:), allocatable :: prev_fjac_sparse !! previous jacobian matrix [sparse] real(wp),dimension(:) ,allocatable :: rhs !! linear system right-hand side real(wp),dimension(:) ,allocatable :: p !! search direction real(wp),dimension(:) ,allocatable :: xold !! previous value of `x` real(wp),dimension(:) ,allocatable :: prev_fvec !! previous function vector real(wp),dimension(:,:),allocatable :: prev_fjac !! previous jacobian matrix real(wp),dimension(:,:),allocatable :: delf !! used for Broyden (rank 2 for `matmul`) real(wp),dimension(:,:),allocatable :: delx !! used for Broyden (rank 2 for `matmul`) logical :: user_stop !! user stop button flag integer :: info !! status flag from the [[linear_solver]] integer :: iter !! iteration counter real(wp) :: f !! magnitude of `fvec` real(wp) :: fold !! previous value of `f` integer :: n_uphill !! number of steps taken in the "uphill" direction !! (where `f` is increasing) real(wp) :: delxmag2 !! used for Broyden logical :: recompute_jac !! if using Broyden, and we want to call the user !! jacobian routine instead integer :: broyden_counter !! number of times the broyden update has been used integer :: alloc_stat !! allocation status flag type(lsqr_solver_ez) :: sparse_solver !! sparse LSQR solver class type(lusol_settings) :: lusol_options integer :: i !! counter integer,dimension(:),allocatable :: idx, index_array !! for sparse indexing character(len=10) :: i_str !! string version of `i` for row string real(wp) :: normA, condA, normr, normAr, normx !! for LSMR integer :: itn !! for LSMR if (me%istat<0) return ! class was not initialized properly if (.not. associated(me%func)) then call me%set_status(istat = -10, string = 'Error: function routine is not associated') return end if if (me%sparsity_mode==NLESOLVER_SPARSITY_DENSE .and. .not. associated(me%grad)) then call me%set_status(istat = -11, string = 'Error: gradient routine is not associated') return end if if (me%sparsity_mode>NLESOLVER_SPARSITY_DENSE .and. .not. associated(me%grad_sparse)) then call me%set_status(istat = -11, string = 'Error: gradient routine is not associated') return end if ! initialize: iter = 0 n_uphill = 0 recompute_jac = .false. broyden_counter = 0 ! allocate the arrays: alloc_stat = 0 if (alloc_stat==0) allocate(fvec (me%m) , stat=alloc_stat) if (me%sparsity_mode==NLESOLVER_SPARSITY_DENSE) then ! dense if (alloc_stat==0) allocate(fjac (me%m,me%n) , stat=alloc_stat) else ! sparse if (alloc_stat==0) allocate(fjac_sparse (me%n_nonzeros) , stat=alloc_stat) if (me%use_broyden .and. alloc_stat==0) allocate(prev_fjac_sparse(me%n_nonzeros) , stat=alloc_stat) end if if (alloc_stat==0) allocate(rhs (me%m) , stat=alloc_stat) if (alloc_stat==0) allocate(p (me%n) , stat=alloc_stat) if (alloc_stat==0) allocate(xold (me%n) , stat=alloc_stat) if (me%use_broyden) then ! only need these for broyden: if (alloc_stat==0) allocate(prev_fvec(me%m) , stat=alloc_stat) if (me%sparsity_mode/=NLESOLVER_SPARSITY_DENSE .and. alloc_stat==0) then allocate(prev_fjac(me%m,me%n) , stat=alloc_stat) index_array = [(i, i=1,me%n_nonzeros)] ! an array to index the sparse jacobian elements end if if (alloc_stat==0) allocate(delf (me%m,1) , stat=alloc_stat) if (alloc_stat==0) allocate(delx (me%n,1) , stat=alloc_stat) end if if (alloc_stat/=0) then call me%set_status(istat = -9, string = 'Error: Out of memory') return else me%istat = -998 me%message = 'Unknown error' end if ! evaluate the function: call me%func(x,fvec) f = norm2(fvec) ! check to see if initial guess is a root: if (f <= me%tol) then call me%set_status(istat = 1, string = 'Required accuracy achieved') else ! main iteration loop: do iter = 1,me%max_iter ! Export the current iteration: if (associated(me%export_iteration)) call me%export_iteration(x,fvec,iter) ! Check for user stop: if (associated(me%user_input_check)) then call me%user_input_check(user_stop) if (user_stop) then call me%set_status(istat = 4, string = 'Stopped by the user') exit end if end if if (me%use_broyden .and. .not. recompute_jac) then if (iter==1) then ! always compute Jacobian on the first iteration !call me%grad(x,fjac) select case (me%sparsity_mode) case (NLESOLVER_SPARSITY_DENSE) call me%grad(x,fjac) case default ! sparse modes call me%grad_sparse(x,fjac_sparse) end select broyden_counter = 0 else ! and use Broyden update to estimate Jacobian ! for subsequent iterations. ! note: fvec was computed the last iteration delx(:,1) = x - xold delf(:,1) = fvec - prev_fvec if (me%sparsity_mode==NLESOLVER_SPARSITY_DENSE) then ! dense delxmag2 = dot_product(delx(:,1),delx(:,1)) if (delxmag2 < eps) then call me%set_status(istat = -8, & string = 'Error: Divide by zero when computing Broyden update') exit end if ! Jacobian estimate: ! This is the equation from wikipedia : https://en.wikipedia.org/wiki/Broyden%27s_method fjac = prev_fjac + matmul((delf-matmul(prev_fjac,delx)), transpose(delx)) / delxmag2 ! see also: eqn 4 in Schubert, which is different ... else ! using a sparse option ! Just a row-by-row version of the fjac equation above. ! see L.K. Schubert, 1970 do i = 1, me%m idx = pack(index_array, mask = me%irow==i) ! the nonzero elements in this row if (size(idx)==0) cycle ! no non-zeros in this row associate (dx => delx(me%icol(idx),:)) ! nonzero x vec for this row delxmag2 = dot_product(dx(:,1),dx(:,1)) ! only those x's for this row if (delxmag2 < eps) then write(i_str,'(I10)') i call me%set_status(istat = -8, & string = 'Error: Divide by zero when computing sparse Broyden update for row '//& trim(adjustl(i_str))) exit end if fjac_sparse(idx) = prev_fjac_sparse(idx) + & matmul((delf(i,1)-matmul(prev_fjac_sparse(idx),dx)), transpose(dx)) / delxmag2 end associate end do end if broyden_counter = broyden_counter + 1 end if prev_fvec = fvec if (me%sparsity_mode==NLESOLVER_SPARSITY_DENSE) then prev_fjac = fjac else prev_fjac_sparse = fjac_sparse end if else ! compute the jacobian: select case (me%sparsity_mode) case (NLESOLVER_SPARSITY_DENSE) call me%grad(x,fjac) case default ! sparse call me%grad_sparse(x,fjac_sparse) end select recompute_jac = .false. ! for broyden broyden_counter = 0 end if xold = x fold = f ! compute the search direction p by solving linear system: rhs = -fvec ! RHS of the linear system select case (me%sparsity_mode) case (NLESOLVER_SPARSITY_DENSE) ! use dense solver call linear_solver(me%m,me%n,fjac,rhs,p,info) case (NLESOLVER_SPARSITY_LSQR) ! initialize the LSQR sparse solver call sparse_solver%initialize(me%m,me%n,fjac_sparse,me%irow,me%icol,& atol = me%atol, & btol = me%btol, & conlim = me%conlim, & itnlim = me%itnlim, & nout = me%nout) call sparse_solver%solve(rhs,me%damp,p,info) ! solve the linear system ! check convergence: select case (info) case(4) call me%set_status(istat = -1004, & string = 'LSQR Error: The system appears to be ill-conditioned. istop =', i=info) exit case(5) call me%set_status(istat = -1005, & string = 'LSQR Error: The iteration limit was reached. istop =', i=info) exit case default info = 0 end select case (NLESOLVER_SPARSITY_LUSOL) ! use lusol solver lusol_options%method = me%lusol_method call solve(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info,& settings=lusol_options) case (NLESOLVER_SPARSITY_LSMR) ! use LSMR solver: !me%conlim = 1.0_wp/(100*epsilon(1.0_wp)) call lsmr_ez ( me%m, me%n, me%irow, me%icol, fjac_sparse, rhs, me%damp, & me%atol, me%btol, me%conlim, me%itnlim, me%localSize, me%nout, & p, info, itn, normA, condA, normr, normAr, normx ) ! check convergence: select case (info) case(4) call me%set_status(istat = -1004, & string = 'LSMR Error: The system appears to be ill-conditioned. istop =', i=info) exit case(5) call me%set_status(istat = -1005, & string = 'LSMR Error: The iteration limit was reached. istop =', i=info) exit case default info = 0 end select case (NLESOLVER_SPARSITY_CUSTOM_SPARSE) if (associated(me%custom_solver_sparse)) then ! a user-provided sparse solver: call me%custom_solver_sparse(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info) else call me%set_status(istat = -1006, & string = 'Error: The custom_solver_sparse procedure has not been set.') exit end if case default error stop 'invalid sparsity_mode' end select ! check for errors: if (info /= 0) then call me%set_status(istat = -6, string = 'Error solving linear system. info =', i=info) exit else ! next step, using the specified method: call me%linesearch(xold,p,x,f,fvec,fjac,fjac_sparse) ! keep track of the number of steps in the "uphill" direction: if (f>fold) then n_uphill = n_uphill + 1 else n_uphill = 0 end if ! check for stopping conditions if (f <= me%tol) then call me%set_status(istat = 1, string = 'Required accuracy achieved') exit elseif ( maxval(abs(x-xold)) <= me%tolx ) then call me%set_status(istat = 2, string = 'Solution cannot be improved') exit elseif (iter == me%max_iter) then call me%set_status(istat = 3, string = 'Maximum number of iterations reached') exit elseif (n_uphill > me%n_uphill_max) then call me%set_status(istat = 5, string = 'Too many steps in the uphill direction') exit elseif (me%use_broyden) then ! If delxmag2 is too small when using broyden, just ! call the user-supplied jacobian function to avoid ! a divide by zero on the next step. This should ! normally only happen when the solution is almost converged. if (norm2(x-xold)**2 <= eps) then recompute_jac = .true. else if (me%broyden_update_n>0) then ! Note that we also recompute if we have taken an uphill step if (broyden_counter==me%broyden_update_n .or. n_uphill>0) then ! time to recompute the full jacobian recompute_jac = .true. end if end if end if endif end if end do !end of iterations loop end if !Export the last iteration: iter = iter + 1 if (associated(me%export_iteration)) call me%export_iteration(x,fvec,iter) end subroutine nlesolver_solver