1C -------------------------------------------------------------------- 2C SUNDIALS Copyright Start 3C Copyright (c) 2002-2020, 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 Diagonal ODE example. Stiff case, with BDF/SPGMR, diagonal 13C preconditioner. Solved with preconditioning on left, then with 14C preconditioning on right. 15C -------------------------------------------------------------------- 16C 17C Include MPI-Fortran header file for MPI_COMM_WORLD, MPI types. 18 19 IMPLICIT NONE 20C 21 INCLUDE "mpif.h" 22C 23C The following declaration specification should match C type long int. 24 INTEGER*8 NLOCAL, NEQ, I, IOUT(25), IPAR(2) 25 PARAMETER (NLOCAL=10) 26C 27 INTEGER*4 LNST, LNFE, LNSETUP, LNNI, LNCF, LNETF, LNPE, LNLI, LNPS 28 INTEGER*4 LNCFL, NOUT, MYPE, NPES, IER, METH, IATOL 29 INTEGER*4 ITASK, IPRE, IGS, JOUT 30 INTEGER*8 NST, NFE, NPSET, NPE, NPS, NNI, NLI 31 INTEGER*8 NCFL, NETF, NCFN 32 DOUBLE PRECISION Y(1024), ROUT(10), RPAR(1) 33 DOUBLE PRECISION ATOL, DTOUT, T, ALPHA, RTOL, TOUT, ERMAX, ERRI 34 DOUBLE PRECISION GERMAX, AVDIM 35C 36 DATA ATOL/1.0D-10/, RTOL/1.0D-5/, DTOUT/0.1D0/, NOUT/10/ 37 DATA LNST/3/, LNFE/4/, LNETF/5/, LNCF/6/, LNNI/7/, LNSETUP/8/, 38 1 LNPE/20/, LNLI/22/, LNPS/21/, LNCFL/23/ 39C 40C Get NPES and MYPE. Requires initialization of MPI. 41 CALL MPI_INIT(IER) 42 IF (IER .NE. 0) THEN 43 WRITE(6,5) IER 44 5 FORMAT(///' MPI_ERROR: MPI_INIT returned IER = ', I5) 45 STOP 46 ENDIF 47 CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPES, IER) 48 IF (IER .NE. 0) THEN 49 WRITE(6,6) IER 50 6 FORMAT(///' MPI_ERROR: MPI_COMM_SIZE returned IER = ', I5) 51 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 52 STOP 53 ENDIF 54 CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYPE, IER) 55 IF (IER .NE. 0) THEN 56 WRITE(6,7) IER 57 7 FORMAT(///' MPI_ERROR: MPI_COMM_RANK returned IER = ', I5) 58 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 59 STOP 60 ENDIF 61C 62C Set input arguments. 63 NEQ = NPES * NLOCAL 64 T = 0.0D0 65 METH = 2 66 IATOL = 1 67 ITASK = 1 68 IPRE = 1 69 IGS = 1 70C Set parameter alpha 71 ALPHA = 10.0D0 72C 73C Load IPAR and RPAR 74 IPAR(1) = NLOCAL 75 IPAR(2) = MYPE 76 RPAR(1) = ALPHA 77C 78C Do remaining initializations for first case: IPRE = 1 (prec. on left). 79C 80 DO 10 I = 1, NLOCAL 81 10 Y(I) = 1.0D0 82C 83 IF (MYPE .EQ. 0) THEN 84 WRITE(6,11) NEQ, ALPHA 85 11 FORMAT('Diagonal test problem:'//' NEQ = ', I3, 86 1 ' parameter alpha = ', F8.3) 87 WRITE(6,12) 88 12 FORMAT(' ydot_i = -alpha*i * y_i (i = 1,...,NEQ)') 89 WRITE(6,13) RTOL, ATOL 90 13 FORMAT(' RTOL, ATOL = ', 2E10.1) 91 WRITE(6,14) 92 14 FORMAT(' Method is BDF/NEWTON/SPGMR'/ 93 1 ' Diagonal preconditioner uses approximate Jacobian') 94 WRITE(6,15) NPES 95 15 FORMAT(' Number of processors = ', I3) 96 WRITE(6,16) 97 16 FORMAT(//'Preconditioning on left'/) 98 ENDIF 99C 100 CALL FNVINITP(MPI_COMM_WORLD, 1, NLOCAL, NEQ, IER) 101C 102 IF (IER .NE. 0) THEN 103 WRITE(6,20) IER 104 20 FORMAT(///' SUNDIALS_ERROR: FNVINITP returned IER = ', I5) 105 CALL MPI_FINALIZE(IER) 106 STOP 107 ENDIF 108C 109C initialize SPGMR linear solver module 110 call FSUNSPGMRINIT(1, IPRE, 0, IER) 111 IF (IER .NE. 0) THEN 112 WRITE(6,25) IER 113 25 FORMAT(///' SUNDIALS_ERROR: FSUNSPGMRINIT IER = ', I5) 114 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 115 STOP 116 ENDIF 117C 118 call FSUNSPGMRSETGSTYPE(1, IGS, IER) 119 IF (IER .NE. 0) THEN 120 WRITE(6,27) IER 121 27 FORMAT(///' SUNDIALS_ERROR: FSUNSPGMRSETGSTYPE IER = ', I5) 122 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 123 STOP 124 ENDIF 125C 126 CALL FCVMALLOC(T, Y, METH, IATOL, RTOL, ATOL, 127 1 IOUT, ROUT, IPAR, RPAR, IER) 128C 129 IF (IER .NE. 0) THEN 130 WRITE(6,30) IER 131 30 FORMAT(///' SUNDIALS_ERROR: FCVMALLOC returned IER = ', I5) 132 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 133 STOP 134 ENDIF 135C 136C attach linear solver module to CVLs interface 137 CALL FCVLSINIT(IER) 138 IF (IER .NE. 0) THEN 139 WRITE(6,32) IER 140 32 FORMAT(///' SUNDIALS_ERROR: FCVLSINIT returned IER = ', I5) 141 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 142 STOP 143 ENDIF 144C 145C attach preconditioner to CVLs interface 146 CALL FCVLSSETPREC(1, IER) 147C 148C Loop through tout values, call solver, print output, test for failure. 149 TOUT = DTOUT 150 DO 70 JOUT = 1, NOUT 151C 152 CALL FCVODE(TOUT, T, Y, ITASK, IER) 153C 154 IF (MYPE .EQ. 0) WRITE(6,40) T, IOUT(LNST), IOUT(LNFE) 155 40 FORMAT(' t = ', E10.2, 5X, 'no. steps = ', I5, 156 & ' no. f-s = ', I5) 157C 158 IF (IER .NE. 0) THEN 159 WRITE(6,60) IER, IOUT(15) 160 60 FORMAT(///' SUNDIALS_ERROR: FCVODE returned IER = ', I5, /, 161 & ' Linear Solver returned IER = ', I5) 162 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 163 STOP 164 ENDIF 165C 166 TOUT = TOUT + DTOUT 167 70 CONTINUE 168C 169C Get max. absolute error in the local vector. 170 ERMAX = 0.0D0 171 DO 75 I = 1, NLOCAL 172 ERRI = Y(I) - EXP(-ALPHA * (MYPE * NLOCAL + I) * T) 173 75 ERMAX = MAX(ERMAX, ABS(ERRI)) 174C Get global max. error from MPI_REDUCE call. 175 CALL MPI_REDUCE(ERMAX, GERMAX, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 176 1 0, MPI_COMM_WORLD, IER) 177 IF (IER .NE. 0) THEN 178 WRITE(6,80) IER 179 80 FORMAT(///' MPI_ERROR: MPI_REDUCE returned IER = ', I5) 180 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 181 STOP 182 ENDIF 183 IF (MYPE .EQ. 0) WRITE(6,85) GERMAX 184 85 FORMAT(/'Max. absolute error is ', E10.2/) 185C 186C Print final statistics. 187 NST = IOUT(LNST) 188 NFE = IOUT(LNFE) 189 NPSET = IOUT(LNSETUP) 190 NPE = IOUT(LNPE) 191 NPS = IOUT(LNPS) 192 NNI = IOUT(LNNI) 193 NLI = IOUT(LNLI) 194 AVDIM = DBLE(NLI) / DBLE(NNI) 195 NCFN = IOUT(LNCF) 196 NCFL = IOUT(LNCFL) 197 NETF = IOUT(LNETF) 198 IF (MYPE .EQ. 0) 199 1 WRITE (6,90) NST, NFE, NPSET, NPE, NPS, NNI, NLI, AVDIM, NCFN, 200 & NCFL, NETF 201 90 FORMAT(/'Final statistics:'// 202 & ' number of steps = ', I5, 5X, 203 & 'number of f evals. =', I5/ 204 & ' number of prec. setups = ', I5/ 205 & ' number of prec. evals. = ', I5, 5X, 206 & 'number of prec. solves = ', I5/ 207 & ' number of nonl. iters. = ', I5, 5X, 208 & 'number of lin. iters. = ', I5/ 209 & ' average Krylov subspace dimension (NLI/NNI) = ', F8.4/ 210 & ' number of conv. failures.. nonlinear = ', I3, 211 & ' linear = ', I3/ 212 & ' number of error test failures = ', I3) 213C 214C Re-initialize to run second case: IPRE = 2 (prec. on right). 215 IPRE = 2 216 T = 0.0D0 217 DO 110 I = 1, NLOCAL 218 110 Y(I) = 1.0D0 219C 220 IF (MYPE .EQ. 0) WRITE(6,111) 221 111 FORMAT(//60('-')///'Preconditioning on right'/) 222C 223 CALL FCVREINIT(T, Y, IATOL, RTOL, ATOL, IER) 224C 225 IF (IER .NE. 0) THEN 226 WRITE(6,130) IER 227 130 FORMAT(///' SUNDIALS_ERROR: FCVREINIT returned IER = ', I5) 228 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 229 STOP 230 ENDIF 231C 232 CALL FSUNSPGMRSETPRECTYPE(1, IPRE, IER) 233 IF (IER .NE. 0) THEN 234 WRITE(6,140) IER 235 140 FORMAT(///' SUNDIALS_ERROR: FSUNSPGMRSETPRECTYPE IER = ',I5) 236 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 237 STOP 238 ENDIF 239C 240C Loop through tout values, call solver, print output, test for failure. 241 TOUT = DTOUT 242 DO 170 JOUT = 1, NOUT 243C 244 CALL FCVODE(TOUT, T, Y, ITASK, IER) 245C 246 IF (MYPE .EQ. 0) WRITE(6,40) T, IOUT(LNST), IOUT(LNFE) 247C 248 IF (IER .NE. 0) THEN 249 WRITE(6,60) IER 250 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 251 STOP 252 ENDIF 253C 254 TOUT = TOUT + DTOUT 255 170 CONTINUE 256C 257C Get max. absolute error in the local vector. 258 ERMAX = 0.0D0 259 DO 175 I = 1, NLOCAL 260 ERRI = Y(I) - EXP(-ALPHA * (MYPE * NLOCAL + I) * T) 261 175 ERMAX = MAX(ERMAX, ABS(ERRI)) 262C Get global max. error from MPI_REDUCE call. 263 CALL MPI_REDUCE(ERMAX, GERMAX, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 264 1 0, MPI_COMM_WORLD, IER) 265 IF (IER .NE. 0) THEN 266 WRITE(6,80) IER 267 CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER) 268 STOP 269 ENDIF 270 IF (MYPE .EQ. 0) WRITE(6,85) GERMAX 271C 272C Print final statistics. 273 NST = IOUT(LNST) 274 NFE = IOUT(LNFE) 275 NPSET = IOUT(LNSETUP) 276 NPE = IOUT(LNPE) 277 NPS = IOUT(LNPS) 278 NNI = IOUT(LNNI) 279 NLI = IOUT(LNLI) 280 AVDIM = DBLE(NLI) / DBLE(NNI) 281 NCFN = IOUT(LNCF) 282 NCFL = IOUT(LNCFL) 283 NETF = IOUT(LNETF) 284 IF (MYPE .EQ. 0) 285 1 WRITE (6,90) NST, NFE, NPSET, NPE, NPS, NNI, NLI, AVDIM, NCFN, 286 & NCFL, NETF 287C 288C Free the memory and finalize MPI. 289 CALL FCVFREE 290 CALL MPI_FINALIZE(IER) 291 IF (IER .NE. 0) THEN 292 WRITE(6,195) IER 293 195 FORMAT(///' MPI_ERROR: MPI_FINALIZE returned IER = ', I5) 294 STOP 295 ENDIF 296C 297 STOP 298 END 299C 300C ------------------------------------------------------------------------ 301C 302 SUBROUTINE FCVFUN(T, Y, YDOT, IPAR, RPAR, IER) 303C Routine for right-hand side function f 304 IMPLICIT NONE 305C 306C The following declaration specification should match C type long int. 307 INTEGER*8 IPAR(*), MYPE, I, NLOCAL 308 INTEGER*4 IER 309 DOUBLE PRECISION T, Y(*), YDOT(*), RPAR(*) 310 DOUBLE PRECISION ALPHA 311C 312 NLOCAL = IPAR(1) 313 MYPE = IPAR(2) 314 ALPHA = RPAR(1) 315C 316 DO I = 1, NLOCAL 317 YDOT(I) = -ALPHA * (MYPE * NLOCAL + I) * Y(I) 318 ENDDO 319C 320 IER = 0 321C 322 RETURN 323 END 324C 325C ------------------------------------------------------------------------ 326C 327 SUBROUTINE FCVPSOL(T, Y, FY, R, Z, GAMMA, DELTA, LR, 328 & IPAR, RPAR, IER) 329C Routine to solve preconditioner linear system 330C This routine uses a diagonal preconditioner P = I - gamma*J, 331C where J is a diagonal approximation to the true Jacobian, given by: 332C J = diag(0, 0, 0, -4*alpha, ..., -N*alpha). 333C The vector r is copied to z, and the inverse of P (restricted to the 334C local vector segment) is applied to the vector z. 335 IMPLICIT NONE 336C 337C The following declaration specification should match C type long int. 338 INTEGER*8 IPAR(*), NLOCAL, I, MYPE, ISTART, IBASE 339 INTEGER*4 IER, LR 340 DOUBLE PRECISION T, Y(*), FY(*), R(*), Z(*) 341 DOUBLE PRECISION GAMMA, DELTA, RPAR(*) 342 DOUBLE PRECISION PSUBI, ALPHA 343C 344 NLOCAL = IPAR(1) 345 MYPE = IPAR(2) 346 ALPHA = RPAR(1) 347C 348 DO I = 1, NLOCAL 349 Z(I) = R(I) 350 ENDDO 351C 352 IBASE = MYPE * NLOCAL 353 ISTART = MAX(1, 4 - IBASE) 354 DO I = ISTART, NLOCAL 355 PSUBI = 1.0D0 + GAMMA * ALPHA * (IBASE + I) 356 Z(I) = Z(I) / PSUBI 357 ENDDO 358C 359 RETURN 360 END 361C 362C ------------------------------------------------------------------------ 363C 364 SUBROUTINE FCVPSET(T, Y, FY, JOK, JCUR, GAMMA, H, 365 & IPAR, RPAR, IER) 366C Empty subroutine. Not needed for the preconditioner, but required 367C by the FCVODE module. 368 INTEGER*4 IER, JOK, JCUR 369 DOUBLE PRECISION T, Y(*), FY(*), GAMMA, H 370C The following declaration specification should match C type long int 371 INTEGER*8 IPAR(*) 372 DOUBLE PRECISION RPAR(*) 373 374 IER = 0 375 RETURN 376 END 377