!+ Subroutine to perform SVD decomposition of a matrix.
!
! Method:
!	This routine is a direct implementation of that given in
!	"Numerical recipes" by D. E. Knuth.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE SVD_DECOMPOSE(IERR
     &   , A, N_ROW, N_COLUMN
     &   , ROW_DIMENSION, COLUMN_DIMENSION
     &   , W, V, RV1
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision_svd.tcf'
      INCLUDE 'incl/svd_parameters.acf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS
      INTEGER 	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER 	!, INTENT(IN)
     &     N_ROW
!             NUMBER OF ROWS
     &   , N_COLUMN
!             NUMBER OF COLUMNS
     &   , ROW_DIMENSION
!             FIRST DIMENSION
     &   , COLUMN_DIMENSION
!             SECOND DIMENSION
      REAL*8 	!, INTENT(INOUT)
     &     A(ROW_DIMENSION, COLUMN_DIMENSION)
!             SVD MATRIX: OVERWRITTEN
      REAL*8 	!, INTENT(OUT)
     &     W(ROW_DIMENSION)
!             DIAGONAL TERMS
     &   , V(COLUMN_DIMENSION, COLUMN_DIMENSION)
!             SECOND MATRIX
     &   , RV1(COLUMN_DIMENSION)
!             WORKING ARRAY
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     I
!             LOOP VARIABLE
     &   , J
!             "
     &   , K
!             "
     &   , L
!             "
     &   , ITERATION
!             ITERATION NUMBER
     &   , NM
!             TEMPORARY SCALAR
     &   , JJ
!             "
      REAL*8
     &     SCALE
!             THESE VARIABLES ARE JUST
     &   , ANORM
!             PART OF THE ALGORITHM
     &   , G
!             AND CNNOT MEANINGFULLY
     &   , F
!             BE DESCRIBEED
     &   , H
     &   , Y
     &   , Z
     &   , S
     &   , C
     &   , X
!
!
!     CHECK THAT THERE ARE NOT TOO MANY POINTS IN THE FIT.
      IF (N_COLUMN.GT.N_ROW) THEN
         WRITE(IU_ERR, '(A)')
     &      'SVD ATTEMPTED WITH TOO MANY UNKNOWNS.'
         IERR=I_ABORT_CALCULATION
         RETURN
      ENDIF
!
!     HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM.
      G=0.0D+00
      SCALE=0.0D+00
      ANORM=0.0D+00
      DO I=1, N_COLUMN
         L=I+1
         RV1(I)=SCALE*G
         G=0.0D+00
         S=0.0D+00
         SCALE=0.0D+00
         IF (I.LE.N_ROW) THEN
            DO K=I, N_ROW
               SCALE=SCALE+ABS(A(K,I))
            ENDDO
            IF (SCALE.GT.TOL_SVD) THEN
               DO K=I, N_ROW
                  A(K, I)=A(K, I)/SCALE
                  S=S+A(K, I)*A(K, I)
               ENDDO
               F=A(I, I)
               G=-SIGN(SQRT(S), F)
               H=F*G-S
               A(I, I)=F-G
               IF (I.NE.N_COLUMN) THEN
                  DO J=L, N_COLUMN
                     S=0.0D+00
                     DO K=I, N_ROW
                        S=S+A(K, I)*A(K, J)
                     ENDDO
                     F=S/H
                     DO K=I, N_ROW
                        A(K, J)=A(K, J)+F*A(K, I)
                     ENDDO
                  ENDDO
               ENDIF
               DO K=I, N_ROW
                  A(K, I)=SCALE*A(K, I)
               ENDDO
            ENDIF
         ENDIF
!
         W(I)=SCALE*G
         G=0.0D+00
         S=0.0D+00
         SCALE=0.0D+00
         IF ( (I.LE.N_ROW).AND.(I.NE.N_COLUMN) ) THEN
            DO K=L, N_COLUMN
               SCALE=SCALE+ABS(A(I, K))
            ENDDO
            IF (SCALE.GT.TOL_SVD) THEN
               DO K=L, N_COLUMN
                  A(I, K)=A(I, K)/SCALE
                  S=S+A(I, K)*A(I, K)
               ENDDO
               F=A(I, L)
               G=-SIGN(SQRT(S), F)
               H=F*G-S
               A(I, L)=F-G
               DO K=L, N_COLUMN
                  RV1(K)=A(I, K)/H
               ENDDO
               IF (I.NE.N_ROW) THEN
                  DO J=L, N_ROW
                     S=0.0D+00
                     DO K=L, N_COLUMN
                        S=S+A(J, K)*A(I, K)
                     ENDDO
                     DO K=L, N_COLUMN
                        A(J, K)=A(J, K)+S*RV1(K)
                     ENDDO
                  ENDDO
               ENDIF
               DO K=L, N_COLUMN
                  A(I, K)=SCALE*A(I, K)
               ENDDO
            ENDIF
         ENDIF
         ANORM=MAX(ANORM, (ABS(W(I))+ABS(RV1(I))))
      ENDDO
!
!     ACCUMULATE THE RIGHT-HAND SIDES OF THE TRANSFORMATIONS.
      DO I=N_COLUMN, 1, -1
         IF (I.LT.N_COLUMN) THEN
            IF (ABS(G).GT.TOL_SVD) THEN
               DO J=L, N_COLUMN
                  V(J, I)=(A(I, J)/A(I, L))/G
               ENDDO
               DO J=L, N_COLUMN
                  S=0.0D+00
                  DO K=L, N_COLUMN
                     S=S+A(I, K)*V(K, J)
                  ENDDO
                  DO K=L, N_COLUMN
                     V(K, J)=V(K, J)+S*V(K, I)
                  ENDDO
               ENDDO
            ENDIF
            DO J=L, N_COLUMN
               V(I, J)=0.0D+00
               V(J, I)=0.0D+00
            ENDDO
         ENDIF
         V(I, I)=1.0D+00
         G=RV1(I)
         L=I
      ENDDO
!
!     ACCUMULATE THE LEFT-HAND SIDES OF THE TRANSFORMATIONS.
      DO I=N_COLUMN, 1, -1
         L=I+1
         G=W(I)
         IF (I.LT.N_COLUMN) THEN
            DO J=L, N_COLUMN
               A(I, J)=0.0D+00
            ENDDO
         ENDIF
	 IF (ABS(G).GT.TOL_SVD) THEN
            G=1.0D+00/G
            IF (I.NE.N_COLUMN) THEN
               DO J=L, N_COLUMN
                  S=0.0D+00
                  DO K=L, N_ROW
                     S=S+A(K, I)*A(K, J)
                  ENDDO
                  F=(S/A(I, I))*G
                  DO K=I, N_ROW
                     A(K, J)=A(K, J)+F*A(K, I)
                  ENDDO
               ENDDO
            ENDIF
            DO J=I, N_ROW
               A(J, I)=A(J, I)*G
            ENDDO
         ELSE
            DO J=I, N_ROW
               A(J, I)=0.0D+00
            ENDDO
         ENDIF
         A(I, I)=A(I, I)+1.0D+00
      ENDDO
!
!     DIAGONALIZE THE BIDIAGONAL FORM.
      DO K=N_COLUMN, 1, -1
         DO ITERATION=1, MAX_SVD_ITERATION
            DO L=K, 1, -1
               NM=L-1
               IF (ABS(RV1(L)).LT.TOL_SVD) GOTO 2
               IF (ABS(W(NM)).LT.TOL_SVD) GOTO 1
            ENDDO
!
1           C=0.0D+00
            S=1.0D+00
            DO I=L, K
               F=S*RV1(I)
               IF (ABS(F).LT.TOL_SVD) THEN
                  G=W(I)
                  H=SQRT(F*F+G*G)
                  W(I)=H
                  H=1.0D+00/H
                  C=G*H
                  S=-F*H
                  DO J=1, N_ROW
                     Y=A(J, NM)
                     Z=A(J, I)
                     A(J, NM)=(Y*C)+(Z*S)
                     A(J, I)=-(Y*S)+(Z*C)
                  ENDDO
               ENDIF
            ENDDO
!
2           Z=W(K)
            IF (L.EQ.K) THEN
               IF (Z.LT.0.0D+00) THEN
                  W(K)=-Z
                  DO J=1, N_COLUMN
                     V(J, K)=-V(J, K)
                  ENDDO
               ENDIF
               GOTO 3
            ENDIF
!
            IF (ITERATION.EQ.MAX_SVD_ITERATION) THEN
               WRITE(IU_MONITOR, '(/A)')
     &            'CONVERGENCE OF SVD NOT ACHIEVED:TOO MANY ITERATIONS.'
            ENDIF
            X=W(L)
            NM=K-1
            Y=W(NM)
            G=RV1(NM)
            H=RV1(K)
            F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0D+00*H*Y)
            G=SQRT(F*F+1.0D+00)
            F=((X-Z)*(X+Z)+H*((Y/(F+SIGN(G, F)))-H))/X
