home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / eispack-1.0-src.tgz / tar.out / contrib / eispack / ex / rstest.f < prev    next >
Text File  |  1996-09-28  |  13KB  |  392 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL SYMMETRIC
  3. C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RSSUMARY).
  6. C
  7. C     THE DIMENSION OF  A  AND  Z  SHOULD BE  NM  BY  NM.
  8. C     THE DIMENSION OF  W,D,E,E2,IND,RV1,RV2,RV3,RV4,RV5,RV6,
  9. C     W1, AND  W2  SHOULD BE  NM.
  10. C     THE DIMENSION OF  AHOLD  SHOULD BE  NM  BY  NM.
  11. C     HERE NM = 20.
  12. C
  13.       DOUBLE PRECISION A(20,20),Z(20,20),AHOLD(20,20),W(20),D(20),E(20),
  14.      X        E2( 20),RV1( 20),RV2( 20),RV3( 20),RV4( 20),RV5( 20),
  15.      X        RV6( 20),W1( 20),W2( 20),TCRIT( 8),EPSLON,RESDUL,MAXEIG,
  16.      X        MAXDIF,U,LB,UB,EPS1,FWORK(1600),DFL
  17.       REAL  XUB,XLB
  18.       INTEGER  IND( 20),IERR( 6),ERROR
  19.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  20. C
  21.       OPEN(UNIT=IREAD1,FILE='FILE35')
  22.       OPEN(UNIT=IREADC,FILE='FILE36')
  23.       REWIND IREAD1
  24.       REWIND IREADC
  25. C
  26.       NM = 20
  27.       LCOUNT = 0
  28.       WRITE(IWRITE,1)
  29.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  30.      XTATISTICS//1H ,95(1H-)/96H ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB
  31.      X      UB    M  IMTQLV   TSTURM   BISECT  M1 NO  TRIDIB  /1H ,
  32.      X95(1H-)//48H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
  33.      X95H UNDER 'TQL2   TQLRAT' ARE THREE NUMBERS.  THE FIRST NUMBER, AN
  34.      X INTEGER, IS THE ABSOLUTE SUM OF/
  35.      X61H THE ERROR FLAGS RETURNED SEPARATELY FROM  TQL2  AND  TQLRAT.,
  36.      X33H THE SECOND NUMBER IS THE MEASURE/
  37.      X62H OF PERFORMANCE BASED UPON THE RESIDUAL COMPUTED FOR THE  TQL2,
  38.      X25H  PATH.  THE THIRD NUMBER        /
  39.      X62H MEASURES THE AGREEMENT OF THE EIGENVALUES FROM THE  TQL2  AND,
  40.      X16H  TQLRAT  PATHS.  //
  41.      X95H UNDER 'IMTQL2 IMTQL1' ARE THREE NUMBERS WITH MEANING LIKE THOS
  42.      XE UNDER  'TQL2   TQLRAT'.       //
  43.      X95H UNDER 'LB' AND 'UB' ARE THE INPUT VARIABLES SPECIFYING THE INT
  44.      XERVAL TO  BISECT  AND  TSTURM.  //
  45.      X61H UNDER 'M' IS THE NUMBER OF EIGENVALUES DETERMINED BY  BISECT,
  46.      X29H AND  TSTURM  THAT LIE IN THE    /18H INTERVAL (LB,UB).//
  47.      X95H UNDER EACH OF 'IMTQLV', 'TSTURM', 'BISECT', AND 'TRIDIB' ARE T
  48.      XWO NUMBERS.  THE FIRST NUMBER,       )
  49.       WRITE(IWRITE,2)
  50.     2 FORMAT(
  51.      X95H AN INTEGER, IS THE ABSOLUTE SUM OF THE ERROR FLAGS RETURNED FR
  52.      XOM THE RESPECTIVE PATH.         /
  53.      X95H THE SECOND NUMBER IS THE MEASURE OF PERFORMANCE BASED UPON THE
  54.      X RESIDUAL COMPUTED FOR THE PATH.//
  55.      X95H UNDER 'M1' AND 'NO' ARE THE VARIABLES SPECIFYING THE LOWER BOU
  56.      XNDARY INDEX AND THE NUMBER       /
  57.      X28H OF EIGENVALUES TO  TRIDIB.   //
  58.      X62H -1.0  AS THE MEASURE OF PERFORMANCE IS PRINTED IF AN ERROR IN,
  59.      X27H THE CORRESPONDING PATH HAS        /
  60.      X47H PREVENTED THE COMPUTATION OF THE EIGENVECTORS. //
  61.      X62H THE  TQL2    PATH USES THE EISPACK CODES   TRED2-TQL2  ,     /
  62.      X38H AS CALLED FROM DRIVER SUBROUTINE  RS. /
  63.      X62H THE  TQLRAT  PATH USES THE EISPACK CODES   TRED1-TQLRAT,     /
  64.      X38H AS CALLED FROM DRIVER SUBROUTINE  RS. /
  65.      X62H THE  IMTQL2  PATH USES THE EISPACK CODES   TRED2-IMTQL2.     /
  66.      X62H THE  IMTQL1  PATH USES THE EISPACK CODES   TRED1-IMTQL1.     /
  67.      X63H THE  IMTQLV  PATH USES THE EISPACK CODES   TRED1-IMTQLV-TINVIT
  68.      X ,8H-TRBAK , /
  69.      X38H AS CALLED FROM DRIVER SUBROUTINE RSM. )
  70.       WRITE(IWRITE,3)
  71.     3 FORMAT(
  72.      X64H THE  TSTURM  PATH USES THE EISPACK CODES   TRED1-TSTURM-TRBAK1
  73.      X.  /
  74.      X63H THE  BISECT  PATH USES THE EISPACK CODES   TRED1-BISECT-TINVIT
  75.      X ,8H-TRBAK1. /
  76.      X63H THE  TRIDIB  PATH USES THE EISPACK CODES   TRED1-TRIDIB-TINVIT
  77.      X ,8H-TRBAK1. /)
  78.       WRITE(IWRITE,15)
  79.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  80.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  81.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  82.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  83.      X    33H FOR REAL SYMMETRIC MATRICES.        /
  84.      X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
  85.      X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
  86.    10 CALL RMATIN(NM,N,A,AHOLD,0)
  87.       READ(IREADC,50) MM,LB,UB,M11,NO
  88.    50 FORMAT(I4,2D24.16,2(4X,I4))
  89. C
  90. C     MM,LB,UB,M11,  AND  NO  ARE READ FROM SYSIN AFTER THE MATRIX IS
  91. C     GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE MAXIMUM
  92. C     NUMBER OF EIGENVALUES AND THE BOUNDS FOR THE INTERVAL WHICH IS
  93. C     TO BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE LOWER
  94. C     BOUNDARY INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
  95. C
  96.       DO  230  ICALL = 1,10
  97.          IF( ICALL .NE. 1 )  CALL  RMATIN(NM,N,A,AHOLD,1)
  98. C
  99. C     IF  TQLRAT  PATH (LABEL 80) IS TAKEN THEN  TQL2  PATH (LABEL 70)
  100. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  101. C     MEANINGFUL.
  102. C     IF  IMTQL1  PATH (LABEL 85) IS TAKEN THEN  IMTQL2  PATH (LABEL 75)
  103. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  104. C     MEANINGFUL.
  105. C     IF  TQL2  (IMTQL2)  PATH FAILS, THEN  TQLRAT  (IMTQL1)  PATH IS
  106. C     OMITTED AND PRINTOUT FLAGGED WITH  -1.0.
  107. C
  108.          GO TO  (70,75,80,85,89,90,95,230,110,230),  ICALL
  109. C
  110. C     RSWZ  USING  TQL2
  111. C     INVOKED FROM DRIVER SUBROUTINE  RS.
  112. C
  113.    70    ICT = 1
  114.          CALL  RS(NM,N,A,W,1,Z,E,E,ERROR)
  115.          IERR(ICT) = ERROR
  116.          M = ERROR - 1
  117.          IF( ERROR .NE. 0 ) GO TO 190
  118.          DO  73  I=1,N
  119.    73      W1(I) = W(I)
  120.          M = N
  121.          GO TO  190
  122. C
  123. C     RSWZ  USING  IMTQL2
  124. C
  125.    75    ICT = 2
  126.          CALL   TRED2(NM,N,A,W,E,Z)
  127.          CALL  IMTQL2(NM,N,W,E,Z,ERROR)
  128.          IERR(ICT) = ERROR
  129.          M = ERROR - 1
  130.          IF( ERROR .NE. 0 ) GO TO 190
  131.          DO 76 I = 1,N
  132.             W2(I) = W(I)
  133.    76    CONTINUE
  134.          M = N
  135.          GO TO  190
  136. C
  137. C     RSW  USING  TQLRAT
  138. C     INVOKED FROM DRIVER SUBROUTINE  RS.
  139. C
  140.    80    ICT = 7
  141.          IF( IERR(1) .NE. 0 ) GO TO 200
  142.          CALL  RS(NM,N,A,W,0,A,E,E2,ERROR)
  143.          IERR(1) = ERROR
  144.          IF( ERROR .NE. 0 ) GO TO 200
  145.          MAXEIG = 0.0D0
  146.          MAXDIF = 0.0D0
  147.          DO 81 I = 1,N
  148.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  149.             U = DABS(W1(I) - W(I))
  150.             IF( U .GT. MAXDIF ) MAXDIF = U
  151.    81    CONTINUE
  152.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  153.          DFL = 10 * N
  154.          TCRIT(7) = MAXDIF/EPSLON(MAXEIG*DFL)
  155.          GO TO  230
  156. C
  157. C     RSW  USING  IMTQL1
  158. C
  159.    85    ICT = 8
  160.          IF( IERR(2) .NE. 0 ) GO TO 200
  161.          CALL   TRED1(NM,N,A,W,E,E)
  162.          CALL  IMTQL1(N,W,E,ERROR)
  163.          IERR(2) = ERROR
  164.          MAXEIG = 0.0D0
  165.          MAXDIF = 0.0D0
  166.          DO 86 I = 1,N
  167.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  168.             U = DABS(W2(I) - W(I))
  169.             IF( U .GT. MAXDIF ) MAXDIF = U
  170.    86    CONTINUE
  171.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  172.          DFL = 10 * N
  173.          TCRIT(8) = MAXDIF/EPSLON(MAXEIG*DFL)
  174.          GO TO  230
  175. C
  176. C     RSW1Z  ( USAGE HERE COMPUTES ALL THE EIGENVECTORS )
  177. C     INVOKED FROM DRIVER SUBROUTINE  RSM.
  178. C
  179.    89    ICT = 3
  180.          M = N
  181.          CALL  RSM(NM,N,A,W,M,Z,FWORK,IND,ERROR)
  182.          IERR(ICT) = IABS(ERROR)
  183.          GO TO 190
  184. C
  185. C     RS1W1Z  USING  TSTURM
  186. C
  187.    90    ICT = 4
  188.          EPS1 = 0.0D0
  189.          CALL   TRED1(NM,N,A,D,E,E2)
  190.          CALL  TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,ERROR,
  191.      X                RV1,RV2,RV3,RV4,RV5,RV6)
  192.          IERR(ICT) = ERROR
  193.          XLB = LB
  194.          XUB = UB
  195.          IF( ERROR .EQ. 3*N + 1 ) GO TO 200
  196.          IF( ERROR .GT. 4*N ) M = ERROR - 4*N - 1
  197.          CALL  TRBAK1(NM,N,A,E,M,Z)
  198.          GO TO  190
  199. C
  200. C     RS1W1Z  USING  BISECT  AND  TINVIT
  201. C
  202.    95    ICT = 5
  203.          EPS1 = 0.0D0
  204.          CALL   TRED1(NM,N,A,D,E,E2)
  205.          CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV4,RV5)
  206.          IERR(ICT) = ERROR
  207.          MBISCT = M
  208.          XLB = LB
  209.          XUB = UB
  210.          IF( ERROR .NE. 0 ) GO TO 200
  211.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  212.          IERR(ICT) = IABS(ERROR)
  213.          CALL  TRBAK1(NM,N,A,E,M,Z)
  214.          GO TO  190
  215. C
  216. C     RS1W1Z  USING  TRIDIB  AND  TINVIT
  217. C
  218.   110    ICT = 6
  219.          EPS1 = 0.0D0
  220.          CALL   TRED1(NM,N,A,D,E,E2)
  221.          CALL  TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,NO,W,IND,ERROR,RV4,RV5)
  222.          IERR(ICT) = ERROR
  223.          IF( ERROR .NE. 0 ) GO TO 200
  224.          M = NO
  225.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  226.          IERR(ICT) = IABS(ERROR)
  227.          CALL  TRBAK1(NM,N,A,E,M,Z)
  228. C
  229.   190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
  230.          CALL  RSWZR(NM,N,M,AHOLD,W,Z,RV1,RESDUL)
  231.          DFL = 10 * N
  232.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  233.          GO TO 230
  234.   200    TCRIT(ICT) = -1.0D0
  235.   230 CONTINUE
  236. C
  237.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  238.       LCOUNT = LCOUNT + 1
  239.       WRITE(IWRITE,240) N,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2), 
  240.      X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,5),
  241.      X             M11,NO,IERR(6),TCRIT(6)
  242.   240 FORMAT(I4,2(I3,2F6.3),2(1PE8.0),I3,3(I3,0PF6.3),3I3,F6.3)
  243.       GO TO  10
  244.       END
  245.       SUBROUTINE RSWZR(NM,N,M,A,W,Z,NORM,RESDUL)
  246. C
  247.       DOUBLE PRECISION NORM(M),W(M),A(NM,N),Z(NM,M),NORMA,TNORM,S,SUM,
  248.      X       SUMZ, SUMA, RESDUL
  249. C
  250. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  251. C     A*Z-Z*DIAG(W)  WHERE  A  IS A REAL SYMMETRIC MATRIX,  W  IS
  252. C     A VECTOR WHICH CONTAINS  M  EIGENVALUES OF  A, AND  Z
  253. C     IS AN ARRAY WHICH CONTAINS THE  M  CORRESPONDING EIGENVECTORS OF
  254. C     OF  A.  ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
  255. C
  256. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RSWZR).
  257. C
  258. C     INPUT.
  259. C
  260. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  261. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  262. C
  263. C        M IS THE NUMBER OF EIGENVECTORS FOR WHICH RESIDUALS ARE
  264. C           DESIRED;
  265. C
  266. C        N IS THE ORDER OF THE MATRIX  A;
  267. C
  268. C        A(NM,N) IS A REAL SYMMETRIC MATRIX.  ONLY THE FULL UPPER
  269. C           TRIANGLE NEED BE SUPPLIED;
  270. C
  271. C        Z(NM,M) IS A REAL MATRIX WHOSE FIRST  M  COLUMNS CONTAIN THE
  272. C           APPROXIMATE EIGENVECTORS OF  A;
  273. C
  274. C        W(M) IS A VECTOR WHOSE FIRST  M   COMPONENTS CONTAIN THE
  275. C           APPROXIMATE EIGENVALUES OF  A.  W(I) IS ASSOCIATED WITH THE
  276. C           I-TH COLUMN OF  Z.
  277. C
  278. C     OUTPUT.
  279. C
  280. C        Z(NM,M) IS AN ARRAY WHICH CONTAINS  M  NORMALIZED
  281. C           APPROXIMATE EIGENVECTORS OF  A.  THE EIGENVECTORS ARE
  282. C           NORMALIZED USING THE 1-NORM IN SUCH A WAY THAT THE FIRST
  283. C           ELEMENT WHOSE MAGNITUDE IS LARGER THAN THE NORM OF THE
  284. C           EIGENVECTOR DIVIDED BY  N  IS POSITIVE;
  285. C
  286. C        A(NM,N) IS ALTERED BY MAKING THE LOWER TRIANGLE OF  A
  287. C           SYMMETRIC WITH THE UPPER TRIANGLE OF  A;
  288. C
  289. C        NORM(M) IS AN ARRAY SUCH THAT FOR EACH  K,
  290. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  291. C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
  292. C
  293. C        RESDUL IS THE REAL NUMBER
  294. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  295. C
  296. C     ----------------------------------------------------------------
  297. C
  298.       RESDUL = 0.0D0
  299.       IF( M .EQ. 0 ) RETURN
  300.       NORMA = 0.0D0
  301. C
  302.       DO 40 I=1,N
  303.          SUMA = 0.0D0
  304.          IF(I .EQ. 1) GO TO 20
  305. C
  306.          DO 10 L=2,I
  307.             A(I,L-1) = A(L-1,I)
  308.    10    CONTINUE
  309. C
  310.    20    DO 30 L=1,N
  311.    30       SUMA = SUMA + DABS(A(I,L))
  312. C
  313.    40    NORMA = DMAX1(NORMA,SUMA)
  314. C
  315.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  316. C
  317.       DO 100 I=1,M
  318.          S = 0.0D0
  319.          SUMZ = 0.0D0
  320. C
  321.          DO 60 L=1,N
  322.             SUMZ = SUMZ + DABS(Z(L,I))
  323.             SUM = -W(I)*Z(L,I)
  324. C
  325.             DO 50 K=1,N
  326.    50          SUM = SUM + A(L,K)*Z(K,I)
  327. C
  328.    60       S = S + DABS(SUM)
  329. C
  330.          NORM(I) = SUMZ
  331.          IF( SUMZ .EQ. 0.0D0 )  GO TO  100
  332. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  333. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  334. C                  LARGER THAN !!Z(I)!!/N..........
  335.          DO 70 L=1,N
  336.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 80
  337.    70       CONTINUE
  338. C
  339.    80    TNORM = DSIGN(NORM(I),Z(L,I))
  340. C
  341.          DO 90 L=1,N
  342.    90       Z(L,I) = Z(L,I)/TNORM
  343. C
  344.          NORM(I) = S/(NORM(I)*NORMA)
  345.   100    RESDUL = DMAX1(NORM(I), RESDUL)
  346. C
  347.       RETURN
  348.       END
  349.       SUBROUTINE RMATIN(NM,N,A,AHOLD,INITIL)
  350. C
  351. C     THIS INPUT SUBROUTINE READS A REAL MATRIX FROM SYSIN OF
  352. C     ORDER N.
  353. C     TO GENERATE THE MATRIX  A  INITIALLY,  INITIL  IS TO BE 0.
  354. C     TO REGENERATE THE MATRIX  A  FOR THE PURPOSE OF THE RESIDUAL
  355. C     CALCULATION,  INITIL  IS TO BE  1.
  356. C
  357. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RGREADI).
  358. C
  359.       DOUBLE PRECISION A(NM,NM),AHOLD(NM,NM)
  360.       INTEGER  IA( 20)
  361.       DATA IREADA/1/,IWRITE/6/
  362. C
  363.       IF( INITIL .EQ. 1 )  GO TO  30
  364.       READ(IREADA,5)N, M
  365.     5 FORMAT(I6,6X,I6)
  366.       IF( N .EQ. 0 )  GO TO  70
  367.       IF (M .NE. 1) GO TO 14
  368.       DO  13  I = 1,N
  369.          READ(IREADA,16) (IA(J), J=I,N)
  370.          DO  12  J = I,N
  371.            A(I,J) = IA(J)
  372.    12      A(J,I) = A(I,J)
  373.    13 CONTINUE
  374.       GO TO 19
  375.    14 DO  17  I = 1,N
  376.          READ(IREADA,16) (IA(J), J=1,N)
  377.    16    FORMAT(6I12)
  378.          DO  17  J = 1,N
  379.    17      A(I,J) = IA(J)
  380.    19 DO  20  I = 1,N
  381.          DO  20  J = 1,N
  382.    20      AHOLD(I,J) = A(I,J)
  383.       RETURN
  384.    30 DO  40  I = 1,N
  385.          DO  40  J = 1,N
  386.    40      A(I,J) = AHOLD(I,J)
  387.       RETURN
  388.    70 WRITE(IWRITE,80)
  389.    80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RGREADI). /1H1)
  390.       STOP
  391.       END
  392.