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