C     ALGORITHM 581, COLLECTED ALGORITHMS FROM ACM. THIS WORK
C     PUBLISHED IN TRANS. MATH. SOFTWARE, 8(1), PP. 84-88.
 
C***********************************************************************
C
C       THIS FILE CONTAINS PROGRAMS AND SUBROUTINES THAT ARE RELATED
C       TO THE PAPER 'AN IMPROVED ALGORITHM FOR COMPUTING THE SINGULAR
C       VALUE DECOMPOSITION' BY T.F. CHAN, WHICH WAS SUBMITTED FOR
C       PUBLICATIONS IN ACM TOMS.
C
C       THIS FILE CONTAINS ONE (1) DRIVER, SIX  (6) MAIN SUBROUTINES,
C       SOME BLAS ROUTINES, A ROUTINE FOR COMPUTING THE MACHINE RELATIVE
C       PRECISION AND THE OUTPUT PRODUCED BY THE DRIVER.
C       THEY ARE LISTED BELOW IN THE ORDER THAT THEY APPEAR IN THE
C       FILE:
C
C           1) DRIVER FOR THE ROUTINE HYBSVD TO COMPUTE THE SVD
C              OF A TEST MATRIX AND ITS TRANSPOSE.
C
C           2) ROUTINE HYBSVD ... HYBRID ALGORITHM ROUTINE.
C
C           3) ROUTINE MGNSVD ... HYBRID ALGORITHM, ASSUMES M .GE. N.
C
C           4) ROUTINE GRSVD  ... GOLUB-REINSCH ALGORITHM ROUTINE
C
C           5) BLAS ROUTINE   ... SSWAP, FOR SWAPPING VECTORS
C
C           6) ROUTINE SRELPR ... COMPUTES MACHINE PRECISION
C
C           7) ROUTINE SSVDC  ... LINPACK SVD ROUTINE
C
C           8) BLAS ROUTINES  ... SAXPY, SDOT, SNRM2, SROT, SROTG,
C                                 SSCAL.
C
C           9) ROUTINE SVD    ... EISPACK SVD ROUTINE
C
C          10) ROUTINE MINFIT ... EISPACK MINFIT ROUTINE
C
C          11) OUTPUT PRODUCED BY THE DRIVER IN 1).
C
C       THE ROUTINES 2) TO 6) CONSTITUTE THE PACKAGE THAT IMPLEMENTS
C       THE NEW HYBRID ALGORITHM AND CAN BE USED BY THEMSELVES.
C
C       PLEASE ADDRESS COMMENTS AND SUGGESTIONS TO:
C
C           TONY CHAN
C           COMPUTER SCIENCE DEPT., YALE UNIV.,
C           BOX 2158, YALE STATION,
C           NEW HAVEN, CT 06520.
C
C***********************************************************************
C
C       DRIVER FOR HYBSVD, A ROUTINE FOR COMPUTING THE SVD OF A MATRIX.
C       RESULTS ARE COMPARED TO THOSE PRODUCED BY SSVDC FROM LINPACK,
C       SVD AND MINFIT FROM EISPACK.
C       PREPARED FOR ACM TRANS. ON MATH. SOFTWARE.
C
C       MATRIX TESTED:  A(I,J) = 1/(I+J),  1ST RHS(I) = 1
C               T                          2ND RHS(I) = I
C       (A AND A )
C
C       PROGRAMMED BY
C           TONY CHAN
C           COMPUTER SCIENCE DEPT., YALE UNIV.,
C           BOX 2158, YALE STATION,
C           NEW HAVEN, CT 06520.
C
C       LAST REVISION.. JANUARY, 1982.
C
      REAL AE(15,15), AH(10,10), AL(20,20), WORK(10), E(10)
      REAL UH(11,10), VH(12,20), Z(13,20), BH(14,5), WH(10)
      REAL UE(15,10), VE(15,10), BE(15,5), WE(10)
      REAL UL(18,10), VL(19,10), WL(10)
      LOGICAL MATU, MATV
      INTEGER NUH, NVH, NZ, NBH
      INTEGER NUE, NVE, NBE
      INTEGER NUL, NVL
      INTEGER NAE, NAH, NAL, M, N, IRHS, IERR, NOUT, JOB
      DATA NAH /10/, NUH /11/, NVH /12/, NZ /13/, NBH /14/
      DATA NUE /15/, NVE /15/, NBE /15/, NAE /15/
      DATA NUL /18/, NVL /19/, NAL /20/
C
C     THIS IS OUTPUT UNIT NUMBER.
      DATA NOUT /6/
C
C     FIRST CASE... A AND RHS B.
C
      M = 8
      N = 4
      IRHS = 2
C
      DO 50 IU=1,2
        DO 40 IV=1,2
          IF (IU.EQ.2) MATU = .TRUE.
          IF (IU.EQ.1) MATU = .FALSE.
          IF (IV.EQ.2) MATV = .TRUE.
          IF (IV.EQ.1) MATV = .FALSE.
C
C       JOB PARAMETER NEEDED BY SSVDC.
          JOB = (IU-1)*20 + IV - 1
C
C       SET UP MATRIX
C
          DO 20 I=1,M
            BH(I,1) = 1.0
            BH(I,2) = FLOAT(I)
            BE(I,1) = 1.0
            BE(I,2) = FLOAT(I)
            DO 10 J=1,N
              AH(I,J) = 1./FLOAT(I+J)
              AE(I,J) = AH(I,J)
              AL(I,J) = AH(I,J)
   10       CONTINUE
   20     CONTINUE
C
          WRITE (NOUT,99999)
99999     FORMAT (/4H A ./)
          DO 30 I=1,M
            WRITE (NOUT,99998) (AH(I,J),J=1,N)
99998       FORMAT (5E15.7)
   30     CONTINUE
C
          WRITE (NOUT,99997)
99997     FORMAT (//45H ******************* HYBSVD *****************//)
          CALL HYBSVD(NAH, NUH, NVH, NZ, NBH, M, N, AH, WH, MATU, UH,
     *     MATV, VH, Z, BH, IRHS, IERR, WORK)
          CALL PRINTS(WH, UH, VH, IERR, BH, NUH, NVH, NBH, MATU, MATV,
     *     M, N, NOUT, IRHS)
C
          WRITE (NOUT,99996)
99996     FORMAT (//45H ******************* EISPAK *****************//)
          CALL SVD(NUE, M, N, AE, WE, MATU, UE, MATV, VE, IERR, WORK)
          CALL MINFIT(NUE, M, N, AE, WE, IRHS, BE, IERR, WORK)
          CALL PRINTS(WE, UE, VE, IERR, BE, NUE, NVE, NBE, MATU, MATV,
     *     M, N, NOUT, IRHS)
C
          WRITE (NOUT,99995)
99995     FORMAT (//45H ******************* SSVDC  *****************//)
          CALL SSVDC(AL, NAL, M, N, WL, E, UL, NUL, VL, NVL, WORK, JOB,
     *     IERR)
          CALL PRINTS(WL, UL, VL, IERR, E, NUL, NVL, NBH, MATU, MATV,
     *     M, N, NOUT, IRHS)
C
          CALL DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
     *     NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU,
     *     MATV)
   40   CONTINUE
   50 CONTINUE
C                      T
C     SECOND CASE...  A
C
      M = 4
      N = 8
      IRHS = 2
C
      DO 100 IU=1,2
        DO 90 IV=1,2
          IF (IU.EQ.2) MATU = .TRUE.
          IF (IU.EQ.1) MATU = .FALSE.
          IF (IV.EQ.2) MATV = .TRUE.
          IF (IV.EQ.1) MATV = .FALSE.
C
C       JOB PARAMETER NEEDED BY SSVDC.
          JOB = (IU-1)*20 + IV - 1
C
C       SET UP MATRIX
C
          DO 70 I=1,M
            BH(I,1) = 1.
            BH(I,2) = FLOAT(I)
            BE(I,1) = 1.
            BE(I,2) = FLOAT(I)
            DO 60 J=1,N
              AH(I,J) = 1./FLOAT(I+J)
              AE(I,J) = AH(I,J)
              AL(I,J) = AH(I,J)
   60       CONTINUE
   70     CONTINUE
C
          WRITE (NOUT,99999)
          DO 80 I=1,M
            WRITE (NOUT,99998) (AH(I,J),J=1,N)
   80     CONTINUE
C
          WRITE (NOUT,99997)
          CALL HYBSVD(NAH, NUH, NVH, NZ, NBH, M, N, AH, WH, MATU, UH,
     *     MATV, VH, Z, BH, IRHS, IERR, WORK)
          CALL PRINTS(WH, UH, VH, IERR, BH, NUH, NVH, NBH, MATU, MATV,
     *     M, N, NOUT, IRHS)
C
          WRITE (NOUT,99996)
          CALL SVD(NUE, M, N, AE, WE, MATU, UE, MATV, VE, IERR, WORK)
          CALL MINFIT(NUE, M, N, AE, WE, IRHS, BE, IERR, WORK)
          CALL PRINTS(WE, UE, VE, IERR, BE, NUE, NVE, NBE, MATU, MATV,
     *     M, N, NOUT, IRHS)
C
          WRITE (NOUT,99995)
          CALL SSVDC(AL, NAL, M, N, WL, E, UL, NUL, VL, NVL, WORK, JOB,
     *     IERR)
          CALL PRINTS(WL, UL, VL, IERR, E, NUL, NVL, NBH, MATU, MATV,
     *     M, N, NOUT, IRHS)
          CALL DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
     *     NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU,
     *     MATV)
   90   CONTINUE
  100 CONTINUE
      STOP
      END
      SUBROUTINE PRINTS(W, U, V, IERR, B, NAU, NV, NB, MATU, MATV, M,
     * N, NOUT, IRHS)
C
C       PRINTS OUT THE SINGULAR VALUES AND THE MATRICES U, V, B.
C
      INTEGER NAU, NV, IERR, M, N, NOUT, IRHS
      REAL U(NAU,1), V(NV,1), W(1), B(NB,1)
      LOGICAL MATU, MATV
      MINMN = MIN0(M,N)
C
      WRITE (NOUT,99999) MATU, MATV
99999 FORMAT (7H MATU =, L3, 7H MATV =, L3/)
      WRITE (NOUT,99998) IERR
99998 FORMAT (7H IERR =, I5/)
C
      WRITE (NOUT,99997)
99997 FORMAT (17H SINGULAR VALUES./)
      DO 10 I=1,MINMN
        WRITE (NOUT,99995) W(I)
   10 CONTINUE
C
      IF (.NOT.MATU) GO TO 30
      WRITE (NOUT,99996)
99996 FORMAT (/4H U ./)
99995 FORMAT (E15.8)
      DO 20 I=1,M
        WRITE (NOUT,99993) (U(I,J),J=1,MINMN)
   20 CONTINUE
C
   30 IF (.NOT.MATV) GO TO 50
      WRITE (NOUT,99994)
99994 FORMAT (/4H V ./)
99993 FORMAT (5E15.7)
      DO 40 I=1,N
        WRITE (NOUT,99993) (V(I,J),J=1,MINMN)
   40 CONTINUE
C
   50 IF (M.LT.N .OR. IRHS.LT.1) GO TO 70
      WRITE (NOUT,99992)
99992 FORMAT (/4H B ./)
      DO 60 I=1,MINMN
        WRITE (NOUT,99993) (B(I,J),J=1,IRHS)
   60 CONTINUE
C
   70 CONTINUE
      RETURN
      END
      SUBROUTINE DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
     * NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU, MATV)
C
C     COMPUTES DIFFERENCE BETWEEN W, U, V AND B AS COMPUTED BY
C     EISPACK SVD AND MINFIT, HYBSVD, AND LINPACK SSVDC.
C
      REAL WE(1), WH(1), WL(1)
      REAL UE(NUE,1), UH(NUH,1), UL(NUL,1)
      REAL VE(NVE,1), VH(NVH,1), VL(NVL,1)
      REAL BE(NBE,1), BH(NBH,1)
      INTEGER NOUT, IRHS, M, N
      LOGICAL MATU, MATV
