!     Last change:  DOS  27 Jul 2000    8:47 pm
!                                       *** copyright 2000 ***
! *** filename bidiag.f95               *** John F. Monahan **
!                                       **********************
      PROGRAM PBIDIAG
!  DEMONSTRATION PROGRAM BIDIAG -- BIDIAGONALIZES MATRIX
!   EXPNDB  EXPANDS ORTHOGONAL MATRIX FROM COMPACT STORAGE
!  FIRST STEP IN COMPUTING SINGULAR VALUE DECOMPOSITION
!
      IMPLICIT NONE
      INTEGER, PARAMETER :: MMAX = 8    ! BIGGEST PROBLEM
      INTEGER, PARAMETER :: NMAX = 4    ! IS 8 BY 4
      REAL, DIMENSION(MMAX,NMAX) :: A,ACOPY
      REAL, DIMENSION(MMAX,MMAX) :: V
      REAL, DIMENSION(NMAX,NMAX) :: U
      REAL, DIMENSION(MMAX) :: PIV
      REAL, DIMENSION(NMAX) :: D,F,PIU
      REAL S
      INTEGER M,N,I,J,K,L,ICASE
!                   INTERFACE BLOCK
      INTERFACE
      SUBROUTINE BIDIAG(A,M,N,D,F,PIU,PIV)
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, DIMENSION(:), INTENT(OUT) :: D,F,PIU
      REAL, DIMENSION(:), INTENT(OUT) :: PIV
      INTEGER, INTENT(IN) :: M,N
      END SUBROUTINE BIDIAG
      SUBROUTINE EXPNDB(A,M,N,PIU,PIV,U,V)
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN) :: PIU
      REAL, DIMENSION(:), INTENT(IN) :: PIV
      REAL, DIMENSION(:,:), INTENT(OUT) :: U
      REAL, DIMENSION(:,:), INTENT(OUT) :: V
      INTEGER, INTENT(IN) :: M,N
      END SUBROUTINE EXPNDB
      END INTERFACE
!
  20  FORMAT(2I4)                                                       
  21  FORMAT(5F6.2)                                                     
  22  FORMAT(8F10.6)                                                    
!
  25  FORMAT(/' ORIGINAL MATRIX WITH DIMENSIONS M=',I2,'  N=',I2)
  26  FORMAT(' BIDIAGONALIZED, FIRST DIAGONAL, THEN SUPER ')
  27  FORMAT(' COMPACT STORAGE OF TRANSFORMATIONS')
  28  FORMAT(' IS V''A U BIDIAGONAL ? ')
  29  FORMAT(' IS V ORTHOGONAL, V''V = I ?')
  30  FORMAT(' IS U ORTHOGONAL, U''U = I ?')
  31  FORMAT(' ORTHOGONAL MATRIX V')
  32  FORMAT(' ORTHOGONAL MATRIX U')
!
      OPEN( UNIT=5, FILE='svdex.dat')
      OPEN( UNIT=6, FILE='bidiag.out')
!
      DO ICASE = 1,3          ! THREE PROBLEMS
!
      READ(5,20) M,N
!                      WRITE OUT THE MATRIX
      WRITE(6,25) M,N 
      DO I = 1,M
      READ(5,21) (A(I,J),J=1,N)
      WRITE(6,21) (A(I,J),J=1,N)                                        
!                      MAKE A COPY
      DO J = 1,N
      ACOPY(I,J) = A(I,J)                                               
      END DO  ! LOOP ON J
      END DO  ! LOOP ON I
!                      BIDIAGONALIZE
      CALL BIDIAG(A,M,N,D,F,PIU,PIV)                                    
!
      WRITE(6,26)                                                    
      WRITE(6,22) (D(I),I=1,N)                                          
      WRITE(6,22) (F(I),I=1,N)                                          
      WRITE(6,27)                                                    
      DO I = 1,M
      WRITE(6,22) (A(I,J),J=1,N)                                        
      END DO  ! LOOP ON I
!                      EXPAND FROM COMPACT STORAGE
      CALL EXPNDB(A,M,N,PIU,PIV,U,V)                                    
!                   MULTIPLY V'A U TO CHECK IF IT'S DIAGONAL
      WRITE(6,28)  
      DO I = 1,M
      DO J = 1,N
      S = 0.                                                            
      DO K = 1,M
      DO L = 1,N
      S = S + V(K,I)*ACOPY(K,L)*U(L,J)                                  
      END DO  ! LOOP ON L
      END DO  ! LOOP ON K
      A(I,J) = S
      END DO  ! LOOP ON J
      WRITE(6,22) (A(I,J),J=1,N)
      END DO  ! LOOP ON I
