home *** CD-ROM | disk | FTP | other *** search
- 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
