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

   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 !     USE KPP_ROOT_Precision
  22
  23       INTEGER i,incX,incY,M,MP1,N
  24       KPP_REAL X(N),Y(N)
  25
  26       IF (N.LE.0) RETURN
  27
  28       M = MOD(N,8)
  29       IF( M .NE. 0 ) THEN
  30         DO i = 1,M
  31           Y(i) = X(i)
  32         END DO
  33         IF( N .LT. 8 ) RETURN
  34       END IF
  35       MP1 = M+1
  36       DO i = MP1,N,8
  37         Y(i) = X(i)
  38         Y(i + 1) = X(i + 1)
  39         Y(i + 2) = X(i + 2)
  40         Y(i + 3) = X(i + 3)
  41         Y(i + 4) = X(i + 4)
  42         Y(i + 5) = X(i + 5)
  43         Y(i + 6) = X(i + 6)
  44         Y(i + 7) = X(i + 7)
  45       END DO
  46
  47       END SUBROUTINE WCOPY
  48
  49
  50 !--------------------------------------------------------------
  51       SUBROUTINE WAXPY(N,Alpha,X,incX,Y,incY)
  52 !--------------------------------------------------------------
  53 !     constant times a vector plus a vector: y <- y + Alpha*x
  54 !     only for incX=incY=1
  55 !     after BLAS
  56 !     replace this by the function from the optimized BLAS implementation:
  57 !         CALL SAXPY(N,Alpha,X,1,Y,1) or  CALL DAXPY(N,Alpha,X,1,Y,1)
  58 !--------------------------------------------------------------
  59 !      USE KPP_ROOT_Precision
  60
  61       INTEGER i,incX,incY,M,MP1,N
  62       KPP_REAL X(N),Y(N),Alpha
  63       KPP_REAL ZERO
  64       PARAMETER( ZERO = 0.0_dp )
  65
  66       IF (Alpha .EQ. ZERO) RETURN
  67       IF (N .LE. 0) RETURN
  68
  69       M = MOD(N,4)
  70       IF( M .NE. 0 ) THEN
  71         DO i = 1,M
  72           Y(i) = Y(i) + Alpha*X(i)
  73         END DO
  74         IF( N .LT. 4 ) RETURN
  75       END IF
  76       MP1 = M + 1
  77       DO i = MP1,N,4
  78         Y(i) = Y(i) + Alpha*X(i)
  79         Y(i + 1) = Y(i + 1) + Alpha*X(i + 1)
  80         Y(i + 2) = Y(i + 2) + Alpha*X(i + 2)
  81         Y(i + 3) = Y(i + 3) + Alpha*X(i + 3)
  82       END DO
  83
  84       END SUBROUTINE WAXPY
  85
  86
  87
  88 !--------------------------------------------------------------
  89       SUBROUTINE WSCAL(N,Alpha,X,incX)
  90 !--------------------------------------------------------------
  91 !     constant times a vector: x(1:N) <- Alpha*x(1:N)
  92 !     only for incX=incY=1
  93 !     after BLAS
  94 !     replace this by the function from the optimized BLAS implementation:
  95 !         CALL SSCAL(N,Alpha,X,1) or  CALL DSCAL(N,Alpha,X,1)
  96 !--------------------------------------------------------------
  97 !      USE KPP_ROOT_Precision
  98
  99       INTEGER i,incX,M,MP1,N
 100       KPP_REAL X(N),Alpha
 101       KPP_REAL ZERO, ONE
 102       PARAMETER( ZERO = 0.0_dp )
 103       PARAMETER( ONE  = 1.0_dp )
 104
 105       IF (Alpha .EQ. ONE) RETURN
 106       IF (N .LE. 0) RETURN
 107
 108       M = MOD(N,5)
 109       IF( M .NE. 0 ) THEN
 110         IF (Alpha .EQ. (-ONE)) THEN
 111           DO i = 1,M
 112             X(i) = -X(i)
 113           END DO
 114         ELSEIF (Alpha .EQ. ZERO) THEN
 115           DO i = 1,M
 116             X(i) = ZERO
 117           END DO
 118         ELSE
 119           DO i = 1,M
 120             X(i) = Alpha*X(i)
 121           END DO
 122         END IF
 123         IF( N .LT. 5 ) RETURN
 124       END IF
 125       MP1 = M + 1
 126       IF (Alpha .EQ. (-ONE)) THEN
 127         DO i = MP1,N,5
 128           X(i)     = -X(i)
 129           X(i + 1) = -X(i + 1)
 130           X(i + 2) = -X(i + 2)
 131           X(i + 3) = -X(i + 3)
 132           X(i + 4) = -X(i + 4)
 133         END DO
 134       ELSEIF (Alpha .EQ. ZERO) THEN
 135         DO i = MP1,N,5
 136           X(i)     = ZERO
 137           X(i + 1) = ZERO
 138           X(i + 2) = ZERO
 139           X(i + 3) = ZERO
 140           X(i + 4) = ZERO
 141         END DO
 142       ELSE
 143         DO i = MP1,N,5
 144           X(i)     = Alpha*X(i)
 145           X(i + 1) = Alpha*X(i + 1)
 146           X(i + 2) = Alpha*X(i + 2)
 147           X(i + 3) = Alpha*X(i + 3)
 148           X(i + 4) = Alpha*X(i + 4)
 149         END DO
 150       END IF
 151
 152       END SUBROUTINE WSCAL
 153
 154 !--------------------------------------------------------------
 155       KPP_REAL FUNCTION WLAMCH( C )
 156 !--------------------------------------------------------------
 157 !     returns epsilon machine
 158 !     after LAPACK
 159 !     replace this by the function from the optimized LAPACK implementation:
 160 !          CALL SLAMCH('E') or CALL DLAMCH('E')
 161 !--------------------------------------------------------------
 162 !      USE KPP_ROOT_Precision
 163
 164       CHARACTER C
 165       INTEGER   i
 166       KPP_REAL  ONE, HALF, Eps, Sum
 167       PARAMETER (ONE  = 1.0_dp)
 168       PARAMETER (HALF = 0.5_dp)
 169       LOGICAL   First
 170       SAVE     First, Eps
 171       DATA     First /.TRUE./
 172
 173       IF (First) THEN
 174         First = .FALSE.
 175         Eps = HALF**(16)
 176         DO i = 17, 80
 177           Eps = Eps*HALF
 178           CALL WLAMCH_ADD(ONE,Eps,Sum)
 179           IF (Sum.LE.ONE) GOTO 10
 180         END DO
 181         PRINT*,'ERROR IN WLAMCH. EPS < ',Eps
 182         RETURN
 183 10      Eps = Eps*2
 184         i = i-1
 185       END IF
 186
 187       WLAMCH = Eps
 188
 189       END FUNCTION WLAMCH
 190
 191       SUBROUTINE WLAMCH_ADD( A, B, Sum )
 192 !      USE KPP_ROOT_Precision
 193
 194       KPP_REAL A, B, Sum
 195       Sum = A + B
 196
 197       END SUBROUTINE WLAMCH_ADD
 198 !--------------------------------------------------------------
 199
 200
 201 !--------------------------------------------------------------
 202       SUBROUTINE SET2ZERO(N,Y)
 203 !--------------------------------------------------------------
 204 !     copies zeros into the vector y:  y <- 0
 205 !     after BLAS
 206 !--------------------------------------------------------------
 207
 208       INTEGER ::  i,M,MP1,N
 209       KPP_REAL ::  Y(N)
 210       KPP_REAL, PARAMETER :: ZERO = 0.0d0
 211
 212       IF (N.LE.0) RETURN
 213
 214       M = MOD(N,8)
 215       IF( M .NE. 0 ) THEN
 216         DO i = 1,M
 217           Y(i) = ZERO
 218         END DO
 219         IF( N .LT. 8 ) RETURN
 220       END IF
 221       MP1 = M+1
 222       DO i = MP1,N,8
 223         Y(i)     = ZERO
 224         Y(i + 1) = ZERO
 225         Y(i + 2) = ZERO
 226         Y(i + 3) = ZERO
 227         Y(i + 4) = ZERO
 228         Y(i + 5) = ZERO
 229         Y(i + 6) = ZERO
 230         Y(i + 7) = ZERO
 231       END DO
 232
 233       END SUBROUTINE SET2ZERO
 234
 235
 236 !--------------------------------------------------------------
 237       KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY)
 238 !--------------------------------------------------------------
 239 !     dot produce: wdot = x(1:N)*y(1:N)
 240 !     only for incX=incY=1
 241 !     after BLAS
 242 !     replace this by the function from the optimized BLAS implementation:
 243 !         CALL SDOT(N,X,1,Y,1) or  CALL DDOT(N,X,1,Y,1)
 244 !--------------------------------------------------------------
 245 !      USE messy_mecca_kpp_Precision
 246 !--------------------------------------------------------------
 247       IMPLICIT NONE
 248       INTEGER :: N, incX, incY
 249       KPP_REAL :: DX(N), DY(N)
 250
 251       INTEGER :: i, IX, IY, M, MP1, NS
 252
 253       WDOT = 0.0D0
 254       IF (N .LE. 0) RETURN
 255       IF (incX .EQ. incY) IF (incX-1) 5,20,60
 256 !
 257 !     Code for unequal or nonpositive increments.
 258 !
 259     5 IX = 1
 260       IY = 1
 261       IF (incX .LT. 0) IX = (-N+1)*incX + 1
 262       IF (incY .LT. 0) IY = (-N+1)*incY + 1
 263       DO i = 1,N
 264         WDOT = WDOT + DX(IX)*DY(IY)
 265         IX = IX + incX
 266         IY = IY + incY
 267       END DO
 268       RETURN
 269 !
 270 !     Code for both increments equal to 1.
 271 !
 272 !     Clean-up loop so remaining vector length is a multiple of 5.
 273 !
 274    20 M = MOD(N,5)
 275       IF (M .EQ. 0) GO TO 40
 276       DO i = 1,M
 277          WDOT = WDOT + DX(i)*DY(i)
 278       END DO
 279       IF (N .LT. 5) RETURN
 280    40 MP1 = M + 1
 281       DO i = MP1,N,5
 282           WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) +  &
 283                    DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4)
 284       END DO
 285       RETURN
 286 !
 287 !     Code for equal, positive, non-unit increments.
 288 !
 289    60 NS = N*incX
 290       DO i = 1,NS,incX
 291         WDOT = WDOT + DX(i)*DY(i)
 292       END DO
 293
 294       END FUNCTION WDOT