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
|
*DECK ALROOT
SUBROUTINE ALROOT(A,M,ROOTS,LFAIL)
*
*-----------------------------------------------------------------------
*
*Purpose:
* find the roots of a polynomial.
*
*Copyright:
* Copyright (C) 1993 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): A. Hebert
*
*Parameters: input
* A polynomial coefficients.
* M polynomial order.
*
*Parameters: output
* ROOTS complex roots.
* LFAIL flag set to .true. in case of failure.
*
*-----------------------------------------------------------------------
*
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER M
DOUBLE PRECISION A(M+1)
COMPLEX*16 ROOTS(M)
LOGICAL LFAIL
*----
* LOCAL VARIABLES
*----
COMPLEX*16 AAA,BBB,SQ1,TEST,SQRTM3
PARAMETER (SQRTM3=(0.0,1.73205080756888))
DOUBLE PRECISION EPS
PARAMETER (EPS=1.D-6,MAXM=101)
COMPLEX*16 AD(MAXM),X,B,C
*
LFAIL=.FALSE.
IF(M+1.GT.MAXM) CALL XABORT('ALROOT: INSUFFICIENT STORAGE.')
IF(A(M+1).EQ.0.0D0) CALL XABORT('ALROOT: INVALID COEFFICIENT.')
IF(M.EQ.1) THEN
ROOTS(1)=-A(1)/A(2)
ELSE IF(M.EQ.2) THEN
CQ=A(2)/A(3)
DQ=A(1)/A(3)
AAA=CQ*CQ-4.0D0*DQ
AAA=SQRT(AAA)
ROOTS(1)=-0.5D0*(CQ+AAA)
ROOTS(2)=-0.5D0*(CQ-AAA)
ELSE IF(M.EQ.3) THEN
IF(A(1).EQ.0.0) THEN
CQ=A(3)/A(4)
DQ=A(2)/A(4)
AAA=CQ*CQ-4.0D0*DQ
AAA=SQRT(AAA)
ROOTS(1)=0.0
ROOTS(2)=-0.5D0*(CQ+AAA)
ROOTS(3)=-0.5D0*(CQ-AAA)
ELSE
BQ=A(3)/A(4)
CQ=A(2)/A(4)
DQ=A(1)/A(4)
AA=(3.0D0*CQ-BQ**2)/3.0D0
BB=(2.0D0*BQ**3-9.0D0*BQ*CQ+27.0D0*DQ)/27.0D0
SQ1=BB**2/4.0D0+AA**3/27.0D0
TEST=BB/2.0D0-SQRT(SQ1)
IF(DBLE(TEST).EQ.0.0) THEN
AAA=0.0D0
ELSE IF(DBLE(TEST).GT.0.0) THEN
AAA=-(TEST)**(1.0D0/3.0D0)
ELSE
AAA=(-TEST)**(1.0D0/3.0D0)
ENDIF
TEST=BB/2.0D0+SQRT(SQ1)
IF(DBLE(TEST).EQ.0.0) THEN
BBB=0.0D0
ELSE IF(DBLE(TEST).GT.0.0) THEN
BBB=-(TEST)**(1.0D0/3.0D0)
ELSE
BBB=(-TEST)**(1.0D0/3.0D0)
ENDIF
ROOTS(1)=AAA+BBB-BQ/3.0D0
ROOTS(2)=-(AAA+BBB)/2.0D0+(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
ROOTS(3)=-(AAA+BBB)/2.0D0-(AAA-BBB)*SQRTM3/2.0D0-BQ/3.0D0
ENDIF
ELSE IF(M.EQ.4) THEN
CALL ALQUAR(A,ROOTS)
ELSE
DO 10 J=1,M+1
AD(J)=CMPLX(A(J),0.D0,KIND=KIND(AD))
10 CONTINUE
DO 25 J=M,1,-1
X=CMPLX(0.D0,0.D0,KIND=KIND(X))
CALL ALGUER(AD,J,X,ITS,LFAIL)
IF(LFAIL) RETURN
IF(ABS(DIMAG(X)).LE.2.D0*EPS**2*ABS(DBLE(X)))
1 X=CMPLX(DBLE(X),0.D0,KIND=KIND(X))
ROOTS(J)=X
B=AD(J+1)
DO 20 JJ=J,1,-1
C=AD(JJ)
AD(JJ)=B
B=X*B+C
20 CONTINUE
25 CONTINUE
DO 30 J=1,M+1
AD(J)=CMPLX(A(J),0.D0,KIND=KIND(AD))
30 CONTINUE
DO 40 J=1,M
CALL ALGUER(AD,M,ROOTS(J),ITS,LFAIL)
IF(LFAIL) RETURN
40 CONTINUE
ENDIF
*
DO 70 J=2,M
X=ROOTS(J)
DO 50 I=J-1,1,-1
IF(DBLE(ROOTS(I)).LE.DBLE(X)) GOTO 60
ROOTS(I+1)=ROOTS(I)
50 CONTINUE
I=0
60 ROOTS(I+1)=X
70 CONTINUE
RETURN
END
|
