1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
|
SUBROUTINE DRV000(NENTRY,HENTRY,IENTRY,JENTRY,KENTRY)
*
*-----------------------------------------------------------------------
*
* MODULE TO FIND A ZERO USING BRENT'S METHOD WITH
* GIVEN FUNCTION VALUES.
*
* INPUT/OUTPUT PARAMETERS:
* NENTRY : NUMBER OF LINKED LISTS AND FILES USED BY THE MODULE.
* HENTRY : CHARACTER*12 NAME OF EACH LINKED LIST OR FILE.
* IENTRY : =0 CLE-2000 VARIABLE; =1 LINKED LIST; =2 XSM FILE;
* =3 SEQUENTIAL BINARY FILE; =4 SEQUENTIAL ASCII FILE;
* =5 DIRECT ACCESS FILE.
* JENTRY : =0 THE LINKED LIST OR FILE IS CREATED.
* =1 THE LINKED LIST OR FILE IS OPEN FOR MODIFICATIONS;
* =2 THE LINKED LIST OR FILE IS OPEN IN READ-ONLY MODE.
* KENTRY : =FILE UNIT NUMBER; =LINKED LIST ADDRESS OTHERWISE.
* DIMENSION HENTRY(NENTRY),IENTRY(NENTRY),JENTRY(NENTRY),
* KENTRY(NENTRY)
*
*-------------------------------------- AUTHOR: R. ROY ; 29/11/94 ---
*
USE GANLIB
IMPLICIT NONE
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NENTRY
CHARACTER HENTRY(NENTRY)*12
INTEGER IENTRY(NENTRY), JENTRY(NENTRY)
TYPE(C_PTR) KENTRY(NENTRY)
*----
* LOCAL VARIABLES
*----
INTEGER NITMA, ITYP, I
CHARACTER TEXT12*12, SGNTUR*12
LOGICAL LSTART, LCONV
DOUBLE PRECISION DFLOTT
INTEGER ITER,ITMAX,IPRT, ISGNTR(3),ITMD
INTEGER ITERV(3), ICONV
REAL A,B,C, D,E, FA,FB,FC, P,Q,R,S, TOL1,XM,TOL
REAL X(3), DE(2), Y(3), PQRS(4), ATOL(3)
REAL FLOTT
REAL EPM,TOLD,Z0,ZH,Z1,Z2,Z3,ZBESTM
TYPE(C_PTR) IPL0
PARAMETER (EPM=3.E-8,TOLD=1.E-5,ITMD=100)
PARAMETER (Z0=0.0,ZH=0.5,Z1=1.0,Z2=2.0,Z3=3.0)
*----
* PARAMETER VALIDATION.
*----
IF(NENTRY.NE.1) CALL XABORT('DRV000: ONLY ONE ENTRY EXPECTED.')
TEXT12=HENTRY(1)
IF(IENTRY(1).NE.1) CALL XABORT('DRV000: LHS L_0 OBJECT EXPECTED ('
> //TEXT12//').')
IF((JENTRY(1).NE.0).AND.(JENTRY(1).NE.1)) CALL XABORT('DRV000: LH'
> //'S L_0 OBJECT IN CREATE OR MODIFICATION MODE EXPECTED.')
*
LSTART= JENTRY(1).EQ.0
IPL0= KENTRY(1)
IF( LSTART )THEN
* DEFINE ALL TEMP VARIABLES
D= 0.0
E= 0.0
P= 0.0
Q= 0.0
R= 0.0
S= 0.0
*
LCONV= .FALSE.
TOL= TOLD
ITMAX= ITMD
ITER= 0
IPRT= 0
10 CONTINUE
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF( ITYP.NE.3 )
> CALL XABORT('DRV000: KEYWORDS *TOL*, *POINT*... EXPECTED.')
IF( TEXT12.EQ.'TOL' )THEN
CALL REDGET(ITYP,NITMA,TOL,TEXT12,DFLOTT)
IF( ITYP.NE.2 )
> CALL XABORT('DRV000: A REAL TOLERANCE *TOL* IS EXPECTED.')
IF( TOL.LT.1.E-7 )
> CALL XABORT('DRV000: TOLERANCE .LT. 1.E-7.')
GO TO 10
ELSEIF( TEXT12.EQ.'ITMAX' )THEN
CALL REDGET(ITYP,ITMAX,FLOTT,TEXT12,DFLOTT)
IF( ITYP.NE.1 )
> CALL XABORT('DRV000: AN INTEGER *ITMAX* IS EXPECTED.')
GO TO 10
ELSEIF( TEXT12.EQ.'DEBUG' )THEN
IPRT= 1
GO TO 10
ELSEIF( TEXT12.EQ.'POINT' )THEN
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'X')
> CALL XABORT('DRV000: *X* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,A ,TEXT12,DFLOTT)
IF(ITYP.NE.2 ) CALL XABORT('DRV000: AFTER *X*,'
> // ' A REAL IS EXPECTED.')
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'Y')
> CALL XABORT('DRV000: *Y* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,FA ,TEXT12,DFLOTT)
IF(ITYP.NE.2 ) CALL XABORT('DRV000: AFTER *Y*,'
> // ' A REAL IS EXPECTED.')
ELSE
CALL XABORT('DRV000: KEYWORDS *TOL* OR *POINT* EXPECTED.')
ENDIF
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'POINT')
> CALL XABORT('DRV000: ONCE MORE, *POINT* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'X')
> CALL XABORT('DRV000: *X* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,B ,TEXT12,DFLOTT)
IF(ITYP.NE.2 ) CALL XABORT('DRV000: AFTER *X*,'
> // ' A REAL IS EXPECTED.')
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'Y')
> CALL XABORT('DRV000: *Y* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,FB ,TEXT12,DFLOTT)
IF(ITYP.NE.2 ) CALL XABORT('DRV000: AFTER *Y*,'
> // ' A REAL IS EXPECTED.')
SGNTUR='L_0'
READ(SGNTUR,'(3A4)') ISGNTR
CALL LCMSIX(IPL0,' ',0)
*
* PUT SIGNATURE
CALL LCMPUT(IPL0,'SIGNATURE',3,3,ISGNTR)
*
* PUT CONVERGENCE FLAG
ICONV=-1
CALL LCMPUT(IPL0,'ICONV',1,1,ICONV)
ELSE
*
CALL LCMSIX(IPL0,' ',0)
*
* VERIFY SIGNATURE
CALL LCMGET(IPL0,'SIGNATURE',ISGNTR)
WRITE(SGNTUR,'(3A4)') (ISGNTR(I),I=1,3)
IF(SGNTUR.NE.'L_0')
> CALL XABORT('DRV000: L_0 OBJECT IS EXPECTED')
*
CALL LCMGET(IPL0,'ICONV',ICONV)
LCONV= ICONV.EQ.+1
*
* NOTIFY USER IF ALREADY CONVERVED
IF( LCONV )
> CALL XABORT('DRV000: PROCESS IS ALREADY CONVERGED')
*
* GET L_0 OBJECT VALUES
CALL LCMGET(IPL0,'X',X)
A=X(1)
B=X(2)
C=X(3)
CALL LCMGET(IPL0,'DE',DE)
D=DE(1)
E=DE(2)
CALL LCMGET(IPL0,'Y',Y)
FA=Y(1)
FB=Y(2)
FC=Y(3)
CALL LCMGET(IPL0,'PQRS',PQRS)
P=PQRS(1)
Q=PQRS(2)
R=PQRS(3)
S=PQRS(4)
CALL LCMGET(IPL0,'ATOL',ATOL)
TOL=ATOL(1)
XM=ATOL(2)
TOL1=ATOL(3)
CALL LCMGET(IPL0,'ITERV',ITERV)
ITER=ITERV(1)
ITMAX=ITERV(2)
IPRT=ITERV(3)
*
* GET NEW *Y* VALUE
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.'Y')
> CALL XABORT('DRV000: *Y* KEYWORD EXPECTED.')
CALL REDGET(ITYP,NITMA,FB ,TEXT12,DFLOTT)
IF(ITYP.NE.2 )
> CALL XABORT('DRV000: AFTER *Y*, A REAL IS EXPECTED.')
ENDIF
*
* METHOD: BRENT'S METHOD FOR FINDING ZEROS.
* FIRST INTERVAL MUST BE BRACKETED: FA*FB < 0.
*
* INPUT: ITER= NUMBER OF ITERATIONS (0 AT START)
* TOL= TOLERANCE FOR ZERO FINDING
* (A,FA)= FIRST POINT
* (B,FB)= SECOND POINT
*
* OUTPUT: ZBESTM= ESTIMATION OF NEXT ZERO
*
IF( ITER.EQ.0 )THEN
IF((FA.GT.Z0.AND.FB.GT.Z0).OR.(FA.LT.Z0.AND.FB.LT.Z0))
> CALL XABORT(' DRV000: ROOT MUST BE BRACKETED')
C=B
FC=FB
ENDIF
IF((FB.GT.Z0.AND.FC.GT.Z0).OR.(FB.LT.Z0.AND.FC.LT.Z0))THEN
C=A
FC=FA
D=B-A
E=D
ENDIF
IF(ABS(FC).LT.ABS(FB)) THEN
A=B
B=C
C=A
FA=FB
FB=FC
FC=FA
ENDIF
TOL1=Z2*EPM*ABS(B)+ ZH*TOL
XM=ZH*(C-B)
IF(ABS(XM).LE.TOL1 .OR. FB.EQ.Z0)THEN
ZBESTM=B
LCONV= .TRUE.
GO TO 20
ENDIF
IF(ABS(E).GE.TOL1 .AND. ABS(FA).GT.ABS(FB)) THEN
S=FB/FA
IF(A.EQ.C) THEN
P=Z2*XM*S
Q=Z1-S
ELSE
Q=FA/FC
R=FB/FC
P=S*(Z2*XM*Q*(Q-R)-(B-A)*(R-Z1))
Q=(Q-Z1)*(R-Z1)*(S-Z1)
ENDIF
IF(P.GT.Z0) Q=-Q
P=ABS(P)
IF(Z2*P .LT. MIN(Z3*XM*Q-ABS(TOL1*Q),ABS(E*Q))) THEN
E=D
D=P/Q
ELSE
D=XM
E=D
ENDIF
ELSE
D=XM
E=D
ENDIF
A=B
FA=FB
IF(ABS(D) .GT. TOL1) THEN
B=B+D
ELSE
B=B+SIGN(TOL1,XM)
ENDIF
ZBESTM=B
ITER= ITER + 1
IF( ITER.GT.ITMAX )
> CALL XABORT('DRV000: MAX NUMBER OF ITERATIONS REACHED.')
*
20 CONTINUE
*
* PUT L_0 OBJECT VALUES
X(1)=A
X(2)=B
X(3)=C
CALL LCMPUT(IPL0,'X',3,2,X)
DE(1)=D
DE(2)=E
CALL LCMPUT(IPL0,'DE',2,2,DE)
Y(1)=FA
Y(2)=FB
Y(3)=FC
CALL LCMPUT(IPL0,'Y',3,2,Y)
PQRS(1)=P
PQRS(2)=Q
PQRS(3)=R
PQRS(4)=S
CALL LCMPUT(IPL0,'PQRS',4,2,PQRS)
ATOL(1)=TOL
ATOL(2)=XM
ATOL(3)=TOL1
CALL LCMPUT(IPL0,'ATOL',3,2,ATOL)
ITERV(1)=ITER
ITERV(2)=ITMAX
ITERV(3)=IPRT
CALL LCMPUT(IPL0,'ITERV',3,1,ITERV)
IF( LCONV )THEN
ICONV=+1
ELSE
ICONV=-1
ENDIF
*
* SAVE CONVERGENCE FLAG
CALL LCMPUT(IPL0,'ICONV',1,1,ICONV)
CALL LCMSIX(IPL0,' ',0)
IF( IPRT.EQ.1 )THEN
WRITE(6,*) 'DEBUG: A=', A,' B=', B,' C=', C
WRITE(6,*) 'DEBUG: FA=',FA,' FB=',FB,' FC=',FC
WRITE(6,*) 'DEBUG: D=', D,' E=', E
WRITE(6,*) 'DEBUG: P=', P,' Q=', Q,' R=',R,' S=',S
WRITE(6,*) 'DEBUG: TOL1=',TOL1,' XM=',XM,' TOL=',TOL
WRITE(6,*) 'DEBUG: ITER=',ITER,' ITMAX=',ITMAX
WRITE(6,*) 'DEBUG: ICONV=',ICONV
ENDIF
*
* NOW, RETURN BACK LOGICAL VALUE AND ZERO ESTIMATE
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
ITYP= -ITYP
IF( ITYP.NE.5 )THEN
CALL XABORT('DRV000: MUST WRITE LOGICAL FLAG INTO >>.<<')
ENDIF
CALL REDPUT(ITYP,ICONV,FLOTT,TEXT12,DFLOTT)
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
ITYP= -ITYP
IF( ITYP.NE.2 )THEN
CALL XABORT('DRV000: MUST WRITE REAL ZERO INTO >>.<<')
ENDIF
CALL REDPUT(ITYP,NITMA,ZBESTM,TEXT12,DFLOTT)
*
CALL REDGET(ITYP,NITMA,FLOTT,TEXT12,DFLOTT)
IF(ITYP.NE.3.OR.TEXT12.NE.';')
> CALL XABORT('DRV000: *;* IS EXPECTED FOR ENDING THE SENTENCE.')
RETURN
END
|