C
C     SORT SINGULAR VALUES AND EXCHANGE COLUMNS OF U AND V FOR EISPACK.
C     SELECTION SORT MINIMIZES SWAPPING OF U AND V.
C
      IF (N.LE.1) GO TO 40
      NM1 = N - 1
      DO 30 I=1,NM1
C...    FIND INDEX OF MAXIMUM SINGULAR VALUE
        ID = I
        IP1 = I + 1
        DO 10 J=IP1,N
          IF (WE(J).GT.WE(ID)) ID = J
   10   CONTINUE
        IF (ID.EQ.I) GO TO 30
C...    SWAP SINGULAR VALUES AND VECTORS
        T = WE(I)
        WE(I) = WE(ID)
        WE(ID) = T
        IF (MATV) CALL SSWAP(N, VE(1,I), 1, VE(1,ID), 1)
        IF (MATU) CALL SSWAP(M, UE(1,I), 1, UE(1,ID), 1)
        IF (IRHS.LT.1) GO TO 30
        DO 20 KRHS=1,IRHS
          T = BE(I,KRHS)
          BE(I,KRHS) = BE(ID,KRHS)
          BE(ID,KRHS) = T
   20   CONTINUE
   30 CONTINUE
C
   40 CONTINUE
      DWHE = 0.
      DWHL = 0.
      DWEL = 0.
      DUHE = 0.
      DUHL = 0.
      DUEL = 0.
      DVHE = 0.
      DVHL = 0.
      DVEL = 0.
      DBHE = 0.
      MINMN = MIN0(M,N)
C
      DO 50 I=1,MINMN
        DWHE = DWHE + (WH(I)-WE(I))**2
        DWHL = DWHL + (WH(I)-WL(I))**2
        DWEL = DWEL + (WE(I)-WL(I))**2
   50 CONTINUE
      DWHE = SQRT(DWHE)
      DWHL = SQRT(DWHL)
      DWEL = SQRT(DWEL)
      WRITE (NOUT,99999) DWHE, DWHL, DWEL
99999 FORMAT (/35H  2-NORM ( HYBSVD W - EISPAK W ) = , 1PE15.7/6H  2-NO,
     * 29HRM ( HYBSVD W - LINPAK W ) = , 1PE15.7/19H  2-NORM ( EISPAK W,
     * 16H - LINPAK W ) = , 1PE15.7/)
C
C     DIFFERENCE OF ABSOLUTE VALUES BECAUSE OF POSSIBLE SIGN DIFFERENCE.
C
      IF (.NOT.MATU) GO TO 80
      DO 70 J=1,MINMN
        DO 60 I=1,M
          DUHE = DUHE + (ABS(UH(I,J))-ABS(UE(I,J)))**2
          DUHL = DUHL + (ABS(UH(I,J))-ABS(UL(I,J)))**2
          DUEL = DUEL + (ABS(UE(I,J))-ABS(UL(I,J)))**2
   60   CONTINUE
   70 CONTINUE
      DUHE = SQRT(DUHE)
      DUHL = SQRT(DUHL)
      DUEL = SQRT(DUEL)
      WRITE (NOUT,99998) DUHE, DUHL, DUEL
99998 FORMAT (/35H  2-NORM ( HYBSVD U - EISPAK U ) = , 1PE15.7/6H  2-NO,
     * 29HRM ( HYBSVD U - LINPAK U ) = , 1PE15.7/19H  2-NORM ( EISPAK U,
     * 16H - LINPAK U ) = , 1PE15.7/)
C
   80 IF (.NOT.MATV) GO TO 110
      DO 100 J=1,MINMN
        DO 90 I=1,N
          DVHE = DVHE + (ABS(VH(I,J))-ABS(VE(I,J)))**2
          DVHL = DVHL + (ABS(VH(I,J))-ABS(VL(I,J)))**2
          DVEL = DVEL + (ABS(VE(I,J))-ABS(VL(I,J)))**2
   90   CONTINUE
  100 CONTINUE
      DVHE = SQRT(DVHE)
      DVHL = SQRT(DVHL)
      DVEL = SQRT(DVEL)
      WRITE (NOUT,99997) DVHE, DVHL, DVEL
99997 FORMAT (/35H  2-NORM ( HYBSVD V - EISPAK V ) = , 1PE15.7/6H  2-NO,
     * 29HRM ( HYBSVD V - LINPAK V ) = , 1PE15.7/19H  2-NORM ( EISPAK V,
     * 16H - LINPAK V ) = , 1PE15.7/)
C
  110 IF (M.LT.N .OR. IRHS.LT.1) GO TO 140
      DO 130 J=1,IRHS
        DO 120 I=1,MINMN
          DBHE = DBHE + (ABS(BE(I,J))-ABS(BH(I,J)))**2
  120   CONTINUE
  130 CONTINUE
      DBHE = SQRT(DBHE)
      WRITE (NOUT,99996) DBHE
99996 FORMAT (/35H  2-NORM ( HYBSVD B - MINFIT B ) = , 1PE15.7/)
C
  140 CONTINUE
      RETURN
      END
      SUBROUTINE HYBSVD(NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV,
     * V, Z, B, IRHS, IERR, RV1)
      INTEGER NA, NU, NV, NZ, M, N, IRHS, IERR, MIN0
      REAL A(NA,1), W(1), U(NU,1), V(NV,1), Z(NZ,1), B(NB,IRHS), RV1(1)
      LOGICAL MATU, MATV
C
C     THIS ROUTINE IS A MODIFICATION OF THE GOLUB-REINSCH PROCEDURE (1)
C                                                           T
C     FOR COMPUTING THE SINGULAR VALUE DECOMPOSITION A = UWV  OF A
C     REAL M BY N RECTANGULAR MATRIX. U IS M BY MIN(M,N) CONTAINING
C     THE LEFT SINGULAR VECTORS, W IS A MIN(M,N) BY MIN(M,N) DIAGONAL
C     MATRIX CONTAINING THE SINGULAR VALUES, AND V IS N BY MIN(M,N)
C     CONTAINING THE RIGHT SINGULAR VECTORS.
C
C     THE ALGORITHM IMPLEMENTED IN THIS
C     ROUTINE HAS A HYBRID NATURE.  WHEN M IS APPROXIMATELY EQUAL TO N,
C     THE GOLUB-REINSCH ALGORITHM IS USED, BUT WHEN EITHER OF THE RATIOS
C     M/N OR N/M IS GREATER THAN ABOUT 2,
C     A MODIFIED VERSION OF THE GOLUB-REINSCH
C     ALGORITHM IS USED.  THIS MODIFIED ALGORITHM FIRST TRANSFORMS A
C                                                                T
C     INTO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS L
C     AND THEN USES THE GOLUB-REINSCH ALGORITHM TO FIND THE SINGULAR
C     VALUE DECOMPOSITION OF THE RESULTING UPPER TRIANGULAR MATRIX R.
C     WHEN U IS NEEDED EXPLICITLY IN THE CASE M.GE.N (OR V IN THE CASE
C     M.LT.N), AN EXTRA ARRAY Z (OF SIZE AT LEAST
C     MIN(M,N)**2) IS NEEDED, BUT OTHERWISE Z IS NOT REFERENCED
C     AND NO EXTRA STORAGE IS REQUIRED.  THIS HYBRID METHOD
C     SHOULD BE MORE EFFICIENT THAN THE GOLUB-REINSCH ALGORITHM WHEN
C     M/N OR N/M IS LARGE.  FOR DETAILS, SEE (2).
C
C     WHEN M .GE. N,
C     HYBSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST
C     SQUARES SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X=B.
C     IF M .LT. N (I.E. FOR UNDERDETERMINED SYSTEMS), THE RHS B
C     IS NOT PROCESSED.
C
C     NOTICE THAT THE SINGULAR VALUE DECOMPOSITION OF A MATRIX
C     IS UNIQUE ONLY UP TO THE SIGN OF THE CORRESPONDING COLUMNS
C     OF U AND V.
C
C     THIS ROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER (3) FOR
C     ADHERENCE TO A LARGE, CAREFULLY DEFINED, PORTABLE SUBSET OF
C     AMERICAN NATIONAL STANDARD FORTRAN CALLED PFORT.
C
C     REFERENCES:
C
C     (1) GOLUB,G.H. AND REINSCH,C. (1970) 'SINGULAR VALUE
C         DECOMPOSITION AND LEAST SQUARES SOLUTIONS,'
C         NUMER. MATH. 14,403-420, 1970.
C
C     (2) CHAN,T.F. (1982) 'AN IMPROVED ALGORITHM FOR COMPUTING
C         THE SINGULAR VALUE DECOMPOSITION,' ACM TOMS, VOL.8,
C         NO. 1, MARCH, 1982.
C
C     (3) RYDER,B.G. (1974) 'THE PFORT VERIFIER,' SOFTWARE -
C         PRACTICE AND EXPERIENCE, VOL.4, 359-377, 1974.
C
C     ON INPUT:
C
C        NA MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER A AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NA MUST BE AT LEAST
C          AS LARGE AS M.
C
C        NU MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY U AS DECLARED IN THE CALLING PROGRAM DIMENSION
C          STATEMENT. NU MUST BE AT LEAST AS LARGE AS M.
C
C        NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT. NV MUST BE AT LEAST AS LARGE AS N.
C
C        NZ MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER Z AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NZ MUST BE AT LEAST
C          AS LARGE AS MIN(M,N).
C
C        NB MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER B AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT. NB MUST BE AT LEAST AS LARGE AS M.
C
C        M IS THE NUMBER OF ROWS OF A (AND U).
C
C        N IS THE NUMBER OF COLUMNS OF A (AND NUMBER OF ROWS OF V).
C
C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED.
C
C        B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED
C         LINEAR SYSTEM A*X=B. IF IRHS .GT. 0 AND M .GE. N,
C         THEN ON OUTPUT, THE FIRST N COMPONENTS OF THESE IRHS COLUMNS
C                       T
C         WILL CONTAIN U B. THUS, TO COMPUTE THE MINIMAL LENGTH LEAST
C                                               +
C         SQUARES SOLUTION, ONE MUST COMPUTE V*W  TIMES THE COLUMNS OF
C                   +                        +
C         B, WHERE W  IS A DIAGONAL MATRIX, W (I)=0 IF W(I) IS
C         NEGLIGIBLE, OTHERWISE IS 1/W(I). IF IRHS=0 OR M.LT.N,
C         B IS NOT REFERENCED.
C
C        IRHS IS THE NUMBER OF RIGHT-HAND-SIDES OF THE OVERDETERMINED
C         SYSTEM A*X=B. IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULAR
C         VALUE DECOMPOSITION OF A IS DESIRED.
C
C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C        WHEN HYBSVD IS USED TO COMPUTE THE MINIMAL LENGTH LEAST
C        SQUARES SOLUTION TO AN OVERDETERMINED SYSTEM, MATU SHOULD
C        BE SET TO .FALSE. , AND MATV SHOULD BE SET TO .TRUE.  .
C
C     ON OUTPUT:
C
C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V).
C
C        W CONTAINS THE (NON-NEGATIVE) SINGULAR VALUES OF A (THE
C          DIAGONAL ELEMENTS OF W).  THEY ARE SORTED IN DESCENDING
C          ORDER.  IF AN ERROR EXIT IS MADE, THE SINGULAR VALUES
C          SHOULD BE CORRECT AND SORTED FOR INDICES IERR+1,...,MIN(M,N).
C
C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  IF MATU IS
C          FALSE, THEN U IS EITHER USED AS A TEMPORARY STORAGE (IF
C          M .GE. N) OR NOT REFERENCED (IF M .LT. N).
C          U MAY COINCIDE WITH A IN THE CALLING SEQUENCE.
C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT.
C
C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
C          MATV HAS BEEN SET TO .TRUE.  IF MATV IS
C          FALSE, THEN V IS EITHER USED AS A TEMPORARY STORAGE (IF
C          M .LT. N) OR NOT REFERENCED (IF M .GE. N).
C          IF M .GE. N, V MAY ALSO COINCIDE WITH A.  IF AN ERROR
C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
C          CORRECT SINGULAR VALUES SHOULD BE CORRECT.
C
C        Z CONTAINS THE MATRIX X IN THE SINGULAR VALUE DECOMPOSITION
C                  T
C          OF R=XSY,  IF THE MODIFIED ALGORITHM IS USED. IF THE
C          GOLUB-REINSCH PROCEDURE IS USED, THEN IT IS NOT REFERENCED.
C          IF MATU HAS BEEN SET TO .FALSE. IN THE CASE M.GE.N (OR
C          MATV SET TO .FALSE. IN THE CASE M.LT.N), THEN Z IS NOT
C          REFERENCED AND NO EXTRA STORAGE IS REQUIRED.
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS.
C          -1         IF IRHS .LT. 0 .
C          -2         IF M .LT. 1 .OR. N .LT. 1
C          -3         IF NA .LT. M .OR. NU .LT. M .OR. NB .LT. M.
C          -4         IF NV .LT. N .
C          -5         IF NZ .LT. MIN(M,N).
C
C        RV1 IS A TEMPORARY STORAGE ARRAY OF LENGTH AT LEAST MIN(M,N).
C
C     PROGRAMMED BY : TONY CHAN
C                     BOX 2158, YALE STATION,
C                     COMPUTER SCIENCE DEPT, YALE UNIV.,
C                     NEW HAVEN, CT 06520.
C     LAST MODIFIED : JANUARY, 1982.
C
C     HYBSVD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
C       INTERNAL  GRSVD, MGNSVD, SRELPR
C       FORTRAN   MIN0,ABS,SQRT,FLOAT,SIGN,AMAX1
C       BLAS      SSWAP
C
C     -----------------------------------------------------------------
C     ERROR CHECK.
C
      IERR = 0
      IF (IRHS.GE.0) GO TO 10
      IERR = -1
      RETURN
   10 IF (M.GE.1 .AND. N.GE.1) GO TO 20
      IERR = -2
      RETURN
   20 IF (NA.GE.M .AND. NU.GE.M .AND. NB.GE.M) GO TO 30
      IERR = -3
      RETURN
   30 IF (NV.GE.N) GO TO 40
      IERR = -4
      RETURN
   40 IF (NZ.GE.MIN0(M,N)) GO TO 50
      IERR = -5
      RETURN
   50 CONTINUE
