
      SUBROUTINE SwapMatVals(x,y)
      REAL x,y,temp
        temp = x
        x = y
        y = temp
      END !SUB swap

      SUBROUTINE CopyMatrixVectors (xin, xout, n, r)
      INCLUDE 'STDHDR.FOR'

      REAL xin(0:maxr,0:maxc),xout(0:maxc)
      INTEGER n,r,i

      DO i = 0, n - 1
         xout(i) = xin(r, i)
      END DO

      END !SUB




      SUBROUTINE Copyto2DVector (xin, xout, n, r)
      INCLUDE 'STDHDR.FOR'
      REAL xout(0:maxr,0:maxc),xin(0:maxc)
      INTEGER n,r,i

      DO i = 0, n - 1
         xout(r, i) = xin(i)
      END DO

      END !SUB



      SUBROUTINE CyclicJacobi (a, n, eigenvalues, v, count, success)
      INCLUDE 'STDHDR.FOR'
      REAL a(0:maxr,0:maxc), eigenvalues(0:maxc)
      REAL v(0:maxr,0:maxc), tolerance, nearzero
      REAL atemp1(0:maxc), atemp2(0: maxc), sumsq
      INTEGER Maxcount, p, q
      INTEGER n,count
      LOGICAL success

      PARAMETER (nearzero = 1E-08)   ! {change to 1.0e-15 FOR 8087 system}
      PARAMETER (tolerance = 1E-08)

      Maxcount = count
      count = 0
      DO p = 0, n - 1
         DO q = 0, n - 1
            v(p, q) = 0.0
         END DO
         v(p, p) = 1.0
      END DO

      DO WHILE (.NOT. success .OR. count .LE. Maxcount)
        count = count + 1
        sumsq = 0.0
        DO p = 0, n - 1
           sumsq = sumsq + a(p, p) ** 2
        END DO
        success = .TRUE.
        DO p = 0, n - 2
           DO q = p + 1, n - 1
              IF (ABS(a(p, q) / sumsq) .GT. tolerance ) THEN
                success = .FALSE.
                CALL CopyMatrixVectors(a, atemp1, n, p)
                CALL CopyMatrixVectors(a, atemp2, n, q)
                CALL JacobiRotation(v, a, p, q, atemp1, atemp2, n)
                CALL Copyto2DVector(atemp1, a, n, p)
                CALL Copyto2DVector(atemp2, a, n, q)
              END IF
          END DO
        END DO
      END DO
      IF (success) THEN
        DO p = 0, n - 1
        eigenvalues(p) = a(p, p)
        END DO
      END IF
      END !SUB cycjacob



      SUBROUTINE JacobiRotation (v, m, p, q, mp, mq, n)
      INCLUDE 'STDHDR.FOR'
      REAL v(0:maxr,0:maxc),m(0:maxr,0:maxc)
      REAL mp(0:maxc), mq(0:maxc), tolerance, nearzero
      REAL g, h, mpp, mqq, mpq, mpj, mqj, tau
      REAL tempr, angle, c, s, cs, cSquared, sSquared
      REAL PiOver4, PiOver2
      INTEGER n,p,q, j

      PARAMETER (nearzero = 1E-08)   ! {change to 1.0e-15 FOR 8087 system}
      PARAMETER (tolerance = 1E-08)

      PiOver4 = pi / 4.0
      PiOver2 = pi / 2.0

      mpp = mp(p)
      mpq = mp(q)
      mqq = mq(q)
      IF (ABS(mpp - mqq) .GT. nearzero) THEN
         angle = ATAN(2.0 * mpq / (mpp - mqq)) / 2.0
        IF (angle .GT. PiOver4) angle = angle - PiOver2
      ELSE
       IF (mpq .GT. 0 ) THEN
           angle = PiOver4
        ELSE
          angle = -PiOver4
        END IF
      END IF
      c = COS(angle)
      s = SIN(angle)
      cSquared = c ** 2
      sSquared = s ** 2
      cs = c * s
      tempr = 2.0 * mpq * cs
      mp(p) = mpp * cSquared + tempr + mqq * sSquared
      mq(q) = mpp * sSquared - tempr + mqq * cSquared
      mp(q) = 0.0
      mq(p) = 0.0
      tau = s / (1.0 + c)
      DO j = 0, n - 1
         IF ( j .NE. p .AND. j .NE. q) THEN
             mpj = mp(j) * c + mq(j) * s
             mqj = mq(j) * c - mp(j) * s
             mp(j) = mpj
             m(j, p) = mpj
             mq(j) = mqj
             m(j, q) = mqj
         END IF
      END DO
      DO j = 0, n - 1
         g = v(j, p)
         h = v(j, q)
         v(j, p) = c * g + s * h
         v(j, q) = -s * g + c * h
      END DO

      END !SUB JacobRot



      SUBROUTINE MatAdd (m1, m2, numrow, numcol, m3)
      INCLUDE 'STDHDR.FOR'
      REAL m1(0:maxr,0:maxc), m2(0:maxr,0:maxc), m3(0:maxr,0:maxc)
      INTEGER numrow,numcol,i,j


      DO i = 0, numrow - 1
         DO j = 0, numcol - 1
            m3(i, j) = m1(i, j) + m2(i, j)
         END DO
      END DO
      END !SUB MatAdd



      SUBROUTINE MatDeter (m1, numrow, det)
      INCLUDE 'STDHDR.FOR'
      REAL m1(0:maxr,0:maxc)
      INTEGER i,j,k,f, lx, nxt,iErr

      REAL bignum, eval, temp1, temp2, rac
      REAL tempmat[ALLOCATABLE](:,:)
      ALLOCATE(tempmat(0:maxc,0:maxc),STAT=iErr)
      CALL CheckMem(iErr)

      DO i = 0, numrow - 1
         DO j = 0, numrow - 1
            tempmat(i, j) = m1(i, j)
         END DO
      END DO

      f = 1
      DO i = 0, numrow - 1
        bignum = 0
        DO k = i, numrow - 1
          eval = ABS(tempmat(k, i))
          IF (eval - bignum .GT. 0) THEN
            bignum = eval
            lx = k
          END IF
        END DO
        IF (i - lx .NE. 0 ) f = -f
        DO j = 0, numrow - 1
          temp1 = tempmat(i, j)
          tempmat(i, j) = tempmat(lx, j)
          tempmat(lx, j) = temp1
        END DO
        pivot = tempmat(i, i)
        nxt = i + 1
        DO j = nxt, numrow - 1
          rac = tempmat(j, i) / pivot
          DO k = i, numrow - 1
            tempmat(j, k) = tempmat(j, k) - rac * tempmat(i, k)
          END DO
        END DO
      END DO
      temp2 = 1
      DO i = 0, numrow - 1
         temp2 = temp2 * tempmat(i, i)
      END DO
      det = temp2 * f
      DEALLOCATE(tempmat,STAT=iErr)
      CALL CheckDealloc(iErr)
      END !SUB       matdet



      SUBROUTINE MatInvert (matdata, numcol, det, invary)
      INCLUDE 'STDHDR.FOR'
      REAL matdata(0:maxr,0:maxc), invary(0:maxr,0:maxc)
      REAL det,piv, t, leval
      INTEGER numcol, i, j, k, l, lerow, lecol, l1, iErr
      INTEGER pivlst[ALLOCATABLE](:,:)
      LOGICAL pivchk[ALLOCATABLE](:)

      ALLOCATE(pivlst(0:maxc,0:1),STAT=iErr)
      CALL CheckMem(iErr)
      ALLOCATE(pivchk(0:maxc),STAT=iErr)
      CALL CheckMem(iErr)

      det = 1
      DO i = 0, numcol - 1
         pivchk(i) = .FALSE.
         DO j = 0, numcol - 1
           invary(i, j) = matdata(i, j)
         END DO
      END DO

      DO i = 0, numcol - 1
        leval = 0.0
        DO j = 0, numcol - 1
           IF (.NOT. pivchk(j)) THEN
             DO k = 0, numcol - 1
               IF (.NOT. pivchk(k)) THEN
                 IF ( ABS(invary(j, k)) .GT. leval ) THEN
                   lerow = j
                   lecol = k
                   leval = ABS(invary(j, k))
                 END IF
               END IF
             END DO
           END IF
        END DO
        pivchk(lecol) = .TRUE.
        pivlst(i, 0) = lerow
        pivlst(i, 1) = lecol
        IF (lerow .NE. lecol) THEN
          det = -det
          DO l = 0, numcol - 1
            call SwapMatVals( invary(lerow, l), invary(lecol, l))
          END DO
        END IF

        piv = invary(lecol, lecol)
        det = det * piv
        IF (det .GT. 1E30 )THEN
          det = 1      ! {prevents numeric overflow}
        END IF
        invary(lecol, lecol) = 1.0
        DO l = 0, numcol - 1
          invary(lecol, l) = invary(lecol, l) / piv
        END DO
        DO l1 = 0, numcol - 1
          IF (l1 .NE. lecol) THEN
            t = invary(l1, lecol)
            invary(l1, lecol) = 0
            DO l = 0, numcol - 1
              invary(l1, l) = invary(l1, l) - invary(lecol, l) * t
            END DO
          END IF
        END DO
      END DO

      DO i = 0, numcol - 1
        l = numcol - i - 1
        IF (pivlst(l, 0) .NE. pivlst(l, 1)) THEN
          lerow = pivlst(l, 0)
          lecol = pivlst(l, 1)
          DO k = 0, numcol-1
            call SWAPMatVals( invary(k, lerow), invary(k, lecol))
          END DO
       END IF
      END DO
      DEALLOCATE(pivlst,STAT=iErr)
      CALL CheckDealloc(iErr)
      DEALLOCATE(pivchk,STAT=iErr)
      CALL CheckDealloc(iErr)
      END !SUB MatInvert



      SUBROUTINE MatPrint (m1, m, n, matname)
      INCLUDE 'STDHDR.FOR'
      REAL m1(0:maxr,0:maxc)
      INTEGER m,n
      INTEGER i,j
      CHARACTER* (*) matname

      WRITE (*, 10)  matname
