C C THIS DRIVER TESTS EISPACK FOR THE CLASS OF REAL MATRICES C EXHIBITING THE USE OF EISPACK TO FIND THE SINGULAR VALUES C AND THE SOLUTION TO THE EQUATION A*X = B . THIS DRIVER C SUMMARIZES THE FIGURES OF MERIT FOR ALL PATHS. C C THIS DRIVER IS CATALOGUED AS EISPDRV4(RLSUMARY). C C THE DIMENSION OF A,U, AND V SHOULD BE NM BY NM. C THE DIMENSION OF SIGMA AND RV1 SHOULD BE NM. C THE DIMENSION OF AHOLD SHOULD BE NM BY NM. C HERE NM = 20. C REAL A( 20, 20),U( 20, 20),V( 20, 20), X SIGMA( 20),RV1( 20),AHOLD( 20, 20),TCRIT( 2), X RESDUL,EPSLON,DFL INTEGER IERR( 2),ERROR DATA IREAD1/1/,IWRITE/6/ C OPEN(UNIT=IREAD1,FILE='FILE48') REWIND IREAD1 C NM = 20 LCOUNT = 0 WRITE(IWRITE,1) 1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S XTATISTICS//1H ,95(1H-)/ 31H M N IP SVD MINFIT / X1H ,95(1H-)// X53H UNDER 'M' IS THE ROW DIMENSION OF EACH TEST MATRIX. // X49H UNDER 'N' IS THE COLUMN DIMENSION OF THE MATRIX. // X93H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX. (A XS PRESENTLY CONSTITUTED, / X58H THE IDENTITY MATRIX IS GENERATED AS THE CONSTANT MATRIX.) // X93H UNDER 'SVD' ARE TWO NUMBERS. THE FIRST NUMBER, AN INTEGER, IS X THE ERROR FLAG RETURNED / X93H FROM SVD. THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BA XSED UPON A RESIDUAL. // X93H UNDER 'MINFIT' ARE TWO NUMBERS. THE FIRST NUMBER, AN INTEGER, X IS THE ERROR FLAG RETURNED / X93H FROM MINFIT. THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE X BASED UPON A RESIDUAL. // X62H -1.0 AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN, X27H THE CORRESPONDING PATH HAS / X50H PREVENTED THE COMPUTATION OF THE SINGULAR VALUES. /) WRITE(IWRITE,15) 15 FORMAT(1X,21HS.P. VERSION 04/15/83 ) 5 FORMAT( 53H1 TABULATION OF THE ERROR FLAG ERROR AND THE , X26HMEASURE OF PERFORMANCE FOR /5X, 21HTHE EISPACK CODES. , X60H THIS RUN DISPLAYS THESE STATISTICS FOR REAL LINEAR SYSTEMS.// X33H M N IP SVD MINFIT ) 10 CALL RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,0) C DO 230 ICALL = 1,2 IF( ICALL .EQ. 1 ) GO TO 90 CALL RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1) GO TO 100 C C RL USING SVD C 90 ICT = 1 CALL SVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,ERROR,RV1) IERR(ICT) = ERROR IF( IERR(ICT) .NE. 0 ) GO TO 225 GO TO 226 C C RL USING MINFIT C 100 ICT = 2 IP = M DO 102 I = 1,M DO 101 J = 1,IP V(I,J) = 0.0E0 101 CONTINUE V(I,I) = 1.0E0 102 CONTINUE CALL MINFIT(NM,M,N,A,SIGMA,IP,V,ERROR,RV1) IERR(ICT) = ERROR IF( IERR(ICT) .NE. 0 ) GO TO 225 DO 104 I = 1,M DO 103 J = 1,N U(I,J) = V(J,I) 103 CONTINUE 104 CONTINUE DO 106 I = 1,N DO 105 J = 1,N V(I,J) = A(I,J) 105 CONTINUE 106 CONTINUE GO TO 226 225 TCRIT(ICT) = -1.0E0 GO TO 230 C 226 CALL RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1) CALL RLRES(NM,M,N,A,U,SIGMA,V,RV1,RESDUL) DFL = 10 * MAX0(M,N) TCRIT(ICT) = RESDUL/EPSLON(DFL) 230 CONTINUE C IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5) LCOUNT = LCOUNT + 1 WRITE(IWRITE,520) M,N,IP,IERR(1),TCRIT(1),IERR(2),TCRIT(2) 520 FORMAT(3I4,I4,F8.3,I4,F8.3) GO TO 10 END SUBROUTINE RLRES(NM,M,N,A,U,S,V,NORM,RESDUL) C REAL A(NM,N),U(NM,N),V(NM,N),S(N),NORM(N) REAL SUM1,SUM2,X,RESDUL,NORMA,SUMA,SUMUV C C THIS SUBROUTINE PERFORMS A RESIDUAL TEST FOR THE DECOMPOSITION C T C OF A MATRIX A INTO USV . C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RLRES). C C INPUT. C C NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT; C C M IS THE ROW DIMENSION OF THE MATRICES A AND U. C C N IS THE COLUMN DIMENSION OF THE MATRICES A AND U , AND C THE ORDER OF V; C C A IS A REAL GENERAL MATRIX OF DIMEMSION M BY N; C C U IS A MATRIX WITH ORTHOGONAL COLUMNS OF DIMENSION M BY N; C C S IS A DIAGONAL MATRIX STORED AS A VECTOR OF DIMENSION N; C C V IS AN ORTHOGONAL MATRIX OF DIMENSION N BY N. C C OUTPUT. C C NORM IS AN ARRAY SUCH THAT FOR EACH K, C T C NORM(K) = !!A*V(K)-S(K)*U(K)!! + !!A *U(K)-S(K)*V(K)!! C --------------------------------------------- C !!A!!*(!!U(K)!! + !!V(K)!!) C C RESDUL IS THE LARGEST ELEMENT OF NORM. C C ------------------------------------------------------------------ C RESDUL = 0.0E0 NORMA = 0.0E0 C DO 20 I = 1,M SUMA = 0.0E0 C DO 10 J = 1,N SUMA = SUMA + ABS(A(I,J)) 10 CONTINUE C NORMA = AMAX1(NORMA,SUMA) 20 CONTINUE C IF( NORMA .EQ. 0.0E0 ) NORMA = 1.0E0 C DO 70 K = 1,N SUMUV = 0.0E0 SUM1 = 0.0E0 C DO 40 I = 1,M SUMUV = SUMUV + ABS(U(I,K)) X = -S(K)*U(I,K) C DO 30 J = 1,N X = X + A(I,J)*V(J,K) 30 CONTINUE C SUM1 = SUM1 + ABS(X) 40 CONTINUE C SUM2 = 0.0E0 C DO 60 I = 1,N SUMUV = SUMUV + ABS(V(I,K)) X = -S(K)*V(I,K) C DO 50 J = 1,M X = X + A(J,I)*U(J,K) 50 CONTINUE C SUM2 = SUM2 + ABS(X) 60 CONTINUE C NORM(K) = (SUM1 + SUM2)/(SUMUV*NORMA) RESDUL = AMAX1(NORM(K),RESDUL) 70 CONTINUE C RETURN END SUBROUTINE RMATIN(NM,NN,M,N,A,AHOLD,BHOLD,AORB,INITIL) C C THIS INPUT SUBROUTINE READS A REAL MATRIX, WITH DIMENSION C M BY N, FROM SYSIN. C TO GENERATE THE MATRIX A (OR B) INITIALLY, INITIL IS C SET TO 0. TO REGENERATE THE MATRIX A (OR B) FOR THE C PURPOSE OF THE RESIDUAL CALCULATION, INITIL IS SET TO 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RLREADI). C INTEGER IA( 20),NM,NN,M,N,AORB,INITIL REAL A(NM,NN),AHOLD(NM,NN),BHOLD(NM,NM) DATA IREADA/1/,IWRITE/6/ C IF( INITIL .EQ. 1 ) GO TO 90 READ(IREADA,5) M,N 5 FORMAT(I6,6X,I6) IF( M .EQ. 0 ) GO TO 170 DO 20 I = 1,M READ(IREADA,10) (IA(J),J=1,N) 10 FORMAT(6I12) DO 15 J = 1,N A(I,J) = IA(J) 15 CONTINUE 20 CONTINUE IF( AORB .EQ. 1 ) GO TO 50 DO 40 J = 1,N DO 30 I = 1,M AHOLD(I,J) = A(I,J) 30 CONTINUE 40 CONTINUE GO TO 80 50 DO 70 J = 1,N DO 60 I = 1,M BHOLD(I,J) = A(I,J) 60 CONTINUE 70 CONTINUE 80 RETURN C 90 IF( AORB .EQ. 1 ) GO TO 120 DO 110 J = 1,N DO 100 I = 1,M A(I,J) = AHOLD(I,J) 100 CONTINUE 110 CONTINUE GO TO 150 120 DO 140 J = 1,N DO 130 I = 1,M A(I,J) = BHOLD(I,J) 130 CONTINUE 140 CONTINUE 150 RETURN C 170 WRITE(IWRITE,175) 175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1) STOP END