← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to routines/control/optml2.f

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
routines/control/optml2.f
 
1
      subroutine optml2(feq,jacl2,neq,q,nch,w,iw)
 
2
C!but
 
3
C      Routine de recherche de minimum du probleme d'approximation L2
 
4
C       par lsoda ( Lsoda = routine de resolution d'equa diff )
 
5
C!liste d'appel
 
6
C     subroutine optml2(feq,jacl2,neq,q,nch,w,iw)
 
7
C
 
8
C     external feq,jacl2
 
9
C     double precision q(*),w(*)
 
10
C     integer nch,iw(*)
 
11
C
 
12
C     Entrees :
 
13
C     - feq est la subroutine qui calcule le gradient,
 
14
C        oppose de la derivee premiere de la fonction phi.
 
15
c     - neq. tableau entier de taille 3+(npara+1)*(npara+2)
 
16
c         neq(1)=nq est le degre effectif du polynome  q.
 
17
c         neq(2)=ng est le nombre de coefficient de fourier
 
18
c         neq(3)=dgmax degre maximum pour q (l'adresse des coeff de fourier dans
 
19
c               q est neq(3)+2
 
20
C     - neq est le degre du polynome q
 
21
c     - tq. tableau reel de taille au moins
 
22
c               6+dgmax+5*nq+6*ng+nq*ng+nq**2*(ng+1)
 
23
c         tq(1:nq+1) est le tableau des coefficients du polynome q.
 
24
c         tq(dgmax+2:dgmax+ng+2) est le tableau des coefficients
 
25
c                      de fourier
 
26
c         tq(dgmax+ng+3:) est un tableau de travail de taille au moins
 
27
c                         5+5*nq+5*ng+nq*ng+nq**2*(ng+1)
 
28
C     - nch est l indice (valant 1 ou 2) qui classifie l
 
29
C       appel comme etant soit celui de la recherche et de la
 
30
C       localisation d un minimum local, soit de la
 
31
C       confirmation d un minimum local.
 
32
C
 
33
C     Sorties :
 
34
C     - neq est toujours le degre du polynome q (il peut  avoir varie).
 
35
C     - q est le polynome (ou plutot le tableau contenant
 
36
C         ses coefficients) qui resulte de la recherche ,il peut
 
37
C         etre du meme degre que le polynome initial mais aussi
 
38
C         de degre inferieur dans le cas d'une sortie de face.
 
39
C
 
40
C     Tableau de travail
 
41
C     - w de taille 25+26*nq+ng+nq**2
 
42
C     - iw de taille 20+nq
 
43
C!
 
44
c     Copyright INRIA
 
45
      implicit double precision (a-h,o-y)
 
46
      dimension q(*), w(*), iw(*), xx(1)
 
47
      integer neq(*)
 
48
      double precision x,phi0,phi,gnrm
 
49
C
 
50
      external feq, jacl2
 
51
      common /temps/ t
 
52
      common /comall/ nall1 /sortie/ io,info,ll
 
53
      common /no2f/ gnrm
 
54
c
 
55
      nq=neq(1)
 
56
      ng=neq(2)
 
57
      ltg=1+neq(3)
 
58
C
 
59
c     taille des tableaux de travail necessaires a lsode
 
60
      lrw = nq**2 + 9*nq + 22  
 
61
      liw = 20+nq
 
62
 
 
63
C     decoupage du tableau de travail w
 
64
      lqi = 1
 
65
      lqdot = lqi + nq + 1
 
66
      latol = lqdot + nq
 
67
      lrtol = latol + nq
 
68
      lwork = lrtol + nq
 
69
      lfree = lwork + 24+22*nq+ng+nq**2
 
70
c
 
71
      lw = lwork + lrw
 
72
 
 
73
C     decoupage du tableau de travail iw 
 
74
      liww=1
 
75
      lifree=liww+liw
 
76
C
 
77
      nqbac = nq
 
78
C
 
79
C     --- Initialisation de lsode ------------------------
 
80
C
 
81
      if (nch .eq. 1) t = 0.0d+0
 
82
      t0 = t
 
83
      tt = 0.10d+0
 
84
      tout = t0 + tt
 
85
      itol = 4
 
86
C
 
87
      if (nq .lt. 7) then
 
88
        ntol = int((nq-1)/3) + 5
 
89
      else
 
90
        ntol = int((nq-7)/2) + 7
 
91
      endif
 
92
      call dset(nq,10.0d+0**(-(ntol)),w(lrtol),1)
 
93
      call dset(nq,10.0d+0**(-(ntol+2)),w(latol),1)
 
94
C
 
95
      itask = 1
 
96
      if (nch .eq. 1) istate = 1
 
97
      if (nch .eq. 2) istate = 3
 
98
      iopt = 0
 
99
      mf = 21
 
100
C
 
101
C     --- Initialisation du nombre maximal d'iteration ---
 
102
C
 
103
      if (nch .eq. 1) then
 
104
        if (nq .le. 11) then
 
105
          nlsode = 11 + 2*(nq-1)
 
106
        else
 
107
          nlsode = 29
 
108
        endif
 
109
      else
 
110
        nlsode = 19
 
111
      endif
 
112
      ilcom = 0
 
113
      ipass = 0
 
114
C
 
115
C     --- Appel  de lsode --------------------------------
 
116
C
 
117
 210  do 290 i = 1,nlsode
 
118
C
 
119
 220    ilcom = ilcom + 1
 
120
C
 
121
C     -- Reinitialisation de la Tolerance --
 
122
C
 
123
        if (ilcom.eq.2 .and. nch.eq.1) then
 
124
          call dset(nq,1.0d-05,w(lrtol),1)
 
125
          call dset(nq,1.0d-07,w(latol),1)
 
126
          istate = 3
 
127
        elseif (ilcom.eq.2 .and. nch.eq.2) then
 
128
          w(lrtol) = 1.0d-08
 
129
          w(latol) = 1.0d-10
 
130
          w(lrtol+1) = 1.0d-07
 
131
          w(latol+1) = 1.0d-09
 
132
          w(lrtol+nq-1) = 1.0d-05
 
133
          w(latol+nq-1) = 1.0d-07
 
134
          do 240 j = 2,nq-2
 
135
            w(lrtol+j) = 1.0d-06
 
136
            w(latol+j) = 1.0d-08
 
137
 240      continue
 
138
          istate = 3
 
139
        endif
 
140
C
 
141
C     --------------------------------------
 
142
C
 
143
        call dcopy(nq+1,q,1,w(lqi),1)
 
144
        ti = t
 
145
        touti = tout
 
146
C
 
147
        if (info .gt. 1) call outl2(30,nq,nq,q,xx,t,tout)
 
148
C
 
149
 
 
150
        call lsode(feq,neq,q,t,tout,itol,w(lrtol),w(latol),itask,
 
151
     &             istate,iopt,w(lwork),lrw,iw(liww),liw,jacl2,mf)
 
152
C
 
153
        call front(nq,q,nbout,w(lw))
 
154
C
 
155
        call feq(neq,t,q,w(lqdot))
 
156
        dnorm0 = dnrm2(nq,w(lqdot),1)
 
157
        if (info .gt. 1) call outl2(31,nq,nbout,q,dnorm0,t,tout)
 
158
C
 
159
C     -- test pour degre1 -----------
 
160
        if (nall1.gt.0 .and. nq.eq.1 .and. nbout.gt.0) return
 
161
C
 
162
C
 
163
C     -- Istate de lsode ------------
 
164
C
 
165
        if (istate .eq. -5) then
 
166
          if (info .gt. 0) call outl2(32,nq,nq,xx,xx,x,x)
 
167
          call dscal(nq,0.10d+0,w(lrtol),1)
 
168
          call dscal(nq,0.10d+0,w(latol),1)
 
169
          if (t .eq. 0.0d+0) istate = 1
 
170
          if (t .ne. 0.0d+0) istate = 3
 
171
          ilcom = 0
 
172
          goto 220
 
173
        endif
 
174
C
 
175
        if (istate .eq. -6) then
 
176
C     echec de l'integration appel avec de nouvelles tolerances
 
177
          if (info .gt. 0) call outl2(33,nq,nq,xx,xx,x,x)
 
178
          if (info .gt. 1)
 
179
     &      call outl2(34,nq,itol,w(latol),w(lrtol),t,tout)
 
180
          iopt = 0
 
181
          itol = 4
 
182
          call dset(nq,0.10d-05,w(lrtol),1)
 
183
          call dset(nq,0.10d-05,w(latol),1)
 
184
          if (info .gt. 1) call outl2(35,nq,itol,w(latol),w(lrtol),x,x)
 
185
          if (info .gt. 0) call outl2(36,nq,nq,xx,xx,x,x)
 
186
          istate = 3
 
187
          if (t .ne. tout) goto 220
 
188
        endif
 
189
C
 
190
        if (istate.lt.-1 .and. istate.ne.-6 .and. istate.ne.-5) then
 
191
          if (info .gt. 0) call outl2(37,nq,iopt,xx,xx,x,x)
 
192
          nch = 15
 
193
          return
 
194
        endif
 
195
C
 
196
C     -------------------------------
 
197
C
 
198
C     -- Sortie de face -------------
 
199
C
 
200
        if (nbout.gt.0 .and. nbout.ne.99) then
 
201
          call domout(neq,q,w(lqi),nbout,ti,t,itol,w(lrtol),
 
202
     &          w(latol),itask,istate,iopt,w(lwork),lrw,iw(liww),liw,
 
203
     &          jacl2,mf,job)
 
204
          nq=neq(1)
 
205
          if (job .eq. -1) then
 
206
C     anomalie dans la recherche de l'intersection
 
207
            nch = 16
 
208
            return
 
209
          endif
 
210
          if (job .eq. 1) then
 
211
            nch = nq - nqbac
 
212
            return
 
213
          endif
 
214
        endif
 
215
C
 
216
C     -------------------------------
 
217
C
 
218
C     -- test sur le gradient -------
 
219
C
 
220
        epstop = (1.0d-06)**nch
 
221
        call feq(neq,t,q,w(lqdot))
 
222
        dnorm0 = dnrm2(nq,w(lqdot),1)
 
223
        if (dnorm0 .lt. epstop) goto 299
 
224
C
 
225
C     -------------------------------
 
226
C
 
227
C     -- Istate de lsode (suite) ----
 
228
C
 
229
        if (istate.eq.-1 .and. t.ne.tout) then
 
230
          if (info .gt. 0) call outl2(38,nq,nq,xx,xx,x,x)
 
231
          istate = 2
 
232
          goto 220
 
233
        endif
 
234
C
 
235
C     -------------------------------
 
236
C
 
237
        tt = sqrt(10.0d+0) * tt
 
238
        tout = t0 + tt
 
239
C
 
240
 290  continue
 
241
C
 
242
      if (nch.eq.2 .and. dnorm0.gt.(1.0d-06)) then
 
243
        ipass = ipass + 1
 
244
        if (ipass .lt. 5) then
 
245
          if (info .gt. 0) then
 
246
             call lq(nq,q,w(lw),q(ltg),ng)
 
247
             x=sqrt(gnrm)
 
248
             call dscal(nq,x,w(lw),1)
 
249
             call outl2(14,nq,nq,q,w(lw),x,x)
 
250
 
 
251
             phi0= abs(phi(q,nq,q(ltg),ng,w(lw)))
 
252
             call feq(neq,t,q,w(lqdot))
 
253
             call outl2(17,nq,nq,q,w(lqdot),phi0,x)
 
254
          endif
 
255
          goto 210
 
256
        else
 
257
          if (info .gt. 0) call outl2(39,nq,nq,xx,xx,x,x)
 
258
          nch = 17
 
259
          return
 
260
        endif
 
261
      endif
 
262
C
 
263
 299  return
 
264
C
 
265
      end
 
266