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