!     Last change:  DOS  27 Jul 2000    8:54 pm
!                                       *** copyright 2000 ***
! *** filename orthit.f95               *** John F. Monahan **
!                                       **********************
      PROGRAM ORTHIT
!  DEMONSTRATION PROGRAM ORTHIT
!  TEST OF ORTHOGONAL ITERATION AND MODIFIED GRAM-SCHMIDT
!
      IMPLICIT NONE
      REAL, DIMENSION(10,10) :: A
      REAL, DIMENSION(10,5) :: Z,Q
      REAL, DIMENSION(5,5) :: R
      REAL S
      INTEGER ICASE,I,J,K,P,N,NIT,KIT
!                   INTERFACE
      INTERFACE
      SUBROUTINE MGSQR(B,N,P,Q,R)
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: B,Q
      REAL, DIMENSION(:,:), INTENT(OUT) :: R
      INTEGER, INTENT(IN) :: N,P
      END SUBROUTINE MGSQR
      END INTERFACE
!                   INPUT FORMATS
  25  FORMAT(3I4)                                                       
  26  FORMAT(8F9.4)                                                     
!                   OUTPUT FORMATS
  20  FORMAT(4X,'ITERATION NUMBER',I4)                                  
  21  FORMAT(8F10.6)
  27  FORMAT(/' PROBLEM NUMBER',I4/'DIMENSION, NUMBER, ITERATIONS',3I4)
  28  FORMAT(' Q MATRIX FROM ITERATION',I4)
  29  FORMAT(' R MATRIX FROM ITERATION',I4)
!                   FIVE EXAMPLES IN 'orthit.dat'
      OPEN( UNIT=5, FILE='orthit.dat')
      OPEN( UNIT=6, FILE='orthit.out')
!                   FIRST ONE IS STATIONARY VECTOR -- SLOW CONVERGENCE
!                     STATIONARY PROBS (36, 12, 24, 12, 18, 10, 9)/121
!                   REMEMBER NORM IS EUCLID, NOT L1
!
      DO ICASE = 1,5
!                   READ SIZE, HOW MANY VECTORS, NIT, THEN MATRIX
      READ(5,25) N,P,NIT
      WRITE(6,27) ICASE,N,P,NIT
!                   WRITE OUT MATRIX
      DO I = 1,N
      READ(5,26) (A(I,J),J=1,N)
      WRITE(6,26) (A(I,J),J=1,N)                                        
      DO J = 1,P    ! COPY FIRST P COLUMNS TO START
      Q(I,J) = A(I,J)                                                   
      END DO  ! LOOP ON J
      END DO  ! LOOP ON I
!                   START ITERATION
      DO KIT = 1,NIT
      WRITE(6,20) KIT                                                   
!                   MULTIPLY
      DO J = 1,P
      DO I = 1,N
      S = 0.                                                            
      DO K = 1,N
      S = S + A(I,K)*Q(K,J)                                             
      END DO  ! LOOP ON K
      Z(I,J) = S
      END DO  ! LOOP ON I
      END DO  ! LOOP ON J
!                   DO QR FACTORIZATION
      CALL MGSQR(Z,N,P,Q,R)
!                   WRITE OUT RESULTS  - SOMETIMES
      IF( MOD(KIT,5) .EQ. 0 ) THEN
      WRITE(6,28) KIT 
      DO I = 1,N
      WRITE(6,21) (Q(I,J),J=1,P)
      END DO  ! LOOP ON I
      WRITE(6,29) KIT
      DO I = 1,P
      WRITE(6,21) (R(I,J),J=1,P)
      END DO  ! LOOP ON I
                              END IF
      END DO  ! LOOP ON KIT
!
      END DO  ! LOOP ON ICASE
      STOP
      END PROGRAM ORTHIT
      SUBROUTINE MGSQR(B,N,P,Q,R)
!  COMPUTES MODIFIED GRAM-SCHMIDT FACTORIZATION OF B = Q R
!  WHERE Q HAS P ORTHONORMAL COLUMNS AND R IS UPPER TRIANGULAR
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: B,Q
      REAL, DIMENSION(:,:), INTENT(OUT) :: R
      INTEGER, INTENT(IN) :: N,P
      REAL S
      INTEGER J,K,JP1,I
!                   MODIFIED GRAM-SCHMIDT -- ONLY FIRST IP COLS
      DO J = 1,P
!                   COMPUTE NORM OF COL J
      S = 0.                                                            
      DO K = 1,N
      S = S + B(K,J)*B(K,J)                                             
      END DO  ! LOOP ON K
      S = SQRT(S)
      R(J,J) = S                                                        
!                   NORMALIZE
      DO K = 1,N
      Q(K,J) = B(K,J) / S                                               
      END DO  ! LOOP ON K
!                   ARE WE DONE ?
      IF ( J .NE. P ) THEN
!                   SUBTRACT FROM OTHER COLUMNS
      JP1 = J + 1                                                       
      DO I = JP1,P
      S = 0.                                                            
      DO K = 1,N
      S = S + Q(K,J)*B(K,I)                                             
      END DO  ! LOOP ON K
      R(J,I) = S
      DO K = 1,N
      B(K,I) = B(K,I) - R(J,I)*Q(K,J)
      END DO  ! LOOP ON K
      END DO  ! LOOP ON I
                      END IF
!                   DONE WITH COL J
!
!               DONE WITH ALL IP COLUMNS
      END DO  ! LOOP ON J
      RETURN
      END SUBROUTINE MGSQR
! *** end of filename orthit.f95        *********************