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
|
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
c Copyright (C) 1986-2008 - INRIA - Francois DELBECQUE
c Copyright (C) 1990-2008 - INRIA - Serge STEER
c
c This file must be used under the terms of the CeCILL.
c This source file is licensed as described in the file COPYING, which
c you should have received as part of this distribution. The terms
c are also available at
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
subroutine sfact1(b,n,w,maxit,ierr)
c!but
c on cherche une factorisation spectrale d'un polynome a donne
c par : a*(a(1/z) = b(n)*z**-n+....+b(0)+ ... +b(n)*z**n
c!liste d'appel
c subroutine sfact1(b,n,w,maxit,ierr)
c
c double precision b(n+1),w(6*(n+1))
c integer n,maxit,ierr
c
c b : contient les coeffs b(n),b(n-1),....,b(0)
c apres execution b contient les coeff du resultat
c n : degre de a
c w : tableau de travail de taille 6*(n+1)
c maxit : nombre maxi d'iterations admis
c ierr : indicateur d'erreur
c si ierr=0 : ok
c si ierr<0 : convergence a 10**ierr pres
c si ierr=1 : non convergence
c si ierr=2 : b(0) est negatif ou nul.
c
c!origine
c V Kucera : Discrete Linear control (john Wiley& Sons) 1979
c!
dimension b(n+1),w(*)
double precision b,b0,w,a0,temp,b00,s,dlamch,eps,best
c
eps=dlamch('p')*10.0d+0
c
lb=n+1
ierr=0
c
lomeg=1
lalpha=lomeg+lb
lro=lalpha+lb
leta=lro+lb
lbold=leta+lb
lambda=lbold+lb
lsave=lambda+lb
c
call dcopy(lb,b,-1,w(lbold),1)
call dcopy(lb,w(lbold),1,b,1)
b00=w(lbold)
if(b00.le.0.0d+0) goto 91
b0=sqrt(b00)
do 1 j=1,lb
w(lalpha-1+j)=b(j)/b0
1 continue
c
do 40 i=1,maxit
c
call dcopy(lb,w(lbold),1,b,1)
call dcopy(lb,w(lalpha),1,w(lomeg),1)
c premier tableau
do 15 k=1,lb-1
call dcopy(lb-k+1,w(lalpha),-1,w(lro),1)
w(lambda+k-1)=w(lalpha+lb-k)/w(lro+lb-k)
do 11 j=1,lb-k
w(lalpha-1+j)=w(lalpha-1+j)-w(lambda+k-1)*w(lro+j-1)
11 continue
a0=w(lalpha)
c calcul de eta
w(leta+lb-k)=2.0d+0*b(lb-k+1)/a0
if (k.lt.lb-1) then
do 12 j=2,lb-k
b(j)=b(j)-0.50d+0*w(leta+lb-k)*w(lalpha+lb-k-j+1)
12 continue
endif
15 continue
w(leta)=b(1)/w(lalpha)
c deuxieme tableau
do 21 k=lb-1,1,-1
call dcopy(lb-k+1,w(leta),-1,b,1)
do 19 j=1,lb-k+1
w(leta+j-1)=w(leta+j-1)-w(lambda+k-1)*b(j)
19 continue
21 continue
s=0.0d+0
do 22 j=1,lb
w(lalpha-1+j)=0.50d+0*(w(leta+j-1)+w(lomeg+j-1))
s=s+w(lalpha-1+j)*w(lalpha-1+j)
22 continue
c
c test de convergence
c calcul de l'erreur element par elements
c$$$ call dcopy(lb,w(lbold),-1,b,1)
c$$$ do 900 iii=0,n
c$$$ x=b(iii+1)-ddot(iii+1,w(lalpha),1,w(lalpha+n-iii),1)
c$$$ write(6,'(10x,e10.3,'','',e10.3,'';'')') x,b(1+iii)
c$$$ 900 continue
temp=abs(s-b00)/b00
if(temp.le.eps) goto 50
if(i.eq.1) best=temp
if(temp.lt.best) then
call dcopy(lb,w(lalpha),1,w(lsave),1)
best=temp
endif
40 continue
goto 90
c
50 do 51 i=1,lb
b(i)=w(lalpha-1+i)
51 continue
return
c
c retours en erreur
c
90 continue
if(best.le.1.d-3) then
call dcopy(lb,w(lsave),1,b,1)
ierr=nint(log10(best))
else
c non convergence
ierr=1
endif
return
91 continue
c b00 est negatif ou nul
ierr=2
return
end
|
