[2696] | 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 | |
---|
| 60 | INTEGER :: i,incX,incY,M,MP1,N |
---|
| 61 | KPP_REAL :: X(N),Y(N),Alpha |
---|
| 62 | KPP_REAL, PARAMETER :: ZERO = 0.0_dp |
---|
| 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 | |
---|
| 82 | END SUBROUTINE WAXPY |
---|
| 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, PARAMETER :: ZERO=0.0_dp, ONE=1.0_dp |
---|
| 99 | |
---|
| 100 | IF (Alpha .EQ. ONE) RETURN |
---|
| 101 | IF (N .LE. 0) RETURN |
---|
| 102 | |
---|
| 103 | M = MOD(N,5) |
---|
| 104 | IF( M .NE. 0 ) THEN |
---|
| 105 | IF (Alpha .EQ. (-ONE)) THEN |
---|
| 106 | DO i = 1,M |
---|
| 107 | X(i) = -X(i) |
---|
| 108 | END DO |
---|
| 109 | ELSEIF (Alpha .EQ. ZERO) THEN |
---|
| 110 | DO i = 1,M |
---|
| 111 | X(i) = ZERO |
---|
| 112 | END DO |
---|
| 113 | ELSE |
---|
| 114 | DO i = 1,M |
---|
| 115 | X(i) = Alpha*X(i) |
---|
| 116 | END DO |
---|
| 117 | END IF |
---|
| 118 | IF( N .LT. 5 ) RETURN |
---|
| 119 | END IF |
---|
| 120 | MP1 = M + 1 |
---|
| 121 | IF (Alpha .EQ. (-ONE)) THEN |
---|
| 122 | DO i = MP1,N,5 |
---|
| 123 | X(i) = -X(i) |
---|
| 124 | X(i + 1) = -X(i + 1) |
---|
| 125 | X(i + 2) = -X(i + 2) |
---|
| 126 | X(i + 3) = -X(i + 3) |
---|
| 127 | X(i + 4) = -X(i + 4) |
---|
| 128 | END DO |
---|
| 129 | ELSEIF (Alpha .EQ. ZERO) THEN |
---|
| 130 | DO i = MP1,N,5 |
---|
| 131 | X(i) = ZERO |
---|
| 132 | X(i + 1) = ZERO |
---|
| 133 | X(i + 2) = ZERO |
---|
| 134 | X(i + 3) = ZERO |
---|
| 135 | X(i + 4) = ZERO |
---|
| 136 | END DO |
---|
| 137 | ELSE |
---|
| 138 | DO i = MP1,N,5 |
---|
| 139 | X(i) = Alpha*X(i) |
---|
| 140 | X(i + 1) = Alpha*X(i + 1) |
---|
| 141 | X(i + 2) = Alpha*X(i + 2) |
---|
| 142 | X(i + 3) = Alpha*X(i + 3) |
---|
| 143 | X(i + 4) = Alpha*X(i + 4) |
---|
| 144 | END DO |
---|
| 145 | END IF |
---|
| 146 | |
---|
| 147 | END SUBROUTINE WSCAL |
---|
| 148 | |
---|
| 149 | !-------------------------------------------------------------- |
---|
| 150 | KPP_REAL FUNCTION WLAMCH( C ) |
---|
| 151 | !-------------------------------------------------------------- |
---|
| 152 | ! returns epsilon machine |
---|
| 153 | ! after LAPACK |
---|
| 154 | ! replace this by the function from the optimized LAPACK implementation: |
---|
| 155 | ! CALL SLAMCH('E') or CALL DLAMCH('E') |
---|
| 156 | !-------------------------------------------------------------- |
---|
| 157 | ! USE KPP_ROOT_Precision |
---|
| 158 | |
---|
| 159 | CHARACTER :: C |
---|
| 160 | INTEGER :: i |
---|
| 161 | KPP_REAL, SAVE :: Eps |
---|
| 162 | KPP_REAL :: Suma |
---|
| 163 | KPP_REAL, PARAMETER :: ONE=1.0_dp, HALF=0.5_dp |
---|
| 164 | LOGICAL, SAVE :: First=.TRUE. |
---|
| 165 | |
---|
| 166 | IF (First) THEN |
---|
| 167 | First = .FALSE. |
---|
| 168 | Eps = HALF**(16) |
---|
| 169 | DO i = 17, 80 |
---|
| 170 | Eps = Eps*HALF |
---|
| 171 | CALL WLAMCH_ADD(ONE,Eps,Suma) |
---|
| 172 | IF (Suma.LE.ONE) GOTO 10 |
---|
| 173 | END DO |
---|
| 174 | PRINT*,'ERROR IN WLAMCH. EPS < ',Eps |
---|
| 175 | RETURN |
---|
| 176 | 10 Eps = Eps*2 |
---|
| 177 | i = i-1 |
---|
| 178 | END IF |
---|
| 179 | |
---|
| 180 | WLAMCH = Eps |
---|
| 181 | |
---|
| 182 | END FUNCTION WLAMCH |
---|
| 183 | |
---|
| 184 | SUBROUTINE WLAMCH_ADD( A, B, Suma ) |
---|
| 185 | ! USE KPP_ROOT_Precision |
---|
| 186 | |
---|
| 187 | KPP_REAL A, B, Suma |
---|
| 188 | Suma = A + B |
---|
| 189 | |
---|
| 190 | END SUBROUTINE WLAMCH_ADD |
---|
| 191 | !-------------------------------------------------------------- |
---|
| 192 | |
---|
| 193 | |
---|
| 194 | !-------------------------------------------------------------- |
---|
| 195 | SUBROUTINE SET2ZERO(N,Y) |
---|
| 196 | !-------------------------------------------------------------- |
---|
| 197 | ! copies zeros into the vector y: y <- 0 |
---|
| 198 | ! after BLAS |
---|
| 199 | !-------------------------------------------------------------- |
---|
| 200 | |
---|
| 201 | INTEGER :: i,M,MP1,N |
---|
| 202 | KPP_REAL :: Y(N) |
---|
| 203 | KPP_REAL, PARAMETER :: ZERO = 0.0d0 |
---|
| 204 | |
---|
| 205 | IF (N.LE.0) RETURN |
---|
| 206 | |
---|
| 207 | M = MOD(N,8) |
---|
| 208 | IF( M .NE. 0 ) THEN |
---|
| 209 | DO i = 1,M |
---|
| 210 | Y(i) = ZERO |
---|
| 211 | END DO |
---|
| 212 | IF( N .LT. 8 ) RETURN |
---|
| 213 | END IF |
---|
| 214 | MP1 = M+1 |
---|
| 215 | DO i = MP1,N,8 |
---|
| 216 | Y(i) = ZERO |
---|
| 217 | Y(i + 1) = ZERO |
---|
| 218 | Y(i + 2) = ZERO |
---|
| 219 | Y(i + 3) = ZERO |
---|
| 220 | Y(i + 4) = ZERO |
---|
| 221 | Y(i + 5) = ZERO |
---|
| 222 | Y(i + 6) = ZERO |
---|
| 223 | Y(i + 7) = ZERO |
---|
| 224 | END DO |
---|
| 225 | |
---|
| 226 | END SUBROUTINE SET2ZERO |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | !-------------------------------------------------------------- |
---|
| 230 | KPP_REAL FUNCTION WDOT (N, DX, incX, DY, incY) |
---|
| 231 | !-------------------------------------------------------------- |
---|
| 232 | ! dot produce: wdot = x(1:N)*y(1:N) |
---|
| 233 | ! only for incX=incY=1 |
---|
| 234 | ! after BLAS |
---|
| 235 | ! replace this by the function from the optimized BLAS implementation: |
---|
| 236 | ! CALL SDOT(N,X,1,Y,1) or CALL DDOT(N,X,1,Y,1) |
---|
| 237 | !-------------------------------------------------------------- |
---|
| 238 | ! USE messy_mecca_kpp_Precision |
---|
| 239 | !-------------------------------------------------------------- |
---|
| 240 | IMPLICIT NONE |
---|
| 241 | INTEGER :: N, incX, incY |
---|
| 242 | KPP_REAL :: DX(N), DY(N) |
---|
| 243 | |
---|
| 244 | INTEGER :: i, IX, IY, M, MP1, NS |
---|
| 245 | |
---|
| 246 | WDOT = 0.0D0 |
---|
| 247 | IF (N .LE. 0) RETURN |
---|
| 248 | IF (incX .EQ. incY) IF (incX-1) 5,20,60 |
---|
| 249 | ! |
---|
| 250 | ! Code for unequal or nonpositive increments. |
---|
| 251 | ! |
---|
| 252 | 5 IX = 1 |
---|
| 253 | IY = 1 |
---|
| 254 | IF (incX .LT. 0) IX = (-N+1)*incX + 1 |
---|
| 255 | IF (incY .LT. 0) IY = (-N+1)*incY + 1 |
---|
| 256 | DO i = 1,N |
---|
| 257 | WDOT = WDOT + DX(IX)*DY(IY) |
---|
| 258 | IX = IX + incX |
---|
| 259 | IY = IY + incY |
---|
| 260 | END DO |
---|
| 261 | RETURN |
---|
| 262 | ! |
---|
| 263 | ! Code for both increments equal to 1. |
---|
| 264 | ! |
---|
| 265 | ! Clean-up loop so remaining vector length is a multiple of 5. |
---|
| 266 | ! |
---|
| 267 | 20 M = MOD(N,5) |
---|
| 268 | IF (M .EQ. 0) GO TO 40 |
---|
| 269 | DO i = 1,M |
---|
| 270 | WDOT = WDOT + DX(i)*DY(i) |
---|
| 271 | END DO |
---|
| 272 | IF (N .LT. 5) RETURN |
---|
| 273 | 40 MP1 = M + 1 |
---|
| 274 | DO i = MP1,N,5 |
---|
| 275 | WDOT = WDOT + DX(i)*DY(i) + DX(i+1)*DY(i+1) + DX(i+2)*DY(i+2) + & |
---|
| 276 | DX(i+3)*DY(i+3) + DX(i+4)*DY(i+4) |
---|
| 277 | END DO |
---|
| 278 | RETURN |
---|
| 279 | ! |
---|
| 280 | ! Code for equal, positive, non-unit increments. |
---|
| 281 | ! |
---|
| 282 | 60 NS = N*incX |
---|
| 283 | DO i = 1,NS,incX |
---|
| 284 | WDOT = WDOT + DX(i)*DY(i) |
---|
| 285 | END DO |
---|
| 286 | |
---|
| 287 | END FUNCTION WDOT |
---|
| 288 | |
---|
| 289 | |
---|
| 290 | !-------------------------------------------------------------- |
---|
| 291 | SUBROUTINE WADD(N,X,Y,Z) |
---|
| 292 | !-------------------------------------------------------------- |
---|
| 293 | ! adds two vectors: z <- x + y |
---|
| 294 | ! BLAS - like |
---|
| 295 | !-------------------------------------------------------------- |
---|
| 296 | ! USE KPP_ROOT_Precision |
---|
| 297 | |
---|
| 298 | INTEGER :: i, M, MP1, N |
---|
| 299 | KPP_REAL :: X(N),Y(N),Z(N) |
---|
| 300 | |
---|
| 301 | IF (N.LE.0) RETURN |
---|
| 302 | |
---|
| 303 | M = MOD(N,5) |
---|
| 304 | IF( M /= 0 ) THEN |
---|
| 305 | DO i = 1,M |
---|
| 306 | Z(i) = X(i) + Y(i) |
---|
| 307 | END DO |
---|
| 308 | IF( N < 5 ) RETURN |
---|
| 309 | END IF |
---|
| 310 | MP1 = M+1 |
---|
| 311 | DO i = MP1,N,5 |
---|
| 312 | Z(i) = X(i) + Y(i) |
---|
| 313 | Z(i + 1) = X(i + 1) + Y(i + 1) |
---|
| 314 | Z(i + 2) = X(i + 2) + Y(i + 2) |
---|
| 315 | Z(i + 3) = X(i + 3) + Y(i + 3) |
---|
| 316 | Z(i + 4) = X(i + 4) + Y(i + 4) |
---|
| 317 | END DO |
---|
| 318 | |
---|
| 319 | END SUBROUTINE WADD |
---|
| 320 | |
---|
| 321 | |
---|
| 322 | |
---|
| 323 | !-------------------------------------------------------------- |
---|
| 324 | SUBROUTINE WGEFA(N,A,Ipvt,info) |
---|
| 325 | !-------------------------------------------------------------- |
---|
| 326 | ! WGEFA FACTORS THE MATRIX A (N,N) BY |
---|
| 327 | ! GAUSS ELIMINATION WITH PARTIAL PIVOTING |
---|
| 328 | ! LINPACK - LIKE |
---|
| 329 | !-------------------------------------------------------------- |
---|
| 330 | ! |
---|
| 331 | INTEGER :: N,Ipvt(N),info |
---|
| 332 | KPP_REAL :: A(N,N) |
---|
| 333 | KPP_REAL :: t, dmax, da |
---|
| 334 | INTEGER :: j,k,l |
---|
| 335 | KPP_REAL, PARAMETER :: ZERO = 0.0, ONE = 1.0 |
---|
| 336 | |
---|
| 337 | info = 0 |
---|
| 338 | |
---|
| 339 | size: IF (n > 1) THEN |
---|
| 340 | |
---|
| 341 | col: DO k = 1, n-1 |
---|
| 342 | |
---|
| 343 | ! find l = pivot index |
---|
| 344 | ! l = idamax(n-k+1,A(k,k),1) + k - 1 |
---|
| 345 | l = k; dmax = abs(A(k,k)) |
---|
| 346 | DO j = k+1,n |
---|
| 347 | da = ABS(A(j,k)) |
---|
| 348 | IF (da > dmax) THEN |
---|
| 349 | l = j; dmax = da |
---|
| 350 | END IF |
---|
| 351 | END DO |
---|
| 352 | Ipvt(k) = l |
---|
| 353 | |
---|
| 354 | ! zero pivot implies this column already triangularized |
---|
| 355 | IF (ABS(A(l,k)) < TINY(ZERO)) THEN |
---|
| 356 | info = k |
---|
| 357 | return |
---|
| 358 | ELSE |
---|
| 359 | IF (l /= k) THEN |
---|
| 360 | t = A(l,k); A(l,k) = A(k,k); A(k,k) = t |
---|
| 361 | END IF |
---|
| 362 | t = -ONE/A(k,k) |
---|
| 363 | CALL WSCAL(n-k,t,A(k+1,k),1) |
---|
| 364 | DO j = k+1, n |
---|
| 365 | t = A(l,j) |
---|
| 366 | IF (l /= k) THEN |
---|
| 367 | A(l,j) = A(k,j); A(k,j) = t |
---|
| 368 | END IF |
---|
| 369 | CALL WAXPY(n-k,t,A(k+1,k),1,A(k+1,j),1) |
---|
| 370 | END DO |
---|
| 371 | END IF |
---|
| 372 | |
---|
| 373 | END DO col |
---|
| 374 | |
---|
| 375 | END IF size |
---|
| 376 | |
---|
| 377 | Ipvt(N) = N |
---|
| 378 | IF (ABS(A(N,N)) == ZERO) info = N |
---|
| 379 | |
---|
| 380 | END SUBROUTINE WGEFA |
---|
| 381 | |
---|
| 382 | |
---|
| 383 | !-------------------------------------------------------------- |
---|
| 384 | SUBROUTINE WGESL(Trans,N,A,Ipvt,b) |
---|
| 385 | !-------------------------------------------------------------- |
---|
| 386 | ! WGESL solves the system |
---|
| 387 | ! a * x = b or trans(a) * x = b |
---|
| 388 | ! using the factors computed by WGEFA. |
---|
| 389 | ! |
---|
| 390 | ! Trans = 'N' to solve A*x = b , |
---|
| 391 | ! = 'T' to solve transpose(A)*x = b |
---|
| 392 | ! LINPACK - LIKE |
---|
| 393 | !-------------------------------------------------------------- |
---|
| 394 | |
---|
| 395 | INTEGER :: N,Ipvt(N) |
---|
| 396 | CHARACTER :: trans |
---|
| 397 | KPP_REAL :: A(N,N),b(N) |
---|
| 398 | KPP_REAL :: t |
---|
| 399 | INTEGER :: k,kb,l |
---|
| 400 | |
---|
| 401 | |
---|
| 402 | SELECT CASE (Trans) |
---|
| 403 | |
---|
| 404 | CASE ('n','N') ! Solve A * x = b |
---|
| 405 | |
---|
| 406 | ! first solve L*y = b |
---|
| 407 | IF (n >= 2) THEN |
---|
| 408 | DO k = 1, n-1 |
---|
| 409 | l = Ipvt(k) |
---|
| 410 | t = b(l) |
---|
| 411 | IF (l /= k) THEN |
---|
| 412 | b(l) = b(k) |
---|
| 413 | b(k) = t |
---|
| 414 | END IF |
---|
| 415 | CALL WAXPY(n-k,t,a(k+1,k),1,b(k+1),1) |
---|
| 416 | END DO |
---|
| 417 | END IF |
---|
| 418 | ! now solve U*x = y |
---|
| 419 | DO kb = 1, n |
---|
| 420 | k = n + 1 - kb |
---|
| 421 | b(k) = b(k)/a(k,k) |
---|
| 422 | t = -b(k) |
---|
| 423 | CALL WAXPY(k-1,t,a(1,k),1,b(1),1) |
---|
| 424 | END DO |
---|
| 425 | |
---|
| 426 | CASE ('t','T') ! Solve transpose(A) * x = b |
---|
| 427 | |
---|
| 428 | ! first solve trans(U)*y = b |
---|
| 429 | DO k = 1, n |
---|
| 430 | t = WDOT(k-1,a(1,k),1,b(1),1) |
---|
| 431 | b(k) = (b(k) - t)/a(k,k) |
---|
| 432 | END DO |
---|
| 433 | ! now solve trans(L)*x = y |
---|
| 434 | IF (n >= 2) THEN |
---|
| 435 | DO kb = 1, n-1 |
---|
| 436 | k = n - kb |
---|
| 437 | b(k) = b(k) + WDOT(n-k,a(k+1,k),1,b(k+1),1) |
---|
| 438 | l = Ipvt(k) |
---|
| 439 | IF (l /= k) THEN |
---|
| 440 | t = b(l); b(l) = b(k); b(k) = t |
---|
| 441 | END IF |
---|
| 442 | END DO |
---|
| 443 | END IF |
---|
| 444 | |
---|
| 445 | END SELECT |
---|
| 446 | |
---|
| 447 | END SUBROUTINE WGESL |
---|