home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / prgramer / adaptor / examples / purdue / prob15.fcm < prev   
Encoding:
Text File  |  1993-06-26  |  6.0 KB  |  239 lines
  1.       PROGRAM PROB15
  2. C
  3. C     PROBLEM 15
  4. C
  5. C  REFERENCE:  PROBLEMS TO TEST PARALLEL AND VECTOR LANGUAGES
  6. C              CSD-TR 516, COMPUTER SCIENCE, PURDUE UNIVERSITY
  7. C              JOHN R. RICE, MAY 1, 1985
  8. C
  9. C              REVISED BY JOHN R. RICE AND J. JING, OCT. 1, 1990
  10. C
  11. C
  12. C      *************************************************
  13. C      *      Adapted for FORTRAN D benchmarking       *
  14. C      *    by  T. HAUPT  (haupt@sccs.npac.syr.edu)    *
  15. C      *                                               *
  16. C      *    Northeast Parallel Architectures Center    *
  17. C      *   at Syracuse University, Syracuse, NY, USA   *
  18. C      *************************************************
  19. C
  20. C
  21. C       VERSION SIMD/CM2-1.00
  22. C       ==================================================
  23. C
  24.       INCLUDE '/usr/include/cm/paris-configuration-fort.h'
  25.       INTEGER KASES
  26.       PARAMETER (KASES=2)
  27.       INTEGER N(KASES)
  28. cmf$  layout N(:serial)
  29.       DATA N /8192,16384/
  30.       INTEGER NK,NPTS
  31.       PARAMETER(NPTS=10)
  32.       REAL ERRMAX(NPTS,2),DECAY(NPTS,2)
  33. cmf$  layout errmax (:serial,:serial)
  34. cmf$  layout decay (:serial,:serial)
  35.  
  36.       PRINT *, 'PROBLEM 15 started'
  37.       DO 50 K = 1, KASES
  38.  
  39.       CALL CM_TIMER_CLEAR(0)
  40.       CALL CM_TIMER_START(0)
  41.       DO MANY=1,10
  42.       NK=N(K)
  43.       CALL DOIT(NK,ERRMAX,DECAY)
  44.       ENDDO
  45.       CALL CM_TIMER_STOP(0)
  46.  
  47.       DO 170 KASE=1,2
  48.          IF (KASE.EQ.1) THEN
  49.             PRINT 110,NK
  50.   110       FORMAT (4X,'PROBLEM 15 WITH NK=',I8/)
  51.             PRINT 120
  52.   120       FORMAT (9X,'EQUALLY SPACED POINTS'/)
  53.             PRINT 130
  54.   130       FORMAT (4X,'N      MAX ERROR      DECAY EXPONENT')
  55.          ELSE
  56.             PRINT 140
  57.   140       FORMAT (/9X,'CHEBYSHEV SPACED POINTS'/)
  58.             PRINT 130
  59.          ENDIF
  60.          DO 160 I = 2, NPTS
  61.             PRINT 150,I,ERRMAX(I,KASE),DECAY(I,KASE)
  62.   150       FORMAT (1X,I4,4X,E12.4,3X,F11.2)
  63.   160    CONTINUE
  64. C
  65. C             PRINT OUT OF INTERPOLATION POINT ARRAY XPTS(J,N,KASE)
  66. C             SUPPRESSED HERE TO REDUCE AMOUNT OF OUTPUT
  67. C
  68.   170 CONTINUE
  69.  
  70.  
  71.       CALL CM_TIMER_PRINT(0)
  72.  
  73.    50 CONTINUE
  74.       STOP
  75.       END
  76.  
  77.       SUBROUTINE DOIT(NK,ERRMAX,DECAY)
  78.       INTEGER NK,NPTS
  79.       PARAMETER(NPTS=10)
  80.       REAL ERRMAX(NPTS,2),DECAY(NPTS,2)
  81. cmf$  layout errmax(:serial,:serial)
  82. cmf$  layout decay(:serial,:serial)
  83.       INTEGER NMAX,I,K,KNOTSMAX
  84.       PARAMETER (NMAX=2*NPTS,KNOTSMAX=NPTS+2)
  85.       REAL XPTS(NPTS),COEF(NMAX),H
  86. cmf$  layout xpts(:serial)
  87. cmf$  layout coef(:serial)
  88.       REAL TKNOTS(KNOTSMAX)
  89. cmf$  layout tknots(:serial)
  90.       REAL LOGERR,OLDERR,INTERP,RATIO,A,B
  91.       DATA A,B / -1.,1. /,PI / 3.141592654 /
  92.       INTEGER ERRTYP1,ERRTYP2
  93.       DATA ERRTYP1,ERRTYP2 /31,32/
  94.       REAL ERR1
  95.  
  96.       DO KASE = 1, 2
  97. C
  98. C        BEGIN             'NUMBER OF POINTS USED'
  99. C
  100.          DO 80 N = 2, NPTS
  101.  
  102. C
  103. C           BEGIN                'CASES OF POINT DISTRIBUTION'
  104. C
  105.             IF (KASE.EQ.1) THEN
  106.                H = (B-A)/(N-1)
  107.                DO 10 J = 1, N
  108.                   XPTS(J) = A+(J-1)*H
  109.    10          CONTINUE
  110.             ELSE
  111.                DO 20 J = 1, N
  112.                   XPTS(J) = (A+B)/2+(A-B)/2*COS((2*J-1)*PI/(2*N))
  113.    20          CONTINUE
  114.                RATIO = (B-A)/(XPTS(N)-XPTS(1))
  115.                DO 30 J = 1, N
  116.                   XPTS(J) = (XPTS(J)-(A+B)/2)*
  117.      *                                RATIO+(A+B)/2
  118.    30          CONTINUE
  119.             ENDIF
  120.  
  121.             DO 40 J = 1, N
  122.                call F(XPTS(J), COEF(2*J-1))
  123.    40       CONTINUE
  124.             DO 50 J = 1, N
  125.                call FPRIME (XPTS(J), COEF(2*J))
  126.    50       CONTINUE
  127.  
  128.             DO 60 J = 1, N
  129.                TKNOTS(J+1) = XPTS(J)
  130.    60       CONTINUE
  131. C
  132. C           ADD THE DUMMY POINTS AT EACH END OF TKNOTS ARRAY
  133. C
  134.             TKNOTS(1) = A-0.1*(ABS(A)+1.)
  135.             TKNOTS(N+2) = B+0.1*(ABS(B)+1.)
  136.             KNOTS = N+2
  137.  
  138.             CALL INTERPOLATE(N,NK,A,B,COEF,TKNOTS,ERRMAX(N,KASE))
  139.  
  140.  
  141.    80    CONTINUE
  142.  
  143.          DO 90 N = 3, NPTS
  144.             OLDERR = ALOG(ERRMAX(N-1,KASE))
  145.             LOGERR = ALOG(ERRMAX(N,KASE))
  146.             DECAY(N,KASE) = (LOGERR-OLDERR)/ALOG(FLOAT(N)/FLOAT(N-1))
  147.    90    CONTINUE
  148.  
  149.       ENDDO
  150.  
  151.  
  152. c     RETURN
  153.       END
  154.  
  155.       SUBROUTINE F (X,R)
  156.       REAL X, R
  157.       R = 1./(X*X+25.0)
  158.       END
  159.  
  160.       SUBROUTINE FPRIME (X,R)
  161.       REAL X, R
  162.       R = -2.0*X/(X*X+25.)**2
  163.       END
  164.  
  165.       SUBROUTINE FUN(X,NK,F)
  166.       INTEGER NK
  167.       REAL, ARRAY(NK) :: X,F
  168.       F=1./(X*X+25.0)
  169.       END
  170.  
  171.       SUBROUTINE INTERPOLATE(N,NK,A,B,COEF,T,ERROR)
  172.       INTEGER NPTS,NMAX,KNOTSMAX
  173.       PARAMETER (NPTS=10)
  174.       PARAMETER (NMAX=2*NPTS,KNOTSMAX=NPTS+2)
  175.       REAL COEF(NMAX),T(KNOTSMAX)
  176. cmf$  layout coef(:serial)
  177. cmf$  layout t(:serial)
  178.       REAL ERROR,A,B
  179.       INTEGER N,NK,J,IKNOT
  180.       REAL, ARRAY(NK) :: X,F,INTERP,DTDX,DTDX2,DX,COEF2
  181.       REAL, ARRAY(NMAX,NK) :: HCUBIC
  182. CMF$  LAYOUT HCUBIC(:SERIAL,:NEWS)
  183.  
  184.       X=A+[0:NK-1]*(B-A)/(NK-1)
  185.       CALL FUN(X,NK,F)
  186.       COEF2=0
  187.  
  188.       DO J=1,2*N
  189.        DTDX=0
  190.        DTDX2=0
  191.        DX=0
  192.  
  193.        IKNOT=(J+3)/2
  194.        IF(MOD(J,2).EQ.1) THEN
  195.  
  196.          WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT)))
  197.              DTDX=(T(IKNOT+1)-X)/(T(IKNOT+1)-T(IKNOT))
  198.          ENDWHERE
  199.  
  200.          WHERE((X.LE.T(IKNOT)) .AND. (X.GT.T(IKNOT-1)))
  201.              DTDX=(X-T(IKNOT-1))/(T(IKNOT)-T(IKNOT-1))
  202.          ENDWHERE
  203.  
  204.          WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT-1)))
  205.              HCUBIC(J,1:NK)=(3.0-2.0*DTDX(1:NK))*DTDX(1:NK)**2
  206.          ELSEWHERE
  207.            HCUBIC(J,1:NK)=0
  208.          ENDWHERE
  209.  
  210.       ELSE
  211.  
  212.          WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT)))
  213.            DX=X-T(IKNOT)
  214.            DTDX2=((X-T(IKNOT+1))/(T(IKNOT)-T(IKNOT+1)))**2
  215.          ENDWHERE
  216.  
  217.          WHERE((X.LE.T(IKNOT)) .AND. (X.GT.T(IKNOT-1)))
  218.            DX=X-T(IKNOT)
  219.            DTDX2=((X-T(IKNOT-1))/(T(IKNOT)-T(IKNOT-1)))**2
  220.  
  221.          ENDWHERE
  222.  
  223.          WHERE((X.LT.T(IKNOT+1)) .AND. (X.GT.T(IKNOT-1)))
  224.            HCUBIC(J,1:NK)=DTDX2(1:NK)*DX(1:NK)
  225.          ELSEWHERE
  226.            HCUBIC(J,1:NK)=0
  227.          ENDWHERE
  228.  
  229.        ENDIF
  230.        COEF2(1:NK)=COEF(J)*HCUBIC(J,1:NK)+COEF2(1:NK)
  231.       ENDDO
  232.  
  233. C      INTERP=SUM(COEF2,DIM=1)
  234.       INTERP=COEF2
  235.       ERROR=MAXVAL(ABS(F-INTERP))
  236.  
  237.       END
  238.  
  239.