!!$=head1 NAME !!$ !!$broyden_mod - on step of broyden_method !!$ !!$=head1 SYNOPSIS !!$ !!$INIT: C !!$ !!$OPERATOR:C !!$ !!$=head1 PARAMETERS !!$ !!$=over 4 !!$ !!$=item forget - logical !!$ !!$ Wheter or not to forget previous step !!$ !!$=item x - C !!$ !!$ Model !!$ !!$=item g - C !!$ !!$ Gradient !!$ !!$back !!$ !!$=head1 DESCRIPTION !!$ !!$One step of brodyen method !!$ !!$=head1 head1 SEE ALSO !!$ !!$L, L !!$ !!=head1 $head1 LIBRARY !!$ !!$B !!$=cut module broyden_mod use ddot_mod use llist_mod implicit none type(list), private :: steps integer, private :: nmem contains subroutine broyden_init (n) integer :: n steps = llist_init () nmem = n end subroutine broyden_init subroutine broyden_close () call llist_close (steps) end subroutine broyden_close function broyden (forget, x, grad) result (stat) integer :: stat real, dimension (:) :: x, grad logical :: forget double precision :: alpha, beta double precision, parameter :: eps = 1.d-20 real, dimension (:), pointer :: g, gi integer :: i, n allocate (g (size (x))) ; g = grad alpha = 0.d0 ; n = llist_depth (steps) call llist_rewind (steps ) do i = 1, n call llist_down (steps , gi , beta ) g = g + alpha * gi alpha = ddot (g, gi) / beta end do g = g / (1. - alpha) alpha = ddot (g, g) if (alpha > eps) then x = x + g if (forget ) then call llist_close (steps) steps = llist_init () forget = .false. deallocate (g) else call llist_add (steps, g, alpha) end if else deallocate (g) end if if ((n > 2 .and. alpha > 100. * beta) .or. n > nmem) then call llist_close (steps) steps = llist_init () forget = .true. end if stat = 0 end function broyden end module broyden_mod