C
C     FIRST COPIES A INTO EITHER U OR V ACCORDING TO WHETHER
C     M .GE. N OR M .LT. N, AND THEN CALLS SUBROUTINE MGNSVD
C     WHICH ASSUMES THAT NUMBER OF ROWS .GE. NUMBER OF COLUMNS.
C
      IF (M.LT.N) GO TO 80
C
C       M .GE. N  CASE.
C
      DO 70 I=1,M
        DO 60 J=1,N
          U(I,J) = A(I,J)
   60   CONTINUE
   70 CONTINUE
C
      CALL MGNSVD(NU, NV, NZ, NB, M, N, W, MATU, U, MATV, V, Z, B,
     * IRHS, IERR, RV1)
      RETURN
C
   80 CONTINUE
C                              T
C       M .LT. N CASE. COPIES A  INTO V.
C
      DO 100 I=1,M
        DO 90 J=1,N
          V(J,I) = A(I,J)
   90   CONTINUE
  100 CONTINUE
      CALL MGNSVD(NV, NU, NZ, NB, N, M, W, MATV, V, MATU, U, Z, B, 0,
     * IERR, RV1)
      RETURN
      END
C
      SUBROUTINE MGNSVD(NU, NV, NZ, NB, M, N, W, MATU, U, MATV, V, Z,
     * B, IRHS, IERR, RV1)
C
C     THE DESCRIPTION OF SUBROUTINE MGNSVD IS ALMOST IDENTICAL
C     TO THAT FOR SUBROUTINE HYBSVD ABOVE, WITH THE EXCEPTION
C     THAT MGNSVD ASSUMES M .GE. N.
C     IT ALSO ASSUMES THAT A COPY OF THE MATRIX A IS IN THE ARRAY U.
C
      INTEGER NU, NV, NZ, M, N, IRHS, IERR, IP1, I, J, K, IM1, IBACK
      REAL W(1), U(NU,1), V(NV,1), Z(NZ,1), B(NB,IRHS), RV1(1)
      REAL XOVRPT, C, R, G, SCALE, SIGN, ABS, SQRT, F, S, H
      REAL FLOAT
      LOGICAL MATU, MATV
C
C     SET VALUE FOR C. THE VALUE FOR C DEPENDS ON THE RELATIVE
C     EFFICIENCY OF FLOATING POINT MULTIPLICATIONS, FLOATING POINT
C     ADDITIONS AND TWO-DIMENSIONAL ARRAY INDEXINGS ON THE
C     COMPUTER WHERE THIS SUBROUTINE IS TO BE RUN.  C SHOULD
C     USUALLY BE BETWEEN 2 AND 4.  FOR DETAILS ON CHOOSING C, SEE
C     (2).  THE ALGORITHM IS NOT SENSITIVE TO THE VALUE OF C
C     ACTUALLY USED AS LONG AS C IS BETWEEN 2 AND 4.
C
      C = 4.0
C
C     DETERMINE CROSS-OVER POINT
C
      IF (MATU .AND. MATV) XOVRPT = (C+7./3.)/C
      IF (MATU .AND. .NOT.MATV) XOVRPT = (C+7./3.)/C
      IF (.NOT.MATU .AND. MATV) XOVRPT = 5./3.
      IF (.NOT.MATU .AND. .NOT.MATV) XOVRPT = 5./3.
C
C     DETERMINE WHETHER TO USE GOLUB-REINSCH OR THE MODIFIED
C     ALGORITHM.
C
      R = FLOAT(M)/FLOAT(N)
      IF (R.GE.XOVRPT) GO TO 10
C
C     USE GOLUB-REINSCH PROCEDURE
C
      CALL GRSVD(NU, NV, NB, M, N, W, MATU, U, MATV, V, B, IRHS, IERR,
     * RV1)
      GO TO 330
C
C     USE MODIFIED ALGORITHM
C
   10 CONTINUE
C
C     TRIANGULARIZE U BY HOUSEHOLDER TRANSFORMATIONS, USING
C     W AND RV1 AS TEMPORARY STORAGE.
C
      DO 110 I=1,N
        G = 0.0
        S = 0.0
        SCALE = 0.0
C
C         PERFORM SCALING OF COLUMNS TO AVOID UNNECSSARY OVERFLOW
C         OR UNDERFLOW
C
        DO 20 K=I,M
          SCALE = SCALE + ABS(U(K,I))
   20   CONTINUE
        IF (SCALE.EQ.0.0) GO TO 110
        DO 30 K=I,M
          U(K,I) = U(K,I)/SCALE
          S = S + U(K,I)**2
   30   CONTINUE
C
C         THE VECTOR E OF THE HOUSEHOLDER TRANSFORMATION I + EE'/H
C         WILL BE STORED IN COLUMN I OF U. THE TRANSFORMED ELEMENT
C         U(I,I) WILL BE STORED IN W(I) AND THE SCALAR H IN
C         RV1(I).
C
        F = U(I,I)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        U(I,I) = F - G
        RV1(I) = H
        W(I) = SCALE*G
C
        IF (I.EQ.N) GO TO 70
C
C         APPLY TRANSFORMATIONS TO REMAINING COLUMNS OF A
C
        IP1 = I + 1
        DO 60 J=IP1,N
          S = 0.0
          DO 40 K=I,M
            S = S + U(K,I)*U(K,J)
   40     CONTINUE
          F = S/H
          DO 50 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
   50     CONTINUE
   60   CONTINUE
C
C         APPLY TRANSFORMATIONS TO COLUMNS OF B IF IRHS .GT. 0
C
   70   IF (IRHS.EQ.0) GO TO 110
        DO 100 J=1,IRHS
          S = 0.0
          DO 80 K=I,M
            S = S + U(K,I)*B(K,J)
   80     CONTINUE
          F = S/H
          DO 90 K=I,M
            B(K,J) = B(K,J) + F*U(K,I)
   90     CONTINUE
  100   CONTINUE
  110 CONTINUE
C
C     COPY R INTO Z IF MATU = .TRUE.
C
      IF (.NOT.MATU) GO TO 290
      DO 130 I=1,N
        DO 120 J=I,N
          Z(J,I) = 0.0
          Z(I,J) = U(I,J)
  120   CONTINUE
        Z(I,I) = W(I)
  130 CONTINUE
C
C     ACCUMULATE HOUSEHOLDER TRANSFORMATIONS IN U
C
      DO 240 IBACK=1,N
        I = N - IBACK + 1
        IP1 = I + 1
        G = W(I)
        H = RV1(I)
        IF (I.EQ.N) GO TO 150
C
        DO 140 J=IP1,N
          U(I,J) = 0.0
  140   CONTINUE
C
  150   IF (H.EQ.0.0) GO TO 210
        IF (I.EQ.N) GO TO 190
C
        DO 180 J=IP1,N
          S = 0.0
          DO 160 K=IP1,M
            S = S + U(K,I)*U(K,J)
  160     CONTINUE
          F = S/H
          DO 170 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
  170     CONTINUE
  180   CONTINUE
C
  190   S = U(I,I)/H
        DO 200 J=I,M
          U(J,I) = U(J,I)*S
  200   CONTINUE
        GO TO 230
C
  210   DO 220 J=I,M
          U(J,I) = 0.0
  220   CONTINUE
  230   U(I,I) = U(I,I) + 1.0
  240 CONTINUE
C
C     COMPUTE SVD OF R (WHICH IS STORED IN Z)
C
      CALL GRSVD(NZ, NV, NB, N, N, W, MATU, Z, MATV, V, B, IRHS, IERR,
     * RV1)
C
C                                      T
C     FORM L*X TO OBTAIN U (WHERE R=XWY ). X IS RETURNED IN Z
C     BY GRSVD. THE MATRIX MULTIPLY IS DONE ONE ROW AT A TIME,
C     USING RV1 AS SCRATCH SPACE.
C
      DO 280 I=1,M
        DO 260 J=1,N
          S = 0.0
          DO 250 K=1,N
            S = S + U(I,K)*Z(K,J)
  250     CONTINUE
          RV1(J) = S
  260   CONTINUE
        DO 270 J=1,N
          U(I,J) = RV1(J)
  270   CONTINUE
  280 CONTINUE
      GO TO 330
C
C     FORM R IN U BY ZEROING THE LOWER TRIANGULAR PART OF R IN U
C
  290 IF (N.EQ.1) GO TO 320
      DO 310 I=2,N
        IM1 = I - 1
        DO 300 J=1,IM1
          U(I,J) = 0.0
  300   CONTINUE
        U(I,I) = W(I)
  310 CONTINUE
  320 U(1,1) = W(1)
C
      CALL GRSVD(NU, NV, NB, N, N, W, MATU, U, MATV, V, B, IRHS, IERR,
     * RV1)
  330 CONTINUE
      IERRP1 = IERR + 1
      IF (IERR.LT.0 .OR. N.LE.1 .OR. IERRP1.EQ.N) RETURN
C
C     SORT SINGULAR VALUES AND EXCHANGE COLUMNS OF U AND V ACCORDINGLY.
C     SELECTION SORT MINIMIZES SWAPPING OF U AND V.
C
      NM1 = N - 1
      DO 360 I=IERRP1,NM1