!                        NOW CHECK TO SEE IF U AND V ARE ORTHOGONAL
      WRITE(6,29) 
      DO I = 1,M
      DO J = 1,M
      S = 0.                                                            
      DO K = 1,M
      S = S + V(I,K)*V(J,K)
      END DO  ! LOOP ON K
      D(J) = S                                                          
      END DO  ! LOOP ON J
      WRITE(6,22) (D(J),J=1,M)
      END DO  ! LOOP ON I
!
      WRITE(6,30) 
      DO I = 1,N
      DO J = 1,N
      S = 0.                                                            
      DO K = 1,N
      S = S + U(I,K)*U(J,K)
      END DO  ! LOOP ON K
      D(J) = S
      END DO  ! LOOP ON J
      WRITE(6,22) (D(J),J=1,N)
      END DO  ! LOOP ON I
!
      WRITE(6,31) 
      DO I = 1,M
      WRITE(6,22) (V(I,J),J=1,M)                                        
      END DO  ! LOOP ON I
      WRITE(6,32)
      DO I = 1,N
      WRITE(6,22) (U(I,J),J=1,N)                                        
      END DO  ! LOOP ON I
!                             END OF THIS PROBLEM
      END DO  ! LOOP ON ICASE
      STOP
      END PROGRAM PBIDIAG
      SUBROUTINE BIDIAG(A,M,N,D,F,PIU,PIV)
!  BIDIAGONALIZES MATRIX A (M BY N) USING HOUSEHOLDER TRANSFORMATIONS
!  U (N BY N) AND V (M BY M) ORTHOGONAL MATRICES SO THAT
!         T
!       V   A  U  =  B    ***  A IS M BY N WITH M .GE. N  ***
!
!  WHERE B IS BIDIAGONAL WITH DIAGONAL ELEMENTS STORED IN D AND
!  SUPERDIAGONAL ELEMENTS STORED IN F
!
!  MATRICES U AND V STORED IN COMPACT FORM OVERWRITING A
!  PI VALUES OF HOUSEHOLDER TRANSFORMATIONS STORED IN PIU, PIV
!  CALL EXPNDB TO EXPAND FROM COMPACT FORM TO FULL MATRICES U AND V
!
!  ADAPTED FROM GOLUB - REINSCH ALGORITHM
!  J F MONAHAN (JULY,1986) DEPT OF STATISTICS, N C S U, RALEIGH NC USA
!                       CORRECTED JULY, 1987
!      RECODED  SEPTEMBER 1999, APRIL 2000  FOR FORTRAN 95
!
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN OUT) :: A
      REAL, DIMENSION(:), INTENT(OUT) :: D,F,PIU
      REAL, DIMENSION(:), INTENT(OUT) :: PIV
      INTEGER, INTENT(IN) :: M,N
      REAL ETA,S
      INTEGER I,J,K,KP1
!
!                   DO N STEPS ON LEFT, N-2 STEPS ON RIGHT
!                    *** BECAUSE M .GT. N ***
      DO K = 1,N
      KP1 = K + 1
!                  DO V TRANSFORMATION, PUT ZEROS IN COLUMN KC
      ETA = 0
      DO I = K,M
      ETA = MAX(ETA, ABS(A(I,K)) )
      END DO  ! LOOP ON I
      IF( (ETA .EQ. 0.) .OR. (K .EQ. M) ) THEN
!                 IF ALL ZERO THEN DO NO TRANSFORMATION
      PIV(K) = 0.
      D(K) = A(K,K)
                                          ELSE
      S = 0.
      DO I = K,M
      A(I,K) = A(I,K)/ETA
      S = S + A(I,K)*A(I,K)
      END DO  ! LOOP ON I
      S = SIGN( SQRT(S), A(K,K) )
      A(K,K) = A(K,K) + S
      D(K) = -S*ETA
      PIV(K) = S*A(K,K)
!  PI VALUE FOR TFM STORED, VECTOR IN COLUMN K (K,M)
!
      IF( K .LT. N ) THEN
!                APPLY TRANSFORM TO REMAINDER OF MATRIX
      DO J = KP1,N
      S = 0.                                                            
      DO I = K,M
      S = S + A(I,K)*A(I,J)
      END DO  ! LOOP ON I
      S = S / PIV(K)
      DO I = K,M
      A(I,J) = A(I,J) - S*A(I,K)
      END DO  ! LOOP ON I
      END DO  ! LOOP ON J
