LCOV - code coverage report
Current view: top level - math - qsf.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 42 89 47.2 %
Date: 2019-09-08 04:53:50 Functions: 1 2 50.0 %

          Line data    Source code
       1             :       MODULE m_qsf
       2             : 
       3             : !      interface qsf
       4             : !        module procedure qsfReal, qsfComplex
       5             : !      end interface
       6             : 
       7             :       CONTAINS
       8             :       
       9        2278 :       SUBROUTINE qsf(h,y,z,ndim,isave)
      10             : c     ..................................................................
      11             : c     subroutine qsf
      12             : c     purpose
      13             : c        to integrate an equidistant table of function values.
      14             : c     description of parameters
      15             : c         h     - the increment of arguement values.
      16             : c         y     - the input vector of function values.  y must be
      17             : c                 dimensioned ndim.
      18             : c         z     - if isave = 0,z is the final value of the integration.
      19             : c                 if isave not 0,z is the resulting vector of integrals
      20             : c                 for each y value. in this case it must be dimensioned
      21             : c                 ndim.
      22             : c         ndim  - the number of points in the table.
      23             : c         isave - determines whether a single value or multiple results
      24             : c                 will be returned
      25             : c                  isave =  0  - only the final value of the integra-
      26             : c                                tion will be returned.
      27             : c                  isave ne 0  - the results of integration at each
      28             : c                                point in the table will be returned.
      29             : c     remarks
      30             : c        no action in case ndim less than 3.
      31             : c     subroutines and function subprograms required
      32             : c        none
      33             : c     method
      34             : c        integration is done by means of simpsons rule together with
      35             : c        newtons 3/8 rule or a combination of these two rules. trunca-
      36             : c        tion error is of order h**5 if ndim greater than 3. if ndim = 3
      37             : c        truncation error is of order h**4.
      38             : c     reference
      39             : c        (1) f.b.hildebrand, introduction to numerical analysis,
      40             : c            mcgraw-hill, new york/toronto/london, 1956, pp.71-76.
      41             : c        (2) r. zurmuehl, praktische mathematic fuer ingenieure und
      42             : c            physiker, springer, berlin/goettingen/heidelberg, 1963,
      43             : c            pp. 214-221.
      44             : c     ..................................................................
      45             : C     .. Scalar Arguments ..
      46             : 
      47             :       IMPLICIT NONE
      48             :       REAL h
      49             :       INTEGER isave,ndim
      50             : C     ..
      51             : C     .. Array Arguments ..
      52             :       REAL y(ndim),z(*)
      53             : C     ..
      54             : C     .. Local Scalars ..
      55             :       REAL aux,aux1,aux2,ht,sum1,sum2,val
      56             :       INTEGER i
      57             :       LOGICAL save
      58             : C     ..
      59        2278 :       save = isave .NE. 0
      60        2278 :       ht = h / 3.0
      61        2278 :       IF ( ndim == 5 ) THEN
      62             :         GOTO 140
      63        2278 :       ELSEIF ( ndim < 5 ) THEN
      64             :         GOTO 130
      65             :       ENDIF
      66             : c     ndim is greater than 5. preparations of integration loop
      67        2278 :       sum1 = y(2) + y(2)
      68        2278 :       sum1 = sum1 + sum1
      69        2278 :       sum1 = ht* (y(1)+sum1+y(3))
      70        2278 :       aux1 = y(4) + y(4)
      71        2278 :       aux1 = aux1 + aux1
      72        2278 :       aux1 = sum1 + ht* (y(3)+aux1+y(5))
      73        2278 :       aux2 = ht* (y(1)+3.875e0* (y(2)+y(5))+2.625e0* (y(3)+y(4))+y(6))
      74        2278 :       IF (.NOT.save) GO TO 30
      75        2250 :    20 z(1) = 0.e0
      76        2250 :       sum2 = y(5) + y(5)
      77        2250 :       sum2 = sum2 + sum2
      78        2250 :       sum2 = aux2 - ht* (y(4)+sum2+y(6))
      79        2250 :       aux = y(3) + y(3)
      80        2250 :       aux = aux + aux
      81        2250 :       z(2) = sum2 - ht* (y(2)+aux+y(4))
      82        2250 :       z(3) = sum1
      83        2278 :       z(4) = sum2
      84             :    30 CONTINUE
      85        2278 :       IF ( ndim <= 6 ) THEN
      86             :          GOTO 70
      87             :       ENDIF
      88             : c     integration loop
      89      112194 :    40 DO 60 i = 7,ndim,2
      90      109916 :          sum1 = aux1
      91      109916 :          sum2 = aux2
      92      109916 :          aux1 = y(i-1) + y(i-1)
      93      109916 :          aux1 = aux1 + aux1
      94      109916 :          aux1 = sum1 + ht* (y(i-2)+aux1+y(i))
      95      109916 :          IF (save) z(i-2) = sum1
      96      109916 :          IF ( i >= ndim ) THEN
      97             :            GOTO 100
      98             :          ENDIF
      99      109916 :          aux2 = y(i) + y(i)
     100      109916 :          aux2 = aux2 + aux2
     101      109916 :          aux2 = sum2 + ht* (y(i-1)+aux2+y(i+1))
     102      109916 :          IF (save) z(i-1) = sum2
     103        2278 :    60 CONTINUE
     104        2278 :    70 IF (.NOT.save) GO TO 90
     105        2250 :    80 z(ndim-1) = aux1
     106        2250 :       z(ndim) = aux2
     107        2250 :       RETURN
     108          28 :    90 z(1) = aux2
     109          28 :       RETURN
     110           0 :   100 IF (.NOT.save) GO TO 120
     111           0 :   110 z(ndim-1) = sum2
     112           0 :       z(ndim) = aux1
     113           0 :       RETURN
     114           0 :   120 z(1) = aux1
     115           0 :       RETURN
     116             : c     end of integration loop
     117             :   130 CONTINUE
     118           0 :       IF ( ndim == 3 ) THEN
     119             :         GOTO 210
     120           0 :       ELSEIF ( ndim < 3 ) THEN
     121             :         GOTO 230
     122             :       ENDIF
     123             : c     ndim is equal to 4 or 5
     124           0 :   140 sum2 = 1.125e0*ht* (y(1)+y(2)+y(2)+y(2)+y(3)+y(3)+y(3)+y(4))
     125           0 :       sum1 = y(2) + y(2)
     126           0 :       sum1 = sum1 + sum1
     127           0 :       sum1 = ht* (y(1)+sum1+y(3))
     128           0 :       IF (.NOT.save) GO TO 160
     129           0 :   150 z(1) = 0.e0
     130           0 :       aux1 = y(3) + y(3)
     131           0 :       aux1 = aux1 + aux1
     132           0 :       z(2) = sum2 - ht* (y(2)+aux1+y(4))
     133             :   160 CONTINUE
     134           0 :       IF ( ndim < 5 ) THEN
     135             :         GOTO 190
     136             :       ENDIF
     137           0 :   170 aux1 = y(4) + y(4)
     138           0 :       aux1 = aux1 + aux1
     139           0 :       val = sum1 + ht* (y(3)+aux1+y(5))
     140           0 :       IF (save) GO TO 180
     141           0 :       z(1) = val
     142           0 :       RETURN
     143           0 :   180 z(5) = val
     144           0 :       GO TO 200
     145           0 :   190 IF (save) GO TO 200
     146           0 :       z(1) = sum2
     147           0 :       RETURN
     148           0 :   200 z(3) = sum1
     149           0 :       z(4) = sum2
     150           0 :       RETURN
     151             : c     ndim is equal to 3
     152           0 :   210 sum2 = y(2) + y(2)
     153           0 :       sum2 = sum2 + sum2
     154           0 :       val = ht* (y(1)+sum2+y(3))
     155           0 :       IF (save) GO TO 220
     156           0 :       z(1) = val
     157           0 :       RETURN
     158           0 :   220 z(1) = 0.e0
     159           0 :       z(2) = ht* (1.25e0*y(1)+y(2)+y(2)-.25e0*y(3))
     160           0 :       z(3) = val
     161             :   230 RETURN
     162             :       END SUBROUTINE qsf
     163             : 
     164           0 :       subroutine qsfComplex( h, y, z, ndim, isave )
     165             : 
     166             :         implicit none
     167             : 
     168             :         REAL h
     169             :         INTEGER isave,ndim
     170             :         COMPLEX y(ndim),z(*)
     171             : 
     172           0 :         real zr(ndim), zi(ndim)
     173             : 
     174           0 :         call qsf( h,  real(y), zr, ndim, isave )
     175           0 :         call qsf( h, aimag(y), zi, ndim, isave )
     176           0 :         z(1:ndim) = cmplx( zr(1:ndim), zi(1:ndim) )
     177             : 
     178           0 :       end subroutine qsfComplex
     179             : 
     180             :       END MODULE m_qsf

Generated by: LCOV version 1.13