C...    FIND INDEX OF MAXIMUM SINGULAR VALUE
        ID = I
        IP1 = I + 1
        DO 340 J=IP1,N
          IF (W(J).GT.W(ID)) ID = J
  340   CONTINUE
        IF (ID.EQ.I) GO TO 360
C...    SWAP SINGULAR VALUES AND VECTORS
        T = W(I)
        W(I) = W(ID)
        W(ID) = T
        IF (MATV) CALL SSWAP(N, V(1,I), 1, V(1,ID), 1)
        IF (MATU) CALL SSWAP(M, U(1,I), 1, U(1,ID), 1)
        IF (IRHS.LT.1) GO TO 360
        DO 350 KRHS=1,IRHS
          T = B(I,KRHS)
          B(I,KRHS) = B(ID,KRHS)
          B(ID,KRHS) = T
  350   CONTINUE
  360 CONTINUE
      RETURN
C     ************** LAST CARD OF HYBSVD *****************
      END
      SUBROUTINE GRSVD(NU, NV, NB, M, N, W, MATU, U, MATV, V, B, IRHS,
     * IERR, RV1)
C
      INTEGER I, J, K, L, M, N, II, I1, KK, K1, LL, L1, MN, NU, NV, NB,
     * ITS, IERR, IRHS
      REAL W(1), U(NU,1), V(NV,1), B(NB,IRHS), RV1(1)
      REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, SRELPR, DUMMY
      REAL SQRT, AMAX1, ABS, SIGN
      LOGICAL MATU, MATV
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD,
C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
C
C     THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION
C          T
C     A=USV  OF A REAL M BY N RECTANGULAR MATRIX.  HOUSEHOLDER
C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
C     GRSVD ASSUMES THAT A COPY OF THE MATRIX A IS IN THE ARRAY U. IT
C     ALSO ASSUMES M .GE. N.  IF M .LT. N, THEN COMPUTE THE SINGULAR
C                             T       T    T             T
C     VALUE DECOMPOSITION OF A .  IF A =UWV  , THEN A=VWU  .
C
C     GRSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST SQUARES
C     SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X=B.
C
C     ON INPUT-
C
C        NU MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS U AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NU MUST BE AT LEAST
C          AS LARGE AS M,
C
C        NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
C          ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NV MUST BE AT LEAST AS LARGE AS N,
C
C        NB MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS B AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NB MUST BE AT LEAST
C          AS LARGE AS M,
C
C        M IS THE NUMBER OF ROWS OF A (AND U),
C
C        N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V,
C
C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED,
C
C        B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED
C          LINEAR SYSTEM A*X=B.  IF IRHS .GT. 0,  THEN ON OUTPUT,
C          THE FIRST N COMPONENTS OF THESE IRHS COLUMNS OF B
C                        T
C          WILL CONTAIN U B.  THUS, TO COMPUTE THE MINIMAL LENGTH LEAST
C                                                +
C          SQUARES SOLUTION, ONE MUST COMPUTE V*W  TIMES THE COLUMNS OF
C                    +                        +
C          B, WHERE W  IS A DIAGONAL MATRIX, W (I)=0 IF W(I) IS
C          NEGLIGIBLE, OTHERWISE IS 1/W(I).  IF IRHS=0, B MAY COINCIDE
C          WITH A OR U AND WILL NOT BE REFERENCED,
C
C        IRHS IS THE NUMBER OF RIGHT-HAND-SIDES OF THE OVERDETERMINED
C          SYSTEM A*X=B.  IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULA
C          VALUE DECOMPOSITION OF A IS DESIRED,
C
C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE,
C
C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C     ON OUTPUT-
C
C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
C          FOR INDICES IERR+1,IERR+2,...,N,
C
C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  OTHERWISE
C          U IS USED AS A TEMPORARY ARRAY.
C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
C
C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
C          MATV HAS BEEN SET TO .TRUE.  OTHERWISE V IS NOT REFERENCED.
C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO
C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS,
C          -1         IF IRHS .LT. 0 ,
C          -2         IF M .LT. N ,
C          -3         IF NU .LT. M .OR. NB .LT. M,
C          -4         IF NV .LT. N .
C
C        RV1 IS A TEMPORARY STORAGE ARRAY.
C
C        THIS SUBROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER
C        (RYDER, B.G. 'THE PFORT VERIFIER', SOFTWARE - PRACTICE AND
C        EXPERIENCE, VOL.4, 359-377, 1974) FOR ADHERENCE TO A LARGE,
C        CAREFULLY DEFINED, PORTABLE SUBSET OF AMERICAN NATIONAL STANDAR
C        FORTRAN CALLED PFORT.
C
C        ORIGINAL VERSION OF THIS CODE IS SUBROUTINE SVD IN RELEASE 2 OF
C        EISPACK.
C
C        MODIFIED BY TONY F. CHAN,
C                    COMP. SCI. DEPT, YALE UNIV.,
C                    BOX 2158, YALE STATION,
C                    CT 06520
C        LAST MODIFIED : JANUARY, 1982.
C
C     ------------------------------------------------------------------
C
C     ********** SRELPR IS A MACHINE-DEPENDENT FUNCTION SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
C
      IERR = 0
      IF (IRHS.GE.0) GO TO 10
      IERR = -1
      RETURN
   10 IF (M.GE.N) GO TO 20
      IERR = -2
      RETURN
   20 IF (NU.GE.M .AND. NB.GE.M) GO TO 30
      IERR = -3
      RETURN
   30 IF (NV.GE.N) GO TO 40
      IERR = -4
      RETURN
   40 CONTINUE
C
C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
      G = 0.0
      SCALE = 0.0
      X = 0.0
