!   Die folgende Reihe von Fortran-Routinen   EIGSYS
!                                             REDUC
!                                             TRED2
!                                             TQL2

!   dient zur numerischen Loesung des "erweiterten Eigenwertproblems"

!                      (A - lambda*B) x = 0    :

!       A = symmetrische Matrix                    [n x n]
!       B = symmetrische, positiv-definite Matrix  [n x n]

!       lambda = zu berechnende Eigenwerte
!       x      = zu berechnende Eigenvektoren

!   Die Programme REDUC, TRED2 und TQL2  entstammen einer sehr
!   angesehenen Quelle, naemlich der Programmbibliothek EISPACK 
!   (Smith et al 1974).

!   Die Routine EIGSYS (Sormann 2007) als erstes der folgenden Programme
!   fasst die Teilprogramme zusammen und ermoeglicht eine sehr
!   komfortable Verwendung dieses Programmsystems, da nur mehr
!   das eine Programm EIGSYS aufgerufen werden muss.

!   Der Aufruf von EIGSYS aus einem Fortran-Programm lautet 
!   --------------------------------------------------------
!       call   eigsys(nm,n,amat,bmat,eigwrt,eigvek)
!   --------------------------------------------------------

!   und enthaelt die folgenden Parameter:

!   INPUT:  INTEGER   nm   = Wert des ersten Index in der Speicher-
!                            platzvereinbarung fuer amat und bmat                            
!           INTEGER   n    = Ordnung der Matrizen
!           DOUBLE PRECISION    amat = Matrix A des Eigenwertproblems
!           DOUBLE PRECISION    bmat = Matrix B des Eigenwertproblems

!   OUTPUT: DOUBLE PRECISION    eigwrt  = Vektor mit allen n Eigenwerten
!           DOUBLE PRECISION    eigvek  = [n x n] Matrix, wobei die i-te 
!                               Spalte den i-ten Eigenvektor enthaelt.

!   WICHTIG: Es wird von den Programmen nicht getestet, ob die gegebenen
!            Matrizen A und B wirklich symmetrisch sind; dies wird dem
!            Benutzer ueberlassen.
!            Nur die Eigenschaft der "positiv-Definitheit" der Matrix B
!            wird im Programm REDUC gecheckt: wenn B nicht positiv-definit
!            ist, kommt eine Fehlermitteilung.


!   Was die Funktion der einzelnen Routinen betrifft, verweise ich auf die
!   Fachliteratur bzw. auf die Unterlagen zu meinen Lehrveranstaltungen.
!   Hier nur eine ganz kurze Charakteristik:

!   REDUC   fuehrt das erweiterte symmetrische Eigenwertproblem auf ein
!           symmetrisches Standard-Eigenwertproblem zurueck.

!   TRED2   wandelt die Matrix A mittels einer sog. Householder-Transformation
!           in eine tridiagonale Matrix um.

!   TQL2    berechnet die Eigenwerte und Eigenvektoren einer tridiagonalen
!           Matrix mittels des QL-Algorithmus.

! H. Sormann,  6-6-2007	    




        subroutine eigsys(nm,n,amat,bmat,eigwrt,eigvek)

! H. Sormann   6-6-2007
! last update  6-6-2007

	implicit none
	integer nm,n,i,j,k,ierr
        double precision  amat(nm,n),bmat(nm,n),eigwrt(n),eigvek(nm,n)
        double precision  e(n),summe,z(nm,n)

	call reduc(nm,n,amat,bmat,ierr)
!	if(ierr.eq.0)print *,'OK in REDUC'
	if(ierr.ne.0)print *,'Problem in REDUC: BMAT not pos.def.'

        call tred2(nm,n,amat,eigwrt,e,z)
!   	print *,'BACK FROM TRED2'

        call tql2(nm,n,eigwrt,e,z,ierr)
!        if(ierr.eq.0)print *,'OK in TQL2'
        if(ierr.ne.0)print *,'eigenvalue nr ',ierr,' did not converge'

! Umrechnung der Eigenvektoren (Spalten von Z) in die Eigenvektoren
! des urspruenglichen erweiterten Eigenwertproblems:

        do 70 k=1,n
          eigvek(n,k)=z(n,k)/bmat(n,n)
          do 80 i=n-1,1,-1
            summe=z(i,k)
            do 90 j=i+1,n
              summe=summe - bmat(j,i)*eigvek(j,k)
90          continue
            eigvek(i,k)=summe/bmat(i,i)
