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 / rbltest.f < prev    next >
Text File  |  1996-09-28  |  9KB  |  295 lines
  1. C
  2. C     THIS DRIVER TESTS  EISPACK  FOR THE CLASS OF REAL
  3. C     BAND MATRICES EXHIBITING THE USE OF  EISPACK  TO FIND THE
  4. C     SOLUTION TO THE EQUATION  A*X = B .  THIS DRIVER SUMMARIZES
  5. C     THE FIGURES OF MERIT FOR THE PATH.
  6. C
  7. C     THIS DRIVER IS CATALOGUED AS EISPDRV4(RBLSUMAR).
  8. C
  9. C     THE DIMENSION OF  A  AND  AHOLD  SHOULD BE  NM  BY  MBB.
  10. C     THE DIMENSION OF  X,B,  AND  BHOLD  SHOULD BE  NM  BY  NM.
  11. C     THE DIMENSION OF  RESDUL  SHOULD BE  NM.
  12. C     THE DIMENSION OF  RV  SHOULD BE  NM*MBB.
  13. C     THE DIMENSION OF  RV1  AND  RV6  SHOULD BE  NM.
  14. C     HERE  NM = 20  AND  MBB = 39.
  15. C
  16. C
  17.       DOUBLE PRECISION A( 20, 39),X( 20, 20),B( 20, 20),RV( 780),
  18.      X        RESDUL( 20),RV1( 20),RV6( 20),TCRIT( 1),EPSLON,
  19.      X        RESMAX,E21,DFL
  20.       DOUBLE PRECISION AHOLD( 20, 39),BHOLD( 20, 20)
  21.       INTEGER IERR( 1),ERROR
  22.       DATA IREAD1/1/,IREADC/5/,IWRITE/6/
  23. C
  24.       OPEN(UNIT=IREAD1,FILE='FILE49')
  25.       OPEN(UNIT=IREADC,FILE='FILE50')
  26.       REWIND IREAD1
  27.       REWIND IREADC
  28. C
  29.       NM = 20
  30.       MBB = 39
  31.       LCOUNT = 0
  32.       WRITE(IWRITE,1)
  33.     1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S
  34.      XTATISTICS//1H ,95(1H-)/22H ORDER B.W. IP   BANDV /
  35.      X1H ,95(1H-)//
  36.      X49H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. //
  37.      X52H UNDER 'B.W.' IS THE FULL BAND WIDTH OF THE MATRIX. //
  38.      X59H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX. //
  39.      X62H UNDER 'BANDV' ARE TWO NUMBERS.  THE FIRST NUMBER, AN INTEGER,,
  40.      X29H IS THE ABSOLUTE VALUE OF THE  /
  41.      X92H ERROR FLAG RETURNED FROM  BANDV.  THE SECOND NUMBER IS THE MEA
  42.      XSURE OF PERFORMANCE BASED     /
  43.      X46H UPON THE RESIDUAL COMPUTED FROM THE SOLUTION. /)
  44.       WRITE(IWRITE,15)
  45.    15 FORMAT(1X,21HD.P. VERSION 04/15/83 )
  46.     5 FORMAT( 80H1       TABULATION OF THE ERROR FLAG  ERROR  AND THE ME
  47.      XASURE OF PERFORMANCE FOR /5X,13HTHE  EISPACK ,
  48.      X73H CODES.  THIS RUN DISPLAYS THESE STATISTICS FOR REAL BAND LINEA
  49.      XR SYSTEMS. //
  50.      X22H ORDER B.W. IP   BANDV )
  51.    10 CALL  RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,0)
  52.       CALL  RMATIN(NM,NM,N,IP,X,AHOLD,BHOLD,1,0)
  53.       READ(IREADC,20) IE21
  54.    20 FORMAT(I6)
  55.       E21 = IE21
  56. C
  57. C     RBL  USING  BANDV
  58. C
  59.       ICT = 1
  60.       DO 905 I = 1,N
  61.          RV1(I) = 0.0D0
  62.   905 CONTINUE
  63.       CALL  BANDV(NM,N,MBW,A,E21,IP,RV1,X,ERROR,780,RV,RV6)
  64.       IERR(ICT) = IABS(ERROR)
  65.       CALL  RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,1)
  66.       CALL  RMATIN(NM,NM,N,IP,B,AHOLD,BHOLD,1,1)
  67.       CALL  RBLRES(NM,N,MBW,IP,A,B,X,RESDUL,E21)
  68.       RESMAX = 0.0D0
  69.       DO 935 I = 1,IP
  70.          IF( RESDUL(I) .GT. RESMAX ) RESMAX = RESDUL(I)
  71.   935 CONTINUE
  72.       DFL = N
  73.       TCRIT(ICT) = RESMAX/EPSLON(DFL)
  74. C
  75.       IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5)
  76.       LCOUNT = LCOUNT + 1
  77.       MBF = MBW
  78.       IF( E21 .EQ. 1.0D0 ) MBF = MBW*2 - 1
  79.       WRITE(IWRITE,520) N,MBF,IP,IERR(1),TCRIT(1)
  80.   520 FORMAT(1X,I3,2I5,I4,F6.2)
  81.       GO TO 10
  82.       END
  83.       SUBROUTINE RBLRES(NM,N,MBW,IP,A,B,X,NORM,E21)
  84. C
  85.       DOUBLE PRECISION A(NM,MBW),B(NM,IP),X(NM,IP),NORM(N)
  86.       DOUBLE PRECISION NORMA,NORMX,NORMR,SUM,SUMA,E21
  87. C
  88. C     THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX
  89. C     A*X-B  WHERE  A  IS A REAL BAND MATRIX, WHICH MAY BE SYMMETRIC,
  90. C     X  IS A REAL  N  BY  IP  MATRIX, AND   B IS A REAL
  91. C     N  BY  IP  MATRIX.
  92. C
  93. C     THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RBLRES).
  94. C
  95. C     INPUT.
  96. C        NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS
  97. C           AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT;
  98. C
  99. C        N IS THE ROW DIMENSION OF THE MATRIX  A;
  100. C
  101. C        MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE
  102. C           BAND MATRIX.  IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF)
  103. C           BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT
  104. C           DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO
  105. C           SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE
  106. C           MATRIX.  IF THE MATRIX IS NON-SYMMETRIC IT MUST HAVE THE SAM
  107. C           NUMBER OF ADJACENT DIAGONALS ABOVE THE MAIN DIAGONAL AS
  108. C           BELOW.  IN THIS CASE  MBW = 2*MB-1;
  109. C
  110. C        IP  IS THE COLUMN DIMENSION OF  B  AND  X;
  111. C
  112. C        A(N,MBW) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE
  113. C           SUBDIAGONALS AND DIAGONAL OF THE SYMMETRIC BAND
  114. C           MATRIX. IF  A  IS NON-SYMMETRIC IT ALSO CONTAINS
  115. C           THE SUPER-DIAGONALS;
  116. C
  117. C        B  IS A REAL FULL  N  BY  IP  MATRIX;
  118. C
  119. C        X  IS A REAL FULL  N  BY  IP  MATRIX;
  120. C
  121. C        E21 TELLS IF THE MATRIX IS SYMMETRIC.  IF E21 EQUALS ONE THE
  122. C           MATRIX IS SYMMETRIC AND IF E21 EQUALS MINUS ONE THE
  123. C           MATRIX IS NON-SYMMETRIC.
  124. C
  125. C     OUTPUT.
  126. C
  127. C        NORM(N)  IS AN ARRAY SUCH THAT FOR EACH K,
  128. C               NORM(K) = !!A*X(K)-B(K)!!/(!!A!!*!!X(K)!!) .
  129. C
  130. C     ------------------------------------------------------------------
  131. C
  132.       IF( E21 .EQ. 1.0D0 ) GO TO 75
  133.       MBWT = MBW
  134.       MB = MBW/2 + 1
  135.       NORMA = 0.0D0
  136. C
  137.       DO 35 I = 1,N
  138.          L = MAX0(MB - I + 1,1)
  139.          M = MIN0(MB - I + N,MBWT)
  140.          SUMA = 0.0D0
  141. C
  142.          DO 34 J = L,M
  143.             SUMA = SUMA + DABS(A(I,J))
  144.    34    CONTINUE
  145. C
  146.          NORMA = DMAX1(NORMA,SUMA)
  147.    35 CONTINUE
  148. C
  149.       DO 70 IPP = 1,IP
  150.          NORMR = 0.0D0
  151.          NORMX = 0.0D0
  152.          DO 50 I = 1,N
  153.             L = MAX0(MB - I + 1,1)
  154.             M = MIN0(MB - I + N,MBWT)
  155.             K = MAX0(-MB + 1 + I,1)
  156.             SUM = -B(I,IPP)
  157.             SUMA = 0.0D0
  158. C
  159.             DO 40 J = L,M
  160.                SUM = SUM + A(I,J)*X(K,IPP)
  161.                K = K + 1
  162.                SUMA = SUMA + DABS(A(I,J))
  163.    40       CONTINUE
  164. C
  165.             NORMR = NORMR + DABS(SUM)
  166.             NORMX = NORMX + DABS(X(I,IPP))
  167.    50    CONTINUE
  168. C
  169.          IF( NORMX .EQ. 0.0D0 ) NORMX = 1.0D0
  170.          IF( NORMA .EQ. 0.0D0 ) NORMA = 1.0D0
  171.          NORM(IPP) = NORMR/(NORMA*NORMX)
  172.    70 CONTINUE
  173. C
  174.       RETURN
  175. C
  176.    75 MB = MBW
  177.       NORMA = 0.0D0
  178.       MB1 = MB - 1
  179. C
  180.       DO 110 I = 1,N
  181.          SUMA = 0.0D0
  182.          IF( I .EQ. 1 ) GO TO 90
  183.          LSTART = MAX0(1,MB + 1 - I)
  184. C
  185.          DO 80 L = LSTART,MB1
  186.             SUMA = SUMA + DABS(A(I,L))
  187.    80    CONTINUE
  188. C
  189.    90    LSTOP = MIN0(MB,N + 1 - I)
  190.          J = I
  191. C
  192.          DO 100 L = 1,LSTOP
  193.             L1 = MB + 1 - L
  194.             SUMA = SUMA + DABS(A(J,L1))
  195.             J = J + 1
  196.   100    CONTINUE
  197. C
  198.          NORMA = DMAX1(NORMA,SUMA)
  199.   110 CONTINUE
  200. C
  201.       DO 170 I = 1,IP
  202.          NORMR = 0.0D0
  203.          NORMX = 0.0D0
  204. C
  205.          DO 150 L = 1,N
  206.             SUM = -B(L,I)
  207.             J = MAX0(0,L - MB)
  208.             IF( L .EQ. 1 ) GO TO 130
  209.             KSTART = MAX0(1,MB + 1 - L)
  210. C
  211.             DO 120 K = KSTART,MB1
  212.                J = J + 1
  213.                SUM = SUM + A(L,K)*X(J,I)
  214.   120       CONTINUE
  215. C
  216.   130       KSTOP = MIN0(MB,N + 1 - L)
  217. C
  218.             DO 140 K = 1,KSTOP
  219.                J = J + 1
  220.                K1 = MB + 1 - K
  221.                SUM = SUM + A(J,K1)*X(J,I)
  222.   140       CONTINUE
  223. C
  224.             NORMR = NORMR + DABS(SUM)
  225.   150    CONTINUE
  226. C
  227.          DO 160 K = 1,N
  228.             NORMX = NORMX + DABS(X(K,I))
  229.   160    CONTINUE
  230. C
  231.          IF( NORMX .EQ. 0.0D0 ) NORMX = 1.0D0
  232.          IF( NORMA .EQ. 0.0D0 ) NORMA = 1.0D0
  233.          NORM(I) = NORMR/(NORMA*NORMX)
  234.   170 CONTINUE
  235.       RETURN
  236. C
  237.       END
  238.       SUBROUTINE RMATIN(NM,NN,M,N,A,AHOLD,BHOLD,AORB,INITIL)
  239. C
  240. C     THIS INPUT SUBROUTINE READS A REAL MATRIX, WITH DIMENSION
  241. C     M  BY  N, FROM SYSIN.
  242. C     TO GENERATE THE MATRIX  A (OR  B)  INITIALLY,  INITIL  IS
  243. C     SET TO 0.  TO REGENERATE THE MATRIX  A (OR  B)  FOR THE
  244. C     PURPOSE OF THE RESIDUAL CALCULATION,  INITIL  IS SET TO 1.
  245. C
  246. C     THIS ROUTINE IS CATALOGUED AS  EISPDRV4(RLREADI).
  247. C
  248.       INTEGER IA( 20),NM,NN,M,N,AORB,INITIL
  249.       DOUBLE PRECISION A(NM,NN),AHOLD(NM,NN),BHOLD(NM,NM)
  250.       DATA IREADA/1/,IWRITE/6/
  251. C
  252.       IF( INITIL .EQ. 1 ) GO TO 90
  253.       READ(IREADA,5) M,N
  254.     5 FORMAT(I6,6X,I6)
  255.       IF( M .EQ. 0 ) GO TO 170
  256.       DO 20 I = 1,M
  257.          READ(IREADA,10) (IA(J),J=1,N)
  258.    10    FORMAT(6I12)
  259.          DO 15 J = 1,N
  260.             A(I,J) = IA(J)
  261.    15    CONTINUE
  262.    20 CONTINUE
  263.       IF( AORB .EQ. 1 ) GO TO 50
  264.       DO 40 J = 1,N
  265.          DO 30 I = 1,M
  266.             AHOLD(I,J) = A(I,J)
  267.    30    CONTINUE
  268.    40 CONTINUE
  269.       GO TO 80
  270.    50 DO 70 J = 1,N
  271.          DO 60 I = 1,M
  272.             BHOLD(I,J) = A(I,J)
  273.    60    CONTINUE
  274.    70 CONTINUE
  275.    80 RETURN
  276. C
  277.    90 IF( AORB .EQ. 1 ) GO TO 120
  278.       DO 110 J = 1,N
  279.          DO 100 I = 1,M
  280.             A(I,J) = AHOLD(I,J)
  281.   100    CONTINUE
  282.   110 CONTINUE
  283.       GO TO 150
  284.   120 DO 140 J = 1,N
  285.          DO 130 I = 1,M
  286.             A(I,J) = BHOLD(I,J)
  287.   130    CONTINUE
  288.   140 CONTINUE
  289.   150 RETURN
  290. C
  291.   170 WRITE(IWRITE,175)
  292.   175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1)
  293.       STOP
  294.       END
  295.