C
      DO 260 I=1,N
        L = I + 1
        RV1(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
C
C     COMPUTE LEFT TRANSFORMATIONS THAT ZERO THE SUBDIAGONAL ELEMENTS
C     OF THE I-TH COLUMN.
C
        DO 50 K=I,M
          SCALE = SCALE + ABS(U(K,I))
   50   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 160
C
        DO 60 K=I,M
          U(K,I) = U(K,I)/SCALE
          S = S + U(K,I)**2
   60   CONTINUE
C
        F = U(I,I)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        U(I,I) = F - G
        IF (I.EQ.N) GO TO 100
C
C     APPLY LEFT TRANSFORMATIONS TO REMAINING COLUMNS OF A.
C
        DO 90 J=L,N
          S = 0.0
C
          DO 70 K=I,M
            S = S + U(K,I)*U(K,J)
   70     CONTINUE
C
          F = S/H
C
          DO 80 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
   80     CONTINUE
   90   CONTINUE
C
C      APPLY LEFT TRANSFORMATIONS TO THE COLUMNS OF B IF IRHS .GT. 0
C
  100   IF (IRHS.EQ.0) GO TO 140
        DO 130 J=1,IRHS
          S = 0.0
          DO 110 K=I,M
            S = S + U(K,I)*B(K,J)
  110     CONTINUE
          F = S/H
          DO 120 K=I,M
            B(K,J) = B(K,J) + F*U(K,I)
  120     CONTINUE
  130   CONTINUE
C
C     COMPUTE RIGHT TRANSFORMATIONS.
C
  140   DO 150 K=I,M
          U(K,I) = SCALE*U(K,I)
  150   CONTINUE
C
  160   W(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
        IF (I.GT.M .OR. I.EQ.N) GO TO 250
C
        DO 170 K=L,N
          SCALE = SCALE + ABS(U(I,K))
  170   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 250
C
        DO 180 K=L,N
          U(I,K) = U(I,K)/SCALE
          S = S + U(I,K)**2
  180   CONTINUE
C
        F = U(I,L)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        U(I,L) = F - G
C
        DO 190 K=L,N
          RV1(K) = U(I,K)/H
  190   CONTINUE
C
        IF (I.EQ.M) GO TO 230
C
        DO 220 J=L,M
          S = 0.0
C
          DO 200 K=L,N
            S = S + U(J,K)*U(I,K)
  200     CONTINUE
C
          DO 210 K=L,N
            U(J,K) = U(J,K) + S*RV1(K)
  210     CONTINUE
  220   CONTINUE
C
  230   DO 240 K=L,N
          U(I,K) = SCALE*U(I,K)
  240   CONTINUE
C
  250   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  260 CONTINUE
C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS **********
      IF (.NOT.MATV) GO TO 350
C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
      DO 340 II=1,N
        I = N + 1 - II
        IF (I.EQ.N) GO TO 330
        IF (G.EQ.0.0) GO TO 310
C
        DO 270 J=L,N
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
          V(J,I) = (U(I,J)/U(I,L))/G
  270   CONTINUE
C
        DO 300 J=L,N
          S = 0.0
C
          DO 280 K=L,N
            S = S + U(I,K)*V(K,J)
  280     CONTINUE
C
          DO 290 K=L,N
            V(K,J) = V(K,J) + S*V(K,I)
  290     CONTINUE
  300   CONTINUE
C
  310   DO 320 J=L,N
          V(I,J) = 0.0
          V(J,I) = 0.0
  320   CONTINUE
C
  330   V(I,I) = 1.0
        G = RV1(I)
        L = I
  340 CONTINUE
C     ********** ACCUMULATION OF LEFT-HAND TRANSFORMATIONS **********
  350 IF (.NOT.MATU) GO TO 470
C     **********FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- **********
      MN = N
      IF (M.LT.N) MN = M
C
      DO 460 II=1,MN
        I = MN + 1 - II
        L = I + 1
        G = W(I)
        IF (I.EQ.N) GO TO 370
C
        DO 360 J=L,N
          U(I,J) = 0.0
  360   CONTINUE
C
  370   IF (G.EQ.0.0) GO TO 430
        IF (I.EQ.MN) GO TO 410
C
        DO 400 J=L,N
          S = 0.0
C
          DO 380 K=L,M
            S = S + U(K,I)*U(K,J)
  380     CONTINUE
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
          F = (S/U(I,I))/G
C
          DO 390 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
  390     CONTINUE
  400   CONTINUE
C
  410   DO 420 J=I,M
          U(J,I) = U(J,I)/G
  420   CONTINUE
C
        GO TO 450
C
  430   DO 440 J=I,M
          U(J,I) = 0.0
  440   CONTINUE
C
  450   U(I,I) = U(I,I) + 1.0
  460 CONTINUE
C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  470 EPS = SRELPR(DUMMY)*X
C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
      DO 650 KK=1,N
        K1 = N - KK
        K = K1 + 1
        ITS = 0
C     ********** TEST FOR SPLITTING.
C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  480   DO 490 LL=1,K
          L1 = K - LL
          L = L1 + 1
          IF (ABS(RV1(L)).LE.EPS) GO TO 550
C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
          IF (ABS(W(L1)).LE.EPS) GO TO 500
  490   CONTINUE
C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  500   C = 0.0
        S = 1.0
C
        DO 540 I=L,K
          F = S*RV1(I)
          RV1(I) = C*RV1(I)
          IF (ABS(F).LE.EPS) GO TO 550
          G = W(I)
          H = SQRT(F*F+G*G)
          W(I) = H
          C = G/H
          S = -F/H
C
C     APPLY LEFT TRANSFORMATIONS TO B IF IRHS .GT. 0
C
          IF (IRHS.EQ.0) GO TO 520
          DO 510 J=1,IRHS
            Y = B(L1,J)
            Z = B(I,J)
            B(L1,J) = Y*C + Z*S
            B(I,J) = -Y*S + Z*C
  510     CONTINUE
  520     CONTINUE
C
          IF (.NOT.MATU) GO TO 540
C
          DO 530 J=1,M
            Y = U(J,L1)
            Z = U(J,I)
            U(J,L1) = Y*C + Z*S
            U(J,I) = -Y*S + Z*C
  530     CONTINUE
C
  540   CONTINUE
C     ********** TEST FOR CONVERGENCE **********
  550   Z = W(K)
        IF (L.EQ.K) GO TO 630
C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
        IF (ITS.EQ.30) GO TO 660
        ITS = ITS + 1
        X = W(L)
        Y = W(K1)
        G = RV1(K1)
        H = RV1(K)
        F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
        G = SQRT(F*F+1.0)
        F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C     ********** NEXT QR TRANSFORMATION **********
        C = 1.0
        S = 1.0
C
        DO 620 I1=L,K1
          I = I1 + 1
          G = RV1(I)
          Y = W(I)
          H = S*G
          G = C*G
          Z = SQRT(F*F+H*H)
          RV1(I1) = Z
          C = F/Z
          S = H/Z
          F = X*C + G*S
          G = -X*S + G*C
          H = Y*S
          Y = Y*C
          IF (.NOT.MATV) GO TO 570
C
          DO 560 J=1,N
            X = V(J,I1)
            Z = V(J,I)
            V(J,I1) = X*C + Z*S
            V(J,I) = -X*S + Z*C
  560     CONTINUE
C
  570     Z = SQRT(F*F+H*H)
          W(I1) = Z
C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
          IF (Z.EQ.0.0) GO TO 580
          C = F/Z
          S = H/Z
  580     F = C*G + S*Y
          X = -S*G + C*Y
C
C     APPLY LEFT TRANSFORMATIONS TO B IF IRHS .GT. 0
C
          IF (IRHS.EQ.0) GO TO 600
          DO 590 J=1,IRHS
            Y = B(I1,J)
            Z = B(I,J)
            B(I1,J) = Y*C + Z*S
            B(I,J) = -Y*S + Z*C
  590     CONTINUE
  600     CONTINUE
C
          IF (.NOT.MATU) GO TO 620
C
          DO 610 J=1,M
            Y = U(J,I1)
            Z = U(J,I)
            U(J,I1) = Y*C + Z*S
            U(J,I) = -Y*S + Z*C
  610     CONTINUE
C
  620   CONTINUE
C
        RV1(L) = 0.0
        RV1(K) = F
        W(K) = X
        GO TO 480
C     ********** CONVERGENCE **********
  630   IF (Z.GE.0.0) GO TO 650
C     ********** W(K) IS MADE NON-NEGATIVE **********
        W(K) = -Z
        IF (.NOT.MATV) GO TO 650
C
        DO 640 J=1,N
          V(J,K) = -V(J,K)
  640   CONTINUE
C
  650 CONTINUE
C
      GO TO 670
C     ********** SET ERROR -- NO CONVERGENCE TO A
C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  660 IERR = K
  670 RETURN
C     ********** LAST CARD OF GRSVD **********
      END
      SUBROUTINE SSWAP(N, SX, INCX, SY, INCY)
C
C     INTERCHANGES TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1), SY(1), STEMP
      INTEGER I, INCX, INCY, IX, IY, M, MP1, N
C
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I=1,N
        STEMP = SX(IX)
        SX(IX) = SY(IY)
        SY(IY) = STEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C       CLEAN-UP LOOP
C
   20 M = MOD(N,3)
      IF (M.EQ.0) GO TO 40
      DO 30 I=1,M
        STEMP = SX(I)
        SX(I) = SY(I)
        SY(I) = STEMP
   30 CONTINUE
      IF (N.LT.3) RETURN
   40 MP1 = M + 1
      DO 50 I=MP1,N,3
        STEMP = SX(I)
        SX(I) = SY(I)
        SY(I) = STEMP
        STEMP = SX(I+1)
        SX(I+1) = SY(I+1)
        SY(I+1) = STEMP
        STEMP = SX(I+2)
        SX(I+2) = SY(I+2)
        SY(I+2) = STEMP
   50 CONTINUE
      RETURN
      END
      REAL FUNCTION SRELPR(DUMMY)
      REAL DUMMY
C
C     SRELPR COMPUTES THE RELATIVE PRECISION OF THE FLOATING POINT
C     ARITHMETIC OF THE MACHINE.
C
C     IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES,
C     THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS.
C     ASSUME THE COMPUTER HAS
C
C        B = BASE OF ARITHMETIC
C        T = NUMBER OF BASE  B  DIGITS
C
C     THEN
C
C        SRELPR = B**(1-T)
C
      REAL S
C
      SRELPR = 1.0
   10 SRELPR = SRELPR/2.0
      S = 1.0 + SRELPR
      IF (S.GT.1.0) GO TO 10
      SRELPR = 2.0*SRELPR
      RETURN
      END
      SUBROUTINE SSVDC(X, LDX, N, P, S, E, U, LDU, V, LDV, WORK, JOB,
     * INFO)
      INTEGER LDX, N, P, LDU, LDV, JOB, INFO
      REAL X(LDX,1), S(1), E(1), U(LDU,1), V(LDV,1), WORK(1)
C
C
C     SSVDC IS A SUBROUTINE TO REDUCE A REAL NXP MATRIX X BY
C     ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
C
C     ON ENTRY
C
C         X         REAL(LDX,P), WHERE LDX.GE.N.
C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
C                   DESTROYED BY SSVDC.
C
C         LDX       INTEGER.
C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
C
C         N         INTEGER.
C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
C
C         P         INTEGER.
C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
C
C         LDU       INTEGER.
C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
C                   (SEE BELOW).
C
C         LDV       INTEGER.
C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
C                   (SEE BELOW).
C
C         WORK      REAL(N).
C                   WORK IS A SCRATCH ARRAY.
C
C         JOB       INTEGER.
C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
C                   WITH THE FOLLOWING MEANING
C
C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
C                                  VECTORS.
C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
C                                  IN U.
C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
C                                  VECTORS IN U.
C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
C                                  VECTORS.
C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
C                                  IN V.
C
C     ON RETURN
C
C         S         REAL(MM), WHERE MM=MIN(N+1,P).
C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
C                   ORDER OF MAGNITUDE.
C
C         E         REAL(P).
C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
C                   DISCUSSION OF INFO FOR EXCEPTIONS.
C
C         U         REAL(LDU,K), WHERE LDU.GE.N.  IF JOBA.EQ.1 THEN
C                                   K.EQ.N, IF JOBA.GE.2 THEN
C                                   K.EQ.MIN(N,P).
C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
C                   IN THE SUBROUTINE CALL.
C
C         V         REAL(LDV,P), WHERE LDV.GE.P.
C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
C                   THEN V MAY BE IDENTIFIED WITH X IN THE
C                   SUBROUTINE CALL.
C
C         INFO      INTEGER.
C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
C                   VALUES OF X AND B ARE THE SAME.
C
C     LINPACK. THIS VERSION DATED 03/19/79 .
C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
C
C     ***** USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
C
C     EXTERNAL SROT
C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2,SROTG
C     FORTRAN ABS,AMAX1,MAX0,MIN0,MOD,SQRT
C
C     INTERNAL VARIABLES
C
      INTEGER I, ITER, J, JOBU, K, KASE, KK, L, LL, LLS, LM1, LP1, LS,
     * LU, M, MAXIT, MM, MM1, MP1, NCT, NCTP1, NCU, NRT, NRTP1
      REAL SDOT, T, R
      REAL B, C, CS, EL, EMM1, F, G, SNRM2, SCALE, SHIFT, SL, SM, SN,
     * SMM1, T1, TEST, ZTEST
      LOGICAL WANTU, WANTV
C
C
C     SET THE MAXIMUM NUMBER OF ITERATIONS.
C
      MAXIT = 30
C
C     DETERMINE WHAT IS TO BE COMPUTED.
C
      WANTU = .FALSE.
      WANTV = .FALSE.
      JOBU = MOD(JOB,100)/10
      NCU = N
      IF (JOBU.GT.1) NCU = MIN0(N,P)
      IF (JOBU.NE.0) WANTU = .TRUE.
      IF (MOD(JOB,10).NE.0) WANTV = .TRUE.
C
C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
C
      INFO = 0
      NCT = MIN0(N-1,P)
      NRT = MAX0(0,MIN0(P-2,N))
      LU = MAX0(NCT,NRT)
      IF (LU.LT.1) GO TO 170
      DO 160 L=1,LU
        LP1 = L + 1
        IF (L.GT.NCT) GO TO 20
C
C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
C           PLACE THE L-TH DIAGONAL IN S(L).
C
        S(L) = SNRM2(N-L+1,X(L,L),1)
        IF (S(L).EQ.0.0E0) GO TO 10
        IF (X(L,L).NE.0.0E0) S(L) = SIGN(S(L),X(L,L))
        CALL SSCAL(N-L+1, 1.0E0/S(L), X(L,L), 1)
        X(L,L) = 1.0E0 + X(L,L)
   10   CONTINUE
        S(L) = -S(L)
   20   CONTINUE
        IF (P.LT.LP1) GO TO 50
        DO 40 J=LP1,P
          IF (L.GT.NCT) GO TO 30
          IF (S(L).EQ.0.0E0) GO TO 30
C
C              APPLY THE TRANSFORMATION.
C
          T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
          CALL SAXPY(N-L+1, T, X(L,L), 1, X(L,J), 1)
   30     CONTINUE
C
C           PLACE THE L-TH ROW OF X INTO  E FOR THE
C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
C
          E(J) = X(L,J)
   40   CONTINUE
   50   CONTINUE
        IF (.NOT.WANTU .OR. L.GT.NCT) GO TO 70
C
C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
C           MULTIPLICATION.
C
        DO 60 I=L,N
          U(I,L) = X(I,L)
   60   CONTINUE
   70   CONTINUE
        IF (L.GT.NRT) GO TO 150
C
C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
C           L-TH SUPER-DIAGONAL IN E(L).
C
        E(L) = SNRM2(P-L,E(LP1),1)
        IF (E(L).EQ.0.0E0) GO TO 80
        IF (E(LP1).NE.0.0E0) E(L) = SIGN(E(L),E(LP1))
        CALL SSCAL(P-L, 1.0E0/E(L), E(LP1), 1)
        E(LP1) = 1.0E0 + E(LP1)
   80   CONTINUE
        E(L) = -E(L)
        IF (LP1.GT.N .OR. E(L).EQ.0.0E0) GO TO 120
C
C              APPLY THE TRANSFORMATION.
C
        DO 90 I=LP1,N
          WORK(I) = 0.0E0
   90   CONTINUE
        DO 100 J=LP1,P
          CALL SAXPY(N-L, E(J), X(LP1,J), 1, WORK(LP1), 1)
  100   CONTINUE
        DO 110 J=LP1,P
          CALL SAXPY(N-L, -E(J)/E(LP1), WORK(LP1), 1, X(LP1,J), 1)
  110   CONTINUE
  120   CONTINUE
        IF (.NOT.WANTV) GO TO 140
C
C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
C              BACK MULTIPLICATION.
C
        DO 130 I=LP1,P
          V(I,L) = E(I)
  130   CONTINUE
  140   CONTINUE
  150   CONTINUE
  160 CONTINUE
  170 CONTINUE
C
C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
C
      M = MIN0(P,N+1)
      NCTP1 = NCT + 1
      NRTP1 = NRT + 1
      IF (NCT.LT.P) S(NCTP1) = X(NCTP1,NCTP1)
      IF (N.LT.M) S(M) = 0.0E0
      IF (NRTP1.LT.M) E(NRTP1) = X(NRTP1,M)
      E(M) = 0.0E0
C
C     IF REQUIRED, GENERATE U.
C
      IF (.NOT.WANTU) GO TO 300
      IF (NCU.LT.NCTP1) GO TO 200
      DO 190 J=NCTP1,NCU
        DO 180 I=1,N
          U(I,J) = 0.0E0
  180   CONTINUE
        U(J,J) = 1.0E0
  190 CONTINUE
  200 CONTINUE
      IF (NCT.LT.1) GO TO 290
      DO 280 LL=1,NCT
        L = NCT - LL + 1
        IF (S(L).EQ.0.0E0) GO TO 250
        LP1 = L + 1
        IF (NCU.LT.LP1) GO TO 220
        DO 210 J=LP1,NCU
          T = -SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
          CALL SAXPY(N-L+1, T, U(L,L), 1, U(L,J), 1)
  210   CONTINUE
  220   CONTINUE
        CALL SSCAL(N-L+1, -1.0E0, U(L,L), 1)
        U(L,L) = 1.0E0 + U(L,L)
        LM1 = L - 1
        IF (LM1.LT.1) GO TO 240
        DO 230 I=1,LM1
          U(I,L) = 0.0E0
  230   CONTINUE
  240   CONTINUE
        GO TO 270
  250   CONTINUE
        DO 260 I=1,N
          U(I,L) = 0.0E0
  260   CONTINUE
        U(L,L) = 1.0E0
  270   CONTINUE
  280 CONTINUE
  290 CONTINUE
  300 CONTINUE
C
C     IF IT IS REQUIRED, GENERATE V.
C
      IF (.NOT.WANTV) GO TO 350
      DO 340 LL=1,P
        L = P - LL + 1
        LP1 = L + 1
        IF (L.GT.NRT) GO TO 320
        IF (E(L).EQ.0.0E0) GO TO 320
        DO 310 J=LP1,P
          T = -SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
          CALL SAXPY(P-L, T, V(LP1,L), 1, V(LP1,J), 1)
  310   CONTINUE
  320   CONTINUE
        DO 330 I=1,P
          V(I,L) = 0.0E0
  330   CONTINUE
        V(L,L) = 1.0E0
  340 CONTINUE
  350 CONTINUE
C
C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
C
      MM = M
      ITER = 0
  360 CONTINUE
C
C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
C
C     ...EXIT
      IF (M.EQ.0) GO TO 620
C
C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
C        FLAG AND RETURN.
C
      IF (ITER.LT.MAXIT) GO TO 370
      INFO = M
C     ......EXIT
      GO TO 620
  370 CONTINUE
C
C        THIS SECTION OF THE PROGRAM INSPECTS FOR
C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
C
C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
C
      DO 390 LL=1,M
        L = M - LL
C        ...EXIT
        IF (L.EQ.0) GO TO 400
        TEST = ABS(S(L)) + ABS(S(L+1))
        ZTEST = TEST + ABS(E(L))
        IF (ZTEST.NE.TEST) GO TO 380
        E(L) = 0.0E0
C        ......EXIT
        GO TO 400
  380   CONTINUE
  390 CONTINUE
  400 CONTINUE
      IF (L.NE.M-1) GO TO 410
      KASE = 4
      GO TO 480
  410 CONTINUE
      LP1 = L + 1
      MP1 = M + 1
      DO 430 LLS=LP1,MP1
        LS = M - LLS + LP1
C           ...EXIT
        IF (LS.EQ.L) GO TO 440
        TEST = 0.0E0
        IF (LS.NE.M) TEST = TEST + ABS(E(LS))
        IF (LS.NE.L+1) TEST = TEST + ABS(E(LS-1))
        ZTEST = TEST + ABS(S(LS))
        IF (ZTEST.NE.TEST) GO TO 420
        S(LS) = 0.0E0
C           ......EXIT
        GO TO 440
  420   CONTINUE
  430 CONTINUE
  440 CONTINUE
      IF (LS.NE.L) GO TO 450
      KASE = 3
      GO TO 470
  450 CONTINUE
      IF (LS.NE.M) GO TO 460
      KASE = 1
      GO TO 470
  460 CONTINUE
      KASE = 2
      L = LS
  470 CONTINUE
  480 CONTINUE
      L = L + 1
C
C        PERFORM THE TASK INDICATED BY KASE.
C
      GO TO (490, 520, 540, 570), KASE
C
C        DEFLATE NEGLIGIBLE S(M).
C
  490 CONTINUE
      MM1 = M - 1
      F = E(M-1)
      E(M-1) = 0.0E0
      DO 510 KK=L,MM1
        K = MM1 - KK + L
        T1 = S(K)
        CALL SROTG(T1, F, CS, SN)
        S(K) = T1
        IF (K.EQ.L) GO TO 500
        F = -SN*E(K-1)
        E(K-1) = CS*E(K-1)
  500   CONTINUE
        IF (WANTV) CALL SROT(P, V(1,K), 1, V(1,M), 1, CS, SN)
  510 CONTINUE
      GO TO 610
C
C        SPLIT AT NEGLIGIBLE S(L).
C
  520 CONTINUE
      F = E(L-1)
      E(L-1) = 0.0E0
      DO 530 K=L,M
        T1 = S(K)
        CALL SROTG(T1, F, CS, SN)
        S(K) = T1
        F = -SN*E(K)
        E(K) = CS*E(K)
        IF (WANTU) CALL SROT(N, U(1,K), 1, U(1,L-1), 1, CS, SN)
  530 CONTINUE
      GO TO 610
C
C        PERFORM ONE QR STEP.
C
  540 CONTINUE
C
C           CALCULATE THE SHIFT.
C
      SCALE = AMAX1(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),ABS(E(L)
     * ))
      SM = S(M)/SCALE
      SMM1 = S(M-1)/SCALE
      EMM1 = E(M-1)/SCALE
      SL = S(L)/SCALE
      EL = E(L)/SCALE
      B = ((SMM1+SM)*(SMM1-SM)+EMM1**2)/2.0E0
      C = (SM*EMM1)**2
      SHIFT = 0.0E0
      IF (B.EQ.0.0E0 .AND. C.EQ.0.0E0) GO TO 550
      SHIFT = SQRT(B**2+C)
      IF (B.LT.0.0E0) SHIFT = -SHIFT
      SHIFT = C/(B+SHIFT)
  550 CONTINUE
      F = (SL+SM)*(SL-SM) - SHIFT
      G = SL*EL
