chem/KPP/kpp/kpp-2.1/util/blas.f

   1 !--------------------------------------------------------------
   2 !
   3 ! BLAS/LAPACK-like subroutines used by the integration algorithms
   4 ! It is recommended to replace them by calls to the optimized
   5 !      BLAS/LAPACK library for your machine
   6 !
   7 !  (C) Adrian Sandu, Aug. 2004
   8 !      Virginia Polytechnic Institute and State University
   9 !--------------------------------------------------------------
  10
  11
  12 !--------------------------------------------------------------
  13       SUBROUTINE  WCOPY(N,X,incX,Y,incY)
  14 !--------------------------------------------------------------
  15 !     copies a vector, x, to a vector, y:  y <- x
  16 !     only for incX=incY=1
  17 !     after BLAS
  18 !     replace this by the function from the optimized BLAS implementation:
  19 !         CALL  SCOPY(N,X,1,Y,1)   or   CALL  DCOPY(N,X,1,Y,1)
  20 !--------------------------------------------------------------
  21
  22       INTEGER i,incX,incY,M,MP1,N
  23       KPP_REAL X(N),Y(N)
  24
  25       IF (N.LE.0) RETURN
  26
  27       M = MOD(N,8)
  28       IF( M .NE. 0 ) THEN
  29         DO i = 1,M
  30           Y(i) = X(i)
  31         END DO
  32         IF( N .LT. 8 ) RETURN
  33       END IF
  34       MP1 = M+1
  35       DO i = MP1,N,8
  36         Y(i) = X(i)
  37         Y(i + 1) = X(i + 1)
  38         Y(i + 2) = X(i + 2)
  39         Y(i + 3) = X(i + 3)
  40         Y(i + 4) = X(i + 4)
  41         Y(i + 5) = X(i + 5)
  42         Y(i + 6) = X(i + 6)
  43         Y(i + 7) = X(i + 7)
  44       END DO
  45       RETURN
  46       END
  47
  48
  49 !--------------------------------------------------------------
  50       SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
  51 !--------------------------------------------------------------
  52 !     constant times a vector plus a vector: y <- y + Alpha*x
  53 !     only for incX=incY=1
  54 !     after BLAS
  55 !     replace this by the function from the optimized BLAS implementation:
  56 !         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
  57 !--------------------------------------------------------------
  58
  59       INTEGER i,incX,incY,M,MP1,N
  60       KPP_REAL X(N),Y(N),Alpha
  61       KPP_REAL ZERO
  62       PARAMETER( ZERO = 0.0d0 )
  63
  64       IF (Alpha .EQ. ZERO) RETURN
  65       IF (N .LE. 0) RETURN
  66
  67       M = MOD(N,4)
  68       IF( M .NE. 0 ) THEN
  69         DO i = 1,M
  70           Y(i) = Y(i) + Alpha*X(i)
  71         END DO
  72         IF( N .LT. 4 ) RETURN
  73       END IF
  74       MP1 = M + 1
  75       DO i = MP1,N,4
  76         Y(i) = Y(i) + Alpha*X(i)
  77         Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
  78         Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
  79         Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
  80       END DO
  81       RETURN
  82       END
  83
  84
  85
  86 !--------------------------------------------------------------
  87       SUBROUTINE WSCAL(N,Alpha,X,incX)
  88 !--------------------------------------------------------------
  89 !     constant times a vector: x(1:N) <- Alpha*x(1:N)
  90 !     only for incX=incY=1
  91 !     after BLAS
  92 !     replace this by the function from the optimized BLAS implementation:
  93 !         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
  94 !--------------------------------------------------------------
  95
  96       INTEGER i,incX,M,MP1,N
  97       KPP_REAL X(N),Alpha
  98       KPP_REAL ZERO, ONE
  99       PARAMETER( ZERO = 0.0d0 )
 100       PARAMETER( ONE  = 1.0d0 )
 101
 102       IF (Alpha .EQ. ONE) RETURN
 103       IF (N .LE. 0) RETURN
 104
 105       M = MOD(N,5)
 106       IF( M .NE. 0 ) THEN
 107         IF (Alpha .EQ. (-ONE)) THEN
 108           DO i = 1,M
 109             X(i) = -X(i)
 110           END DO
 111         ELSEIF (Alpha .EQ. ZERO) THEN
 112           DO i = 1,M
 113             X(i) = ZERO
 114           END DO
 115         ELSE
 116           DO i = 1,M
 117             X(i) = Alpha*X(i)
 118           END DO
 119         END IF
 120         IF( N .LT. 5 ) RETURN
 121       END IF
 122       MP1 = M + 1
 123       IF (Alpha .EQ. (-ONE)) THEN
 124         DO i = MP1,N,5
 125           X(i)     = -X(i)
 126           X(i + 1) = -X(i + 1)
 127           X(i + 2) = -X(i + 2)
 128           X(i + 3) = -X(i + 3)
 129           X(i + 4) = -X(i + 4)
 130         END DO
 131       ELSEIF (Alpha .EQ. ZERO) THEN
 132         DO i = MP1,N,5
 133           X(i)     = ZERO
 134           X(i + 1) = ZERO
 135           X(i + 2) = ZERO
 136           X(i + 3) = ZERO
 137           X(i + 4) = ZERO
 138         END DO
 139       ELSE
 140         DO i = MP1,N,5
 141           X(i)     = Alpha*X(i)
 142           X(i + 1) = Alpha*X(i + 1)
 143           X(i + 2) = Alpha*X(i + 2)
 144           X(i + 3) = Alpha*X(i + 3)
 145           X(i + 4) = Alpha*X(i + 4)
 146         END DO
 147       END IF
 148       RETURN
 149       END
 150
 151 !--------------------------------------------------------------
 152       KPP_REAL FUNCTION WLAMCH( C )
 153 !--------------------------------------------------------------
 154 !     returns epsilon machine
 155 !     after LAPACK
 156 !     replace this by the function from the optimized LAPACK implementation:
 157 !          CALL SLAMCH('E') or CALL DLAMCH('E')
 158 !--------------------------------------------------------------
 159
 160       CHARACTER C
 161       INTEGER   i
 162       KPP_REAL  ONE, HALF, Eps, Sum
 163       PARAMETER (ONE  = 1.0d0)
 164       PARAMETER (HALF = 0.5d0)
 165       LOGICAL   First
 166       SAVE     First, Eps
 167       DATA     First /.TRUE./
 168
 169       IF (First) THEN
 170         First = .FALSE.
 171         Eps = HALF**(16)
 172         DO i = 17, 80
 173           Eps = Eps*HALF
 174           CALL WLAMCH_ADD(ONE,Eps,Sum)
 175           IF (Sum.LE.ONE) GOTO 10
 176         END DO
 177         PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
 178         RETURN
 179 10      Eps = Eps*2
 180         i = i-1
 181       END IF
 182
 183       WLAMCH = Eps
 184
 185       RETURN
 186       END
 187
 188       SUBROUTINE WLAMCH_ADD( A, B, Sum )
 189       KPP_REAL A, B, Sum
 190       Sum = A + B
 191       RETURN
 192       END
 193 !--------------------------------------------------------------