1C     ----------------------------------------------------------------
2C     Programmer(s): Daniel R. Reynolds @ SMU
3C     ----------------------------------------------------------------
4C     SUNDIALS Copyright Start
5C     Copyright (c) 2002-2020, Lawrence Livermore National Security
6C     and Southern Methodist University.
7C     All rights reserved.
8C
9C     See the top-level LICENSE and NOTICE files for details.
10C
11C     SPDX-License-Identifier: BSD-3-Clause
12C     SUNDIALS Copyright End
13C     ----------------------------------------------------------------
14C     Diagonal ODE example. Nonstiff case: alpha = 10/NEQ.
15C     ----------------------------------------------------------------
16C
17C     Include MPI-Fortran header file for MPI_COMM_WORLD, MPI types.
18C
19      IMPLICIT NONE
20C
21      INCLUDE "mpif.h"
22      INTEGER*8 NLOCAL
23      PARAMETER (NLOCAL=2)
24C
25      INTEGER*4 IER, MYPE, NPES, NOUT, LNST, LNST_ATT, LNFE, LNFI, LNNI
26      INTEGER*4 LNCF, LNETF, METH, IATOL, ITASK, JOUT
27      INTEGER*8 NEQ, I, NST, NST_ATT, NFE, NFI, NNI, NCFN, NETF
28      INTEGER*8 IOUT(35), IPAR(2)
29      DOUBLE PRECISION Y(128), ROUT(6), RPAR(1)
30      DOUBLE PRECISION ATOL, RTOL, DTOUT, T, ALPHA, TOUT
31      DOUBLE PRECISION ERMAX, ERRI, GERMAX
32C
33      DATA ATOL/1.0D-10/, RTOL/1.0D-5/, DTOUT/0.1D0/, NOUT/10/
34      DATA LNST/3/, LNST_ATT/6/, LNFE/7/, LNFI/8/, LNNI/11/, LNCF/12/,
35     1     LNETF/10/
36C
37C     Get NPES and MYPE.  Requires initialization of MPI.
38      CALL MPI_INIT(IER)
39      IF (IER .NE. 0) THEN
40         WRITE(6,5) IER
41 5       FORMAT(///' MPI_ERROR: MPI_INIT returned IER = ', I5)
42         STOP
43      ENDIF
44      CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NPES, IER)
45      IF (IER .NE. 0) THEN
46         WRITE(6,6) IER
47 6       FORMAT(///' MPI_ERROR: MPI_COMM_SIZE returned IER = ', I5)
48         CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER)
49         STOP
50      ENDIF
51      CALL MPI_COMM_RANK(MPI_COMM_WORLD, MYPE, IER)
52      IF (IER .NE. 0) THEN
53         WRITE(6,7) IER
54 7       FORMAT(///' MPI_ERROR: MPI_COMM_RANK returned IER = ', I5)
55         CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER)
56         STOP
57      ENDIF
58C
59C     Set input arguments.
60      NEQ = NPES * NLOCAL
61      T = 0.0D0
62      METH = 1
63      IATOL = 1
64      ITASK = 1
65c     Set parameter ALPHA
66      ALPHA  = 10.0D0 / NEQ
67C
68C     Load IPAR and RPAR
69      IPAR(1) = NLOCAL
70      IPAR(2) = MYPE
71      RPAR(1) = ALPHA
72C
73      DO 10 I = 1, NLOCAL
74  10    Y(I) = 1.0D0
75C
76      IF (MYPE .EQ. 0) THEN
77         WRITE(6,11) NEQ, NLOCAL, ALPHA
78  11     FORMAT('Diagonal test problem:'//' NEQ = ', I3, /
79     1        ' NLOCAL = ', I3, /
80     2        ' parameter alpha = ', F8.3)
81         WRITE(6,12)
82  12     FORMAT(' ydot_i = -alpha*i * y_i (i = 1,...,NEQ)')
83         WRITE(6,13) RTOL, ATOL
84  13     FORMAT(' RTOL, ATOL = ', 2E10.1)
85         WRITE(6,14)
86  14     FORMAT(' Method is ERK')
87         WRITE(6,15) NPES
88  15     FORMAT(' Number of processors = ', I3//)
89      ENDIF
90C
91      CALL FNVINITP(MPI_COMM_WORLD, 4, NLOCAL, NEQ, IER)
92C
93      IF (IER .NE. 0) THEN
94         WRITE(6,20) IER
95  20     FORMAT(///' SUNDIALS_ERROR: FNVINITP returned IER = ', I5)
96         CALL MPI_FINALIZE(IER)
97         STOP
98      ENDIF
99C
100      CALL FARKMALLOC(T, Y, METH, IATOL, RTOL, ATOL,
101     1                IOUT, ROUT, IPAR, RPAR, IER)
102C
103      IF (IER .NE. 0) THEN
104         WRITE(6,30) IER
105  30     FORMAT(///' SUNDIALS_ERROR: FARKMALLOC returned IER = ', I5)
106         CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER)
107         STOP
108      ENDIF
109C
110C     Loop through tout values, call solver, print output, test for failure.
111      TOUT = DTOUT
112      DO 70 JOUT = 1, NOUT
113C
114        CALL FARKODE(TOUT, T, Y, ITASK, IER)
115C
116        IF (MYPE .EQ. 0) WRITE(6,40) T, IOUT(LNST), IOUT(LNST_ATT),
117     &       IOUT(LNFE), IOUT(LNFI)
118  40    FORMAT(' t = ', D10.2, 5X, 'steps = ', I5,
119     &         '  (attempted = ', I5, '),  fe = ', I5,
120     &         '  fi = ', I5)
121C
122        IF (IER .NE. 0) THEN
123           WRITE(6,60) IER, IOUT(16)
124  60       FORMAT(///' SUNDIALS_ERROR: FARKODE returned IER = ', I5, /,
125     &            '                 Linear Solver returned IER = ', I5)
126           CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER)
127           STOP
128        ENDIF
129C
130        TOUT = TOUT + DTOUT
131  70    CONTINUE
132C
133C     Get max. absolute error in the local vector.
134      ERMAX = 0.0D0
135      DO 75 I = 1, NLOCAL
136        ERRI  = Y(I) - EXP(-ALPHA * (MYPE * NLOCAL + I) * T)
137        ERMAX = MAX(ERMAX, ABS(ERRI))
138  75    CONTINUE
139C     Get global max. error from MPI_REDUCE call.
140      CALL MPI_REDUCE(ERMAX, GERMAX, 1, MPI_DOUBLE_PRECISION, MPI_MAX,
141     1                0, MPI_COMM_WORLD, IER)
142      IF (IER .NE. 0) THEN
143         WRITE(6,80) IER
144  80     FORMAT(///' MPI_ERROR: MPI_REDUCE returned IER = ', I5)
145         CALL MPI_ABORT(MPI_COMM_WORLD, 1, IER)
146         STOP
147      ENDIF
148      IF (MYPE .EQ. 0) WRITE(6,85) GERMAX
149  85  FORMAT(/'Max. absolute error is ', E10.2/)
150C
151C     Print final statistics.
152      NST = IOUT(LNST)
153      NST_ATT = IOUT(LNST_ATT)
154      NFE = IOUT(LNFE)
155      NFI = IOUT(LNFI)
156      NNI = IOUT(LNNI)
157      NCFN = IOUT(LNCF)
158      NETF = IOUT(LNETF)
159      IF (MYPE .EQ. 0) WRITE (6,90) NST, NST_ATT, NFE, NFI, NNI, NCFN,
160     &     NETF
161  90  FORMAT(/'Final statistics:'//
162     &       ' number of steps = ', I5/
163     &       ' number of step attempts = ', I5/
164     &       ' number of fe evals. = ', I5/
165     &       ' number of fi evals. = ', I5/
166     &       ' number of nonlinear iters. = ', I5/
167     &       ' number of nonlinear conv. failures = ', I3/
168     &       ' number of error test failures = ', I3)
169C
170C     Free the memory and finalize MPI.
171      CALL FARKFREE
172      CALL MPI_FINALIZE(IER)
173      IF (IER .NE. 0) THEN
174         WRITE(6,95) IER
175 95      FORMAT(///' MPI_ERROR: MPI_FINALIZE returned IER = ', I5)
176         STOP
177      ENDIF
178C
179      STOP
180      END
181C
182C     ------------------------------------------------------------------------
183C
184      SUBROUTINE FARKIFUN(T, Y, YDOT, IPAR, RPAR, IER)
185C     Routine for right-hand side function fi
186C
187      IMPLICIT NONE
188C
189      INTEGER*4 IER
190      INTEGER*8 IPAR(*)
191      DOUBLE PRECISION T, Y(*), YDOT(*), RPAR(*)
192C
193      INTEGER*8 NLOCAL, I
194C
195      NLOCAL = IPAR(1)
196C
197      DO I = 1, NLOCAL
198         YDOT(I) = 0.D0
199      ENDDO
200C
201      IER = 0
202C
203      RETURN
204      END
205C     ------------------------------------------------------------------------
206C
207      SUBROUTINE FARKEFUN(T, Y, YDOT, IPAR, RPAR, IER)
208C     Routine for right-hand side function fe
209C
210      IMPLICIT NONE
211C
212      INTEGER*4 IER
213      INTEGER*8 IPAR(*)
214      DOUBLE PRECISION T, Y(*), YDOT(*), RPAR(*)
215C
216      INTEGER*8 MYPE
217      INTEGER*8 NLOCAL, I
218      DOUBLE PRECISION ALPHA
219C
220      NLOCAL = IPAR(1)
221      MYPE = IPAR(2)
222      ALPHA = RPAR(1)
223C
224      DO I = 1, NLOCAL
225         YDOT(I) = -ALPHA * (MYPE * NLOCAL + I) * Y(I)
226      ENDDO
227C
228      IER = 0
229C
230      RETURN
231      END
232