C
C           CHASE ZEROS.
C
      MM1 = M - 1
      DO 560 K=L,MM1
        CALL SROTG(F, G, CS, SN)
        IF (K.NE.L) E(K-1) = F
        F = CS*S(K) + SN*E(K)
        E(K) = CS*E(K) - SN*S(K)
        G = SN*S(K+1)
        S(K+1) = CS*S(K+1)
        IF (WANTV) CALL SROT(P, V(1,K), 1, V(1,K+1), 1, CS, SN)
        CALL SROTG(F, G, CS, SN)
        S(K) = F
        F = CS*E(K) + SN*S(K+1)
        S(K+1) = -SN*E(K) + CS*S(K+1)
        G = SN*E(K+1)
        E(K+1) = CS*E(K+1)
        IF (WANTU .AND. K.LT.N) CALL SROT(N, U(1,K), 1, U(1,K+1), 1,
     *   CS, SN)
  560 CONTINUE
      E(M-1) = F
      ITER = ITER + 1
      GO TO 610
C
C        CONVERGENCE.
C
  570 CONTINUE
C
C           MAKE THE SINGULAR VALUE  POSITIVE.
C
      IF (S(L).GE.0.0E0) GO TO 580
      S(L) = -S(L)
      IF (WANTV) CALL SSCAL(P, -1.0E0, V(1,L), 1)
  580 CONTINUE
C
C           ORDER THE SINGULAR VALUE.
C
  590 IF (L.EQ.MM) GO TO 600
C           ...EXIT
      IF (S(L).GE.S(L+1)) GO TO 600
      T = S(L)
      S(L) = S(L+1)
      S(L+1) = T
      IF (WANTV .AND. L.LT.P) CALL SSWAP(P, V(1,L), 1, V(1,L+1), 1)
      IF (WANTU .AND. L.LT.N) CALL SSWAP(N, U(1,L), 1, U(1,L+1), 1)
      L = L + 1
      GO TO 590
  600 CONTINUE
      ITER = 0
      M = M - 1
  610 CONTINUE
      GO TO 360
  620 CONTINUE
      RETURN
      END
      REAL FUNCTION SASUM(N, SX, INCX)
C
C     TAKES THE SUM OF THE ABSOLUTE VALUES.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1), STEMP
      INTEGER I, INCX, M, MP1, N, NINCX
C
      SASUM = 0.0E0
      STEMP = 0.0E0
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1) GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I=1,NINCX,INCX
        STEMP = STEMP + ABS(SX(I))
   10 CONTINUE
      SASUM = STEMP
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,6)
      IF (M.EQ.0) GO TO 40
      DO 30 I=1,M
        STEMP = STEMP + ABS(SX(I))
   30 CONTINUE
      IF (N.LT.6) GO TO 60
   40 MP1 = M + 1
      DO 50 I=MP1,N,6
        STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) + ABS(SX(I+2)) +
     *   ABS(SX(I+3)) + ABS(SX(I+4)) + ABS(SX(I+5))
   50 CONTINUE
   60 SASUM = STEMP
      RETURN
      END
      SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY)
C
C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
C     USES UNROLLED LOOP FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1), SY(1), SA
      INTEGER I, INCX, INCY, IX, IY, M, MP1, N
C
      IF (N.LE.0) RETURN
      IF (SA.EQ.0.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I=1,N
        SY(IY) = SY(IY) + SA*SX(IX)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,4)
      IF (M.EQ.0) GO TO 40
      DO 30 I=1,M
        SY(I) = SY(I) + SA*SX(I)
   30 CONTINUE
      IF (N.LT.4) RETURN
   40 MP1 = M + 1
      DO 50 I=MP1,N,4
        SY(I) = SY(I) + SA*SX(I)
        SY(I+1) = SY(I+1) + SA*SX(I+1)
        SY(I+2) = SY(I+2) + SA*SX(I+2)
        SY(I+3) = SY(I+3) + SA*SX(I+3)
   50 CONTINUE
      RETURN
      END
      REAL FUNCTION SDOT(N, SX, INCX, SY, INCY)
C
C     FORMS THE DOT PRODUCT OF TWO VECTORS.
C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1), SY(1), STEMP
      INTEGER I, INCX, INCY, IX, IY, M, MP1, N
C
      STEMP = 0.0E0
      SDOT = 0.0E0
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
C
C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
C          NOT EQUAL TO 1
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I=1,N
        STEMP = STEMP + SX(IX)*SY(IY)
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      SDOT = STEMP
      RETURN
C
C        CODE FOR BOTH INCREMENTS EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF (M.EQ.0) GO TO 40
      DO 30 I=1,M
        STEMP = STEMP + SX(I)*SY(I)
   30 CONTINUE
      IF (N.LT.5) GO TO 60
   40 MP1 = M + 1
      DO 50 I=MP1,N,5
        STEMP = STEMP + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2)
     *   + SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4)
   50 CONTINUE
   60 SDOT = STEMP
      RETURN
      END
      REAL FUNCTION SNRM2(N, SX, INCX)
      INTEGER NEXT
      REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
      DATA ZERO, ONE /0.0E0,1.0E0/
C
C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE
C     INCREMENT INCX .
C     IF    N .LE. 0 RETURN WITH RESULT = 0.
C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
C
C           C.L.LAWSON, 1978 JAN 08
C
C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
C     HOPEFULLY APPLICABLE TO ALL MACHINES.
C         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
C         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
C     WHERE
C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
C
C     BRIEF OUTLINE OF ALGORITHM..
C
C     PHASE 1    SCANS ZERO COMPONENTS.
C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
C
C     VALUES FOR CUTLO AND CUTHI..
C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
C                   UNIVAC AND DEC AT 2**(-103)
C                   THUS CUTLO = 2**(-51) = 4.44089E-16
C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
C                   THUS CUTHI = 2**(63.5) = 1.30438E19
C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
      DATA CUTLO, CUTHI /4.441E-16,1.304E19/
