summaryrefslogtreecommitdiff
path: root/Ganlib/src/DRV000.f
diff options
context:
space:
mode:
Diffstat (limited to 'Ganlib/src/DRV000.f')
-rw-r--r--Ganlib/src/DRV000.f323
1 files changed, 323 insertions, 0 deletions
diff --git a/Ganlib/src/DRV000.f b/Ganlib/src/DRV000.f
new file mode 100644
index 0000000..7e19563
--- /dev/null
+++ b/Ganlib/src/DRV000.f
@@ -0,0 +1,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
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:16:02 +0000