home *** CD-ROM | disk | FTP | other *** search
/ High Voltage Shareware / high1.zip / high1 / DIR9 / ACMALG01.ZIP / ACM581.FOR < prev    next >
Text File  |  1991-02-18  |  86KB  |  2,950 lines
  1. C     ALGORITHM 581, COLLECTED ALGORITHMS FROM ACM. THIS WORK
  2. C     PUBLISHED IN TRANS. MATH. SOFTWARE, 8(1), PP. 84-88.
  3.  
  4. C***********************************************************************
  5. C
  6. C       THIS FILE CONTAINS PROGRAMS AND SUBROUTINES THAT ARE RELATED
  7. C       TO THE PAPER 'AN IMPROVED ALGORITHM FOR COMPUTING THE SINGULAR
  8. C       VALUE DECOMPOSITION' BY T.F. CHAN, WHICH WAS SUBMITTED FOR
  9. C       PUBLICATIONS IN ACM TOMS.
  10. C
  11. C       THIS FILE CONTAINS ONE (1) DRIVER, SIX  (6) MAIN SUBROUTINES,
  12. C       SOME BLAS ROUTINES, A ROUTINE FOR COMPUTING THE MACHINE RELATIVE
  13. C       PRECISION AND THE OUTPUT PRODUCED BY THE DRIVER.
  14. C       THEY ARE LISTED BELOW IN THE ORDER THAT THEY APPEAR IN THE
  15. C       FILE:
  16. C
  17. C           1) DRIVER FOR THE ROUTINE HYBSVD TO COMPUTE THE SVD
  18. C              OF A TEST MATRIX AND ITS TRANSPOSE.
  19. C
  20. C           2) ROUTINE HYBSVD ... HYBRID ALGORITHM ROUTINE.
  21. C
  22. C           3) ROUTINE MGNSVD ... HYBRID ALGORITHM, ASSUMES M .GE. N.
  23. C
  24. C           4) ROUTINE GRSVD  ... GOLUB-REINSCH ALGORITHM ROUTINE
  25. C
  26. C           5) BLAS ROUTINE   ... SSWAP, FOR SWAPPING VECTORS
  27. C
  28. C           6) ROUTINE SRELPR ... COMPUTES MACHINE PRECISION
  29. C
  30. C           7) ROUTINE SSVDC  ... LINPACK SVD ROUTINE
  31. C
  32. C           8) BLAS ROUTINES  ... SAXPY, SDOT, SNRM2, SROT, SROTG,
  33. C                                 SSCAL.
  34. C
  35. C           9) ROUTINE SVD    ... EISPACK SVD ROUTINE
  36. C
  37. C          10) ROUTINE MINFIT ... EISPACK MINFIT ROUTINE
  38. C
  39. C          11) OUTPUT PRODUCED BY THE DRIVER IN 1).
  40. C
  41. C       THE ROUTINES 2) TO 6) CONSTITUTE THE PACKAGE THAT IMPLEMENTS
  42. C       THE NEW HYBRID ALGORITHM AND CAN BE USED BY THEMSELVES.
  43. C
  44. C       PLEASE ADDRESS COMMENTS AND SUGGESTIONS TO:
  45. C
  46. C           TONY CHAN
  47. C           COMPUTER SCIENCE DEPT., YALE UNIV.,
  48. C           BOX 2158, YALE STATION,
  49. C           NEW HAVEN, CT 06520.
  50. C
  51. C***********************************************************************
  52. C
  53. C       DRIVER FOR HYBSVD, A ROUTINE FOR COMPUTING THE SVD OF A MATRIX.
  54. C       RESULTS ARE COMPARED TO THOSE PRODUCED BY SSVDC FROM LINPACK,
  55. C       SVD AND MINFIT FROM EISPACK.
  56. C       PREPARED FOR ACM TRANS. ON MATH. SOFTWARE.
  57. C
  58. C       MATRIX TESTED:  A(I,J) = 1/(I+J),  1ST RHS(I) = 1
  59. C               T                          2ND RHS(I) = I
  60. C       (A AND A )
  61. C
  62. C       PROGRAMMED BY
  63. C           TONY CHAN
  64. C           COMPUTER SCIENCE DEPT., YALE UNIV.,
  65. C           BOX 2158, YALE STATION,
  66. C           NEW HAVEN, CT 06520.
  67. C
  68. C       LAST REVISION.. JANUARY, 1982.
  69. C
  70.       REAL AE(15,15), AH(10,10), AL(20,20), WORK(10), E(10)
  71.       REAL UH(11,10), VH(12,20), Z(13,20), BH(14,5), WH(10)
  72.       REAL UE(15,10), VE(15,10), BE(15,5), WE(10)
  73.       REAL UL(18,10), VL(19,10), WL(10)
  74.       LOGICAL MATU, MATV
  75.       INTEGER NUH, NVH, NZ, NBH
  76.       INTEGER NUE, NVE, NBE
  77.       INTEGER NUL, NVL
  78.       INTEGER NAE, NAH, NAL, M, N, IRHS, IERR, NOUT, JOB
  79.       DATA NAH /10/, NUH /11/, NVH /12/, NZ /13/, NBH /14/
  80.       DATA NUE /15/, NVE /15/, NBE /15/, NAE /15/
  81.       DATA NUL /18/, NVL /19/, NAL /20/
  82. C
  83. C     THIS IS OUTPUT UNIT NUMBER.
  84.       DATA NOUT /6/
  85. C
  86. C     FIRST CASE... A AND RHS B.
  87. C
  88.       M = 8
  89.       N = 4
  90.       IRHS = 2
  91. C
  92.       DO 50 IU=1,2
  93.         DO 40 IV=1,2
  94.           IF (IU.EQ.2) MATU = .TRUE.
  95.           IF (IU.EQ.1) MATU = .FALSE.
  96.           IF (IV.EQ.2) MATV = .TRUE.
  97.           IF (IV.EQ.1) MATV = .FALSE.
  98. C
  99. C       JOB PARAMETER NEEDED BY SSVDC.
  100.           JOB = (IU-1)*20 + IV - 1
  101. C
  102. C       SET UP MATRIX
  103. C
  104.           DO 20 I=1,M
  105.             BH(I,1) = 1.0
  106.             BH(I,2) = FLOAT(I)
  107.             BE(I,1) = 1.0
  108.             BE(I,2) = FLOAT(I)
  109.             DO 10 J=1,N
  110.               AH(I,J) = 1./FLOAT(I+J)
  111.               AE(I,J) = AH(I,J)
  112.               AL(I,J) = AH(I,J)
  113.    10       CONTINUE
  114.    20     CONTINUE
  115. C
  116.           WRITE (NOUT,99999)
  117. 99999     FORMAT (/4H A ./)
  118.           DO 30 I=1,M
  119.             WRITE (NOUT,99998) (AH(I,J),J=1,N)
  120. 99998       FORMAT (5E15.7)
  121.    30     CONTINUE
  122. C
  123.           WRITE (NOUT,99997)
  124. 99997     FORMAT (//45H ******************* HYBSVD *****************//)
  125.           CALL HYBSVD(NAH, NUH, NVH, NZ, NBH, M, N, AH, WH, MATU, UH,
  126.      *     MATV, VH, Z, BH, IRHS, IERR, WORK)
  127.           CALL PRINTS(WH, UH, VH, IERR, BH, NUH, NVH, NBH, MATU, MATV,
  128.      *     M, N, NOUT, IRHS)
  129. C
  130.           WRITE (NOUT,99996)
  131. 99996     FORMAT (//45H ******************* EISPAK *****************//)
  132.           CALL SVD(NUE, M, N, AE, WE, MATU, UE, MATV, VE, IERR, WORK)
  133.           CALL MINFIT(NUE, M, N, AE, WE, IRHS, BE, IERR, WORK)
  134.           CALL PRINTS(WE, UE, VE, IERR, BE, NUE, NVE, NBE, MATU, MATV,
  135.      *     M, N, NOUT, IRHS)
  136. C
  137.           WRITE (NOUT,99995)
  138. 99995     FORMAT (//45H ******************* SSVDC  *****************//)
  139.           CALL SSVDC(AL, NAL, M, N, WL, E, UL, NUL, VL, NVL, WORK, JOB,
  140.      *     IERR)
  141.           CALL PRINTS(WL, UL, VL, IERR, E, NUL, NVL, NBH, MATU, MATV,
  142.      *     M, N, NOUT, IRHS)
  143. C
  144.           CALL DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
  145.      *     NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU,
  146.      *     MATV)
  147.    40   CONTINUE
  148.    50 CONTINUE
  149. C                      T
  150. C     SECOND CASE...  A
  151. C
  152.       M = 4
  153.       N = 8
  154.       IRHS = 2
  155. C
  156.       DO 100 IU=1,2
  157.         DO 90 IV=1,2
  158.           IF (IU.EQ.2) MATU = .TRUE.
  159.           IF (IU.EQ.1) MATU = .FALSE.
  160.           IF (IV.EQ.2) MATV = .TRUE.
  161.           IF (IV.EQ.1) MATV = .FALSE.
  162. C
  163. C       JOB PARAMETER NEEDED BY SSVDC.
  164.           JOB = (IU-1)*20 + IV - 1
  165. C
  166. C       SET UP MATRIX
  167. C
  168.           DO 70 I=1,M
  169.             BH(I,1) = 1.
  170.             BH(I,2) = FLOAT(I)
  171.             BE(I,1) = 1.
  172.             BE(I,2) = FLOAT(I)
  173.             DO 60 J=1,N
  174.               AH(I,J) = 1./FLOAT(I+J)
  175.               AE(I,J) = AH(I,J)
  176.               AL(I,J) = AH(I,J)
  177.    60       CONTINUE
  178.    70     CONTINUE
  179. C
  180.           WRITE (NOUT,99999)
  181.           DO 80 I=1,M
  182.             WRITE (NOUT,99998) (AH(I,J),J=1,N)
  183.    80     CONTINUE
  184. C
  185.           WRITE (NOUT,99997)
  186.           CALL HYBSVD(NAH, NUH, NVH, NZ, NBH, M, N, AH, WH, MATU, UH,
  187.      *     MATV, VH, Z, BH, IRHS, IERR, WORK)
  188.           CALL PRINTS(WH, UH, VH, IERR, BH, NUH, NVH, NBH, MATU, MATV,
  189.      *     M, N, NOUT, IRHS)
  190. C
  191.           WRITE (NOUT,99996)
  192.           CALL SVD(NUE, M, N, AE, WE, MATU, UE, MATV, VE, IERR, WORK)
  193.           CALL MINFIT(NUE, M, N, AE, WE, IRHS, BE, IERR, WORK)
  194.           CALL PRINTS(WE, UE, VE, IERR, BE, NUE, NVE, NBE, MATU, MATV,
  195.      *     M, N, NOUT, IRHS)
  196. C
  197.           WRITE (NOUT,99995)
  198.           CALL SSVDC(AL, NAL, M, N, WL, E, UL, NUL, VL, NVL, WORK, JOB,
  199.      *     IERR)
  200.           CALL PRINTS(WL, UL, VL, IERR, E, NUL, NVL, NBH, MATU, MATV,
  201.      *     M, N, NOUT, IRHS)
  202.           CALL DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
  203.      *     NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU,
  204.      *     MATV)
  205.    90   CONTINUE
  206.   100 CONTINUE
  207.       STOP
  208.       END
  209.       SUBROUTINE PRINTS(W, U, V, IERR, B, NAU, NV, NB, MATU, MATV, M,
  210.      * N, NOUT, IRHS)
  211. C
  212. C       PRINTS OUT THE SINGULAR VALUES AND THE MATRICES U, V, B.
  213. C
  214.       INTEGER NAU, NV, IERR, M, N, NOUT, IRHS
  215.       REAL U(NAU,1), V(NV,1), W(1), B(NB,1)
  216.       LOGICAL MATU, MATV
  217.       MINMN = MIN0(M,N)
  218. C
  219.       WRITE (NOUT,99999) MATU, MATV
  220. 99999 FORMAT (7H MATU =, L3, 7H MATV =, L3/)
  221.       WRITE (NOUT,99998) IERR
  222. 99998 FORMAT (7H IERR =, I5/)
  223. C
  224.       WRITE (NOUT,99997)
  225. 99997 FORMAT (17H SINGULAR VALUES./)
  226.       DO 10 I=1,MINMN
  227.         WRITE (NOUT,99995) W(I)
  228.    10 CONTINUE
  229. C
  230.       IF (.NOT.MATU) GO TO 30
  231.       WRITE (NOUT,99996)
  232. 99996 FORMAT (/4H U ./)
  233. 99995 FORMAT (E15.8)
  234.       DO 20 I=1,M
  235.         WRITE (NOUT,99993) (U(I,J),J=1,MINMN)
  236.    20 CONTINUE
  237. C
  238.    30 IF (.NOT.MATV) GO TO 50
  239.       WRITE (NOUT,99994)
  240. 99994 FORMAT (/4H V ./)
  241. 99993 FORMAT (5E15.7)
  242.       DO 40 I=1,N
  243.         WRITE (NOUT,99993) (V(I,J),J=1,MINMN)
  244.    40 CONTINUE
  245. C
  246.    50 IF (M.LT.N .OR. IRHS.LT.1) GO TO 70
  247.       WRITE (NOUT,99992)
  248. 99992 FORMAT (/4H B ./)
  249.       DO 60 I=1,MINMN
  250.         WRITE (NOUT,99993) (B(I,J),J=1,IRHS)
  251.    60 CONTINUE
  252. C
  253.    70 CONTINUE
  254.       RETURN
  255.       END
  256.       SUBROUTINE DIFF(WE, WH, WL, UE, UH, UL, VE, VH, VL, BH, BE, NUE,
  257.      * NUH, NUL, NVE, NVH, NVL, NBH, NBE, M, N, NOUT, IRHS, MATU, MATV)
  258. C
  259. C     COMPUTES DIFFERENCE BETWEEN W, U, V AND B AS COMPUTED BY
  260. C     EISPACK SVD AND MINFIT, HYBSVD, AND LINPACK SSVDC.
  261. C
  262.       REAL WE(1), WH(1), WL(1)
  263.       REAL UE(NUE,1), UH(NUH,1), UL(NUL,1)
  264.       REAL VE(NVE,1), VH(NVH,1), VL(NVL,1)
  265.       REAL BE(NBE,1), BH(NBH,1)
  266.       INTEGER NOUT, IRHS, M, N
  267.       LOGICAL MATU, MATV
  268. C
  269. C     SORT SINGULAR VALUES AND EXCHANGE COLUMNS OF U AND V FOR EISPACK.
  270. C     SELECTION SORT MINIMIZES SWAPPING OF U AND V.
  271. C
  272.       IF (N.LE.1) GO TO 40
  273.       NM1 = N - 1
  274.       DO 30 I=1,NM1
  275. C...    FIND INDEX OF MAXIMUM SINGULAR VALUE
  276.         ID = I
  277.         IP1 = I + 1
  278.         DO 10 J=IP1,N
  279.           IF (WE(J).GT.WE(ID)) ID = J
  280.    10   CONTINUE
  281.         IF (ID.EQ.I) GO TO 30
  282. C...    SWAP SINGULAR VALUES AND VECTORS
  283.         T = WE(I)
  284.         WE(I) = WE(ID)
  285.         WE(ID) = T
  286.         IF (MATV) CALL SSWAP(N, VE(1,I), 1, VE(1,ID), 1)
  287.         IF (MATU) CALL SSWAP(M, UE(1,I), 1, UE(1,ID), 1)
  288.         IF (IRHS.LT.1) GO TO 30
  289.         DO 20 KRHS=1,IRHS
  290.           T = BE(I,KRHS)
  291.           BE(I,KRHS) = BE(ID,KRHS)
  292.           BE(ID,KRHS) = T
  293.    20   CONTINUE
  294.    30 CONTINUE
  295. C
  296.    40 CONTINUE
  297.       DWHE = 0.
  298.       DWHL = 0.
  299.       DWEL = 0.
  300.       DUHE = 0.
  301.       DUHL = 0.
  302.       DUEL = 0.
  303.       DVHE = 0.
  304.       DVHL = 0.
  305.       DVEL = 0.
  306.       DBHE = 0.
  307.       MINMN = MIN0(M,N)
  308. C
  309.       DO 50 I=1,MINMN
  310.         DWHE = DWHE + (WH(I)-WE(I))**2
  311.         DWHL = DWHL + (WH(I)-WL(I))**2
  312.         DWEL = DWEL + (WE(I)-WL(I))**2
  313.    50 CONTINUE
  314.       DWHE = SQRT(DWHE)
  315.       DWHL = SQRT(DWHL)
  316.       DWEL = SQRT(DWEL)
  317.       WRITE (NOUT,99999) DWHE, DWHL, DWEL
  318. 99999 FORMAT (/35H  2-NORM ( HYBSVD W - EISPAK W ) = , 1PE15.7/6H  2-NO,
  319.      * 29HRM ( HYBSVD W - LINPAK W ) = , 1PE15.7/19H  2-NORM ( EISPAK W,
  320.      * 16H - LINPAK W ) = , 1PE15.7/)
  321. C
  322. C     DIFFERENCE OF ABSOLUTE VALUES BECAUSE OF POSSIBLE SIGN DIFFERENCE.
  323. C
  324.       IF (.NOT.MATU) GO TO 80
  325.       DO 70 J=1,MINMN
  326.         DO 60 I=1,M
  327.           DUHE = DUHE + (ABS(UH(I,J))-ABS(UE(I,J)))**2
  328.           DUHL = DUHL + (ABS(UH(I,J))-ABS(UL(I,J)))**2
  329.           DUEL = DUEL + (ABS(UE(I,J))-ABS(UL(I,J)))**2
  330.    60   CONTINUE
  331.    70 CONTINUE
  332.       DUHE = SQRT(DUHE)
  333.       DUHL = SQRT(DUHL)
  334.       DUEL = SQRT(DUEL)
  335.       WRITE (NOUT,99998) DUHE, DUHL, DUEL
  336. 99998 FORMAT (/35H  2-NORM ( HYBSVD U - EISPAK U ) = , 1PE15.7/6H  2-NO,
  337.      * 29HRM ( HYBSVD U - LINPAK U ) = , 1PE15.7/19H  2-NORM ( EISPAK U,
  338.      * 16H - LINPAK U ) = , 1PE15.7/)
  339. C
  340.    80 IF (.NOT.MATV) GO TO 110
  341.       DO 100 J=1,MINMN
  342.         DO 90 I=1,N
  343.           DVHE = DVHE + (ABS(VH(I,J))-ABS(VE(I,J)))**2
  344.           DVHL = DVHL + (ABS(VH(I,J))-ABS(VL(I,J)))**2
  345.           DVEL = DVEL + (ABS(VE(I,J))-ABS(VL(I,J)))**2
  346.    90   CONTINUE
  347.   100 CONTINUE
  348.       DVHE = SQRT(DVHE)
  349.       DVHL = SQRT(DVHL)
  350.       DVEL = SQRT(DVEL)
  351.       WRITE (NOUT,99997) DVHE, DVHL, DVEL
  352. 99997 FORMAT (/35H  2-NORM ( HYBSVD V - EISPAK V ) = , 1PE15.7/6H  2-NO,
  353.      * 29HRM ( HYBSVD V - LINPAK V ) = , 1PE15.7/19H  2-NORM ( EISPAK V,
  354.      * 16H - LINPAK V ) = , 1PE15.7/)
  355. C
  356.   110 IF (M.LT.N .OR. IRHS.LT.1) GO TO 140
  357.       DO 130 J=1,IRHS
  358.         DO 120 I=1,MINMN
  359.           DBHE = DBHE + (ABS(BE(I,J))-ABS(BH(I,J)))**2
  360.   120   CONTINUE
  361.   130 CONTINUE
  362.       DBHE = SQRT(DBHE)
  363.       WRITE (NOUT,99996) DBHE
  364. 99996 FORMAT (/35H  2-NORM ( HYBSVD B - MINFIT B ) = , 1PE15.7/)
  365. C
  366.   140 CONTINUE
  367.       RETURN
  368.       END
  369.       SUBROUTINE HYBSVD(NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV,
  370.      * V, Z, B, IRHS, IERR, RV1)
  371.       INTEGER NA, NU, NV, NZ, M, N, IRHS, IERR, MIN0
  372.       REAL A(NA,1), W(1), U(NU,1), V(NV,1), Z(NZ,1), B(NB,IRHS), RV1(1)
  373.       LOGICAL MATU, MATV
  374. C
  375. C     THIS ROUTINE IS A MODIFICATION OF THE GOLUB-REINSCH PROCEDURE (1)
  376. C                                                           T
  377. C     FOR COMPUTING THE SINGULAR VALUE DECOMPOSITION A = UWV  OF A
  378. C     REAL M BY N RECTANGULAR MATRIX. U IS M BY MIN(M,N) CONTAINING
  379. C     THE LEFT SINGULAR VECTORS, W IS A MIN(M,N) BY MIN(M,N) DIAGONAL
  380. C     MATRIX CONTAINING THE SINGULAR VALUES, AND V IS N BY MIN(M,N)
  381. C     CONTAINING THE RIGHT SINGULAR VECTORS.
  382. C
  383. C     THE ALGORITHM IMPLEMENTED IN THIS
  384. C     ROUTINE HAS A HYBRID NATURE.  WHEN M IS APPROXIMATELY EQUAL TO N,
  385. C     THE GOLUB-REINSCH ALGORITHM IS USED, BUT WHEN EITHER OF THE RATIOS
  386. C     M/N OR N/M IS GREATER THAN ABOUT 2,
  387. C     A MODIFIED VERSION OF THE GOLUB-REINSCH
  388. C     ALGORITHM IS USED.  THIS MODIFIED ALGORITHM FIRST TRANSFORMS A
  389. C                                                                T
  390. C     INTO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS L
  391. C     AND THEN USES THE GOLUB-REINSCH ALGORITHM TO FIND THE SINGULAR
  392. C     VALUE DECOMPOSITION OF THE RESULTING UPPER TRIANGULAR MATRIX R.
  393. C     WHEN U IS NEEDED EXPLICITLY IN THE CASE M.GE.N (OR V IN THE CASE
  394. C     M.LT.N), AN EXTRA ARRAY Z (OF SIZE AT LEAST
  395. C     MIN(M,N)**2) IS NEEDED, BUT OTHERWISE Z IS NOT REFERENCED
  396. C     AND NO EXTRA STORAGE IS REQUIRED.  THIS HYBRID METHOD
  397. C     SHOULD BE MORE EFFICIENT THAN THE GOLUB-REINSCH ALGORITHM WHEN
  398. C     M/N OR N/M IS LARGE.  FOR DETAILS, SEE (2).
  399. C
  400. C     WHEN M .GE. N,
  401. C     HYBSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST
  402. C     SQUARES SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X=B.
  403. C     IF M .LT. N (I.E. FOR UNDERDETERMINED SYSTEMS), THE RHS B
  404. C     IS NOT PROCESSED.
  405. C
  406. C     NOTICE THAT THE SINGULAR VALUE DECOMPOSITION OF A MATRIX
  407. C     IS UNIQUE ONLY UP TO THE SIGN OF THE CORRESPONDING COLUMNS
  408. C     OF U AND V.
  409. C
  410. C     THIS ROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER (3) FOR
  411. C     ADHERENCE TO A LARGE, CAREFULLY DEFINED, PORTABLE SUBSET OF
  412. C     AMERICAN NATIONAL STANDARD FORTRAN CALLED PFORT.
  413. C
  414. C     REFERENCES:
  415. C
  416. C     (1) GOLUB,G.H. AND REINSCH,C. (1970) 'SINGULAR VALUE
  417. C         DECOMPOSITION AND LEAST SQUARES SOLUTIONS,'
  418. C         NUMER. MATH. 14,403-420, 1970.
  419. C
  420. C     (2) CHAN,T.F. (1982) 'AN IMPROVED ALGORITHM FOR COMPUTING
  421. C         THE SINGULAR VALUE DECOMPOSITION,' ACM TOMS, VOL.8,
  422. C         NO. 1, MARCH, 1982.
  423. C
  424. C     (3) RYDER,B.G. (1974) 'THE PFORT VERIFIER,' SOFTWARE -
  425. C         PRACTICE AND EXPERIENCE, VOL.4, 359-377, 1974.
  426. C
  427. C     ON INPUT:
  428. C
  429. C        NA MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  430. C          ARRAY PARAMETER A AS DECLARED IN THE CALLING PROGRAM
  431. C          DIMENSION STATEMENT.  NOTE THAT NA MUST BE AT LEAST
  432. C          AS LARGE AS M.
  433. C
  434. C        NU MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  435. C          ARRAY U AS DECLARED IN THE CALLING PROGRAM DIMENSION
  436. C          STATEMENT. NU MUST BE AT LEAST AS LARGE AS M.
  437. C
  438. C        NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  439. C          ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM
  440. C          DIMENSION STATEMENT. NV MUST BE AT LEAST AS LARGE AS N.
  441. C
  442. C        NZ MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  443. C          ARRAY PARAMETER Z AS DECLARED IN THE CALLING PROGRAM
  444. C          DIMENSION STATEMENT.  NOTE THAT NZ MUST BE AT LEAST
  445. C          AS LARGE AS MIN(M,N).
  446. C
  447. C        NB MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  448. C          ARRAY PARAMETER B AS DECLARED IN THE CALLING PROGRAM
  449. C          DIMENSION STATEMENT. NB MUST BE AT LEAST AS LARGE AS M.
  450. C
  451. C        M IS THE NUMBER OF ROWS OF A (AND U).
  452. C
  453. C        N IS THE NUMBER OF COLUMNS OF A (AND NUMBER OF ROWS OF V).
  454. C
  455. C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED.
  456. C
  457. C        B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED
  458. C         LINEAR SYSTEM A*X=B. IF IRHS .GT. 0 AND M .GE. N,
  459. C         THEN ON OUTPUT, THE FIRST N COMPONENTS OF THESE IRHS COLUMNS
  460. C                       T
  461. C         WILL CONTAIN U B. THUS, TO COMPUTE THE MINIMAL LENGTH LEAST
  462. C                                               +
  463. C         SQUARES SOLUTION, ONE MUST COMPUTE V*W  TIMES THE COLUMNS OF
  464. C                   +                        +
  465. C         B, WHERE W  IS A DIAGONAL MATRIX, W (I)=0 IF W(I) IS
  466. C         NEGLIGIBLE, OTHERWISE IS 1/W(I). IF IRHS=0 OR M.LT.N,
  467. C         B IS NOT REFERENCED.
  468. C
  469. C        IRHS IS THE NUMBER OF RIGHT-HAND-SIDES OF THE OVERDETERMINED
  470. C         SYSTEM A*X=B. IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULAR
  471. C         VALUE DECOMPOSITION OF A IS DESIRED.
  472. C
  473. C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
  474. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  475. C
  476. C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
  477. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  478. C
  479. C        WHEN HYBSVD IS USED TO COMPUTE THE MINIMAL LENGTH LEAST
  480. C        SQUARES SOLUTION TO AN OVERDETERMINED SYSTEM, MATU SHOULD
  481. C        BE SET TO .FALSE. , AND MATV SHOULD BE SET TO .TRUE.  .
  482. C
  483. C     ON OUTPUT:
  484. C
  485. C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V).
  486. C
  487. C        W CONTAINS THE (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  488. C          DIAGONAL ELEMENTS OF W).  THEY ARE SORTED IN DESCENDING
  489. C          ORDER.  IF AN ERROR EXIT IS MADE, THE SINGULAR VALUES
  490. C          SHOULD BE CORRECT AND SORTED FOR INDICES IERR+1,...,MIN(M,N).
  491. C
  492. C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
  493. C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  IF MATU IS
  494. C          FALSE, THEN U IS EITHER USED AS A TEMPORARY STORAGE (IF
  495. C          M .GE. N) OR NOT REFERENCED (IF M .LT. N).
  496. C          U MAY COINCIDE WITH A IN THE CALLING SEQUENCE.
  497. C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
  498. C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT.
  499. C
  500. C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
  501. C          MATV HAS BEEN SET TO .TRUE.  IF MATV IS
  502. C          FALSE, THEN V IS EITHER USED AS A TEMPORARY STORAGE (IF
  503. C          M .LT. N) OR NOT REFERENCED (IF M .GE. N).
  504. C          IF M .GE. N, V MAY ALSO COINCIDE WITH A.  IF AN ERROR
  505. C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
  506. C          CORRECT SINGULAR VALUES SHOULD BE CORRECT.
  507. C
  508. C        Z CONTAINS THE MATRIX X IN THE SINGULAR VALUE DECOMPOSITION
  509. C                  T
  510. C          OF R=XSY,  IF THE MODIFIED ALGORITHM IS USED. IF THE
  511. C          GOLUB-REINSCH PROCEDURE IS USED, THEN IT IS NOT REFERENCED.
  512. C          IF MATU HAS BEEN SET TO .FALSE. IN THE CASE M.GE.N (OR
  513. C          MATV SET TO .FALSE. IN THE CASE M.LT.N), THEN Z IS NOT
  514. C          REFERENCED AND NO EXTRA STORAGE IS REQUIRED.
  515. C
  516. C        IERR IS SET TO
  517. C          ZERO       FOR NORMAL RETURN,
  518. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  519. C                     DETERMINED AFTER 30 ITERATIONS.
  520. C          -1         IF IRHS .LT. 0 .
  521. C          -2         IF M .LT. 1 .OR. N .LT. 1
  522. C          -3         IF NA .LT. M .OR. NU .LT. M .OR. NB .LT. M.
  523. C          -4         IF NV .LT. N .
  524. C          -5         IF NZ .LT. MIN(M,N).
  525. C
  526. C        RV1 IS A TEMPORARY STORAGE ARRAY OF LENGTH AT LEAST MIN(M,N).
  527. C
  528. C     PROGRAMMED BY : TONY CHAN
  529. C                     BOX 2158, YALE STATION,
  530. C                     COMPUTER SCIENCE DEPT, YALE UNIV.,
  531. C                     NEW HAVEN, CT 06520.
  532. C     LAST MODIFIED : JANUARY, 1982.
  533. C
  534. C     HYBSVD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES.
  535. C       INTERNAL  GRSVD, MGNSVD, SRELPR
  536. C       FORTRAN   MIN0,ABS,SQRT,FLOAT,SIGN,AMAX1
  537. C       BLAS      SSWAP
  538. C
  539. C     -----------------------------------------------------------------
  540. C     ERROR CHECK.
  541. C
  542.       IERR = 0
  543.       IF (IRHS.GE.0) GO TO 10
  544.       IERR = -1
  545.       RETURN
  546.    10 IF (M.GE.1 .AND. N.GE.1) GO TO 20
  547.       IERR = -2
  548.       RETURN
  549.    20 IF (NA.GE.M .AND. NU.GE.M .AND. NB.GE.M) GO TO 30
  550.       IERR = -3
  551.       RETURN
  552.    30 IF (NV.GE.N) GO TO 40
  553.       IERR = -4
  554.       RETURN
  555.    40 IF (NZ.GE.MIN0(M,N)) GO TO 50
  556.       IERR = -5
  557.       RETURN
  558.    50 CONTINUE
  559. C
  560. C     FIRST COPIES A INTO EITHER U OR V ACCORDING TO WHETHER
  561. C     M .GE. N OR M .LT. N, AND THEN CALLS SUBROUTINE MGNSVD
  562. C     WHICH ASSUMES THAT NUMBER OF ROWS .GE. NUMBER OF COLUMNS.
  563. C
  564.       IF (M.LT.N) GO TO 80
  565. C
  566. C       M .GE. N  CASE.
  567. C
  568.       DO 70 I=1,M
  569.         DO 60 J=1,N
  570.           U(I,J) = A(I,J)
  571.    60   CONTINUE
  572.    70 CONTINUE
  573. C
  574.       CALL MGNSVD(NU, NV, NZ, NB, M, N, W, MATU, U, MATV, V, Z, B,
  575.      * IRHS, IERR, RV1)
  576.       RETURN
  577. C
  578.    80 CONTINUE
  579. C                              T
  580. C       M .LT. N CASE. COPIES A  INTO V.
  581. C
  582.       DO 100 I=1,M
  583.         DO 90 J=1,N
  584.           V(J,I) = A(I,J)
  585.    90   CONTINUE
  586.   100 CONTINUE
  587.       CALL MGNSVD(NV, NU, NZ, NB, N, M, W, MATV, V, MATU, U, Z, B, 0,
  588.      * IERR, RV1)
  589.       RETURN
  590.       END
  591. C
  592.       SUBROUTINE MGNSVD(NU, NV, NZ, NB, M, N, W, MATU, U, MATV, V, Z,
  593.      * B, IRHS, IERR, RV1)
  594. C
  595. C     THE DESCRIPTION OF SUBROUTINE MGNSVD IS ALMOST IDENTICAL
  596. C     TO THAT FOR SUBROUTINE HYBSVD ABOVE, WITH THE EXCEPTION
  597. C     THAT MGNSVD ASSUMES M .GE. N.
  598. C     IT ALSO ASSUMES THAT A COPY OF THE MATRIX A IS IN THE ARRAY U.
  599. C
  600.       INTEGER NU, NV, NZ, M, N, IRHS, IERR, IP1, I, J, K, IM1, IBACK
  601.       REAL W(1), U(NU,1), V(NV,1), Z(NZ,1), B(NB,IRHS), RV1(1)
  602.       REAL XOVRPT, C, R, G, SCALE, SIGN, ABS, SQRT, F, S, H
  603.       REAL FLOAT
  604.       LOGICAL MATU, MATV
  605. C
  606. C     SET VALUE FOR C. THE VALUE FOR C DEPENDS ON THE RELATIVE
  607. C     EFFICIENCY OF FLOATING POINT MULTIPLICATIONS, FLOATING POINT
  608. C     ADDITIONS AND TWO-DIMENSIONAL ARRAY INDEXINGS ON THE
  609. C     COMPUTER WHERE THIS SUBROUTINE IS TO BE RUN.  C SHOULD
  610. C     USUALLY BE BETWEEN 2 AND 4.  FOR DETAILS ON CHOOSING C, SEE
  611. C     (2).  THE ALGORITHM IS NOT SENSITIVE TO THE VALUE OF C
  612. C     ACTUALLY USED AS LONG AS C IS BETWEEN 2 AND 4.
  613. C
  614.       C = 4.0
  615. C
  616. C     DETERMINE CROSS-OVER POINT
  617. C
  618.       IF (MATU .AND. MATV) XOVRPT = (C+7./3.)/C
  619.       IF (MATU .AND. .NOT.MATV) XOVRPT = (C+7./3.)/C
  620.       IF (.NOT.MATU .AND. MATV) XOVRPT = 5./3.
  621.       IF (.NOT.MATU .AND. .NOT.MATV) XOVRPT = 5./3.
  622. C
  623. C     DETERMINE WHETHER TO USE GOLUB-REINSCH OR THE MODIFIED
  624. C     ALGORITHM.
  625. C
  626.       R = FLOAT(M)/FLOAT(N)
  627.       IF (R.GE.XOVRPT) GO TO 10
  628. C
  629. C     USE GOLUB-REINSCH PROCEDURE
  630. C
  631.       CALL GRSVD(NU, NV, NB, M, N, W, MATU, U, MATV, V, B, IRHS, IERR,
  632.      * RV1)
  633.       GO TO 330
  634. C
  635. C     USE MODIFIED ALGORITHM
  636. C
  637.    10 CONTINUE
  638. C
  639. C     TRIANGULARIZE U BY HOUSEHOLDER TRANSFORMATIONS, USING
  640. C     W AND RV1 AS TEMPORARY STORAGE.
  641. C
  642.       DO 110 I=1,N
  643.         G = 0.0
  644.         S = 0.0
  645.         SCALE = 0.0
  646. C
  647. C         PERFORM SCALING OF COLUMNS TO AVOID UNNECSSARY OVERFLOW
  648. C         OR UNDERFLOW
  649. C
  650.         DO 20 K=I,M
  651.           SCALE = SCALE + ABS(U(K,I))
  652.    20   CONTINUE
  653.         IF (SCALE.EQ.0.0) GO TO 110
  654.         DO 30 K=I,M
  655.           U(K,I) = U(K,I)/SCALE
  656.           S = S + U(K,I)**2
  657.    30   CONTINUE
  658. C
  659. C         THE VECTOR E OF THE HOUSEHOLDER TRANSFORMATION I + EE'/H
  660. C         WILL BE STORED IN COLUMN I OF U. THE TRANSFORMED ELEMENT
  661. C         U(I,I) WILL BE STORED IN W(I) AND THE SCALAR H IN
  662. C         RV1(I).
  663. C
  664.         F = U(I,I)
  665.         G = -SIGN(SQRT(S),F)
  666.         H = F*G - S
  667.         U(I,I) = F - G
  668.         RV1(I) = H
  669.         W(I) = SCALE*G
  670. C
  671.         IF (I.EQ.N) GO TO 70
  672. C
  673. C         APPLY TRANSFORMATIONS TO REMAINING COLUMNS OF A
  674. C
  675.         IP1 = I + 1
  676.         DO 60 J=IP1,N
  677.           S = 0.0
  678.           DO 40 K=I,M
  679.             S = S + U(K,I)*U(K,J)
  680.    40     CONTINUE
  681.           F = S/H
  682.           DO 50 K=I,M
  683.             U(K,J) = U(K,J) + F*U(K,I)
  684.    50     CONTINUE
  685.    60   CONTINUE
  686. C
  687. C         APPLY TRANSFORMATIONS TO COLUMNS OF B IF IRHS .GT. 0
  688. C
  689.    70   IF (IRHS.EQ.0) GO TO 110
  690.         DO 100 J=1,IRHS
  691.           S = 0.0
  692.           DO 80 K=I,M
  693.             S = S + U(K,I)*B(K,J)
  694.    80     CONTINUE
  695.           F = S/H
  696.           DO 90 K=I,M
  697.             B(K,J) = B(K,J) + F*U(K,I)
  698.    90     CONTINUE
  699.   100   CONTINUE
  700.   110 CONTINUE
  701. C
  702. C     COPY R INTO Z IF MATU = .TRUE.
  703. C
  704.       IF (.NOT.MATU) GO TO 290
  705.       DO 130 I=1,N
  706.         DO 120 J=I,N
  707.           Z(J,I) = 0.0
  708.           Z(I,J) = U(I,J)
  709.   120   CONTINUE
  710.         Z(I,I) = W(I)
  711.   130 CONTINUE
  712. C
  713. C     ACCUMULATE HOUSEHOLDER TRANSFORMATIONS IN U
  714. C
  715.       DO 240 IBACK=1,N
  716.         I = N - IBACK + 1
  717.         IP1 = I + 1
  718.         G = W(I)
  719.         H = RV1(I)
  720.         IF (I.EQ.N) GO TO 150
  721. C
  722.         DO 140 J=IP1,N
  723.           U(I,J) = 0.0
  724.   140   CONTINUE
  725. C
  726.   150   IF (H.EQ.0.0) GO TO 210
  727.         IF (I.EQ.N) GO TO 190
  728. C
  729.         DO 180 J=IP1,N
  730.           S = 0.0
  731.           DO 160 K=IP1,M
  732.             S = S + U(K,I)*U(K,J)
  733.   160     CONTINUE
  734.           F = S/H
  735.           DO 170 K=I,M
  736.             U(K,J) = U(K,J) + F*U(K,I)
  737.   170     CONTINUE
  738.   180   CONTINUE
  739. C
  740.   190   S = U(I,I)/H
  741.         DO 200 J=I,M
  742.           U(J,I) = U(J,I)*S
  743.   200   CONTINUE
  744.         GO TO 230
  745. C
  746.   210   DO 220 J=I,M
  747.           U(J,I) = 0.0
  748.   220   CONTINUE
  749.   230   U(I,I) = U(I,I) + 1.0
  750.   240 CONTINUE
  751. C
  752. C     COMPUTE SVD OF R (WHICH IS STORED IN Z)
  753. C
  754.       CALL GRSVD(NZ, NV, NB, N, N, W, MATU, Z, MATV, V, B, IRHS, IERR,
  755.      * RV1)
  756. C
  757. C                                      T
  758. C     FORM L*X TO OBTAIN U (WHERE R=XWY ). X IS RETURNED IN Z
  759. C     BY GRSVD. THE MATRIX MULTIPLY IS DONE ONE ROW AT A TIME,
  760. C     USING RV1 AS SCRATCH SPACE.
  761. C
  762.       DO 280 I=1,M
  763.         DO 260 J=1,N
  764.           S = 0.0
  765.           DO 250 K=1,N
  766.             S = S + U(I,K)*Z(K,J)
  767.   250     CONTINUE
  768.           RV1(J) = S
  769.   260   CONTINUE
  770.         DO 270 J=1,N
  771.           U(I,J) = RV1(J)
  772.   270   CONTINUE
  773.   280 CONTINUE
  774.       GO TO 330
  775. C
  776. C     FORM R IN U BY ZEROING THE LOWER TRIANGULAR PART OF R IN U
  777. C
  778.   290 IF (N.EQ.1) GO TO 320
  779.       DO 310 I=2,N
  780.         IM1 = I - 1
  781.         DO 300 J=1,IM1
  782.           U(I,J) = 0.0
  783.   300   CONTINUE
  784.         U(I,I) = W(I)
  785.   310 CONTINUE
  786.   320 U(1,1) = W(1)
  787. C
  788.       CALL GRSVD(NU, NV, NB, N, N, W, MATU, U, MATV, V, B, IRHS, IERR,
  789.      * RV1)
  790.   330 CONTINUE
  791.       IERRP1 = IERR + 1
  792.       IF (IERR.LT.0 .OR. N.LE.1 .OR. IERRP1.EQ.N) RETURN
  793. C
  794. C     SORT SINGULAR VALUES AND EXCHANGE COLUMNS OF U AND V ACCORDINGLY.
  795. C     SELECTION SORT MINIMIZES SWAPPING OF U AND V.
  796. C
  797.       NM1 = N - 1
  798.       DO 360 I=IERRP1,NM1
  799. C...    FIND INDEX OF MAXIMUM SINGULAR VALUE
  800.         ID = I
  801.         IP1 = I + 1
  802.         DO 340 J=IP1,N
  803.           IF (W(J).GT.W(ID)) ID = J
  804.   340   CONTINUE
  805.         IF (ID.EQ.I) GO TO 360
  806. C...    SWAP SINGULAR VALUES AND VECTORS
  807.         T = W(I)
  808.         W(I) = W(ID)
  809.         W(ID) = T
  810.         IF (MATV) CALL SSWAP(N, V(1,I), 1, V(1,ID), 1)
  811.         IF (MATU) CALL SSWAP(M, U(1,I), 1, U(1,ID), 1)
  812.         IF (IRHS.LT.1) GO TO 360
  813.         DO 350 KRHS=1,IRHS
  814.           T = B(I,KRHS)
  815.           B(I,KRHS) = B(ID,KRHS)
  816.           B(ID,KRHS) = T
  817.   350   CONTINUE
  818.   360 CONTINUE
  819.       RETURN
  820. C     ************** LAST CARD OF HYBSVD *****************
  821.       END
  822.       SUBROUTINE GRSVD(NU, NV, NB, M, N, W, MATU, U, MATV, V, B, IRHS,
  823.      * IERR, RV1)
  824. C
  825.       INTEGER I, J, K, L, M, N, II, I1, KK, K1, LL, L1, MN, NU, NV, NB,
  826.      * ITS, IERR, IRHS
  827.       REAL W(1), U(NU,1), V(NV,1), B(NB,IRHS), RV1(1)
  828.       REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, SRELPR, DUMMY
  829.       REAL SQRT, AMAX1, ABS, SIGN
  830.       LOGICAL MATU, MATV
  831. C
  832. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD,
  833. C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
  834. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
  835. C
  836. C     THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION
  837. C          T
  838. C     A=USV  OF A REAL M BY N RECTANGULAR MATRIX.  HOUSEHOLDER
  839. C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
  840. C     GRSVD ASSUMES THAT A COPY OF THE MATRIX A IS IN THE ARRAY U. IT
  841. C     ALSO ASSUMES M .GE. N.  IF M .LT. N, THEN COMPUTE THE SINGULAR
  842. C                             T       T    T             T
  843. C     VALUE DECOMPOSITION OF A .  IF A =UWV  , THEN A=VWU  .
  844. C
  845. C     GRSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST SQUARES
  846. C     SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X=B.
  847. C
  848. C     ON INPUT-
  849. C
  850. C        NU MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  851. C          ARRAY PARAMETERS U AS DECLARED IN THE CALLING PROGRAM
  852. C          DIMENSION STATEMENT.  NOTE THAT NU MUST BE AT LEAST
  853. C          AS LARGE AS M,
  854. C
  855. C        NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  856. C          ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM
  857. C          DIMENSION STATEMENT.  NV MUST BE AT LEAST AS LARGE AS N,
  858. C
  859. C        NB MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  860. C          ARRAY PARAMETERS B AS DECLARED IN THE CALLING PROGRAM
  861. C          DIMENSION STATEMENT.  NOTE THAT NB MUST BE AT LEAST
  862. C          AS LARGE AS M,
  863. C
  864. C        M IS THE NUMBER OF ROWS OF A (AND U),
  865. C
  866. C        N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V,
  867. C
  868. C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED,
  869. C
  870. C        B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED
  871. C          LINEAR SYSTEM A*X=B.  IF IRHS .GT. 0,  THEN ON OUTPUT,
  872. C          THE FIRST N COMPONENTS OF THESE IRHS COLUMNS OF B
  873. C                        T
  874. C          WILL CONTAIN U B.  THUS, TO COMPUTE THE MINIMAL LENGTH LEAST
  875. C                                                +
  876. C          SQUARES SOLUTION, ONE MUST COMPUTE V*W  TIMES THE COLUMNS OF
  877. C                    +                        +
  878. C          B, WHERE W  IS A DIAGONAL MATRIX, W (I)=0 IF W(I) IS
  879. C          NEGLIGIBLE, OTHERWISE IS 1/W(I).  IF IRHS=0, B MAY COINCIDE
  880. C          WITH A OR U AND WILL NOT BE REFERENCED,
  881. C
  882. C        IRHS IS THE NUMBER OF RIGHT-HAND-SIDES OF THE OVERDETERMINED
  883. C          SYSTEM A*X=B.  IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULA
  884. C          VALUE DECOMPOSITION OF A IS DESIRED,
  885. C
  886. C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
  887. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE,
  888. C
  889. C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
  890. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  891. C
  892. C     ON OUTPUT-
  893. C
  894. C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  895. C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
  896. C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
  897. C          FOR INDICES IERR+1,IERR+2,...,N,
  898. C
  899. C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
  900. C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  OTHERWISE
  901. C          U IS USED AS A TEMPORARY ARRAY.
  902. C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
  903. C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
  904. C
  905. C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
  906. C          MATV HAS BEEN SET TO .TRUE.  OTHERWISE V IS NOT REFERENCED.
  907. C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO
  908. C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
  909. C
  910. C        IERR IS SET TO
  911. C          ZERO       FOR NORMAL RETURN,
  912. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  913. C                     DETERMINED AFTER 30 ITERATIONS,
  914. C          -1         IF IRHS .LT. 0 ,
  915. C          -2         IF M .LT. N ,
  916. C          -3         IF NU .LT. M .OR. NB .LT. M,
  917. C          -4         IF NV .LT. N .
  918. C
  919. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  920. C
  921. C        THIS SUBROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER
  922. C        (RYDER, B.G. 'THE PFORT VERIFIER', SOFTWARE - PRACTICE AND
  923. C        EXPERIENCE, VOL.4, 359-377, 1974) FOR ADHERENCE TO A LARGE,
  924. C        CAREFULLY DEFINED, PORTABLE SUBSET OF AMERICAN NATIONAL STANDAR
  925. C        FORTRAN CALLED PFORT.
  926. C
  927. C        ORIGINAL VERSION OF THIS CODE IS SUBROUTINE SVD IN RELEASE 2 OF
  928. C        EISPACK.
  929. C
  930. C        MODIFIED BY TONY F. CHAN,
  931. C                    COMP. SCI. DEPT, YALE UNIV.,
  932. C                    BOX 2158, YALE STATION,
  933. C                    CT 06520
  934. C        LAST MODIFIED : JANUARY, 1982.
  935. C
  936. C     ------------------------------------------------------------------
  937. C
  938. C     ********** SRELPR IS A MACHINE-DEPENDENT FUNCTION SPECIFYING
  939. C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
  940. C
  941. C                **********
  942. C
  943.       IERR = 0
  944.       IF (IRHS.GE.0) GO TO 10
  945.       IERR = -1
  946.       RETURN
  947.    10 IF (M.GE.N) GO TO 20
  948.       IERR = -2
  949.       RETURN
  950.    20 IF (NU.GE.M .AND. NB.GE.M) GO TO 30
  951.       IERR = -3
  952.       RETURN
  953.    30 IF (NV.GE.N) GO TO 40
  954.       IERR = -4
  955.       RETURN
  956.    40 CONTINUE
  957. C
  958. C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
  959.       G = 0.0
  960.       SCALE = 0.0
  961.       X = 0.0
  962. C
  963.       DO 260 I=1,N
  964.         L = I + 1
  965.         RV1(I) = SCALE*G
  966.         G = 0.0
  967.         S = 0.0
  968.         SCALE = 0.0
  969. C
  970. C     COMPUTE LEFT TRANSFORMATIONS THAT ZERO THE SUBDIAGONAL ELEMENTS
  971. C     OF THE I-TH COLUMN.
  972. C
  973.         DO 50 K=I,M
  974.           SCALE = SCALE + ABS(U(K,I))
  975.    50   CONTINUE
  976. C
  977.         IF (SCALE.EQ.0.0) GO TO 160
  978. C
  979.         DO 60 K=I,M
  980.           U(K,I) = U(K,I)/SCALE
  981.           S = S + U(K,I)**2
  982.    60   CONTINUE
  983. C
  984.         F = U(I,I)
  985.         G = -SIGN(SQRT(S),F)
  986.         H = F*G - S
  987.         U(I,I) = F - G
  988.         IF (I.EQ.N) GO TO 100
  989. C
  990. C     APPLY LEFT TRANSFORMATIONS TO REMAINING COLUMNS OF A.
  991. C
  992.         DO 90 J=L,N
  993.           S = 0.0
  994. C
  995.           DO 70 K=I,M
  996.             S = S + U(K,I)*U(K,J)
  997.    70     CONTINUE
  998. C
  999.           F = S/H
  1000. C
  1001.           DO 80 K=I,M
  1002.             U(K,J) = U(K,J) + F*U(K,I)
  1003.    80     CONTINUE
  1004.    90   CONTINUE
  1005. C
  1006. C      APPLY LEFT TRANSFORMATIONS TO THE COLUMNS OF B IF IRHS .GT. 0
  1007. C
  1008.   100   IF (IRHS.EQ.0) GO TO 140
  1009.         DO 130 J=1,IRHS
  1010.           S = 0.0
  1011.           DO 110 K=I,M
  1012.             S = S + U(K,I)*B(K,J)
  1013.   110     CONTINUE
  1014.           F = S/H
  1015.           DO 120 K=I,M
  1016.             B(K,J) = B(K,J) + F*U(K,I)
  1017.   120     CONTINUE
  1018.   130   CONTINUE
  1019. C
  1020. C     COMPUTE RIGHT TRANSFORMATIONS.
  1021. C
  1022.   140   DO 150 K=I,M
  1023.           U(K,I) = SCALE*U(K,I)
  1024.   150   CONTINUE
  1025. C
  1026.   160   W(I) = SCALE*G
  1027.         G = 0.0
  1028.         S = 0.0
  1029.         SCALE = 0.0
  1030.         IF (I.GT.M .OR. I.EQ.N) GO TO 250
  1031. C
  1032.         DO 170 K=L,N
  1033.           SCALE = SCALE + ABS(U(I,K))
  1034.   170   CONTINUE
  1035. C
  1036.         IF (SCALE.EQ.0.0) GO TO 250
  1037. C
  1038.         DO 180 K=L,N
  1039.           U(I,K) = U(I,K)/SCALE
  1040.           S = S + U(I,K)**2
  1041.   180   CONTINUE
  1042. C
  1043.         F = U(I,L)
  1044.         G = -SIGN(SQRT(S),F)
  1045.         H = F*G - S
  1046.         U(I,L) = F - G
  1047. C
  1048.         DO 190 K=L,N
  1049.           RV1(K) = U(I,K)/H
  1050.   190   CONTINUE
  1051. C
  1052.         IF (I.EQ.M) GO TO 230
  1053. C
  1054.         DO 220 J=L,M
  1055.           S = 0.0
  1056. C
  1057.           DO 200 K=L,N
  1058.             S = S + U(J,K)*U(I,K)
  1059.   200     CONTINUE
  1060. C
  1061.           DO 210 K=L,N
  1062.             U(J,K) = U(J,K) + S*RV1(K)
  1063.   210     CONTINUE
  1064.   220   CONTINUE
  1065. C
  1066.   230   DO 240 K=L,N
  1067.           U(I,K) = SCALE*U(I,K)
  1068.   240   CONTINUE
  1069. C
  1070.   250   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  1071.   260 CONTINUE
  1072. C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS **********
  1073.       IF (.NOT.MATV) GO TO 350
  1074. C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
  1075.       DO 340 II=1,N
  1076.         I = N + 1 - II
  1077.         IF (I.EQ.N) GO TO 330
  1078.         IF (G.EQ.0.0) GO TO 310
  1079. C
  1080.         DO 270 J=L,N
  1081. C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
  1082.           V(J,I) = (U(I,J)/U(I,L))/G
  1083.   270   CONTINUE
  1084. C
  1085.         DO 300 J=L,N
  1086.           S = 0.0
  1087. C
  1088.           DO 280 K=L,N
  1089.             S = S + U(I,K)*V(K,J)
  1090.   280     CONTINUE
  1091. C
  1092.           DO 290 K=L,N
  1093.             V(K,J) = V(K,J) + S*V(K,I)
  1094.   290     CONTINUE
  1095.   300   CONTINUE
  1096. C
  1097.   310   DO 320 J=L,N
  1098.           V(I,J) = 0.0
  1099.           V(J,I) = 0.0
  1100.   320   CONTINUE
  1101. C
  1102.   330   V(I,I) = 1.0
  1103.         G = RV1(I)
  1104.         L = I
  1105.   340 CONTINUE
  1106. C     ********** ACCUMULATION OF LEFT-HAND TRANSFORMATIONS **********
  1107.   350 IF (.NOT.MATU) GO TO 470
  1108. C     **********FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- **********
  1109.       MN = N
  1110.       IF (M.LT.N) MN = M
  1111. C
  1112.       DO 460 II=1,MN
  1113.         I = MN + 1 - II
  1114.         L = I + 1
  1115.         G = W(I)
  1116.         IF (I.EQ.N) GO TO 370
  1117. C
  1118.         DO 360 J=L,N
  1119.           U(I,J) = 0.0
  1120.   360   CONTINUE
  1121. C
  1122.   370   IF (G.EQ.0.0) GO TO 430
  1123.         IF (I.EQ.MN) GO TO 410
  1124. C
  1125.         DO 400 J=L,N
  1126.           S = 0.0
  1127. C
  1128.           DO 380 K=L,M
  1129.             S = S + U(K,I)*U(K,J)
  1130.   380     CONTINUE
  1131. C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
  1132.           F = (S/U(I,I))/G
  1133. C
  1134.           DO 390 K=I,M
  1135.             U(K,J) = U(K,J) + F*U(K,I)
  1136.   390     CONTINUE
  1137.   400   CONTINUE
  1138. C
  1139.   410   DO 420 J=I,M
  1140.           U(J,I) = U(J,I)/G
  1141.   420   CONTINUE
  1142. C
  1143.         GO TO 450
  1144. C
  1145.   430   DO 440 J=I,M
  1146.           U(J,I) = 0.0
  1147.   440   CONTINUE
  1148. C
  1149.   450   U(I,I) = U(I,I) + 1.0
  1150.   460 CONTINUE
  1151. C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  1152.   470 EPS = SRELPR(DUMMY)*X
  1153. C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
  1154.       DO 650 KK=1,N
  1155.         K1 = N - KK
  1156.         K = K1 + 1
  1157.         ITS = 0
  1158. C     ********** TEST FOR SPLITTING.
  1159. C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  1160.   480   DO 490 LL=1,K
  1161.           L1 = K - LL
  1162.           L = L1 + 1
  1163.           IF (ABS(RV1(L)).LE.EPS) GO TO 550
  1164. C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
  1165. C                THROUGH THE BOTTOM OF THE LOOP **********
  1166.           IF (ABS(W(L1)).LE.EPS) GO TO 500
  1167.   490   CONTINUE
  1168. C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  1169.   500   C = 0.0
  1170.         S = 1.0
  1171. C
  1172.         DO 540 I=L,K
  1173.           F = S*RV1(I)
  1174.           RV1(I) = C*RV1(I)
  1175.           IF (ABS(F).LE.EPS) GO TO 550
  1176.           G = W(I)
  1177.           H = SQRT(F*F+G*G)
  1178.           W(I) = H
  1179.           C = G/H
  1180.           S = -F/H
  1181. C
  1182. C     APPLY LEFT TRANSFORMATIONS TO B IF IRHS .GT. 0
  1183. C
  1184.           IF (IRHS.EQ.0) GO TO 520
  1185.           DO 510 J=1,IRHS
  1186.             Y = B(L1,J)
  1187.             Z = B(I,J)
  1188.             B(L1,J) = Y*C + Z*S
  1189.             B(I,J) = -Y*S + Z*C
  1190.   510     CONTINUE
  1191.   520     CONTINUE
  1192. C
  1193.           IF (.NOT.MATU) GO TO 540
  1194. C
  1195.           DO 530 J=1,M
  1196.             Y = U(J,L1)
  1197.             Z = U(J,I)
  1198.             U(J,L1) = Y*C + Z*S
  1199.             U(J,I) = -Y*S + Z*C
  1200.   530     CONTINUE
  1201. C
  1202.   540   CONTINUE
  1203. C     ********** TEST FOR CONVERGENCE **********
  1204.   550   Z = W(K)
  1205.         IF (L.EQ.K) GO TO 630
  1206. C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
  1207.         IF (ITS.EQ.30) GO TO 660
  1208.         ITS = ITS + 1
  1209.         X = W(L)
  1210.         Y = W(K1)
  1211.         G = RV1(K1)
  1212.         H = RV1(K)
  1213.         F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
  1214.         G = SQRT(F*F+1.0)
  1215.         F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
  1216. C     ********** NEXT QR TRANSFORMATION **********
  1217.         C = 1.0
  1218.         S = 1.0
  1219. C
  1220.         DO 620 I1=L,K1
  1221.           I = I1 + 1
  1222.           G = RV1(I)
  1223.           Y = W(I)
  1224.           H = S*G
  1225.           G = C*G
  1226.           Z = SQRT(F*F+H*H)
  1227.           RV1(I1) = Z
  1228.           C = F/Z
  1229.           S = H/Z
  1230.           F = X*C + G*S
  1231.           G = -X*S + G*C
  1232.           H = Y*S
  1233.           Y = Y*C
  1234.           IF (.NOT.MATV) GO TO 570
  1235. C
  1236.           DO 560 J=1,N
  1237.             X = V(J,I1)
  1238.             Z = V(J,I)
  1239.             V(J,I1) = X*C + Z*S
  1240.             V(J,I) = -X*S + Z*C
  1241.   560     CONTINUE
  1242. C
  1243.   570     Z = SQRT(F*F+H*H)
  1244.           W(I1) = Z
  1245. C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
  1246.           IF (Z.EQ.0.0) GO TO 580
  1247.           C = F/Z
  1248.           S = H/Z
  1249.   580     F = C*G + S*Y
  1250.           X = -S*G + C*Y
  1251. C
  1252. C     APPLY LEFT TRANSFORMATIONS TO B IF IRHS .GT. 0
  1253. C
  1254.           IF (IRHS.EQ.0) GO TO 600
  1255.           DO 590 J=1,IRHS
  1256.             Y = B(I1,J)
  1257.             Z = B(I,J)
  1258.             B(I1,J) = Y*C + Z*S
  1259.             B(I,J) = -Y*S + Z*C
  1260.   590     CONTINUE
  1261.   600     CONTINUE
  1262. C
  1263.           IF (.NOT.MATU) GO TO 620
  1264. C
  1265.           DO 610 J=1,M
  1266.             Y = U(J,I1)
  1267.             Z = U(J,I)
  1268.             U(J,I1) = Y*C + Z*S
  1269.             U(J,I) = -Y*S + Z*C
  1270.   610     CONTINUE
  1271. C
  1272.   620   CONTINUE
  1273. C
  1274.         RV1(L) = 0.0
  1275.         RV1(K) = F
  1276.         W(K) = X
  1277.         GO TO 480
  1278. C     ********** CONVERGENCE **********
  1279.   630   IF (Z.GE.0.0) GO TO 650
  1280. C     ********** W(K) IS MADE NON-NEGATIVE **********
  1281.         W(K) = -Z
  1282.         IF (.NOT.MATV) GO TO 650
  1283. C
  1284.         DO 640 J=1,N
  1285.           V(J,K) = -V(J,K)
  1286.   640   CONTINUE
  1287. C
  1288.   650 CONTINUE
  1289. C
  1290.       GO TO 670
  1291. C     ********** SET ERROR -- NO CONVERGENCE TO A
  1292. C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  1293.   660 IERR = K
  1294.   670 RETURN
  1295. C     ********** LAST CARD OF GRSVD **********
  1296.       END
  1297.       SUBROUTINE SSWAP(N, SX, INCX, SY, INCY)
  1298. C
  1299. C     INTERCHANGES TWO VECTORS.
  1300. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO 1.
  1301. C     JACK DONGARRA, LINPACK, 3/11/78.
  1302. C
  1303.       REAL SX(1), SY(1), STEMP
  1304.       INTEGER I, INCX, INCY, IX, IY, M, MP1, N
  1305. C
  1306.       IF (N.LE.0) RETURN
  1307.       IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
  1308. C
  1309. C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
  1310. C         TO 1
  1311. C
  1312.       IX = 1
  1313.       IY = 1
  1314.       IF (INCX.LT.0) IX = (-N+1)*INCX + 1
  1315.       IF (INCY.LT.0) IY = (-N+1)*INCY + 1
  1316.       DO 10 I=1,N
  1317.         STEMP = SX(IX)
  1318.         SX(IX) = SY(IY)
  1319.         SY(IY) = STEMP
  1320.         IX = IX + INCX
  1321.         IY = IY + INCY
  1322.    10 CONTINUE
  1323.       RETURN
  1324. C
  1325. C       CODE FOR BOTH INCREMENTS EQUAL TO 1
  1326. C
  1327. C
  1328. C       CLEAN-UP LOOP
  1329. C
  1330.    20 M = MOD(N,3)
  1331.       IF (M.EQ.0) GO TO 40
  1332.       DO 30 I=1,M
  1333.         STEMP = SX(I)
  1334.         SX(I) = SY(I)
  1335.         SY(I) = STEMP
  1336.    30 CONTINUE
  1337.       IF (N.LT.3) RETURN
  1338.    40 MP1 = M + 1
  1339.       DO 50 I=MP1,N,3
  1340.         STEMP = SX(I)
  1341.         SX(I) = SY(I)
  1342.         SY(I) = STEMP
  1343.         STEMP = SX(I+1)
  1344.         SX(I+1) = SY(I+1)
  1345.         SY(I+1) = STEMP
  1346.         STEMP = SX(I+2)
  1347.         SX(I+2) = SY(I+2)
  1348.         SY(I+2) = STEMP
  1349.    50 CONTINUE
  1350.       RETURN
  1351.       END
  1352.       REAL FUNCTION SRELPR(DUMMY)
  1353.       REAL DUMMY
  1354. C
  1355. C     SRELPR COMPUTES THE RELATIVE PRECISION OF THE FLOATING POINT
  1356. C     ARITHMETIC OF THE MACHINE.
  1357. C
  1358. C     IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES,
  1359. C     THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS.
  1360. C     ASSUME THE COMPUTER HAS
  1361. C
  1362. C        B = BASE OF ARITHMETIC
  1363. C        T = NUMBER OF BASE  B  DIGITS
  1364. C
  1365. C     THEN
  1366. C
  1367. C        SRELPR = B**(1-T)
  1368. C
  1369.       REAL S
  1370. C
  1371.       SRELPR = 1.0
  1372.    10 SRELPR = SRELPR/2.0
  1373.       S = 1.0 + SRELPR
  1374.       IF (S.GT.1.0) GO TO 10
  1375.       SRELPR = 2.0*SRELPR
  1376.       RETURN
  1377.       END
  1378.       SUBROUTINE SSVDC(X, LDX, N, P, S, E, U, LDU, V, LDV, WORK, JOB,
  1379.      * INFO)
  1380.       INTEGER LDX, N, P, LDU, LDV, JOB, INFO
  1381.       REAL X(LDX,1), S(1), E(1), U(LDU,1), V(LDV,1), WORK(1)
  1382. C
  1383. C
  1384. C     SSVDC IS A SUBROUTINE TO REDUCE A REAL NXP MATRIX X BY
  1385. C     ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
  1386. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
  1387. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
  1388. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
  1389. C
  1390. C     ON ENTRY
  1391. C
  1392. C         X         REAL(LDX,P), WHERE LDX.GE.N.
  1393. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
  1394. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
  1395. C                   DESTROYED BY SSVDC.
  1396. C
  1397. C         LDX       INTEGER.
  1398. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  1399. C
  1400. C         N         INTEGER.
  1401. C                   N IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  1402. C
  1403. C         P         INTEGER.
  1404. C                   P IS THE NUMBER OF ROWS OF THE MATRIX X.
  1405. C
  1406. C         LDU       INTEGER.
  1407. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
  1408. C                   (SEE BELOW).
  1409. C
  1410. C         LDV       INTEGER.
  1411. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
  1412. C                   (SEE BELOW).
  1413. C
  1414. C         WORK      REAL(N).
  1415. C                   WORK IS A SCRATCH ARRAY.
  1416. C
  1417. C         JOB       INTEGER.
  1418. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
  1419. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
  1420. C                   WITH THE FOLLOWING MEANING
  1421. C
  1422. C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
  1423. C                                  VECTORS.
  1424. C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
  1425. C                                  IN U.
  1426. C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
  1427. C                                  VECTORS IN U.
  1428. C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
  1429. C                                  VECTORS.
  1430. C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
  1431. C                                  IN V.
  1432. C
  1433. C     ON RETURN
  1434. C
  1435. C         S         REAL(MM), WHERE MM=MIN(N+1,P).
  1436. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
  1437. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
  1438. C                   ORDER OF MAGNITUDE.
  1439. C
  1440. C         E         REAL(P).
  1441. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
  1442. C                   DISCUSSION OF INFO FOR EXCEPTIONS.
  1443. C
  1444. C         U         REAL(LDU,K), WHERE LDU.GE.N.  IF JOBA.EQ.1 THEN
  1445. C                                   K.EQ.N, IF JOBA.GE.2 THEN
  1446. C                                   K.EQ.MIN(N,P).
  1447. C                   U CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  1448. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
  1449. C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
  1450. C                   IN THE SUBROUTINE CALL.
  1451. C
  1452. C         V         REAL(LDV,P), WHERE LDV.GE.P.
  1453. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  1454. C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
  1455. C                   THEN V MAY BE IDENTIFIED WITH X IN THE
  1456. C                   SUBROUTINE CALL.
  1457. C
  1458. C         INFO      INTEGER.
  1459. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
  1460. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
  1461. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
  1462. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
  1463. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
  1464. C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
  1465. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
  1466. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
  1467. C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
  1468. C                   VALUES OF X AND B ARE THE SAME.
  1469. C
  1470. C     LINPACK. THIS VERSION DATED 03/19/79 .
  1471. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  1472. C
  1473. C     ***** USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  1474. C
  1475. C     EXTERNAL SROT
  1476. C     BLAS SAXPY,SDOT,SSCAL,SSWAP,SNRM2,SROTG
  1477. C     FORTRAN ABS,AMAX1,MAX0,MIN0,MOD,SQRT
  1478. C
  1479. C     INTERNAL VARIABLES
  1480. C
  1481.       INTEGER I, ITER, J, JOBU, K, KASE, KK, L, LL, LLS, LM1, LP1, LS,
  1482.      * LU, M, MAXIT, MM, MM1, MP1, NCT, NCTP1, NCU, NRT, NRTP1
  1483.       REAL SDOT, T, R
  1484.       REAL B, C, CS, EL, EMM1, F, G, SNRM2, SCALE, SHIFT, SL, SM, SN,
  1485.      * SMM1, T1, TEST, ZTEST
  1486.       LOGICAL WANTU, WANTV
  1487. C
  1488. C
  1489. C     SET THE MAXIMUM NUMBER OF ITERATIONS.
  1490. C
  1491.       MAXIT = 30
  1492. C
  1493. C     DETERMINE WHAT IS TO BE COMPUTED.
  1494. C
  1495.       WANTU = .FALSE.
  1496.       WANTV = .FALSE.
  1497.       JOBU = MOD(JOB,100)/10
  1498.       NCU = N
  1499.       IF (JOBU.GT.1) NCU = MIN0(N,P)
  1500.       IF (JOBU.NE.0) WANTU = .TRUE.
  1501.       IF (MOD(JOB,10).NE.0) WANTV = .TRUE.
  1502. C
  1503. C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
  1504. C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
  1505. C
  1506.       INFO = 0
  1507.       NCT = MIN0(N-1,P)
  1508.       NRT = MAX0(0,MIN0(P-2,N))
  1509.       LU = MAX0(NCT,NRT)
  1510.       IF (LU.LT.1) GO TO 170
  1511.       DO 160 L=1,LU
  1512.         LP1 = L + 1
  1513.         IF (L.GT.NCT) GO TO 20
  1514. C
  1515. C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
  1516. C           PLACE THE L-TH DIAGONAL IN S(L).
  1517. C
  1518.         S(L) = SNRM2(N-L+1,X(L,L),1)
  1519.         IF (S(L).EQ.0.0E0) GO TO 10
  1520.         IF (X(L,L).NE.0.0E0) S(L) = SIGN(S(L),X(L,L))
  1521.         CALL SSCAL(N-L+1, 1.0E0/S(L), X(L,L), 1)
  1522.         X(L,L) = 1.0E0 + X(L,L)
  1523.    10   CONTINUE
  1524.         S(L) = -S(L)
  1525.    20   CONTINUE
  1526.         IF (P.LT.LP1) GO TO 50
  1527.         DO 40 J=LP1,P
  1528.           IF (L.GT.NCT) GO TO 30
  1529.           IF (S(L).EQ.0.0E0) GO TO 30
  1530. C
  1531. C              APPLY THE TRANSFORMATION.
  1532. C
  1533.           T = -SDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
  1534.           CALL SAXPY(N-L+1, T, X(L,L), 1, X(L,J), 1)
  1535.    30     CONTINUE
  1536. C
  1537. C           PLACE THE L-TH ROW OF X INTO  E FOR THE
  1538. C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
  1539. C
  1540.           E(J) = X(L,J)
  1541.    40   CONTINUE
  1542.    50   CONTINUE
  1543.         IF (.NOT.WANTU .OR. L.GT.NCT) GO TO 70
  1544. C
  1545. C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
  1546. C           MULTIPLICATION.
  1547. C
  1548.         DO 60 I=L,N
  1549.           U(I,L) = X(I,L)
  1550.    60   CONTINUE
  1551.    70   CONTINUE
  1552.         IF (L.GT.NRT) GO TO 150
  1553. C
  1554. C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
  1555. C           L-TH SUPER-DIAGONAL IN E(L).
  1556. C
  1557.         E(L) = SNRM2(P-L,E(LP1),1)
  1558.         IF (E(L).EQ.0.0E0) GO TO 80
  1559.         IF (E(LP1).NE.0.0E0) E(L) = SIGN(E(L),E(LP1))
  1560.         CALL SSCAL(P-L, 1.0E0/E(L), E(LP1), 1)
  1561.         E(LP1) = 1.0E0 + E(LP1)
  1562.    80   CONTINUE
  1563.         E(L) = -E(L)
  1564.         IF (LP1.GT.N .OR. E(L).EQ.0.0E0) GO TO 120
  1565. C
  1566. C              APPLY THE TRANSFORMATION.
  1567. C
  1568.         DO 90 I=LP1,N
  1569.           WORK(I) = 0.0E0
  1570.    90   CONTINUE
  1571.         DO 100 J=LP1,P
  1572.           CALL SAXPY(N-L, E(J), X(LP1,J), 1, WORK(LP1), 1)
  1573.   100   CONTINUE
  1574.         DO 110 J=LP1,P
  1575.           CALL SAXPY(N-L, -E(J)/E(LP1), WORK(LP1), 1, X(LP1,J), 1)
  1576.   110   CONTINUE
  1577.   120   CONTINUE
  1578.         IF (.NOT.WANTV) GO TO 140
  1579. C
  1580. C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
  1581. C              BACK MULTIPLICATION.
  1582. C
  1583.         DO 130 I=LP1,P
  1584.           V(I,L) = E(I)
  1585.   130   CONTINUE
  1586.   140   CONTINUE
  1587.   150   CONTINUE
  1588.   160 CONTINUE
  1589.   170 CONTINUE
  1590. C
  1591. C     SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M.
  1592. C
  1593.       M = MIN0(P,N+1)
  1594.       NCTP1 = NCT + 1
  1595.       NRTP1 = NRT + 1
  1596.       IF (NCT.LT.P) S(NCTP1) = X(NCTP1,NCTP1)
  1597.       IF (N.LT.M) S(M) = 0.0E0
  1598.       IF (NRTP1.LT.M) E(NRTP1) = X(NRTP1,M)
  1599.       E(M) = 0.0E0
  1600. C
  1601. C     IF REQUIRED, GENERATE U.
  1602. C
  1603.       IF (.NOT.WANTU) GO TO 300
  1604.       IF (NCU.LT.NCTP1) GO TO 200
  1605.       DO 190 J=NCTP1,NCU
  1606.         DO 180 I=1,N
  1607.           U(I,J) = 0.0E0
  1608.   180   CONTINUE
  1609.         U(J,J) = 1.0E0
  1610.   190 CONTINUE
  1611.   200 CONTINUE
  1612.       IF (NCT.LT.1) GO TO 290
  1613.       DO 280 LL=1,NCT
  1614.         L = NCT - LL + 1
  1615.         IF (S(L).EQ.0.0E0) GO TO 250
  1616.         LP1 = L + 1
  1617.         IF (NCU.LT.LP1) GO TO 220
  1618.         DO 210 J=LP1,NCU
  1619.           T = -SDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
  1620.           CALL SAXPY(N-L+1, T, U(L,L), 1, U(L,J), 1)
  1621.   210   CONTINUE
  1622.   220   CONTINUE
  1623.         CALL SSCAL(N-L+1, -1.0E0, U(L,L), 1)
  1624.         U(L,L) = 1.0E0 + U(L,L)
  1625.         LM1 = L - 1
  1626.         IF (LM1.LT.1) GO TO 240
  1627.         DO 230 I=1,LM1
  1628.           U(I,L) = 0.0E0
  1629.   230   CONTINUE
  1630.   240   CONTINUE
  1631.         GO TO 270
  1632.   250   CONTINUE
  1633.         DO 260 I=1,N
  1634.           U(I,L) = 0.0E0
  1635.   260   CONTINUE
  1636.         U(L,L) = 1.0E0
  1637.   270   CONTINUE
  1638.   280 CONTINUE
  1639.   290 CONTINUE
  1640.   300 CONTINUE
  1641. C
  1642. C     IF IT IS REQUIRED, GENERATE V.
  1643. C
  1644.       IF (.NOT.WANTV) GO TO 350
  1645.       DO 340 LL=1,P
  1646.         L = P - LL + 1
  1647.         LP1 = L + 1
  1648.         IF (L.GT.NRT) GO TO 320
  1649.         IF (E(L).EQ.0.0E0) GO TO 320
  1650.         DO 310 J=LP1,P
  1651.           T = -SDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
  1652.           CALL SAXPY(P-L, T, V(LP1,L), 1, V(LP1,J), 1)
  1653.   310   CONTINUE
  1654.   320   CONTINUE
  1655.         DO 330 I=1,P
  1656.           V(I,L) = 0.0E0
  1657.   330   CONTINUE
  1658.         V(L,L) = 1.0E0
  1659.   340 CONTINUE
  1660.   350 CONTINUE
  1661. C
  1662. C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
  1663. C
  1664.       MM = M
  1665.       ITER = 0
  1666.   360 CONTINUE
  1667. C
  1668. C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
  1669. C
  1670. C     ...EXIT
  1671.       IF (M.EQ.0) GO TO 620
  1672. C
  1673. C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
  1674. C        FLAG AND RETURN.
  1675. C
  1676.       IF (ITER.LT.MAXIT) GO TO 370
  1677.       INFO = M
  1678. C     ......EXIT
  1679.       GO TO 620
  1680.   370 CONTINUE
  1681. C
  1682. C        THIS SECTION OF THE PROGRAM INSPECTS FOR
  1683. C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
  1684. C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
  1685. C
  1686. C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
  1687. C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
  1688. C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
  1689. C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
  1690. C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
  1691. C
  1692.       DO 390 LL=1,M
  1693.         L = M - LL
  1694. C        ...EXIT
  1695.         IF (L.EQ.0) GO TO 400
  1696.         TEST = ABS(S(L)) + ABS(S(L+1))
  1697.         ZTEST = TEST + ABS(E(L))
  1698.         IF (ZTEST.NE.TEST) GO TO 380
  1699.         E(L) = 0.0E0
  1700. C        ......EXIT
  1701.         GO TO 400
  1702.   380   CONTINUE
  1703.   390 CONTINUE
  1704.   400 CONTINUE
  1705.       IF (L.NE.M-1) GO TO 410
  1706.       KASE = 4
  1707.       GO TO 480
  1708.   410 CONTINUE
  1709.       LP1 = L + 1
  1710.       MP1 = M + 1
  1711.       DO 430 LLS=LP1,MP1
  1712.         LS = M - LLS + LP1
  1713. C           ...EXIT
  1714.         IF (LS.EQ.L) GO TO 440
  1715.         TEST = 0.0E0
  1716.         IF (LS.NE.M) TEST = TEST + ABS(E(LS))
  1717.         IF (LS.NE.L+1) TEST = TEST + ABS(E(LS-1))
  1718.         ZTEST = TEST + ABS(S(LS))
  1719.         IF (ZTEST.NE.TEST) GO TO 420
  1720.         S(LS) = 0.0E0
  1721. C           ......EXIT
  1722.         GO TO 440
  1723.   420   CONTINUE
  1724.   430 CONTINUE
  1725.   440 CONTINUE
  1726.       IF (LS.NE.L) GO TO 450
  1727.       KASE = 3
  1728.       GO TO 470
  1729.   450 CONTINUE
  1730.       IF (LS.NE.M) GO TO 460
  1731.       KASE = 1
  1732.       GO TO 470
  1733.   460 CONTINUE
  1734.       KASE = 2
  1735.       L = LS
  1736.   470 CONTINUE
  1737.   480 CONTINUE
  1738.       L = L + 1
  1739. C
  1740. C        PERFORM THE TASK INDICATED BY KASE.
  1741. C
  1742.       GO TO (490, 520, 540, 570), KASE
  1743. C
  1744. C        DEFLATE NEGLIGIBLE S(M).
  1745. C
  1746.   490 CONTINUE
  1747.       MM1 = M - 1
  1748.       F = E(M-1)
  1749.       E(M-1) = 0.0E0
  1750.       DO 510 KK=L,MM1
  1751.         K = MM1 - KK + L
  1752.         T1 = S(K)
  1753.         CALL SROTG(T1, F, CS, SN)
  1754.         S(K) = T1
  1755.         IF (K.EQ.L) GO TO 500
  1756.         F = -SN*E(K-1)
  1757.         E(K-1) = CS*E(K-1)
  1758.   500   CONTINUE
  1759.         IF (WANTV) CALL SROT(P, V(1,K), 1, V(1,M), 1, CS, SN)
  1760.   510 CONTINUE
  1761.       GO TO 610
  1762. C
  1763. C        SPLIT AT NEGLIGIBLE S(L).
  1764. C
  1765.   520 CONTINUE
  1766.       F = E(L-1)
  1767.       E(L-1) = 0.0E0
  1768.       DO 530 K=L,M
  1769.         T1 = S(K)
  1770.         CALL SROTG(T1, F, CS, SN)
  1771.         S(K) = T1
  1772.         F = -SN*E(K)
  1773.         E(K) = CS*E(K)
  1774.         IF (WANTU) CALL SROT(N, U(1,K), 1, U(1,L-1), 1, CS, SN)
  1775.   530 CONTINUE
  1776.       GO TO 610
  1777. C
  1778. C        PERFORM ONE QR STEP.
  1779. C
  1780.   540 CONTINUE
  1781. C
  1782. C           CALCULATE THE SHIFT.
  1783. C
  1784.       SCALE = AMAX1(ABS(S(M)),ABS(S(M-1)),ABS(E(M-1)),ABS(S(L)),ABS(E(L)
  1785.      * ))
  1786.       SM = S(M)/SCALE
  1787.       SMM1 = S(M-1)/SCALE
  1788.       EMM1 = E(M-1)/SCALE
  1789.       SL = S(L)/SCALE
  1790.       EL = E(L)/SCALE
  1791.       B = ((SMM1+SM)*(SMM1-SM)+EMM1**2)/2.0E0
  1792.       C = (SM*EMM1)**2
  1793.       SHIFT = 0.0E0
  1794.       IF (B.EQ.0.0E0 .AND. C.EQ.0.0E0) GO TO 550
  1795.       SHIFT = SQRT(B**2+C)
  1796.       IF (B.LT.0.0E0) SHIFT = -SHIFT
  1797.       SHIFT = C/(B+SHIFT)
  1798.   550 CONTINUE
  1799.       F = (SL+SM)*(SL-SM) - SHIFT
  1800.       G = SL*EL
  1801. C
  1802. C           CHASE ZEROS.
  1803. C
  1804.       MM1 = M - 1
  1805.       DO 560 K=L,MM1
  1806.         CALL SROTG(F, G, CS, SN)
  1807.         IF (K.NE.L) E(K-1) = F
  1808.         F = CS*S(K) + SN*E(K)
  1809.         E(K) = CS*E(K) - SN*S(K)
  1810.         G = SN*S(K+1)
  1811.         S(K+1) = CS*S(K+1)
  1812.         IF (WANTV) CALL SROT(P, V(1,K), 1, V(1,K+1), 1, CS, SN)
  1813.         CALL SROTG(F, G, CS, SN)
  1814.         S(K) = F
  1815.         F = CS*E(K) + SN*S(K+1)
  1816.         S(K+1) = -SN*E(K) + CS*S(K+1)
  1817.         G = SN*E(K+1)
  1818.         E(K+1) = CS*E(K+1)
  1819.         IF (WANTU .AND. K.LT.N) CALL SROT(N, U(1,K), 1, U(1,K+1), 1,
  1820.      *   CS, SN)
  1821.   560 CONTINUE
  1822.       E(M-1) = F
  1823.       ITER = ITER + 1
  1824.       GO TO 610
  1825. C
  1826. C        CONVERGENCE.
  1827. C
  1828.   570 CONTINUE
  1829. C
  1830. C           MAKE THE SINGULAR VALUE  POSITIVE.
  1831. C
  1832.       IF (S(L).GE.0.0E0) GO TO 580
  1833.       S(L) = -S(L)
  1834.       IF (WANTV) CALL SSCAL(P, -1.0E0, V(1,L), 1)
  1835.   580 CONTINUE
  1836. C
  1837. C           ORDER THE SINGULAR VALUE.
  1838. C
  1839.   590 IF (L.EQ.MM) GO TO 600
  1840. C           ...EXIT
  1841.       IF (S(L).GE.S(L+1)) GO TO 600
  1842.       T = S(L)
  1843.       S(L) = S(L+1)
  1844.       S(L+1) = T
  1845.       IF (WANTV .AND. L.LT.P) CALL SSWAP(P, V(1,L), 1, V(1,L+1), 1)
  1846.       IF (WANTU .AND. L.LT.N) CALL SSWAP(N, U(1,L), 1, U(1,L+1), 1)
  1847.       L = L + 1
  1848.       GO TO 590
  1849.   600 CONTINUE
  1850.       ITER = 0
  1851.       M = M - 1
  1852.   610 CONTINUE
  1853.       GO TO 360
  1854.   620 CONTINUE
  1855.       RETURN
  1856.       END
  1857.       REAL FUNCTION SASUM(N, SX, INCX)
  1858. C
  1859. C     TAKES THE SUM OF THE ABSOLUTE VALUES.
  1860. C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
  1861. C     JACK DONGARRA, LINPACK, 3/11/78.
  1862. C
  1863.       REAL SX(1), STEMP
  1864.       INTEGER I, INCX, M, MP1, N, NINCX
  1865. C
  1866.       SASUM = 0.0E0
  1867.       STEMP = 0.0E0
  1868.       IF (N.LE.0) RETURN
  1869.       IF (INCX.EQ.1) GO TO 20
  1870. C
  1871. C        CODE FOR INCREMENT NOT EQUAL TO 1
  1872. C
  1873.       NINCX = N*INCX
  1874.       DO 10 I=1,NINCX,INCX
  1875.         STEMP = STEMP + ABS(SX(I))
  1876.    10 CONTINUE
  1877.       SASUM = STEMP
  1878.       RETURN
  1879. C
  1880. C        CODE FOR INCREMENT EQUAL TO 1
  1881. C
  1882. C
  1883. C        CLEAN-UP LOOP
  1884. C
  1885.    20 M = MOD(N,6)
  1886.       IF (M.EQ.0) GO TO 40
  1887.       DO 30 I=1,M
  1888.         STEMP = STEMP + ABS(SX(I))
  1889.    30 CONTINUE
  1890.       IF (N.LT.6) GO TO 60
  1891.    40 MP1 = M + 1
  1892.       DO 50 I=MP1,N,6
  1893.         STEMP = STEMP + ABS(SX(I)) + ABS(SX(I+1)) + ABS(SX(I+2)) +
  1894.      *   ABS(SX(I+3)) + ABS(SX(I+4)) + ABS(SX(I+5))
  1895.    50 CONTINUE
  1896.    60 SASUM = STEMP
  1897.       RETURN
  1898.       END
  1899.       SUBROUTINE SAXPY(N, SA, SX, INCX, SY, INCY)
  1900. C
  1901. C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
  1902. C     USES UNROLLED LOOP FOR INCREMENTS EQUAL TO ONE.
  1903. C     JACK DONGARRA, LINPACK, 3/11/78.
  1904. C
  1905.       REAL SX(1), SY(1), SA
  1906.       INTEGER I, INCX, INCY, IX, IY, M, MP1, N
  1907. C
  1908.       IF (N.LE.0) RETURN
  1909.       IF (SA.EQ.0.0) RETURN
  1910.       IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
  1911. C
  1912. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  1913. C          NOT EQUAL TO 1
  1914. C
  1915.       IX = 1
  1916.       IY = 1
  1917.       IF (INCX.LT.0) IX = (-N+1)*INCX + 1
  1918.       IF (INCY.LT.0) IY = (-N+1)*INCY + 1
  1919.       DO 10 I=1,N
  1920.         SY(IY) = SY(IY) + SA*SX(IX)
  1921.         IX = IX + INCX
  1922.         IY = IY + INCY
  1923.    10 CONTINUE
  1924.       RETURN
  1925. C
  1926. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  1927. C
  1928. C
  1929. C        CLEAN-UP LOOP
  1930. C
  1931.    20 M = MOD(N,4)
  1932.       IF (M.EQ.0) GO TO 40
  1933.       DO 30 I=1,M
  1934.         SY(I) = SY(I) + SA*SX(I)
  1935.    30 CONTINUE
  1936.       IF (N.LT.4) RETURN
  1937.    40 MP1 = M + 1
  1938.       DO 50 I=MP1,N,4
  1939.         SY(I) = SY(I) + SA*SX(I)
  1940.         SY(I+1) = SY(I+1) + SA*SX(I+1)
  1941.         SY(I+2) = SY(I+2) + SA*SX(I+2)
  1942.         SY(I+3) = SY(I+3) + SA*SX(I+3)
  1943.    50 CONTINUE
  1944.       RETURN
  1945.       END
  1946.       REAL FUNCTION SDOT(N, SX, INCX, SY, INCY)
  1947. C
  1948. C     FORMS THE DOT PRODUCT OF TWO VECTORS.
  1949. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  1950. C     JACK DONGARRA, LINPACK, 3/11/78.
  1951. C
  1952.       REAL SX(1), SY(1), STEMP
  1953.       INTEGER I, INCX, INCY, IX, IY, M, MP1, N
  1954. C
  1955.       STEMP = 0.0E0
  1956.       SDOT = 0.0E0
  1957.       IF (N.LE.0) RETURN
  1958.       IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
  1959. C
  1960. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  1961. C          NOT EQUAL TO 1
  1962. C
  1963.       IX = 1
  1964.       IY = 1
  1965.       IF (INCX.LT.0) IX = (-N+1)*INCX + 1
  1966.       IF (INCY.LT.0) IY = (-N+1)*INCY + 1
  1967.       DO 10 I=1,N
  1968.         STEMP = STEMP + SX(IX)*SY(IY)
  1969.         IX = IX + INCX
  1970.         IY = IY + INCY
  1971.    10 CONTINUE
  1972.       SDOT = STEMP
  1973.       RETURN
  1974. C
  1975. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  1976. C
  1977. C
  1978. C        CLEAN-UP LOOP
  1979. C
  1980.    20 M = MOD(N,5)
  1981.       IF (M.EQ.0) GO TO 40
  1982.       DO 30 I=1,M
  1983.         STEMP = STEMP + SX(I)*SY(I)
  1984.    30 CONTINUE
  1985.       IF (N.LT.5) GO TO 60
  1986.    40 MP1 = M + 1
  1987.       DO 50 I=MP1,N,5
  1988.         STEMP = STEMP + SX(I)*SY(I) + SX(I+1)*SY(I+1) + SX(I+2)*SY(I+2)
  1989.      *   + SX(I+3)*SY(I+3) + SX(I+4)*SY(I+4)
  1990.    50 CONTINUE
  1991.    60 SDOT = STEMP
  1992.       RETURN
  1993.       END
  1994.       REAL FUNCTION SNRM2(N, SX, INCX)
  1995.       INTEGER NEXT
  1996.       REAL SX(1), CUTLO, CUTHI, HITEST, SUM, XMAX, ZERO, ONE
  1997.       DATA ZERO, ONE /0.0E0,1.0E0/
  1998. C
  1999. C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN SX() WITH STORAGE
  2000. C     INCREMENT INCX .
  2001. C     IF    N .LE. 0 RETURN WITH RESULT = 0.
  2002. C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
  2003. C
  2004. C           C.L.LAWSON, 1978 JAN 08
  2005. C
  2006. C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
  2007. C     HOPEFULLY APPLICABLE TO ALL MACHINES.
  2008. C         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
  2009. C         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
  2010. C     WHERE
  2011. C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
  2012. C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
  2013. C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
  2014. C
  2015. C     BRIEF OUTLINE OF ALGORITHM..
  2016. C
  2017. C     PHASE 1    SCANS ZERO COMPONENTS.
  2018. C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
  2019. C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
  2020. C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
  2021. C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
  2022. C
  2023. C     VALUES FOR CUTLO AND CUTHI..
  2024. C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
  2025. C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
  2026. C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
  2027. C                   UNIVAC AND DEC AT 2**(-103)
  2028. C                   THUS CUTLO = 2**(-51) = 4.44089E-16
  2029. C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
  2030. C                   THUS CUTHI = 2**(63.5) = 1.30438E19
  2031. C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
  2032. C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
  2033. C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
  2034. C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
  2035. C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
  2036.       DATA CUTLO, CUTHI /4.441E-16,1.304E19/
  2037. C
  2038.       IF (N.GT.0) GO TO 10
  2039.       SNRM2 = ZERO
  2040.       GO TO 140
  2041. C
  2042.    10 ASSIGN 30 TO NEXT
  2043.       SUM = ZERO
  2044.       NN = N*INCX
  2045. C                                                 BEGIN MAIN LOOP
  2046.       I = 1
  2047.    20 GO TO NEXT, (30, 40, 70, 80)
  2048.    30 IF (ABS(SX(I)).GT.CUTLO) GO TO 110
  2049.       ASSIGN 40 TO NEXT
  2050.       XMAX = ZERO
  2051. C
  2052. C                        PHASE 1.  SUM IS ZERO
  2053. C
  2054.    40 IF (SX(I).EQ.ZERO) GO TO 130
  2055.       IF (ABS(SX(I)).GT.CUTLO) GO TO 110
  2056. C
  2057. C                                PREPARE FOR PHASE 2.
  2058.       ASSIGN 70 TO NEXT
  2059.       GO TO 60
  2060. C
  2061. C                                PREPARE FOR PHASE 4.
  2062. C
  2063.    50 I = J
  2064.       ASSIGN 80 TO NEXT
  2065.       SUM = (SUM/SX(I))/SX(I)
  2066.    60 XMAX = ABS(SX(I))
  2067.       GO TO 90
  2068. C
  2069. C                   PHASE 2.  SUM IS SMALL.
  2070. C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
  2071. C
  2072.    70 IF (ABS(SX(I)).GT.CUTLO) GO TO 100
  2073. C
  2074. C                     COMMON CODE FOR PHASES 2 AND 4.
  2075. C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
  2076. C
  2077.    80 IF (ABS(SX(I)).LE.XMAX) GO TO 90
  2078.       SUM = ONE + SUM*(XMAX/SX(I))**2
  2079.       XMAX = ABS(SX(I))
  2080.       GO TO 130
  2081. C
  2082.    90 SUM = SUM + (SX(I)/XMAX)**2
  2083.       GO TO 130
  2084. C
  2085. C
  2086. C                  PREPARE FOR PHASE 3.
  2087. C
  2088.   100 SUM = (SUM*XMAX)*XMAX
  2089. C
  2090. C
  2091. C     FOR REAL OR D.P. SET HITEST = CUTHI/N
  2092. C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
  2093. C
  2094.   110 HITEST = CUTHI/FLOAT(N)
  2095. C
  2096. C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
  2097. C
  2098.       DO 120 J=I,NN,INCX
  2099.         IF (ABS(SX(J)).GE.HITEST) GO TO 50
  2100.         SUM = SUM + SX(J)**2
  2101.   120 CONTINUE
  2102.       SNRM2 = SQRT(SUM)
  2103.       GO TO 140
  2104. C
  2105.   130 CONTINUE
  2106.       I = I + INCX
  2107.       IF (I.LE.NN) GO TO 20
  2108. C
  2109. C              END OF MAIN LOOP.
  2110. C
  2111. C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
  2112. C
  2113.       SNRM2 = XMAX*SQRT(SUM)
  2114.   140 CONTINUE
  2115.       RETURN
  2116.       END
  2117.       SUBROUTINE SROT(N, SX, INCX, SY, INCY, C, S)
  2118. C
  2119. C     APPLIES A PLANE ROTATION.
  2120. C     JACK DONGARRA, LINPACK, 3/11/78.
  2121. C
  2122.       REAL SX(1), SY(1), STEMP, C, S
  2123.       INTEGER I, INCX, INCY, IX, IY, N
  2124. C
  2125.       IF (N.LE.0) RETURN
  2126.       IF (INCX.EQ.1 .AND. INCY.EQ.1) GO TO 20
  2127. C
  2128. C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
  2129. C         TO 1
  2130. C
  2131.       IX = 1
  2132.       IY = 1
  2133.       IF (INCX.LT.0) IX = (-N+1)*INCX + 1
  2134.       IF (INCY.LT.0) IY = (-N+1)*INCY + 1
  2135.       DO 10 I=1,N
  2136.         STEMP = C*SX(IX) + S*SY(IY)
  2137.         SY(IY) = C*SY(IY) - S*SX(IX)
  2138.         SX(IX) = STEMP
  2139.         IX = IX + INCX
  2140.         IY = IY + INCY
  2141.    10 CONTINUE
  2142.       RETURN
  2143. C
  2144. C       CODE FOR BOTH INCREMENTS EQUAL TO 1
  2145. C
  2146.    20 DO 30 I=1,N
  2147.         STEMP = C*SX(I) + S*SY(I)
  2148.         SY(I) = C*SY(I) - S*SX(I)
  2149.         SX(I) = STEMP
  2150.    30 CONTINUE
  2151.       RETURN
  2152.       END
  2153.       SUBROUTINE SROTG(SA, SB, C, S)
  2154. C
  2155. C     CONSTRUCT GIVENS PLANE ROTATION.
  2156. C     JACK DONGARRA, LINPACK, 3/11/78.
  2157. C
  2158.       REAL SA, SB, C, S, ROE, SCALE, R, Z
  2159. C
  2160.       ROE = SB
  2161.       IF (ABS(SA).GT.ABS(SB)) ROE = SA
  2162.       SCALE = ABS(SA) + ABS(SB)
  2163.       IF (SCALE.NE.0.0) GO TO 10
  2164.       C = 1.0
  2165.       S = 0.0
  2166.       R = 0.0
  2167.       GO TO 20
  2168.    10 R = SCALE*SQRT((SA/SCALE)**2+(SB/SCALE)**2)
  2169.       R = SIGN(1.0,ROE)*R
  2170.       C = SA/R
  2171.       S = SB/R
  2172.    20 Z = 1.0
  2173.       IF (ABS(SA).GT.ABS(SB)) Z = S
  2174.       IF (ABS(SB).GE.ABS(SA) .AND. C.NE.0.0) Z = 1.0/C
  2175.       SA = R
  2176.       SB = Z
  2177.       RETURN
  2178.       END
  2179.       SUBROUTINE SSCAL(N, SA, SX, INCX)
  2180. C
  2181. C     SCALES A VECTOR BY A CONSTANT.
  2182. C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO 1.
  2183. C     JACK DONGARRA, LINPACK, 3/11/78.
  2184. C
  2185.       REAL SA, SX(1)
  2186.       INTEGER I, INCX, M, MP1, N, NINCX
  2187. C
  2188.       IF (N.LE.0) RETURN
  2189.       IF (INCX.EQ.1) GO TO 20
  2190. C
  2191. C        CODE FOR INCREMENT NOT EQUAL TO 1
  2192. C
  2193.       NINCX = N*INCX
  2194.       DO 10 I=1,NINCX,INCX
  2195.         SX(I) = SA*SX(I)
  2196.    10 CONTINUE
  2197.       RETURN
  2198. C
  2199. C        CODE FOR INCREMENT EQUAL TO 1
  2200. C
  2201. C
  2202. C        CLEAN-UP LOOP
  2203. C
  2204.    20 M = MOD(N,5)
  2205.       IF (M.EQ.0) GO TO 40
  2206.       DO 30 I=1,M
  2207.         SX(I) = SA*SX(I)
  2208.    30 CONTINUE
  2209.       IF (N.LT.5) RETURN
  2210.    40 MP1 = M + 1
  2211.       DO 50 I=MP1,N,5
  2212.         SX(I) = SA*SX(I)
  2213.         SX(I+1) = SA*SX(I+1)
  2214.         SX(I+2) = SA*SX(I+2)
  2215.         SX(I+3) = SA*SX(I+3)
  2216.         SX(I+4) = SA*SX(I+4)
  2217.    50 CONTINUE
  2218.       RETURN
  2219.       END
  2220.       SUBROUTINE SVD(NM, M, N, A, W, MATU, U, MATV, V, IERR, RV1)
  2221. C
  2222.       INTEGER I, J, K, L, M, N, II, I1, KK, K1, LL, L1, MN, NM, ITS,
  2223.      * IERR
  2224.       REAL A(NM,N), W(N), U(NM,N), V(NM,N), RV1(N)
  2225.       REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, MACHEP
  2226.       REAL SQRT, AMAX1, ABS, SIGN
  2227.       LOGICAL MATU, MATV
  2228. C
  2229. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE SVD,
  2230. C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
  2231. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
  2232. C
  2233. C     THIS SUBROUTINE DETERMINES THE SINGULAR VALUE DECOMPOSITION
  2234. C          T
  2235. C     A=USV  OF A REAL M BY N RECTANGULAR MATRIX.  HOUSEHOLDER
  2236. C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
  2237. C
  2238. C     ON INPUT-
  2239. C
  2240. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2241. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2242. C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
  2243. C          AS LARGE AS THE MAXIMUM OF M AND N,
  2244. C
  2245. C        M IS THE NUMBER OF ROWS OF A (AND U),
  2246. C
  2247. C        N IS THE NUMBER OF COLUMNS OF A (AND U) AND THE ORDER OF V,
  2248. C
  2249. C        A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED,
  2250. C
  2251. C        MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE
  2252. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE,
  2253. C
  2254. C        MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE
  2255. C          DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE.
  2256. C
  2257. C     ON OUTPUT-
  2258. C
  2259. C        A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V),
  2260. C
  2261. C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  2262. C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
  2263. C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
  2264. C          FOR INDICES IERR+1,IERR+2,...,N,
  2265. C
  2266. C        U CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE
  2267. C          DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE.  OTHERWISE
  2268. C          U IS USED AS A TEMPORARY ARRAY.  U MAY COINCIDE WITH A.
  2269. C          IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING
  2270. C          TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
  2271. C
  2272. C        V CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF
  2273. C          MATV HAS BEEN SET TO .TRUE.  OTHERWISE V IS NOT REFERENCED.
  2274. C          V MAY ALSO COINCIDE WITH A IF U IS NOT NEEDED.  IF AN ERROR
  2275. C          EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF
  2276. C          CORRECT SINGULAR VALUES SHOULD BE CORRECT,
  2277. C
  2278. C        IERR IS SET TO
  2279. C          ZERO       FOR NORMAL RETURN,
  2280. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  2281. C                     DETERMINED AFTER 30 ITERATIONS,
  2282. C
  2283. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  2284. C
  2285. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2286. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2287. C
  2288. C     ------------------------------------------------------------------
  2289. C
  2290. C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
  2291. C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
  2292. C
  2293. C                **********
  2294.       MACHEP = 2.**(-26)
  2295. C
  2296.       IERR = 0
  2297. C
  2298.       DO 20 I=1,M
  2299. C
  2300.         DO 10 J=1,N
  2301.           U(I,J) = A(I,J)
  2302.    10   CONTINUE
  2303.    20 CONTINUE
  2304. C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
  2305.       G = 0.0
  2306.       SCALE = 0.0
  2307.       X = 0.0
  2308. C
  2309.       DO 200 I=1,N
  2310.         L = I + 1
  2311.         RV1(I) = SCALE*G
  2312.         G = 0.0
  2313.         S = 0.0
  2314.         SCALE = 0.0
  2315.         IF (I.GT.M) GO TO 100
  2316. C
  2317.         DO 30 K=I,M
  2318.           SCALE = SCALE + ABS(U(K,I))
  2319.    30   CONTINUE
  2320. C
  2321.         IF (SCALE.EQ.0.0) GO TO 100
  2322. C
  2323.         DO 40 K=I,M
  2324.           U(K,I) = U(K,I)/SCALE
  2325.           S = S + U(K,I)**2
  2326.    40   CONTINUE
  2327. C
  2328.         F = U(I,I)
  2329.         G = -SIGN(SQRT(S),F)
  2330.         H = F*G - S
  2331.         U(I,I) = F - G
  2332.         IF (I.EQ.N) GO TO 80
  2333. C
  2334.         DO 70 J=L,N
  2335.           S = 0.0
  2336. C
  2337.           DO 50 K=I,M
  2338.             S = S + U(K,I)*U(K,J)
  2339.    50     CONTINUE
  2340. C
  2341.           F = S/H
  2342. C
  2343.           DO 60 K=I,M
  2344.             U(K,J) = U(K,J) + F*U(K,I)
  2345.    60     CONTINUE
  2346.    70   CONTINUE
  2347. C
  2348.    80   DO 90 K=I,M
  2349.           U(K,I) = SCALE*U(K,I)
  2350.    90   CONTINUE
  2351. C
  2352.   100   W(I) = SCALE*G
  2353.         G = 0.0
  2354.         S = 0.0
  2355.         SCALE = 0.0
  2356.         IF (I.GT.M .OR. I.EQ.N) GO TO 190
  2357. C
  2358.         DO 110 K=L,N
  2359.           SCALE = SCALE + ABS(U(I,K))
  2360.   110   CONTINUE
  2361. C
  2362.         IF (SCALE.EQ.0.0) GO TO 190
  2363. C
  2364.         DO 120 K=L,N
  2365.           U(I,K) = U(I,K)/SCALE
  2366.           S = S + U(I,K)**2
  2367.   120   CONTINUE
  2368. C
  2369.         F = U(I,L)
  2370.         G = -SIGN(SQRT(S),F)
  2371.         H = F*G - S
  2372.         U(I,L) = F - G
  2373. C
  2374.         DO 130 K=L,N
  2375.           RV1(K) = U(I,K)/H
  2376.   130   CONTINUE
  2377. C
  2378.         IF (I.EQ.M) GO TO 170
  2379. C
  2380.         DO 160 J=L,M
  2381.           S = 0.0
  2382. C
  2383.           DO 140 K=L,N
  2384.             S = S + U(J,K)*U(I,K)
  2385.   140     CONTINUE
  2386. C
  2387.           DO 150 K=L,N
  2388.             U(J,K) = U(J,K) + S*RV1(K)
  2389.   150     CONTINUE
  2390.   160   CONTINUE
  2391. C
  2392.   170   DO 180 K=L,N
  2393.           U(I,K) = SCALE*U(I,K)
  2394.   180   CONTINUE
  2395. C
  2396.   190   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  2397.   200 CONTINUE
  2398. C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS **********
  2399.       IF (.NOT.MATV) GO TO 290
  2400. C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
  2401.       DO 280 II=1,N
  2402.         I = N + 1 - II
  2403.         IF (I.EQ.N) GO TO 270
  2404.         IF (G.EQ.0.0) GO TO 250
  2405. C
  2406.         DO 210 J=L,N
  2407. C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
  2408.           V(J,I) = (U(I,J)/U(I,L))/G
  2409.   210   CONTINUE
  2410. C
  2411.         DO 240 J=L,N
  2412.           S = 0.0
  2413. C
  2414.           DO 220 K=L,N
  2415.             S = S + U(I,K)*V(K,J)
  2416.   220     CONTINUE
  2417. C
  2418.           DO 230 K=L,N
  2419.             V(K,J) = V(K,J) + S*V(K,I)
  2420.   230     CONTINUE
  2421.   240   CONTINUE
  2422. C
  2423.   250   DO 260 J=L,N
  2424.           V(I,J) = 0.0
  2425.           V(J,I) = 0.0
  2426.   260   CONTINUE
  2427. C
  2428.   270   V(I,I) = 1.0
  2429.         G = RV1(I)
  2430.         L = I
  2431.   280 CONTINUE
  2432. C     ********** ACCUMULATION OF LEFT-HAND TRANSFORMATIONS **********
  2433.   290 IF (.NOT.MATU) GO TO 410
  2434. C     **********FOR I=MIN(M,N) STEP -1 UNTIL 1 DO -- **********
  2435.       MN = N
  2436.       IF (M.LT.N) MN = M
  2437. C
  2438.       DO 400 II=1,MN
  2439.         I = MN + 1 - II
  2440.         L = I + 1
  2441.         G = W(I)
  2442.         IF (I.EQ.N) GO TO 310
  2443. C
  2444.         DO 300 J=L,N
  2445.           U(I,J) = 0.0
  2446.   300   CONTINUE
  2447. C
  2448.   310   IF (G.EQ.0.0) GO TO 370
  2449.         IF (I.EQ.MN) GO TO 350
  2450. C
  2451.         DO 340 J=L,N
  2452.           S = 0.0
  2453. C
  2454.           DO 320 K=L,M
  2455.             S = S + U(K,I)*U(K,J)
  2456.   320     CONTINUE
  2457. C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
  2458.           F = (S/U(I,I))/G
  2459. C
  2460.           DO 330 K=I,M
  2461.             U(K,J) = U(K,J) + F*U(K,I)
  2462.   330     CONTINUE
  2463.   340   CONTINUE
  2464. C
  2465.   350   DO 360 J=I,M
  2466.           U(J,I) = U(J,I)/G
  2467.   360   CONTINUE
  2468. C
  2469.         GO TO 390
  2470. C
  2471.   370   DO 380 J=I,M
  2472.           U(J,I) = 0.0
  2473.   380   CONTINUE
  2474. C
  2475.   390   U(I,I) = U(I,I) + 1.0
  2476.   400 CONTINUE
  2477. C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  2478.   410 EPS = MACHEP*X
  2479. C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
  2480.       DO 550 KK=1,N
  2481.         K1 = N - KK
  2482.         K = K1 + 1
  2483.         ITS = 0
  2484. C     ********** TEST FOR SPLITTING.
  2485. C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  2486.   420   DO 430 LL=1,K
  2487.           L1 = K - LL
  2488.           L = L1 + 1
  2489.           IF (ABS(RV1(L)).LE.EPS) GO TO 470
  2490. C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
  2491. C                THROUGH THE BOTTOM OF THE LOOP **********
  2492.           IF (ABS(W(L1)).LE.EPS) GO TO 440
  2493.   430   CONTINUE
  2494. C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  2495.   440   C = 0.0
  2496.         S = 1.0
  2497. C
  2498.         DO 460 I=L,K
  2499.           F = S*RV1(I)
  2500.           RV1(I) = C*RV1(I)
  2501.           IF (ABS(F).LE.EPS) GO TO 470
  2502.           G = W(I)
  2503.           H = SQRT(F*F+G*G)
  2504.           W(I) = H
  2505.           C = G/H
  2506.           S = -F/H
  2507.           IF (.NOT.MATU) GO TO 460
  2508. C
  2509.           DO 450 J=1,M
  2510.             Y = U(J,L1)
  2511.             Z = U(J,I)
  2512.             U(J,L1) = Y*C + Z*S
  2513.             U(J,I) = -Y*S + Z*C
  2514.   450     CONTINUE
  2515. C
  2516.   460   CONTINUE
  2517. C     ********** TEST FOR CONVERGENCE **********
  2518.   470   Z = W(K)
  2519.         IF (L.EQ.K) GO TO 530
  2520. C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
  2521.         IF (ITS.EQ.30) GO TO 560
  2522.         ITS = ITS + 1
  2523.         X = W(L)
  2524.         Y = W(K1)
  2525.         G = RV1(K1)
  2526.         H = RV1(K)
  2527.         F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
  2528.         G = SQRT(F*F+1.0)
  2529.         F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
  2530. C     ********** NEXT QR TRANSFORMATION **********
  2531.         C = 1.0
  2532.         S = 1.0
  2533. C
  2534.         DO 520 I1=L,K1
  2535.           I = I1 + 1
  2536.           G = RV1(I)
  2537.           Y = W(I)
  2538.           H = S*G
  2539.           G = C*G
  2540.           Z = SQRT(F*F+H*H)
  2541.           RV1(I1) = Z
  2542.           C = F/Z
  2543.           S = H/Z
  2544.           F = X*C + G*S
  2545.           G = -X*S + G*C
  2546.           H = Y*S
  2547.           Y = Y*C
  2548.           IF (.NOT.MATV) GO TO 490
  2549. C
  2550.           DO 480 J=1,N
  2551.             X = V(J,I1)
  2552.             Z = V(J,I)
  2553.             V(J,I1) = X*C + Z*S
  2554.             V(J,I) = -X*S + Z*C
  2555.   480     CONTINUE
  2556. C
  2557.   490     Z = SQRT(F*F+H*H)
  2558.           W(I1) = Z
  2559. C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
  2560.           IF (Z.EQ.0.0) GO TO 500
  2561.           C = F/Z
  2562.           S = H/Z
  2563.   500     F = C*G + S*Y
  2564.           X = -S*G + C*Y
  2565.           IF (.NOT.MATU) GO TO 520
  2566. C
  2567.           DO 510 J=1,M
  2568.             Y = U(J,I1)
  2569.             Z = U(J,I)
  2570.             U(J,I1) = Y*C + Z*S
  2571.             U(J,I) = -Y*S + Z*C
  2572.   510     CONTINUE
  2573. C
  2574.   520   CONTINUE
  2575. C
  2576.         RV1(L) = 0.0
  2577.         RV1(K) = F
  2578.         W(K) = X
  2579.         GO TO 420
  2580. C     ********** CONVERGENCE **********
  2581.   530   IF (Z.GE.0.0) GO TO 550
  2582. C     ********** W(K) IS MADE NON-NEGATIVE **********
  2583.         W(K) = -Z
  2584.         IF (.NOT.MATV) GO TO 550
  2585. C
  2586.         DO 540 J=1,N
  2587.           V(J,K) = -V(J,K)
  2588.   540   CONTINUE
  2589. C
  2590.   550 CONTINUE
  2591. C
  2592.       GO TO 570
  2593. C     ********** SET ERROR -- NO CONVERGENCE TO A
  2594. C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  2595.   560 IERR = K
  2596.   570 RETURN
  2597. C     ********** LAST CARD OF SVD **********
  2598.       END
  2599. C
  2600. C     ------------------------------------------------------------------
  2601. C
  2602.       SUBROUTINE MINFIT(NM, M, N, A, W, IP, B, IERR, RV1)
  2603. C
  2604.       INTEGER I, J, K, L, M, N, II, IP, I1, KK, K1, LL, L1, M1, NM,
  2605.      * ITS, IERR
  2606.       REAL A(NM,N), W(N), B(NM,IP), RV1(N)
  2607.       REAL C, F, G, H, S, X, Y, Z, EPS, SCALE, MACHEP
  2608.       REAL SQRT, AMAX1, ABS, SIGN
  2609. C
  2610. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT,
  2611. C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.
  2612. C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).
  2613. C
  2614. C     THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR
  2615. C                                                        T
  2616. C     SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV  OF A REAL
  2617. C                                         T
  2618. C     M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U.  HOUSEHOLDER
  2619. C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.
  2620. C
  2621. C     ON INPUT-
  2622. C
  2623. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2624. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2625. C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST
  2626. C          AS LARGE AS THE MAXIMUM OF M AND N,
  2627. C
  2628. C        M IS THE NUMBER OF ROWS OF A AND B,
  2629. C
  2630. C        N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V,
  2631. C
  2632. C        A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM,
  2633. C
  2634. C        IP IS THE NUMBER OF COLUMNS OF B.  IP CAN BE ZERO,
  2635. C
  2636. C        B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM
  2637. C          IF IP IS NOT ZERO.  OTHERWISE B IS NOT REFERENCED.
  2638. C
  2639. C     ON OUTPUT-
  2640. C
  2641. C        A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE
  2642. C          DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS.  IF AN
  2643. C          ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO
  2644. C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT,
  2645. C
  2646. C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE
  2647. C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN
  2648. C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT
  2649. C          FOR INDICES IERR+1,IERR+2,...,N,
  2650. C
  2651. C                                   T
  2652. C        B HAS BEEN OVERWRITTEN BY U B.  IF AN ERROR EXIT IS MADE,
  2653. C                       T
  2654. C          THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT
  2655. C          SINGULAR VALUES SHOULD BE CORRECT,
  2656. C
  2657. C        IERR IS SET TO
  2658. C          ZERO       FOR NORMAL RETURN,
  2659. C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN
  2660. C                     DETERMINED AFTER 30 ITERATIONS,
  2661. C
  2662. C        RV1 IS A TEMPORARY STORAGE ARRAY.
  2663. C
  2664. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2665. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2666. C
  2667. C     ------------------------------------------------------------------
  2668. C
  2669. C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
  2670. C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
  2671. C
  2672. C                **********
  2673.       MACHEP = 2.**(-26)
  2674. C
  2675.       IERR = 0
  2676. C     ********** HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM **********
  2677.       G = 0.0
  2678.       SCALE = 0.0
  2679.       X = 0.0
  2680. C
  2681.       DO 220 I=1,N
  2682.         L = I + 1
  2683.         RV1(I) = SCALE*G
  2684.         G = 0.0
  2685.         S = 0.0
  2686.         SCALE = 0.0
  2687.         IF (I.GT.M) GO TO 120
  2688. C
  2689.         DO 10 K=I,M
  2690.           SCALE = SCALE + ABS(A(K,I))
  2691.    10   CONTINUE
  2692. C
  2693.         IF (SCALE.EQ.0.0) GO TO 120
  2694. C
  2695.         DO 20 K=I,M
  2696.           A(K,I) = A(K,I)/SCALE
  2697.           S = S + A(K,I)**2
  2698.    20   CONTINUE
  2699. C
  2700.         F = A(I,I)
  2701.         G = -SIGN(SQRT(S),F)
  2702.         H = F*G - S
  2703.         A(I,I) = F - G
  2704.         IF (I.EQ.N) GO TO 60
  2705. C
  2706.         DO 50 J=L,N
  2707.           S = 0.0
  2708. C
  2709.           DO 30 K=I,M
  2710.             S = S + A(K,I)*A(K,J)
  2711.    30     CONTINUE
  2712. C
  2713.           F = S/H
  2714. C
  2715.           DO 40 K=I,M
  2716.             A(K,J) = A(K,J) + F*A(K,I)
  2717.    40     CONTINUE
  2718.    50   CONTINUE
  2719. C
  2720.    60   IF (IP.EQ.0) GO TO 100
  2721. C
  2722.         DO 90 J=1,IP
  2723.           S = 0.0
  2724. C
  2725.           DO 70 K=I,M
  2726.             S = S + A(K,I)*B(K,J)
  2727.    70     CONTINUE
  2728. C
  2729.           F = S/H
  2730. C
  2731.           DO 80 K=I,M
  2732.             B(K,J) = B(K,J) + F*A(K,I)
  2733.    80     CONTINUE
  2734.    90   CONTINUE
  2735. C
  2736.   100   DO 110 K=I,M
  2737.           A(K,I) = SCALE*A(K,I)
  2738.   110   CONTINUE
  2739. C
  2740.   120   W(I) = SCALE*G
  2741.         G = 0.0
  2742.         S = 0.0
  2743.         SCALE = 0.0
  2744.         IF (I.GT.M .OR. I.EQ.N) GO TO 210
  2745. C
  2746.         DO 130 K=L,N
  2747.           SCALE = SCALE + ABS(A(I,K))
  2748.   130   CONTINUE
  2749. C
  2750.         IF (SCALE.EQ.0.0) GO TO 210
  2751. C
  2752.         DO 140 K=L,N
  2753.           A(I,K) = A(I,K)/SCALE
  2754.           S = S + A(I,K)**2
  2755.   140   CONTINUE
  2756. C
  2757.         F = A(I,L)
  2758.         G = -SIGN(SQRT(S),F)
  2759.         H = F*G - S
  2760.         A(I,L) = F - G
  2761. C
  2762.         DO 150 K=L,N
  2763.           RV1(K) = A(I,K)/H
  2764.   150   CONTINUE
  2765. C
  2766.         IF (I.EQ.M) GO TO 190
  2767. C
  2768.         DO 180 J=L,M
  2769.           S = 0.0
  2770. C
  2771.           DO 160 K=L,N
  2772.             S = S + A(J,K)*A(I,K)
  2773.   160     CONTINUE
  2774. C
  2775.           DO 170 K=L,N
  2776.             A(J,K) = A(J,K) + S*RV1(K)
  2777.   170     CONTINUE
  2778.   180   CONTINUE
  2779. C
  2780.   190   DO 200 K=L,N
  2781.           A(I,K) = SCALE*A(I,K)
  2782.   200   CONTINUE
  2783. C
  2784.   210   X = AMAX1(X,ABS(W(I))+ABS(RV1(I)))
  2785.   220 CONTINUE
  2786. C     ********** ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS.
  2787. C                FOR I=N STEP -1 UNTIL 1 DO -- **********
  2788.       DO 300 II=1,N
  2789.         I = N + 1 - II
  2790.         IF (I.EQ.N) GO TO 290
  2791.         IF (G.EQ.0.0) GO TO 270
  2792. C
  2793.         DO 230 J=L,N
  2794. C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
  2795.           A(J,I) = (A(I,J)/A(I,L))/G
  2796.   230   CONTINUE
  2797. C
  2798.         DO 260 J=L,N
  2799.           S = 0.0
  2800. C
  2801.           DO 240 K=L,N
  2802.             S = S + A(I,K)*A(K,J)
  2803.   240     CONTINUE
  2804. C
  2805.           DO 250 K=L,N
  2806.             A(K,J) = A(K,J) + S*A(K,I)
  2807.   250     CONTINUE
  2808.   260   CONTINUE
  2809. C
  2810.   270   DO 280 J=L,N
  2811.           A(I,J) = 0.0
  2812.           A(J,I) = 0.0
  2813.   280   CONTINUE
  2814. C
  2815.   290   A(I,I) = 1.0
  2816.         G = RV1(I)
  2817.         L = I
  2818.   300 CONTINUE
  2819. C
  2820.       IF (M.GE.N .OR. IP.EQ.0) GO TO 330
  2821.       M1 = M + 1
  2822. C
  2823.       DO 320 I=M1,N
  2824. C
  2825.         DO 310 J=1,IP
  2826.           B(I,J) = 0.0
  2827.   310   CONTINUE
  2828.   320 CONTINUE
  2829. C     ********** DIAGONALIZATION OF THE BIDIAGONAL FORM **********
  2830.   330 EPS = MACHEP*X
  2831. C     ********** FOR K=N STEP -1 UNTIL 1 DO -- **********
  2832.       DO 460 KK=1,N
  2833.         K1 = N - KK
  2834.         K = K1 + 1
  2835.         ITS = 0
  2836. C     ********** TEST FOR SPLITTING.
  2837. C                FOR L=K STEP -1 UNTIL 1 DO -- **********
  2838.   340   DO 350 LL=1,K
  2839.           L1 = K - LL
  2840.           L = L1 + 1
  2841.           IF (ABS(RV1(L)).LE.EPS) GO TO 390
  2842. C     ********** RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT
  2843. C                THROUGH THE BOTTOM OF THE LOOP **********
  2844.           IF (ABS(W(L1)).LE.EPS) GO TO 360
  2845.   350   CONTINUE
  2846. C     ********** CANCELLATION OF RV1(L) IF L GREATER THAN 1 **********
  2847.   360   C = 0.0
  2848.         S = 1.0
  2849. C
  2850.         DO 380 I=L,K
  2851.           F = S*RV1(I)
  2852.           RV1(I) = C*RV1(I)
  2853.           IF (ABS(F).LE.EPS) GO TO 390
  2854.           G = W(I)
  2855.           H = SQRT(F*F+G*G)
  2856.           W(I) = H
  2857.           C = G/H
  2858.           S = -F/H
  2859.           IF (IP.EQ.0) GO TO 380
  2860. C
  2861.           DO 370 J=1,IP
  2862.             Y = B(L1,J)
  2863.             Z = B(I,J)
  2864.             B(L1,J) = Y*C + Z*S
  2865.             B(I,J) = -Y*S + Z*C
  2866.   370     CONTINUE
  2867. C
  2868.   380   CONTINUE
  2869. C     ********** TEST FOR CONVERGENCE **********
  2870.   390   Z = W(K)
  2871.         IF (L.EQ.K) GO TO 440
  2872. C     ********** SHIFT FROM BOTTOM 2 BY 2 MINOR **********
  2873.         IF (ITS.EQ.30) GO TO 470
  2874.         ITS = ITS + 1
  2875.         X = W(L)
  2876.         Y = W(K1)
  2877.         G = RV1(K1)
  2878.         H = RV1(K)
  2879.         F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)
  2880.         G = SQRT(F*F+1.0)
  2881.         F = ((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X
  2882. C     ********** NEXT QR TRANSFORMATION **********
  2883.         C = 1.0
  2884.         S = 1.0
  2885. C
  2886.         DO 430 I1=L,K1
  2887.           I = I1 + 1
  2888.           G = RV1(I)
  2889.           Y = W(I)
  2890.           H = S*G
  2891.           G = C*G
  2892.           Z = SQRT(F*F+H*H)
  2893.           RV1(I1) = Z
  2894.           C = F/Z
  2895.           S = H/Z
  2896.           F = X*C + G*S
  2897.           G = -X*S + G*C
  2898.           H = Y*S
  2899.           Y = Y*C
  2900. C
  2901.           DO 400 J=1,N
  2902.             X = A(J,I1)
  2903.             Z = A(J,I)
  2904.             A(J,I1) = X*C + Z*S
  2905.             A(J,I) = -X*S + Z*C
  2906.   400     CONTINUE
  2907. C
  2908.           Z = SQRT(F*F+H*H)
  2909.           W(I1) = Z
  2910. C     ********** ROTATION CAN BE ARBITRARY IF Z IS ZERO **********
  2911.           IF (Z.EQ.0.0) GO TO 410
  2912.           C = F/Z
  2913.           S = H/Z
  2914.   410     F = C*G + S*Y
  2915.           X = -S*G + C*Y
  2916.           IF (IP.EQ.0) GO TO 430
  2917. C
  2918.           DO 420 J=1,IP
  2919.             Y = B(I1,J)
  2920.             Z = B(I,J)
  2921.             B(I1,J) = Y*C + Z*S
  2922.             B(I,J) = -Y*S + Z*C
  2923.   420     CONTINUE
  2924. C
  2925.   430   CONTINUE
  2926. C
  2927.         RV1(L) = 0.0
  2928.         RV1(K) = F
  2929.         W(K) = X
  2930.         GO TO 340
  2931. C     ********** CONVERGENCE **********
  2932.   440   IF (Z.GE.0.0) GO TO 460
  2933. C     ********** W(K) IS MADE NON-NEGATIVE **********
  2934.         W(K) = -Z
  2935. C
  2936.         DO 450 J=1,N
  2937.           A(J,K) = -A(J,K)
  2938.   450   CONTINUE
  2939. C
  2940.   460 CONTINUE
  2941. C
  2942.       GO TO 480
  2943. C     ********** SET ERROR -- NO CONVERGENCE TO A
  2944. C                SINGULAR VALUE AFTER 30 ITERATIONS **********
  2945.   470 IERR = K
  2946.   480 RETURN
  2947. C     ********** LAST CARD OF MINFIT **********
  2948.       END
  2949.  
  2950.