program main c*********************************************************************72 c cc MAIN is the main program for LAWSON_PRB3. c C DEMONSTRATE THE USE OF SUBROUTINE SVDRS TO COMPUTE THE C SINGULAR VALUE DECOMPOSITION OF A MATRIX, A, AND SOLVE A LEAST C SQUARES PROBLEM, A*X=B. c c Modified: c c 19 October 2008 c c Author: c c Charles Lawson, Richard Hanson c c Reference: c c Charles Lawson, Richard Hanson, c Solving Least Squares Problems, c SIAM, 1995, c ISBN: 0898713560, c LC: QA275.L38. c C THE S.V.D. A= U*(S,0)**T*V**T IS C COMPUTED SO THAT.. C (1) U**T*B REPLACES THE M+1 ST. COL. OF B. C C (2) U**T REPLACES THE M BY M IDENTITY IN C THE FIRST M COLS. OF B. C C (3) V REPLACES THE FIRST N ROWS AND COLS. C OF A. C C (4) THE DIAGONAL ENTRIES OF THE S MATRIX C REPLACE THE FIRST N ENTRIES OF THE ARRAY S. C C THE ARRAY S( ) MUST BE DIMENSIONED AT LEAST 3*N . C integer MDA, MDB, MX parameter(MDA = 8, MDB = 8, MX = 8) integer I, J, KP1, KRANK, L, M, MINMN, MN1, MN2, N, NOISE double precision A(MDA, MX), AA(MDA, MX), ANOISE, B(MDB, MDB+1) double precision DN, DUMMY, FLOATN, GEN, ONE, RHO double precision S(MX), SM, SMALL3, SMALL4, SRSMSQ double precision T, TAU, TEN, WORK(2*MX), X(MX), ZERO parameter(ONE = 1.0d0, SMALL3 = 1.0d-3, SMALL4 = 1.0d-4) parameter(TEN = 10.0d0) parameter(ZERO = 0.0d0) 170 format (/////' M N'/1X,2I4) 180 format (/' SINGULAR VALUES OF MATRIX') 190 format (/' TRANSFORMED RIGHT SIDE, U**T*B') 200 format (/ * ' ESTIMATED PARAMETERS, X=A**(+)*B, COMPUTED BY ''SVDRS''') 210 format (/' RESIDUAL VECTOR LENGTH = ',E12.4) 220 format (/ * ' FROBENIUS NORM(A-U*(S,0)**T*V**T)'/ * ' --------------------------------- = ',g12.4/ * ' SQRT(N) * SPECTRAL NORM OF A') 230 format (2X,I5,E16.8,I5,E16.8,I5,E16.8) 240 format (/ * ' PROG3. THIS PROGRAM DEMONSTRATES THE ALGORITHM SVDRS.'/) 250 format (/ * ' THE RELATIVE NOISE LEVEL OF THE GENERATED DATA WILL BE', * g12.4/ * ' THE RELATIVE TOLERANCE FOR PSEUDORANK DETERMINATION IS RHO =', * g12.4) 260 format (/' ABSOLUTE PSEUDORANK TOLERANCE, TAU = ',g12.4/ * ' PSEUDORANK = ',I5) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB3' DO 166 NOISE=1,2 if(NOISE .eq. 1) then ANOISE= ZERO RHO= SMALL3 else ANOISE= SMALL4 RHO= TEN * ANOISE endif write (*,240) write (*,250) ANOISE,RHO c C INITIALIZE DATA GENERATION FUNCTION C DUMMY=GEN(-ONE) DO 164 MN1=1,6,5 MN2=MN1+2 DO 162 M=MN1,MN2 DO 160 N=MN1,MN2 write (*,170) M,N DO 35 I=1,M DO 20 J=1,M 20 B(I,J)= ZERO B(I,I)= ONE DO 30 J=1,N A(I,J)=GEN(ANOISE) AA(I,J)=A(I,J) 30 continue 35 continue DO 40 I=1,M 40 B(I,M+1)=GEN(ANOISE) C C The arrays are now filled in.. C Compute the singular value decomposition. C CALL SVDRS (A,MDA,M,N,B,MDB,M+1,S, WORK) write (*,180) write (*,230) (I,S(I),I=1,N) write (*,190) write (*,230) (I,B(I,M+1),I=1,M) C C TEST FOR DISPARITY OF RATIO OF SINGULAR VALUES. C KRANK=N TAU=RHO*S(1) DO 50 I=1,N IF (S(I).LE.TAU) GO TO 60 50 CONTINUE GO TO 70 60 KRANK=I-1 70 continue write (*,260) TAU,KRANK c C COMPUTE SOLUTION VECTOR ASSUMING PSEUDORANK IS KRANK C DO 80 I=1,KRANK 80 B(I,M+1)=B(I,M+1)/S(I) DO 100 I=1,N SM= ZERO DO 90 J=1,KRANK 90 SM=SM+A(I,J)*B(J,M+1) 100 X(I)=SM c C COMPUTE PREDICTED NORM OF RESIDUAL VECTOR. C SRSMSQ= ZERO IF (KRANK .ne. M) then KP1=KRANK+1 DO 110 I=KP1,M 110 SRSMSQ=SRSMSQ+B(I,M+1)**2 SRSMSQ=sqrt(SRSMSQ) endif write (*,200) write (*,230) (I,X(I),I=1,N) write (*,210) SRSMSQ C COMPUTE THE FROBENIUS NORM OF A**T- V*(S,0)*U**T. C C COMPUTE V*S FIRST. C MINMN=min(M,N) DO 130 J=1,MINMN DO 130 I=1,N A(I,J)=A(I,J)*S(J) 130 continue 135 continue DN= ZERO DO 150 J=1,M DO 145 I=1,N SM= ZERO DO 140 L=1,MINMN 140 SM=SM+A(I,L)*B(L,J) c C COMPUTED DIFFERENCE OF (I,J) TH ENTRY C OF A**T-V*(S,0)*U**T. C T=AA(J,I)-SM DN=DN+T**2 145 continue 150 continue FLOATN = N DN=sqrt(DN)/(sqrt(FLOATN)*S(1)) write (*,220) DN 160 continue 162 continue 164 continue 166 continue write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LAWSON_PRB3:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end