!
!           MAKE THE QR-TRANSFORMATION.
            C=1.0D+00
            S=1.0D+00
            DO J=L, NM
               I=J+1
               G=RV1(I)
               Y=W(I)
               H=S*G
               G=C*G
               Z=SQRT(F*F+H*H)
               RV1(J)=Z
               C=F/Z
               S=H/Z
               F=(X*C)+(G*S)
               G=-(X*S)+(G*C)
               H=Y*S
               Y=Y*C
               DO JJ=1, N_COLUMN
                  X=V(JJ, J)
                  Z=V(JJ, I)
                  V(JJ, J)=(X*C)+(Z*S)
                  V(JJ, I)=-(X*S)+(Z*C)
               ENDDO
               Z=SQRT(F*F+H*H)
               W(J)=Z
               IF (Z.GT.TOL_SVD) THEN
                  Z=1.0D+00/Z
                  C=F*Z
                  S=H*Z
               ENDIF
               F=(C*G)+(S*Y)
               X=-(S*G)+(C*Y)
               DO JJ=1, N_ROW
                  Y=A(JJ, J)
                  Z=A(JJ, I)
                  A(JJ, J)=(Y*C)+(Z*S)
                  A(JJ, I)=-(Y*S)+(Z*C)
               ENDDO
            ENDDO
            RV1(L)=0.0D+00
            RV1(K)=F
            W(K)=X
         ENDDO
3     CONTINUE
      ENDDO
!
!
      RETURN
      END
