module arnoldi_mod use ddot_mod use llist_mod implicit none contains subroutine arnoldi (grads, g, go) type (list) :: grads real, dimension (:), pointer :: g logical, intent (in) :: go real, dimension (:), pointer :: gi integer :: i, n double precision :: beta double precision, dimension (:), allocatable :: alpha double precision, parameter :: eps = 1.d-20 n = llist_depth (grads) allocate (alpha (n)) call llist_rewind (grads) do i = 1, n call llist_down (grads, gi, beta) alpha (i) = ddot (g, gi) / beta end do call llist_rewind (grads) do i = 1, n call llist_down (grads, gi, beta) g = g - alpha (i) * gi end do deallocate (alpha) beta = ddot (g, g) if (beta > eps .and. go) call llist_add (grads, g, beta) end subroutine arnoldi subroutine arnoldi_diag (grads, res) type (list) :: grads real, dimension (:), intent (out) :: res integer :: i real, dimension (:), pointer :: g double precision :: gn res = 0. call llist_rewind (grads) do i = 1, llist_depth (grads) call llist_down (grads, g, gn) res = res + dprod (g, g) / gn end do end subroutine arnoldi_diag subroutine arnoldi_mat (grads, res) type (list) :: grads real, dimension (:,:), intent (out) :: res integer :: i, j real, dimension (:), pointer :: g double precision :: gn res = 0. call llist_rewind (grads) do i = 1, llist_depth (grads) call llist_down (grads, g, gn) do j = 1, size (g) res (:,j) = res (:,j) + dprod (g (j), g) / gn end do end do end subroutine arnoldi_mat end module arnoldi_mod