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) use lsqr_module, only: lsqr_solver_ez use lusol_ez_module, only: solve, lusol_settings 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 :: 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 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==1 .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>1 .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==1) 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) 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 (alloc_stat==0) allocate(prev_fjac(me%m,me%n) , stat=alloc_stat) 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) 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 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: fjac = prev_fjac + & matmul((delf-matmul(prev_fjac,delx)),& transpose(delx))/delxmag2 broyden_counter = broyden_counter + 1 end if prev_fjac = fjac prev_fvec = fvec else ! compute the jacobian: select case (me%sparsity_mode) case (1) call me%grad(x,fjac) case (2:) 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 (1) ! use dense solver call linear_solver(me%m,me%n,fjac,rhs,p,info) case (2) ! 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 (3) ! use lusol solver lusol_options%method = me%lusol_method ! 0 => TPP: Threshold Partial Pivoting. ! 1 => TRP: Threshold Rook Pivoting. ! 2 => TCP: Threshold Complete Pivoting. call solve(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info,& settings=lusol_options) 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