C
      IF (N.GT.0) GO TO 10
      SNRM2 = ZERO
      GO TO 140
C
   10 ASSIGN 30 TO NEXT
      SUM = ZERO
      NN = N*INCX
C                                                 BEGIN MAIN LOOP
      I = 1
   20 GO TO NEXT, (30, 40, 70, 80)
   30 IF (ABS(SX(I)).GT.CUTLO) GO TO 110
      ASSIGN 40 TO NEXT
      XMAX = ZERO
C
C                        PHASE 1.  SUM IS ZERO
C
   40 IF (SX(I).EQ.ZERO) GO TO 130
      IF (ABS(SX(I)).GT.CUTLO) GO TO 110
C
C                                PREPARE FOR PHASE 2.
      ASSIGN 70 TO NEXT
      GO TO 60
C
C                                PREPARE FOR PHASE 4.
C
   50 I = J
      ASSIGN 80 TO NEXT
      SUM = (SUM/SX(I))/SX(I)
   60 XMAX = ABS(SX(I))
      GO TO 90
C
C                   PHASE 2.  SUM IS SMALL.
C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
C
   70 IF (ABS(SX(I)).GT.CUTLO) GO TO 100
C
C                     COMMON CODE FOR PHASES 2 AND 4.
C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
C
   80 IF (ABS(SX(I)).LE.XMAX) GO TO 90
      SUM = ONE + SUM*(XMAX/SX(I))**2
      XMAX = ABS(SX(I))
      GO TO 130
C
   90 SUM = SUM + (SX(I)/XMAX)**2
      GO TO 130
C
C
C                  PREPARE FOR PHASE 3.
C
  100 SUM = (SUM*XMAX)*XMAX
C
C
C     FOR REAL OR D.P. SET HITEST = CUTHI/N
C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
C
  110 HITEST = CUTHI/FLOAT(N)
C
C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
C
      DO 120 J=I,NN,INCX
        IF (ABS(SX(J)).GE.HITEST) GO TO 50
        SUM = SUM + SX(J)**2
  120 CONTINUE
      SNRM2 = SQRT(SUM)
      GO TO 140
C
  130 CONTINUE
      I = I + INCX
      IF (I.LE.NN) GO TO 20
C
C              END OF MAIN LOOP.
C
C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
C
      SNRM2 = XMAX*SQRT(SUM)
  140 CONTINUE
      RETURN
      END
      SUBROUTINE SROT(N, SX, INCX, SY, INCY, C, S)
C
C     APPLIES A PLANE ROTATION.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SX(1), SY(1), STEMP, C, S
      INTEGER I, INCX, INCY, IX, IY, N
C
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
C
C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
C         TO 1
C
      IX = 1
      IY = 1
      IF (INCX.LT.0) IX = (-N+1)*INCX + 1
      IF (INCY.LT.0) IY = (-N+1)*INCY + 1
      DO 10 I=1,N
        STEMP = C*SX(IX) + S*SY(IY)
        SY(IY) = C*SY(IY) - S*SX(IX)
        SX(IX) = STEMP
        IX = IX + INCX
        IY = IY + INCY
   10 CONTINUE
      RETURN
C
C       CODE FOR BOTH INCREMENTS EQUAL TO 1
C
   20 DO 30 I=1,N
        STEMP = C*SX(I) + S*SY(I)
        SY(I) = C*SY(I) - S*SX(I)
        SX(I) = STEMP
   30 CONTINUE
      RETURN
      END
      SUBROUTINE SROTG(SA, SB, C, S)
C
C     CONSTRUCT GIVENS PLANE ROTATION.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SA, SB, C, S, ROE, SCALE, R, Z
C
      ROE = SB
      IF (ABS(SA).GT.ABS(SB)) ROE = SA
      SCALE = ABS(SA) + ABS(SB)
      IF (SCALE.NE.0.0) GO TO 10
      C = 1.0
      S = 0.0
      R = 0.0
      GO TO 20
   10 R = SCALE*SQRT((SA/SCALE)**2+(SB/SCALE)**2)
      R = SIGN(1.0,ROE)*R
      C = SA/R
      S = SB/R
   20 Z = 1.0
      IF (ABS(SA).GT.ABS(SB)) Z = S
      IF (ABS(SB).GE.ABS(SA) .AND. C.NE.0.0) Z = 1.0/C
      SA = R
      SB = Z
      RETURN
      END
      SUBROUTINE SSCAL(N, SA, SX, INCX)
C
C     SCALES A VECTOR BY A CONSTANT.
C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1.
C     JACK DONGARRA, LINPACK, 3/11/78.
C
      REAL SA, SX(1)
      INTEGER I, INCX, M, MP1, N, NINCX
C
      IF (N.LE.0) RETURN
      IF (INCX.EQ.1) GO TO 20
C
C        CODE FOR INCREMENT NOT EQUAL TO 1
C
      NINCX = N*INCX
      DO 10 I=1,NINCX,INCX
        SX(I) = SA*SX(I)
   10 CONTINUE
      RETURN
C
C        CODE FOR INCREMENT EQUAL TO 1
C
C
C        CLEAN-UP LOOP
C
   20 M = MOD(N,5)
      IF (M.EQ.0) GO TO 40
      DO 30 I=1,M
        SX(I) = SA*SX(I)
   30 CONTINUE
      IF (N.LT.5) RETURN
   40 MP1 = M + 1
      DO 50 I=MP1,N,5
        SX(I) = SA*SX(I)
        SX(I+1) = SA*SX(I+1)
        SX(I+2) = SA*SX(I+2)
        SX(I+3) = SA*SX(I+3)
        SX(I+4) = SA*SX(I+4)
   50 CONTINUE
      RETURN
      END
      SUBROUTINE SVD(NM, M, N, A, W, MATU, U, MATV, V, IERR, RV1)
C
      INTEGER I, J, K, L, M, N, II, I1, KK, K1, LL, L1, MN, NM, ITS,
     * IERR
      REAL A(NM,N), W(N), U(NM,N), V(NM,N), RV1(N)
      REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, MACHEP
      REAL SQRT, AMAX1, ABS, SIGN
      LOGICAL MATU, MATV
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD,
C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
C
C     THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION
C          T
C     A=USV  OF A REAL M BY N RECTANGULAR MATRIX.  HOUSEHOLDER
C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
C          AS LARGE AS THE MAXIMUM OF M AND N,
C
C        M IS THE NUMBER OF ROWS OF A (AND U),
C
C        N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V,
C
C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED,
C
C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE,
C
C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
C
C     ON OUTPUT-
C
C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V),
C
C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
C          FOR INDICES IERR+1,IERR+2,...,N,
C
C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  OTHERWISE
C          U IS USED AS A TEMPORARY ARRAY.  U MAY COINCIDE WITH A.
C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
C
C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
C          MATV HAS BEEN SET TO .TRUE.  OTHERWISE V IS NOT REFERENCED.
C          V MAY ALSO COINCIDE WITH A IF U IS NOT NEEDED.  IF AN ERROR
C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
C          CORRECT SINGULAR VALUES SHOULD BE CORRECT,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS,
C
C        RV1 IS A TEMPORARY STORAGE ARRAY.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
      MACHEP = 2.**(-26)
C
      IERR = 0
C
      DO 20 I=1,M
C
        DO 10 J=1,N
          U(I,J) = A(I,J)
   10   CONTINUE
   20 CONTINUE
C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
      G = 0.0
      SCALE = 0.0
      X = 0.0
