      PROGRAM CalcStats

      INCLUDE 'STDHDR.FOR'
      REAL mat1(0: maxr, 0: maxc),mat2(0: maxr, 0: maxc)
      REAL mat3(0: maxr, 0: maxc)
      REAL eigenvalues(0: maxc)
      REAL a(0:maxr,0:maxc), ev(0:maxr, 0:maxc)
      CHARACTER message * 80
      INTEGER n,count
      LOGICAL success

      ! matmath demo main

      mat1(0, 0) = 13.0
      mat1(0, 1) = -8.0
      mat1(0, 2) = -3.0
      mat1(1, 0) = -8.0
      mat1(1, 1) = 10.0
      mat1(1, 2) = -1.0
      mat1(2, 0) = -3.0
      mat1(2, 1) = -1.0
      mat1(2, 2) = 11.0

      mat2(0, 0) = 1.0
      mat2(0, 1) = 5.0
      mat2(0, 2) = 1.0
      mat2(1, 0) = 3.0
      mat2(1, 1) = 6.0
      mat2(1, 2) = 3.0
      mat2(2, 0) = 0.0
      mat2(2, 1) = 1.0
      mat2(2, 2) = -9.0

      CALL MatProd(mat1, mat2, 3, 3, 3, mat3)
      message = 'Matrix product mat1 x mat2 = '
      CALL MatPrint(mat3, 3, 3, message)
      WRITE (*,*) 'Press carriage return to advance to next procedure'
      READ (*,*)

      vm = 11
      CALL MatScalarProd(mat1, vm, 3, 3, mat3)
      message = 'Matrix scalar product mat1 x 11.000  = '
      CALL MatPrint(mat3, 3, 3, message)
      WRITE(*,*) 'Press carriage return to advance to next procedure'
      READ (*,*)


      CALL MatAdd(mat1, mat2, 3, 3, mat3)
      message = 'Matrix add mat1 + mat 2 = '
      CALL MatPrint(mat3, 3, 3, message)
      WRITE(*,*) 'Press carriage return to advance to next procedure'
      READ (*,*)


      CALL MatTranspose(mat1, 3, 3, mat3)
      message = 'Matrix transpose of mat1  = '
      CALL MatPrint(mat3, 3, 3, message)
      WRITE(*,*) 'Press carriage return to advance to next procedure'
      READ (*,*)


      CALL MatDeter(mat1, 3, det)
      WRITE (*,5) 'Determinant of matrix mat1 = ', det
5     FORMAT (A28, F7.2)
      WRITE(*,*)
      WRITE(*,*) 'Press carriage return to advance to next procedure'
      READ (*,*)

      CALL MatInvert(mat1, 3, det, mat3)
      message = 'Inverse of matrix mat1  = '
      CALL MatPrint(mat3, 3, 3, message)
      WRITE(*,*)'Press carriage return to begin Jacobi Demo '
      READ (*,*)

      n = 4

       a(0, 0) = 5.0
       a(0, 1) = 4.0
       a(0, 2) = 1.0
       a(0, 3) = 1.0
       a(1, 0) = 4.0
       a(1, 1) = 5.0
       a(1, 2) = 1.0
       a(1, 3) = 1.0
       a(2, 0) = 1.0
       a(2, 1) = 1.0
       a(2, 2) = 4.0
       a(2, 3) = 2.0
       a(3, 0) = 1.0
       a(3, 1) = 1.0
       a(3, 2) = 2.0
       a(3, 3) = 4.0


      count = 10
      CALL CyclicJacobi(a, n, eigenvalues, ev, count, success)
      CALL NormalizeEigenVectors(ev, n)
      WRITE (*,*) 'Success ', success
      DO i = 0, n - 1
        WRITE(*,7) 'Eigenvalue ', i, ' =  ', eigenvalues(i)
        WRITE(*,8) 'Eigenvector ', i, ' = '
        DO j = 0, n - 1
          WRITE(* ,10)  ev(j, i)
        END DO
        WRITE (*,*)
      END DO

      WRITE(*,*) 'count ', count
      READ (*,*)
7     FORMAT (A12, I2, A3, F7.2)
8     FORMAT (1X, A12, I2, A3 \)
10    FORMAT (F11.6\)



      END
