Re: Help with translation: Fortran77 => Ada95

[Gary Scott <scottg@flash.net>]

Hi,

I guess I wasn't intending to start a competition to see how few
statements this could be written in.  I simply noted that the non-blank,
non-comment portion of the posted example was about twice as long as the
original, no big deal.  I also noted that significant improvements could
be made to the original.  Here's a minimalist update which only
eliminates the GOTOs.  At least it's more readable to me, if no one else
(not 100% guaranteed, certainly there's still lots of room for
improvement, a few opportunities for  whole or partial array
operations):

!
!        Input:  n - - the number of points in the Gaussian quadrature
!                      formula; type integer
!                alpha,beta - - arrays of dimension  n  to be filled
!                      with the values of  alpha(k-1), beta(k-1), k=1,2,
!                      ...,n
!                eps - the relative accuracy desired in the nodes
!                      and weights
!
!        Output: zero- array of dimension  n  containing the Gaussian
!                      nodes (in increasing order)  zero(k)=x(k), k=1,2,
!                      ...,n
!                weight - array of dimension  n  containing the
!                      Gaussian weights  weight(k)=w(k), k=1,2,...,n
!                ierr- an error flag equal to  0  on normal return,
!                      equal to  i  if the QR algorithm does not
!                      converge within 30 iterations on evaluating the
!                      i-th eigenvalue, equal to  -1  if  n  is not in
!                      range, and equal to  -2  if one of the beta's is
!                      negative.
!
! The array  e  is needed for working space.
!
subroutine gauss(n, alpha, beta, eps, zero, weight, ierr)

   implicit none

   integer, intent(in)    :: n                    !Automatic extent
   integer, intent(out)   :: ierr                 !Return code
   integer                :: ii, i, j, k, l, m    !Loop counters

   real, intent(in)       :: alpha(n), beta(n), eps
   real, intent(out)      :: zero(n), weight(n)
   real                   :: e(n), b, c, f, g, p, r, s, mml 
!Intermediates and workspaces

   if (n < 1) then  !Invalid extent
      ierr = -1
      return
   end if

   ierr = 0; zero = 0.0; weight = 0.0; e = 0.0 !Initialize

   zero(1) = alpha(1)

   if (beta(1) < 0.0) then
      ierr = -2
      return
   end if

   weight(1) = beta(1)

   if (n == 1) return

   weight(1) = 1.0

   do k = 2, n
      zero(k) = alpha(k)

      if (beta(k) < 0.0) then
        ierr = -2
        return
      end if

      e(k - 1) = sqrt( beta(k) )
      weight(k) = 0.0
   end do

   j = 0

   do l = 1, n

      do m = l, n           ! Look for a small subdiagonal element
         if (m == n) exit
         if (abs( e(m) ) <= eps * (abs( zero(m) ) + abs( zero(m + 1) )))
exit
      end do

      p = zero(l)

      if (m == l) then
         j = 0
         cycle
      end if

      if (j == 30) then !No convergence after 30 interations
         ierr = 1
         return
      end if

      j = j + 1
!
! Form shift
!
      g = (zero(l + 1) - p) / (2.0 * e(l))

      r = sqrt(g * g + 1.0)
      g = zero(m) - p + e(l) / (g + sign(r,g))

      s = 1.0
      c = 1.0
      p = 0.0

      mml = m - l
!
! For i=m-1 step -1 until l do ...
!
      do ii = 1, mml

         i = m - ii
         f = s * e(i)
         b = c * e(i)

         if (abs(f) < abs(g)) then
            s = f / g
            r = sqrt(s * s + 1.0)
            e(i + 1) = g * r
            c = 1.0 / r
            s = s * c
         else
            c = g / f
            r = sqrt(c * c + 1.0)
            e(i + 1) = f * r
            s = 1.0 / r
            c = c * s
         end if

         g = zero(i + 1) - p
         r = (zero(i) - g) * s + 2.0 * c * b
         p = s * r
         zero(i + 1) = g + p
         g = c * r - b
!
! Form first component of vector.
!
         f = weight(i + 1)
         weight(i + 1) = s * weight(i) + c * f
         weight(i) = c * weight(i) - s * f
      end do

      zero(l) = zero(l) - p
      e(l) = g
      e(m) = 0.0

   end do
!
! Order eigenvalues and eigenvectors.
!
   do ii = 2, n
      i = ii - 1
      k = i
      p = zero(i)

      do j = ii, n
         if (zero(j) >= p) cycle
         k = j
         p = zero(j)
      end do

      if (k >= i) cycle

      zero(k) = zero(i)
      zero(i) = p
      p = weight(i)
      weight(i) = weight(k)
      weight(k) = p
   end do

   do k = 1, n
      weight(k) = beta(1) * weight(k) * weight(k)
   end do

   return

end subroutine




tmoran@acm.org wrote:
> 
> >Interesting that it's almost exactly twice as long as the Fortran.
>   Would you prefer this version, shorter than the Fortran one?  I dropped
> the comments that just showed the original Fortran code, tightened up
> things that fit nicely on one line and didn't really need multiple lines,
> pretty printed, and did a few other "personal preference" changes.
> 

<snip>



-- 

Gary Scott
mailto:scottg@flash.net

mailto:webmaster@fortranlib.com
http://www.fortranlib.com

Support the GNU Fortran G95 Project:  http://g95.sourceforge.net