C     PROG3                                                             PR300100
C     C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUN 12 PR300200
C     TO APPEAR IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974PR300300
C          DEMONSTRATE THE USE OF SUBROUTINE   SVDRS  TO COMPUTE THE    PR300400
C     SINGULAR VALUE DECOMPOSITION OF A MATRIX, A, AND SOLVE A LEAST    PR300500
C     SQUARES PROBLEM,  A*X=B.                                          PR300600
C                                                                       PR300700
C     THE S.V.D.  A= U*(S,0)**T*V**T IS                                 PR300800
C     COMPUTED SO THAT..                                                PR300900
C     (1)  U**T*B REPLACES THE M+1 ST. COL. OF  B.                      PR301000
C                                                                       PR301100
C     (2)  U**T   REPLACES THE M BY M IDENTITY IN                       PR301200
C     THE FIRST M COLS. OF  B.                                          PR301300
C                                                                       PR301400
C     (3) V REPLACES THE FIRST N ROWS AND COLS.                         PR301500
C     OF A.                                                             PR301600
C                                                                       PR301700
C     (4) THE DIAGONAL ENTRIES OF THE S MATRIX                          PR301800
C     REPLACE THE FIRST N ENTRIES OF THE ARRAY S.                       PR301900
C                                                                       PR302000
C     THE ARRAY S( ) MUST BE DIMENSIONED AT LEAST 3*N .                 PR302100
C                                                                       PR302200
      DIMENSION A(8,8),B(8,9),S(24),X(8),AA(8,8)                        PR302300
      REAL GEN,ANOISE                                                   PR302400
      DOUBLE PRECISION SM                                               PR302500
      DATA MDA,MDB/8,8/                                                 PR302600
C                                                                       PR302700
      DO 150 NOISE=1,2                                                  PR302800
      ANOISE = 0.                                                       PR302900
      RHO = 1.E-3                                                       PR303000
      IF(NOISE .EQ. 1) GO TO 5                                          PR303100
      ANOISE = 1.E-4                                                    PR303200
      RHO = 10. * ANOISE                                                PR303300
    5 CONTINUE                                                          PR303400
      WRITE(6,230)                                                      PR303500
      WRITE(6,240) ANOISE,RHO                                           PR303600
C     INITIALIZE DATA GENERATION FUNCTION                               PR303700
C    ..                                                                 PR303800
      DUMMY= GEN(-1.)                                                   PR303900
C                                                                       PR304000
           DO 150 MN1=1,6,5                                             PR304100
           MN2=MN1+2                                                    PR304200
                DO 150 M=MN1,MN2                                        PR304300
                DO 150 N=MN1,MN2                                        PR304400
                WRITE (6,160) M,N                                       PR304500
                     DO 20 I=1,M                                        PR304600
                          DO 10 J=1,M                                   PR304700
   10                     B(I,J)=0.                                     PR304800
                     B(I,I)=1.                                          PR304900
                          DO 20 J=1,N                                   PR305000
                          A(I,J)= GEN(ANOISE)                           PR305100
   20                     AA(I,J)=A(I,J)                                PR305200
                     DO 30 I=1,M                                        PR305300
   30                B(I,M+1)= GEN(ANOISE)                              PR305400
C                                                                       PR305500
C     THE ARRAYS ARE NOW FILLED IN..                                    PR305600
C     COMPUTE THE S.V.D.                                                PR305700
C     ******************************************************            PR305800
                CALL SVDRS (A,MDA,M,N,B,MDB,M+1,S)                      PR305900
C     ******************************************************            PR306000
C                                                                       PR306100
                WRITE (6,170)                                           PR306200
                WRITE (6,220) (I,S(I),I=1,N)                            PR306300
                WRITE (6,180)                                           PR306400
                WRITE (6,220) (I,B(I,M+1),I=1,M)                        PR306500
C                                                                       PR306600
C     TEST FOR DISPARITY OF RATIO OF SINGULAR VALUES.                   PR306700
C    ..                                                                 PR306800
                KRANK=N                                                 PR306900
                TAU=RHO*S(1)                                            PR307000
                     DO 40 I=1,N                                        PR307100
                     IF (S(I).LE.TAU) GO TO 50                          PR307200
   40                CONTINUE                                           PR307300
                GO TO 55                                                PR307400
   50           KRANK=I-1                                               PR307500
   55 WRITE(6,250) TAU, KRANK                                           PR307600