10    FORMAT (1X, A80)
      DO i = 0, m - 1
         DO j = 0, n - 1
           WRITE (*,20) m1(i, j)
20         FORMAT ( 5x,  F9.3 \)
         END DO
         WRITE (*,*)
      END DO
      END !SUB matPrint



      SUBROUTINE MatProd (m1, m2, l, m, n, m3)
      INCLUDE 'STDHDR.FOR'
      REAL m1(0:maxr,0:maxc), m2(0:maxr,0:maxc), m3(0:maxr,0:maxc)
      INTEGER l,m,n,i,j,k

      DO i = 0, l - 1
         DO j = 0, n - 1
            m3(i, j) = 0
            DO k = 0, m - 1
               m3(i, j) = m3(i, j) + m1(i, k) * m2(k, j)
            END DO
         END DO
      END DO
      END !SUB MatProd



      SUBROUTINE MatScalarProd (m1, sval, numrow, numcol, m2)
      INCLUDE 'STDHDR.FOR'
      REAL m1,m2,sval
      DIMENSION m1(0:maxr,0:maxc), m2(0:maxr,0:maxc)
      INTEGER numrow,numcol
      INTEGER i,j

      DO i = 0, numrow - 1
         DO j = 0, numcol - 1
            m2(i, j) = sval * m1(i, j)
         END DO
      END DO
      END !SUB  MatScalarProduct



      SUBROUTINE MatTranspose (m1, numrow, numcol, m2)
      INCLUDE 'STDHDR.FOR'
      REAL m1,m2
      DIMENSION m1(0:maxr,0:maxc), m2(0:maxr,0:maxc)
      INTEGER numrow,numcol
      INTEGER i,j

      DO i = 0, numrow - 1
         DO j = 0, numcol - 1
            m2(j, i) = m1(i, j)
         END DO
      END DO
      END !SUB MatTranspose



      SUBROUTINE NormalizeEigenVectors (a, n)
      INCLUDE 'STDHDR.FOR'
      REAL a, max
      DIMENSION a(0:maxr,0:maxc)
      INTEGER n, i, j

      DO j = 0, n - 1
        max = a(0, j)
        DO i = 1, n - 1
           IF (ABS(a(i, j)) .GT. ABS(max)) max = a(i, j)
        END DO
        DO i = 0, n - 1
          a(i, j) = a(i, j) / max
        END DO
      END DO
      END !SUB NormEigVec

