1C -------------------------------------------------------------------- 2C SUNDIALS Copyright Start 3C Copyright (c) 2002-2021, Lawrence Livermore National Security 4C and Southern Methodist University. 5C All rights reserved. 6C 7C See the top-level LICENSE and NOTICE files for details. 8C 9C SPDX-License-Identifier: BSD-3-Clause 10C SUNDIALS Copyright End 11C -------------------------------------------------------------------- 12C FCVODE Example Problem: 2D kinetics-transport, precond. Krylov 13C solver. 14C 15C An ODE system is generated from the following 2-species diurnal 16C kinetics advection-diffusion PDE system in 2 space dimensions: 17C 18C dc(i)/dt = Kh*(d/dx)**2 c(i) + V*dc(i)/dx + (d/dy)(Kv(y)*dc(i)/dy) 19C + Ri(c1,c2,t) for i = 1,2, where 20C R1(c1,c2,t) = -q1*c1*c3 - q2*c1*c2 + 2*q3(t)*c3 + q4(t)*c2 , 21C R2(c1,c2,t) = q1*c1*c3 - q2*c1*c2 - q4(t)*c2 , 22C Kv(y) = Kv0*exp(y/5) , 23C Kh, V, Kv0, q1, q2, and c3 are constants, and q3(t) and q4(t) 24C vary diurnally. 25C 26C The problem is posed on the square 27C 0 .le. x .le. 20, 30 .le. y .le. 50 (all in km), 28C with homogeneous Neumann boundary conditions, and for time t 29C in 0 .le. t .le. 86400 sec (1 day). 30C 31C The PDE system is treated by central differences on a uniform 32C 10 x 10 mesh, with simple polynomial initial profiles. 33C The problem is solved with CVODE, with the BDF/GMRES method and 34C the block-diagonal part of the Jacobian as a left 35C preconditioner. 36C 37C Note: this program includes modified (64-bit integer) versions 38C of the dense linear solver routines DGEFA and DGESL from LINPACK, 39C and BLAS routines DCOPY and DSCAL. 40C 41C The second and third dimensions of U here must match the values 42C of MX and MY, for consistency with the output statements below. 43C -------------------------------------------------------------------- 44C 45 IMPLICIT NONE 46C 47 INTEGER*4 MX, MY 48 PARAMETER (MX=10, MY=10) 49 INTEGER*4 LENIPAR, LENRPAR 50 PARAMETER (LENIPAR=6+2*MX*MY, LENRPAR=12+8*MX*MY) 51C 52 INTEGER*4 METH,IATOL,ITASK,IER,LNCFL,LNPS 53 INTEGER*4 LNST,LNFE,LNSETUP,LNNI,LNCF,LQ,LH,LNPE,LNLI,LNETF 54 INTEGER*4 JOUT,JPRETYPE,IGSTYPE,MAXL 55C The following declaration specification should match C type long int. 56 INTEGER*8 NEQ, IOUT(25), IPAR(LENIPAR) 57 INTEGER*4 NST,NFE,NPSET,NPE,NPS,NNI,NETF 58 INTEGER*4 NLI,NCFN,NCFL 59 DOUBLE PRECISION ATOL,AVDIM,T,TOUT,TWOHR,RTOL,FLOOR,DELT 60 DOUBLE PRECISION U(2,MX,MY),ROUT(10),RPAR(LENRPAR) 61C 62 DATA TWOHR/7200.0D0/, RTOL/1.0D-5/, FLOOR/100.0D0/, 63 & JPRETYPE/1/, IGSTYPE/1/, MAXL/0/, DELT/0.0D0/ 64 DATA LNST/3/, LNFE/4/, LNETF/5/, LNCF/6/, LNNI/7/, LNSETUP/8/, 65 & LQ/9/, LNPE/20/, LNLI/22/, LNPS/21/, LNCFL/23/ 66 DATA LH/2/ 67C 68C Load problem constants into IPAR, RPAR, and set initial values 69 CALL INITKX(MX, MY, U, IPAR, RPAR) 70C 71C Set other input arguments. 72 NEQ = 2*MX*MY 73 T = 0.0D0 74 METH = 2 75 IATOL = 1 76 ATOL = RTOL * FLOOR 77 ITASK = 1 78C 79 WRITE(6,10) NEQ 80 10 FORMAT('Krylov example problem:'// 81 & ' Kinetics-transport, NEQ = ', I4/) 82C 83C Initialize vector specification 84 CALL FNVINITS(1, NEQ, IER) 85 IF (IER .NE. 0) THEN 86 WRITE(6,20) IER 87 20 FORMAT(///' SUNDIALS_ERROR: FNVINITS returned IER = ', I5) 88 STOP 89 ENDIF 90C 91C Initialize SPGMR linear solver module 92 call FSUNSPGMRINIT(1, JPRETYPE, MAXL, IER) 93 IF (IER .NE. 0) THEN 94 WRITE(6,25) IER 95 25 FORMAT(///' SUNDIALS_ERROR: FSUNSPGMRINIT IER = ', I5) 96 STOP 97 ENDIF 98 call FSUNSPGMRSETGSTYPE(1, IGSTYPE, IER) 99 IF (IER .NE. 0) THEN 100 WRITE(6,27) IER 101 27 FORMAT(///' SUNDIALS_ERROR: FSUNSPGMRSETGSTYPE IER = ', I5) 102 STOP 103 ENDIF 104C 105C Initialize CVODE 106 CALL FCVMALLOC(T, U, METH, IATOL, RTOL, ATOL, 107 & IOUT, ROUT, IPAR, RPAR, IER) 108 IF (IER .NE. 0) THEN 109 WRITE(6,30) IER 110 30 FORMAT(///' SUNDIALS_ERROR: FCVMALLOC returned IER = ', I5) 111 STOP 112 ENDIF 113C 114C attach linear solver module to CVLs interface 115 CALL FCVLSINIT(IER) 116 IF (IER .NE. 0) THEN 117 WRITE(6,40) IER 118 40 FORMAT(///' SUNDIALS_ERROR: FCVLSINIT returned IER = ',I5) 119 CALL FCVFREE 120 STOP 121 ENDIF 122C 123C attach preconditioner to CVLs interface 124 CALL FCVLSSETPREC(1, IER) 125C 126C Loop over output points, call FCVODE, print sample solution values. 127 TOUT = TWOHR 128 DO JOUT = 1, 12 129C 130 CALL FCVODE(TOUT, T, U, ITASK, IER) 131C 132 WRITE(6,50) T, IOUT(LNST), IOUT(LQ), ROUT(LH) 133 50 FORMAT(/' t = ', E11.3, 3X, 'nst = ', I5, 134 & ' q = ', I2, ' h = ', E14.6) 135 WRITE(6,55) U(1,1,1), U(1,5,5), U(1,10,10), 136 & U(2,1,1), U(2,5,5), U(2,10,10) 137 55 FORMAT(' c1 (bot.left/middle/top rt.) = ', 3E14.6/ 138 & ' c2 (bot.left/middle/top rt.) = ', 3E14.6) 139C 140 IF (IER .NE. 0) THEN 141 WRITE(6,60) IER, IOUT(15) 142 60 FORMAT(///' SUNDIALS_ERROR: FCVODE returned IER = ', I5, /, 143 & ' Linear Solver returned IER = ', I5) 144 CALL FCVFREE 145 STOP 146 ENDIF 147C 148 TOUT = TOUT + TWOHR 149C 150 ENDDO 151 152C Print final statistics. 153 NST = IOUT(LNST) 154 NFE = IOUT(LNFE) 155 NPSET = IOUT(LNSETUP) 156 NPE = IOUT(LNPE) 157 NPS = IOUT(LNPS) 158 NNI = IOUT(LNNI) 159 NLI = IOUT(LNLI) 160 AVDIM = DBLE(NLI) / DBLE(NNI) 161 NCFN = IOUT(LNCF) 162 NCFL = IOUT(LNCFL) 163 NETF = IOUT(LNETF) 164 WRITE(6,80) NST, NFE, NPSET, NPE, NPS, NNI, NLI, AVDIM, NCFN, 165 & NCFL, NETF 166 80 FORMAT(//'Final statistics:'// 167 & ' number of steps = ', I5, 5X, 168 & ' number of f evals. =', I5/ 169 & ' number of prec. setups = ', I5/ 170 & ' number of prec. evals. = ', I5, 5X, 171 & ' number of prec. solves = ', I5/ 172 & ' number of nonl. iters. = ', I5, 5X, 173 & ' number of lin. iters. = ', I5/ 174 & ' average Krylov subspace dimension (NLI/NNI) = ', E14.6/ 175 & ' number of conv. failures.. nonlinear = ', I3, 176 & ' linear = ', I3/ 177 & ' number of error test failures = ', I3) 178C 179 CALL FCVFREE 180C 181 STOP 182 END 183 184C ---------------------------------------------------------------- 185 186 SUBROUTINE INITKX(MX, MY, U0, IPAR, RPAR) 187C Routine to set problem constants and initial values 188C 189 IMPLICIT NONE 190C 191 INTEGER*4 MX, MY 192C The following declaration specification should match C type long int. 193 INTEGER*8 IPAR(*) 194 DOUBLE PRECISION RPAR(*) 195C 196 INTEGER*4 MM, JY, JX, P_IPP, P_BD, P_P 197 DOUBLE PRECISION U0 198 DIMENSION U0(2,MX,MY) 199 DOUBLE PRECISION Q1, Q2, Q3, Q4, A3, A4, OM, C3, DY, HDCO 200 DOUBLE PRECISION VDCO, HACO, X, Y 201 DOUBLE PRECISION CX, CY, DKH, DKV0, DX, HALFDA, PI, VEL 202C 203 DATA DKH/4.0D-6/, VEL/0.001D0/, DKV0/1.0D-8/, HALFDA/4.32D4/, 204 1 PI/3.1415926535898D0/ 205C 206C Problem constants 207 MM = MX * MY 208 Q1 = 1.63D-16 209 Q2 = 4.66D-16 210 Q3 = 0.0D0 211 Q4 = 0.0D0 212 A3 = 22.62D0 213 A4 = 7.601D0 214 OM = PI / HALFDA 215 C3 = 3.7D16 216 DX = 20.0D0 / (MX - 1.0D0) 217 DY = 20.0D0 / (MY - 1.0D0) 218 HDCO = DKH / DX**2 219 HACO = VEL / (2.0D0 * DX) 220 VDCO = (1.0D0 / DY**2) * DKV0 221C 222C Load constants in IPAR and RPAR 223 IPAR(1) = MX 224 IPAR(2) = MY 225 IPAR(3) = MM 226C 227 RPAR(1) = Q1 228 RPAR(2) = Q2 229 RPAR(3) = Q3 230 RPAR(4) = Q4 231 RPAR(5) = A3 232 RPAR(6) = A4 233 RPAR(7) = OM 234 RPAR(8) = C3 235 RPAR(9) = DY 236 RPAR(10) = HDCO 237 RPAR(11) = VDCO 238 RPAR(12) = HACO 239C 240C Pointers into IPAR and RPAR 241 P_IPP = 7 242 P_BD = 13 243 P_P = P_BD + 4*MM 244C 245 IPAR(4) = P_IPP 246 IPAR(5) = P_BD 247 IPAR(6) = P_P 248C 249C Set initial profiles. 250 DO JY = 1, MY 251 Y = 30.0D0 + (JY - 1.0D0) * DY 252 CY = (0.1D0 * (Y - 40.0D0))**2 253 CY = 1.0D0 - CY + 0.5D0 * CY**2 254 DO JX = 1, MX 255 X = (JX - 1.0D0) * DX 256 CX = (0.1D0 * (X - 10.0D0))**2 257 CX = 1.0D0 - CX + 0.5D0 * CX**2 258 U0(1,JX,JY) = 1.0D6 * CX * CY 259 U0(2,JX,JY) = 1.0D12 * CX * CY 260 ENDDO 261 ENDDO 262C 263 RETURN 264 END 265 266C ---------------------------------------------------------------- 267 268 SUBROUTINE FCVFUN(T, U, UDOT, IPAR, RPAR, IER) 269C Routine for right-hand side function f 270C 271 IMPLICIT NONE 272C 273 DOUBLE PRECISION T, U(2,*), UDOT(2,*), RPAR(*) 274C The following declaration specification should match C type long int. 275 INTEGER*8 IPAR(*) 276 INTEGER*4 IER 277C 278 INTEGER*4 ILEFT, IRIGHT 279 INTEGER*4 JX, JY, MX, MY, MM, IBLOK0, IBLOK, IDN, IUP 280 DOUBLE PRECISION Q1, Q2, Q3, Q4, A3, A4, OM, C3, DY, HDCO 281 DOUBLE PRECISION VDCO, HACO 282 DOUBLE PRECISION C1, C2, C1DN, C2DN, C1UP, C2UP, C1LT, C2LT 283 DOUBLE PRECISION C1RT, C2RT, CYDN, CYUP, HORD1, HORD2, HORAD1 284 DOUBLE PRECISION HORAD2, QQ1, QQ2, QQ3, QQ4, RKIN1, RKIN2, S 285 DOUBLE PRECISION VERTD1, VERTD2, YDN, YUP 286C 287C Extract constants from IPAR and RPAR 288 MX = IPAR(1) 289 MY = IPAR(2) 290 MM = IPAR(3) 291C 292 Q1 = RPAR(1) 293 Q2 = RPAR(2) 294 Q3 = RPAR(3) 295 Q4 = RPAR(4) 296 A3 = RPAR(5) 297 A4 = RPAR(6) 298 OM = RPAR(7) 299 C3 = RPAR(8) 300 DY = RPAR(9) 301 HDCO = RPAR(10) 302 VDCO = RPAR(11) 303 HACO = RPAR(12) 304C 305C Set diurnal rate coefficients. 306 S = SIN(OM * T) 307 IF (S .GT. 0.0D0) THEN 308 Q3 = EXP(-A3 / S) 309 Q4 = EXP(-A4 / S) 310 ELSE 311 Q3 = 0.0D0 312 Q4 = 0.0D0 313 ENDIF 314 RPAR(3) = Q3 315 RPAR(4) = Q4 316C 317C Loop over all grid points. 318 DO JY = 1, MY 319 YDN = 30.0D0 + (JY - 1.5D0) * DY 320 YUP = YDN + DY 321 CYDN = VDCO * EXP(0.2D0 * YDN) 322 CYUP = VDCO * EXP(0.2D0 * YUP) 323 IBLOK0 = (JY - 1) * MX 324 IDN = -MX 325 IF (JY .EQ. 1) IDN = MX 326 IUP = MX 327 IF (JY .EQ. MY) IUP = -MX 328 DO JX = 1, MX 329 IBLOK = IBLOK0 + JX 330 C1 = U(1,IBLOK) 331 C2 = U(2,IBLOK) 332C Set kinetic rate terms. 333 QQ1 = Q1 * C1 * C3 334 QQ2 = Q2 * C1 * C2 335 QQ3 = Q3 * C3 336 QQ4 = Q4 * C2 337 RKIN1 = -QQ1 - QQ2 + 2.0D0 * QQ3 + QQ4 338 RKIN2 = QQ1 - QQ2 - QQ4 339C Set vertical diffusion terms. 340 C1DN = U(1,IBLOK + IDN) 341 C2DN = U(2,IBLOK + IDN) 342 C1UP = U(1,IBLOK + IUP) 343 C2UP = U(2,IBLOK + IUP) 344 VERTD1 = CYUP * (C1UP - C1) - CYDN * (C1 - C1DN) 345 VERTD2 = CYUP * (C2UP - C2) - CYDN * (C2 - C2DN) 346C Set horizontal diffusion and advection terms. 347 ILEFT = -1 348 IF (JX .EQ. 1) ILEFT = 1 349 IRIGHT = 1 350 IF (JX .EQ. MX) IRIGHT = -1 351 C1LT = U(1,IBLOK + ILEFT) 352 C2LT = U(2,IBLOK + ILEFT) 353 C1RT = U(1,IBLOK + IRIGHT) 354 C2RT = U(2,IBLOK + IRIGHT) 355 HORD1 = HDCO * (C1RT - 2.0D0 * C1 + C1LT) 356 HORD2 = HDCO * (C2RT - 2.0D0 * C2 + C2LT) 357 HORAD1 = HACO * (C1RT - C1LT) 358 HORAD2 = HACO * (C2RT - C2LT) 359C Load all terms into UDOT. 360 UDOT(1,IBLOK) = VERTD1 + HORD1 + HORAD1 + RKIN1 361 UDOT(2,IBLOK) = VERTD2 + HORD2 + HORAD2 + RKIN2 362 ENDDO 363 ENDDO 364C 365 IER = 0 366C 367 RETURN 368 END 369 370C ---------------------------------------------------------------- 371 372 SUBROUTINE FCVPSET(T, U, FU, JOK, JCUR, GAMMA, H, 373 & IPAR, RPAR, IER) 374C Routine to set and preprocess block-diagonal preconditioner. 375C Note: The dimensions in /BDJ/ below assume at most 100 mesh points. 376C 377 IMPLICIT NONE 378C 379 INTEGER*4 IER, JOK, JCUR 380 DOUBLE PRECISION T, U(2,*), FU(*), GAMMA, H 381C The following declaration specification should match C type long int. 382 INTEGER*8 IPAR(*) 383 DOUBLE PRECISION RPAR(*) 384C 385 INTEGER*4 MX, MY, MM, P_IPP, P_BD, P_P 386 DOUBLE PRECISION Q1, Q2, Q3, Q4, C3, DY, HDCO, VDCO 387C 388 IER = 0 389C 390C Extract constants from IPAR and RPAR 391 MX = IPAR(1) 392 MY = IPAR(2) 393 MM = IPAR(3) 394C 395 Q1 = RPAR(1) 396 Q2 = RPAR(2) 397 Q3 = RPAR(3) 398 Q4 = RPAR(4) 399 C3 = RPAR(8) 400 DY = RPAR(9) 401 HDCO = RPAR(10) 402 VDCO = RPAR(11) 403C 404C Extract pointers into IPAR and RPAR 405 P_IPP = IPAR(4) 406 P_BD = IPAR(5) 407 P_P = IPAR(6) 408C 409C If needed, recompute BD 410C 411 IF (JOK .EQ. 1) THEN 412C JOK = 1. Reuse saved BD 413 JCUR = 0 414 ELSE 415C JOK = 0. Compute diagonal Jacobian blocks. 416C (using q4 value computed on last FCVFUN call). 417 CALL PREC_JAC(MX, MY, MM, U, RPAR(P_BD), 418 & Q1, Q2, Q3, Q4, C3, DY, HDCO, VDCO) 419 JCUR = 1 420 ENDIF 421C 422C Copy BD to P 423 CALL DCOPY(4*MM, RPAR(P_BD), 1, RPAR(P_P), 1) 424C 425C Scale P by -GAMMA 426 CALL DSCAL(4*MM, -GAMMA, RPAR(P_P), 1) 427C 428C Perform LU decomposition 429 CALL PREC_LU(MM, RPAR(P_P), IPAR(P_IPP), IER) 430C 431 RETURN 432 END 433 434C ---------------------------------------------------------------- 435 436 SUBROUTINE FCVPSOL(T, U, FU, R, Z, GAMMA, DELTA, LR, 437 & IPAR, RPAR, IER) 438C Routine to solve preconditioner linear system. 439C 440 IMPLICIT NONE 441C 442 INTEGER*4 IER, LR 443C The following declaration specification should match C type long int. 444 INTEGER*8 IPAR(*) 445 DOUBLE PRECISION T, U(*), FU(*), R(*), Z(2,*) 446 DOUBLE PRECISION GAMMA, DELTA, RPAR(*) 447C 448 INTEGER*4 MM, P_IPP, P_P 449C 450 IER = 0 451C 452C Extract constants from IPAR and RPAR 453 MM = IPAR(3) 454C 455C Extract pointers into IPAR and RPAR 456 P_IPP = IPAR(4) 457 P_P = IPAR(6) 458C 459C Copy RHS into Z 460 CALL DCOPY(2*MM, R, 1, Z, 1) 461C 462C Solve the block-diagonal system Px = r using LU factors stored in P 463C and pivot data in IPP, and return the solution in Z. 464 CALL PREC_SOL(MM, RPAR(P_P), IPAR(P_IPP), Z) 465 466 RETURN 467 END 468 469C ---------------------------------------------------------------- 470 471 SUBROUTINE PREC_JAC(MX, MY, MM, U, BD, 472 & Q1, Q2, Q3, Q4, C3, DY, HDCO, VDCO) 473C Routine to compute diagonal Jacobian blocks 474C 475 IMPLICIT NONE 476C 477 INTEGER*4 MX, MY, MM 478 DOUBLE PRECISION U(2,*), BD(2,2,MM) 479 DOUBLE PRECISION Q1, Q2, Q3, Q4, C3, DY, HDCO, VDCO 480C 481 INTEGER*4 JY, JX, IBLOK, IBLOK0 482 DOUBLE PRECISION C1, C2, CYDN, CYUP, DIAG, YDN, YUP 483C 484 DO JY = 1, MY 485 YDN = 30.0D0 + (JY - 1.5D0) * DY 486 YUP = YDN + DY 487 CYDN = VDCO * EXP(0.2D0 * YDN) 488 CYUP = VDCO * EXP(0.2D0 * YUP) 489 DIAG = -(CYDN + CYUP + 2.0D0 * HDCO) 490 IBLOK0 = (JY - 1) * MX 491 DO JX = 1, MX 492 IBLOK = IBLOK0 + JX 493 C1 = U(1,IBLOK) 494 C2 = U(2,IBLOK) 495 BD(1,1,IBLOK) = (-Q1 * C3 - Q2 * C2) + DIAG 496 BD(1,2,IBLOK) = -Q2 * C1 + Q4 497 BD(2,1,IBLOK) = Q1 * C3 - Q2 * C2 498 BD(2,2,IBLOK) = (-Q2 * C1 - Q4) + DIAG 499 ENDDO 500 ENDDO 501 502 RETURN 503 END 504 505C ---------------------------------------------------------------- 506 507 SUBROUTINE PREC_LU(MM, P, IPP, IER) 508C Routine to perform LU decomposition on (P+I) 509C 510 IMPLICIT NONE 511C 512 INTEGER*4 IER 513 INTEGER*4 MM 514 INTEGER*8 IPP(2,MM) 515 DOUBLE PRECISION P(2,2,MM) 516C 517 INTEGER*4 I 518C 519C Add identity matrix and do LU decompositions on blocks, in place. 520 DO I = 1, MM 521 P(1,1,I) = P(1,1,I) + 1.0D0 522 P(2,2,I) = P(2,2,I) + 1.0D0 523 CALL DGEFA(P(1,1,I), 2, 2, IPP(1,I), IER) 524 IF (IER .NE. 0) RETURN 525 ENDDO 526C 527 RETURN 528 END 529 530C ---------------------------------------------------------------- 531 532 SUBROUTINE PREC_SOL(MM, P, IPP, Z) 533C Routine for backsolve 534C 535 IMPLICIT NONE 536C 537 INTEGER*4 MM 538 INTEGER*8 IPP(2,MM) 539 DOUBLE PRECISION P(2,2,MM), Z(2,MM) 540C 541 INTEGER*4 I 542C 543 DO I = 1, MM 544 CALL DGESL(P(1,1,I), 2, 2, IPP(1,I), Z(1,I), 0) 545 ENDDO 546 547 RETURN 548 END 549 550C ---------------------------------------------------------------- 551 552 subroutine dgefa(a, lda, n, ipvt, info) 553c 554 implicit none 555c 556 integer info, idamax, j, k, kp1, l, nm1, n 557 integer*4 lda 558 integer*8 ipvt(1) 559 double precision a(lda,1), t 560c 561c dgefa factors a double precision matrix by gaussian elimination. 562c 563c dgefa is usually called by dgeco, but it can be called 564c directly with a saving in time if rcond is not needed. 565c (time for dgeco) = (1 + 9/n)*(time for dgefa) . 566c 567c on entry 568c 569c a double precision(lda, n) 570c the matrix to be factored. 571c 572c lda integer 573c the leading dimension of the array a . 574c 575c n integer 576c the order of the matrix a . 577c 578c on return 579c 580c a an upper triangular matrix and the multipliers 581c which were used to obtain it. 582c the factorization can be written a = l*u where 583c l is a product of permutation and unit lower 584c triangular matrices and u is upper triangular. 585c 586c ipvt integer(n) 587c an integer vector of pivot indices. 588c 589c info integer 590c = 0 normal value. 591c = k if u(k,k) .eq. 0.0 . this is not an error 592c condition for this subroutine, but it does 593c indicate that dgesl or dgedi will divide by zero 594c if called. employ rcond in dgeco for a reliable 595c indication of singularity. 596c 597c linpack. this version dated 08/14/78 . 598c cleve moler, university of new mexico, argonne national lab. 599c 600c subroutines and functions 601c 602c blas daxpy,dscal,idamax 603c 604c internal variables 605c 606c gaussian elimination with partial pivoting 607c 608 info = 0 609 nm1 = n - 1 610 if (nm1 .lt. 1) go to 70 611 do 60 k = 1, nm1 612 kp1 = k + 1 613c 614c find l = pivot index 615c 616 l = idamax(n - k + 1, a(k,k), 1) + k - 1 617 ipvt(k) = l 618c 619c zero pivot implies this column already triangularized 620c 621 if (a(l,k) .eq. 0.0d0) go to 40 622c 623c interchange if necessary 624c 625 if (l .eq. k) go to 10 626 t = a(l,k) 627 a(l,k) = a(k,k) 628 a(k,k) = t 629 10 continue 630c 631c compute multipliers 632c 633 t = -1.0d0 / a(k,k) 634 call dscal(n - k, t, a(k + 1,k), 1) 635c 636c row elimination with column indexing 637c 638 do 30 j = kp1, n 639 t = a(l,j) 640 if (l .eq. k) go to 20 641 a(l,j) = a(k,j) 642 a(k,j) = t 643 20 continue 644 call daxpy(n - k, t, a(k + 1,k), 1, a(k + 1,j), 1) 645 30 continue 646 go to 50 647 40 continue 648 info = k 649 50 continue 650 60 continue 651 70 continue 652 ipvt(n) = n 653 if (a(n,n) .eq. 0.0d0) info = n 654 return 655 end 656 657C ---------------------------------------------------------------- 658 659 subroutine dgesl(a, lda, n, ipvt, b, job) 660c 661 implicit none 662c 663 integer lda, n, job, k, kb, l, nm1 664 integer*8 ipvt(1) 665 double precision a(lda,1), b(1), ddot, t 666c 667c dgesl solves the double precision system 668c a * x = b or trans(a) * x = b 669c using the factors computed by dgeco or dgefa. 670c 671c on entry 672c 673c a double precision(lda, n) 674c the output from dgeco or dgefa. 675c 676c lda integer 677c the leading dimension of the array a . 678c 679c n integer 680c the order of the matrix a . 681c 682c ipvt integer(n) 683c the pivot vector from dgeco or dgefa. 684c 685c b double precision(n) 686c the right hand side vector. 687c 688c job integer 689c = 0 to solve a*x = b , 690c = nonzero to solve trans(a)*x = b where 691c trans(a) is the transpose. 692c 693c on return 694c 695c b the solution vector x . 696c 697c error condition 698c 699c a division by zero will occur if the input factor contains a 700c zero on the diagonal. technically this indicates singularity 701c but it is often caused by improper arguments or improper 702c setting of lda . it will not occur if the subroutines are 703c called correctly and if dgeco has set rcond .gt. 0.0 704c or dgefa has set info .eq. 0 . 705c 706c to compute inverse(a) * c where c is a matrix 707c with p columns 708c call dgeco(a,lda,n,ipvt,rcond,z) 709c if (rcond is too small) go to ... 710c do 10 j = 1, p 711c call dgesl(a,lda,n,ipvt,c(1,j),0) 712c 10 continue 713c 714c linpack. this version dated 08/14/78 . 715c cleve moler, university of new mexico, argonne national lab. 716c 717c subroutines and functions 718c 719c blas daxpy,ddot 720c 721c internal variables 722c 723 nm1 = n - 1 724 if (job .ne. 0) go to 50 725c 726c job = 0 , solve a * x = b 727c first solve l*y = b 728c 729 if (nm1 .lt. 1) go to 30 730 do 20 k = 1, nm1 731 l = ipvt(k) 732 t = b(l) 733 if (l .eq. k) go to 10 734 b(l) = b(k) 735 b(k) = t 736 10 continue 737 call daxpy(n - k, t, a(k + 1,k), 1, b(k + 1), 1) 738 20 continue 739 30 continue 740c 741c now solve u*x = y 742c 743 do 40 kb = 1, n 744 k = n + 1 - kb 745 b(k) = b(k) / a(k,k) 746 t = -b(k) 747 call daxpy(k - 1, t, a(1,k), 1, b(1), 1) 748 40 continue 749 go to 100 750 50 continue 751c 752c job = nonzero, solve trans(a) * x = b 753c first solve trans(u)*y = b 754c 755 do 60 k = 1, n 756 t = ddot(k - 1, a(1,k), 1, b(1), 1) 757 b(k) = (b(k) - t) / a(k,k) 758 60 continue 759c 760c now solve trans(l)*x = y 761c 762 if (nm1 .lt. 1) go to 90 763 do 80 kb = 1, nm1 764 k = n - kb 765 b(k) = b(k) + ddot(n - k, a(k + 1,k), 1, b(k + 1), 1) 766 l = ipvt(k) 767 if (l .eq. k) go to 70 768 t = b(l) 769 b(l) = b(k) 770 b(k) = t 771 70 continue 772 80 continue 773 90 continue 774 100 continue 775 return 776 end 777 778C ---------------------------------------------------------------- 779 780 subroutine daxpy(n, da, dx, incx, dy, incy) 781c 782c constant times a vector plus a vector. 783c uses unrolled loops for increments equal to one. 784c jack dongarra, linpack, 3/11/78. 785c 786 implicit none 787c 788 integer i, incx, incy, ix, iy, m, mp1 789 integer*4 n 790 double precision dx(1), dy(1), da 791c 792 if (n .le. 0) return 793 if (da .eq. 0.0d0) return 794 if (incx .eq. 1 .and. incy .eq. 1) go to 20 795c 796c code for unequal increments or equal increments 797c not equal to 1 798c 799 ix = 1 800 iy = 1 801 if (incx .lt. 0) ix = (-n + 1) * incx + 1 802 if (incy .lt. 0) iy = (-n + 1) * incy + 1 803 do 10 i = 1, n 804 dy(iy) = dy(iy) + da * dx(ix) 805 ix = ix + incx 806 iy = iy + incy 807 10 continue 808 return 809c 810c code for both increments equal to 1 811c 812c 813c clean-up loop 814c 815 20 m = mod(n, 4) 816 if ( m .eq. 0 ) go to 40 817 do 30 i = 1, m 818 dy(i) = dy(i) + da * dx(i) 819 30 continue 820 if ( n .lt. 4 ) return 821 40 mp1 = m + 1 822 do 50 i = mp1, n, 4 823 dy(i) = dy(i) + da * dx(i) 824 dy(i + 1) = dy(i + 1) + da * dx(i + 1) 825 dy(i + 2) = dy(i + 2) + da * dx(i + 2) 826 dy(i + 3) = dy(i + 3) + da * dx(i + 3) 827 50 continue 828 return 829 end 830C ---------------------------------------------------------------- 831 832 subroutine dscal(n, da, dx, incx) 833c 834c scales a vector by a constant. 835c uses unrolled loops for increment equal to one. 836c jack dongarra, linpack, 3/11/78. 837c 838 implicit none 839c 840 integer i, incx, m, mp1, nincx 841 integer*4 n 842 double precision da, dx(1) 843c 844 if (n.le.0) return 845 if (incx .eq. 1) go to 20 846c 847c code for increment not equal to 1 848c 849 nincx = n * incx 850 do 10 i = 1, nincx, incx 851 dx(i) = da * dx(i) 852 10 continue 853 return 854c 855c code for increment equal to 1 856c 857c 858c clean-up loop 859c 860 20 m = mod(n, 5) 861 if ( m .eq. 0 ) go to 40 862 do 30 i = 1, m 863 dx(i) = da * dx(i) 864 30 continue 865 if ( n .lt. 5 ) return 866 40 mp1 = m + 1 867 do 50 i = mp1, n, 5 868 dx(i) = da * dx(i) 869 dx(i + 1) = da * dx(i + 1) 870 dx(i + 2) = da * dx(i + 2) 871 dx(i + 3) = da * dx(i + 3) 872 dx(i + 4) = da * dx(i + 4) 873 50 continue 874 return 875 end 876 877C ---------------------------------------------------------------- 878 879 double precision function ddot(n, dx, incx, dy, incy) 880c 881c forms the dot product of two vectors. 882c uses unrolled loops for increments equal to one. 883c jack dongarra, linpack, 3/11/78. 884c 885 implicit none 886c 887 integer i, incx, incy, ix, iy, m, mp1 888 integer*4 n 889 double precision dx(1), dy(1), dtemp 890c 891 ddot = 0.0d0 892 dtemp = 0.0d0 893 if (n .le. 0) return 894 if (incx .eq. 1 .and. incy .eq. 1) go to 20 895c 896c code for unequal increments or equal increments 897c not equal to 1 898c 899 ix = 1 900 iy = 1 901 if (incx .lt. 0) ix = (-n + 1) * incx + 1 902 if (incy .lt. 0) iy = (-n + 1) * incy + 1 903 do 10 i = 1, n 904 dtemp = dtemp + dx(ix) * dy(iy) 905 ix = ix + incx 906 iy = iy + incy 907 10 continue 908 ddot = dtemp 909 return 910c 911c code for both increments equal to 1 912c 913c 914c clean-up loop 915c 916 20 m = mod(n, 5) 917 if ( m .eq. 0 ) go to 40 918 do 30 i = 1,m 919 dtemp = dtemp + dx(i) * dy(i) 920 30 continue 921 if ( n .lt. 5 ) go to 60 922 40 mp1 = m + 1 923 do 50 i = mp1, n, 5 924 dtemp = dtemp + dx(i) * dy(i) + dx(i + 1) * dy(i + 1) + 925 * dx(i + 2) * dy(i + 2) + dx(i + 3) * dy(i + 3) + 926 * dx(i + 4) * dy(i + 4) 927 50 continue 928 60 ddot = dtemp 929 return 930 end 931 932C ---------------------------------------------------------------- 933 934 integer function idamax(n, dx, incx) 935c 936c finds the index of element having max. absolute value. 937c jack dongarra, linpack, 3/11/78. 938c 939 implicit none 940c 941 integer i, incx, ix 942 integer*4 n 943 double precision dx(1), dmax 944c 945 idamax = 0 946 if (n .lt. 1) return 947 idamax = 1 948 if (n .eq. 1) return 949 if (incx .eq. 1) go to 20 950c 951c code for increment not equal to 1 952c 953 ix = 1 954 dmax = abs(dx(1)) 955 ix = ix + incx 956 do 10 i = 2, n 957 if (abs(dx(ix)) .le. dmax) go to 5 958 idamax = i 959 dmax = abs(dx(ix)) 960 5 ix = ix + incx 961 10 continue 962 return 963c 964c code for increment equal to 1 965c 966 20 dmax = abs(dx(1)) 967 do 30 i = 2, n 968 if (abs(dx(i)) .le. dmax) go to 30 969 idamax = i 970 dmax = abs(dx(i)) 971 30 continue 972 return 973 end 974 975C ---------------------------------------------------------------- 976 977 subroutine dcopy(n, dx, incx, dy, incy) 978c 979c copies a vector, x, to a vector, y. 980c uses unrolled loops for increments equal to one. 981c jack dongarra, linpack, 3/11/78. 982c 983 implicit none 984c 985 integer i, incx, incy, ix, iy, m, mp1 986 integer*4 n 987 double precision dx(1), dy(1) 988c 989 if (n .le. 0) return 990 if (incx .eq. 1 .and. incy .eq. 1) go to 20 991c 992c code for unequal increments or equal increments 993c not equal to 1 994c 995 ix = 1 996 iy = 1 997 if (incx .lt. 0) ix = (-n + 1) * incx + 1 998 if (incy .lt. 0) iy = (-n + 1) * incy + 1 999 do 10 i = 1, n 1000 dy(iy) = dx(ix) 1001 ix = ix + incx 1002 iy = iy + incy 1003 10 continue 1004 return 1005c 1006c code for both increments equal to 1 1007c 1008c 1009c clean-up loop 1010c 1011 20 m = mod(n, 7) 1012 if ( m .eq. 0 ) go to 40 1013 do 30 i = 1, m 1014 dy(i) = dx(i) 1015 30 continue 1016 if ( n .lt. 7 ) return 1017 40 mp1 = m + 1 1018 do 50 i = mp1, n, 7 1019 dy(i) = dx(i) 1020 dy(i + 1) = dx(i + 1) 1021 dy(i + 2) = dx(i + 2) 1022 dy(i + 3) = dx(i + 3) 1023 dy(i + 4) = dx(i + 4) 1024 dy(i + 5) = dx(i + 5) 1025 dy(i + 6) = dx(i + 6) 1026 50 continue 1027 return 1028 end 1029