qsf.f 4.65 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157
      MODULE m_qsf
      CONTAINS
      SUBROUTINE qsf(h,y,z,ndim,isave)
c     ..................................................................
c     subroutine qsf
c     purpose
c        to integrate an equidistant table of function values.
c     description of parameters
c         h     - the increment of arguement values.
c         y     - the input vector of function values.  y must be
c                 dimensioned ndim.
c         z     - if isave = 0,z is the final value of the integration.
c                 if isave not 0,z is the resulting vector of integrals
c                 for each y value. in this case it must be dimensioned
c                 ndim.
c         ndim  - the number of points in the table.
c         isave - determines whether a single value or multiple results
c                 will be returned
c                  isave =  0  - only the final value of the integra-
c                                tion will be returned.
c                  isave ne 0  - the results of integration at each
c                                point in the table will be returned.
c     remarks
c        no action in case ndim less than 3.
c     subroutines and function subprograms required
c        none
c     method
c        integration is done by means of simpsons rule together with
c        newtons 3/8 rule or a combination of these two rules. trunca-
c        tion error is of order h**5 if ndim greater than 3. if ndim = 3
c        truncation error is of order h**4.
c     reference
c        (1) f.b.hildebrand, introduction to numerical analysis,
c            mcgraw-hill, new york/toronto/london, 1956, pp.71-76.
c        (2) r. zurmuehl, praktische mathematic fuer ingenieure und
c            physiker, springer, berlin/goettingen/heidelberg, 1963,
c            pp. 214-221.
c     ..................................................................
C     .. Scalar Arguments ..

      IMPLICIT NONE
      REAL h
      INTEGER isave,ndim
C     ..
C     .. Array Arguments ..
      REAL y(ndim),z(*)
C     ..
C     .. Local Scalars ..
      REAL aux,aux1,aux2,ht,sum1,sum2,val
      INTEGER i
      LOGICAL save
C     ..
      save = isave .NE. 0
      ht = h / 3.0
      IF ( ndim == 5 ) THEN
        GOTO 140
      ELSEIF ( ndim < 5 ) THEN
        GOTO 130
      ENDIF
c     ndim is greater than 5. preparations of integration loop
      sum1 = y(2) + y(2)
      sum1 = sum1 + sum1
      sum1 = ht* (y(1)+sum1+y(3))
      aux1 = y(4) + y(4)
      aux1 = aux1 + aux1
      aux1 = sum1 + ht* (y(3)+aux1+y(5))
      aux2 = ht* (y(1)+3.875e0* (y(2)+y(5))+2.625e0* (y(3)+y(4))+y(6))
      IF (.NOT.save) GO TO 30
   20 z(1) = 0.e0
      sum2 = y(5) + y(5)
      sum2 = sum2 + sum2
      sum2 = aux2 - ht* (y(4)+sum2+y(6))
      aux = y(3) + y(3)
      aux = aux + aux
      z(2) = sum2 - ht* (y(2)+aux+y(4))
      z(3) = sum1
      z(4) = sum2
   30 CONTINUE
      IF ( ndim <= 6 ) THEN
         GOTO 70
      ENDIF
c     integration loop
   40 DO 60 i = 7,ndim,2
         sum1 = aux1
         sum2 = aux2
         aux1 = y(i-1) + y(i-1)
         aux1 = aux1 + aux1
         aux1 = sum1 + ht* (y(i-2)+aux1+y(i))
         IF (save) z(i-2) = sum1
         IF ( i >= ndim ) THEN
           GOTO 100
         ENDIF
         aux2 = y(i) + y(i)
         aux2 = aux2 + aux2
         aux2 = sum2 + ht* (y(i-1)+aux2+y(i+1))
         IF (save) z(i-1) = sum2
   60 CONTINUE
   70 IF (.NOT.save) GO TO 90
   80 z(ndim-1) = aux1
      z(ndim) = aux2
      RETURN
   90 z(1) = aux2
      RETURN
  100 IF (.NOT.save) GO TO 120
  110 z(ndim-1) = sum2
      z(ndim) = aux1
      RETURN
  120 z(1) = aux1
      RETURN
c     end of integration loop
  130 CONTINUE
      IF ( ndim == 3 ) THEN
        GOTO 210
      ELSEIF ( ndim < 3 ) THEN
        GOTO 230
      ENDIF
c     ndim is equal to 4 or 5
  140 sum2 = 1.125e0*ht* (y(1)+y(2)+y(2)+y(2)+y(3)+y(3)+y(3)+y(4))
      sum1 = y(2) + y(2)
      sum1 = sum1 + sum1
      sum1 = ht* (y(1)+sum1+y(3))
      IF (.NOT.save) GO TO 160
  150 z(1) = 0.e0
      aux1 = y(3) + y(3)
      aux1 = aux1 + aux1
      z(2) = sum2 - ht* (y(2)+aux1+y(4))
  160 CONTINUE
      IF ( ndim < 5 ) THEN
        GOTO 190
      ENDIF
  170 aux1 = y(4) + y(4)
      aux1 = aux1 + aux1
      val = sum1 + ht* (y(3)+aux1+y(5))
      IF (save) GO TO 180
      z(1) = val
      RETURN
  180 z(5) = val
      GO TO 200
  190 IF (save) GO TO 200
      z(1) = sum2
      RETURN
  200 z(3) = sum1
      z(4) = sum2
      RETURN
c     ndim is equal to 3
  210 sum2 = y(2) + y(2)
      sum2 = sum2 + sum2
      val = ht* (y(1)+sum2+y(3))
      IF (save) GO TO 220
      z(1) = val
      RETURN
  220 z(1) = 0.e0
      z(2) = ht* (1.25e0*y(1)+y(2)+y(2)-.25e0*y(3))
      z(3) = val
  230 RETURN
      END SUBROUTINE qsf
      END MODULE m_qsf