80        continue
70      continue

	return
	end



      subroutine reduc(nm,n,a,b,ierr)
!
      integer i,j,k,n,i1,j1,nm,nn,ierr
      double precision a(nm,n),b(nm,n),dl(n)
      double precision x,y
!
!     this subroutine is a translation of the algol procedure reduc1,
!     num. math. 11, 99-110(1968) by martin and wilkinson.
!     handbook for auto. comp., vol.ii-linear algebra, 303-314(1971).
!
!     this subroutine reduces the generalized symmetric eigenproblem
!     ax=(lambda)bx, where b is positive definite, to the standard
!     symmetric eigenproblem using the cholesky factorization of b.
!
!     on input
!
!        nm must be set to the row dimension of two-dimensional
!          array parameters as declared in the calling program
!          dimension statement.
!
!        n is the order of the matrices a and b.  if the cholesky
!          factor l of b is already available, n should be prefixed
!          with a minus sign.
!
!        a and b contain the real symmetric input matrices.  only the
!          full upper triangles of the matrices need be supplied.  if
!          n is negative, the strict lower triangle of b contains,
!          instead, the strict lower triangle of its cholesky factor l.
!
!        dl contains, if n is negative, the diagonal elements of l.
!
!     on output
!
!        a contains in its full lower triangle the full lower triangle
!          of the symmetric matrix derived from the reduction to the
!          standard form.  the strict upper triangle of a is unaltered.
!
!        b contains in its strict lower triangle the strict lower
!          triangle of its cholesky factor l.  the full upper
!          triangle of b is unaltered.
!
!        dl contains the diagonal elements of l.
!
!        ierr is set to
!          zero       for normal return,
!          7*n+1      if b is not positive definite.
!
!     questions and comments should be directed to burton s. garbow,
!     mathematics and computer science div, argonne national laboratory
!
!     this version dated august 1983.
!
!     ------------------------------------------------------------------
!
      ierr = 0
      nn = iabs(n)
!      if (n .lt. 0) go to 100
!     .......... form l in the arrays b and dl ..........
      do 80 i = 1, n
         i1 = i - 1
!
         do 80 j = i, n
            x = b(i,j)
            if (i .eq. 1) go to 40
!
            do 20 k = 1, i1
   20       x = x - b(i,k) * b(j,k)
!
   40       if (j .ne. i) go to 60
            if (x .le. 0.0d0) go to 1000
            y = dsqrt(x)
            dl(i) = y
            go to 80
   60       b(j,i) = x / y
   80 continue
!     .......... form the transpose of the upper triangle of inv(l)*a
!                in the lower triangle of the array a ..........
!  100 continue
       do 200 i = 1, nn
         i1 = i - 1
         y = dl(i)
!
         do 200 j = i, nn
            x = a(i,j)
            if (i .eq. 1) go to 180
!
            do 160 k = 1, i1
  160       x = x - b(i,k) * a(j,k)
!
  180       a(j,i) = x / y
  200 continue
!     .......... pre-multiply by inv(l) and overwrite ..........
      do 300 j = 1, nn
         j1 = j - 1
!
         do 300 i = j, nn
            x = a(i,j)
            if (i .eq. j) go to 240
            i1 = i - 1
!
            do 220 k = j, i1
  220       x = x - a(k,j) * b(i,k)
!
  240       if (j .eq. 1) go to 280
!
            do 260 k = 1, j1
  260       x = x - a(j,k) * b(i,k)
!
  280       a(i,j) = x / dl(i)
  300 continue
!
! Sormann 6-6-2007: beim Return wird als  b  die Cholesky-Matrix L
!                   zurueckgegeben:
            do 400 i=1,n
              do 410 j=i,n
                b(i,j)=0.0
410           continue
              b(i,i)=dl(i)
400         continue

      go to 1001
!     .......... set error -- b is not positive definite ..........
 1000 ierr = 7 * n + 1
 1001 return
      end



      subroutine tred2(nm,n,a,d,e,z)
!
      integer i,j,k,l,n,ii,nm,jp1
      double precision a(nm,n),d(n),e(n),z(nm,n)
      double precision f,g,h,hh,scale
