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
|
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!
c Copyright INRIA
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
|
