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
|
subroutine feq(neq,t,tq,tqdot)
c!but
c Etablir la valeur de l'oppose du gradient au point q
c!liste d'appel
c subroutine feq(neq,t,tq,tqdot)
c - neq. tableau entier de taille 3+(nq+1)*(nq+2)
c neq(1)=nq est le degre effectif du polynome tq (ou q).
c neq(2)=ng est le nombre de coefficient de fourier
c neq(3)=dgmax degre maximum pour q (l'adresse des coeff de fourier dans
c tq est neq(3)+2
c - t . variable parametrique necessaire a l'execution de
c la routine lsoda .
c - tq. tableau reel de taille au moins
c 3+dgmax+nq+2*ng
c tq(1:nq+1) est le tableau des coefficients du polynome q.
c tq(dgmax+2:dgmax+ng+2) est le tableau des coefficients
c de fourier
c tq(dgmax+ng+3:) est un tableau de travail de taille au moins
c nq+ng+1
c Sortie :
c - tqdot . tableau contenant les opposes des coordonnees du
c gradient de la fonction PHI au point q
c!Remarque
c la structure particuliere pour neq et tq est liee au fait que feq peut
c etre appele comme un external de lsode
c!
implicit double precision (a-h,o-y)
dimension tq(*),tqdot(*)
dimension neq(*)
c
nq=neq(1)
ng=neq(2)
c
c decoupage du tableau tq
itq=1
itg=itq+neq(3)+1
iw=itg+ng+1
ifree=iw+1+nq+ng
call feq1(nq,t,tq,tq(itg),ng,tqdot,tq(iw))
return
end
subroutine feqn(neq,t,tq,tqdot)
c!but
c Etablir la valeur du gradient au point q
c!liste d'appel
c subroutine feqn(neq,t,tq,tqdot)
c - neq. tableau entier de taille 3+(nq+1)*(nq+2)
c neq(1)=nq est le degre effectif du polynome tq (ou q).
c neq(2)=ng est le nombre de coefficient de fourier
c neq(3)=dgmax degre maximum pour q (l'adresse des coeff de fourier dans
c tq est neq(3)+2
c - t . variable parametrique necessaire a l'execution de
c la routine lsoda .
c - tq. tableau reel de taille au moins
c 3+dgmax+nq+2*ng
c tq(1:nq+1) est le tableau des coefficients du polynome q.
c tq(dgmax+2:dgmax+ng+2) est le tableau des coefficients
c de fourier
c tq(dgmax+ng+3:) est un tableau de travail de taille au moins
c nq+ng+1
c Sortie :
c - tqdot . tableau contenant les opposes des coordonnees du
c gradient de la fonction PHI au point q
c!Remarque
c la structure particuliere pour neq et tq est liee au fait que feq peut
c etre appele comme un external de lsode
c!
implicit double precision (a-h,o-y)
dimension tq(*),tqdot(*)
dimension neq(*)
c
nq=neq(1)
ng=neq(2)
c
c decoupage du tableau tq
itq=1
itg=itq+neq(3)+1
iw=itg+ng+1
ifree=iw+1+nq+ng
call feq1(nq,t,tq,tq(itg),ng,tqdot,tq(iw))
do 10 i=1,nq
tqdot(i)=-tqdot(i)
10 continue
return
end
subroutine feq1(nq,t,tq,tg,ng,tqdot,tr)
implicit double precision (a-h,o-y)
dimension tq(nq+1),tqdot(nq),tg(*)
dimension tr(nq+ng+1)
c
do 199 i=1,nq
c
c -- calcul du terme general --
c
if (i.eq.1) then
call lq(nq,tq,tr,tg,ng)
c . tlq =tr(1:nq); tvq =tr(nq+1:nq+ng+1)
ltlq=1
ltvq=nq+1
c
c division de tvq par q
call dpodiv(tr(ltvq),tq,ng,nq)
nv=ng-nq
else
ichoix=1
call mzdivq(ichoix,nv,tr(ltvq),nq,tq)
endif
c
c calcul de tvq~ sur place
nr=nq-1
call tild(nr,tr(ltvq),tr)
call calsca(nq,tq,tr,y0,tg,ng)
c
c -- conclusion --
c
tqdot(i)=-2.0d+0*y0
c
199 continue
c write(6,'(''tqdot='',5(e10.3,2x))') (tqdot(i),i=1,nq)
c
return
end
|