!
!                DONE WITH THIS TRANSFORMATION ON COLUMN K
!
                      END IF ! ( KC .LE. N )
                                     END IF ! ( (ETA .EQ. 0.)
                                            ! OR (K .EQ. M) )
      IF( K .LE. N-2 ) THEN
!               NOW TRANSFORM ON RIGHT AS U -- ZEROS IN ROWS
      ETA = 0.
      DO I = KP1,N
      ETA = MAX( ETA, ABS(A(K,I)) )
      END DO  ! LOOP ON I
      IF( ETA .EQ. 0. ) THEN
!             IF RELEVANT PART OF ROW IS ALL ZERO, THEN NO TFM
      PIU(K) = 0.
      F(K) = 0.
                        ELSE
      S = 0.
      DO I = KP1,N
      A(K,I) = A(K,I) / ETA
      S = S +A(K,I)*A(K,I)
      END DO  ! LOOP ON I
      S = SIGN( SQRT(S), A(K,KP1) )
      A(K,KP1) = A(K,KP1) + S
      PIU(K) = S*A(K,KP1)
      F(K) = - S * ETA
!              NOW APPLY TRANSFORMATION TO REST OF MATRIX
      DO J = KP1,M
      S = 0.                                                            
      DO I = KP1,N
      S = S + A(K,I)*A(J,I)
      END DO  ! LOOP ON I
      S = S / PIU(K)
      DO I = KP1,N
      A(J,I) = A(J,I) - S * A(K,I)
      END DO  ! LOOP ON I
      END DO  ! LOOP ON J
!
                        END IF ! ( ETA .EQ. 0. )
                       END IF ! ( K .LE. N-2 )
      END DO  ! LOOP ON K
      F(N-1) = A(N-1,N)
      RETURN                                                            
      END SUBROUTINE BIDIAG
      SUBROUTINE EXPNDB(A,M,N,PIU,PIV,U,V)
!  TAKES MATRIX A HOLDING HOUSEHOLDER TRANSFORMATIONS FROM BIDIAG
!  IN COMPACT FORM, WITH CONSTANT PI IN PIU AND PIV
!  AND CREATES MATRICES U AND V
!
      IMPLICIT NONE
      REAL, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN) :: PIU
      REAL, DIMENSION(:), INTENT(IN) :: PIV
      REAL, DIMENSION(:,:), INTENT(OUT) :: U
      REAL, DIMENSION(:,:), INTENT(OUT) :: V
      INTEGER, INTENT(IN) :: M,N
      REAL S
      INTEGER I,J,K,KP1,KVMAX,KUMAX
!                   CREATE IDENTITY IN V
      DO I = 1,M
      DO J = 1,M
      V(I,J) = 0.                                                       
      END DO  ! LOOP ON J
      V(I,I) = 1.
      END DO  ! LOOP ON I
!                   CREATE IDENTITY IN U
      DO I = 1,N
      DO J = 1,N
      U(I,J) = 0.
      END DO  ! LOOP ON J
      U(I,I) = 1.
      END DO  ! LOOP ON I
!                   HOW MANY TRANSFORMATIONS ?
      KVMAX = MIN(M-1,N)
      KUMAX = MIN(M,N-2)
!
      DO K = 1,KVMAX
!                IF NO TRANSFORMATION SKIP IT ALL
      IF( PIV(K) .NE. 0. ) THEN
      DO I = 1,M
      S = 0.                                                            
      DO J = K,M
      S = S + V(I,J)*A(J,K)                                             
      END DO  ! LOOP ON J
      S = S / PIV(K)
      DO J = K,M
      V(I,J) = V(I,J) - S * A(J,K)                                      
      END DO  ! LOOP ON J
      END DO  ! LOOP ON I
                           END IF ! ( PIV(K) .NE. 0. )
      END DO  ! LOOP ON K -- TRANSFORMATIONS
!                   NOW DO THE U TRANSFORMATIONS
      IF( KUMAX .LE. 0 ) RETURN
      DO K = 1,KUMAX
      IF(PIU(K) .NE. 0. ) THEN
      KP1 = K + 1                                                       
      DO I = 1,N
      S = 0.                                                            
      DO J = KP1,N
      S = S + U(I,J)*A(K,J)                                             
      END DO  ! LOOP ON J
      S = S / PIU(K)
      DO J = KP1,N
      U(I,J) = U(I,J) - S * A(K,J)                                      
      END DO  ! LOOP ON J
      END DO  ! LOOP ON I
                          END IF ! (PIU(K) .NE. 0.)
      END DO  ! LOOP ON K  -- TRANSFORMATIONS
      RETURN
      END SUBROUTINE EXPNDB
! *** end of filename bidiag.f95        *********************