home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / tests / rltest.f < prev    next >
Text File  |  1996-09-28  |  7KB  |  252 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL MATRICES
  3. C     EXHIBITING THE USE OF  EISPACK  TO FIND THE SINGULAR VALUES
  4. C     AND THE SOLUTION TO THE EQUATION  A*X = B .  THIS DRIVER
  5. C     SUMMARIZES THE FIGURES OF MERIT FOR ALL PATHS.
  6. C
  7. C     THIS DRIVER IS CATALOGUED AS EISPDRV4(RLSUMARY).
  8. C
  9. C     THE DIMENSION OF  A,U,  AND  V  SHOULD BE  NM  BY  NM.
  10. C     THE DIMENSION OF  SIGMA  AND  RV1  SHOULD BE  NM.
  11. C     THE DIMENSION OF  AHOLD  SHOULD BE  NM  BY  NM.
  12. C     HERE  NM = 20.
  13. C
  14.       DOUBLE PRECISION A( 20, 20),U( 20, 20),V( 20, 20),
  15.      X        SIGMA( 20),RV1( 20),AHOLD( 20, 20),TCRIT( 2),
  16.      X        RESDUL,EPSLON,DFL
  17.       INTEGER IERR( 2),ERROR
  18.       DATA IREAD1/1/,IWRITE/6/
  19. C
  20.       OPEN(UNIT=IREAD1,FILE='FILE48')
  21.       REWIND IREAD1
  22. C
  23.       NM = 20
  24.       LCOUNT = 0
  25.       WRITE(IWRITE,1)
  26.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  27.      XTATISTICS//1H ,95(1H-)/  31H M   N  IP   SVD         MINFIT /
  28.      X1H ,95(1H-)//
  29.      X53H UNDER 'M' IS THE ROW DIMENSION OF EACH TEST MATRIX. //
  30.      X49H UNDER 'N' IS THE COLUMN DIMENSION OF THE MATRIX. //
  31.      X93H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX.  (A
  32.      XS PRESENTLY CONSTITUTED,      /
  33.      X58H THE IDENTITY MATRIX IS GENERATED AS THE CONSTANT MATRIX.) //
  34.      X93H UNDER 'SVD' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER, IS
  35.      X THE ERROR FLAG RETURNED      /
  36.      X93H FROM  SVD.  THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BA
  37.      XSED UPON A RESIDUAL.          //
  38.      X93H UNDER 'MINFIT' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,
  39.      X IS THE ERROR FLAG RETURNED   /
  40.      X93H FROM  MINFIT.  THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE
  41.      X BASED UPON A RESIDUAL.       //
  42.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  43.      X27H THE CORRESPONDING PATH HAS        /
  44.      X50H PREVENTED THE COMPUTATION OF THE SINGULAR VALUES.  /)
  45.       WRITE(IWRITE,15)
  46.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  47.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  48.      X26HMEASURE OF PERFORMANCE FOR /5X, 21HTHE  EISPACK  CODES. ,
  49.      X60H THIS RUN DISPLAYS THESE STATISTICS FOR REAL LINEAR SYSTEMS.//
  50.      X33H   M   N  IP   SVD         MINFIT  )
  51.    10 CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,0)
  52. C
  53.       DO 230 ICALL = 1,2
  54.          IF( ICALL .EQ. 1 ) GO TO 90
  55.          CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1)
  56.          GO TO 100
  57. C
  58. C     RL  USING  SVD
  59. C
  60.    90    ICT = 1
  61.          CALL  SVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,ERROR,RV1)
  62.          IERR(ICT) = ERROR
  63.          IF( IERR(ICT) .NE. 0 ) GO TO 225
  64.          GO TO 226
  65. C
  66. C     RL  USING  MINFIT
  67. C
  68.   100    ICT = 2
  69.          IP = M
  70.          DO 102 I = 1,M
  71.             DO 101 J = 1,IP
  72.                V(I,J) = 0.0D0
  73.   101       CONTINUE
  74.             V(I,I) = 1.0D0
  75.   102    CONTINUE
  76.          CALL  MINFIT(NM,M,N,A,SIGMA,IP,V,ERROR,RV1)
  77.          IERR(ICT) = ERROR
  78.          IF( IERR(ICT) .NE. 0 ) GO TO 225
  79.          DO 104 I = 1,M
  80.             DO 103 J = 1,N
  81.                U(I,J) = V(J,I)
  82.   103       CONTINUE
  83.   104    CONTINUE
  84.          DO 106 I = 1,N
  85.             DO 105 J = 1,N
  86.                V(I,J) = A(I,J)
  87.   105       CONTINUE
  88.   106    CONTINUE
  89.          GO TO 226
  90.   225    TCRIT(ICT) = -1.0D0
  91.          GO TO 230
  92. C
  93.   226    CALL  RMATIN(NM,NM,M,N,A,AHOLD,AHOLD,0,1)
  94.          CALL  RLRES(NM,M,N,A,U,SIGMA,V,RV1,RESDUL)
  95.          DFL = 10 * MAX0(M,N)
  96.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  97.   230 CONTINUE
  98. C
  99.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  100.       LCOUNT = LCOUNT + 1
  101.       WRITE(IWRITE,520) M,N,IP,IERR(1),TCRIT(1),IERR(2),TCRIT(2)
  102.   520 FORMAT(3I4,I4,F8.3,I4,F8.3)
  103.       GO TO 10
  104.       END
  105.       SUBROUTINE RLRES(NM,M,N,A,U,S,V,NORM,RESDUL)
  106. C
  107.       DOUBLE PRECISION A(NM,N),U(NM,N),V(NM,N),S(N),NORM(N)
  108.       DOUBLE PRECISION SUM1,SUM2,X,RESDUL,NORMA,SUMA,SUMUV
  109. C
  110. C     THIS SUBROUTINE PERFORMS A RESIDUAL TEST FOR THE DECOMPOSITION
  111. C                              T
  112. C     OF A MATRIX  A  INTO  USV .
  113. C
  114. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RLRES).
  115. C
  116. C     INPUT.
  117. C
  118. C       NM  IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  119. C         AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  120. C
  121. C       M  IS THE ROW DIMENSION OF THE MATRICES  A  AND  U.
  122. C
  123. C       N  IS THE COLUMN DIMENSION OF THE MATRICES  A  AND U , AND
  124. C         THE ORDER OF  V;
  125. C
  126. C       A  IS A REAL GENERAL MATRIX OF DIMEMSION  M  BY  N;
  127. C
  128. C       U  IS A MATRIX WITH ORTHOGONAL COLUMNS OF DIMENSION  M  BY  N;
  129. C
  130. C       S  IS A DIAGONAL MATRIX STORED AS A VECTOR OF DIMENSION  N;
  131. C
  132. C       V  IS AN ORTHOGONAL MATRIX OF DIMENSION  N  BY  N.
  133. C
  134. C     OUTPUT.
  135. C
  136. C       NORM  IS AN ARRAY SUCH THAT FOR EACH K,
  137. C                                              T
  138. C         NORM(K) = !!A*V(K)-S(K)*U(K)!!  + !!A *U(K)-S(K)*V(K)!!
  139. C                   ---------------------------------------------
  140. C                           !!A!!*(!!U(K)!! + !!V(K)!!)
  141. C
  142. C       RESDUL  IS THE LARGEST ELEMENT OF  NORM.
  143. C
  144. C     ------------------------------------------------------------------
  145. C
  146.       RESDUL = 0.0D0
  147.       NORMA = 0.0D0
  148. C
  149.       DO 20 I = 1,M
  150.          SUMA = 0.0D0
  151. C
  152.          DO 10 J = 1,N
  153.             SUMA = SUMA + DABS(A(I,J))
  154.    10    CONTINUE
  155. C
  156.          NORMA = DMAX1(NORMA,SUMA)
  157.    20 CONTINUE
  158. C
  159.       IF( NORMA .EQ. 0.0D0 ) NORMA = 1.0D0
  160. C
  161.       DO 70 K = 1,N
  162.          SUMUV = 0.0D0
  163.          SUM1 = 0.0D0
  164. C
  165.          DO 40 I = 1,M
  166.             SUMUV = SUMUV + DABS(U(I,K))
  167.             X = -S(K)*U(I,K)
  168. C
  169.             DO 30 J = 1,N
  170.                X = X + A(I,J)*V(J,K)
  171.    30       CONTINUE
  172. C
  173.             SUM1 = SUM1 + DABS(X)
  174.    40    CONTINUE
  175. C
  176.          SUM2 = 0.0D0
  177. C
  178.          DO 60 I = 1,N
  179.             SUMUV = SUMUV + DABS(V(I,K))
  180.             X = -S(K)*V(I,K)
  181. C
  182.             DO 50 J = 1,M
  183.                X = X + A(J,I)*U(J,K)
  184.    50       CONTINUE
  185. C
  186.             SUM2 = SUM2 + DABS(X)
  187.    60    CONTINUE
  188. C
  189.          NORM(K) = (SUM1 + SUM2)/(SUMUV*NORMA)
  190.          RESDUL = DMAX1(NORM(K),RESDUL)
  191.    70 CONTINUE
  192. C
  193.       RETURN
  194.       END
  195.       SUBROUTINE RMATIN(NM,NN,M,N,A,AHOLD,BHOLD,AORB,INITIL)
  196. C
  197. C     THIS INPUT SUBROUTINE READS A REAL MATRIX, WITH DIMENSION
  198. C     M  BY  N, FROM SYSIN.
  199. C     TO GENERATE THE MATRIX  A (OR  B)  INITIALLY,  INITIL  IS
  200. C     SET TO 0.  TO REGENERATE THE MATRIX  A (OR  B)  FOR THE
  201. C     PURPOSE OF THE RESIDUAL CALCULATION,  INITIL  IS SET TO 1.
  202. C
  203. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RLREADI).
  204. C
  205.       INTEGER IA( 20),NM,NN,M,N,AORB,INITIL
  206.       DOUBLE PRECISION A(NM,NN),AHOLD(NM,NN),BHOLD(NM,NM)
  207.       DATA IREADA/1/,IWRITE/6/
  208. C
  209.       IF( INITIL .EQ. 1 ) GO TO 90
  210.       READ(IREADA,5) M,N
  211.     5 FORMAT(I6,6X,I6)
  212.       IF( M .EQ. 0 ) GO TO 170
  213.       DO 20 I = 1,M
  214.          READ(IREADA,10) (IA(J),J=1,N)
  215.    10    FORMAT(6I12)
  216.          DO 15 J = 1,N
  217.             A(I,J) = IA(J)
  218.    15    CONTINUE
  219.    20 CONTINUE
  220.       IF( AORB .EQ. 1 ) GO TO 50
  221.       DO 40 J = 1,N
  222.          DO 30 I = 1,M
  223.             AHOLD(I,J) = A(I,J)
  224.    30    CONTINUE
  225.    40 CONTINUE
  226.       GO TO 80
  227.    50 DO 70 J = 1,N
  228.          DO 60 I = 1,M
  229.             BHOLD(I,J) = A(I,J)
  230.    60    CONTINUE
  231.    70 CONTINUE
  232.    80 RETURN
  233. C
  234.    90 IF( AORB .EQ. 1 ) GO TO 120
  235.       DO 110 J = 1,N
  236.          DO 100 I = 1,M
  237.             A(I,J) = AHOLD(I,J)
  238.   100    CONTINUE
  239.   110 CONTINUE
  240.       GO TO 150
  241.   120 DO 140 J = 1,N
  242.          DO 130 I = 1,M
  243.             A(I,J) = BHOLD(I,J)
  244.   130    CONTINUE
  245.   140 CONTINUE
  246.   150 RETURN
  247. C
  248.   170 WRITE(IWRITE,175)
  249.   175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1)
  250.       STOP
  251.       END
  252.