!
!     this subroutine is a translation of the algol procedure tred2,
!     num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
!     handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
!
!     this subroutine reduces a real symmetric matrix to a
!     symmetric tridiagonal matrix using and accumulating
!     orthogonal similarity transformations.
!
!     on input
!
!        nm must be set to the row dimension of two-dimensional
!          array parameters as declared in the calling program
!          dimension statement.
!
!        n is the order of the matrix.
!
!        a contains the real symmetric input matrix.  only the
!          lower triangle of the matrix need be supplied.
!
!     on output
!
!        d contains the diagonal elements of the tridiagonal matrix.
!
!        e contains the subdiagonal elements of the tridiagonal
!          matrix in its last n-1 positions.  e(1) is set to zero.
!
!        z contains the orthogonal transformation matrix
!          produced in the reduction.
!
!        a and z may coincide.  if distinct, a is unaltered.
!
!     questions and comments should be directed to burton s. garbow,
!     mathematics and computer science div, argonne national laboratory
!
!     this version dated august 1983.
!
!     ------------------------------------------------------------------
!
      do 100 i = 1, n
!
         do 80 j = i, n
   80    z(j,i) = a(j,i)
!
         d(i) = a(n,i)
  100 continue
!
      if (n .eq. 1) go to 510
!     .......... for i=n step -1 until 2 do -- ..........
      do 300 ii = 2, n
         i = n + 2 - ii
         l = i - 1
         h = 0.0d0
         scale = 0.0d0
         if (l .lt. 2) go to 130
!     .......... scale row (algol tol then not needed) ..........
         do 120 k = 1, l
  120    scale = scale + dabs(d(k))
!
         if (scale .ne. 0.0d0) go to 140
  130    e(i) = d(l)
!
         do 135 j = 1, l
            d(j) = z(l,j)
            z(i,j) = 0.0d0
            z(j,i) = 0.0d0
  135    continue
!
         go to 290
!
  140    do 150 k = 1, l
            d(k) = d(k) / scale
            h = h + d(k) * d(k)
  150    continue
!
         f = d(l)
         g = -dsign(dsqrt(h),f)
         e(i) = scale * g
         h = h - f * g
         d(l) = f - g
!     .......... form a*u ..........
         do 170 j = 1, l
  170    e(j) = 0.0d0
!
         do 240 j = 1, l
            f = d(j)
            z(j,i) = f
            g = e(j) + z(j,j) * f
            jp1 = j + 1
            if (l .lt. jp1) go to 220
!
            do 200 k = jp1, l
               g = g + z(k,j) * d(k)
               e(k) = e(k) + z(k,j) * f
  200       continue
!
  220       e(j) = g
  240    continue
!     .......... form p ..........
         f = 0.0d0
!
         do 245 j = 1, l
            e(j) = e(j) / h
            f = f + e(j) * d(j)
  245    continue
!
         hh = f / (h + h)
!     .......... form q ..........
         do 250 j = 1, l
  250    e(j) = e(j) - hh * d(j)
!     .......... form reduced a ..........
         do 280 j = 1, l
            f = d(j)
            g = e(j)
!
            do 260 k = j, l
  260       z(k,j) = z(k,j) - f * e(k) - g * d(k)
!
            d(j) = z(l,j)
            z(i,j) = 0.0d0
  280    continue
!
  290    d(i) = h
  300 continue
!     .......... accumulation of transformation matrices ..........
      do 500 i = 2, n
         l = i - 1
         z(n,l) = z(l,l)
         z(l,l) = 1.0d0
         h = d(i)
         if (h .eq. 0.0d0) go to 380
!
         do 330 k = 1, l
  330    d(k) = z(k,i) / h
!
         do 360 j = 1, l
            g = 0.0d0
!
            do 340 k = 1, l
  340       g = g + z(k,i) * z(k,j)
!
            do 360 k = 1, l
               z(k,j) = z(k,j) - g * d(k)
  360    continue
!
  380    do 400 k = 1, l
  400    z(k,i) = 0.0d0
!
  500 continue
!
  510 do 520 i = 1, n
         d(i) = z(n,i)
         z(n,i) = 0.0d0
  520 continue
!
      z(n,n) = 1.0d0
      e(1) = 0.0d0
      return
      end



      subroutine tql2(nm,n,d,e,z,ierr)
!
      integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr
      double precision d(n),e(n),z(nm,n)
      double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2
!     double precision pythag
!
!     this subroutine is a translation of the algol procedure tql2,
!     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and
!     wilkinson.
!     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971).
!
!     this subroutine finds the eigenvalues and eigenvectors
!     of a symmetric tridiagonal matrix by the ql method.
!     the eigenvectors of a full symmetric matrix can also
!     be found if  tred2  has been used to reduce this
!     full matrix to tridiagonal form.
!
!     on input
!
!        nm must be set to the row dimension of two-dimensional
!          array parameters as declared in the calling program
!          dimension statement.
!
!        n is the order of the matrix.
!
!        d contains the diagonal elements of the input matrix.
!
!        e contains the subdiagonal elements of the input matrix
!          in its last n-1 positions.  e(1) is arbitrary.
!
!        z contains the transformation matrix produced in the
!          reduction by  tred2, if performed.  if the eigenvectors
!          of the tridiagonal matrix are desired, z must contain
!          the identity matrix.
!
!      on output
!
!        d contains the eigenvalues in ascending order.  if an
!          error exit is made, the eigenvalues are correct but
!          unordered for indices 1,2,...,ierr-1.
!
!        e has been destroyed.
!
!        z contains orthonormal eigenvectors of the symmetric
!          tridiagonal (or full) matrix.  if an error exit is made,
!          z contains the eigenvectors associated with the stored
!          eigenvalues.
!
!        ierr is set to
!          zero       for normal return,
!          j          if the j-th eigenvalue has not been
!                     determined after 30 iterations.
!
!     calls pythag for  dsqrt(a*a + b*b) .
!
!     questions and comments should be directed to burton s. garbow,
!     mathematics and computer science div, argonne national laboratory
!
!     this version dated august 1983.
!
!     ------------------------------------------------------------------
!
      ierr = 0
      if (n .eq. 1) go to 1001
!
      do 100 i = 2, n
  100 e(i-1) = e(i)
!
      f = 0.0d0
      tst1 = 0.0d0
      e(n) = 0.0d0
!
      do 240 l = 1, n
         j = 0
         h = dabs(d(l)) + dabs(e(l))
         if (tst1 .lt. h) tst1 = h
!     .......... look for small sub-diagonal element ..........
         do 110 m = l, n
            tst2 = tst1 + dabs(e(m))
            if (tst2 .eq. tst1) go to 120
!     .......... e(n) is always zero, so there is no exit
!                through the bottom of the loop ..........
  110    continue
!
  120    if (m .eq. l) go to 220
  130    if (j .eq. 30) go to 1000
         j = j + 1
!     .......... form shift ..........
         l1 = l + 1
         l2 = l1 + 1
         g = d(l)
         p = (d(l1) - g) / (2.0d0 * e(l))
!         r = pythag(p,1.0d0)
         r=sqrt(p*p+1.d0)
         d(l) = e(l) / (p + dsign(r,p))
         d(l1) = e(l) * (p + dsign(r,p))
         dl1 = d(l1)
         h = g - d(l)
         if (l2 .gt. n) go to 145
!
         do 140 i = l2, n
  140    d(i) = d(i) - h
!
  145    f = f + h
!     .......... ql transformation ..........
         p = d(m)
         c = 1.0d0
         c2 = c
         el1 = e(l1)
         s = 0.0d0
         mml = m - l
!     .......... for i=m-1 step -1 until l do -- ..........
         do 200 ii = 1, mml
            c3 = c2
            c2 = c
            s2 = s
            i = m - ii
            g = c * e(i)
            h = c * p
!            r = pythag(p,e(i))
            r=sqrt(p*p+e(i)*e(i))
            e(i+1) = s * r
            s = e(i) / r
            c = p / r
            p = c * d(i) - s * g
            d(i+1) = h + s * (c * g + s * d(i))
!     .......... form vector ..........
            do 180 k = 1, n
               h = z(k,i+1)
               z(k,i+1) = s * z(k,i) + c * h
               z(k,i) = c * z(k,i) - s * h
  180       continue
!
  200    continue
!
         p = -s * s2 * c3 * el1 * e(l) / dl1
         e(l) = s * p
         d(l) = c * p
         tst2 = tst1 + dabs(e(l))
         if (tst2 .gt. tst1) go to 130
  220    d(l) = d(l) + f
  240 continue
!     .......... order eigenvalues and eigenvectors ..........
      do 300 ii = 2, n
         i = ii - 1
         k = i
         p = d(i)
!
         do 260 j = ii, n
            if (d(j) .ge. p) go to 260
            k = j
            p = d(j)
  260    continue
!
         if (k .eq. i) go to 300
         d(k) = d(i)
         d(i) = p
!
         do 280 j = 1, n
            p = z(j,i)
            z(j,i) = z(j,k)
            z(j,k) = p
  280    continue
!
  300 continue
!
      go to 1001
!     .......... set error -- no convergence to an
!                eigenvalue after 30 iterations ..........
 1000 ierr = l
 1001 return
      end
