1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2c Copyright (C) INRIA - M Cardelli, L Baratchart  INRIA sophia-Antipolis 1989, S Steer
3c
4c Copyright (C) 2012 - 2016 - Scilab Enterprises
5c
6c This file is hereby licensed under the terms of the GNU GPL v2.0,
7c pursuant to article 5.3.4 of the CeCILL v.2.1.
8c This file was originally licensed under the terms of the CeCILL v2.1,
9c and continues to be available under such terms.
10c For more information, see the COPYING file which you should have received
11c along with this program.
12
13      subroutine arl2(f,nf,num,tq,dgmin,dgmax,errl2,w,iw,inf,ierr,ilog)
14C!but
15C     Cette procedure a pour but de gerer l'execution dans
16C     le cas ou un unique polynome approximant est desire
17C!liste d'appel
18C      subroutine arl2(f,nf,num,tq,dgmin,dgmax,errl2,w,
19C     $     inf,ierr,ilog)
20C
21C     double precision tq(dgmax+1),f(nf),num(dgmax)
22C     double precision w(*)
23C     integer dgmin,dgmax,dginit,info,ierr,iw(*)
24C
25C     Entree :
26C     dgmin. est le degre du polynome de depart quand il est
27C        fourni, (vaux 0 s'il ne l'est pas).
28C     dginit. est le premier degre pour lequel aura lieu la
29C        recherche.
30C     dgmax. est le degre desire du dernier approximant
31C     tq. est le tableau contenant le polynome qui peut etre
32C        fourni comme point de depart par l'utilisateur.
33C
34C     Sortie :
35C     tq. contient la solution obtenu de degre dgmax.
36C     num. contient les coefficients du numerateur optimal
37C     errl2. contient l'erreur L2 pour l'optimum retourne
38C     ierr. contient l'information sur le deroulement du programme
39C          ierr=0 : ok
40C          ierr=3 : boucle indesirable sur 2 ordres
41C          ierr=4 : plantage lsode
42C          ierr=5 : plantage dans recherche de l'intersection avec une face
43C
44C  tableau de travail
45C     w: dimension:  32+32*dgmax+7*ng+dgmax*ng+dgmax**2*(ng+2)
46C     iw : dimension  29+dgmax**2+4*dgmax
47C!sous programme appeles
48C  optml2,feq,jacl2,outl2,lq,phi (arl2)
49C  dcopy,dnrm2,dscal,dpmul1
50C!organigramme
51C arl2
52C    optml2
53C          outl2
54C          feq
55C          domout
56C                 onface
57C                        rootgp
58C                        feq
59C                        outl2
60C                 outl2
61C                 phi
62C                 lsode
63C                 front
64C                 watfac
65C          front
66C          lsode
67C                 feq
68C                 jacl2
69C                        hessl2
70C                              lq
71C    outl2
72C          feq
73C          phi
74C          lq
75C          jacl2
76C    phi
77C          lq
78C          calsca
79C    feq
80C          lq
81C          calsca
82C    lq
83C!
84c     Copyright INRIA
85      integer dgmin,dgmax,dginit,info,ierr,iw(*)
86      double precision tq(dgmax+1),f(nf),num(dgmax),w(*),x
87C
88      double precision errl2,xx(1)
89      double precision tps(2),tms(2),dnrm2,sqrt,phi,gnrm,phi0
90      integer dg,dgback,dgr
91      external feq, jacl2
92cDEC$ ATTRIBUTES DLLIMPORT:: /sortie/
93cDEC$ ATTRIBUTES DLLIMPORT:: /no2f/
94      common /sortie/ io,info,ll
95      common /no2f/ gnrm
96C
97c     taille des tableaux de travail necessaires a lsode
98      lrw = dgmax**2 + 9*dgmax + 22
99      liw = 20+dgmax
100
101C     decoupage du tableau de travail w
102      ncoeff=nf
103      ng=nf-1
104      ltq   = 1
105      ltg   = ltq+dgmax+1
106      lwode = ltg+ng+1
107      ltr   = lwode+5+5*dgmax+5*ng+dgmax*ng+dgmax**2*(ng+1)
108      lfree = ltr + 25+26*dgmax+ng+dgmax**2
109
110c     les lrw elements de w suivant w(ltr) ne doivent pas etre modifies
111c     d'un appel de optml2 a l'autre
112      lw=ltr+lrw
113C
114C     decoupage du tableau de travail iw
115      liwode = 1
116      liww   = liwode+4+(dgmax+1)*(dgmax+2)
117      lifree = liww+20+dgmax
118      iw(liwode+1)=ng
119      iw(liwode+2)=dgmax
120      ll   = 80
121      info = inf
122      io   = ilog
123C
124C test validite des arguments
125C
126      if (dgmin .gt. 0) then
127        dginit = dgmin
128        call dcopy(dgmin+1,tq,1,w(ltq),1)
129      else
130        w(ltq) = 1.d0
131        dginit = 1
132      endif
133C
134      dgr=dginit
135      ierr = 0
136      ntest1 = -1
137C
138      ng = nf - 1
139      call dcopy(nf,f,1,w(ltg),1)
140      gnrm = dnrm2(nf,f,1)
141      call dscal(nf,1.0d+0/gnrm,w(ltg),1)
142      gnrm = gnrm**2
143C
144      tps(1) = 1.0d+0
145      tps(2) = 1.0d+0
146      tms(1) = -1.0d+0
147      tms(2) = 1.0d+0
148C
149C     ---- Boucle de calcul ---------------------------------------------
150C
151      do 500 nnn = dginit,dgmax
152C
153        ifaceo = 0
154C
155        if (nnn .eq. dginit) then
156          if (dgmin .gt. 0) then
157            dg = dginit
158            goto 230
159          else
160            dg = dginit - 1
161          endif
162        endif
163C
164 200    dg = dg + 1
165C
166C     -- Initialisation du nouveau point de depart --
167C     (dans l'espace de dimension dg , Hyperespace superieur
168C     d'une dimension par rapport au precedent ).
169C
170        if (ntest1 .eq. 1) then
171          call dpmul1(w(ltq),dg-1,tps,1,w(ltr))
172          call dcopy(dg+1,w(ltr),1,w(ltq),1)
173        elseif (ntest1 .eq. -1) then
174          call dpmul1(w(ltq),dg-1,tms,1,w(ltr))
175          call dcopy(dg+1,w(ltr),1,w(ltq),1)
176        endif
177C
178C     ------------------------
179C
180 230    dgback = dg
181C
182        if (info .gt. 1) call outl2(20,dg,dgback,xx,xx,x,x)
183C
184        nch = 1
185        iw(liwode)=dg
186        call optml2(feq,jacl2,iw(liwode),w(ltq),nch,w(ltr),iw)
187        dg=iw(liwode)
188        if (info .gt. 1) then
189             call lq(dg,w(ltq),w(lw),w(ltg),ng)
190             x=sqrt(gnrm)
191             call dscal(dg,x,w(lw),1)
192             call outl2(nch,dg,dg,w(ltq),w(lw),x,x)
193
194             phi0= abs(phi(w(ltq),dg,w(ltg),ng,w(lw)))
195             lqdot=lw
196             call feq(iw(liwode),t,w(ltq),w(lqdot))
197             call outl2(17,dg,dg,w(ltq),w(lqdot),phi0,x)
198        endif
199
200        if (nch .ge. 15) then
201          if(nch.eq.17) then
202             call dcopy(dg+1,w(ltq),1,tq,1)
203             dgr=dg
204             goto 231
205          endif
206          ierr = 4 + nch - 15
207          goto 510
208        endif
209C
210        if (nch .lt. 0) then
211          ifaceo = ifaceo + 1
212          ntest1 = (-1) * ntest1
213          if (dg .eq. 0) goto 200
214          goto 230
215        endif
216C
217        if (info .gt. 1) call outl2(21,dg,dg,xx,xx,x,x)
218        nch = 2
219        iw(liwode)=dg
220        call optml2(feq,jacl2,iw(liwode),w(ltq),nch,w(ltr),iw)
221        if (info .gt. 0) then
222             call lq(dg,w(ltq),w(lw),w(ltg),ng)
223             x=sqrt(gnrm)
224             call dscal(dg,x,w(lw),1)
225             call outl2(nch,dg,dg,w(ltq),w(lw),x,x)
226
227             phi0= abs(phi(w(ltq),dg,w(ltg),ng,w(lw)))
228             lqdot=lw
229             call feq(iw(liwode),t,w(ltq),w(lqdot))
230             call outl2(17,dg,dg,w(ltq),w(lqdot),phi0,x)
231        endif
232        if (nch .ge. 15) then
233          if(nch.eq.17) then
234             call dcopy(dg+1,w(ltq),1,tq,1)
235             dgr=dg
236             goto 231
237          endif
238          ierr = 4 + nch - 15
239          goto 510
240        endif
241C
242        if (nch .lt. 0) then
243          ifaceo = ifaceo + 1
244          ntest1 = (-1) * ntest1
245          if (dg .eq. 0) goto 200
246          goto 230
247        endif
248C
249C
250 231    if (ifaceo .eq. 8) then
251          if (info .ge. 0) call outl2(22,dg,dg,xx,xx,x,x)
252          ierr = 3
253          goto 510
254        endif
255C
256        if (dg .lt. nnn) goto 200
257        call dcopy(dg+1,w(ltq),1,tq,1)
258        dgr=dg
259C
260 500  continue
261C
262C Fin de la recherche Optimale
263C numerateur optimal
264 510  gnrm = sqrt(gnrm)
265      call lq(dgr,tq,w(ltr),w(ltg),ng)
266      call dcopy(dgr,w(ltr),1,num,1)
267      call dscal(dgr,gnrm,num,1)
268C     Le gradient de la fonction critere y vaut :-tqdot
269C     call feq(dg,t,w(ltq),tqdot)
270C     valeur du critere
271      lw = ltg+ncoeff+1
272      errl2 = sqrt(phi(tq,dgr,w(ltg),ng,w(lw))) * gnrm
273      dgmax=dgr
274C
275      return
276      end
277
278