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