C     COMPUTE SOLUTION VECTOR ASSUMING PSEUDORANK IS KRANK              PR307700
C     ..                                                                PR307800
   60                DO 70 I=1,KRANK                                    PR307900
   70                B(I,M+1)=B(I,M+1)/S(I)                             PR308000
                     DO 90 I=1,N                                        PR308100
                     SM=0.D0                                            PR308200
                          DO 80 J=1,KRANK                               PR308300
   80                     SM=SM+A(I,J)*DBLE(B(J,M+1))                   PR308400
   90                X(I)=SM                                            PR308500
C     COMPUTE PREDICTED NORM OF RESIDUAL VECTOR.                        PR308600
C     ..                                                                PR308700
                SRSMSQ=0.                                               PR308800
                IF (KRANK.EQ.M) GO TO 110                               PR308900
                KP1=KRANK+1                                             PR309000
                     DO 100 I=KP1,M                                     PR309100
  100                SRSMSQ=SRSMSQ+B(I,M+1)**2                          PR309200
                SRSMSQ=SQRT(SRSMSQ)                                     PR309300
C                                                                       PR309400
  110           CONTINUE                                                PR309500
                WRITE (6,190)                                           PR309600
                WRITE (6,220) (I,X(I),I=1,N)                            PR309700
                WRITE (6,200) SRSMSQ                                    PR309800
C     COMPUTE THE FROBENIUS NORM OF  A**T- V*(S,0)*U**T.                PR309900
C                                                                       PR310000
C     COMPUTE  V*S  FIRST.                                              PR310100
C                                                                       PR310200
                MINMN=MIN0(M,N)                                         PR310250
                     DO 120 J=1,MINMN                                   PR310300
                          DO 120 I=1,N                                  PR310400
  120                     A(I,J)=A(I,J)*S(J)                            PR310500
                DN=0.                                                   PR310600
                     DO 140 J=1,M                                       PR310700
                          DO 140 I=1,N                                  PR310800
                          SM=0.D0                                       PR310900
                               DO 130 L=1,MINMN                         PR311000
  130                          SM=SM+A(I,L)*DBLE(B(L,J))                PR311100
C     COMPUTED DIFFERENCE OF (I,J) TH ENTRY                             PR311200
C     OF  A**T-V*(S,0)*U**T.                                            PR311300
C     ..                                                                PR311400
                          T=AA(J,I)-SM                                  PR311500
  140                     DN=DN+T**2                                    PR311600
                DN=SQRT(DN)/(SQRT(FLOAT(N))*S(1))                       PR311700
                WRITE (6,210) DN                                        PR311800
  150           CONTINUE                                                PR311900
      STOP                                                              PR312000
  160 FORMAT (1H0////9H0   M   N/1X,2I4)                                PR312100
  170 FORMAT (1H0,8X,25HSINGULAR VALUES OF MATRIX)                      PR312200
  180 FORMAT (1H0,8X,30HTRANSFORMED RIGHT SIDE, U**T*B)                 PR312300
  190 FORMAT (1H0,8X,33HESTIMATED PARAMETERS,  X=A**(+)*B,21H COMPUTED BPR312400
     1Y 'SVDRS' )                                                       PR312500
  200 FORMAT (1H0,8X,24HRESIDUAL VECTOR LENGTH =,E12.4)                 PR312600
  210 FORMAT (1H0,8X,43HFROBENIUS NORM (A-U*(S,0)**T*V**T)/(SQRT(N),    PR312700
     *22H*SPECTRAL NORM OF A) =,E12.4)                                  PR312800
  220 FORMAT (9X,I5,E16.8,I5,E16.8,I5,E16.8,I5,E16.8,I5,E16.8)          PR312900
  230 FORMAT(51H1PROG3.     THIS PROGRAM DEMONSTRATES THE ALGORITHM,    PR313000
     * 9H, SVDRS. )                                                     PR313100
  240 FORMAT(55H0THE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILL BE,PR313200
     *E16.4/44H0THE RELATIVE TOLERANCE, RHO, FOR PSEUDORANK,            PR313300
     *17H DETERMINATION IS,E16.4)                                       PR313400
  250 FORMAT(1H0,8X,36HABSOLUTE PSEUDORANK TOLERANCE, TAU =,            PR313500
     *E12.4,10X,12HPSEUDORANK =,I5)                                     PR313600
      END                                                               PR313700
