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 / rttest.f < prev    next >
Text File  |  1996-09-28  |  14KB  |  395 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL TRIDIAGONAL
  3. C     MATRICES SUMMARIZING THE FIGURES OF MERIT FOR ALL PATHS.
  4. C
  5. C     THIS DRIVER IS CATALOGUED AS  EISPDRV4(RTSUMARY).
  6. C
  7. C     THE DIMENSION OF  A  SHOULD BE  NM  BY  3.
  8. C     THE DIMENSION OF  Z  SHOULD BE  NM  BY  NM.
  9. C     THE DIMENSION OF  W,D,E,E2,IND,RV1,RV2,RV3,RV4,RV5,RV6,
  10. C     W1,  AND  W2  SHOULD BE  NM.
  11. C     HERE NM = 20.
  12. C
  13.       DOUBLE PRECISION Z( 20, 20),A( 20, 3),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,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='FILE39')
  22.       OPEN(UNIT=IREADC,FILE='FILE40')
  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.      X34H  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.      X30H  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   FIGI2-TQL2  .     /
  62.      X62H THE  TQLRAT  PATH USES THE EISPACK CODES   FIGI -TQLRAT.     /
  63.      X62H THE  IMTQL2  PATH USES THE EISPACK CODES   FIGI2-IMTQL2,     /
  64.      X38H AS CALLED FROM DRIVER SUBROUTINE  RT. /
  65.      X62H THE  IMTQL1  PATH USES THE EISPACK CODES   FIGI -IMTQL1,     /
  66.      X38H AS CALLED FROM DRIVER SUBROUTINE  RT. /
  67.      X63H THE  IMTQLV  PATH USES THE EISPACK CODES   FIGI -IMTQLV-TINVIT
  68.      X ,8H-BAKVEC.)
  69.       WRITE(IWRITE,3)
  70.     3 FORMAT(
  71.      X64H THE  TSTURM  PATH USES THE EISPACK CODES   FIGI -TSTURM-BAKVEC
  72.      X.  /
  73.      X63H THE  BISECT  PATH USES THE EISPACK CODES   FIGI -BISECT-TINVIT
  74.      X ,8H-BAKVEC. /
  75.      X63H THE  TRIDIB  PATH USES THE EISPACK CODES   FIGI -TRIDIB-TINVIT
  76.      X ,8H-BAKVEC. /)
  77.       WRITE(IWRITE,15)
  78.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  79.     5 FORMAT( 53H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ,
  80.      X    31HMEASURE OF PERFORMANCE  Y  FOR /5X,
  81.      X    56HTHE  EISPACK  CODES.  THIS RUN DISPLAYS THESE STATISTICS ,
  82.      X    33H FOR REAL TRIDIAGONAL MATRICES.      /
  83.      X    55H0ORDER TQL2   TQLRAT IMTQL2 IMTQL1    LB      UB    M   ,
  84.      X    40HIMTQLV   TSTURM   BISECT  M1 NO  TRIDIB )
  85.    10 CALL RMATIN(NM,N,A)
  86.       READ(IREADC,50) MM,LB,UB,M11,NO
  87.    50 FORMAT(I4,2D24.16,2(4X,I4))
  88. C
  89. C     MM,LB,UB,M11,  AND  NO  ARE READ FROM SYSIN AFTER THE MATRIX IS
  90. C     GENERATED.  MM,LB,  AND  UB  SPECIFY TO  BISECT  THE MAXIMUM
  91. C     NUMBER OF EIGENVALUES AND BOUNDS FOR THE INTERVAL WHICH IF TO
  92. C     BE SEARCHED.  M11  AND  NO  SPECIFY TO  TRIDIB  THE LOWER BOUNDARY
  93. C     INDEX AND THE NUMBER OF DESIRED EIGENVALUES.
  94. C
  95.       DO  230  ICALL = 1,10
  96. C
  97. C     IF  TQLRAT  PATH (LABEL 80) IS TAKEN THEN  TQL2  PATH (LABEL 70)
  98. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  99. C     MEANINGFUL.
  100. C     IF  IMTQL1  PATH (LABEL 85) IS TAKEN THEN  IMTQL2  PATH (LABEL 75)
  101. C     MUST ALSO BE TAKEN IN ORDER THAT THE MEASURE OF PERFORMANCE BE
  102. C     MEANINGFUL.
  103. C     IF  TQL2  (IMTQL2)  PATH FAILS, THEN  TQLRAT  (IMTQL1)  PATH IS
  104. C     OMITTED AND PRINTOUT FLAGGED WITH  -1.0.
  105. C
  106.          GO TO  (70,75,80,85,89,90,95,230,110,230),  ICALL
  107. C
  108. C     RTWZ  USING  TQL2
  109. C
  110.    70    ICT = 1
  111.          CALL  FIGI2(NM,N,A,W,E,Z,ERROR)
  112.          IERR(ICT) = ERROR
  113.          IF( ERROR .NE. 0 ) GO TO 200
  114.          CALL  TQL2(NM,N,W,E,Z,ERROR)
  115.          IERR(ICT) = ERROR
  116.          M = ERROR - 1
  117.          IF( ERROR .NE. 0 ) GO TO 190
  118.          DO 71 I = 1,N
  119.             W1(I) = W(I)
  120.    71    CONTINUE
  121.          M = N
  122.          GO TO  190
  123. C
  124. C     RTWZ  USING  IMTQL2
  125. C     INVOKED FROM DRIVER SUBROUTINE  RT.
  126. C
  127.    75    ICT = 2
  128.          CALL  RT(NM,N,A,W,1,Z,E,ERROR)
  129.          IERR(ICT) = ERROR
  130.          IF( ERROR .GT. N )  GO TO  200
  131.          M = ERROR - 1
  132.          IF( ERROR .NE. 0 ) GO TO 190
  133.          DO 78 I = 1,N
  134.             W2(I) = W(I)
  135.    78    CONTINUE
  136.          M = N
  137.          GO TO  190
  138. C
  139. C     RTW  USING  TQLRAT
  140. C
  141.    80    ICT = 7
  142.          IF( IERR(1) .NE. 0 ) GO TO 200
  143.          CALL  FIGI(NM,N,A,W,E,E2,ERROR)
  144.          CALL TQLRAT(N,W,E2,ERROR)
  145.          IERR(1) = ERROR
  146.          IF( ERROR .NE. 0 ) GO TO 200
  147.          MAXEIG = 0.0D0
  148.          MAXDIF = 0.0D0
  149.          DO 81 I = 1,N
  150.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  151.             U = DABS(W1(I) - W(I))
  152.             IF( U .GT. MAXDIF ) MAXDIF = U
  153.    81    CONTINUE
  154.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  155.          DFL = 10 * N
  156.          TCRIT(7) = MAXDIF/EPSLON(MAXEIG*DFL)
  157.          GO TO  230
  158. C
  159. C     RTW  USING  IMTQL1
  160. C     INVOKED FROM DRIVER SUBROUTINE  RT.
  161. C
  162.    85    ICT = 8
  163.          IF( IERR(2) .NE. 0 ) GO TO 200
  164.          CALL  RT(NM,N,A,W,0,Z,E,ERROR)
  165.          IERR(2) = ERROR
  166.          MAXEIG = 0.0D0
  167.          MAXDIF = 0.0D0
  168.          DO 86 I = 1,N
  169.             IF( DABS(W(I)) .GT. MAXEIG ) MAXEIG = DABS(W(I))
  170.             U = DABS(W2(I) - W(I))
  171.             IF( U .GT. MAXDIF ) MAXDIF = U
  172.    86    CONTINUE
  173.          IF( MAXEIG .EQ. 0.0D0 ) MAXEIG = 1.0D0
  174.          DFL = 10 * N
  175.          TCRIT(8) = MAXDIF/EPSLON(MAXEIG*DFL)
  176.          GO TO  230
  177. C
  178. C     RTW1Z  USING  ( USAGE HERE COMPUTES ALL THE EIGENVALUES )
  179. C
  180.    89    ICT = 3
  181.          CALL  FIGI(NM,N,A,D,E,E2,ERROR)
  182.          IERR(ICT) = IABS(ERROR)
  183.          IF( ERROR .NE. 0 ) GO TO 200
  184.          CALL IMTQLV(N,D,E,E2,W,IND,ERROR,RV1)
  185.          IERR(ICT) = ERROR
  186.          M = N
  187.          IF( ERROR .NE. 0 ) M = ERROR - 1
  188.          CALL TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  189.          IERR(ICT) = IERR(ICT) + IABS(ERROR)
  190.          CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
  191.          IERR(ICT) = IERR(ICT) + ERROR
  192.          GO TO 190
  193. C
  194. C     RT1W1Z  USING  TSTURM
  195. C
  196.    90    ICT = 4
  197.          EPS1 = 0.0D0
  198.          CALL  FIGI(NM,N,A,D,E,E2,ERROR)
  199.          IERR(ICT) = IABS(ERROR)
  200.          IF( ERROR .NE. 0 )  GO TO  200
  201.          CALL  TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,ERROR,
  202.      X                RV1,RV2,RV3,RV4,RV5,RV6)
  203.          IERR(ICT) = ERROR
  204.          XLB = LB
  205.          XUB = UB
  206.          IERR(ICT) = ERROR
  207.          IF( ERROR .EQ. 3*N + 1 ) GO TO 200
  208.          IF( ERROR .GT. 4*N ) M = ERROR - 4*N - 1
  209.          CALL  BAKVEC(NM,N,A,E,M,Z,ERROR)
  210.          IERR(ICT) = IERR(ICT) + ERROR
  211.          GO TO  190
  212. C
  213. C     RT1W1Z  USING  BISECT  AND  TINVIT
  214. C
  215.    95    ICT = 5
  216.          EPS1 = 0.0D0
  217.          CALL  FIGI(NM,N,A,D,E,E2,ERROR)
  218.          IERR(ICT) = IABS(ERROR)
  219.          IF( ERROR .NE. 0 )  GO TO  200
  220.          CALL  BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,ERROR,RV4,RV5)
  221.          IERR(ICT) = ERROR
  222.          MBISCT = M
  223.          XLB = LB
  224.          XUB = UB
  225.          IF( ERROR .NE. 0 ) GO TO 200
  226.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  227.          IERR(ICT) = IABS(ERROR)
  228.          CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
  229.          IERR(ICT) = IERR(ICT) + ERROR
  230.          GO TO  190
  231. C
  232. C     RT1W1Z  USING  TRIDIB  AND  TINVIT
  233. C
  234.   110    ICT = 6
  235.          EPS1 = 0.0D0
  236.          CALL  FIGI(NM,N,A,D,E,E2,ERROR)
  237.          IERR(ICT) = IABS(ERROR)
  238.          IF( ERROR .NE. 0 )  GO TO  200
  239.          CALL TRIDIB(N,EPS1,D,E,E2,LB,UB,M11,NO,W,IND,ERROR,RV4,RV5)
  240.          IERR(ICT) = ERROR
  241.          IF( ERROR .NE. 0 )  GO TO  200
  242.          M = NO
  243.          CALL  TINVIT(NM,N,D,E,E2,M,W,IND,Z,ERROR,RV1,RV2,RV3,RV4,RV6)
  244.          IERR(ICT) = IABS(ERROR)
  245.          CALL BAKVEC(NM,N,A,E,M,Z,ERROR)
  246.          IERR(ICT) = IERR(ICT) + ERROR
  247. C
  248.   190    IF( M .EQ. 0 .AND. ERROR .NE. 0 ) GO TO 200
  249.          CALL  RTWZR(NM,N,M,A,W,Z,RV1,RESDUL)
  250.          DFL = 10 * N
  251.          TCRIT(ICT) = RESDUL/EPSLON(DFL)
  252.          GO TO 230
  253.   200    TCRIT(ICT) = -1.0D0
  254.   230 CONTINUE
  255. C
  256.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  257.       LCOUNT = LCOUNT + 1
  258.       WRITE(IWRITE,240) N,IERR(1),TCRIT(1),TCRIT(7),IERR(2),TCRIT(2), 
  259.      X             TCRIT(8),XLB,XUB,MBISCT,(IERR(I),TCRIT(I),I=3,5),
  260.      X             M11,NO,IERR(6),TCRIT(6)
  261.   240 FORMAT(I4,2(I3,2F6.3),2(1PE8.0),I3,3(I3,0PF6.3),3I3,F6.3)
  262.       GO TO  10
  263.       END
  264.       SUBROUTINE RTWZR(NM,N,M,A,W,Z,NORM,RESDUL)
  265. C
  266.       DOUBLE PRECISION NORM(M), W(M), A(NM,3), Z(NM,M), NORMA, TNORM, 
  267.      X       S, SUM, SUMA, SUMZ, RESDUL
  268. C
  269. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  270. C     A*Z-Z*DIAG(W)  WHERE  A  IS A REAL NON-SYMMETRIC TRIDIAGONAL
  271. C     MATRIX,  W    IS A VECTOR WHICH CONTAINS  M  EIGENVALUES OF  A,
  272. C     AND  Z  IS AN ARRAY WHICH CONTAINS THE  M  CORRESPONDING
  273. C     EIGENVECTORS OF  A.  ALL NORMS APPEARING IN THE COMMENTS BELOW
  274. C     ARE 1-NORMS.
  275. C
  276. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RTWZR).
  277. C
  278. C     INPUT.
  279. C
  280. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  281. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  282. C
  283. C        N IS THE ORDER OF THE MATRIX  A;
  284. C
  285. C        M IS THE NUMBERS OF EIGENVECTORS WHOSE RESIDUALS ARE DESIRED;
  286. C
  287. C        A(N,3) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
  288. C           SUBDIAGONAL,DIAGONAL AND SUPERDIAGONAL OF THE SYMMETRIC
  289. C           TRIDIAGONAL MATRIX.  THE SUBDIAGONAL BEGINS AT  A(2,1), THE
  290. C           SUPERDIAGONAL BEGINS AT  A(1,3)  AND  A(1,1)  AND  A(N,3)
  291. C           ARE ARBITRARY;
  292. C
  293. C        W(M) IS A VECTOR WHOSE FIRST  M  COMPONENTS CONTAIN EIGENVALUES
  294. C           OF  A;
  295. C
  296. C        Z(NM,M) IS AN ARRAY WHOSE FIRST  M  COLUMNS CONTAIN THE
  297. C           EIGENVECTORS OF  A  CORRESPONDING TO THE EIGENVALUES IN  W.
  298. C
  299. C     OUTPUT.
  300. C
  301. C        Z(NM,M) IS AN ARRAY WHICH CONTAINS THE NORMALIZED
  302. C           APPROXIMATE EIGENVECTORS OF  A.  THE EIGENVECTORS
  303. C           ARE NORMALIZED USING THE 1-NORM IN SUCH A WAY
  304. C           THAT THE FIRST ELEMENT WHOSE MAGNITUDE IS LARGER
  305. C           THAN THE NORM OF THE EIGENVECTOR DIVIDED BY  N  IS
  306. C           POSITIVE;
  307. C
  308. C        NORM(M) IS AN ARRAY SUCH THAT FOR EACH  K
  309. C           NORM(K) = !!A*Z(K)-Z(K)*W(K)!!/(!!A!!*!!Z(K)!!)
  310. C           WHERE  Z(K)  IS THE K-TH EIGENVECTOR;
  311. C
  312. C        RESDUL IS THE REAL NUMBER
  313. C           !!A*Z-Z*DIAG(W)!!/(!!A!!*!!Z!!).
  314. C
  315. C     ----------------------------------------------------------------
  316. C
  317.       RESDUL = 0.0D0
  318.       IF( M .EQ. 0 ) RETURN
  319.       NORMA = 0.0D0
  320. C
  321.       DO 40 I=1,N
  322.          J = MIN0(N,I+1)
  323.          SUMA = 0.0D0
  324.          LSTART = MAX0(1,I+2-N)
  325.          LSTOP = MIN0(3,I+1)
  326. C
  327.          DO 10 L=LSTART,LSTOP
  328.             SUMA = SUMA + DABS(A(J,L))
  329.    10       J = J-1
  330. C
  331.    40    NORMA = DMAX1(SUMA,NORMA)
  332. C
  333.       IF(NORMA .EQ. 0.0D0) NORMA = 1.0D0
  334. C
  335.       DO 120 I=1,M
  336.          S = 0.0D0
  337.          SUMZ = 0.0D0
  338. C
  339.          DO 80 L=1,N
  340.             SUM = -W(I)*Z(L,I)
  341.             SUMZ = SUMZ + DABS(Z(L,I))
  342.             J = MAX0(0,L-2)
  343.             KSTOP = MIN0(3,N+2-L)
  344.             KSTART = MAX0(1,3-L)
  345. C
  346.             DO 50 K=KSTART,KSTOP
  347.                J = J+1
  348.    50          SUM = SUM + A(L,K)*Z(J,I)
  349. C
  350.    80       S = S + DABS(SUM)
  351. C
  352.          NORM(I) = SUMZ
  353.          IF( SUMZ .EQ. 0.0D0 )  GO TO  120
  354. C        ..........THIS LOOP WILL NEVER BE COMPLETED SINCE THERE
  355. C                  WILL ALWAYS EXIST AN ELEMENT IN THE VECTOR Z(I)
  356. C                  LARGER THAN !!Z(I)!!/N..........
  357.          DO 90 L=1,N
  358.             IF(DABS(Z(L,I)) .GE. NORM(I)/N) GO TO 100
  359.    90       CONTINUE
  360. C
  361.   100    TNORM = DSIGN(NORM(I),Z(L,I))
  362. C
  363.          DO 110 L=1,N
  364.   110       Z(L,I) = Z(L,I)/TNORM
  365. C
  366.          NORM(I) = S/(NORM(I)*NORMA)
  367.   120    RESDUL = DMAX1(NORM(I),RESDUL)
  368. C
  369.       RETURN
  370.       END
  371.       SUBROUTINE RMATIN(NM,N,A)
  372. C
  373. C     THIS INPUT SUBROUTINE READS A REAL TRIDIAGONAL MATRIX
  374. C     FROM SYSIN OF ORDER N.
  375. C
  376. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RTREADI).
  377. C
  378.       DOUBLE PRECISION A(NM,3)
  379.       INTEGER  IA(3)
  380.       DATA IREADA/1/,IWRITE/6/
  381. C
  382.       READ(IREADA,5) N
  383.     5 FORMAT(I6)
  384.       IF( N .EQ. 0 )  GO TO  70
  385.       DO  15  I = 1,N
  386.          READ(IREADA,10) (IA(J),J=1,3)
  387.    10    FORMAT(3I12)
  388.          DO  15  J = 1,3
  389.    15      A(I,J) = IA(J)
  390.       RETURN
  391.    70 WRITE(IWRITE,80)
  392.    80 FORMAT(46H0END OF DATA FOR SUBROUTINE  RMATIN(RTREADI). /1H1)
  393.       STOP
  394.       END
  395.