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