!     Last change:  DOS  27 Jul 2000    8:41 pm
!                                       *** copyright 2000 ***
! *** filename bidiac.f95               *** John F. Monahan **
!                                       **********************
      PROGRAM PBIDIAC
!  DEMONSTRATION PROGRAM BIDIAC -- BIDIAGONALIZES COMPLEX MATRIX
!   EXPNDC  EXPANDS UNITARY MATRIX FROM COMPACT STORAGE
!  FIRST STEP IN COMPUTING COMPLEX SINGULAR VALUE DECOMPOSITION
!
      IMPLICIT NONE
      INTEGER, PARAMETER :: MMAX = 8    ! BIGGEST PROBLEM
      INTEGER, PARAMETER :: NMAX = 4    ! IS 8 BY 4
      COMPLEX, DIMENSION(MMAX,NMAX) :: A,ACOPY
      COMPLEX, DIMENSION(MMAX,MMAX) :: V
      COMPLEX, DIMENSION(NMAX,NMAX) :: U
      REAL, DIMENSION(MMAX) :: PIV
      REAL, DIMENSION(NMAX) :: PIU
      COMPLEX, DIMENSION(NMAX) :: D,F
      COMPLEX S
      INTEGER M,N,I,J,K,L,ICASE
!                   INTERFACE BLOCK
      INTERFACE
      SUBROUTINE BIDIAC(A,M,N,D,F,PIU,PIV)
      complex, DIMENSION(:,:), INTENT(IN OUT) :: A
      complex, DIMENSION(:), INTENT(OUT) :: D,F
      REAL, DIMENSION(:), INTENT(OUT) :: PIU
      REAL, DIMENSION(:), INTENT(OUT) :: PIV
      INTEGER, INTENT(IN) :: M,N
      END SUBROUTINE BIDIAC
      SUBROUTINE EXPNDC(A,M,N,PIU,PIV,U,V)
      complex, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN) :: PIU
      REAL, DIMENSION(:), INTENT(IN) :: PIV
      complex, DIMENSION(:,:), INTENT(OUT) :: U
      complex, DIMENSION(:,:), INTENT(OUT) :: V
      INTEGER, INTENT(IN) :: M,N
      END SUBROUTINE EXPNDC
      END INTERFACE
!
  20  FORMAT(2I4)                                                       
  21  FORMAT(8F6.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(' UNITARY MATRIX V')
  32  FORMAT(' UNITARY MATRIX U')
!
      OPEN( UNIT=5, FILE='csvdex.dat')
      OPEN( UNIT=6, FILE='bidiac.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 BIDIAC(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 EXPNDC(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., 0. )
      DO K = 1,M
      DO L = 1,N
      S = S + conjg( 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 unitary
      WRITE(6,29) 
      DO I = 1,M
      DO J = 1,M
      S = ( 0., 0. )
      DO K = 1,M
      S = S + conjg( 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., 0. )
      DO K = 1,N
      S = S + conjg( 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 PBIDIAC
      SUBROUTINE BIDIAC(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) unitary MATRICES SO THAT
!         H
!       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 EXPNDC 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  MARCH 2000  FOR FORTRAN 95
!      complex form of BIDAIG
!
      IMPLICIT NONE
      complex, DIMENSION(:,:), INTENT(IN OUT) :: A
      complex, DIMENSION(:), INTENT(OUT) :: D,F
      REAL, DIMENSION(:), INTENT(OUT) :: PIU
      REAL, DIMENSION(:), INTENT(OUT) :: PIV
      INTEGER, INTENT(IN) :: M,N
      REAL ETA,AS,C
      complex 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( real(A(I,K))), abs(aimag(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
      AS = ( 0., 0. )
      DO I = K,M
      A(I,K) = A(I,K)/ETA
      AS = AS + A(I,K)*conjg( A(I,K) )
      END DO  ! LOOP ON I
      AS = SQRT(AS)
      c = Cabs( A(K,K) )
      s = AS
      if( c .GT. 0. ) s = A(K,K)*AS/C
      A(K,K) = A(K,K) + S
      D(K) = -S*ETA
      PIV(K) = AS * (AS + C )
!  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 + conjg( 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(real(A(K,I))), abs(aimag(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
      AS = ( 0., 0. )
      DO I = KP1,N
      A(K,I) = A(K,I) / ETA
      AS = AS +A(K,I)*conjg( A(K,I) )
      END DO  ! LOOP ON I
      AS = SQRT(AS)
      C = cabs( A(K,KP1) )
      S = AS
      if( c .gt. 0. ) S = A(K,KP1)*AS/C
      A(K,KP1) = A(K,KP1) + S
      PIU(K) = AS * (AS + C)
      F(K) = - S * ETA
!              NOW APPLY TRANSFORMATION TO REST OF MATRIX
      DO J = KP1,M
      S = ( 0., 0. )
      DO I = KP1,N
      S = S + conjg( 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 BIDIAC
      SUBROUTINE EXPNDC(A,M,N,PIU,PIV,U,V)
!  TAKES MATRIX A HOLDING HOUSEHOLDER TRANSFORMATIONS FROM BIDIAC
!  IN COMPACT FORM, WITH CONSTANT PI IN PIU AND PIV
!  AND CREATES MATRICES U AND V
!
      IMPLICIT NONE
      complex, DIMENSION(:,:), INTENT(IN) :: A
      REAL, DIMENSION(:), INTENT(IN) :: PIU
      REAL, DIMENSION(:), INTENT(IN) :: PIV
      complex, DIMENSION(:,:), INTENT(OUT) :: U
      complex, DIMENSION(:,:), INTENT(OUT) :: V
      INTEGER, INTENT(IN) :: M,N
      complex S
      INTEGER I,J,K,KP1,KVMAX,KUMAX
!                   CREATE IDENTITY IN V
      DO I = 1,M
      DO J = 1,M
      V(I,J) = ( 0., 0. )
      END DO  ! LOOP ON J
      V(I,I) = ( 1., 0. )
      END DO  ! LOOP ON I
!                   CREATE IDENTITY IN U
      DO I = 1,N
      DO J = 1,N
      U(I,J) = ( 0., 0. )
      END DO  ! LOOP ON J
      U(I,I) = ( 1., 0. )
      END DO  ! LOOP ON I
!                   HOW MANY TRANSFORMATIONS ?
      KVMAX = MIN0(M-1,N)                                               
      KUMAX = MIN0(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., 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 * conjg( 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., 0.)
      DO J = KP1,N
      S = S + U(I,J)*conjg( 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 EXPNDC
! *** end of filename bidiac.f95        *********************