!     Last change:  VGC  11 Mar 2003    9:11 pm
!************************************************************************
!*                                                          7/31/92     *
!* MARKOV CHAIN PROGRAM                                                 *
!*                                                                      *
!* by J R Nanney                                                        *
!*                                                    n                 *
!* Computes the n step transition probability matrix P  .               *
!* The 1 step transition probability matrix is put into file P.DAT.     *
!* It should have the number of rows on the first line. The number of   *
!* columns is the same as the number of rows and should not be put      *
!* into the file P.DAT. Then the matrix follows row after row. The file *
!* P.DAT for the FRESHMOUTH/BLECH/ZAPP problem would appear as follows  * 
!* in file P.DAT (without the stars, of course).                        *
!*                                                                      *
!*    3                                                                 *
!*    .8    .1     .1                                                   *
!*    .3    .2     .5                                                   *
!*    .2    .3     .5                                                   *
!*                                                                      *
!*                                                                      *
!* The output matrices will be found in file PN.DAT.                    *
!* Trial runs resulted in underflow errors which are due to numbers     *
!* approaching zero. Ignore the underflow warnings.                     *
!*                                                                      *
!************************************************************************


      PROGRAM MARKOV
      IMPLICIT NONE
      CHARACTER(LEN=1):: CC                                                   
      INTEGER:: STEP, ICC, I, J, K, N
      REAL,DIMENSION(1:18,1:18):: P, PN, C
      LOGICAL:: LL
      ICC = 0

      OPEN (UNIT = 1, FILE = 'P.DAT', STATUS = 'OLD')
      OPEN (UNIT = 2, FILE = 'PN.DAT', STATUS = 'NEW')

! Reads the 1 step transition probability matrix P from file P.DAT

      READ (1,*) N
      DO I = 1, N
         READ(1,*) (P(I,J), J= 1,N)
      END DO

! Initializes PN to an identity matrix

      DO I = 1,N
      DO J = 1,N
         IF (I == J) THEN
            PN(I,J) = 1.00
         ELSE
            PN(I,J) = 0.00
         END IF
      END DO
      END DO

! Multiplies P times PN and stores the product temporarily in C
! which is then read back into PN

      DO STEP = 1, 5000
         DO I = 1, N
         DO J = 1,N
            C(I,J) = 0.0
            DO K = 1,N
               C(I,J) = C(I,J) + P(I,K)*PN(k,J)
            END DO
         END DO
         END DO

         DO I = 1,N
         DO J = 1,N
            PN(I,J) = C(I,J)
         END DO
         END DO

         LL = (STEP <= 10)                                    &
     &        .OR. (STEP <= 100 .AND. MOD(STEP,10) == 0)      &
     &        .OR. (STEP <= 1000 .AND. MOD(STEP,100) == 0)    &
     &        .OR. (STEP <= 10000 .AND. MOD(STEP, 1000) == 0)
         
! For STEP values 1,2,3,...10,20,30,...100,200,300,...1000,2000...
! prints the n step transition probability matrix into file PN.DAT
! Carriage control (CC) keeps a matrix from being split between two pages
! of hardcopy. This might not work depending on paper size, 
! computer used, printer, etc. 
         IF (LL) THEN

            IF (ICC+N+3 .LE. 60) THEN
               CC = ' '
            ELSE 
               CC = char(12)
               ICC = 0
            END IF

            WRITE(2,'(A1, I4, " step transition probability matrix")') CC, STEP

            WRITE(2,*)
            DO I = 1,N
               WRITE (2,'(1X, 18F9.4)') (PN(I,J), J = 1,N)
            END DO
            WRITE(2,*)
            ICC = ICC + N + 3

         END IF
      END DO
      STOP
      END PROGRAM MARKOV