@@ -0,0 +1,20 @@
Loading
1 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c Copyright (C) INRIA
3 +
c 
4 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine proj(n,binf,bsup,x)
14 +
      implicit double precision (a-h,o-z)
15 +
      dimension binf(n),bsup(n),x(n)
16 +
      do i=1,n
17 +
         x(i)=max(binf(i),min(x(i),bsup(i)))
18 +
      end do
19 +
      return
20 +
      end

@@ -0,0 +1,78 @@
Loading
1 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c Copyright (C) INRIA
3 +
c 
4 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine satur (n,x,binf,bsup,d,ttmin,ttsup,topt,tg,td,
14 +
     &                   tmi,icoi,icos,iproj)
15 +
c
16 +
c      subroutine calculant ,dans un intervalle donne, un pas proche
17 +
c      de tmi saturant une contrainte
18 +
c         topt:pas calculer (=0 s'il n'existe pas un tel pas         (s)
19 +
c        ttmin,ttsup:bornes de l'intervalle dans lequel doit
20 +
c         etre topt                                                  (e)
21 +
c        tmi:pas au voisinnage duquel on calcul topt                 (e)
22 +
c        iproj:indicateur de projection                              (e)
23 +
c             =0:on cherche un pas saturant une contrainte dans
24 +
c                 l'intervalle ttmin,ttsup
25 +
c             =1:on cherche un pas dans l'intervalle tg,td et on
26 +
c                le ramene dans l'intervalle ttmin,ttsup
27 +
c                par projection
28 +
c       icos:indice de la variable saturee par la borne superieure
29 +
c       icoi:indice de la variable saturee par la borne inferieure
30 +
c       inf:indicateur pour l initialisation de icoi et icos
31 +
c            =0:la borne superieure est atteinte
32 +
c            =1:la borne superieure est atteinte
33 +
c            =2:le pas est obtenu par projection sur ttmin ttsup
34 +
c
35 +
      implicit double precision(a-h,o-z)
36 +
      integer iproj
37 +
      dimension x(n),binf(n),bsup(n),d(n)
38 +
c
39 +
      icoi=0
40 +
      icos=0
41 +
      ep=tmi
42 +
c
43 +
c        boucle
44 +
      do 70 i=1,n
45 +
      inf=0
46 +
c        calcul du pas saturant la ieme contrainte:tb
47 +
      CRES=d(i)
48 +
      if (CRES .lt. 0) then
49 +
         goto 61
50 +
      elseif (CRES .eq. 0) then
51 +
         goto 70
52 +
      else
53 +
         goto 62
54 +
      endif
55 +
61    tb=(binf(i)-x(i))/d(i)
56 +
      inf=1
57 +
      go to 63
58 +
62    tb=(bsup(i)-x(i))/d(i)
59 +
63    if ((tb.gt.ttsup).or.(tb.lt.ttmin))then
60 +
c        projection de tb sur l intervalle ttmin,ttsup
61 +
        if ((iproj.eq.0).or.(tb.lt.tg).or.(tb.gt.td)) go to 70
62 +
        tb=max(tb,ttmin)
63 +
        tb=min(tb,ttsup)
64 +
        inf=2
65 +
      end if
66 +
c        recherche du pas le plus proche de tmi
67 +
      e=abs(tb-tmi)
68 +
      if (e.ge.ep) go to  70
69 +
      topt=tb
70 +
      ep=e
71 +
c        mise a jour de icoi,icos
72 +
      icoi=0
73 +
      icos=0
74 +
      if (inf.eq.0) icos=i
75 +
      if (inf.eq.1) icoi=i
76 +
70    continue
77 +
      return
78 +
      end

@@ -6,12 +6,11 @@
Loading
6 6
Rcpp::EvalBase *gev = NULL;                  // pointer to abstract base class
7 7
8 8
typedef void (*S2_fp) (int *, int *, double *, double *, double *, int *, float *, double *);
9 -
10 -
extern "C" void n1qn1c_ (S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
9 +
extern "C" void n1qn1_ (S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
11 10
                        int mode[], int niter[], int nsim[], int imp[], double zm[], int izs[], float rzs[], double dzs[]);
12 11
13 -
unsigned int nq1n1c_calls = 0, nq1n1c_grads = 0;
14 -
int nq1n1c_fprint = 0;
12 +
unsigned int n1qn1_calls = 0, n1qn1_grads = 0;
13 +
int n1qn1_fprint = 0;
15 14
static void fwrap(int *ind, int *n, double *x, double *f, double *g, int *ti, float *tr, double *td)
16 15
{
17 16
  int i;
@@ -19,17 +18,17 @@
Loading
19 18
  std::copy(&x[0], &x[0]+*n, &par[0]);
20 19
  
21 20
  if (*ind==2 || *ind==4) {
22 -
    nq1n1c_calls++;
21 +
    n1qn1_calls++;
23 22
    ret = fev->eval(par);
24 -
    if (nq1n1c_fprint){
25 -
      Rprintf("%3d:%#14.8g:", nq1n1c_calls, ret[0]);
23 +
    if (n1qn1_fprint){
24 +
      Rprintf("%3d:%#14.8g:", n1qn1_calls, ret[0]);
26 25
      for (i = 0; i < *n; i++) Rprintf(" %#8g", x[i]);
27 26
      Rprintf("\n");
28 27
    }
29 28
    *f = ret[0];
30 29
  }
31 30
  if (*ind==3 || *ind==4) {
32 -
    nq1n1c_grads++;
31 +
    n1qn1_grads++;
33 32
    ret = gev->eval(par);
34 33
    std::copy(&ret[0],&ret[0]+*n,&g[0]);
35 34
    // for (i = 0; i < *n; i++) g[i] = ret[i];
@@ -47,14 +46,15 @@
Loading
47 46
}
48 47
49 48
50 -
RcppExport SEXP n1qn1c_wrap(
49 +
RcppExport SEXP
50 +
n1qn1_wrap(
51 51
           SEXP fSEXP, SEXP gSEXP, SEXP rhoSEXP, SEXP xSEXP, SEXP epsSEXP, 
52 52
           SEXP nSEXP, SEXP modeSEXP, SEXP niterSEXP, SEXP nsimSEXP, SEXP impSEXP,
53 53
           SEXP nzmSEXP, SEXP zmSEXP, SEXP fprint_sexp) {
54 54
  BEGIN_RCPP
55 -
    nq1n1c_calls=0;
56 -
  nq1n1c_grads=0;
57 -
  nq1n1c_fprint = INTEGER(fprint_sexp)[0];
55 +
    n1qn1_calls=0;
56 +
  n1qn1_grads=0;
57 +
  n1qn1_fprint = INTEGER(fprint_sexp)[0];
58 58
  if (TYPEOF(fSEXP) == EXTPTRSXP){
59 59
    fev = new Rcpp::EvalCompiled(fSEXP, rhoSEXP); // xptr
60 60
  } else {
@@ -86,7 +86,7 @@
Loading
86 86
  eps = REAL(epsSEXP)[0];
87 87
  std::fill(&var[0], &var[0]+n, 0.1);
88 88
  
89 -
  n1qn1c_(fwrap,&n,x,&f,g,var,&eps,
89 +
  n1qn1_(fwrap,&n,x,&f,g,var,&eps,
90 90
         &mode,&niter,&nsim,&imp,zm,izs,rzs,dzs);
91 91
        
92 92
  Rcpp::NumericVector par(n);
@@ -124,8 +124,8 @@
Loading
124 124
                            Rcpp::Named("H") = H,
125 125
                            // Rcpp::Named("zm")=zms,
126 126
                            Rcpp::Named("c.hess") = hess,
127 -
			    Rcpp::Named("n.fn") = nq1n1c_calls,
128 -
			    Rcpp::Named("n.gr") = nq1n1c_grads);
127 +
			    Rcpp::Named("n.fn") = n1qn1_calls,
128 +
			    Rcpp::Named("n.gr") = n1qn1_grads);
129 129
  END_RCPP
130 130
 }
131 131
@@ -229,9 +229,9 @@
Loading
229 229
	  SEXP itmaxSEXP, SEXP epsfSEXP, SEXP epsgSEXP, SEXP epsxSEXP, SEXP nSEXP,
230 230
          SEXP fprint_sexp) {
231 231
  BEGIN_RCPP
232 -
    nq1n1c_calls=0;
233 -
  nq1n1c_grads=0;
234 -
  nq1n1c_fprint = INTEGER(fprint_sexp)[0];
232 +
    n1qn1_calls=0;
233 +
  n1qn1_grads=0;
234 +
  n1qn1_fprint = INTEGER(fprint_sexp)[0];
235 235
    int indqn[1];
236 236
  if (TYPEOF(fSEXP) == EXTPTRSXP){
237 237
    fev = new Rcpp::EvalCompiled(fSEXP, rhoSEXP); // xptr
@@ -331,7 +331,7 @@
Loading
331 331
			    // Rcpp::Named("c.hess") = hess,
332 332
			    Rcpp::Named("code") = indqn[0],
333 333
                            Rcpp::Named("codeDesc") = ret,
334 -
			    Rcpp::Named("n.fn") = nq1n1c_calls,
335 -
			    Rcpp::Named("n.gr") = nq1n1c_grads);
334 +
			    Rcpp::Named("n.fn") = n1qn1_calls,
335 +
			    Rcpp::Named("n.gr") = n1qn1_grads);
336 336
  END_RCPP
337 337
    }

@@ -0,0 +1,183 @@
Loading
1 +
c     removed io and comment out any printing from fortran.
2 +
3 +
c     Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 +
c     Copyright (C) 1986 - INRIA - F. BONNANS
5 +
c 
6 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
7 +
c
8 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
9 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
10 +
c This file was originally licensed under the terms of the CeCILL v2.1,
11 +
c and continues to be available under such terms.
12 +
c For more information, see the COPYING file which you should have received
13 +
c along with this program.
14 +
c
15 +
      subroutine qnbd(indqn,simul,n,x,f,g,zero,
16 +
     &     napmax,itmax,epsf,epsg,epsx,df0,binf,bsup,nfac,
17 +
     & trav,ntrav,itrav,nitrav,izs,rzs,dzs)
18 +
c!but
19 +
c     code de minimisation d une fonction reguliere sous contraintes
20 +
c     de bornes , aux normes modulopt
21 +
c!methode
22 +
c     principe de l algorithme : quasi-newton + projection
23 +
c     details dans le rapport inria n. 242,1983
24 +
c     cette version permet de tester plusieurs variantes de l algorithme
25 +
c     en modifiant certains parametres internes (cf comment dans le code)
26 +
c     taille memoire necessaire de l ordre de n**2/2
27 +
c     pour les problemes de grande taille le code gcbd est mieux adapte
28 +
c
29 +
c!sous programmes appeles
30 +
c     zqnbd   optimiseur effectif
31 +
c     proj    projection
32 +
c     calmaj  mise a jour du hessien 
33 +
c     ajour mise a jour des facteurs de choleski 
34 +
c     rlbd,satur   recherche lineaire avec bornes
35 +
c
36 +
c!liste d'appel
37 +
c
38 +
c     subroutine qnbd(indqn,simul,n,x,f,g,iprint,io,zero,
39 +
c    & napmax,itmax,epsf,epsg,epsx,df0,binf,bsup,nfac,
40 +
c    & trav,ntrav,itrav,nitrav,izs,rzs,dzs)
41 +
c
42 +
c     indqn   indicateur de qnbd                                  es
43 +
c       en entree =1 standard
44 +
c                 =2 dh et indic initialises au debut de trav et itrav
45 +
c                    ifac,f,g initialises
46 +
c       en sortie
47 +
c        si < 0 incapacite de calculer un point meilleur que le point initial
48 +
c        si = 0 arret demande par l utilisateur
49 +
c        si > 0 on fournit un point meilleur que le point de depart
50 +
c        < -10 parametres d entree non convenables
51 +
c        = -6  arret lors du calcul de la direction de descente et iter=1
52 +
c        = -5  arret lors du calcul de l approximation du hessien  iter=1
53 +
c        = -3  anomalie de simul : indic negatif en un point ou
54 +
c              f et g ont ete precedemment calcules
55 +
c        = -2  echec de la recherche lineaire a la premiere iteration
56 +
c        = -1  f non definie au point initial
57 +
c        =  1  arret sur epsg
58 +
c        =  2            epsf
59 +
c        =  3            epsx
60 +
c        =  4            napmax
61 +
c        =  5            itmax
62 +
c        =  6  pente dans la direction opposee au gradient trop petite
63 +
c        =  7  arret lors du calcul de la direction de descente
64 +
c        =  8  arret lors du calcul de l approximation du hessien
65 +
c        = 10  arret par echec de la recherche lineaire , cause non precisee
66 +
c        = 11  idem avec indsim < 0
67 +
c        = 12            un pas trop petit proche d un pas trop grand
68 +
c                        ceci peut resulter d une erreur dans le gradient
69 +
c        = 13            trop grand nombre d appels dans une recherche lineaire
70 +
c     simul voir les normes modulopt
71 +
c     n dim de x                                                 e
72 +
c     binf,bsup bornes inf,sup,de dim n                          e
73 +
c     x variables a optimiser (controle)                          es
74 +
c     f valeur du critere                                         s
75 +
c     g gradient de f                                             s
76 +
c     zero  proche zero machine                                             e
77 +
c     napmax nombre maximum d appels de simul                               e
78 +
c     itmax nombre maximum d iterations de descente               e
79 +
c     itrav vect travail dim nitrav=2n , se decompose en indic et izig
80 +
c     nfac nombre de variables factorisees                  (e si indqn=2)  s
81 +
c     iprint facteur d impression                                              e
82 +
c     varie de 0 (pas d impressions) a 3 (nombreuses impressions)
83 +
c     io numero du fichier de resultats                                     e
84 +
c     epsx vect dim n precision sur x                                       e
85 +
c     epsf critere arret sur f                                              e
86 +
c     epsg arret si sup a norm2(g+)/n                                       e
87 +
c     trav vect travail dim ntrav
88 +
c     il faut ntrav > n(n+1)/2 +6n
89 +
c     df0>0 decroissance f prevue     (prendre 1. par defaut)               e
90 +
c     izs,rzs,dzs : cf normes modulopt                                     es
91 +
c!
92 +
c     indications sur les variables internes a qnbd et zqnbd
93 +
c     izig  sert a la memorisation des contraintes (actif si izag>1)
94 +
c     si i ne change pas d ens on enleve 1 a izig  (positif)
95 +
c     sinon on ajoute izag
96 +
c     factorisation seulement si izig est nul
97 +
c     dh estimation hessien dim n(n+1)/2 rangee en troismorceaux es
98 +
c     indic(i) nouvel indice de l indice i
99 +
c     indic vect dim n ordre de rangement des indices                       es
100 +
c     pas necessaire de l initialiser si indqn=1
101 +
c
102 +
c     parametres de la recherche lineaire
103 +
c     amd,amf param. du test de wolfe .    (.7,.1)
104 +
c     napm nombre max d appels dans la rl  (=15)
105 +
c
106 +
      implicit double precision (a-h,o-z)
107 +
      real rzs(*)
108 +
      double precision dzs(*)
109 +
      dimension binf(n),bsup(n),x(n),g(n),epsx(n)
110 +
      dimension trav(ntrav),itrav(nitrav),izs(*)
111 +
      external simul
112 +
c
113 +
c---- initial printing
114 +
c$$$      if(iprint.ge.1) then
115 +
c$$$         call basout(io_out, io, '')
116 +
c$$$         write(bufstr,1010)
117 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
118 +
c$$$         write(bufstr,750) n,epsg,iprint
119 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
120 +
c$$$         write(bufstr,751) itmax
121 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
122 +
c$$$         write(bufstr,752) napmax
123 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
124 +
c$$$         call basout(io_out ,io ,
125 +
c$$$     $    '------------------------------------------------')
126 +
c$$$1010    format(' *********** qnbd (with bound cstr) ****************')
127 +
c$$$750     format('dimension=',i10,', epsq=',e24.16,
128 +
c$$$     $ ', verbosity level: iprint=',i10)
129 +
c$$$751     format('max number of iterations allowed: iter=',i10)
130 +
c$$$752     format('max number of calls to costf allowed: nap=',i10)
131 +
c$$$      endif
132 +
c
133 +
c
134 +
c     parametres caracteristiques de l algorithme
135 +
c     si les parametres sont nuls l algorithme est celui du rr 242
136 +
c     ig=1 test sur grad(i) pour relach var
137 +
c     in=1 limite le nombre de factorisations par iter a n/10
138 +
c     irel=1 test sur decroissance grad pour rel a iter courante
139 +
c     epsrel taux de decroissance permettant le relachement (cf irit)
140 +
c     iact blocage variables dans ib (gestion contraintes actives)
141 +
c     ieps1 =1 eps1 egal a zero
142 +
c     note eps1 correspond a eps(xk)
143 +
      ig=0
144 +
      in=0
145 +
      irel=1
146 +
      epsrel=.50d+0
147 +
      izag=0
148 +
      iact=1
149 +
      ieps1=0
150 +
c
151 +
c     decoupage du vecteur trav
152 +
      n1=(n*(n+1))/2 +1
153 +
      n2=n1+n
154 +
      n3=n2+n
155 +
      n4=n3+n
156 +
      n5=n4+n-1
157 +
      if(ntrav.lt.n5) then
158 +
c$$$         if(iprint.gt.0) then
159 +
c$$$           write(bufstr,110)ntrav,n5
160 +
c$$$           call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
161 +
c$$$           endif
162 +
c$$$110      format(' qnbd : ntrav=',i8,' devrait valoir ',i8)
163 +
         indqn=-11
164 +
         return
165 +
      endif
166 +
      ni1=n+1
167 +
      if(nitrav.lt.2*n) then
168 +
         ni2=2*n
169 +
c$$$         if(iprint.gt.0) then
170 +
c$$$           write(bufstr,111)nitrav,ni2
171 +
c$$$           call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
172 +
c$$$           endif
173 +
c$$$111      format(' qnbd : nitrav=',i8,'devrait valoir',i8)
174 +
         indqn=-12
175 +
         return
176 +
      endif
177 +
      call zqnbd(indqn,simul,trav(1),n,binf,bsup,x,f,g,zero,napmax,
178 +
     &itmax,itrav,itrav(ni1),nfac,epsx,epsf,epsg,trav(n1),
179 +
     &trav(n2),trav(n3),trav(n4),df0,ig,in,irel,izag,iact,
180 +
     &epsrel,ieps1,izs,rzs,dzs)
181 +
      return
182 +
      end
183 +
      

@@ -0,0 +1,278 @@
Loading
1 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c Copyright (C) INRIA
3 +
c 
4 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine ajour(mode,n,nc,nr,h,w,indi)
14 +
c
15 +
      implicit double precision (a-h,o-z)
16 +
      dimension h(*),w(n),indi(n)
17 +
c
18 +
c     mode = +1 factorise la ligne nc (indices de depart)
19 +
c          = -1 defactorise   '
20 +
c     nr nbre de lignes factorisees
21 +
c     h mat de dim n
22 +
c     w,d vect de travail
23 +
c     indi(i) ligne ou est stockee la ligne i de depart
24 +
c
25 +
      inc=indi(nc)
26 +
      nr1=nr+1
27 +
      nr2=nr-1
28 +
      nrr=n-nr
29 +
      nii=n-inc
30 +
      nkk=nr-inc
31 +
      if(mode.eq.-1)go to 240
32 +
c
33 +
c      addition d'une ligne a l
34 +
c
35 +
c          stockage des elements de la colonne inc dans w
36 +
      nsaut=nii+1
37 +
      nh=inc*(n+1)-inc*(inc+1)/2
38 +
      nw=n
39 +
      if(inc.eq.n) go to 20
40 +
      do i=1,nii
41 +
         w(nw)=h(nh)
42 +
         nw=nw-1
43 +
         nh=nh-1
44 +
      end do
45 +
   20 w(nr1)=h(nh)
46 +
      nh=nh-1
47 +
      if(inc.eq.nr1) go to 60
48 +
      do i=1,inc-nr1
49 +
         nl=nii+i-1
50 +
         if(nl.eq.0) go to 35
51 +
         do j=1,nl
52 +
            h(nh+nsaut)=h(nh)
53 +
            nh=nh-1
54 +
         end do 
55 +
 35      w(nw)=h(nh)
56 +
         nw=nw-1
57 +
         nh=nh-1
58 +
         nsaut=nsaut+1
59 +
      end do
60 +
      do j=1,inc-nr1
61 +
         h(nh+nsaut)=h(nh)
62 +
         nh=nh-1
63 +
      end do
64 +
c
65 +
   60 nw=nw-1
66 +
      nsaut=1
67 +
      if(nr.eq.0) go to 125
68 +
      if(inc.eq.n) go to 80
69 +
      do i=1,nii
70 +
         h(nh+nsaut)=h(nh)
71 +
         nh=nh-1
72 +
      end do
73 +
   80 if(nr.eq.1) go to 110
74 +
      do i=1,nr2
75 +
         w(nw)=h(nh)
76 +
         nw=nw-1
77 +
         nh=nh-1
78 +
         nsaut=nsaut+1
79 +
         if(n.eq.nr1) go to 100
80 +
         do  j=1,n-nr1
81 +
            h(nh+nsaut)=h(nh)
82 +
            nh=nh-1
83 +
         end do
84 +
      end do
85 +
  100 continue
86 +
  110 w(nw)=h(nh)
87 +
      nh=nh-1
88 +
      nsaut=nsaut+1
89 +
      if(inc.eq.nr1) go to 125
90 +
      do i=1,inc-nr1
91 +
         h(nh+nsaut)=h(nh)
92 +
         nh=nh-1
93 +
      end do
94 +
c         mise a jour de l
95 +
  125 if(nr.ne.0) go to 130
96 +
      if(w(1).gt.0.0d+0) go to 220
97 +
      mode=-1
98 +
      return
99 +
  130 if(nr.eq.1) go to 160
100 +
      do i=2,nr
101 +
         ij=i
102 +
         i1=i-1
103 +
         v=w(i)
104 +
         do j=1,i1
105 +
            v=v-h(ij)*w(j)
106 +
            ij=ij+nr-j
107 +
         end do
108 +
         w(i)=v
109 +
      end do
110 +
  160 ij=1
111 +
      v=w(nr1)
112 +
      do i=1,nr
113 +
         wi=w(i)
114 +
         hij=h(ij)
115 +
         v=v-(wi**2)/hij
116 +
         w(i)=wi/hij
117 +
         ij=ij+nr1-i
118 +
      end do
119 +
      if(v.gt.0.0d+0) go to 180
120 +
      mode=-1
121 +
      return
122 +
  180 w(nr1)=v
123 +
c          stockage de w dans h
124 +
      nh=nr*(nr+1)/2
125 +
      nw=nr1
126 +
      nsaut=nw
127 +
      h(nh+nsaut)=w(nw)
128 +
      nw=nw-1
129 +
      nsaut=nsaut-1
130 +
      if(nr.eq.1) go to 220
131 +
      do i=1,nr2
132 +
         h(nh+nsaut)=w(nw)
133 +
         nw=nw-1
134 +
         nsaut=nsaut-1
135 +
         do  j=1,i
136 +
            h(nh+nsaut)=h(nh)
137 +
            nh=nh-1
138 +
         end do
139 +
      end do
140 +
  220 h(nr1)=w(1)
141 +
      if(n.eq.nr1) go to 233
142 +
      nh1=nr*(n+1)-nr*(nr+1)/2+1
143 +
      nw=nr1
144 +
      do i=1,n-nr1
145 +
         h(nh1+i)=w(nw+i)
146 +
      end do
147 +
c          mise a jour de indi
148 +
 233  do i=1,n
149 +
         ii=indi(i)
150 +
         if(ii.le.nr.or.ii.ge.inc) go to 235
151 +
         indi(i)=ii+1
152 +
      end do
153 +
  235 continue
154 +
      nr=nr+1
155 +
      indi(nc)=nr
156 +
      mode=0
157 +
      return
158 +
c
159 +
c      soustraction d'une ligne a l
160 +
c
161 +
c          recherche des composantes de h
162 +
 240  do i=1,nr
163 +
         ik=i
164 +
         ij=inc
165 +
         ii=1
166 +
         ko=min(ik,inc)
167 +
         v=0.0d+0
168 +
         if(ko.eq.1) go to 252
169 +
         do k=1,ko-1
170 +
            nk=nr1-k
171 +
            v=v+h(ij)*h(ik)*h(ii)
172 +
            ij=ij+nk-1
173 +
            ii=ii+nk
174 +
            ik=ik+nk-1
175 +
         end do
176 +
 252     a=1.0d+0
177 +
         b=1.0d+0
178 +
         if(ko.eq.i) go to 253
179 +
         a=h(ik)
180 +
 253     if(ko.eq.inc) go to 260
181 +
         b=h(ij)
182 +
         w(i)=v+a*b*h(ii)
183 +
      end do
184 +
 260  continue
185 +
c          mise a jour de l
186 +
      if(inc.eq.nr) go to 315
187 +
      inc1=inc-1
188 +
      nh=inc1*nr1-inc1*inc/2+2
189 +
      nh1=nh+nkk
190 +
      di=h(nh-1)
191 +
      do j=1,nkk
192 +
         di1=h(nh1)
193 +
         nh1=nh1+1
194 +
         a=h(nh)
195 +
         ai=a*di
196 +
         c=(a**2)*di+di1
197 +
         h(nh)=c
198 +
         nh=nh+1
199 +
         if(j.eq.nkk) go to 315
200 +
         do i=1,nkk-j
201 +
            h1=h(nh)
202 +
            h2=h(nh1)
203 +
            u=ai*h1+h2*di1
204 +
            h(nh)=u/c
205 +
            h(nh1)=-h1+a*h2
206 +
            nh=nh+1
207 +
            nh1=nh1+1
208 +
         end do
209 +
         nh=nh+1
210 +
         di=di*di1/c
211 +
      end do
212 +
c          condensation de la matrice l
213 +
  315 nh=inc+1
214 +
      nsaut=1
215 +
      nj=nr-2
216 +
      if(inc.eq.1) nj=nj+1
217 +
      if(nr.eq.1) go to 440
218 +
      do i=1,nr2
219 +
         do j=1,nj
220 +
            h(nh-nsaut)=h(nh)
221 +
            nh=nh+1
222 +
         end do
223 +
         nsaut=nsaut+1
224 +
         nh=nh+1
225 +
         if(i.eq.inc-1) go to 430
226 +
         nj=nj-1
227 +
         if(nj.eq.0) go to 440
228 +
      end do
229 +
  430 continue
230 +
c          mise a jour de la matrice h
231 +
  440 nh=((nr*nr2)/2)+1
232 +
      nw=1
233 +
      nsaut=nr
234 +
      if(inc.eq.1) go to 470
235 +
      do i=1,inc-1
236 +
         h(nh)=w(nw)
237 +
         nw=nw+1
238 +
         nsaut=nsaut-1
239 +
         if(n.eq.nr) go to 455
240 +
         do j=1,nrr
241 +
            h(nh+j)=h(nh+nsaut+j)
242 +
         end do
243 +
 455     nh=nh+nrr+1
244 +
      end do
245 +
  470 nw=nw+1
246 +
      if(nr.eq.n) go to 485
247 +
      do i=1,nrr
248 +
         w(nr+i)=h(nh+nsaut+i-1)
249 +
      end do
250 +
      nsaut=nsaut+nrr
251 +
  485 if(inc.eq.nr) go to 510
252 +
      do i=1,nkk
253 +
         nsaut=nsaut-1
254 +
         h(nh)=w(nw)
255 +
         nw=nw+1
256 +
         if(nr.eq.n) go to 495
257 +
         do j=1,nrr
258 +
            h(nh+j)=h(nh+nsaut+j)
259 +
         end do
260 +
 495     nh=nh+nrr+1
261 +
      end do
262 +
  510 h(nh)=w(inc)
263 +
      if(nr.eq.n) go to 540
264 +
      do i=1,nrr
265 +
         h(nh+i)=w(nr+i)
266 +
      end do
267 +
c          mise a jour de indi
268 +
 540  do i=1,n
269 +
         ii=indi(i)
270 +
         if(ii.le.inc.or.ii.gt.nr) go to 550
271 +
         indi(i)=ii-1
272 +
      end do
273 +
  550 continue
274 +
      indi(nc)=nr
275 +
      nr=nr-1
276 +
      mode=0
277 +
      return
278 +
      end

@@ -4,23 +4,22 @@
Loading
4 4
#include <R_ext/Rdynload.h>
5 5
6 6
typedef void (*S2_fp) (int *, int *, double *, double *, double *, int *, float *, double *);
7 -
extern void n1qn1c_(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
7 +
extern void n1qn1_(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
8 8
		   int mode[], int niter[], int nsim[], int imp[], double zm[], int izs[], 
9 9
		   float rzs[], double dzs[]);
10 -
void n1qn1cF(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
10 +
void n1qn1F(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
11 11
		   int mode[], int niter[], int nsim[], int imp[], int lp[], double zm[], int izs[], 
12 12
		   float rzs[], double dzs[]) {
13 -
  n1qn1c_(simul, n, x, f, g, var, eps, mode, niter, nsim, imp, zm, izs, rzs, dzs);
13 +
  n1qn1_(simul, n, x, f, g, var, eps, mode, niter, nsim, imp, zm, izs, rzs, dzs);
14 14
}
15 15
16 -
void n1qn1cF2(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
16 +
void n1qn1F2(S2_fp simul, int n[], double x[], double f[], double g[], double var[], double eps[],
17 17
		   int mode[], int niter[], int nsim[], int imp[], double zm[], int izs[], 
18 18
		   float rzs[], double dzs[]) {
19 -
  n1qn1c_(simul, n, x, f, g, var, eps, mode, niter, nsim, imp, zm, izs, rzs, dzs);
19 +
  n1qn1_(simul, n, x, f, g, var, eps, mode, niter, nsim, imp, zm, izs, rzs, dzs);
20 20
}
21 21
/* .C calls */
22 -
extern SEXP n1qn1c_wrap(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
23 -
22 +
extern SEXP n1qn1_wrap(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP);
24 23
extern SEXP qnbd_wrap(SEXP, SEXP, SEXP, SEXP, SEXP, 
25 24
		      SEXP, SEXP, SEXP, SEXP, SEXP, 
26 25
		      SEXP, SEXP, SEXP, SEXP);
@@ -30,16 +29,19 @@
Loading
30 29
		  double* binf, double* binsup, int* nfac, double* trav, int* ntrav, int* itrav, int* nitrav, 
31 30
		  int* izs, float* rzs, double* dzs);
32 31
33 -
void R_init_n1qn1c(DllInfo *dll)
32 +
33 +
static const R_CallMethodDef CallEntries[] = {
34 +
  {"n1qn1_wrap", (DL_FUNC) &n1qn1_wrap, 13},
35 +
  {"qnbd_wrap", (DL_FUNC) &qnbd_wrap, 14},
36 +
  {NULL, NULL, 0}
37 +
};
38 +
39 +
void R_init_n1qn1(DllInfo *dll)
34 40
{
35 -
  R_CallMethodDef callMethods[]  = {
36 -
    {"n1qn1c_wrap", (DL_FUNC) &n1qn1c_wrap, 13},
37 -
    {"qnbd_wrap", (DL_FUNC) &qnbd_wrap, 14},
38 -
    {NULL, NULL, 0}
39 -
  };
40 -
  R_RegisterCCallable("n1qn1c","n1qn1cF", (DL_FUNC) n1qn1cF);
41 -
  R_RegisterCCallable("n1qn1c","n1qn1cF2", (DL_FUNC) n1qn1cF2);
42 -
  R_RegisterCCallable("n1qn1c","qnbdF", (DL_FUNC) qnbd_);
43 -
  R_registerRoutines(dll, NULL, callMethods, NULL, NULL);
44 -
  R_useDynamicSymbols(dll, FALSE);
45 -
} 
41 +
  R_RegisterCCallable("n1qn1","n1qn1F", (DL_FUNC) n1qn1F);
42 +
  R_RegisterCCallable("n1qn1","n1qn1F2", (DL_FUNC) n1qn1F2);
43 +
  R_RegisterCCallable("n1qn1","qnbdF", (DL_FUNC) qnbd_);
44 +
  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
45 +
  R_useDynamicSymbols(dll, TRUE);
46 +
  R_forceSymbols(dll,TRUE);
47 +
}

@@ -0,0 +1,660 @@
Loading
1 +
c     Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c     Copyright (C) INRIA
3 +
c 
4 +
c     Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine zqnbd(indqn,simul,dh,n,binf,bsup,x,f,g,zero,napmax,
14 +
     &     itmax,indic,izig,nfac,epsx,epsf,epsg,x1,x2,g1,dir,df0,
15 +
     &ig,in,irel,izag,iact,epsrel,ieps1,izs,rzs,dzs)
16 +
c
17 +
      implicit double precision (a-h,o-z)
18 +
      real rzs(*)
19 +
      double precision dzs(*)
20 +
      dimension x1(n),x2(n),g1(n),dir(n),epsx(n)
21 +
      dimension binf(n),bsup(n),x(n),g(n),dh(*),indic(n),izig(n),
22 +
     &izs(*)
23 +
      external simul,proj
24 +
c
25 +
c$$$      if(iprint.lt.4)go to 3
26 +
c$$$      write(bufstr,1020)izag,ig,in,irel,iact,epsrel
27 +
c$$$      call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
28 +
c$$$1020  format(' qnbd :  izag,ig,in,irel,iact,epsrel=',5i3,f11.4)
29 +
c$$$c
30 +
c$$$      if(ig.eq.1) then
31 +
c$$$        write(bufstr,110)
32 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
33 +
c$$$        endif
34 +
c$$$110   format(' test sur gradient pour sortie ib')
35 +
c$$$      if(in.eq.1) then 
36 +
c$$$        write(bufstr,111)
37 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
38 +
c$$$        endif
39 +
c$$$111   format(' test sur nombre de defactorisations pour sortie ib')
40 +
c$$$      if(izag.ne.0) then
41 +
c$$$        write(bufstr,112)izag
42 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
43 +
c$$$        endif
44 +
c$$$112   format(' memorisation de variables izag=',i3)
45 +
c$$$      if(irel.eq.1) then
46 +
c$$$        write(bufstr,114)epsrel
47 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
48 +
c$$$        endif
49 +
c$$$114   format(' methode de minimisations incompletes ; epsrel=',d11.4)
50 +
c$$$      if(iact.eq.1) then
51 +
c$$$        write(bufstr,116)
52 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
53 +
c$$$        endif
54 +
c$$$116   format(' blocage des variables dans ib')
55 +
c$$$      if(ieps1.eq.1) then
56 +
c$$$        write(bufstr,118)
57 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
58 +
c$$$        endif
59 +
c$$$118   format(' parametre eps1 nul')
60 +
c$$$      if(ieps1.eq.2) then
61 +
c$$$        write(bufstr,119)
62 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
63 +
c$$$        endif
64 +
c$$$119   format(' parametre eps1 grand')
65 +
c$$$c
66 +
c$$$c     cscal1 utilise pour calculer eps(x) = eps1 cf avant 310
67 +
c$$$      cscal1=1.0d+8
68 +
c$$$      if(ieps1.eq.2) then
69 +
c$$$        write(bufstr,120)cscal1
70 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
71 +
c$$$        endif
72 +
c$$$120   format(' parametre eps1=eps(x) calcule avec cscal1=',d11.4)
73 +
c$$$3     continue
74 +
c
75 +
      difg0=1.0d+0
76 +
      difg1=0.0d+0
77 +
c
78 +
c     eps0 sert a partitionner les variables
79 +
      eps0=0.0d+0
80 +
      do i=1,n
81 +
         izig(i)=0
82 +
         eps0=eps0+epsx(i)
83 +
      end do
84 +
      eps0=10.*eps0/n
85 +
c
86 +
c     section 1  mise en forme de dh
87 +
c     si indqn=1 on init dh a ident puis scal a it=2
88 +
c
89 +
      call proj(n,binf,bsup,x)
90 +
      ndh=n*(n+1)/2
91 +
      if(indqn.eq.1)go to 10
92 +
      if(indqn.eq.2)go to 30
93 +
c     erreur
94 +
c$$$      if(iprint.gt.0) then
95 +
c$$$        write(bufstr,105)indqn
96 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
97 +
c$$$        endif
98 +
c$$$105   format(' qnbd  : valeur non admissible de indqn  ',i5)
99 +
      indqn=-105
100 +
c$$$      if(iprint.gt.0) then
101 +
c$$$        write(bufstr,123)indqn
102 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
103 +
c$$$        endif
104 +
      return
105 +
10    continue
106 +
c     on initialise dh a l identite puis a l iteration 2
107 +
c     on met a l echelle
108 +
      nfac=0
109 +
      do i=1,n
110 +
         indic(i)=i
111 +
      end do
112 +
      do i=1,ndh
113 +
         dh(i)=0.0d+0
114 +
      end do 
115 +
30    continue
116 +
c
117 +
c     section 2  mise a jour dh
118 +
c
119 +
c     iter nombre d iterations de descente
120 +
      iter=0
121 +
      scal=1.0d+0
122 +
      nap=1
123 +
      indsim=4
124 +
      if(indqn.eq.1) call simul(indsim,n,x,f,g,izs,rzs,dzs)
125 +
      if(indsim.le.0)then
126 +
      indqn=-1
127 +
      if(indsim.eq.0)indqn=0
128 +
c$$$      if(iprint.gt.0) then
129 +
c$$$        write(bufstr,123)indqn
130 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
131 +
c$$$        endif
132 +
c$$$123   format(' qnbd : indqn=',i8)
133 +
      return
134 +
      endif
135 +
      if(indqn.ne.1)go to 200
136 +
c     mise a echelle dh
137 +
c     df0 decroissance prevue . si mod quad df0=((dh)-1g,g)/2
138 +
c     et on cherche dh diag de la forme cst/(dx)**2
139 +
c     d ou cst=som((y(i)*(dx))**2))/(2*df0)
140 +
      cof1=0.0d+0
141 +
      do i=1,n
142 +
         cof1=cof1+(g(i)*epsx(i))**2
143 +
      end do
144 +
      cof1=cof1/(2.0d+0*df0)
145 +
      i1=-n
146 +
      do i=1,n
147 +
         i1=i1+n+2-i
148 +
         dh(i1)=(cof1 + zero)/(epsx(i)**2 + zero)
149 +
      end do
150 +
      iconv=0
151 +
200   iter=iter +1
152 +
      if(iter.le.itmax)go to 202
153 +
c$$$      if(iprint.gt.0) then
154 +
c$$$        write(bufstr,1202)
155 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
156 +
c$$$        endif
157 +
c$$$1202  format(' qnbd : maximum d iterations atteint')
158 +
      indqn=5
159 +
c$$$      if(iprint.gt.0) then
160 +
c$$$        write(bufstr,123)indqn
161 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
162 +
c$$$        endif
163 +
      return
164 +
c$$$202   if(iprint.ge.2) then
165 +
c$$$         write(bufstr,1210)iter,f
166 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
167 +
c$$$         endif
168 +
c$$$1210  format(' qnbd : iter=',i3,'  f=',d15.7)
169 +
c     x1,g1 valeurs a l iteration precedente
170 +
202   if(iter.eq.1)go to 300
171 +
      cof1=0.0d+0
172 +
      do i=1,n
173 +
         x1(i)=x(i)-x1(i)
174 +
         g1(i)=g(i)-g1(i)
175 +
         cof1=cof1 + x1(i)*g1(i)
176 +
      end do
177 +
      if(cof1.le.zero)go to 250
178 +
      if(iter.gt.2.or.indqn.ne.1)go to 250
179 +
c     mise a l echelle de dh par methode shanno-phua
180 +
c      dh=(y,y)/(y,s)*id
181 +
      cof2=0.0d+0
182 +
      do i=1,n
183 +
         cof2=cof2 + g1(i)**2
184 +
      end do
185 +
      cof2=cof2/cof1
186 +
c$$$      if(iprint.gt.3) then
187 +
c$$$        write(bufstr,1203)cof2
188 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
189 +
c$$$        endif
190 +
c$$$1203  format(' qnbd : facteur d echelle=',d11.4)
191 +
      dh(1)=cof2
192 +
      i1=1
193 +
      do i=1,nfac
194 +
         i1=i1+n+1-i
195 +
         dh(i1)=cof2
196 +
      end do
197 +
c
198 +
c     scal= (y,s)/(y,y)
199 +
c     scal sert de coeff a g dans le calcul de dir pour i dans ib
200 +
      scal=1.0d+0/cof2
201 +
250   continue
202 +
c
203 +
c     mise a jour dh par methode bfgs (majour) si iter ge 2
204 +
c     dh1=dh +y*yt/(y,s) - dh*s*st*dh/(s,dh*s)
205 +
c     exprimons ds=x1 et y=g1 dans les nouv variables soit x2 et g1
206 +
      do i=1,n
207 +
         i1=indic(i)
208 +
         x2(i1)=g1(i)
209 +
         dir(i1)=x1(i)
210 +
      end do
211 +
      do i=1,n
212 +
         g1(i)=x2(i)
213 +
      end do
214 +
      do i=1,n
215 +
         i1=indic(i)
216 +
         x2(i1)=x1(i)
217 +
      end do
218 +
c     on stocke d abord dh*s dans x2
219 +
c     calcul des nfac premieres variables,en deux fois
220 +
      continue
221 +
      if(nfac.eq.0) go to 2312
222 +
      if(nfac.gt.1) go to 2300
223 +
      dir(1)=dir(1)*dh(1)
224 +
      go to 2312
225 +
2300  continue
226 +
      np=nfac+1
227 +
      ii=1
228 +
      n1=nfac-1
229 +
      do i=1,n1
230 +
         y=dir(i)
231 +
         if(dh(ii).eq.0.0d+0) go to 2302
232 +
         ij=ii
233 +
         ip=i+1
234 +
         do j=ip,nfac
235 +
            ij=ij+1
236 +
            y=y+dir(j)*dh(ij)
237 +
         end do
238 +
 2302    dir(i)=y*dh(ii)
239 +
         ii=ii+np-i
240 +
      end do
241 +
      dir(nfac)=dir(nfac)*dh(ii)
242 +
      do k=1,n1
243 +
         i=nfac-k
244 +
         ii=ii-np+i
245 +
         if(dir(i).eq.0.0d+0) go to 2311
246 +
         ip=i+1
247 +
         ij=ii
248 +
         y=dir(i)
249 +
         do j=ip,nfac
250 +
            ij=ij+1
251 +
            dir(j)=dir(j)+dh(ij)*dir(i)
252 +
         end do
253 +
      end do
254 +
2311  continue
255 +
2312  continue
256 +
      nfac1=nfac+1
257 +
      n2fac=(nfac*nfac1)/2
258 +
      nnfac=n-nfac
259 +
      k=n2fac
260 +
      if(nfac.eq.n)go to 268
261 +
      do i=nfac1,n
262 +
         dir(i)=0.0d+0
263 +
      end do
264 +
      if(nfac.eq.0)go to 265
265 +
      do i=1,nfac
266 +
         do j=nfac1,n
267 +
            k=k+1
268 +
            if(x2(j).eq.0.)go to 260
269 +
            dir(i)= dir(i) + dh(k)*x2(j)
270 +
         end do
271 +
      end do
272 +
260   continue
273 +
c     calcul autres comp de dh*s=d en deux fois
274 +
      k=n2fac
275 +
      do j=1,nfac
276 +
         do i=nfac1,n
277 +
            k=k+1
278 +
            dir(i)=dir(i) + dh(k)*x2(j)
279 +
         end do
280 +
      end do
281 +
265   continue
282 +
      k=n2fac+nfac*nnfac
283 +
      do j=nfac1,n
284 +
         do i=j,n
285 +
            k=k+1
286 +
            if(x2(j).eq.0.)go to 266
287 +
            dir(i)=dir(i) + dh(k)*x2(j)
288 +
         end do
289 +
      end do
290 +
266   continue
291 +
      if(nfac.eq.n-1)go to 268
292 +
      nm1=n-1
293 +
      k=n2fac+nfac*nnfac
294 +
      do 267 i=nfac1,nm1
295 +
      k=k+1
296 +
      i1=i+1
297 +
      do j=i1,n
298 +
         k=k+1
299 +
         if(x2(j).eq.0.)go to 267
300 +
         dir(i)=dir(i)+dh(k)*x2(j)
301 +
      end do
302 +
267   continue
303 +
c     calcul de dh*s fini
304 +
c     calcul sig1 pour 2eme mise a jour
305 +
268   sig1=0.0d+0
306 +
      do i=1,n
307 +
         sig1=sig1+dir(i)*x2(i)
308 +
      end do
309 +
      if(sig1.gt.0.0d+0)go to 272
310 +
c
311 +
c     ******************************************************
312 +
      indqn=8
313 +
      if(iter.eq.1)indqn=-5
314 +
272      sig1=-1.0d+0/sig1
315 +
c     truc powell si (y,s) negatif
316 +
      if(cof1.gt.zero)go to 277
317 +
      teta=-1.0d+0/sig1
318 +
      teta=.8*teta/(teta-cof1)
319 +
      teta1=1.0d+0-teta
320 +
      do i=1,n
321 +
         g1(i)=teta*g1(i)+teta1*dir(i)
322 +
      end do
323 +
      cof1=-.2/sig1
324 +
277   continue
325 +
c
326 +
c      premiere mise a jour de dh
327 +
      sig=1.0d+0/cof1
328 +
      zsig1=1.0d+0/sig1
329 +
      mk=0
330 +
      ir=nfac
331 +
      epsmc=1.0d-9
332 +
      call calmaj(dh,n,g1,sig,x2,ir,mk,epsmc,nfac)
333 +
      if(ir.ne.nfac)go to 280
334 +
      call calmaj(dh,n,dir,sig1,x2,ir,mk,epsmc,nfac)
335 +
      if(ir.ne.nfac)go to 280
336 +
      go to 300
337 +
c$$$280   if(iprint.gt.0) then
338 +
c$$$         write(bufstr,282)
339 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
340 +
c$$$         endif
341 +
c$$$282   format(' qnbd : pb dans appel majour')
342 +
280   indqn=8
343 +
      if(iter.eq.1)indqn=-5
344 +
c$$$      if(iprint.gt.0) then
345 +
c$$$        write(bufstr,123)indqn
346 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
347 +
c$$$  endif
348 +
      return
349 +
 300  continue
350 +
c     
351 +
c     section 3 determination des variables libres et bloquees
352 +
c     
353 +
c     calcul eps1
354 +
c     
355 +
      scal1=scal
356 +
      if(ieps1.eq.1)scal1=0.0d+0
357 +
      if(ieps1.eq.2)scal1=scal*cscal1
358 +
      do i=1,n
359 +
         x1(i)=x(i)-scal1*abs(g(i))*g(i)
360 +
      end do
361 +
      call proj(n,binf,bsup,x1)
362 +
      eps1=0.0d+0
363 +
      do i=1,n
364 +
         eps1=eps1 + abs(x1(i)-x(i))
365 +
      end do
366 +
      eps1=min(eps0,eps1)
367 +
      if(ieps1.eq.1)eps1=0.0d+0
368 +
      if(ieps1.eq.2)eps1=eps1*1.0d+4
369 +
      ifac=0
370 +
      idfac=0
371 +
      k=0
372 +
c     
373 +
c     
374 +
      gr=0.0d+0
375 +
      if(ig.eq.1)gr=0.2*difg/n
376 +
      n3=n
377 +
      if(in.eq.1)n3=n/10
378 +
c     si irit=1 on peut relacher des variables
379 +
      irit=0
380 +
      if(difg1.le.epsrel*difg0)irit=1
381 +
      if(irel.eq.0.or.iter.eq.1)irit=1
382 +
c     
383 +
      tiers=1.0d+0/3.0d+0
384 +
      do k=1,n
385 +
         izig(k)=izig(k)-1
386 +
         if(izig(k).le.0)izig(k)=0
387 +
         bi=binf(k)
388 +
         bs=bsup(k)
389 +
         ic=indic(k)
390 +
         d1=x(k)-bi
391 +
         d2=bs-x(k)
392 +
         dd=(bs-bi)*tiers
393 +
         ep=min(eps1,dd)
394 +
         if(d1.gt.ep)go to 324
395 +
         if(g(k).gt.0.)go to 330
396 +
         go to 335
397 +
 324     if(d2.gt.ep)go to 335
398 +
         if(g(k).gt.0.)go to 335
399 +
         go to 330
400 +
c     on defactorise si necessaire
401 +
 330     continue
402 +
         if(ic.gt.nfac)go to 340
403 +
         idfac=idfac+1
404 +
         mode=-1
405 +
         izig(k)=izig(k) + izag
406 +
         call ajour(mode,n,k,nfac,dh,x2,indic)
407 +
         if(mode.eq.0) go to 340
408 +
         indqn=8
409 +
         if(iter.eq.1)indqn=-5
410 +
         return
411 +
c     on factorise
412 +
 335     continue
413 +
         if(irit.eq.0)go to 340
414 +
         if(ic.le.nfac)go to 340
415 +
         if(izig(k).ge.1)go to 340
416 +
         mode=1
417 +
         if(ifac.ge.n3.and.iter.gt.1)go to 340
418 +
         if(abs(g(k)).le.gr)go to 340
419 +
         ifac=ifac+1
420 +
         call ajour(mode,n,k,nfac,dh,x2,indic)
421 +
         if(mode.eq.0)go to 340
422 +
         indqn=8
423 +
         if(iter.eq.1)indqn=-5
424 +
         return
425 +
      end do
426 +
 340  continue
427 +
c
428 +
c     *********************************************** a voir
429 +
      if(iconv.eq.1)return
430 +
c
431 +
c     section 6 resolution systeme lineaire et expression de dir
432 +
c     on inverse le syst correspondant aux nl premieres composantes
433 +
c     dans le nouveau syst d indices
434 +
      ir=nfac
435 +
      do i=1,n
436 +
         i1=indic(i)
437 +
         x2(i1)=g(i)
438 +
      end do
439 +
      if(ir.lt.nfac) go to 412
440 +
      if(nfac.gt.1) go to 400
441 +
      x2(1)=x2(1)/dh(1)
442 +
      go to 412
443 +
400   continue
444 +
      do i=2,nfac
445 +
         ij=i
446 +
         i1=i-1
447 +
         v=x2(i)
448 +
         do j=1,i1
449 +
            v=v-dh(ij)*x2(j)
450 +
            ij=ij+nfac-j
451 +
         end do
452 +
         x2(i)=v
453 +
         x2(i)=v
454 +
      end do
455 +
      x2(nfac)=x2(nfac)/dh(ij)
456 +
      np=nfac+1
457 +
      do nip=2,nfac
458 +
         i=np-nip
459 +
         ii=ij-nip
460 +
         v=x2(i)/dh(ii)
461 +
         ip=i+1
462 +
         ij=ii
463 +
         do j=ip,nfac
464 +
            ii=ii+1
465 +
            v=v-dh(ii)*x2(j)
466 +
         end do
467 +
         x2(i)=v
468 +
      end do
469 +
412   continue
470 +
      if(ir.eq.nfac)go to 660
471 +
      indqn=7
472 +
      if(iter.eq.1)indqn=-6
473 +
      return
474 +
660   continue
475 +
      do i=1,n
476 +
         i1=indic(i)
477 +
         dir(i)=-g(i)*scal
478 +
         if(i1.le.nfac) dir(i)=-x2(i1)
479 +
      end do
480 +
      continue
481 +
c
482 +
c     gestion contraintes actives (si iact=1)
483 +
      if(iact.ne.1)go to 675
484 +
      do i=1,n
485 +
         if(izig(i).gt.0)dir(i)=0.
486 +
         if(indic(i).gt.nfac)dir(i)=0.0d+0
487 +
      end do
488 +
675   continue
489 +
c
490 +
c     recherche lineaire
491 +
c     conservation de x et g . calcul de dir+ et fpn
492 +
      do i=1,n
493 +
         g1(i)=g(i)
494 +
         x1(i)=x(i)
495 +
      end do
496 +
c     ifp =1 si fpn trop petit. on prend alors d=-g
497 +
      ifp=0
498 +
      fn=f
499 +
709   fpn=0.0d+0
500 +
      do i=1,n
501 +
         if(x(i)-binf(i).le.epsx(i).and.dir(i).lt.0.)dir(i)=0.0d+0
502 +
         if(bsup(i)-x(i).le.epsx(i).and.dir(i).gt.0.)dir(i)=0.0d+0
503 +
         fpn=fpn + g(i)*dir(i)
504 +
      end do
505 +
      if(fpn.gt.0.0d+0) then
506 +
         if(ifp.eq.1) then
507 +
c$$$            if(iprint.gt.0) then
508 +
c$$$              write(bufstr,1705) fpn
509 +
c$$$              call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
510 +
c$$$              endif
511 +
c$$$1705        format(' qnbd : arret fpn non negatif=',d11.4)
512 +
            indqn=6
513 +
            if(iter.eq.1)indqn=-3
514 +
c$$$            if(iprint.gt.0) then
515 +
c$$$              write(bufstr,123)indqn
516 +
c$$$              call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
517 +
c$$$              endif
518 +
            return
519 +
         else
520 +
            ifp=1
521 +
            do i=1,n
522 +
               if(izig(i).gt.0)dir(i)=-scal*g(i)
523 +
            end do
524 +
            irit=1
525 +
            go to 709
526 +
         endif
527 +
      endif
528 +
c     calcul du t initial suivant une idee de fletcher
529 +
      t1=t
530 +
      if(iter.eq.1) diff=df0
531 +
      t=-2.0d+0*diff/fpn
532 +
      if(t.gt.0.30d+0.and.t.lt.3.0d+0)t=1.0d+0
533 +
      if(eps1.lt.eps0)t=1.0d+0
534 +
      if(indqn.eq.2)t=1.0d+0
535 +
      if(iter.gt.1.and.t1.gt.0.010d+0.and.t1.lt.100.0d+0)t=1.0d+0
536 +
      tmax=1.0d+10
537 +
      t=min(t,tmax)
538 +
      t=max(t,10.*zero)
539 +
c     amd,amf tests sur h'(t) et diff
540 +
      amd=.7
541 +
      amf=.1
542 +
      napm=15
543 +
      napm1=nap + napm
544 +
      if(napm1.gt.napmax)napm1=napmax
545 +
      call rlbd(indrl,n,simul,x,binf,bsup,fn,fpn,t,tmax,dir,g,
546 +
     &     tproj,amd,amf,zero,nap,napm1,x2,izs,rzs,dzs)
547 +
      if(indrl.ge.10)then
548 +
         indsim=4
549 +
         nap=nap + 1
550 +
         call simul(indsim,n,x,f,g,izs,rzs,dzs)
551 +
         if(indsim.le.0)then
552 +
            indqn=-3
553 +
            if(indsim.eq.0)indqn=0
554 +
c$$$            if(iprint.gt.0) then
555 +
c$$$              write(bufstr,123)indqn
556 +
c$$$              call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
557 +
c$$$              endif
558 +
            return
559 +
         endif
560 +
      endif
561 +
      if(indrl.le.0)then
562 +
         indqn=10
563 +
         if(indrl.eq.0)indqn=0
564 +
         if(indrl.eq.-3)indqn=13
565 +
         if(indrl.eq.-4)indqn=12
566 +
         if(indrl.le.-1000)indqn=11
567 +
c$$$         if(iprint.gt.0) then
568 +
c$$$           write(bufstr,123)indqn
569 +
c$$$           call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
570 +
c$$$           endif
571 +
         return
572 +
      endif
573 +
c
574 +
      if(nap.lt.napmax)go to 758
575 +
      f=fn
576 +
c$$$      if(iprint.gt.0) then
577 +
c$$$        write(bufstr,755)napmax
578 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
579 +
c$$$        endif
580 +
c$$$755   format(' qnbd : retour cause max appels simul',i9)
581 +
      indqn=4
582 +
c$$$      if(iprint.gt.0) then
583 +
c$$$        write(bufstr,123)indqn
584 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
585 +
c$$$        endif
586 +
      return
587 +
758   continue
588 +
c     section 8 test de convergence
589 +
c
590 +
      do i=1,n
591 +
         if(abs(x(i)-x1(i)).gt.epsx(i))go to 806
592 +
      end do
593 +
      f=fn
594 +
c$$$      if(iprint.gt.0) then
595 +
c$$$        write(bufstr,1805)
596 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
597 +
c$$$        endif
598 +
c$$$1805  format(' qnbd : retour apres convergence de x')
599 +
      indqn=3
600 +
c$$$      if(iprint.gt.0) then
601 +
c$$$        write(bufstr,123)indqn
602 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
603 +
c$$$        endif
604 +
      return
605 +
806   continue
606 +
      difg=0.0d+0
607 +
      do i=1,n
608 +
         aa=g(i)
609 +
         if(x(i)-binf(i).le.epsx(i))aa=min(0.0d+0,aa)
610 +
         if(bsup(i)-x(i).le.epsx(i))aa=max(0.0d+0,aa)
611 +
         difg=difg + aa**2
612 +
      end do
613 +
      difg1=0.0d+0
614 +
      do i=1,n
615 +
         if(indic(i).gt.nfac)go to 820
616 +
         aa=g(i)
617 +
         if(x(i)-binf(i).le.epsx(i))aa=min(0.0d+0,aa)
618 +
         if(bsup(i)-x(i).le.epsx(i))aa=max(0.0d+0,aa)
619 +
         difg1=difg1 + aa**2
620 +
      end do
621 +
820   continue
622 +
      difg1=sqrt(difg1)
623 +
      difg=sqrt(difg)
624 +
      difg=difg/sqrt(real(n))
625 +
      diff=abs(f-fn)
626 +
      df0=-diff
627 +
      if(irit.eq.1)difg0=difg1
628 +
      f=fn
629 +
c$$$      if(iprint.ge.2) then
630 +
c$$$         write(bufstr,860)epsg,difg,epsf,diff,nap
631 +
c$$$         call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
632 +
c$$$         endif
633 +
c$$$860   format(' qnbd : epsg,difg=',2d11.4,'  epsf,diff=',2d11.4
634 +
c$$$     &,'  nap=',i3)
635 +
      if(diff.lt.epsf)then
636 +
      indqn=2
637 +
c$$$      if(iprint.gt.0) then
638 +
c$$$        write(bufstr,1865)diff
639 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
640 +
c$$$        endif
641 +
c$$$1865  format(' qnbd : retour cause decroissance f trop petite=',d11.4)
642 +
c$$$      if(iprint.gt.0) then
643 +
c$$$        write(bufstr,123)indqn
644 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
645 +
c$$$        endif
646 +
      return
647 +
      endif
648 +
      if(difg.gt.epsg)go to 200
649 +
      indqn=1
650 +
c$$$      if(iprint.gt.0) then
651 +
c$$$        write(bufstr,1900)difg
652 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
653 +
c$$$        endif
654 +
c$$$1900  format(' qnbd : retour cause gradient projete petit=',d11.4)
655 +
c$$$      if(iprint.gt.0) then
656 +
c$$$        write(bufstr,123)indqn
657 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
658 +
c$$$        endif
659 +
      return
660 +
      end

@@ -1,6 +1,6 @@
Loading
1 -
##' n1qn1c optimization
1 +
##' n1qn1 optimization
2 2
##'
3 -
##' This is an R port of the nq1n1c optimization procedure in scilab.
3 +
##' This is an R port of the n1qn1 optimization procedure in scilab.
4 4
##'
5 5
##' @param call_eval Objective function
6 6
##' @param call_grad Gradient Function
@@ -35,7 +35,7 @@
Loading
35 35
##' @author C. Lemarechal, Stephen L. Campbell, Jean-Philippe
36 36
##'     Chancelier, Ramine Nikoukhah, Wenping Wang & Matthew L. Fidler
37 37
##' @importFrom Rcpp evalCpp
38 -
##' @useDynLib nq1n1c, .registration=TRUE
38 +
##' @useDynLib n1qn1, .registration=TRUE
39 39
##' @examples
40 40
##'
41 41
##' ## Rosenbrock's banana function
@@ -69,7 +69,7 @@
Loading
69 69
##' nzm=as.integer(n*(n+13L)/2L)
70 70
##' zm=double(nzm)
71 71
##'
72 -
##' (op1 <- nq1n1c(fr, grr, x, imp=3))
72 +
##' (op1 <- n1qn1(fr, grr, x, imp=3))
73 73
##'
74 74
##' ## Note there are 40 function calls and 40 gradient calls in the above optimization
75 75
##'
@@ -82,11 +82,11 @@
Loading
82 82
##'             200.024349)
83 83
##' c.hess <- c(c.hess, rep(0, 24 - length(c.hess)))
84 84
##'
85 -
##' (op2 <- nq1n1c(fr, grr, x,imp=3, zm=c.hess))
85 +
##' (op2 <- n1qn1(fr, grr, x,imp=3, zm=c.hess))
86 86
##'
87 87
##' ## Note with this knowledge, there were only 29 function/gradient calls
88 88
##'
89 -
##' (op3 <- nq1n1c(fr, grr, x, imp=3, zm=op1$c.hess))
89 +
##' (op3 <- n1qn1(fr, grr, x, imp=3, zm=op1$c.hess))
90 90
##'
91 91
##' ## The number of function evaluations is still reduced because the Hessian
92 92
##' ## is closer to what it should be than the initial guess.
@@ -95,47 +95,45 @@
Loading
95 95
##' ## Optimization time.
96 96
##'
97 97
##' @export
98 -
n1qn1c <- function(call_eval, call_grad, vars, environment = parent.frame(1), ...,
99 -
                   epsilon = .Machine$double.eps, max_iterations = 100, nsim = 100,
100 -
                   imp = 0,
101 -
                   invisible = NULL,
102 -
                   zm = NULL, restart = FALSE,
103 -
                   assign = FALSE,
104 -
                   print.functions = FALSE) {
105 -
  if (!is.null(invisible)) {
106 -
    if (invisible == 1) {
107 -
      imp <- 0
108 -
      print.functions <- FALSE
109 -
    } else {
110 -
      print.functions <- TRUE
98 +
n1qn1 <- function(call_eval, call_grad, vars, environment=parent.frame(1), ...,
99 +
                  epsilon=.Machine$double.eps, max_iterations=100, nsim=100,
100 +
                  imp=0,
101 +
                  invisible=NULL,
102 +
                  zm=NULL, restart=FALSE,
103 +
                  assign=FALSE,
104 +
                  print.functions=FALSE){
105 +
    if (!is.null(invisible)){
106 +
        if (invisible == 1){
107 +
            imp <- 0;
108 +
            print.functions <- FALSE
109 +
        } else {
110 +
            print.functions <- TRUE
111 +
        }
112 +
    }
113 +
    if (!missing(max_iterations) && missing(nsim)){
114 +
        nsim <- max_iterations * 10;
111 115
    }
112 -
  }
113 -
  if (!missing(max_iterations) && missing(nsim)) {
114 -
    nsim <- max_iterations * 10
115 -
  }
116 -
  n <- as.integer(length(vars))
117 -
  imp <- as.integer(imp)
118 -
  max_iterations <- as.integer(max_iterations)
119 -
  nsim <- as.integer(nsim)
120 -
  nzm <- as.integer(ceiling(n * (n + 13) / 2))
121 -
  nsim <- as.integer(nsim)
122 -
  epsilon <- as.double(epsilon)
123 -
  if (is.null(zm)) {
124 -
    mode <- 1L
125 -
    zm <- double(nzm)
126 -
  } else {
127 -
    mode <- 2L
128 -
    if (restart) model <- 3L
129 -
    if (length(zm) != nzm) {
130 -
      stop(sprintf("Compressed Hessian not the right length for this problem.  It should be %d.", nzm))
116 +
    n <- as.integer(length(vars));
117 +
    imp <- as.integer(imp);
118 +
    max_iterations <- as.integer(max_iterations)
119 +
    nsim <- as.integer(nsim)
120 +
    nzm <- as.integer(ceiling(n * (n + 13) / 2));
121 +
    nsim <- as.integer(nsim);
122 +
    epsilon <- as.double(epsilon)
123 +
    if (is.null(zm)){
124 +
        mode <- 1L
125 +
        zm <- double(nzm);
126 +
    } else {
127 +
        mode <- 2L
128 +
        if (restart) model <- 3L
129 +
        if (length(zm) != nzm){
130 +
            stop(sprintf("Compressed Hessian not the right length for this problem.  It should be %d.", nzm))
131 +
        }
131 132
    }
132 -
  }
133 -
  ret <- .Call(
134 -
    n1qn1c_wrap, call_eval, call_grad, environment,
135 -
    vars, epsilon, n, mode, max_iterations, nsim, imp, nzm, zm, as.integer(print.functions)
136 -
  )
137 -
  if (assign) environment$c.hess <- ret$hess
138 -
  return(ret)
133 +
    ret <- .Call(n1qn1_wrap, call_eval, call_grad, environment,
134 +
                 vars, epsilon, n, mode, max_iterations, nsim, imp, nzm, zm, as.integer(print.functions));
135 +
    if (assign) environment$c.hess <- ret$hess;
136 +
    return(ret)
139 137
}
140 138
141 139
@@ -156,7 +154,7 @@
Loading
156 154
##' @param epsg Gradient eps for exiting
157 155
##' @param epsx Parameter eps for exiting
158 156
##' @param print.functions
159 -
##' @inheritParams nq1n1c
157 +
##' @inheritParams n1qn1
160 158
##' @export
161 159
##' @examples
162 160
##'
@@ -190,16 +188,14 @@
Loading
190 188
##' op1 <- qnbd(x, fr, grr)
191 189
##'
192 190
##' @export
193 -
qnbd <- function(par, fn, gr, lower = -Inf, upper = Inf, environment = parent.frame(1),
194 -
                 zero = sqrt(.Machine$double.eps / 7e-07), maxFn = 10000L, maxIt = 10000L,
195 -
                 epsf = sqrt(.Machine$double.eps), epsg = sqrt(.Machine$double.eps),
196 -
                 epsx = sqrt(.Machine$double.eps), print.functions = FALSE) {
197 -
  n <- length(par)
198 -
  if (length(lower) == 1) lower <- rep(lower, n)
199 -
  if (length(upper) == 1) upper <- rep(upper, n)
200 -
  ret <- .Call(
201 -
    qnbd_wrap, fn, gr, environment, par, lower, upper, zero, as.integer(maxFn),
202 -
    as.integer(maxIt), epsf, epsg, epsx, as.integer(n), as.integer(print.functions)
203 -
  )
204 -
  return(ret)
191 +
qnbd <- function(par, fn, gr, lower= -Inf, upper=Inf, environment=parent.frame(1),
192 +
                 zero=sqrt(.Machine$double.eps/7e-07), maxFn=10000L, maxIt=10000L,
193 +
                 epsf=sqrt(.Machine$double.eps), epsg=sqrt(.Machine$double.eps),
194 +
                 epsx=sqrt(.Machine$double.eps), print.functions=FALSE){
195 +
    n <- length(par);
196 +
    if (length(lower) == 1) lower <- rep(lower, n);
197 +
    if (length(upper) == 1) upper <- rep(upper, n);
198 +
    ret <- .Call(qnbd_wrap, fn, gr, environment, par, lower, upper, zero, as.integer(maxFn),
199 +
                           as.integer(maxIt), epsf, epsg, epsx, as.integer(n), as.integer(print.functions));
200 +
    return(ret);
205 201
}

@@ -0,0 +1,535 @@
Loading
1 +
c     Modified by Matthew Fidler in 2017 for different outputs to the R console
2 +
      integer function vff(n, g)
3 +
      integer n, i, ret
4 +
      double precision g, x
5 +
      dimension g(n)
6 +
      ret = 0
7 +
      do i=1, n
8 +
         if (abs(g(i)) .ge. huge(x)) then
9 +
            ret = 1
10 +
            goto 7710
11 +
         endif
12 +
      end do 
13 +
 7710    continue
14 +
      vff = ret
15 +
      end
16 +
c     Note this should be thread safe since it doesn't use global variables
17 +
c     Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
18 +
c Copyright (C) 1987 - INRIA - Claude LEMARECHAL
19 +
c 
20 +
c This file must be used under the terms of the CeCILL.
21 +
c This source file is licensed as described in the file COPYING, which
22 +
c you should have received as part of this distribution.  The terms
23 +
c are also available at    
24 +
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
25 +
c
26 +
      subroutine n1qn1 (simul,n,x,f,g,var,eps,
27 +
     1     mode,niter,nsim,imp,zm,izs,rzs,dzs)
28 +
c
29 +
c!but
30 +
c     minimisation d une fonction reguliere sans contraintes
31 +
c!origine
32 +
c     c. lemarechal, inria, 1987
33 +
c     Copyright INRIA
34 +
c!methode
35 +
c     direction de descente calculee par une methode de quasi-newton
36 +
c     recherche lineaire de type wolfe
37 +
c!liste d appel
38 +
c     simul    : point d'entree au module de simulation (cf normes modulopt i)
39 +
c     n1qn1 appelle toujours simul avec indic = 4 ; le module de
40 +
c     simulation doit se presenter sous la forme subroutine simul
41 +
c     (n,x, f, g, izs, rzs, dzs) et e^tre declare en external dans le
42 +
c     programme appelant n1qn1.
43 +
c     n (e)    : nombre de variables dont depend f.
44 +
c     x (e-s)   : vecteur de dimension n ; en entree le point initial ;
45 +
c                 en sortie : le point final calcule par n1qn1.
46 +
c     f (e-s)   : scalaire ; en entree valeur de f en x (initial), en sortie
47 +
c                 valeur de f en x (final).
48 +
c     g (e-s)   : vecteur de dimension n : en entree valeur du gradient en x
49 +
c                 (initial), en sortie valeur du gradient en x (final).
50 +
c     var (e)   : vecteur strictement positif de dimension n. amplitude de la
51 +
c                 modif souhaitee a la premiere iteration sur x(i).une bonne
52 +
c                 valeur est 10% de la difference (en valeur absolue) avec la
53 +
c                 coordonee x(i) optimale
54 +
c     eps (e-s) : en entree scalaire definit la precision du test d'arret.
55 +
c      le programme considere que la convergence est obtenue lorque il lui
56 +
c      est impossible de diminuer f en attribuant a au moins une coordonnee
57 +
c      x(i) une variation superieure a eps*var(i).
58 +
c      en sortie, eps contient le carre de la norme du gradient en x (final).
59 +
c     mode (e)     : definit l approximation initiale du hessien
60 +
c                  =1 n1qn1 l initialise lui-meme
61 +
c                  =2 le hessien est fourni dans zm sous forme compressee (zm
62 +
c                     contient les colonnes de la partie inferieure du hessien)
63 +
c     niter (e-s)  : en entree nombre maximal d'iterations : en sortie nombre
64 +
c                    d'iterations reellement effectuees.
65 +
c     nsim (e-s)  : en entree nombre maximal d'appels a simul (c'est a dire
66 +
c         avec indic = 4). en sortie le nombre de tels appels reellement faits.
67 +
c      imp (e)   : contro^le les messages d'impression :
68 +
c                  0 rien n'est imprime
69 +
c                  = 1 impressions initiales et finales
70 +
c                  = 2 une impression par iteration (nombre d'iterations,
71 +
c                      nombre d'appels a simul, valeur courante de f).
72 +
c                  >=3 informations supplementaires sur les recherches
73 +
c                      lineaires ;
74 +
c                      tres utile pour detecter les erreurs dans le gradient.
75 +
c      lp (e)    : le numero du canal de sortie, i.e. les impressions
76 +
c                  commandees par imp sont faites par write (lp, format).
77 +
c     zm     : memoire de travail pour n1qn1 de   dimension n*(n+13)/2.
78 +
c     izs,rzs,dzs memoires reservees au simulateur (cf doc)
79 +
c
80 +
c!
81 +
      implicit double precision (a-h,o-z)
82 +
      dimension x(n),g(n),var(n),zm(*),izs(*),dzs(*)
83 +
      real rzs(*)
84 +
      external simul
85 +
      nd=1+(n*(n+1))/2
86 +
      nw=nd+n
87 +
      nxa=nw+n
88 +
      nga=nxa+n
89 +
      nxb=nga+n
90 +
      ngb=nxb+n
91 +
      call n1qn1a (simul,n,x,f,g,var,eps,mode,
92 +
     1 niter,nsim,imp,zm,zm(nd),zm(nw),zm(nxa),zm(nga),
93 +
     2 zm(nxb),zm(ngb),izs,rzs,dzs)
94 +
      end
95 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
96 +
c Copyright (C) INRIA
97 +
c
98 +
c This file must be used under the terms of the CeCILL.
99 +
c This source file is licensed as described in the file COPYING, which
100 +
c you should have received as part of this distribution.  The terms
101 +
c are also available at
102 +
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
103 +
c
104 +
      subroutine n1qn1a (simul,n,x,f,g,scale,acc,mode,
105 +
     1     niter,nsim,iprint,h,d,w,xa,ga,xb,gb,izs,rzs,dzs)
106 +
c
107 +
108 +
*     A (very) few modifs by Bruno (14 March 2005): I have translated some output
109 +
*     information in english (but they don't use format instruction
110 +
*     which is put in the second arg of write). Also for the linear
111 +
*     search output information I divide by the direction vector norm
112 +
*     to get the "normalized" directionnal derivative. Note that this is
113 +
*     just for output (the computing code is normally not modified).
114 +
115 +
      implicit double precision (a-h,o-z)
116 +
      dimension x(n),g(n),scale(n),h(*),d(n),w(n),
117 +
     1 xa(n),ga(n),xb(n),gb(n),izs(*),dzs(*)
118 +
      real rzs(*)
119 +
      external simul
120 +
      ! (blas routine) added by Bruno to get
121 +
      ! a better information concerning directionnal derivative
122 +
      integer vff
123 +
c
124 +
c              calcul initial de fonction-gradient
125 +
c
126 +
      indic=4
127 +
      call simul (indic,n,x,f,g,izs,rzs,dzs)
128 +
c     next line added by Serge to avoid Inf and Nan's (04/2007)
129 +
      if ( abs(f) .ge. huge(f) .and. vff(n, g).ne.1) indic=-1
130 +
      if (indic.gt.0) go to 13
131 +
      acc=0.0d+0
132 +
      niter=1
133 +
      nsim=1
134 +
      return
135 +
   13 nfun=1
136 +
      iecri=0
137 +
      itr=0
138 +
      np=n+1
139 +
c                  initialisation du hessien, en fonction de var
140 +
      if (mode.ge.2) go to 60
141 +
   20 c=0.0d+0
142 +
      do i=1,n
143 +
         c=max(c,abs(g(i)*scale(i)))
144 +
      end do
145 +
      if (c.le.0.0d+0) c=1.0d+0
146 +
      k=(n*np)/2
147 +
      do i=1,k
148 +
         h(i)=0.0d+0
149 +
      end do
150 +
      k=1
151 +
      do i=1,n
152 +
         h(k)=0.010d+0*c/(scale(i)*scale(i))
153 +
         k=k+np-i
154 +
      end do
155 +
      go to 100
156 +
c               factorisation du hessien
157 +
   60 if (mode.ge.3) go to 80
158 +
      k=n
159 +
      if(n.gt.1) go to 300
160 +
      if(h(1).gt.0.0d+0) go to 305
161 +
      h(1)=0.0d+0
162 +
      k=0
163 +
      go to 305
164 +
  300 continue
165 +
      np=n+1
166 +
      ii=1
167 +
      do i=2,n
168 +
         hh=h(ii)
169 +
         ni=ii+np-i
170 +
         if(hh.gt.0.0d+0) go to 301
171 +
         h(ii)=0.0d+0
172 +
         k=k-1
173 +
         ii=ni+1
174 +
         go to 304
175 +
 301     continue
176 +
         ip=ii+1
177 +
         ii=ni+1
178 +
         jk=ii
179 +
         do ij=ip,ni
180 +
            v=h(ij)/hh
181 +
            do ik=ij,ni
182 +
               h(jk)=h(jk)-h(ik)*v
183 +
               jk=jk+1
184 +
            end do
185 +
            h(ij)=v
186 +
         end do
187 +
      end do
188 +
  304 continue
189 +
      if(h(ii).gt.0.0d+0) go to 305
190 +
      h(ii)=0.0d+0
191 +
      k=k-1
192 +
 305  continue
193 +
c
194 +
      if (k.ge.n) go to 100
195 +
 70   go to 20
196 +
c                verification que la diagonale est positive
197 +
   80 k=1
198 +
      do i=1,n
199 +
         if (h(k).le.0.0d+0) go to 70
200 +
         k=k+np-i
201 +
      end do
202 +
c                quelques initialisations
203 +
  100 dff=0.0d+0
204 +
  110 fa=f
205 +
      isfv=1
206 +
      do i=1,n
207 +
         xa(i)=x(i)
208 +
         ga(i)=g(i)
209 +
      end do
210 +
c                   iteration
211 +
  130 itr=itr+1
212 +
      ial=0
213 +
      if (itr.gt.niter) go to 250
214 +
      iecri=iecri+1
215 +
      if (iecri.ne.-iprint) go to 140
216 +
      iecri=0
217 +
      indic=1
218 +
      call simul(indic,n,x,f,g,izs,rzs,dzs)
219 +
c     error in user function
220 +
      if(indic.eq.0) goto 250
221 +
c               calcul de la direction de recherche
222 +
  140 do i=1,n
223 +
         d(i)=-ga(i)
224 +
      end do
225 +
      w(1)=d(1)
226 +
      if(n.gt.1)go to 400
227 +
      d(1)=d(1)/h(1)
228 +
      go to 412
229 +
  400 continue
230 +
      do i=2,n
231 +
         ij=i
232 +
         i1=i-1
233 +
         v=d(i)
234 +
         do j=1,i1
235 +
            v=v-h(ij)*d(j)
236 +
            ij=ij+n-j
237 +
         end do
238 +
         w(i)=v
239 +
         d(i)=v
240 +
      end do
241 +
      d(n)=d(n)/h(ij)
242 +
      np=n+1
243 +
      do nip=2,n
244 +
         i=np-nip
245 +
         ii=ij-nip
246 +
         v=d(i)/h(ii)
247 +
         ip=i+1
248 +
         ij=ii
249 +
         do j=ip,n
250 +
            ii=ii+1
251 +
            v=v-h(ii)*d(j)
252 +
         end do
253 +
         d(i)=v
254 +
      end do 
255 +
  412 continue
256 +
c               calcul du pas minimum
257 +
c               et de la derivee directionnelle initiale
258 +
      c=0.0d+0
259 +
      dga=0.0d+0
260 +
      do i=1,n
261 +
         c=max(c,abs(d(i)/scale(i)))
262 +
         dga=dga+ga(i)*d(i)
263 +
      end do
264 +
c               test si la direction est de descente
265 +
      if (dga.ge.0.0d+0) go to 240
266 +
c               initialisation du pas
267 +
      stmin=0.0d+0
268 +
      stepbd=0.0d+0
269 +
      steplb=acc/c
270 +
      fmin=fa
271 +
      gmin=dga
272 +
      step=1.0d+0
273 +
      if (dff.le.0.0d+0) step=min(step,1.0d+0/c)
274 +
      if (dff.gt.0.0d+0) step=min(step,(dff+dff)/(-dga))
275 +
276 +
  170 c=stmin+step
277 +
      if (nfun.ge.nsim) go to 250
278 +
      nfun=nfun+1
279 +
c              calcul de fonction-gradient
280 +
      do i=1,n
281 +
         xb(i)=xa(i)+c*d(i)
282 +
      end do
283 +
      indic=4
284 +
      call simul (indic,n,xb,fb,gb,izs,rzs,dzs)
285 +
c     next line added by Serge to avoid Inf and Nan's (04/2007)
286 +
      if ( abs(fb) .ge. huge(fb) .and. vff(n,gb).ne.1) indic=-1
287 +
c              test sur indic
288 +
      if (indic.gt.0) goto 185
289 +
      if (indic.lt.0) goto 183
290 +
      do i=1,n
291 +
         x(i)=xb(i)
292 +
         g(i)=gb(i)
293 +
      end do
294 +
      go to 250
295 +
  183 stepbd=step
296 +
      ial=1
297 +
      step=step/10.0d+0
298 +
      if (stepbd.gt.steplb) goto 170
299 +
      goto 240
300 +
c             stockage si c'est la plus petite valeur
301 +
  185 isfv=min(2,isfv)
302 +
      if (fb.gt.f) go to 220
303 +
      if (fb.lt.f) go to 200
304 +
      gl1=0.0d+0
305 +
      gl2=0.0d+0
306 +
      do i=1,n
307 +
         gl1=gl1+(scale(i)*g(i))**2
308 +
         gl2=gl2+(scale(i)*gb(i))**2
309 +
      end do
310 +
      if (gl2.ge.gl1) go to 220
311 +
  200 isfv=3
312 +
      f=fb
313 +
      do i=1,n
314 +
         x(i)=xb(i)
315 +
         g(i)=gb(i)
316 +
      end do
317 +
c               calcul de la derivee directionnelle
318 +
  220 dgb=0.0d+0
319 +
      do i=1,n
320 +
         dgb=dgb+gb(i)*d(i)
321 +
      end do
322 +
      if (iprint.lt.3) goto 231
323 +
      s=fb-fa
324 +
  231 if (fb-fa.le.0.10d+0*c*dga) go to 280
325 +
      ial=0
326 +
c               iteration terminee si le pas est minimum
327 +
      if (step.gt.steplb) go to 270
328 +
  240 if (isfv.ge.2) go to 110
329 +
c               ici, tout est termine
330 +
 250  acc=0.0d+0
331 +
      do i=1,n
332 +
         acc=acc+g(i)*g(i)
333 +
      end do
334 +
      niter=itr
335 +
      nsim=nfun
336 +
      return
337 +
c               interpolation cubique
338 +
  270 stepbd=step
339 +
      c=gmin+dgb-3.0d+0*(fb-fmin)/step
340 +
      if(c.eq.0.0d+0) goto 250
341 +
      cc=abs(c)-gmin*(dgb/abs(c))
342 +
      cc=sqrt(abs(c))*sqrt(max(0.0d+0,cc))
343 +
      c=(c-gmin+cc)/(dgb-gmin+cc+cc)
344 +
      step=step*max(0.10d+0,c)
345 +
      go to 170
346 +
c               ceci est un pas de descente
347 +
  280 if (ial.eq.0) goto 285
348 +
      if (stepbd.gt.steplb) go to 285
349 +
      go to 240
350 +
  285 stepbd=stepbd-step
351 +
      stmin=c
352 +
      fmin=fb
353 +
      gmin=dgb
354 +
c               extrapolation
355 +
      step=9.0d+0*stmin
356 +
      if (stepbd.gt.0.0d+0) step=0.50d+0*stepbd
357 +
      c=dga+3.0d+0*dgb-4.0d+0*(fb-fa)/stmin
358 +
      if (c.gt.0.0d+0) step=min(step,stmin*max(1.0d+0,-dgb/c))
359 +
      if (dgb.lt.0.70d+0*dga) go to 170
360 +
c                 recherche lineaire terminee, test de convergence
361 +
      isfv=4-isfv
362 +
      if (stmin+step.le.steplb) go to 240
363 +
c                 formule de bfgs
364 +
      ir=-n
365 +
      do i=1,n
366 +
         xa(i)=xb(i)
367 +
         xb(i)=ga(i)
368 +
         d(i)=gb(i)-ga(i)
369 +
         ga(i)=gb(i)
370 +
      end do
371 +
      call majour(h,xb,w,n,1.0d+0/dga,ir,1,0.0d+0)
372 +
      ir=-ir
373 +
      call majour(h,d,d,n,1.0d+0/(stmin*(dgb-dga)),ir,1,0.0d+0)
374 +
c ww edits
375 +
c     write(*,*) (h(kk), kk=1,(n*(n+1))/2)
376 +
c                  test du rang de la nouvelle matrice
377 +
      if (ir.lt.n) go to 250
378 +
c                  nouvelle iteration
379 +
      dff=fa-fb
380 +
      fa=fb
381 +
      go to 130
382 +
      end
383 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
384 +
c Copyright (C) INRIA
385 +
c 
386 +
c This file must be used under the terms of the CeCILL.
387 +
c This source file is licensed as described in the file COPYING, which
388 +
c you should have received as part of this distribution.  The terms
389 +
c are also available at    
390 +
c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
391 +
c
392 +
      subroutine majour(hm,hd,dd,n,hno,ir,indic,eps)
393 +
c
394 +
      implicit double precision (a-h,o-z)
395 +
      dimension hm(*),hd(n),dd(n)
396 +
      if(n.eq.1) go to 100
397 +
c
398 +
      np=n+1
399 +
      if(hno.gt.0.0d+0) go to 99
400 +
c
401 +
      if(hno.eq.0.0d+0) go to 999
402 +
      if(ir.eq.0) go to 999
403 +
      hon=1.0d+0/hno
404 +
      ll=1
405 +
      if(indic.eq.0) go to 1
406 +
c
407 +
      do i=1,n
408 +
         if(hm(ll).eq.0.0d+0) go to 2 
409 +
         hon=hon+dd(i)**2/hm(ll)
410 +
         ll=ll+np-i
411 +
      end do
412 +
 2    continue
413 +
      go to 3
414 +
c
415 +
 1    continue
416 +
      do i=1,n
417 +
         dd(i)=hd(i)
418 +
      end do
419 +
      do i=1,n
420 +
         iplus=i+1
421 +
         del=dd(i)
422 +
         if(hm(ll).gt.0.0d+0) go to 6
423 +
         dd(i)=0.0d+0
424 +
         ll=ll+np-i
425 +
         go to 5
426 +
 6       continue
427 +
         hon=hon+del**2/hm(ll)
428 +
         if(i.eq.n) go to 7
429 +
         do j=iplus,n
430 +
            ll=ll+1
431 +
            dd(j)=dd(j)-del*hm(ll)
432 +
         end do
433 +
 7       ll=ll+1
434 +
      end do
435 +
 5    continue
436 +
c
437 +
 3    continue
438 +
      if(ir.le.0) go to 9
439 +
      if(hon.gt.0.0d+0) go to 10
440 +
      if ((indic-1) .le. 0) then 
441 +
         goto 99
442 +
      else
443 +
         goto 11
444 +
      endif
445 +
 9    continue
446 +
      hon=0.0d+0
447 +
      ir=-ir-1
448 +
      go to 11
449 +
 10   hon=eps/hno
450 +
      if(eps.eq.0.0d+0)ir=ir-1
451 +
 11   continue
452 +
      mm=1
453 +
      honm=hon
454 +
      do  i=1,n
455 +
         j=np-i
456 +
         ll=ll-i
457 +
         if(hm(ll).ne.0.0d+0) honm=hon-dd(j)**2/hm(ll)
458 +
         dd(j)=hon
459 +
         hon=honm
460 +
      end do
461 +
      go to 13
462 +
c
463 +
 99   continue
464 +
      mm=0
465 +
      honm=1.0d+0/hno
466 +
 13   continue
467 +
      ll=1
468 +
c
469 +
      do 98 i=1,n
470 +
         iplus=i+1
471 +
         del=hd(i)
472 +
         if(hm(ll).gt.0.0d+0) go to 14
473 +
         if(ir.gt.0) go to 15
474 +
         if(hno.lt.0.0d+0) go to 15
475 +
         if(del.eq.0.0d+0) go to 15
476 +
         ir=1-ir
477 +
         hm(ll)=del**2/honm
478 +
         if(i.eq.n) go to 999
479 +
         do j=iplus,n
480 +
            ll=ll+1
481 +
            hm(ll)=hd(j)/del
482 +
         end do
483 +
         go to 999
484 +
 15      continue
485 +
         hon=honm
486 +
         ll=ll+np-i
487 +
         go to 98
488 +
 14      continue
489 +
         hml=del/hm(ll)
490 +
         if (mm .le. 0) then 
491 +
            goto 17
492 +
         else
493 +
            goto 18
494 +
         endif
495 +
 17      hon=honm+del*hml
496 +
         go to 19
497 +
 18      hon=dd(i)
498 +
 19      continue
499 +
         r=hon/honm
500 +
         hm(ll)=hm(ll)*r
501 +
         if(r.eq.0.0d+0) go to 20
502 +
         if(i.eq.n)go to 20
503 +
         b=hml/hon
504 +
         if(r.gt.4.0d+0) go to 21
505 +
         do j=iplus,n
506 +
            ll=ll+1
507 +
            hd(j)=hd(j)-del*hm(ll)
508 +
            hm(ll)=hm(ll)+b*hd(j)
509 +
         end do
510 +
         go to 23
511 +
 21      gm=honm/hon
512 +
         do j=iplus,n
513 +
            ll=ll+1
514 +
            y=hm(ll)
515 +
            hm(ll)=b*hd(j)+y*gm
516 +
            hd(j)=hd(j)-del*y
517 +
         end do
518 +
 23     continue
519 +
        honm=hon
520 +
        ll=ll+1
521 +
 98   continue
522 +
c
523 +
 20   continue
524 +
      if(ir.lt.0)ir=-ir
525 +
      go to 999
526 +
 100  continue
527 +
      hm(1)=hm(1)+hno *hd(1)**2
528 +
      ir=1
529 +
      if(hm(1).gt.0.0d+0) go to 999
530 +
      hm(1)=0.0d+0
531 +
      ir=0
532 +
 999  continue
533 +
      return
534 +
      end
535 +
       

@@ -0,0 +1,44 @@
Loading
1 +
c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c Copyright (C) INRIA
3 +
c 
4 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine calmaj(dh,n,g1,sig,w,ir,mk,epsmc,nfac)
14 +
c
15 +
      implicit double precision (a-h,o-z)
16 +
c     subroutine de qnbd
17 +
      dimension dh(*),g1(*),w(*)
18 +
      if(nfac.eq.n) go to 50
19 +
      nfac1=nfac+1
20 +
      nnfac=n-nfac
21 +
      n2fac=(nfac*nfac1)/2
22 +
      do i=1,n
23 +
         w(i)=g1(i)*sig
24 +
      end do
25 +
      k=n2fac
26 +
      if(nfac.eq.0)go to 25
27 +
      do j=1,nfac
28 +
         do i=nfac1,n
29 +
            k=k+1
30 +
            dh(k)=dh(k)+g1(i)*w(j)
31 +
         end do
32 +
      end do
33 +
25    k=n2fac+nfac*nnfac
34 +
      do j=nfac1,n
35 +
         do i=j,n
36 +
            k=k+1
37 +
            dh(k)=dh(k) + g1(i)*w(j)
38 +
         end do
39 +
      end do
40 +
50    ir=nfac
41 +
      if(nfac.eq.0)return
42 +
      call majour(dh,g1,w,nfac,sig,ir,mk,epsmc)
43 +
      return
44 +
      end

@@ -0,0 +1,527 @@
Loading
1 +
c     Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 +
c Copyright (C) 1988 - INRIA - F. BONNANS
3 +
c 
4 +
c Copyright (C) 2012 - 2016 - Scilab Enterprises
5 +
c
6 +
c This file is hereby licensed under the terms of the GNU GPL v2.0,
7 +
c pursuant to article 5.3.4 of the CeCILL v.2.1.
8 +
c This file was originally licensed under the terms of the CeCILL v2.1,
9 +
c and continues to be available under such terms.
10 +
c For more information, see the COPYING file which you should have received
11 +
c along with this program.
12 +
c
13 +
      subroutine rlbd(indrl,n,simul,x,binf,bsup,f,hp,t,tmax,d,gn,
14 +
     &     tproj,amd,amf,zero,nap,napmax,xn,izs,rzs,dzs)
15 +
c
16 +
c!but
17 +
c     subroutine de recherche lineaire pour des problemes avec
18 +
c       contraintes de borne (traitees par projection)
19 +
c       le critere de retour est une extension de celui de wolfe
20 +
c!methode
21 +
c     pour chaque valeur du parametre t , sont calcules le critere
22 +
c     et son gradient.
23 +
c     une phase d extrapolation permet d obtenir un encadrement.
24 +
c     l intervalle est ensuite reduit suivant les cas par une methode
25 +
c     de dichotomie, d interpolation lineaire sur les derivees ou
26 +
c     d interpolation cubique.
27 +
c
28 +
c!impressions
29 +
c      si iprint > 2 , rlbd fournit les impressions suivantes :
30 +
c
31 +
c      la premiere ligne indique :
32 +
c      t      premiere valeur de t fournie en liste d' appel
33 +
c      tproj  plus petit t > 0 pour lequel on bute sur une borne
34 +
c      dh/dt  derivee en zero de h(t)=f(x+t*d)-f(x)
35 +
c      tmax   valeur maximale de t fournie en liste d' appel
36 +
c
37 +
c      lignes suivantes :
38 +
c      chaine de caracteres en debut de ligne : indique comment sera calcule
39 +
c      le pas de la ligne suivante ;
40 +
c      ic : interpolation cubique
41 +
c      s  : saturation d une variable sur une borne
42 +
c      id : interpolation lineaire sur la derivee
43 +
c      e  : extrapolation
44 +
c      d  :interpolation cubique ayant echouee t est calcule par dichotomie
45 +
c      b  :sauvegarde de convergence active
46 +
c
47 +
c!subroutines utilisees
48 +
c     proj et satur (bibl. modulopt)
49 +
c!liste d appel
50 +
c
51 +
c      subroutine rlbd(indrl,n,simul,proj,x,binf,bsup,f,hp,t,tmax,d,gn,
52 +
c     &  tproj,amd,amf,iprint,io,zero,nap,napmax,xn,izs,rzs,dzs)
53 +
c
54 +
c      e;s;e,s:parametres initialises en entree,en sortie,en entree et
55 +
c              en sortie
56 +
c      indrl<0:la recherche lineaire n a pas trouve de meilleur pas(e,s)
57 +
c           =0:arret demande par l'utilisateur dans simul
58 +
c           >0:meilleur pas fourni avec f et g
59 +
c           >9:meilleur pas fourni avec f et sans g
60 +
c           =14:deltat trop petit
61 +
c           =13:nap=napmax
62 +
c           =8:toutes les variables sont saturees
63 +
c           =4:deltat trop petit
64 +
c           =3:nap=napmax
65 +
c           =2:t=tmax
66 +
c           =1:descente serieuse avec t<tmax
67 +
c           =0:arret demande par l'utilisateur
68 +
c           =-3:nap=napmax
69 +
c           =-4:deltat trop petit
70 +
c           =-1000+indic:nap=napmax et indic<0
71 +
c      n:dimension de x                                    (e,s)
72 +
c      simul: subroutine fournissant le critere et le gradient (e)
73 +
c      x:valeur initiale de la variable a optimiser en entree;valeur a
74 +
c        l optimum en sortie.                                     (e,s)
75 +
c      binf,bsup:bornes inf et sup de dimension n                 (e,s)
76 +
c      f:valeur du critere en x                            (e,s)
77 +
c      hp:derivee de f(x+t*d) par rapport a t en 0         (e)
78 +
c      t:pas                                               (e)
79 +
c      tmax:valeur maximal du pas                          (e,s)
80 +
c      d:direction de descente                             (e)
81 +
c      gn: gradient de f en xn                             (e,s)
82 +
c      tproj:plus petit pas saturant une nouvelle contrainte(e,s)
83 +
c      amf,amd:constantes du test de wolfe                 (e)
84 +
c      iprint<=2:pas d'impression                             (e)
85 +
c         >=3:une impression par calcul de simul           (e)
86 +
c      io:numero du fichier resultat                       (e)
87 +
c      zero:proche du zero machine                         (e)
88 +
c      nap:nombre d'appel a simul                          (e)
89 +
c      napmax:nombre maximum d'appel a simul               (e)
90 +
c      xn:tableau de travail de dimension n (=x+t*d)
91 +
c      izs,rzs,dzs:cf norme modulopt                       (e,s)
92 +
c!
93 +
c       parametres de l algorithme
94 +
c      eps1:sauvegarde de conv.dans l interpolation lineaire sur la derivee
95 +
c      eps:sauvegarde de conv.dans la l interpolation par saturation
96 +
c          d une contrainte.
97 +
c      epst:sauvegerde de conv.dans l interpolation cubique
98 +
c      extra,extrp:servent a calculer la limite sur la variation relative
99 +
c      de t dans l extrapolation et l interpolation lineaire sur la derivee
100 +
c      cofder: intervient dans le choix entre les methodes d' interpolation
101 +
c
102 +
c        variables de travail
103 +
c      fn:valeur du critere en xn
104 +
c      hpn:derivee de f(x+t*d) par rapport a t
105 +
c      hpd:valeur de hpn a droite
106 +
c      hpg:valeur de hpn a gauche
107 +
c      td:pas trop grand
108 +
c      tg:pas trop petit
109 +
c      tproj:plus petit pas saturant une contrainte
110 +
c      tmaxp:plus grand pas saturant une contraite
111 +
c      ftd:valeur de f en x+td*d
112 +
c      ftg:valeur de f en x+tg*d
113 +
c      hptd:valeur de hpn en x+td*d
114 +
c      hptg:valeur de hpn en x+tg*d
115 +
c      imax=1:tmax a ete atteint
116 +
c          =0:tmax n a pas ete atteint
117 +
c      icos:indice de la variable saturee par la borne superieure
118 +
c      icoi:indice de la variable saturee par la borne inferieure
119 +
c      ico1:indice de la variable saturee en tmaxp
120 +
c      icop:indice de la variable saturee en tproj
121 +
c
122 +
      implicit double precision(a-h,o-z)
123 +
      external simul,proj
124 +
      character var2*3
125 +
      dimension x(n),xn(n),gn(n),d(n),binf(n),bsup(n),izs(*)
126 +
      double precision dzs(*)
127 +
      real rzs(*)
128 +
c
129 +
      indrl=1
130 +
      eps1=.90d+0
131 +
      eps=.10d+0
132 +
      epst=.10d+0
133 +
      extrp=100.0d+0
134 +
      extra=10.0d+0
135 +
      cofder=100.
136 +
      var2='   '
137 +
c
138 +
      ta1=0.0d+0
139 +
      f0=f
140 +
      fa1=f
141 +
      hpta1=hp
142 +
      imax=0
143 +
      hptg=hp
144 +
      ftg=f
145 +
      tg=0.0d+0
146 +
      td=0.0d+0
147 +
      icos=0
148 +
      icoi=0
149 +
      icop=0
150 +
c
151 +
c     calcul de tproj:plus petit point de discontinuite de h'(t)
152 +
      tproj=0.0d+0
153 +
      do i=1,n
154 +
         CRES=d(i)
155 +
         if (CRES .lt. 0) then
156 +
            goto 4
157 +
         elseif (CRES .eq. 0) then 
158 +
            goto 7
159 +
         else
160 +
            goto 5
161 +
         endif
162 +
 4       t2=(binf(i)-x(i))/d(i)
163 +
         go to 6
164 +
 5       t2=(bsup(i)-x(i))/d(i)
165 +
 6       if (t2.le.0.0d+0) go to 7
166 +
         if (tproj.eq.0.0d+0) tproj=t2
167 +
         if (t2.gt.tproj) go to 7
168 +
         tproj=t2
169 +
         icop=i
170 +
      end do
171 +
7     continue
172 +
c
173 +
c
174 +
c              boucle
175 +
c
176 +
c     calcul de xn,de fn et de gn
177 +
 200  if (nap.ge.napmax) then
178 +
        k=3
179 +
        goto 1000
180 +
      end if
181 +
      do i=1,n
182 +
         xn(i)=x(i)+t*d(i)
183 +
      end do
184 +
      call proj (n,binf,bsup,xn)
185 +
      if (icos.gt.0) xn(icos)=bsup(icos)
186 +
      if (icoi.gt.0) xn(icoi)=binf(icoi)
187 +
      indic=4
188 +
      call simul (indic,n,xn,fn,gn,izs,rzs,dzs)
189 +
      nap=nap+1
190 +
      if (indic.lt.0) then
191 +
c$$$         if (iprint.ge.3) then
192 +
c$$$           write (bufstr,16000) indic,t
193 +
c$$$           call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
194 +
c$$$           endif
195 +
         if (nap.ge.napmax) go to 1000
196 +
         t=tg+(t-tg)/4.0d+0
197 +
         tmax=t
198 +
         imax=1
199 +
         icoi=0
200 +
         icos=0
201 +
         var2='dd '
202 +
         go to 800
203 +
      endif
204 +
      if(indic.eq.0) then
205 +
         indrl=0
206 +
         go to 1010
207 +
      endif
208 +
c
209 +
c      calcul de hpg et hpd
210 +
      hpg=0.0d+0
211 +
      do i=1,n
212 +
         xn(i)=x(i)+t*d(i)
213 +
      end do
214 +
      if (icoi.gt.0) xn(icoi)=bsup(icoi)
215 +
      if (icos.gt.0) xn(icos)=bsup(icos)
216 +
      call proj(n,binf,bsup,xn)
217 +
      do i=1,n
218 +
         xni=xn(i)
219 +
         if(binf(i).lt.xni.and.xni.lt.bsup(i)) hpg=hpg + d(i)*gn(i)
220 +
      end do
221 +
      hpd=hpg
222 +
      if(icoi.gt.0) hpg=hpg + d(icoi)*gn(icoi)
223 +
      if(icos.gt.0) hpg=hpg + d(icos)*gn(icos)
224 +
c
225 +
      icoi=0
226 +
      icos=0
227 +
      if((hpd.ne.0.0d+0).or.(hpg.ne.0.0d+0)) go to 360
228 +
c
229 +
c      la derivee de h est nulle
230 +
c      calcul du pas saturant toutes les bornes:tmaxp
231 +
      tmaxp=0.0d+0
232 +
      ico1=0
233 +
      do i=1,n
234 +
         CRES=d(i)
235 +
         if (CRES .lt. 0) then 
236 +
            goto 310
237 +
         elseif (CRES .eq. 0) then
238 +
            goto 350
239 +
         else
240 +
            goto 320
241 +
         endif
242 +
 310     t2=(binf(i)-x(i))/d(i)
243 +
         go to 330
244 +
 320     t2=(bsup(i)-x(i))/d(i)
245 +
 330     if (t2.le.0.0d+0) go to 350
246 +
         if (tmaxp.eq.0.0d+0) tmaxp=t2
247 +
         if (tmaxp.gt.t2)go to 350
248 +
         tmaxp=t2
249 +
         ico1=i
250 +
      end do
251 +
350   continue
252 +
      if (t.lt.tmaxp) then
253 +
         if(fn.le.f+amf*hp*t) goto 1010
254 +
         t=t/10.0d+0
255 +
         var2='d  '
256 +
         goto 800
257 +
      end if
258 +
      icos=ico1
259 +
      icoi=0
260 +
      if (d(ico1).lt.0.0d+0) then
261 +
      icoi=ico1
262 +
      icos=0
263 +
      end if
264 +
c
265 +
c     toutes les variables sont saturees
266 +
c$$$      if (iprint.ge.3) then
267 +
c$$$        write (bufstr,3330) tmaxp
268 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
269 +
c$$$        endif
270 +
c$$$3330   format ('toutes les variables sont saturees:tmaxp= ',e11.4)
271 +
      t=tmaxp
272 +
      if (fn.lt.f+amf*hp*tmaxp)  then
273 +
      indrl=8
274 +
      goto 1010
275 +
      end if
276 +
      hpg=d(ico1)*gn(ico1)
277 +
      if ((fn.lt.f).and.(hpg.lt.0.0d+0)) then
278 +
      indrl=8
279 +
      goto 1010
280 +
      end if
281 +
360    continue
282 +
c
283 +
c       test de wolfe
284 +
c
285 +
      a=f+amf*hp*t
286 +
      if (fn.gt.a) then
287 +
c      le pas est trop grand
288 +
c       (dans le cas quadratique changer eps1 et extra si td<tproj)
289 +
            td=t
290 +
            t1=t-ta1
291 +
            h1=(fn-fa1)/t1
292 +
            ftd=fn
293 +
            hptd=hpg
294 +
            ta=tg
295 +
            hpn=hptd
296 +
            hpa=hptg
297 +
            fa=ftg
298 +
      else
299 +
            if (hpd.ge.(amd*hp)) go to 1010
300 +
c      le pas est trop petit
301 +
            tg=t
302 +
            t1=t-ta1
303 +
            h1=(fn-fa1)/t1
304 +
            ftg=fn
305 +
            hptg=hpd
306 +
            ta=td
307 +
            hpn=hptg
308 +
            hpa=hptd
309 +
            fa=ftd
310 +
            if (td.eq.0.0d+0) go to 700
311 +
            a1=abs(hptd/hp)
312 +
            if ((a1.gt.cofder).and.(ftd.gt.f).and.(hptg.gt.(.99*hp)))
313 +
     &       then
314 +
                hpta1=hp
315 +
                fa1=f
316 +
                ta1=0.0d+0
317 +
                go to 700
318 +
             end if
319 +
      endif
320 +
      a1=abs(hpn/hp)
321 +
      if ((tg.ne.0.0d+0).or.(fn.le.f).or.(a1.le.cofder).or.
322 +
     & (hpn.lt.0.0d+0)) then
323 +
        if (td.le.tproj) go to 600
324 +
        go to 500
325 +
      end if
326 +
c
327 +
c       calcul du nouveau t
328 +
c
329 +
c      par interpolation lineaire sur la derivee
330 +
c
331 +
      ta1=t
332 +
      fa1=fn
333 +
      div=hp-hptd
334 +
      text=t/10.0d+0
335 +
      if(abs(div).gt.zero) text=t*(hp/div)
336 +
      if (text.gt.tproj) text=t/10.0d+0
337 +
      text=max(text,t/(extrp*extra))
338 +
      t=min(text,t*eps1)
339 +
      ttsup=1.50d+0*t
340 +
      extrp=10.
341 +
      if (tproj.gt.ta1) then
342 +
          var2='id '
343 +
          go to 800
344 +
      end if
345 +
      ttmin=0.70d+0*t
346 +
      tmi=t
347 +
      topt=0.0d+0
348 +
      iproj=0
349 +
      call satur(n,x,binf,bsup,d,ttmin,ttsup,topt,tg,td,tmi,
350 +
     &            icoi,icos,iproj)
351 +
      var2='id '
352 +
      if (topt.ne.0.0d+0) then
353 +
         t=topt
354 +
         var2='ids'
355 +
      end if
356 +
      go to 800
357 +
c
358 +
c      interpolation par saturation d une contrainte
359 +
c
360 +
500   if (td.le.tproj) go to 600
361 +
      topt=0.0d+0
362 +
      iproj=1
363 +
      ta1=t
364 +
      fa1=fn
365 +
      ttmin=tg+eps*(td-tg)
366 +
      ttsup=td-eps*(td-tg)
367 +
      tmi=(td+tg)/2.0d+0
368 +
      call satur(n,x,binf,bsup,d,ttmin,ttsup,topt,tg,td,tmi,
369 +
     &            icoi,icos,iproj)
370 +
      if (topt.eq.0.0d+0) go to 600
371 +
      t=topt
372 +
      var2='s  '
373 +
      if (t.eq.ttsup.or.t.eq.ttmin) var2='sb '
374 +
      go to 800
375 +
c
376 +
c      interpolation cubique
377 +
c
378 +
c      test de securite
379 +
600   if ((td-tg).lt.1.0d+2*zero) then
380 +
        k=4
381 +
        goto 1000
382 +
      end if
383 +
c
384 +
c      calcul du minimum
385 +
      b=1.0d+0
386 +
      p=hpn+hpa-3.0d+0*(fn-fa)/(t-ta)
387 +
      di=p*p-hpn*hpa
388 +
      if (di.lt.0.0d+0) go to 690
389 +
      if ((t-ta).lt.0.0d+0) b=-1
390 +
      div=hpn+p+b*sqrt(di)
391 +
      if (abs(div).le.zero) go to 690
392 +
      r=hpn/div
393 +
      topt=t-r*(t-ta)
394 +
      if ((topt.lt.tg).or.(topt.gt.td))go to 690
395 +
c
396 +
c      sauvegarde de convergence
397 +
      e=epst*(td-tg)
398 +
      var2='ic '
399 +
      if (topt.gt.(td-e)) then
400 +
        topt=td-e
401 +
        var2='icb'
402 +
      end if
403 +
      if (topt.lt.(tg+e)) then
404 +
         topt=tg+e
405 +
         var2='icb'
406 +
      end if
407 +
      ta1=t
408 +
      fa1=fn
409 +
      t=topt
410 +
      goto 800
411 +
690   ta1=t
412 +
      fa1=fn
413 +
      t=0.50d+0*(tg+td)
414 +
      var2='d  '
415 +
      go to 800
416 +
c
417 +
c      extrapolation
418 +
c
419 +
700   if (imax.ge.1) then
420 +
        k=2
421 +
        goto 1000
422 +
      end if
423 +
      text=10.0d+0*t
424 +
      difhp=hptg-hpta1
425 +
      if (difhp.gt.zero)then
426 +
        text=(amd*hp/3.0d+0-hptg)*((tg-ta1)/difhp)+tg
427 +
        if ((td.ne.0.0d+0).and.(text.ge.td)) go to 600
428 +
c        dans le cas quadratique prendre extrp plus grand
429 +
        text=min(text,extra*extrp*t)
430 +
        text=max(text,2.50d+0*t)
431 +
      else
432 +
        text=extra*extrp*t
433 +
      end if
434 +
      ta1=t
435 +
      fa1=fn
436 +
      hpta1=hpn
437 +
      extrp=10.
438 +
      if (text.ge.tmax/2.0d+0) then
439 +
        text=tmax
440 +
        imax=1
441 +
      end if
442 +
      if ((t.lt.tproj).and.(text.gt.tproj)) then
443 +
        t=max(tproj,2.50d+0*t)
444 +
        icoi=0
445 +
        icos=icop
446 +
        if(d(icop).lt.0.0d+0) then
447 +
          icoi=icop
448 +
          icos=0
449 +
        end if
450 +
        var2='es '
451 +
        goto 800
452 +
      end if
453 +
      ttsup=min(1.50d+0*text,tmax)
454 +
      if (ttsup.lt.tproj) go to 785
455 +
      ttmin=2*t
456 +
      iproj=0
457 +
      tmi=text
458 +
      topt=0.0d+0
459 +
      call satur(n,x,binf,bsup,d,ttmin,ttsup,topt,tg,td,tmi,
460 +
     &            icoi,icos,iproj)
461 +
      if (topt.gt.0.0d+0)then
462 +
          t=topt
463 +
          var2='es '
464 +
          go to 800
465 +
      endif
466 +
785   t=text
467 +
      var2='e  '
468 +
800   f11=fn-f
469 +
c$$$      if (iprint.ge.3.and.indic.gt.0) then
470 +
c$$$        write (bufstr,15000)var2,ta1,f11,hpn,h1,t1
471 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
472 +
c$$$        endif
473 +
c
474 +
c      test sur deltat
475 +
      if(abs(ta1-t).ge.1.0d+2*zero) go to 200
476 +
      k=4
477 +
c      calcul de indrl
478 +
1000  if (indic.lt.0) then
479 +
        indrl=13
480 +
        if (tg.eq.0.0d+0) indrl=-1000+indic
481 +
        fn=ftg
482 +
        hpn=hptg
483 +
        t=tg
484 +
        go to 1010
485 +
      endif
486 +
      if (fn.le.ftg) then
487 +
      indrl=k
488 +
      t=tg
489 +
      go to 1010
490 +
      end if
491 +
      if (tg.eq.0.0d+0) then
492 +
       indrl=-1*k
493 +
      go to 1010
494 +
      end if
495 +
      indrl=10+k
496 +
      t=tg
497 +
      fn=ftg
498 +
      hpn=hptg
499 +
c
500 +
c      fin du programme
501 +
1010  f=fn
502 +
      do i=1,n
503 +
         x(i)=x(i)+t*d(i)
504 +
      end do
505 +
      call proj(n,binf,bsup,x)
506 +
      if (icos.gt.0) x(icos)=bsup(icos)
507 +
      if (icoi.gt.0) x(icoi)=binf(icoi)
508 +
c
509 +
      if (indrl.lt.0) then
510 +
        nap=nap+1
511 +
        indic=4
512 +
        call simul (indic,n,x,f,gn,izs,rzs,dzs)
513 +
      endif
514 +
c
515 +
      t1=t-ta1
516 +
      if (t1.eq.0.0d+0) then
517 +
      t1=1.
518 +
      end if
519 +
      h1=(fn-fa1)/t1
520 +
      hp=hpd
521 +
      f0=f-f0
522 +
c$$$      if (iprint.ge.3) then
523 +
c$$$        write (bufstr,15020)t,f0,hpd,h1,t1
524 +
c$$$        call basout(io_out ,io ,bufstr(1:lnblnk(bufstr)))
525 +
c$$$        endif
526 +
      return
527 +
      end
Files Coverage
src 22.93%
R/n1qn1.R 54.84%
Project Totals (12 files) 23.61%
Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading