isamax Function

  function isamax(n,sx,incx)

    !! Routine to finds the index of the array element having
    !! the maximum absolute value
    !! author: Jack Dongarra, Linpack
    !! author: P J Knight, CCFE, Culham Science Centre
    !! n        : input integer : order of the matrix sx
    !! sx(n*incx) : input real array : array being checked
    !! incx     : input integer : interval in storage between sx array elements
    !! This routine finds the array element with the maximum
    !! absolute value, and returns the element index.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    integer :: isamax

    !  Arguments

    integer, intent(in) :: n, incx
    real(dp), dimension(n*incx), intent(in) :: sx

    !  Local variables

    integer :: i,ix
    real(dp) :: smax

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    isamax = 0
    if (n < 1) return

    isamax = 1
    if (n == 1) return

    if (incx /= 1) then

       ix = 1
       if (incx < 0) ix = (-n+1)*incx + 1
       smax = abs(sx(ix))
       ix = ix + incx
       do i = 2,n
          if (abs(sx(ix)) > smax) then
             isamax = i
             smax = abs(sx(ix))
          end if
          ix = ix + incx
       end do

    else

       smax = abs(sx(1))
       do i = 2,n
          if (abs(sx(i)) <= smax) cycle
          isamax = i
          smax = abs(sx(i))
       end do

    end if

  end function isamax