在 Fortran 中查找复根
Finding complex root in Fortran
我发现了一些资源，其中建议使用算法来确定一般分析函数（无多项式）中的复根。例如。 “Muller 求解非线性方程的 Fortran 程序
方法”，由 I. Barrodale 和 K. B. Wilson 编写。不幸的是，该代码仅以 pdf 格式提供，并且转换为 ASCII 非常麻烦……还有 Henning Bach，第 121 卷，第 12 期，第 675-678 页，1969 年 12 月. 它在 ASCII 中可用。字母 1 工作正常但对我的目的来说相当慢。在 Fortran 中是否有计算复数根的标准例程？
谢谢 :)
对于那些有兴趣的人。我从 I. Barrodale 和 K. B. Wilson 那里拿了一张。由于这篇文章是免费在线阅读的，所以我想 post 这里的代码。这是迄今为止我发现的最快和最可靠的。要使用它，请定义一个子程序 calcf(x,f)。 f 是 x 的函数值。 x 和 f 是复数。
了解更多信息 -> http://www.sciencedirect.com/science/article/pii/0771050X78900414
这是双精度版本：
SUBROUTINE ROOTS(KNOWN,NMORE,KM,REALRT,EP1,EP2,GUESS,MAXITS,RTS,NFOUND, FNV, IFAIL)
COMPLEX*16 RTS(KM), FNV(KM), X1, X2, X3, F1, F2, F3, &
X21, X32, F21, F32, F31, F321, B, XNEW, FNEW, DENOM, &
FSAVE, CSQRT, RADICL, RT, DIF, CSTEP, FSLAST
REAL*8 EP1, EP2, ONE, HUNDRD, TWO, FOUR, HALF, STEP, &
ZERO, EP1DEF, EP2DEF
INTEGER GUESS
LOGICAL REALRT
DATA IZERO,IONE,ITWO,ITHREE,ZERO,HALF,ONE,TWO,FOUR,HUNDRD,ITSDEF,EP1DEF,EP2DEF &
/0,1,2,3,0.0d0,0.5d0,1.0d0,2.0d0,4.0d0,100.0d0,100,0.5d-5,1.0d-6/
EP1 = AMAX1(EP1,EP1DEF)
EP2 = AMAX1(EP2,EP2DEF)
MAXITS = MIN0(MAXITS,ITSDEF)
IFAIL = IZERO
NFOUND = IZERO
IF (KNOWN.LT.IONE) GO TO 30
DO 10 I=1,KNOWN
II = I
CALL CALCF(RTS(I), FNV(I))
IF (ABS(FNV(II)).GE.EP2) GO TO 20
10 CONTINUE
20 IFAIL =IONE
GUESS = GUESS + KNOWN - II +IONE
NMORE = NMORE + KNOWN - II +IONE
KNOWN = II - IONE
30 CONTINUE
LOOP1 = KNOWN +IONE
LOOP2 = KNOWN + NMORE
IF (LOOP1.GT.LOOP2 .OR. LOOP1.LT.IONE) GO TO 180
IF (GUESS-NMORE) 40, 70, 60
40 ILO = GUESS + IONE
DO 50 I=ILO,LOOP2
RTS(I) = ZERO
50 CONTINUE
GO TO 70
60 CONTINUE
GUESS = NMORE
70 CONTINUE
STEP = HALF
DO 160 NEW=LOOP1,LOOP2
KOUNT = ITHREE
NEWM1 = NEW - IONE
RT = RTS(NEW)
X1 = RT - STEP
X2 = RT + STEP
X3 = RT
CALL CALCF(X1, F1)
CALL CALCF(X2, F2)
CALL CALCF(X3, F3)
FSLAST = F3
IF (NEW.GT.IONE) CALL TEST(X1, F1, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X2, F2, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X3, F3, FSLAST, RTS,NEWM1, EP2, KOUNT)
F21 = (F2-F1)/(X2-X1)
80 X32 = X3 - X2
F32 = (F3-F2)/X32
F321 = (F32-F21)/(X3-X1)
B = F32 + X32*F321
RADICL = B**ITWO - FOUR*F321*F3
IF (REALRT .AND. REAL(RADICL).LT.ZERO) RADICL = ZERO
RADICL = SQRT(RADICL)
IF (REAL(B)*REAL(RADICL)+AIMAG(B)*AIMAG(RADICL).LT.ZERO) RADICL = -RADICL
DENOM = B + RADICL
IF (ABS(DENOM).NE.ZERO) GO TO 100
IF (ABS(F3).GE.EP2) GO TO 90
XNEW = X3
GO TO 120
90 XNEW = X3 + X32
GO TO 120
100 CSTEP = TWO*F3/DENOM
IF (.NOT.REALRT .OR. ABS(F3).EQ.ZERO .OR. ABS(F32).EQ.ZERO) GO TO 110
CSTEP = F32/ABS(F32)*F3/ABS(F3)*ABS(CSTEP)
110 XNEW = X3 - CSTEP
120 CALL CALCF(XNEW, FNEW)
FSAVE = FNEW
IF (NEW.LE.IONE) GO TO 130
CALL TEST(XNEW, FNEW, FSAVE, RTS, NEWM1, EP2,KOUNT)
130 KOUNT = KOUNT +IONE
IF (KOUNT.GT.MAXITS) GO TO 170
DIF = XNEW - X3
IF (ABS(DIF).LT.EP1*ABS(XNEW)) GO TO 150
IF (ABS(FSAVE).LT.EP2) GO TO 150
IF (REALRT .OR. ABS(FSAVE).LT.HUNDRD*ABS(FSLAST)) GO TO 140
CSTEP = CSTEP*HALF
XNEW = XNEM + CSTEP
GO TO 120
140 X1 = X2
X2 = X3
X3 = XNEW
F1 = F2
F2 = F3
F3 = FNEW
FSLAST = FSAVE
F21 = F32
GO TO 80
150 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
NFOUND = NEW
160 CONTINUE
GO TO 190
170 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
IFAIL = ITHREE
GO TO 190
180 CONTINUE
IFAIL = ITWO
190 CONTINUE
RETURN
END
SUBROUTINE TEST(X, F, FSAVE, RTS, K, EPS, KOUNT)
COMPLEX*16 X, F, RTS(K), D, FSAVE
REAL*8 EPS!, CABS
DATA IONE, PERTB /1,0.01d0/
10 CONTINUE
DO 20 I=1,K
D = X - RTS(I)
IF (ABS(D).LE.EPS) GO TO 30
F = F/D
20 CONTINUE
GO TO 40
30 CONTINUE
X = X + PERTB
CALL CALCF(X, F)
FSAVE = F
KOUNT = KOUNT +IONE
GO TO 10
40 RETURN
END
我已经修改了 ROOTS 子例程，并针对 gfortran 对其进行了调整（例如，REALPART 和 IMAGPART 用于复杂变量）。还更正了 KELLEKAI 发布的源代码中的一个错误，即 XNEW = XNEM + CSTEP，显然应该是 XNEW = XNEW + CSTEP（这就是为什么 IMPLICIT NONE 更可取）。原始评论也包括在内。我测试了具有多个复根的非线性方程的子程序，它工作正常。
C source code from I. Barrodale and K. B. Wilson,
C "A Fortran program for solving a nonlinear equation by Muller's method"
C Journal of Computational and Applied Mathematics, vol. 4, no. 2, 1978
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
SUBROUTINE ROOTS(KNOWN, NMORE, KM, REALRT, EP1, EP2, ROO 10
& GUESS, MAXITS, RTS, NFOUND, FNV, IFAIL) ROO 20
C-----------------------------------------------------------------------
C The user must supply a subprogram as follows
C SUBROUTINE CALCF (X, F)
C where X and F are complex numbers, X is an independent variable
C and F is the corresponding function value.
C-----------------------------------------------------------------------
C
C This subroutine uses Muller's method to find roots (real or complex)
C of a general equation
C
C Description of parameters:
C
C KNOWN an INTEGER denoting the number of roots already known.
C If the function value for one of these "known" roots
C is not "small" the value of known is reset to the last
C index of a root with a small function value, and the
C remaining "known" roots are recalculated.
C NMORE an INTEGER denoting the number of roots to be found
C (i.e. in addition to those already known).
C NMORE is reset to account for the case when KNOWN is reset.
C KM an INTEGER used to define dimensions.
C KM must be set at least (KNOWN + NMORE)
C REALRT a LOGICAL variable with values as follows:
C .TRUE. indicates that only REAL roots are to be found
C (i.e. complex numbers with zero imaginary parts)
C .FALSE. indicates that complex roots are to be found
C EP1 a small positive real number denoting the "first"
C convergence criterion: if at the K-th stage the relative
C change in two successive approximations to a root is less
C than EP1, convergence is said to have taken place
C (i.e. the absolute value of (ROOT(K)-ROOT(K-1))/ROOT(K)
C is less than EP1 the convergence occurs). See below for
C minimum default values
C EP2 a small positive real number denoting the "second"
C convergence criterion: if at anu stage the absolute value
C of the function is less than EP2, convergence is said to
C have taken place. See below for minimum default values.
C EP2 is also used by SUBROUTINE TEST to determine if a
C current approximation is "close" to a root already identified
C GUESS an INTEGER denoting the number of roots for which the user
C is supplying estimates. GUESS is reset to account for the
C case when KNOWN is reset. Also, GUESS must not exceed NMORE.
C MAXITS an INTEGER denoting the maximum number of function calls
C allowed when finding any root. If MAXITS function evaluations
C are performed the subroutine "fails" and return control to
C the user. The roots already found and the last approximation
C to the root causing "failure" are stored in RTS. See below
C for maximum default values.
C RTS a COMPLEX vector of length at least KM:
C On entry if any roots are available they must be in
C positions 1 to KNOWN of RTS.
C On entry, if estimates of the roots to be found are known,
C (i.e. GUESS.GT.ZERO) they must be in positions
C (KNOWN+1) to (KNOWN + GUESS).
C On return the roots of the equation are stored in the first
C NFOUND positions. If "failure" due to MAXITS return, the
C (NFOUND+1) position of RTS contains the last approximation
C to the root which caused the "failure".
C NFOUND an INTEGER denoting the number of roots now known and
C stored in RTS.
C FNV a COMPLEX vector of length at least KM containing the
C function values associated with the vector RTS.
C IFAIL an integer exit code with values as follows:
C 0 indicates all roots were found as requested by the user
C 1 indicates (as above) all roots were found but in addition
C some of the input "known" roots were recalculated
C 2 if the input parameters KNOWN or NMORE satisfy the
C following: KNOWN.LT.0, NMORE.LT.1 control is returned
C immediately to the user
C 3 if MAXITS causes "failure" control is returned
C immediately to the user
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
INTEGER KNOWN, NMORE, KM, GUESS, MAXITS, NFOUND, IFAIL ROO 30
COMPLEX*16 RTS(KM), FNV(KM), X1, X2, X3, F1, F2, F3, ROO 40
& X21, X32, F21, F32, F31, F321, B, XNEW, FNEW, DENOM, ROO 50
& FSAVE, CSQRT, RADICL, RT, DIF, CSTEP, FSLAST ROO 60
REAL*8 EP1, EP2, ONE, HUNDRD, TWO, FOUR, HALF, STEP, ROO 70
& ZERO, EP1DEF, EP2DEF ROO 80
LOGICAL REALRT ROO 90
DATA IZERO, IONE, ITWO, ITHREE, ZERO, HALF, ONE, TWO, ROO 100
& FOUR, HUNDRD, ITSDEF, EP1DEF, EP2DEF ROO 110
& /0,1,2,3,0.0d0,0.5d0,1.0d0,2.0d0,4.0d0,100.0d0,100, ROO 120
& 0.5D-6,1.0D-7/ ROO 130
C
C Initialization
C Check for errors in the input and set default values if necessary
C
EP1 = DMAX1(EP1,EP1DEF)
EP2 = DMAX1(EP2,EP2DEF)
MAXITS = MIN0(MAXITS,ITSDEF)
IFAIL = IZERO
NFOUND = IZERO
IF (KNOWN.LT.IONE) GO TO 30
DO 10 I=1,KNOWN
II = I
CALL CALCF(RTS(I), FNV(I))
C no record is kept of the number of function evaluations
C used here to ckeck the list of known roots
IF (ABS(FNV(II)).GE.EP2) GO TO 20
10 CONTINUE
20 IFAIL =IONE
GUESS = GUESS + KNOWN - II +IONE
NMORE = NMORE + KNOWN - II +IONE
KNOWN = II - IONE
30 CONTINUE
LOOP1 = KNOWN +IONE
LOOP2 = KNOWN + NMORE
IF (LOOP1.GT.LOOP2 .OR. LOOP1.LT.IONE) GO TO 180
IF (GUESS-NMORE) 40, 70, 60
40 ILO = GUESS + IONE
DO 50 I=ILO,LOOP2
RTS(I) = ZERO
50 CONTINUE
GO TO 70
60 CONTINUE
GUESS = NMORE
70 CONTINUE
STEP = HALF
C
C try to find NMORE roots in this major loop
DO 160 NEW=LOOP1,LOOP2
KOUNT = ITHREE
NEWM1 = NEW - IONE
RT = RTS(NEW)
X1 = RT - STEP
X2 = RT + STEP
X3 = RT
C use first estimates
CALL CALCF(X1, F1)
CALL CALCF(X2, F2)
CALL CALCF(X3, F3)
FSLAST = F3
IF (NEW.GT.IONE)
& CALL TEST(X1, F1, FSAVE, RTS, NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE)
& CALL TEST(X2, F2, FSAVE, RTS, NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE)
& CALL TEST(X3, F3, FSLAST, RTS, NEWM1, EP2, KOUNT)
F21 = (F2-F1)/(X2-X1)
C find the roots of the quadratic
80 X32 = X3 - X2
F32 = (F3-F2)/X32
F321 = (F32-F21)/(X3-X1)
B = F32 + X32*F321
RADICL = B**ITWO - FOUR*F321*F3
IF (REALRT .AND. REALPART(RADICL) .LT. ZERO) RADICL = ZERO
RADICL = CDSQRT(RADICL)
IF (
& (REALPART(B)*REALPART(RADICL)+IMAGPART(B)*IMAGPART(RADICL)
& ) .LT. ZERO) RADICL = -RADICL
DENOM = B + RADICL
IF (ABS(DENOM).NE.ZERO) GO TO 100
IF (ABS(F3).GE.EP2) GO TO 90
XNEW = X3
GO TO 120
90 XNEW = X3 + X32
GO TO 120
100 CSTEP = TWO*F3/DENOM
IF (.NOT. REALRT .OR. ABS(F3) .EQ. ZERO .OR.
& ABS(F32) .EQ. ZERO) GO TO 110
CSTEP = F32/ABS(F32)*F3/ABS(F3)*ABS(CSTEP)
110 XNEW = X3 - CSTEP
120 CALL CALCF(XNEW, FNEW)
FSAVE = FNEW
IF (NEW.LE.IONE) GO TO 130
CALL TEST(XNEW, FNEW, FSAVE, RTS, NEWM1, EP2,KOUNT)
130 KOUNT = KOUNT +IONE
C test iterations, convergence and divergence
IF (KOUNT.GT.MAXITS) GO TO 170
DIF = XNEW - X3
IF (ABS(DIF).LT.EP1*ABS(XNEW)) GO TO 150
IF (ABS(FSAVE).LT.EP2) GO TO 150
IF (REALRT .OR. ABS(FSAVE) .LT. HUNDRD*ABS(FSLAST))
& GO TO 140
C divergence "fix"
CSTEP = CSTEP*HALF
XNEW = XNEW + CSTEP
GO TO 120
140 X1 = X2
X2 = X3
X3 = XNEW
F1 = F2
F2 = F3
F3 = FNEW
FSLAST = FSAVE
F21 = F32
GO TO 80
150 CONTINUE
C convergence to one more root
C (the number of function evaluations used in determining the current root
C is available at this stage; output the variable KOUNT if desired)
write(*,*) ' number function evaluations = ', KOUNT
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
NFOUND = NEW
160 CONTINUE
GO TO 190
170 CONTINUE
C MAXITS reached
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
IFAIL = ITHREE
GO TO 190
180 CONTINUE
IFAIL = ITWO
190 CONTINUE
RETURN
END SUBROUTINE ROOTS
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
SUBROUTINE TEST(X, F, FSAVE, RTS, K, EPS, KOUNT) TES 10
C
C This subroutine serves two purposes:
C 1. Deflation of the function is performed,
C i.e. if X is the current approximation to a root
C and RTS is a vector of roots already known, the
C deflated function value is
C F=(...((F.(X-RTS(1)))/(X-RTS(2)))/...)
C 2. If a current approximation to a root is "too close"
C to a root already known, X is perturbed away from
C the known root by 0.01. This device avoids division
C by (essentially) zero in the deflation process, but it
C does not prevent convergence to multiple zeros.
C
COMPLEX*16 X, F, RTS(K), D, FSAVE TES 20
REAL*8 EPS, PERTB TES 30
DATA IONE, PERTB /1, 0.01D0/ TES 40
10 CONTINUE TES 50
DO 20 I=1,K TES 60
D = X - RTS(I) TES 70
IF (ABS(D) .LE. EPS) GO TO 30 TES 80
F = F/D TES 90
20 CONTINUE TES 100
GO TO 40 TES 110
30 CONTINUE TES 120
X = X + PERTB TES 130
CALL CALCF(X, F) TES 140
FSAVE = F TES 150
KOUNT = KOUNT + IONE TES 160
GO TO 10 TES 170
40 RETURN TES 180
END SUBROUTINE TEST TES 190
我发现了一些资源，其中建议使用算法来确定一般分析函数（无多项式）中的复根。例如。 “Muller 求解非线性方程的 Fortran 程序 方法”，由 I. Barrodale 和 K. B. Wilson 编写。不幸的是，该代码仅以 pdf 格式提供，并且转换为 ASCII 非常麻烦……还有 Henning Bach，第 121 卷，第 12 期，第 675-678 页，1969 年 12 月. 它在 ASCII 中可用。字母 1 工作正常但对我的目的来说相当慢。在 Fortran 中是否有计算复数根的标准例程？
谢谢 :)
对于那些有兴趣的人。我从 I. Barrodale 和 K. B. Wilson 那里拿了一张。由于这篇文章是免费在线阅读的，所以我想 post 这里的代码。这是迄今为止我发现的最快和最可靠的。要使用它，请定义一个子程序 calcf(x,f)。 f 是 x 的函数值。 x 和 f 是复数。
了解更多信息 -> http://www.sciencedirect.com/science/article/pii/0771050X78900414
这是双精度版本：
SUBROUTINE ROOTS(KNOWN,NMORE,KM,REALRT,EP1,EP2,GUESS,MAXITS,RTS,NFOUND, FNV, IFAIL)
COMPLEX*16 RTS(KM), FNV(KM), X1, X2, X3, F1, F2, F3, &
X21, X32, F21, F32, F31, F321, B, XNEW, FNEW, DENOM, &
FSAVE, CSQRT, RADICL, RT, DIF, CSTEP, FSLAST
REAL*8 EP1, EP2, ONE, HUNDRD, TWO, FOUR, HALF, STEP, &
ZERO, EP1DEF, EP2DEF
INTEGER GUESS
LOGICAL REALRT
DATA IZERO,IONE,ITWO,ITHREE,ZERO,HALF,ONE,TWO,FOUR,HUNDRD,ITSDEF,EP1DEF,EP2DEF &
/0,1,2,3,0.0d0,0.5d0,1.0d0,2.0d0,4.0d0,100.0d0,100,0.5d-5,1.0d-6/
EP1 = AMAX1(EP1,EP1DEF)
EP2 = AMAX1(EP2,EP2DEF)
MAXITS = MIN0(MAXITS,ITSDEF)
IFAIL = IZERO
NFOUND = IZERO
IF (KNOWN.LT.IONE) GO TO 30
DO 10 I=1,KNOWN
II = I
CALL CALCF(RTS(I), FNV(I))
IF (ABS(FNV(II)).GE.EP2) GO TO 20
10 CONTINUE
20 IFAIL =IONE
GUESS = GUESS + KNOWN - II +IONE
NMORE = NMORE + KNOWN - II +IONE
KNOWN = II - IONE
30 CONTINUE
LOOP1 = KNOWN +IONE
LOOP2 = KNOWN + NMORE
IF (LOOP1.GT.LOOP2 .OR. LOOP1.LT.IONE) GO TO 180
IF (GUESS-NMORE) 40, 70, 60
40 ILO = GUESS + IONE
DO 50 I=ILO,LOOP2
RTS(I) = ZERO
50 CONTINUE
GO TO 70
60 CONTINUE
GUESS = NMORE
70 CONTINUE
STEP = HALF
DO 160 NEW=LOOP1,LOOP2
KOUNT = ITHREE
NEWM1 = NEW - IONE
RT = RTS(NEW)
X1 = RT - STEP
X2 = RT + STEP
X3 = RT
CALL CALCF(X1, F1)
CALL CALCF(X2, F2)
CALL CALCF(X3, F3)
FSLAST = F3
IF (NEW.GT.IONE) CALL TEST(X1, F1, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X2, F2, FSAVE, RTS,NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE) CALL TEST(X3, F3, FSLAST, RTS,NEWM1, EP2, KOUNT)
F21 = (F2-F1)/(X2-X1)
80 X32 = X3 - X2
F32 = (F3-F2)/X32
F321 = (F32-F21)/(X3-X1)
B = F32 + X32*F321
RADICL = B**ITWO - FOUR*F321*F3
IF (REALRT .AND. REAL(RADICL).LT.ZERO) RADICL = ZERO
RADICL = SQRT(RADICL)
IF (REAL(B)*REAL(RADICL)+AIMAG(B)*AIMAG(RADICL).LT.ZERO) RADICL = -RADICL
DENOM = B + RADICL
IF (ABS(DENOM).NE.ZERO) GO TO 100
IF (ABS(F3).GE.EP2) GO TO 90
XNEW = X3
GO TO 120
90 XNEW = X3 + X32
GO TO 120
100 CSTEP = TWO*F3/DENOM
IF (.NOT.REALRT .OR. ABS(F3).EQ.ZERO .OR. ABS(F32).EQ.ZERO) GO TO 110
CSTEP = F32/ABS(F32)*F3/ABS(F3)*ABS(CSTEP)
110 XNEW = X3 - CSTEP
120 CALL CALCF(XNEW, FNEW)
FSAVE = FNEW
IF (NEW.LE.IONE) GO TO 130
CALL TEST(XNEW, FNEW, FSAVE, RTS, NEWM1, EP2,KOUNT)
130 KOUNT = KOUNT +IONE
IF (KOUNT.GT.MAXITS) GO TO 170
DIF = XNEW - X3
IF (ABS(DIF).LT.EP1*ABS(XNEW)) GO TO 150
IF (ABS(FSAVE).LT.EP2) GO TO 150
IF (REALRT .OR. ABS(FSAVE).LT.HUNDRD*ABS(FSLAST)) GO TO 140
CSTEP = CSTEP*HALF
XNEW = XNEM + CSTEP
GO TO 120
140 X1 = X2
X2 = X3
X3 = XNEW
F1 = F2
F2 = F3
F3 = FNEW
FSLAST = FSAVE
F21 = F32
GO TO 80
150 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
NFOUND = NEW
160 CONTINUE
GO TO 190
170 CONTINUE
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
IFAIL = ITHREE
GO TO 190
180 CONTINUE
IFAIL = ITWO
190 CONTINUE
RETURN
END
SUBROUTINE TEST(X, F, FSAVE, RTS, K, EPS, KOUNT)
COMPLEX*16 X, F, RTS(K), D, FSAVE
REAL*8 EPS!, CABS
DATA IONE, PERTB /1,0.01d0/
10 CONTINUE
DO 20 I=1,K
D = X - RTS(I)
IF (ABS(D).LE.EPS) GO TO 30
F = F/D
20 CONTINUE
GO TO 40
30 CONTINUE
X = X + PERTB
CALL CALCF(X, F)
FSAVE = F
KOUNT = KOUNT +IONE
GO TO 10
40 RETURN
END
我已经修改了 ROOTS 子例程，并针对 gfortran 对其进行了调整（例如，REALPART 和 IMAGPART 用于复杂变量）。还更正了 KELLEKAI 发布的源代码中的一个错误，即 XNEW = XNEM + CSTEP，显然应该是 XNEW = XNEW + CSTEP（这就是为什么 IMPLICIT NONE 更可取）。原始评论也包括在内。我测试了具有多个复根的非线性方程的子程序，它工作正常。
C source code from I. Barrodale and K. B. Wilson,
C "A Fortran program for solving a nonlinear equation by Muller's method"
C Journal of Computational and Applied Mathematics, vol. 4, no. 2, 1978
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
SUBROUTINE ROOTS(KNOWN, NMORE, KM, REALRT, EP1, EP2, ROO 10
& GUESS, MAXITS, RTS, NFOUND, FNV, IFAIL) ROO 20
C-----------------------------------------------------------------------
C The user must supply a subprogram as follows
C SUBROUTINE CALCF (X, F)
C where X and F are complex numbers, X is an independent variable
C and F is the corresponding function value.
C-----------------------------------------------------------------------
C
C This subroutine uses Muller's method to find roots (real or complex)
C of a general equation
C
C Description of parameters:
C
C KNOWN an INTEGER denoting the number of roots already known.
C If the function value for one of these "known" roots
C is not "small" the value of known is reset to the last
C index of a root with a small function value, and the
C remaining "known" roots are recalculated.
C NMORE an INTEGER denoting the number of roots to be found
C (i.e. in addition to those already known).
C NMORE is reset to account for the case when KNOWN is reset.
C KM an INTEGER used to define dimensions.
C KM must be set at least (KNOWN + NMORE)
C REALRT a LOGICAL variable with values as follows:
C .TRUE. indicates that only REAL roots are to be found
C (i.e. complex numbers with zero imaginary parts)
C .FALSE. indicates that complex roots are to be found
C EP1 a small positive real number denoting the "first"
C convergence criterion: if at the K-th stage the relative
C change in two successive approximations to a root is less
C than EP1, convergence is said to have taken place
C (i.e. the absolute value of (ROOT(K)-ROOT(K-1))/ROOT(K)
C is less than EP1 the convergence occurs). See below for
C minimum default values
C EP2 a small positive real number denoting the "second"
C convergence criterion: if at anu stage the absolute value
C of the function is less than EP2, convergence is said to
C have taken place. See below for minimum default values.
C EP2 is also used by SUBROUTINE TEST to determine if a
C current approximation is "close" to a root already identified
C GUESS an INTEGER denoting the number of roots for which the user
C is supplying estimates. GUESS is reset to account for the
C case when KNOWN is reset. Also, GUESS must not exceed NMORE.
C MAXITS an INTEGER denoting the maximum number of function calls
C allowed when finding any root. If MAXITS function evaluations
C are performed the subroutine "fails" and return control to
C the user. The roots already found and the last approximation
C to the root causing "failure" are stored in RTS. See below
C for maximum default values.
C RTS a COMPLEX vector of length at least KM:
C On entry if any roots are available they must be in
C positions 1 to KNOWN of RTS.
C On entry, if estimates of the roots to be found are known,
C (i.e. GUESS.GT.ZERO) they must be in positions
C (KNOWN+1) to (KNOWN + GUESS).
C On return the roots of the equation are stored in the first
C NFOUND positions. If "failure" due to MAXITS return, the
C (NFOUND+1) position of RTS contains the last approximation
C to the root which caused the "failure".
C NFOUND an INTEGER denoting the number of roots now known and
C stored in RTS.
C FNV a COMPLEX vector of length at least KM containing the
C function values associated with the vector RTS.
C IFAIL an integer exit code with values as follows:
C 0 indicates all roots were found as requested by the user
C 1 indicates (as above) all roots were found but in addition
C some of the input "known" roots were recalculated
C 2 if the input parameters KNOWN or NMORE satisfy the
C following: KNOWN.LT.0, NMORE.LT.1 control is returned
C immediately to the user
C 3 if MAXITS causes "failure" control is returned
C immediately to the user
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
INTEGER KNOWN, NMORE, KM, GUESS, MAXITS, NFOUND, IFAIL ROO 30
COMPLEX*16 RTS(KM), FNV(KM), X1, X2, X3, F1, F2, F3, ROO 40
& X21, X32, F21, F32, F31, F321, B, XNEW, FNEW, DENOM, ROO 50
& FSAVE, CSQRT, RADICL, RT, DIF, CSTEP, FSLAST ROO 60
REAL*8 EP1, EP2, ONE, HUNDRD, TWO, FOUR, HALF, STEP, ROO 70
& ZERO, EP1DEF, EP2DEF ROO 80
LOGICAL REALRT ROO 90
DATA IZERO, IONE, ITWO, ITHREE, ZERO, HALF, ONE, TWO, ROO 100
& FOUR, HUNDRD, ITSDEF, EP1DEF, EP2DEF ROO 110
& /0,1,2,3,0.0d0,0.5d0,1.0d0,2.0d0,4.0d0,100.0d0,100, ROO 120
& 0.5D-6,1.0D-7/ ROO 130
C
C Initialization
C Check for errors in the input and set default values if necessary
C
EP1 = DMAX1(EP1,EP1DEF)
EP2 = DMAX1(EP2,EP2DEF)
MAXITS = MIN0(MAXITS,ITSDEF)
IFAIL = IZERO
NFOUND = IZERO
IF (KNOWN.LT.IONE) GO TO 30
DO 10 I=1,KNOWN
II = I
CALL CALCF(RTS(I), FNV(I))
C no record is kept of the number of function evaluations
C used here to ckeck the list of known roots
IF (ABS(FNV(II)).GE.EP2) GO TO 20
10 CONTINUE
20 IFAIL =IONE
GUESS = GUESS + KNOWN - II +IONE
NMORE = NMORE + KNOWN - II +IONE
KNOWN = II - IONE
30 CONTINUE
LOOP1 = KNOWN +IONE
LOOP2 = KNOWN + NMORE
IF (LOOP1.GT.LOOP2 .OR. LOOP1.LT.IONE) GO TO 180
IF (GUESS-NMORE) 40, 70, 60
40 ILO = GUESS + IONE
DO 50 I=ILO,LOOP2
RTS(I) = ZERO
50 CONTINUE
GO TO 70
60 CONTINUE
GUESS = NMORE
70 CONTINUE
STEP = HALF
C
C try to find NMORE roots in this major loop
DO 160 NEW=LOOP1,LOOP2
KOUNT = ITHREE
NEWM1 = NEW - IONE
RT = RTS(NEW)
X1 = RT - STEP
X2 = RT + STEP
X3 = RT
C use first estimates
CALL CALCF(X1, F1)
CALL CALCF(X2, F2)
CALL CALCF(X3, F3)
FSLAST = F3
IF (NEW.GT.IONE)
& CALL TEST(X1, F1, FSAVE, RTS, NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE)
& CALL TEST(X2, F2, FSAVE, RTS, NEWM1, EP2, KOUNT)
IF (NEW.GT.IONE)
& CALL TEST(X3, F3, FSLAST, RTS, NEWM1, EP2, KOUNT)
F21 = (F2-F1)/(X2-X1)
C find the roots of the quadratic
80 X32 = X3 - X2
F32 = (F3-F2)/X32
F321 = (F32-F21)/(X3-X1)
B = F32 + X32*F321
RADICL = B**ITWO - FOUR*F321*F3
IF (REALRT .AND. REALPART(RADICL) .LT. ZERO) RADICL = ZERO
RADICL = CDSQRT(RADICL)
IF (
& (REALPART(B)*REALPART(RADICL)+IMAGPART(B)*IMAGPART(RADICL)
& ) .LT. ZERO) RADICL = -RADICL
DENOM = B + RADICL
IF (ABS(DENOM).NE.ZERO) GO TO 100
IF (ABS(F3).GE.EP2) GO TO 90
XNEW = X3
GO TO 120
90 XNEW = X3 + X32
GO TO 120
100 CSTEP = TWO*F3/DENOM
IF (.NOT. REALRT .OR. ABS(F3) .EQ. ZERO .OR.
& ABS(F32) .EQ. ZERO) GO TO 110
CSTEP = F32/ABS(F32)*F3/ABS(F3)*ABS(CSTEP)
110 XNEW = X3 - CSTEP
120 CALL CALCF(XNEW, FNEW)
FSAVE = FNEW
IF (NEW.LE.IONE) GO TO 130
CALL TEST(XNEW, FNEW, FSAVE, RTS, NEWM1, EP2,KOUNT)
130 KOUNT = KOUNT +IONE
C test iterations, convergence and divergence
IF (KOUNT.GT.MAXITS) GO TO 170
DIF = XNEW - X3
IF (ABS(DIF).LT.EP1*ABS(XNEW)) GO TO 150
IF (ABS(FSAVE).LT.EP2) GO TO 150
IF (REALRT .OR. ABS(FSAVE) .LT. HUNDRD*ABS(FSLAST))
& GO TO 140
C divergence "fix"
CSTEP = CSTEP*HALF
XNEW = XNEW + CSTEP
GO TO 120
140 X1 = X2
X2 = X3
X3 = XNEW
F1 = F2
F2 = F3
F3 = FNEW
FSLAST = FSAVE
F21 = F32
GO TO 80
150 CONTINUE
C convergence to one more root
C (the number of function evaluations used in determining the current root
C is available at this stage; output the variable KOUNT if desired)
write(*,*) ' number function evaluations = ', KOUNT
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
NFOUND = NEW
160 CONTINUE
GO TO 190
170 CONTINUE
C MAXITS reached
RTS(NEW) = XNEW
FNV(NEW) = FSAVE
IFAIL = ITHREE
GO TO 190
180 CONTINUE
IFAIL = ITWO
190 CONTINUE
RETURN
END SUBROUTINE ROOTS
cccccCcccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxcccccccccxccCCCCCCCC
SUBROUTINE TEST(X, F, FSAVE, RTS, K, EPS, KOUNT) TES 10
C
C This subroutine serves two purposes:
C 1. Deflation of the function is performed,
C i.e. if X is the current approximation to a root
C and RTS is a vector of roots already known, the
C deflated function value is
C F=(...((F.(X-RTS(1)))/(X-RTS(2)))/...)
C 2. If a current approximation to a root is "too close"
C to a root already known, X is perturbed away from
C the known root by 0.01. This device avoids division
C by (essentially) zero in the deflation process, but it
C does not prevent convergence to multiple zeros.
C
COMPLEX*16 X, F, RTS(K), D, FSAVE TES 20
REAL*8 EPS, PERTB TES 30
DATA IONE, PERTB /1, 0.01D0/ TES 40
10 CONTINUE TES 50
DO 20 I=1,K TES 60
D = X - RTS(I) TES 70
IF (ABS(D) .LE. EPS) GO TO 30 TES 80
F = F/D TES 90
20 CONTINUE TES 100
GO TO 40 TES 110
30 CONTINUE TES 120
X = X + PERTB TES 130
CALL CALCF(X, F) TES 140
FSAVE = F TES 150
KOUNT = KOUNT + IONE TES 160
GO TO 10 TES 170
40 RETURN TES 180
END SUBROUTINE TEST TES 190