SUBROUTINE kaiser(a, nrows, n, eigenv, trace, sume, ier)

!  EIGENVALUES AND VECTORS OF A SYMMETRIC +VE DEFINITE MATRIX,
!  USING KAISER'S METHOD.
!  REFERENCE: KAISER,H.F. 'THE JK METHOD: A PROCEDURE FOR FINDING THE
!  EIGENVALUES OF A REAL SYMMETRIC MATRIX', COMPUT.J., VOL.15, 271-273, 1972.

!  ARGUMENTS:-
!  A       = INPUT, AN ARRAY CONTAINING THE MATRIX
!            OUTPUT, THE COLUMNS OF A CONTAIN THE NORMALIZED EIGENVECTORS
!            OF A.   N.B. A IS OVERWRITTEN !
!  NROWS   = INPUT, THE FIRST DIMENSION OF A IN THE CALLING PROGRAM.
!  N       = INPUT, THE ORDER OF A, I.E. NO. OF COLUMNS.
!            N MUST BE <= NROWS.
!  EIGENV()= OUTPUT, A VECTOR CONTAINING THE ORDERED EIGENVALUES.
!  TRACE   = OUTPUT, THE TRACE OF THE INPUT MATRIX.
!  SUME    = OUTPUT, THE SUM OF THE EIGENVALUES COMPUTED.
!            N.B. ANY SYMMETRIC MATRIX MAY BE INPUT, BUT IF IT IS NOT +VE
!            DEFINITE, THE ABSOLUTE VALUES OF THE EIGENVALUES WILL BE FOUND.
!            IF TRACE = SUME, THEN ALL OF THE EIGENVALUES ARE POSITIVE
!            OR ZERO.   IF SUME > TRACE, THE DIFFERENCE IS TWICE THE SUM OF
!            THE EIGENVALUES WHICH HAVE BEEN GIVEN THE WRONG SIGNS !
!  IER     = OUTPUT, ERROR INDICATOR
!             = 0 NO ERROR
!             = 1 N > NROWS OR N < 1
!             = 2 FAILED TO CONVERGE IN 10 ITERATIONS

!  LATEST REVISION - 6 September 1990
!  Fortran 90 version - 20 November 1998

!*************************************************************************

IMPLICIT NONE
INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

REAL (dp), INTENT(IN OUT) :: a(:,:)
INTEGER, INTENT(IN)       :: nrows
INTEGER, INTENT(IN)       :: n
REAL (dp), INTENT(OUT)    :: eigenv(:)
REAL (dp), INTENT(OUT)    :: trace
REAL (dp), INTENT(OUT)    :: sume
INTEGER, INTENT(OUT)      :: ier

! Local variables

REAL (dp), PARAMETER :: small = 1.0e-12_dp, zero = 0.0_dp, half = 0.5_dp,  &
                        one = 1.0_dp
INTEGER   :: i, iter, j, k, ncount, nn
REAL (dp) :: absp, absq, COS, ctn, eps, halfp, p, q, SIN, ss, TAN, temp, xj, xk

!   CALCULATE CONVERGENCE TOLERANCE, EPS.
!   CALCULATE TRACE.   INITIAL SETTINGS.

ier = 1
IF(n < 1 .OR. n > nrows) RETURN
ier = 0
iter = 0
trace = zero
ss = zero
DO j = 1,n
  trace = trace + a(j,j)
  DO i = 1,n
    ss = ss + a(i,j)**2
  END DO
END DO
sume = zero
eps = small*ss/n
nn = n*(n-1)/2
ncount = nn

!   ORTHOGONALIZE PAIRS OF COLUMNS J & K, K > J.

20 DO j = 1,n-1
  DO k = j+1,n
    
!   CALCULATE PLANAR ROTATION REQUIRED
    
    halfp = zero
    q = zero
    DO i = 1,n
      xj = a(i,j)
      xk = a(i,k)
      halfp = halfp + xj*xk
      q = q + (xj+xk) * (xj-xk)
    END DO
    p = halfp + halfp
    absp = ABS(p)
    
!   If P is very small, the vectors are almost orthogonal.
!   Skip the rotation if Q >= 0 (correct ordering).
    
    IF (absp < eps .AND. q >= zero) THEN
      ncount = ncount - 1
      IF (ncount <= 0) GO TO 160
      CYCLE
    END IF
    
!   Rotation needed.
    
    absq = ABS(q)
    IF(absp <= absq) THEN
      TAN = absp/absq
      COS = one/SQRT(one + TAN*TAN)
      SIN = TAN*COS
    ELSE
      ctn = absq/absp
      SIN = one/SQRT(one + ctn*ctn)
      COS = ctn*SIN
    END IF
    COS = SQRT((one + COS)*half)
    SIN = SIN/(COS + COS)
    IF(q < zero) THEN
      temp = COS
      COS = SIN
      SIN = temp
    END IF
    IF(p < zero) SIN = -SIN
    
!   PERFORM ROTATION
    
    DO i = 1,n
      temp = a(i,j)
      a(i,j) = temp*COS + a(i,k)*SIN
      a(i,k) = -temp*SIN + a(i,k)*COS
    END DO
  END DO
END DO
ncount = nn
iter = iter + 1
IF(iter < 10) GO TO 20
ier = 2

!   CONVERGED, OR GAVE UP AFTER 10 ITERATIONS

160 DO j = 1,n
  temp = SUM( a(1:n,j)**2 )
  eigenv(j) = SQRT(temp)
  sume = sume + eigenv(j)
END DO

!   SCALE COLUMNS TO HAVE UNIT LENGTH

DO j = 1,n
  IF (eigenv(j) > zero) THEN
    temp = one/eigenv(j)
  ELSE
    temp = zero
  END IF
  a(1:n,j) = a(1:n,j)*temp
END DO

RETURN
END SUBROUTINE kaiser