New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trosk.F90 in NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED – NEMO

source: NEMO/branches/2021/dev_r14383_PISCES_NEWDEV_PISCO/src/TOP/PISCES/SED/trosk.F90 @ 15127

Last change on this file since 15127 was 15127, checked in by cetlod, 2 years ago

dev_PISCO : merge with trunk@15119

File size: 51.7 KB
Line 
1MODULE trosk
2# if ! defined key_agrif
3!****************************************************************
4!* NUMERICAL SOLUTION OF A STIFF SYSTEM OF FIRST 0RDER ORDINARY *
5!* DIFFERENTIAL EQUATIONS Y'=F(X,Y) BY ROSENBROCK METHOD.       *
6!* ------------------------------------------------------------ *
7!* SAMPLE RUN:                                                  *
8!* Example #1:                                                  *
9!* (Solve set of differential equations (N=2):                  *
10!*  F(1) = Y(1) * Y(2) + COS(X) - HALF * SIN(TWO * X)           *
11!*  F(2) = Y(1) * Y(1) + Y(2) * Y(2) - (ONE + SIN(X))           *
12!*  Find values of F(1), F(2) at X=1.5).                        *
13!*                                                              *
14!*  SOLUTION AT X=    1.50000000000000                          *
15!*  Y(1) =  0.12359935E+01                                      *
16!*  Y(2) = -0.10494372E+00                                      *
17!*                                                              *
18!*  LAST STEP SIZE =  4.150113101356574E-002                    *
19!*  ERROR CODE =           1                                    *
20!*                                                              *
21!* Example #2:                                                  *
22!* (Solve set of differential equations (N=5):                  *
23!*  F(1) = Y(2)                                                 *
24!*  F(2) = Y(3)                                                 *
25!*  F(3) = Y(4)                                                 *
26!*  F(4) = Y(5)                                                 *
27!*  F(5) = (45.d0 * Y(3) * Y(4) * Y(5) -                        *
28!*          40.d0 * Y(4) * Y(4) * Y(4)) / (NINE * Y(3) * Y(3))  *
29!*  Find values of F(1), F(2), ..., F(5) at X=1.5).             *
30!*                                                              *
31!*  SOLUTION AT X=    1.50000000000000                          *
32!*  Y(1) =  0.43639610E+01                                      *
33!*  Y(2) =  0.40000000E+01                                      *
34!*  Y(3) =  0.28284271E+01                                      *
35!*  Y(4) =  0.14790900E-10                                      *
36!*  Y(5) = -0.37712362E+01                                      *
37!*                                                              *
38!*  LAST STEP SIZE =  3.825256949194526E-003                    *
39!*  ERROR CODE =           1                                    *
40!* ------------------------------------------------------------ *
41!* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77   *
42!*       [BIBLI 18].                                            *
43!*                                                              *
44!*                       F90 Release 1.0 By J-P Moreau, Paris   *
45!*                                (www.jpmoreau.fr)             *
46!****************************************************************
47  USE sed
48  USE sedfunc
49
50  INTEGER, PRIVATE :: JINDIC
51
52  PUBLIC rosk
53
54!define example #1
55  INTERFACE
56       SUBROUTINE JAC(NEQ,JINDEX,X,Y,DFY,LDFY)
57         INTEGER, PARAMETER :: WP = KIND(1.0D0)
58         INTEGER NEQ, JINDEX, LDFY
59         REAL(WP), DIMENSION(LDFY,NEQ) :: DFY
60         INTENT(IN) :: NEQ,JINDEX,X,Y,LDFY
61         INTENT(OUT) :: DFY
62       END SUBROUTINE JAC
63  END INTERFACE
64
65
66  CONTAINS
67
68!**********************************************************************
69SUBROUTINE rosk(ROSM,N,JINDEX,X,Y,XEND,H,        &
70                RTOL,ATOL,ITOL,                  &
71                JAC ,IJAC,MLJAC,MUJAC,           &
72                WORK,LWORK,IWORK,LIWORK,IDID,ISTAT,RSTAT)
73! ---------------------------------------------------------------------
74!     NUMERICAL SOLUTION OF A STIFF (OR DIFFERENTIAL ALGEBRAIC)
75!     SYSTEM OF FIRST 0RDER ORDINARY DIFFERENTIAL EQUATIONS  MY'=F(X,Y).
76!     THIS IS AN EMBEDDED ROSENBROCK METHOD OF ORDER (3)4
77!     (WITH STEP SIZE CONTROL).
78!     C.F. SECTION IV.7
79!
80!     AUTHORS: E. HAIRER AND G. WANNER
81!              UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
82!              CH-1211 GENEVE 24, SWITZERLAND
83!              E-MAIL:  HAIRER@CGEUGE51.BITNET,  WANNER@CGEUGE51.BITNET
84!
85!     THIS CODE IS PART OF THE BOOK:
86!         E. HAIRER AND G. WANNER, SOLVING ORDINARY DIFFERENTIAL
87!         EQUATIONS II. STIFF AND DIFFERENTIAL-ALGEBRAIC PROBLEMS.
88!         SPRINGER SERIES IN COMPUTATIONAL MATHEMATICS,
89!         SPRINGER-VERLAG (1990)
90!
91!     VERSION OF OCTOBER 12, 1990
92!
93!     INPUT PARAMETERS
94!     ----------------
95!     N           DIMENSION OF THE SYSTEM
96!
97!     FCN         NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE
98!                 VALUE OF F(X,Y):
99!                    SUBROUTINE FCN(N,X,Y,F)
100!                    REAL*8 X,Y(N),F(N)
101!                    F(1)=...   ETC.
102!
103!     X           INITIAL X-VALUE
104!
105!     Y(N)        INITIAL VALUES FOR Y
106!
107!     XEND        FINAL X-VALUE (XEND-X MAY BE POSITIVE OR NEGATIVE)
108!
109!     H           INITIAL STEP SIZE GUESS;
110!                 FOR STIFF EQUATIONS WITH INITIAL TRANSIENT,
111!                 H=1.D0/(NORM OF F'), USUALLY 1.D-2 OR 1.D-3, IS GOOD.
112!                 THIS CHOICE IS NOT VERY IMPORTANT, THE CODE QUICKLY
113!                 ADAPTS ITS STEP SIZE. STUDY THE CHOSEN VALUES FOR A FEW
114!                 STEPS IN SUBROUTINE "SOLOUT", WHEN YOU ARE NOT SURE.
115!                 (IF H=0.D0, THE CODE PUTS H=1.D-6).
116!
117!     RTOL,ATOL   RELATIVE AND ABSOLUTE ERROR TOLERANCES. THEY
118!                 CAN BE BOTH SCALARS OR ELSE BOTH VECTORS OF LENGTH N.
119!
120!     ITOL        SWITCH FOR RTOL AND ATOL:
121!                   ITOL=0: BOTH RTOL AND ATOL ARE SCALARS.
122!                     THE CODE KEEPS, ROUGHLY, THE LOCAL ERROR OF
123!                     Y(I) BELOW RTOL*ABS(Y(I))+ATOL
124!                   ITOL=1: BOTH RTOL AND ATOL ARE VECTORS.
125!                     THE CODE KEEPS THE LOCAL ERROR OF Y(I) BELOW
126!                     RTOL(I)*ABS(Y(I))+ATOL(I).
127!
128!     JAC         NAME (EXTERNAL) OF THE SUBROUTINE WHICH COMPUTES
129!                 THE PARTIAL DERIVATIVES OF F(X,Y) WITH RESPECT TO Y
130!                 (THIS ROUTINE IS ONLY CALLED IF IJAC=1; SUPPLY
131!                 A DUMMY SUBROUTINE IN THE CASE IJAC=0).
132!                 FOR IJAC=1, THIS SUBROUTINE MUST HAVE THE FORM:
133!                    SUBROUTINE JAC(N,X,Y,DFY,LDFY)
134!                    REAL*8 X,Y(N),DFY(LDFY,N)
135!                    DFY(1,1)= ...
136!                 LDFY, THE COLUMN-LENGTH OF THE ARRAY, IS
137!                 FURNISHED BY THE CALLING PROGRAM.
138!                 IF (MLJAC.EQ.N) THE JACOBIAN IS SUPPOSED TO
139!                    BE FULL AND THE PARTIAL DERIVATIVES ARE
140!                    STORED IN DFY AS
141!                       DFY(I,J) = PARTIAL F(I) / PARTIAL Y(J)
142!                 ELSE, THE JACOBIAN IS TAKEN AS BANDED AND
143!                    THE PARTIAL DERIVATIVES ARE STORED
144!                    DIAGONAL-WISE AS
145!                       DFY(I-J+MUJAC+1,J) = PARTIAL F(I) / PARTIAL Y(J).
146!
147!     IJAC        SWITCH FOR THE COMPUTATION OF THE JACOBIAN:
148!                    IJAC=0: JACOBIAN IS COMPUTED INTERNALLY BY FINITE
149!                       DIFFERENCES, SUBROUTINE "JAC" IS NEVER CALLED.
150!                    IJAC=1: JACOBIAN IS SUPPLIED BY SUBROUTINE JAC.
151!
152!     MLJAC       SWITCH FOR THE BANDED STRUCTURE OF THE JACOBIAN:
153!                    MLJAC=N: JACOBIAN IS A FULL MATRIX. THE LINEAR
154!                       ALGEBRA IS DONE BY FULL-MATRIX GAUSS-ELIMINATION.
155!                    0<=MLJAC<N: MLJAC IS THE LOWER BANDWITH OF JACOBIAN
156!                       MATRIX (>= NUMBER OF NON-ZERO DIAGONALS BELOW
157!                       THE MAIN DIAGONAL).
158!
159!     MUJAC       UPPER BANDWITH OF JACOBIAN  MATRIX (>= NUMBER OF NON-
160!                 ZERO DIAGONALS ABOVE THE MAIN DIAGONAL).
161!                 NEED NOT BE DEFINED IF MLJAC=N.
162!
163!     WORK        ARRAY OF WORKING SPACE OF LENGTH "LWORK".
164!                 SERVES AS WORKING SPACE FOR ALL VECTORS AND MATRICES.
165!                 "LWORK" MUST BE AT LEAST
166!                             N*(LJAC+LE1+8)+5
167!                 WHERE
168!                    LJAC=N              IF MLJAC=N (FULL JACOBIAN)
169!                    LJAC=MLJAC+MUJAC+1  IF MLJAC<N (BANDED JAC.)
170!                 AND
171!                    LE1=N               IF MLJAC=N (FULL JACOBIAN)
172!                    LE1=2*MLJAC+MUJAC+1 IF MLJAC<N (BANDED JAC.).
173!
174!                 IN THE USUAL CASE WHERE THE JACOBIAN IS FULL AND THE
175!                 MASS-MATRIX IS THE INDENTITY (IMAS=0), THE MINIMUM
176!                 STORAGE REQUIREMENT IS
177!                             LWORK = 2*N*N+8*N+5.
178!
179!     LWORK       DECLARED LENGHT OF ARRAY "WORK".
180!
181!     IWORK       INTEGER WORKING SPACE OF LENGTH "LIWORK".
182!                 "LIWORK" MUST BE AT LEAST N+2.
183!
184!     LIWORK      DECLARED LENGHT OF ARRAY "IWORK".
185!
186! ----------------------------------------------------------------------
187!
188!     SOPHISTICATED SETTING OF PARAMETERS
189!     -----------------------------------
190!              SEVERAL PARAMETERS OF THE CODE ARE TUNED TO MAKE IT WORK
191!              WELL. THEY MAY BE DEFINED BY SETTING WORK(1),..,WORK(5)
192!              AS WELL AS IWORK(1),IWORK(2) DIFFERENT FROM ZERO.
193!              FOR ZERO INPUT, THE CODE CHOOSES DEFAULT VALUES:
194!
195!    IWORK(1)  THIS IS THE MAXIMAL NUMBER OF ALLOWED STEPS.
196!              THE DEFAULT VALUE (FOR IWORK(1)=0) IS 100000.
197!
198!    IWORK(2)  SWITCH FOR THE CHOICE OF THE COEFFICIENTS
199!              IF IWORK(2).EQ.1  METHOD OF SHAMPINE
200!              IF IWORK(2).EQ.2  METHOD GRK4T OF KAPS-RENTROP
201!              IF IWORK(2).EQ.3  METHOD GRK4A OF KAPS-RENTROP
202!              IF IWORK(2).EQ.4  METHOD OF VAN VELDHUIZEN (GAMMA=1/2)
203!              IF IWORK(2).EQ.5  METHOD OF VAN VELDHUIZEN ("D-STABLE")
204!              IF IWORK(2).EQ.6  AN L-STABLE METHOD
205!              THE DEFAULT VALUE (FOR IWORK(2)=0) IS IWORK(2)=2.
206!
207!    WORK(1)   UROUND, THE ROUNDING UNIT, DEFAULT 1.D-16.
208!
209!    WORK(2)   MAXIMAL STEP SIZE, DEFAULT XEND-X.
210!
211!    WORK(3), WORK(4)   PARAMETERS FOR STEP SIZE SELECTION
212!              THE NEW STEP SIZE IS CHOSEN SUBJECT TO THE RESTRICTION
213!                 WORK(3) <= HNEW/HOLD <= WORK(4)
214!              DEFAULT VALUES: WORK(3)=0.2D0, WORK(4)=6.D0
215!
216!    WORK(5)   AVOID THE HUMP: AFTER TWO CONSECUTIVE STEP REJECTIONS
217!              THE STEP SIZE IS MULTIPLIED BY WORK(5)
218!              DEFAULT VALUES: WORK(5)=0.1D0
219!
220!-----------------------------------------------------------------------
221!
222!     OUTPUT PARAMETERS
223!     -----------------
224!     X           X-VALUE WHERE THE SOLUTION IS COMPUTED
225!                 (AFTER SUCCESSFUL RETURN X=XEND)
226!
227!     Y(N)        SOLUTION AT X
228!
229!     H           PREDICTED STEP SIZE OF THE LAST ACCEPTED STEP
230!
231!     IDID        REPORTS ON SUCCESSFULNESS UPON RETURN:
232!                   IDID=1  COMPUTATION SUCCESSFUL,
233!                   IDID=-1 COMPUTATION UNSUCCESSFUL.
234!
235! ---------------------------------------------------------
236! *** *** *** *** *** *** *** *** *** *** *** *** ***
237!          DECLARATIONS
238! *** *** *** *** *** *** *** *** *** *** *** *** ***
239      IMPLICIT REAL(wp) (A-H,O-Z)
240      INTEGER :: JINDEX, ROSM
241      DIMENSION Y(N),ATOL(1),RTOL(1),WORK(LWORK),IWORK(LIWORK)
242      INTEGER , DIMENSION(3) :: ISTAT
243      REAL(wp), DIMENSION(2) :: RSTAT
244      LOGICAL JBAND,ARRET
245      EXTERNAL JAC
246      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
247! --------------------------------------------------------------------
248! --- COMMON STAT CAN BE USED FOR STATISTICS
249! ---    NFCN      NUMBER OF FUNCTION EVALUATIONS (THOSE FOR NUMERICAL
250!                  EVALUATION OF THE JACOBIAN ARE NOT COUNTED)
251! ---    NJAC      NUMBER OF JACOBIAN EVALUATIONS (EITHER ANALYTICALLY
252!                  OR NUMERICALLY)
253! ---    NSTEP     NUMBER OF COMPUTED STEPS
254! ---    NACCPT    NUMBER OF ACCEPTED STEPS
255! ---    NREJCT    NUMBER OF REJECTED STEPS (AFTER AT LEAST ONE STEP
256!                  HAS BEEN ACCEPTED)
257! ---    NDEC      NUMBER OF LU-DECOMPOSITIONS (N-DIMENSIONAL MATRIX)
258! ---    NSOL      NUMBER OF FORWARD-BACKWARD SUBSTITUTIONS
259! --------------------------------------------------------------------
260! *** *** *** *** *** *** ***
261!    SETTING THE PARAMETERS
262! *** *** *** *** *** *** ***
263      NFCN=0
264      NJAC=0
265      NSTEP=0
266      NACCPT=0
267      NREJCT=0
268      NDEC=0
269      NSOL=0
270      ARRET=.FALSE.
271      JINDIC = JINDEX
272! -------- NMAX , THE MAXIMAL NUMBER OF STEPS -----
273      IF(IWORK(1).EQ.0)THEN
274         NMAX=100000
275      ELSE
276         NMAX=IWORK(1)
277         IF(NMAX.LE.0)THEN
278            WRITE(NUMSED,*)' WRONG INPUT IWORK(1)=',IWORK(1)
279            ARRET=.TRUE.
280         END IF
281      END IF
282! -------- METH   COEFFICIENTS OF THE METHOD
283      IF(IWORK(2).EQ.0)THEN
284         METH=2
285      ELSE
286         METH=IWORK(2)
287         IF(METH.LE.0.OR.METH.GE.7)THEN
288            WRITE(NUMSED,*)' CURIOUS INPUT IWORK(2)=',IWORK(2)
289            ARRET=.TRUE.
290         END IF
291      END IF
292! -------- UROUND   SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.D0
293      IF(WORK(1).EQ.0.D0)THEN
294         UROUND=1.D-16
295      ELSE
296         UROUND=WORK(1)
297         IF(UROUND.LE.1.D-14.OR.UROUND.GE.1.D0)THEN
298            WRITE(NUMSED,*)' COEFFICIENTS HAVE 16 DIGITS, UROUND=',WORK(1)
299            ARRET=.TRUE.
300         END IF
301      END IF
302! -------- MAXIMAL STEP SIZE
303      IF(WORK(2).EQ.0.D0)THEN
304         HMAX=XEND-X
305      ELSE
306         HMAX=WORK(2)
307      END IF
308! -------  FAC1,FAC2     PARAMETERS FOR STEP SIZE SELECTION
309      IF(WORK(3).EQ.0.D0)THEN
310         FAC1=5.D0
311      ELSE
312         FAC1=1.D0/WORK(3)
313      END IF
314      IF(WORK(4).EQ.0.D0)THEN
315         FAC2=1.D0/6.0D0
316      ELSE
317         FAC2=1.D0/WORK(4)
318      END IF
319! -------  FACREJ    FOR THE HUMP
320      IF(WORK(5).EQ.0.D0)THEN
321         FACREJ=0.1D0
322      ELSE
323         FACREJ=WORK(5)
324      END IF
325! --------- CHECK IF TOLERANCES ARE O.K.
326      IF (ITOL.EQ.0) THEN
327          IF (ATOL(1).LE.0.D0.OR.RTOL(1).LE.10.D0*UROUND) THEN
328              WRITE (NUMSED,*) ' TOLERANCES ARE TOO SMALL'
329              ARRET=.TRUE.
330          END IF
331      ELSE
332          DO 15 I=1,N
333          IF (ATOL(I).LE.0.D0.OR.RTOL(I).LE.10.D0*UROUND) THEN
334              WRITE (NUMSED,*) ' TOLERANCES(',I,') ARE TOO SMALL'
335              ARRET=.TRUE.
336          END IF
337  15      CONTINUE
338      END IF
339! *** *** *** *** *** *** *** *** *** *** *** *** ***
340!         COMPUTATION OF ARRAY ENTRIES
341! *** *** *** *** *** *** *** *** *** *** *** *** ***
342! ---- AUTONOMOUS, IMPLICIT, BANDED OR NOT ?
343      JBAND=MLJAC.NE.N
344      ARRET=.FALSE.
345! -------- COMPUTATION OF THE ROW-DIMENSIONS OF THE 2-ARRAYS ---
346! -- JACOBIAN
347      IF(JBAND)THEN
348         LDJAC=MLJAC+MUJAC+1
349      ELSE
350         LDJAC=N
351      END IF
352! -- MATRIX E FOR LINEAR ALGEBRA
353      IF(JBAND)THEN
354         LDE=2*MLJAC+MUJAC+1
355      ELSE
356         LDE=N
357      END IF
358! ------- PREPARE THE ENTRY-POINTS FOR THE ARRAYS IN WORK -----
359      IEYNEW=6
360      IEDY1=IEYNEW+N
361      IEDY=IEDY1+N
362      IEAK1=IEDY+N
363      IEAK2=IEAK1+N
364      IEAK3=IEAK2+N
365      IEAK4=IEAK3+N
366      IEFX =IEAK4+N
367      IEJAC=IEFX +N
368      IEMAS=IEJAC+N*LDJAC
369      IEE  =IEMAS
370! ------ TOTAL STORAGE REQUIREMENT -----------
371      ISTORE=IEE+N*LDE-1
372      IF(ISTORE.GT.LWORK)THEN
373         WRITE(NUMSED,*)' INSUFFICIENT STORAGE FOR WORK, MIN. LWORK=',ISTORE
374         ARRET=.TRUE.
375      END IF
376! ------- ENTRY POINTS FOR INTEGER WORKSPACE -----
377      IEIP=3
378! --------- TOTAL REQUIREMENT ---------------
379      ISTORE=IEIP+N-1
380      IF(ISTORE.GT.LIWORK)THEN
381         WRITE(NUMSED,*)' INSUFF. STORAGE FOR IWORK, MIN. LIWORK=',ISTORE
382         ARRET=.TRUE.
383      END IF
384! ------ WHEN A FAIL HAS OCCURED, WE RETURN WITH IDID=-1
385      IF (ARRET) THEN
386         IDID=-1
387         RETURN
388      END IF
389! -------- CALL TO CORE INTEGRATOR ------------
390      IF (ROSM == 1) THEN
391         CALL RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC,        &
392            MLJAC,MUJAC,IDID,                 &
393            NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND,     &
394            LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),      &
395            WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),           &
396            WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT)
397      ELSE IF (ROSM == 2) THEN
398! -------- CALL TO CORE INTEGRATOR ------------
399         CALL RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC,        &
400            MLJAC,MUJAC,IDID,                 &
401            NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND,     &
402            LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),      &
403            WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),           &
404            WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT)
405      ELSE IF (ROSM == 3) THEN
406         CALL RO2COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,IJAC,        &
407            MLJAC,MUJAC,IDID,                 &
408            NMAX,UROUND,METH,FAC1,FAC2,FACREJ,JBAND,     &
409            LDJAC,LDE,WORK(IEYNEW),WORK(IEDY1),WORK(IEDY),      &
410            WORK(IEAK1),WORK(IEAK2),WORK(IEAK3),WORK(IEAK4),           &
411            WORK(IEFX),WORK(IEJAC),WORK(IEE),IWORK(IEIP),RSTAT)
412      ENDIF
413! ----------- RETURN -----------
414
415      ISTAT(1) = NFCN
416      ISTAT(2) = NJAC
417      ISTAT(3) = NSTEP
418
419      RETURN
420
421      END SUBROUTINE rosk
422
423      SUBROUTINE RO4COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,         &
424        IJAC,MLJAC,MUJAC,IDID,               &
425        NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED,       &
426        LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT)
427! ----------------------------------------------------------
428!     CORE INTEGRATOR FOR ROS4
429!     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
430! ----------------------------------------------------------
431! ----------------------------------------------------------
432!         DECLARATIONS
433! ----------------------------------------------------------
434      IMPLICIT REAL*8 (A-H,O-Z)
435      REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), &
436            AK2(N),AK3(N),AK4(N),FX(N),            &
437        FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1)
438      INTEGER IP(N)
439      LOGICAL REJECT,RJECT2,BANDED
440      REAL(wp), DIMENSION(2) :: RSTAT
441      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
442
443! ---- PREPARE BANDWIDTHS -----
444      IF (BANDED) THEN
445          MLE=MLJAC
446          MUE=MUJAC
447          MBJAC=MLJAC+MUJAC+1
448          MDIAG=MLE+MUE+1
449      END IF
450! *** *** *** *** *** *** ***
451!  INITIALISATIONS
452! *** *** *** *** *** *** ***
453      POSNEG=DSIGN(1.D0,XEND-X)
454      IF (METH.EQ.1) CALL SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
455                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
456      IF (METH.EQ.2) CALL GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
457                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
458      IF (METH.EQ.3) CALL GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
459                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
460      IF (METH.EQ.4) CALL VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
461                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
462      IF (METH.EQ.5) CALL VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
463                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
464      IF (METH.EQ.6) CALL LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,  &
465                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
466
467! --- INITIAL PREPARATIONS
468      HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
469      H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN)
470      H=DSIGN(H,POSNEG)
471      REJECT=.FALSE.
472      NSING=0
473      IRTRN=1
474      XXOLD=X
475
476      IF (IRTRN.LT.0) GOTO 79
477! --- BASIC INTEGRATION STEP
478   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
479      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
480          H=HOPT
481          IDID=1
482          RETURN
483      END IF
484      HOPT=H
485      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
486
487      CALL sed_func(N,JINDIC,X,Y,DY1)
488
489      NFCN=NFCN+1
490! *** *** *** *** *** *** ***
491!  COMPUTATION OF THE JACOBIAN
492! *** *** *** *** *** *** ***
493      NJAC=NJAC+1
494      IF (IJAC.EQ.0) THEN
495! --- COMPUTE JACOBIAN MATRIX NUMERICALLY
496          IF (BANDED) THEN
497! --- JACOBIAN IS BANDED
498              MUJACP=MUJAC+1
499              MD=MIN(MBJAC,N)
500              DO 16 K=1,MD
501              J=K
502 12           AK2(J)=Y(J)
503              AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
504              Y(J)=Y(J)+AK3(J)
505              J=J+MD
506              IF (J.LE.N) GOTO 12
507              CALL sed_func(N,JINDIC,X,Y,AK1)
508              J=K
509              LBEG=MAX(1,J-MUJAC)
510 14           LEND=MIN(N,J+MLJAC)
511              Y(J)=AK2(J)
512              MUJACJ=MUJACP-J
513              DO 15 L=LBEG,LEND
514 15           FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J)
515              J=J+MD
516              LBEG=LEND+1
517              IF (J.LE.N) GOTO 14
518 16           CONTINUE
519          ELSE
520! --- JACOBIAN IS FULL
521              DO 18 I=1,N
522              YSAFE=Y(I)
523              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
524              Y(I)=YSAFE+DELT
525              CALL sed_func(N,JINDIC,X,Y,AK1)
526              DO 17 J=1,N
527  17          FJAC(J,I)=(AK1(J)-DY1(J))/DELT
528  18          Y(I)=YSAFE
529              MLJAC=N
530          END IF
531      ELSE
532! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
533          CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC)
534      END IF
535   2  CONTINUE
536! *** *** *** *** *** *** ***
537!  COMPUTE THE STAGES
538! *** *** *** *** *** *** ***
539      NDEC=NDEC+1
540      HC21=C21/H
541      HC31=C31/H
542      HC32=C32/H
543      HC41=C41/H
544      HC42=C42/H
545      HC43=C43/H
546      FAC=1.D0/(H*GAMMA)
547          IF (BANDED) THEN
548! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
549              DO 601 J=1,N
550              I1=MAX0(1,MUJAC+2-J)
551              I2=MIN0(MBJAC,N+MUJAC+1-J)
552              DO 600 I=I1,I2
553  600         E(I+MLE,J)=-FJAC(I,J)
554  601         E(MDIAG,J)=E(MDIAG,J)+FAC
555              CALL DECB(N,LE,E,MLE,MUE,IP,INFO)
556              IF (INFO.NE.0) GOTO 80
557! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
558! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
559! ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
560! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
561              DO 603 I=1,N
562  603         AK1(I)=DY1(I)
563              CALL SOLB(N,LE,E,MLE,MUE,AK1,IP)
564              DO 610 I=1,N
565  610         YNEW(I)=Y(I)+A21*AK1(I)
566              CALL sed_func(N,JINDIC,X,YNEW,DY)
567              DO 611 I=1,N
568  611         AK2(I)=DY(I)+HC21*AK1(I)
569              CALL SOLB(N,LE,E,MLE,MUE,AK2,IP)
570              DO 620 I=1,N
571  620         YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
572              CALL sed_func(N,JINDIC,X,YNEW,DY)
573              DO 621 I=1,N
574  621         AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
575              CALL SOLB(N,LE,E,MLE,MUE,AK3,IP)
576              DO 631 I=1,N
577  631         AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I)
578              CALL SOLB(N,LE,E,MLE,MUE,AK4,IP)
579          ELSE
580! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
581              DO 801 J=1,N
582              DO 800 I=1,N
583  800         E(I,J)=-FJAC(I,J)
584  801         E(J,J)=E(J,J)+FAC
585              CALL DEC(N,LE,E,IP,INFO)
586              IF (INFO.NE.0) GOTO 80
587! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
588! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
589! ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
590! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
591              DO 803 I=1,N
592  803         AK1(I)=DY1(I)
593              CALL SOL(N,LE,E,AK1,IP)
594              DO 810 I=1,N
595  810         YNEW(I)=Y(I)+A21*AK1(I)
596              CALL sed_func(N,JINDIC,X,YNEW,DY)
597              DO 811 I=1,N
598  811         AK2(I)=DY(I)+HC21*AK1(I)
599              CALL SOL(N,LE,E,AK2,IP)
600              DO 820 I=1,N
601  820         YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
602              CALL sed_func(N,JINDIC,X,YNEW,DY)
603              DO 821 I=1,N
604  821         AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
605              CALL SOL(N,LE,E,AK3,IP)
606              DO 831 I=1,N
607  831         AK4(I)=DY(I)+HC41*AK1(I)+HC42*AK2(I)+HC43*AK3(I)
608              CALL SOL(N,LE,E,AK4,IP)
609          END IF
610      NSOL=NSOL+4
611      NFCN=NFCN+2
612! *** *** *** *** *** *** ***
613!  ERROR ESTIMATION
614! *** *** *** *** *** *** ***
615      NSTEP=NSTEP+1
616! ------------ NEW SOLUTION ---------------
617      DO 240 I=1,N
618  240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)+B4*AK4(I)
619! ------------ COMPUTE ERROR ESTIMATION ----------------
620      ERR=0.D0
621      DO 300 I=1,N
622      S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)+E4*AK4(I)
623      IF (ITOL.EQ.0) THEN
624         SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
625      ELSE
626         SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
627      END IF
628  300 ERR=ERR+(S/SK)**2
629      ERR=DSQRT(ERR/N)
630! --- COMPUTATION OF HNEW
631! --- WE REQUIRE .2<=HNEW/H<=6.
632      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0))
633      HNEW=H/FAC
634! *** *** *** *** *** *** ***
635!  IS THE ERROR SMALL ENOUGH ?
636! *** *** *** *** *** *** ***
637      RJECT2 = .TRUE.
638      IF (ERR.LE.1.D0) THEN
639! --- STEP IS ACCEPTED
640         NACCPT=NACCPT+1
641         DO 44 I=1,N
642  44     Y(I)=YNEW(I)
643         XXOLD=X
644         X=X+H
645         RSTAT(2) = H
646         IF (IRTRN.LT.0) GOTO 79
647         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
648         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
649         REJECT=.FALSE.
650         RJECT2=.FALSE.
651         IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0
652         H=HNEW
653         GOTO 1
654      ELSE
655! --- STEP IS REJECTED
656         IF (RJECT2) HNEW=H*FACREJ
657         IF (REJECT) RJECT2=.TRUE.
658         REJECT=.TRUE.
659         H=HNEW
660         IF (NACCPT.GE.1) NREJCT=NREJCT+1
661         GOTO 2
662      END IF
663! --- EXIT
664  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
665      NSING=NSING+1
666      IF (NSING.GE.5) GOTO 79
667      H=H*0.5D0
668      GOTO 2
669  79  WRITE(NUMSED,979)X,H
670 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7)
671      IDID=-1
672      RETURN
673
674      END SUBROUTINE RO4COR
675
676
677      SUBROUTINE RO3COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,         &
678        IJAC,MLJAC,MUJAC,IDID,               &
679        NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED,       &
680        LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT)
681! ----------------------------------------------------------
682!     CORE INTEGRATOR FOR ROS4
683!     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
684! ----------------------------------------------------------
685!         DECLARATIONS
686! ----------------------------------------------------------
687      IMPLICIT REAL*8 (A-H,O-Z)
688      REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), &
689       AK2(N),AK3(N),AK4(N),FX(N),            &
690        FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1)
691      INTEGER IP(N)
692      LOGICAL REJECT,RJECT2,BANDED
693      REAL(wp), DIMENSION(2) :: RSTAT
694      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
695
696! ---- PREPARE BANDWIDTHS -----
697      IF (BANDED) THEN
698          MLE=MLJAC
699          MUE=MUJAC
700          MBJAC=MLJAC+MUJAC+1
701          MDIAG=MLE+MUE+1
702      END IF
703! *** *** *** *** *** *** ***
704!  INITIALISATIONS
705! *** *** *** *** *** *** ***
706      POSNEG=DSIGN(1.D0,XEND-X)
707      CALL ROS3PR(A21,A31,A32,C21,C31,C32,B1,B2,B3,E1,E2,E3,GAMMA)
708
709! --- INITIAL PREPARATIONS
710      HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
711      H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN)
712      H=DSIGN(H,POSNEG)
713      REJECT=.FALSE.
714      NSING=0
715      IRTRN=1
716      XXOLD=X
717
718      IF (IRTRN.LT.0) GOTO 79
719! --- BASIC INTEGRATION STEP
720   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
721      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
722          H=HOPT
723          IDID=1
724          RETURN
725      END IF
726      HOPT=H
727      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
728
729      CALL sed_func(N,JINDIC,X,Y,DY1)
730
731      NFCN=NFCN+1
732! *** *** *** *** *** *** ***
733!  COMPUTATION OF THE JACOBIAN
734! *** *** *** *** *** *** ***
735      NJAC=NJAC+1
736      IF (IJAC.EQ.0) THEN
737! --- COMPUTE JACOBIAN MATRIX NUMERICALLY
738          IF (BANDED) THEN
739! --- JACOBIAN IS BANDED
740              MUJACP=MUJAC+1
741              MD=MIN(MBJAC,N)
742              DO 16 K=1,MD
743              J=K
744 12           AK2(J)=Y(J)
745              AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
746              Y(J)=Y(J)+AK3(J)
747              J=J+MD
748              IF (J.LE.N) GOTO 12
749              CALL sed_func(N,JINDIC,X,Y,AK1)
750              J=K
751              LBEG=MAX(1,J-MUJAC)
752 14           LEND=MIN(N,J+MLJAC)
753              Y(J)=AK2(J)
754              MUJACJ=MUJACP-J
755              DO 15 L=LBEG,LEND
756 15           FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J)
757              J=J+MD
758              LBEG=LEND+1
759              IF (J.LE.N) GOTO 14
760 16           CONTINUE
761          ELSE
762! --- JACOBIAN IS FULL
763              DO 18 I=1,N
764              YSAFE=Y(I)
765              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
766              Y(I)=YSAFE+DELT
767              CALL sed_func(N,JINDIC,X,Y,AK1)
768              DO 17 J=1,N
769  17          FJAC(J,I)=(AK1(J)-DY1(J))/DELT
770  18          Y(I)=YSAFE
771              MLJAC=N
772          END IF
773      ELSE
774! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
775          CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC)
776      END IF
777   2  CONTINUE
778! *** *** *** *** *** *** ***
779!  COMPUTE THE STAGES
780! *** *** *** *** *** *** ***
781      NDEC=NDEC+1
782      HC21=C21/H
783      HC31=C31/H
784      HC32=C32/H
785      FAC=1.D0/(H*GAMMA)
786          IF (BANDED) THEN
787! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
788              DO 601 J=1,N
789              I1=MAX0(1,MUJAC+2-J)
790              I2=MIN0(MBJAC,N+MUJAC+1-J)
791              DO 600 I=I1,I2
792  600         E(I+MLE,J)=-FJAC(I,J)
793  601         E(MDIAG,J)=E(MDIAG,J)+FAC
794              CALL DECB(N,LE,E,MLE,MUE,IP,INFO)
795              IF (INFO.NE.0) GOTO 80
796! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
797! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
798! ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
799! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
800              DO 603 I=1,N
801  603         AK1(I)=DY1(I)
802              CALL SOLB(N,LE,E,MLE,MUE,AK1,IP)
803              DO 610 I=1,N
804  610         YNEW(I)=Y(I)+A21*AK1(I)
805              CALL sed_func(N,JINDIC,X,YNEW,DY)
806              DO 611 I=1,N
807  611         AK2(I)=DY(I)+HC21*AK1(I)
808              CALL SOLB(N,LE,E,MLE,MUE,AK2,IP)
809              DO 620 I=1,N
810  620         YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
811              CALL sed_func(N,JINDIC,X,YNEW,DY)
812              DO 621 I=1,N
813  621         AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
814              CALL SOLB(N,LE,E,MLE,MUE,AK3,IP)
815          ELSE
816! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
817              DO 801 J=1,N
818              DO 800 I=1,N
819  800         E(I,J)=-FJAC(I,J)
820  801         E(J,J)=E(J,J)+FAC
821              CALL DEC(N,LE,E,IP,INFO)
822              IF (INFO.NE.0) GOTO 80
823! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
824! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
825! ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
826! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
827              DO 803 I=1,N
828  803         AK1(I)=DY1(I)
829              CALL SOL(N,LE,E,AK1,IP)
830              DO 810 I=1,N
831  810         YNEW(I)=Y(I)+A21*AK1(I)
832              CALL sed_func(N,JINDIC,X,YNEW,DY)
833              DO 811 I=1,N
834  811         AK2(I)=DY(I)+HC21*AK1(I)
835              CALL SOL(N,LE,E,AK2,IP)
836              DO 820 I=1,N
837  820         YNEW(I)=Y(I)+A31*AK1(I)+A32*AK2(I)
838              CALL sed_func(N,JINDIC,X,YNEW,DY)
839              DO 821 I=1,N
840  821         AK3(I)=DY(I)+HC31*AK1(I)+HC32*AK2(I)
841              CALL SOL(N,LE,E,AK3,IP)
842          END IF
843      NSOL=NSOL+3
844      NFCN=NFCN+2
845! *** *** *** *** *** *** ***
846!  ERROR ESTIMATION
847! *** *** *** *** *** *** ***
848      NSTEP=NSTEP+1
849! ------------ NEW SOLUTION ---------------
850      DO 240 I=1,N
851  240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)+B3*AK3(I)
852! ------------ COMPUTE ERROR ESTIMATION ----------------
853      ERR=0.D0
854      DO 300 I=1,N
855      S=E1*AK1(I)+E2*AK2(I)+E3*AK3(I)
856      IF (ITOL.EQ.0) THEN
857         SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
858      ELSE
859         SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
860      END IF
861  300 ERR=ERR+(S/SK)**2
862      ERR=DSQRT(ERR/N)
863! --- COMPUTATION OF HNEW
864! --- WE REQUIRE .2<=HNEW/H<=6.
865      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.33D0/0.9D0))
866      HNEW=H/FAC
867! *** *** *** *** *** *** ***
868!  IS THE ERROR SMALL ENOUGH ?
869! *** *** *** *** *** *** ***
870      RJECT2 = .TRUE.
871      IF (ERR.LE.1.D0) THEN
872! --- STEP IS ACCEPTED
873         NACCPT=NACCPT+1
874         DO 44 I=1,N
875  44     Y(I)=YNEW(I)
876         XXOLD=X
877         X=X+H
878         RSTAT(2) = H
879         IF (IRTRN.LT.0) GOTO 79
880         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
881         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
882         REJECT=.FALSE.
883         RJECT2=.FALSE.
884         IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0
885         H=HNEW
886         GOTO 1
887      ELSE
888! --- STEP IS REJECTED
889         IF (RJECT2) HNEW=H*FACREJ
890         IF (REJECT) RJECT2=.TRUE.
891         REJECT=.TRUE.
892         H=HNEW
893         IF (NACCPT.GE.1) NREJCT=NREJCT+1
894         GOTO 2
895      END IF
896! --- EXIT
897  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
898      NSING=NSING+1
899      IF (NSING.GE.5) GOTO 79
900      H=H*0.5D0
901      GOTO 2
902  79  WRITE(NUMSED,979)X,H
903 979  FORMAT(' EXIT OF ROS4 AT X=',D16.7,'   H=',D16.7)
904      IDID=-1
905      RETURN
906
907      END SUBROUTINE RO3COR
908
909!
910      SUBROUTINE RO2COR(N,X,Y,XEND,HMAX,H,RTOL,ATOL,ITOL,JAC,         &
911        IJAC,MLJAC,MUJAC,IDID,               &
912        NMAX,UROUND,METH,FAC1,FAC2,FACREJ,BANDED,       &
913        LFJAC,LE,YNEW,DY1,DY,AK1,AK2,AK3,AK4,FX,FJAC,E,IP,RSTAT)
914! ----------------------------------------------------------
915!     CORE INTEGRATOR FOR ROS4
916!     PARAMETERS SAME AS IN ROS4 WITH WORKSPACE ADDED
917! ----------------------------------------------------------
918!         DECLARATIONS
919! ----------------------------------------------------------
920      IMPLICIT REAL*8 (A-H,O-Z)
921      REAL*8 Y(N),YNEW(N),DY1(N),DY(N),AK1(N), &
922            AK2(N),AK3(N),AK4(N),FX(N),            &
923        FJAC(LFJAC,N),E(LE,N),ATOL(1),RTOL(1)
924      INTEGER IP(N)
925      LOGICAL REJECT,RJECT2,BANDED
926      REAL(wp), DIMENSION(2) :: RSTAT
927      COMMON/STAT/NFCN,NJAC,NSTEP,NACCPT,NREJCT,NDEC,NSOL
928
929! ---- PREPARE BANDWIDTHS -----
930      IF (BANDED) THEN
931          MLE=MLJAC
932          MUE=MUJAC
933          MBJAC=MLJAC+MUJAC+1
934          MDIAG=MLE+MUE+1
935      END IF
936! *** *** *** *** *** *** ***
937!  INITIALISATIONS
938! *** *** *** *** *** *** ***
939      POSNEG=DSIGN(1.D0,XEND-X)
940      CALL RO2COEF (A21,C21,B1,B2,E1,E2,GAMMA)
941! --- INITIAL PREPARATIONS
942      HMAXN=DMIN1(DABS(HMAX),DABS(XEND-X))
943      H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAXN)
944      H=DSIGN(H,POSNEG)
945      REJECT=.FALSE.
946      NSING=0
947      IRTRN=1
948      XXOLD=X
949
950      IF (IRTRN.LT.0) GOTO 79
951! --- BASIC INTEGRATION STEP
952   1  IF (NSTEP.GT.NMAX.OR.X+.1D0*H.EQ.X.OR.ABS(H).LE.UROUND) GOTO 79
953      IF ((X-XEND)*POSNEG+UROUND.GT.0.D0) THEN
954          H=HOPT
955          IDID=1
956          RETURN
957      END IF
958      HOPT=H
959      IF ((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X
960
961      CALL sed_func(N,JINDIC,X,Y,DY1)
962
963      NFCN=NFCN+1
964! *** *** *** *** *** *** ***
965!  COMPUTATION OF THE JACOBIAN
966! *** *** *** *** *** *** ***
967      NJAC=NJAC+1
968      IF (IJAC.EQ.0) THEN
969! --- COMPUTE JACOBIAN MATRIX NUMERICALLY
970          IF (BANDED) THEN
971! --- JACOBIAN IS BANDED
972              MUJACP=MUJAC+1
973              MD=MIN(MBJAC,N)
974              DO 16 K=1,MD
975              J=K
976 12           AK2(J)=Y(J)
977              AK3(J)=DSQRT(UROUND*MAX(1.D-5,ABS(Y(J))))
978              Y(J)=Y(J)+AK3(J)
979              J=J+MD
980              IF (J.LE.N) GOTO 12
981              CALL sed_func(N,JINDIC,X,Y,AK1)
982              J=K
983              LBEG=MAX(1,J-MUJAC)
984 14           LEND=MIN(N,J+MLJAC)
985              Y(J)=AK2(J)
986              MUJACJ=MUJACP-J
987              DO 15 L=LBEG,LEND
988 15           FJAC(L+MUJACJ,J)=(AK1(L)-DY1(L))/AK3(J)
989              J=J+MD
990              LBEG=LEND+1
991              IF (J.LE.N) GOTO 14
992 16           CONTINUE
993          ELSE
994! --- JACOBIAN IS FULL
995              DO 18 I=1,N
996              YSAFE=Y(I)
997              DELT=DSQRT(UROUND*MAX(1.D-5,ABS(YSAFE)))
998              Y(I)=YSAFE+DELT
999              CALL sed_func(N,JINDIC,X,Y,AK1)
1000              DO 17 J=1,N
1001  17          FJAC(J,I)=(AK1(J)-DY1(J))/DELT
1002  18          Y(I)=YSAFE
1003              MLJAC=N
1004          END IF
1005      ELSE
1006! --- COMPUTE JACOBIAN MATRIX ANALYTICALLY
1007          CALL JAC(N,JINDIC,X,Y,FJAC,LFJAC)
1008      END IF
1009   2  CONTINUE
1010! *** *** *** *** *** *** ***
1011!  COMPUTE THE STAGES
1012! *** *** *** *** *** *** ***
1013      NDEC=NDEC+1
1014      HC21=C21/H
1015      FAC=1.D0/(H*GAMMA)
1016          IF (BANDED) THEN
1017! --- THE MATRIX E (B=IDENTITY, JACOBIAN A BANDED MATRIX)
1018              DO 601 J=1,N
1019              I1=MAX0(1,MUJAC+2-J)
1020              I2=MIN0(MBJAC,N+MUJAC+1-J)
1021              DO 600 I=I1,I2
1022  600         E(I+MLE,J)=-FJAC(I,J)
1023  601         E(MDIAG,J)=E(MDIAG,J)+FAC
1024              CALL DECB(N,LE,E,MLE,MUE,IP,INFO)
1025              IF (INFO.NE.0) GOTO 80
1026! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
1027! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
1028! ---   2) THE JACOBIAN OF THE PROBLEM IS A BANDED MATRIX
1029! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
1030              DO 603 I=1,N
1031  603         AK1(I)=DY1(I)
1032              CALL SOLB(N,LE,E,MLE,MUE,AK1,IP)
1033              DO 610 I=1,N
1034  610         YNEW(I)=Y(I)+A21*AK1(I)
1035              CALL sed_func(N,JINDIC,X,YNEW,DY)
1036              DO 611 I=1,N
1037  611         AK2(I)=DY(I)+HC21*AK1(I)
1038              CALL SOLB(N,LE,E,MLE,MUE,AK2,IP)
1039          ELSE
1040! --- THE MATRIX E (B=IDENTITY, JACOBIAN A FULL MATRIX)
1041              DO 801 J=1,N
1042              DO 800 I=1,N
1043  800         E(I,J)=-FJAC(I,J)
1044  801         E(J,J)=E(J,J)+FAC
1045              CALL DEC(N,LE,E,IP,INFO)
1046              IF (INFO.NE.0) GOTO 80
1047! --- THIS PART COMPUTES THE STAGES IN THE CASE WHERE
1048! ---   1) THE DIFFERENTIAL EQUATION IS IN EXPLICIT FORM
1049! ---   2) THE JACOBIAN OF THE PROBLEM IS A FULL MATRIX
1050! ---   3) THE DIFFERENTIAL EQUATION IS AUTONOMOUS
1051              DO 803 I=1,N
1052  803         AK1(I)=DY1(I)
1053              CALL SOL(N,LE,E,AK1,IP)
1054              DO 810 I=1,N
1055  810         YNEW(I)=Y(I)+A21*AK1(I)
1056              CALL sed_func(N,JINDIC,X,YNEW,DY)
1057              DO 811 I=1,N
1058  811         AK2(I)=DY(I)+HC21*AK1(I)
1059              CALL SOL(N,LE,E,AK2,IP)
1060          END IF
1061      NSOL=NSOL+2
1062      NFCN=NFCN+1
1063! *** *** *** *** *** *** ***
1064!  ERROR ESTIMATION
1065! *** *** *** *** *** *** ***
1066      NSTEP=NSTEP+1
1067! ------------ NEW SOLUTION ---------------
1068      DO 240 I=1,N
1069  240 YNEW(I)=Y(I)+B1*AK1(I)+B2*AK2(I)
1070! ------------ COMPUTE ERROR ESTIMATION ----------------
1071      ERR=0.D0
1072      DO 300 I=1,N
1073      S=E1*AK1(I)+E2*AK2(I)
1074      IF (ITOL.EQ.0) THEN
1075         SK=ATOL(1)+RTOL(1)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
1076      ELSE
1077         SK=ATOL(I)+RTOL(I)*DMAX1(DABS(Y(I)),DABS(YNEW(I)))
1078      END IF
1079  300 ERR=ERR+(S/SK)**2
1080      ERR=DSQRT(ERR/N)
1081! --- COMPUTATION OF HNEW
1082! --- WE REQUIRE .2<=HNEW/H<=6.
1083      FAC=DMAX1(FAC2,DMIN1(FAC1,(ERR)**.25D0/.9D0))
1084      HNEW=H/FAC
1085! *** *** *** *** *** *** ***
1086!  IS THE ERROR SMALL ENOUGH ?
1087! *** *** *** *** *** *** ***
1088      RJECT2 = .TRUE.
1089      IF (ERR.LE.1.D0) THEN
1090! --- STEP IS ACCEPTED
1091         NACCPT=NACCPT+1
1092         DO 44 I=1,N
1093  44     Y(I)=YNEW(I)
1094         XXOLD=X
1095         X=X+H
1096         RSTAT(2) = H
1097         IF (IRTRN.LT.0) GOTO 79
1098         IF (DABS(HNEW).GT.HMAXN) HNEW=POSNEG*HMAXN
1099         IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H))
1100         REJECT=.FALSE.
1101         RJECT2=.FALSE.
1102         IF (NACCPT == 1) RSTAT(1) = (H+HNEW)/2.0
1103         H=HNEW
1104         GOTO 1
1105      ELSE
1106! --- STEP IS REJECTED
1107         IF (RJECT2) HNEW=H*FACREJ
1108         IF (REJECT) RJECT2=.TRUE.
1109         REJECT=.TRUE.
1110         H=HNEW
1111         IF (NACCPT.GE.1) NREJCT=NREJCT+1
1112         GOTO 2
1113      END IF
1114! --- EXIT
1115  80  WRITE (NUMSED,*) ' MATRIX E IS SINGULAR, INFO = ',INFO
1116      NSING=NSING+1
1117      IF (NSING.GE.5) GOTO 79
1118      H=H*0.5D0
1119      GOTO 2
1120  79  WRITE(NUMSED,979)X,H
1121 979  FORMAT(' EXIT OF ROS2 AT X=',D16.7,'   H=',D16.7)
1122      IDID=-1
1123      RETURN
1124
1125      END SUBROUTINE RO2COR
1126
1127      SUBROUTINE SHAMP (A21,A31,A32,C21,C31,C32,C41,C42,C43,      &
1128                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1129      IMPLICIT REAL*8 (A-H,O-Z)
1130         A21=2.D0
1131         A31=48.D0/25.D0
1132         A32=6.D0/25.D0
1133         C21=-8.D0
1134         C31=372.D0/25.D0
1135         C32=12.D0/5.D0
1136         C41=-112.D0/125.D0
1137         C42=-54.D0/125.D0
1138         C43=-2.D0/5.D0
1139         B1=19.D0/9.D0
1140         B2=1.D0/2.D0
1141         B3=25.D0/108.D0
1142         B4=125.D0/108.D0
1143         E1=17.D0/54.D0
1144         E2=7.D0/36.D0
1145         E3=0.D0
1146         E4=125.D0/108.D0
1147         GAMMA=.5D0
1148      RETURN
1149      END SUBROUTINE SHAMP
1150!
1151      SUBROUTINE GRK4A (A21,A31,A32,C21,C31,C32,C41,C42,C43,     &
1152                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1153      IMPLICIT REAL*8 (A-H,O-Z)
1154       A21= 0.1108860759493671D+01
1155       A31= 0.2377085261983360D+01
1156       A32= 0.1850114988899692D+00
1157       C21=-0.4920188402397641D+01
1158       C31= 0.1055588686048583D+01
1159       C32= 0.3351817267668938D+01
1160       C41= 0.3846869007049313D+01
1161       C42= 0.3427109241268180D+01
1162       C43=-0.2162408848753263D+01
1163       B1= 0.1845683240405840D+01
1164       B2= 0.1369796894360503D+00
1165       B3= 0.7129097783291559D+00
1166       B4= 0.6329113924050632D+00
1167       E1= 0.4831870177201765D-01
1168       E2=-0.6471108651049505D+00
1169       E3= 0.2186876660500240D+00
1170       E4=-0.6329113924050632D+00
1171       GAMMA= 0.3950000000000000D+00
1172      RETURN
1173
1174      END SUBROUTINE GRK4A
1175!
1176      SUBROUTINE GRK4T (A21,A31,A32,C21,C31,C32,C41,C42,C43,     &
1177                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1178      IMPLICIT REAL*8 (A-H,O-Z)
1179       A21= 0.2000000000000000D+01
1180       A31= 0.4524708207373116D+01
1181       A32= 0.4163528788597648D+01
1182       C21=-0.5071675338776316D+01
1183       C31= 0.6020152728650786D+01
1184       C32= 0.1597506846727117D+00
1185       C41=-0.1856343618686113D+01
1186       C42=-0.8505380858179826D+01
1187       C43=-0.2084075136023187D+01
1188       B1= 0.3957503746640777D+01
1189       B2= 0.4624892388363313D+01
1190       B3= 0.6174772638750108D+00
1191       B4= 0.1282612945269037D+01
1192       E1= 0.2302155402932996D+01
1193       E2= 0.3073634485392623D+01
1194       E3=-0.8732808018045032D+00
1195       E4=-0.1282612945269037D+01
1196       GAMMA= 0.2310000000000000D+00
1197      RETURN
1198
1199      END SUBROUTINE GRK4T
1200!
1201      SUBROUTINE VELDS (A21,A31,A32,C21,C31,C32,C41,C42,C43,     &
1202                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1203! --- METHOD GIVEN BY VAN VELDHUIZEN
1204      IMPLICIT REAL*8 (A-H,O-Z)
1205       A21= 0.2000000000000000D+01
1206       A31= 0.1750000000000000D+01
1207       A32= 0.2500000000000000D+00
1208       C21=-0.8000000000000000D+01
1209       C31=-0.8000000000000000D+01
1210       C32=-0.1000000000000000D+01
1211       C41= 0.5000000000000000D+00
1212       C42=-0.5000000000000000D+00
1213       C43= 0.2000000000000000D+01
1214       B1= 0.1333333333333333D+01
1215       B2= 0.6666666666666667D+00
1216       B3=-0.1333333333333333D+01
1217       B4= 0.1333333333333333D+01
1218       E1=-0.3333333333333333D+00
1219       E2=-0.3333333333333333D+00
1220       E3=-0.0000000000000000D+00
1221       E4=-0.1333333333333333D+01
1222       GAMMA= 0.5000000000000000D+00
1223      RETURN
1224      END SUBROUTINE VELDS
1225!
1226      SUBROUTINE VELDD (A21,A31,A32,C21,C31,C32,C41,C42,C43,     &
1227                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1228! --- METHOD GIVEN BY VAN VELDHUIZEN
1229      IMPLICIT REAL*8 (A-H,O-Z)
1230       A21= 0.2000000000000000D+01
1231       A31= 0.4812234362695436D+01
1232       A32= 0.4578146956747842D+01
1233       C21=-0.5333333333333331D+01
1234       C31= 0.6100529678848254D+01
1235       C32= 0.1804736797378427D+01
1236       C41=-0.2540515456634749D+01
1237       C42=-0.9443746328915205D+01
1238       C43=-0.1988471753215993D+01
1239       B1= 0.4289339254654537D+01
1240       B2= 0.5036098482851414D+01
1241       B3= 0.6085736420673917D+00
1242       B4= 0.1355958941201148D+01
1243       E1= 0.2175672787531755D+01
1244       E2= 0.2950911222575741D+01
1245       E3=-0.7859744544887430D+00
1246       E4=-0.1355958941201148D+01
1247       GAMMA= 0.2257081148225682D+00
1248      RETURN
1249      END SUBROUTINE VELDD
1250!
1251      SUBROUTINE LSTAB (A21,A31,A32,C21,C31,C32,C41,C42,C43,     &
1252                B1,B2,B3,B4,E1,E2,E3,E4,GAMMA)
1253! --- AN L-STABLE METHOD
1254      IMPLICIT REAL*8 (A-H,O-Z)
1255       A21= 0.2000000000000000D+01
1256       A31= 0.1867943637803922D+01
1257       A32= 0.2344449711399156D+00
1258       C21=-0.7137615036412310D+01
1259       C31= 0.2580708087951457D+01
1260       C32= 0.6515950076447975D+00
1261       C41=-0.2137148994382534D+01
1262       C42=-0.3214669691237626D+00
1263       C43=-0.6949742501781779D+00
1264       B1= 0.2255570073418735D+01
1265       B2= 0.2870493262186792D+00
1266       B3= 0.4353179431840180D+00
1267       B4= 0.1093502252409163D+01
1268       E1=-0.2815431932141155D+00
1269       E2=-0.7276199124938920D-01
1270       E3=-0.1082196201495311D+00
1271       E4=-0.1093502252409163D+01
1272       GAMMA= 0.5728200000000000D+00
1273      RETURN
1274
1275      END SUBROUTINE LSTAB
1276
1277      SUBROUTINE ROS3PR(A21,A31,A32,C21,C31,C32,B1,B2,B3,E1,E2,E3,GAMMA)
1278      IMPLICIT REAL*8 (A-H,O-Z)
1279         A21=1.0
1280         A31=1.0
1281         A32=0.0
1282         C21=-0.10156171083877702091975600115545E+01
1283         C31= 0.40759956452537699824805835358067E+01
1284         C32= 0.92076794298330791242156818474003E+01
1285         B1= 0.1E+01
1286         B2= 0.61697947043828245592553615689730E+01
1287         B3= -0.42772256543218573326238373806514
1288         E1= 0.5
1289         E2=-0.29079558716805469821718236208017E+01
1290         E3= 0.22354069897811569627360909276199
1291         GAMMA=0.43586652150845899941601945119356
1292      RETURN
1293      END SUBROUTINE ROS3PR
1294
1295      SUBROUTINE RO2COEF (A21,C21,B1,B2,E1,E2,GAMMA)
1296      IMPLICIT REAL*8 (A-H,O-Z)
1297         A21=1.D0
1298         C21=-2.D0
1299         B1=3./2.
1300         B2=1./2.
1301         E1=1./2.
1302         E2=1./2.
1303         GAMMA=1.0 + 1.0/SQRT(2.)
1304      RETURN
1305      END SUBROUTINE RO2COEF
1306!
1307      SUBROUTINE DEC (N, NDIM, A, IP, IER)
1308! VERSION REAL DOUBLE PRECISION
1309      INTEGER N,NDIM,IP,IER,NM1,K,KP1,M,I,J
1310      DOUBLE PRECISION A,T
1311      DIMENSION A(NDIM,N), IP(N)
1312!-----------------------------------------------------------------------
1313!  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION.
1314!  INPUT..
1315!     N = ORDER OF MATRIX.
1316!     NDIM = DECLARED DIMENSION OF ARRAY  A .
1317!     A = MATRIX TO BE TRIANGULARIZED.
1318!  OUTPUT..
1319!     A(I,J), I.LE.J = UPPER TRIANGULAR FACTOR, U .
1320!     A(I,J), I.GT.J = MULTIPLIERS = LOWER TRIANGULAR FACTOR, I - L.
1321!     IP(K), K.LT.N = INDEX OF K-TH PIVOT ROW.
1322!     IP(N) = (-1)**(NUMBER OF INTERCHANGES) OR O .
1323!     IER = 0 IF MATRIX A IS NONSINGULAR, OR K IF FOUND TO BE
1324!           SINGULAR AT STAGE K.
1325!  USE  SOL  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
1326!  DETERM(A) = IP(N)*A(1,1)*A(2,2)*...*A(N,N).
1327!  IF IP(N)=O, A IS SINGULAR, SOL WILL DIVIDE BY ZERO.
1328!
1329!  REFERENCE..
1330!     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
1331!     C.A.C.M. 15 (1972), P. 274.
1332!-----------------------------------------------------------------------
1333      IER = 0
1334      IP(N) = 1
1335      IF (N .EQ. 1) GO TO 70
1336      NM1 = N - 1
1337      DO 60 K = 1,NM1
1338        KP1 = K + 1
1339        M = K
1340        DO 10 I = KP1,N
1341          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
1342 10     CONTINUE
1343        IP(K) = M
1344        T = A(M,K)
1345        IF (M .EQ. K) GO TO 20
1346        IP(N) = -IP(N)
1347        A(M,K) = A(K,K)
1348        A(K,K) = T
1349 20     CONTINUE
1350        IF (T .EQ. 0.D0) GO TO 80
1351        T = 1.D0/T
1352        DO 30 I = KP1,N
1353 30       A(I,K) = -A(I,K)*T
1354        DO 50 J = KP1,N
1355          T = A(M,J)
1356          A(M,J) = A(K,J)
1357          A(K,J) = T
1358          IF (T .EQ. 0.D0) GO TO 45
1359          DO 40 I = KP1,N
1360 40         A(I,J) = A(I,J) + A(I,K)*T
1361 45       CONTINUE
1362 50       CONTINUE
1363 60     CONTINUE
1364 70   K = N
1365      IF (A(N,N) .EQ. 0.D0) GO TO 80
1366      RETURN
1367 80   IER = K
1368      IP(N) = 0
1369      RETURN
1370!----------------------- END OF SUBROUTINE DEC -------------------------
1371      END SUBROUTINE DEC
1372
1373      SUBROUTINE SOL (N, NDIM, A, B, IP)
1374! VERSION REAL DOUBLE PRECISION
1375      INTEGER N,NDIM,IP,NM1,K,KP1,M,I,KB,KM1
1376      DOUBLE PRECISION A,B,T
1377      DIMENSION A(NDIM,N), B(N), IP(N)
1378!-----------------------------------------------------------------------
1379!  SOLUTION OF LINEAR SYSTEM, A*X = B .
1380!  INPUT..
1381!    N = ORDER OF MATRIX.
1382!    NDIM = DECLARED DIMENSION OF ARRAY  A .
1383!    A = TRIANGULARIZED MATRIX OBTAINED FROM DEC.
1384!    B = RIGHT HAND SIDE VECTOR.
1385!    IP = PIVOT VECTOR OBTAINED FROM DEC.
1386!  DO NOT USE IF DEC HAS SET IER <> 0.
1387!  OUTPUT..
1388!    B = SOLUTION VECTOR, X .
1389!-----------------------------------------------------------------------
1390      IF (N .EQ. 1) GO TO 50
1391      NM1 = N - 1
1392      DO 20 K = 1,NM1
1393        KP1 = K + 1
1394        M = IP(K)
1395        T = B(M)
1396        B(M) = B(K)
1397        B(K) = T
1398        DO 10 I = KP1,N
1399 10       B(I) = B(I) + A(I,K)*T
1400 20     CONTINUE
1401      DO 40 KB = 1,NM1
1402        KM1 = N - KB
1403        K = KM1 + 1
1404        B(K) = B(K)/A(K,K)
1405        T = -B(K)
1406        DO 30 I = 1,KM1
1407 30       B(I) = B(I) + A(I,K)*T
1408 40     CONTINUE
1409 50   B(1) = B(1)/A(1,1)
1410      RETURN
1411!----------------------- END OF SUBROUTINE SOL -------------------------
1412      END SUBROUTINE SOL
1413
1414      SUBROUTINE DECB (N, NDIM, A, ML, MU, IP, IER)
1415      REAL*8 A,T
1416      DIMENSION A(NDIM,N), IP(N)
1417!-----------------------------------------------------------------------
1418!  MATRIX TRIANGULARIZATION BY GAUSSIAN ELIMINATION OF A BANDED
1419!  MATRIX WITH LOWER BANDWIDTH ML AND UPPER BANDWIDTH MU
1420!  INPUT..
1421!     N       ORDER OF THE ORIGINAL MATRIX A.
1422!     NDIM    DECLARED DIMENSION OF ARRAY  A.
1423!     A       CONTAINS THE MATRIX IN BAND STORAGE.   THE COLUMNS
1424!                OF THE MATRIX ARE STORED IN THE COLUMNS OF  A  AND
1425!                THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
1426!                ML+1 THROUGH 2*ML+MU+1 OF  A.
1427!     ML      LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1428!     MU      UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1429!  OUTPUT..
1430!     A       AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
1431!                THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
1432!     IP      INDEX VECTOR OF PIVOT INDICES.
1433!     IP(N)   (-1)**(NUMBER OF INTERCHANGES) OR O .
1434!     IER     = 0 IF MATRIX A IS NONSINGULAR, OR  = K IF FOUND TO BE
1435!                SINGULAR AT STAGE K.
1436!  USE  SOLB  TO OBTAIN SOLUTION OF LINEAR SYSTEM.
1437!  DETERM(A) = IP(N)*A(MD,1)*A(MD,2)*...*A(MD,N)  WITH MD=ML+MU+1.
1438!  IF IP(N)=O, A IS SINGULAR, SOLB WILL DIVIDE BY ZERO.
1439!
1440!  REFERENCE..
1441!     THIS IS A MODIFICATION OF
1442!     C. B. MOLER, ALGORITHM 423, LINEAR EQUATION SOLVER,
1443!     C.A.C.M. 15 (1972), P. 274.
1444!-----------------------------------------------------------------------
1445      IER = 0
1446      IP(N) = 1
1447      MD = ML + MU + 1
1448      MD1 = MD + 1
1449      JU = 0
1450      IF (ML .EQ. 0) GO TO 70
1451      IF (N .EQ. 1) GO TO 70
1452      IF (N .LT. MU+2) GO TO 7
1453      DO 5 J = MU+2,N
1454      DO 5 I = 1,ML
1455  5   A(I,J) = 0.D0
1456  7   NM1 = N - 1
1457      DO 60 K = 1,NM1
1458        KP1 = K + 1
1459        M = MD
1460        MDL = MIN(ML,N-K) + MD
1461        DO 10 I = MD1,MDL
1462          IF (DABS(A(I,K)) .GT. DABS(A(M,K))) M = I
1463 10     CONTINUE
1464        IP(K) = M + K - MD
1465        T = A(M,K)
1466        IF (M .EQ. MD) GO TO 20
1467        IP(N) = -IP(N)
1468        A(M,K) = A(MD,K)
1469        A(MD,K) = T
1470 20     CONTINUE
1471        IF (T .EQ. 0.D0) GO TO 80
1472        T = 1.D0/T
1473        DO 30 I = MD1,MDL
1474 30       A(I,K) = -A(I,K)*T
1475        JU = MIN0(MAX0(JU,MU+IP(K)),N)
1476        MM = MD
1477        IF (JU .LT. KP1) GO TO 55
1478        DO 50 J = KP1,JU
1479          M = M - 1
1480          MM = MM - 1
1481          T = A(M,J)
1482          IF (M .EQ. MM) GO TO 35
1483          A(M,J) = A(MM,J)
1484          A(MM,J) = T
1485 35       CONTINUE
1486          IF (T .EQ. 0.D0) GO TO 45
1487          JK = J - K
1488          DO 40 I = MD1,MDL
1489            IJK = I - JK
1490 40         A(IJK,J) = A(IJK,J) + A(I,K)*T
1491 45       CONTINUE
1492 50       CONTINUE
1493 55     CONTINUE
1494 60     CONTINUE
1495 70   K = N
1496      IF (A(MD,N) .EQ. 0.D0) GO TO 80
1497      RETURN
1498 80   IER = K
1499      IP(N) = 0
1500      RETURN
1501!----------------------- END OF SUBROUTINE DECB ------------------------
1502      END SUBROUTINE DECB
1503
1504      SUBROUTINE SOLB (N, NDIM, A, ML, MU, B, IP)
1505      REAL*8 A,B,T
1506      DIMENSION A(NDIM,N), B(N), IP(N)
1507      INTEGER :: IER
1508!-----------------------------------------------------------------------
1509!  SOLUTION OF LINEAR SYSTEM, A*X = B .
1510!  INPUT..
1511!    N      ORDER OF MATRIX A.
1512!    NDIM   DECLARED DIMENSION OF ARRAY  A .
1513!    A      TRIANGULARIZED MATRIX OBTAINED FROM DECB.
1514!    ML     LOWER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1515!    MU     UPPER BANDWIDTH OF A (DIAGONAL IS NOT COUNTED).
1516!    B      RIGHT HAND SIDE VECTOR.
1517!    IP     PIVOT VECTOR OBTAINED FROM DECB.
1518!  DO NOT USE IF DECB HAS SET IER  <> 0.
1519!  OUTPUT..
1520!    B      SOLUTION VECTOR, X .
1521!-----------------------------------------------------------------------
1522      MD = ML + MU + 1
1523      MD1 = MD + 1
1524      MDM = MD - 1
1525      NM1 = N - 1
1526      IF (ML .EQ. 0) GO TO 25
1527      IF (N .EQ. 1) GO TO 50
1528      DO 20 K = 1,NM1
1529        M = IP(K)
1530        T = B(M)
1531        B(M) = B(K)
1532        B(K) = T
1533        MDL = MIN(ML,N-K) + MD
1534        DO 10 I = MD1,MDL
1535          IMD = I + K - MD
1536 10       B(IMD) = B(IMD) + A(I,K)*T
1537 20     CONTINUE
1538 25   CONTINUE
1539      DO 40 KB = 1,NM1
1540        K = N + 1 - KB
1541        B(K) = B(K)/A(MD,K)
1542        T = -B(K)
1543        KMD = MD - K
1544        LM = MAX0(1,KMD+1)
1545        DO 30 I = LM,MDM
1546          IMD = I - KMD
1547 30       B(IMD) = B(IMD) + A(I,K)*T
1548 40     CONTINUE
1549 50   B(1) = B(1)/A(MD,1)
1550      RETURN
1551!----------------------- END OF SUBROUTINE SOLB ------------------------
1552      END SUBROUTINE SOLB
1553#endif
1554END MODULE trosk
Note: See TracBrowser for help on using the repository browser.