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