C
      DO 200 I=1,N
        L = I + 1
        RV1(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
        IF (I.GT.M) GO TO 100
C
        DO 30 K=I,M
          SCALE = SCALE + ABS(U(K,I))
   30   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 100
C
        DO 40 K=I,M
          U(K,I) = U(K,I)/SCALE
          S = S + U(K,I)**2
   40   CONTINUE
C
        F = U(I,I)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        U(I,I) = F - G
        IF (I.EQ.N) GO TO 80
C
        DO 70 J=L,N
          S = 0.0
C
          DO 50 K=I,M
            S = S + U(K,I)*U(K,J)
   50     CONTINUE
C
          F = S/H
C
          DO 60 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
   60     CONTINUE
   70   CONTINUE
C
   80   DO 90 K=I,M
          U(K,I) = SCALE*U(K,I)
   90   CONTINUE
C
  100   W(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
        IF (I.GT.M .OR. I.EQ.N) GO TO 190
C
        DO 110 K=L,N
          SCALE = SCALE + ABS(U(I,K))
  110   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 190
C
        DO 120 K=L,N
          U(I,K) = U(I,K)/SCALE
          S = S + U(I,K)**2
  120   CONTINUE
C
        F = U(I,L)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        U(I,L) = F - G
C
        DO 130 K=L,N
          RV1(K) = U(I,K)/H
  130   CONTINUE
C
        IF (I.EQ.M) GO TO 170
C
        DO 160 J=L,M
          S = 0.0
C
          DO 140 K=L,N
            S = S + U(J,K)*U(I,K)
  140     CONTINUE
C
          DO 150 K=L,N
            U(J,K) = U(J,K) + S*RV1(K)
  150     CONTINUE
  160   CONTINUE
C
  170   DO 180 K=L,N
          U(I,K) = SCALE*U(I,K)
  180   CONTINUE
C
  190   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  200 CONTINUE
C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS **********
      IF (.NOT.MATV) GO TO 290
C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
      DO 280 II=1,N
        I = N + 1 - II
        IF (I.EQ.N) GO TO 270
        IF (G.EQ.0.0) GO TO 250
C
        DO 210 J=L,N
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
          V(J,I) = (U(I,J)/U(I,L))/G
  210   CONTINUE
C
        DO 240 J=L,N
          S = 0.0
C
          DO 220 K=L,N
            S = S + U(I,K)*V(K,J)
  220     CONTINUE
C
          DO 230 K=L,N
            V(K,J) = V(K,J) + S*V(K,I)
  230     CONTINUE
  240   CONTINUE
C
  250   DO 260 J=L,N
          V(I,J) = 0.0
          V(J,I) = 0.0
  260   CONTINUE
C
  270   V(I,I) = 1.0
        G = RV1(I)
        L = I
  280 CONTINUE
C     ********** ACCUMULATION OF LEFT-HAND TRANSFORMATIONS **********
  290 IF (.NOT.MATU) GO TO 410
C     **********FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- **********
      MN = N
      IF (M.LT.N) MN = M
C
      DO 400 II=1,MN
        I = MN + 1 - II
        L = I + 1
        G = W(I)
        IF (I.EQ.N) GO TO 310
C
        DO 300 J=L,N
          U(I,J) = 0.0
  300   CONTINUE
C
  310   IF (G.EQ.0.0) GO TO 370
        IF (I.EQ.MN) GO TO 350
C
        DO 340 J=L,N
          S = 0.0
C
          DO 320 K=L,M
            S = S + U(K,I)*U(K,J)
  320     CONTINUE
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
          F = (S/U(I,I))/G
C
          DO 330 K=I,M
            U(K,J) = U(K,J) + F*U(K,I)
  330     CONTINUE
  340   CONTINUE
C
  350   DO 360 J=I,M
          U(J,I) = U(J,I)/G
  360   CONTINUE
C
        GO TO 390
C
  370   DO 380 J=I,M
          U(J,I) = 0.0
  380   CONTINUE
C
  390   U(I,I) = U(I,I) + 1.0
  400 CONTINUE
C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  410 EPS = MACHEP*X
C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
      DO 550 KK=1,N
        K1 = N - KK
        K = K1 + 1
        ITS = 0
C     ********** TEST FOR SPLITTING.
C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  420   DO 430 LL=1,K
          L1 = K - LL
          L = L1 + 1
          IF (ABS(RV1(L)).LE.EPS) GO TO 470
C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
          IF (ABS(W(L1)).LE.EPS) GO TO 440
  430   CONTINUE
C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  440   C = 0.0
        S = 1.0
C
        DO 460 I=L,K
          F = S*RV1(I)
          RV1(I) = C*RV1(I)
          IF (ABS(F).LE.EPS) GO TO 470
          G = W(I)
          H = SQRT(F*F+G*G)
          W(I) = H
          C = G/H
          S = -F/H
          IF (.NOT.MATU) GO TO 460
C
          DO 450 J=1,M
            Y = U(J,L1)
            Z = U(J,I)
            U(J,L1) = Y*C + Z*S
            U(J,I) = -Y*S + Z*C
  450     CONTINUE
C
  460   CONTINUE
C     ********** TEST FOR CONVERGENCE **********
  470   Z = W(K)
        IF (L.EQ.K) GO TO 530
C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
        IF (ITS.EQ.30) GO TO 560
        ITS = ITS + 1
        X = W(L)
        Y = W(K1)
        G = RV1(K1)
        H = RV1(K)
        F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
        G = SQRT(F*F+1.0)
        F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C     ********** NEXT QR TRANSFORMATION **********
        C = 1.0
        S = 1.0
C
        DO 520 I1=L,K1
          I = I1 + 1
          G = RV1(I)
          Y = W(I)
          H = S*G
          G = C*G
          Z = SQRT(F*F+H*H)
          RV1(I1) = Z
          C = F/Z
          S = H/Z
          F = X*C + G*S
          G = -X*S + G*C
          H = Y*S
          Y = Y*C
          IF (.NOT.MATV) GO TO 490
C
          DO 480 J=1,N
            X = V(J,I1)
            Z = V(J,I)
            V(J,I1) = X*C + Z*S
            V(J,I) = -X*S + Z*C
  480     CONTINUE
C
  490     Z = SQRT(F*F+H*H)
          W(I1) = Z
C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
          IF (Z.EQ.0.0) GO TO 500
          C = F/Z
          S = H/Z
  500     F = C*G + S*Y
          X = -S*G + C*Y
          IF (.NOT.MATU) GO TO 520
C
          DO 510 J=1,M
            Y = U(J,I1)
            Z = U(J,I)
            U(J,I1) = Y*C + Z*S
            U(J,I) = -Y*S + Z*C
  510     CONTINUE
C
  520   CONTINUE
C
        RV1(L) = 0.0
        RV1(K) = F
        W(K) = X
        GO TO 420
C     ********** CONVERGENCE **********
  530   IF (Z.GE.0.0) GO TO 550
C     ********** W(K) IS MADE NON-NEGATIVE **********
        W(K) = -Z
        IF (.NOT.MATV) GO TO 550
C
        DO 540 J=1,N
          V(J,K) = -V(J,K)
  540   CONTINUE
C
  550 CONTINUE
C
      GO TO 570
C     ********** SET ERROR -- NO CONVERGENCE TO A
C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  560 IERR = K
  570 RETURN
C     ********** LAST CARD OF SVD **********
      END
C
C     ------------------------------------------------------------------
C
      SUBROUTINE MINFIT(NM, M, N, A, W, IP, B, IERR, RV1)
C
      INTEGER I, J, K, L, M, N, II, IP, I1, KK, K1, LL, L1, M1, NM,
     * ITS, IERR
      REAL A(NM,N), W(N), B(NM,IP), RV1(N)
      REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, MACHEP
      REAL SQRT, AMAX1, ABS, SIGN
C
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT,
C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
C
C     THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR
C                                                        T
C     SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV  OF A REAL
C                                         T
C     M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U.  HOUSEHOLDER
C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
C
C     ON INPUT-
C
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
C          AS LARGE AS THE MAXIMUM OF M AND N,
C
C        M IS THE NUMBER OF ROWS OF A AND B,
C
C        N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V,
C
C        A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM,
C
C        IP IS THE NUMBER OF COLUMNS OF B.  IP CAN BE ZERO,
C
C        B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM
C          IF IP IS NOT ZERO.  OTHERWISE B IS NOT REFERENCED.
C
C     ON OUTPUT-
C
C        A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE
C          DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS.  IF AN
C          ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO
C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
C
C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
C          FOR INDICES IERR+1,IERR+2,...,N,
C
C                                   T
C        B HAS BEEN OVERWRITTEN BY U B.  IF AN ERROR EXIT IS MADE,
C                       T
C          THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT
C          SINGULAR VALUES SHOULD BE CORRECT,
C
C        IERR IS SET TO
C          ZERO       FOR NORMAL RETURN,
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
C                     DETERMINED AFTER 30 ITERATIONS,
C
C        RV1 IS A TEMPORARY STORAGE ARRAY.
C
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
C
C     ------------------------------------------------------------------
C
C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
C
C                **********
      MACHEP = 2.**(-26)
C
      IERR = 0
C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
      G = 0.0
      SCALE = 0.0
      X = 0.0
C
      DO 220 I=1,N
        L = I + 1
        RV1(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
        IF (I.GT.M) GO TO 120
C
        DO 10 K=I,M
          SCALE = SCALE + ABS(A(K,I))
   10   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 120
C
        DO 20 K=I,M
          A(K,I) = A(K,I)/SCALE
          S = S + A(K,I)**2
   20   CONTINUE
C
        F = A(I,I)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        A(I,I) = F - G
        IF (I.EQ.N) GO TO 60
C
        DO 50 J=L,N
          S = 0.0
C
          DO 30 K=I,M
            S = S + A(K,I)*A(K,J)
   30     CONTINUE
C
          F = S/H
C
          DO 40 K=I,M
            A(K,J) = A(K,J) + F*A(K,I)
   40     CONTINUE
   50   CONTINUE
C
   60   IF (IP.EQ.0) GO TO 100
C
        DO 90 J=1,IP
          S = 0.0
C
          DO 70 K=I,M
            S = S + A(K,I)*B(K,J)
   70     CONTINUE
C
          F = S/H
C
          DO 80 K=I,M
            B(K,J) = B(K,J) + F*A(K,I)
   80     CONTINUE
   90   CONTINUE
C
  100   DO 110 K=I,M
          A(K,I) = SCALE*A(K,I)
  110   CONTINUE
C
  120   W(I) = SCALE*G
        G = 0.0
        S = 0.0
        SCALE = 0.0
        IF (I.GT.M .OR. I.EQ.N) GO TO 210
C
        DO 130 K=L,N
          SCALE = SCALE + ABS(A(I,K))
  130   CONTINUE
C
        IF (SCALE.EQ.0.0) GO TO 210
C
        DO 140 K=L,N
          A(I,K) = A(I,K)/SCALE
          S = S + A(I,K)**2
  140   CONTINUE
C
        F = A(I,L)
        G = -SIGN(SQRT(S),F)
        H = F*G - S
        A(I,L) = F - G
C
        DO 150 K=L,N
          RV1(K) = A(I,K)/H
  150   CONTINUE
C
        IF (I.EQ.M) GO TO 190
C
        DO 180 J=L,M
          S = 0.0
C
          DO 160 K=L,N
            S = S + A(J,K)*A(I,K)
  160     CONTINUE
C
          DO 170 K=L,N
            A(J,K) = A(J,K) + S*RV1(K)
  170     CONTINUE
  180   CONTINUE
C
  190   DO 200 K=L,N
          A(I,K) = SCALE*A(I,K)
  200   CONTINUE
C
  210   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  220 CONTINUE
C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS.
C                FOR I=N STEP -1 UNTIL 1 DO -- **********
      DO 300 II=1,N
        I = N + 1 - II
        IF (I.EQ.N) GO TO 290
        IF (G.EQ.0.0) GO TO 270
C
        DO 230 J=L,N
C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
          A(J,I) = (A(I,J)/A(I,L))/G
  230   CONTINUE
C
        DO 260 J=L,N
          S = 0.0
C
          DO 240 K=L,N
            S = S + A(I,K)*A(K,J)
  240     CONTINUE
C
          DO 250 K=L,N
            A(K,J) = A(K,J) + S*A(K,I)
  250     CONTINUE
  260   CONTINUE
C
  270   DO 280 J=L,N
          A(I,J) = 0.0
          A(J,I) = 0.0
  280   CONTINUE
C
  290   A(I,I) = 1.0
        G = RV1(I)
        L = I
  300 CONTINUE
C
      IF (M.GE.N .OR. IP.EQ.0) GO TO 330
      M1 = M + 1
C
      DO 320 I=M1,N
C
        DO 310 J=1,IP
          B(I,J) = 0.0
  310   CONTINUE
  320 CONTINUE
C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  330 EPS = MACHEP*X
C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
      DO 460 KK=1,N
        K1 = N - KK
        K = K1 + 1
        ITS = 0
C     ********** TEST FOR SPLITTING.
C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  340   DO 350 LL=1,K
          L1 = K - LL
          L = L1 + 1
          IF (ABS(RV1(L)).LE.EPS) GO TO 390
C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
C                THROUGH THE BOTTOM OF THE LOOP **********
          IF (ABS(W(L1)).LE.EPS) GO TO 360
  350   CONTINUE
C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  360   C = 0.0
        S = 1.0
C
        DO 380 I=L,K
          F = S*RV1(I)
          RV1(I) = C*RV1(I)
          IF (ABS(F).LE.EPS) GO TO 390
          G = W(I)
          H = SQRT(F*F+G*G)
          W(I) = H
          C = G/H
          S = -F/H
          IF (IP.EQ.0) GO TO 380
C
          DO 370 J=1,IP
            Y = B(L1,J)
            Z = B(I,J)
            B(L1,J) = Y*C + Z*S
            B(I,J) = -Y*S + Z*C
  370     CONTINUE
C
  380   CONTINUE
C     ********** TEST FOR CONVERGENCE **********
  390   Z = W(K)
        IF (L.EQ.K) GO TO 440
C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
        IF (ITS.EQ.30) GO TO 470
        ITS = ITS + 1
        X = W(L)
        Y = W(K1)
        G = RV1(K1)
        H = RV1(K)
        F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
        G = SQRT(F*F+1.0)
        F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
C     ********** NEXT QR TRANSFORMATION **********
        C = 1.0
        S = 1.0
C
        DO 430 I1=L,K1
          I = I1 + 1
          G = RV1(I)
          Y = W(I)
          H = S*G
          G = C*G
          Z = SQRT(F*F+H*H)
          RV1(I1) = Z
          C = F/Z
          S = H/Z
          F = X*C + G*S
          G = -X*S + G*C
          H = Y*S
          Y = Y*C
C
          DO 400 J=1,N
            X = A(J,I1)
            Z = A(J,I)
            A(J,I1) = X*C + Z*S
            A(J,I) = -X*S + Z*C
  400     CONTINUE
C
          Z = SQRT(F*F+H*H)
          W(I1) = Z
C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
          IF (Z.EQ.0.0) GO TO 410
          C = F/Z
          S = H/Z
  410     F = C*G + S*Y
          X = -S*G + C*Y
          IF (IP.EQ.0) GO TO 430
C
          DO 420 J=1,IP
            Y = B(I1,J)
            Z = B(I,J)
            B(I1,J) = Y*C + Z*S
            B(I,J) = -Y*S + Z*C
  420     CONTINUE
C
  430   CONTINUE
C
        RV1(L) = 0.0
        RV1(K) = F
        W(K) = X
        GO TO 340
C     ********** CONVERGENCE **********
  440   IF (Z.GE.0.0) GO TO 460
C     ********** W(K) IS MADE NON-NEGATIVE **********
        W(K) = -Z
C
        DO 450 J=1,N
          A(J,K) = -A(J,K)
  450   CONTINUE
C
  460 CONTINUE
C
      GO TO 480
C     ********** SET ERROR -- NO CONVERGENCE TO A
C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  470 IERR = K
  480 RETURN
C     ********** LAST CARD OF MINFIT **********
      END
 
