COPASI API  4.16.103
CTruncatedNewton.cpp
Go to the documentation of this file.
1 // Copyright (C) 2010 - 2013 by Pedro Mendes, Virginia Tech Intellectual
2 // Properties, Inc., University of Heidelberg, and The University
3 // of Manchester.
4 // All rights reserved.
5 
6 // Copyright (C) 2008 - 2009 by Pedro Mendes, Virginia Tech Intellectual
7 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
8 // and The University of Manchester.
9 // All rights reserved.
10 
11 // Copyright (C) 2007 by Pedro Mendes, Virginia Tech Intellectual
12 // Properties, Inc. and EML Research, gGmbH.
13 // All rights reserved.
14 
15 /* tn.f -- translated by f2c (version 19970805).
16  You must link the resulting object file with the libraries:
17  -lf2c -lm (in that order)
18  */
19 
20 #include <cmath>
21 
22 #include "copasi.h"
23 #include "CTruncatedNewton.h"
24 
25 #include "lapack/blaswrap.h"
26 #include "lapack/lapackwrap.h"
27 
29  C_FLOAT64 *smax);
30 int lsout_(C_INT *, C_INT *, C_FLOAT64 *,
34 int initp3_(C_FLOAT64 *, C_FLOAT64 *, C_INT *,
36  C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *);
37 int ssbfgs_(C_INT *, C_FLOAT64 *, C_FLOAT64 *,
40 int mslv_(C_FLOAT64 *, C_FLOAT64 *, C_INT *,
43  C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *);
44 C_FLOAT64 mchpr1_(void);
45 int negvec_(C_INT *, C_FLOAT64 *);
46 int ztime_(C_INT *, C_FLOAT64 *, C_INT *);
47 int ndia3_(C_INT *, C_FLOAT64 *, C_FLOAT64 *,
48  C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *);
49 int stpmax_(C_FLOAT64 *, C_FLOAT64 *,
50  C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *,
51  C_FLOAT64 *, C_FLOAT64 *);
52 int cnvtst_(C_INT *, C_FLOAT64 *,
56  C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *);
57 
58 int chkucp_(C_INT *, C_INT *, C_INT *,
62  C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *);
64  C_INT *, C_INT *, C_INT *, C_INT *, C_INT *);
65 int crash_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *);
66 int dxpy_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *);
67 int modz_(C_INT *, C_FLOAT64 *, C_FLOAT64 *,
69  C_FLOAT64 *);
70 
71 C_FLOAT64 mchpr1_(void);
72 
73 #define subscr_1 (mpsubscr_->_1)
74 #define subscr_2 (mpsubscr_->_2)
75 #define subscr_3 (mpsubscr_->_3)
76 
77 /* Common Block Declarations */
78 #define TRUE_ (1)
79 #define FALSE_ (0)
80 
81 /* Table of constant values */
82 
83 static C_INT c__1 = 1;
84 static C_INT c_false = FALSE_;
85 static C_INT c_true = TRUE_;
86 static C_FLOAT64 c_b246 = .6666;
87 
89 {
90  mpsubscr_ = new subscr_;
91 }
92 
94 {
95  if (mpsubscr_ != NULL) {delete mpsubscr_; mpsubscr_ = NULL;}
96 }
97 
98 /* %% TRUNCATED-NEWTON METHOD: SUBROUTINES */
99 /* FOR OTHER MACHINES, MODIFY ROUTINE MCHPR1 (MACHINE EPSILON) */
100 /* WRITTEN BY: STEPHEN G. NASH */
101 /* OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. */
102 /* GEORGE MASON UNIVERSITY */
103 /* FAIRFAX, VA 22030 */
104 /* ****************************************************************** */
105 /* Subroutine */ int CTruncatedNewton::tn_(C_INT *ierror, C_INT *n, C_FLOAT64 *x,
106  C_FLOAT64 *f, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *lw, FTruncatedNewton *sfun)
107 {
108  /* Format strings */
109  /* static char fmt_800[] = "(//,\002 ERROR CODE =\002,i3)";
110  static char fmt_810[] = "(//,\002 OPTIMAL FUNCTION VALUE = \002,1pd22.15)"
111  ;
112  static char fmt_820[] = "(10x,\002CURRENT SOLUTION IS (AT MOST 10 COMPON"
113  "ENTS)\002,/,14x,\002I\002,11x,\002X(I)\002)";
114  static char fmt_830[] = "(10x,i5,2x,1pd22.15)";*/
115 
116  /* Local variables */
117  C_FLOAT64 xtol;
118  C_INT maxit;
119  C_FLOAT64 accrcy;
120  C_INT maxfun, msglvl;
121  C_FLOAT64 stepmx, eta;
122 
123  /* Fortran I/O blocks */ /*
124  cilist io___8 = {0, 6, 0, fmt_800, 0 };
125  cilist io___9 = {0, 6, 0, fmt_810, 0 };
126  cilist io___10 = {0, 6, 0, fmt_820, 0 };
127  cilist io___12 = {0, 6, 0, fmt_830, 0 };*/
128 
129  /* THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM */
130 
131  /* MINIMIZE F(X) */
132  /* X */
133 
134  /* WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS */
135  /* A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA */
136  /* THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), */
137  /* PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES
138  */
139  /* NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A */
140  /* GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW.
141  */
142  /* IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS */
143  /* ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. */
144 
145  /* SUBROUTINE PARAMETERS: */
146 
147  /* IERROR - (C_INT) ERROR CODE */
148  /* (0 => NORMAL RETURN) */
149  /* (2 => MORE THAN MAXFUN EVALUATIONS) */
150  /* (3 => LINE SEARCH FAILED TO FIND */
151  /* (LOWER POINT (MAY NOT BE SERIOUS) */
152  /* (-1 => ERROR IN INPUT PARAMETERS) */
153  /* N - (C_INT) NUMBER OF VARIABLES */
154  /* X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL */
155  /* ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. */
156  /* G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL */
157  /* VALUE OF THE GRADIENT */
158  /* F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE */
159  /* OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE */
160  /* OF THE OBJECTIVE FUNCTION AT THE SOLUTION */
161  /* W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N */
162  /* LW - (C_INT) THE DECLARED DIMENSION OF W */
163  /* SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION */
164  /* AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE */
165  /* THE CALLING SEQUENCE */
166  /* SUBROUTINE SFUN (N, X, F, G) */
167  /* C_INT N */
168  /* DOUBLE PRECISION X(N), G(N), F */
169 
170  /* THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE */
171  /* LMQN. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE */
172  /* OF THIS ALGORITHM SHOULD CALL LMQN DIRECTLY. */
173 
174  /* ----------------------------------------------------------------------
175  */
176  /* THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON */
177  /* ALGORITHM. THE PARAMETERS ARE: */
178 
179  /* ETA - SEVERITY OF THE LINESEARCH */
180  /* MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS */
181  /* XTOL - DESIRED ACCURACY FOR THE SOLUTION X* */
182  /* STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH */
183  /* ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES */
184  /* MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT */
185  /* 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. */
186  /* MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP */
187 
188  /* SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE */
189 
190  /* Parameter adjustments */
191  --g;
192  --x;
193  --w;
194 
195  /* Function Body */
196  maxit = *n / 2;
197 
198  if (maxit > 50)
199  {
200  maxit = 50;
201  }
202 
203  if (maxit <= 0)
204  {
205  maxit = 1;
206  }
207 
208  msglvl = 0;
209  maxfun = *n * 150;
210  eta = .25;
211  stepmx = 10.;
212  accrcy = mchpr1_() * 100.;
213  xtol = sqrt(accrcy);
214 
215  /* MINIMIZE THE FUNCTION */
216 
217  CTruncatedNewton::lmqn_(ierror, n, &x[1], f, &g[1], &w[1], lw, sfun, &msglvl, &maxit,
218  &maxfun, &eta, &stepmx, &accrcy, &xtol);
219 
220  /* PRINT THE RESULTS */
221  //remove the printing
222  /* if (*ierror != 0) {
223  s_wsfe(&io___8);
224  do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(C_INT));
225  e_wsfe();
226  }
227  s_wsfe(&io___9);
228  do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(C_FLOAT64));
229  e_wsfe();
230  if (msglvl < 1) {
231  return 0;
232  }
233  s_wsfe(&io___10);
234  e_wsfe();
235  nmax = 10;
236  if (*n < nmax) {
237  nmax = *n;
238  }
239  s_wsfe(&io___12);
240  i__1 = nmax;
241  for (i__ = 1; i__ <= i__1; ++i__) {
242  do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(C_INT));
243  do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(C_FLOAT64));
244  }
245  e_wsfe();*/
246  return 0;
247 } /* tn_ */
248 
249 /* Subroutine */ int CTruncatedNewton::tnbc_(C_INT *ierror, C_INT *n, C_FLOAT64 *x,
250  C_FLOAT64 *f, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *lw, FTruncatedNewton * sfun,
251  C_FLOAT64 *low, C_FLOAT64 *up, C_INT *ipivot)
252 {
253  /* Format strings */
254  /* char fmt_800[] = "(//,\002 ERROR CODE =\002,i3)";
255  char fmt_810[] = "(//,\002 OPTIMAL FUNCTION VALUE = \002,1pd22.15)"
256  ;
257  char fmt_820[] = "(10x,\002CURRENT SOLUTION IS (AT MOST 10 COMPON"
258  "ENTS)\002,/,14x,\002I\002,11x,\002X(I)\002)";
259  char fmt_830[] = "(10x,i5,2x,1pd22.15)";*/
260 
261  /* System generated locals */
262  C_INT i__1;
263 
264  // C_INT s_wsfe(cilist *), do_fio(C_INT *, char *, ftnlen), e_wsfe(void);
265 
266  /* Local variables */
267  C_INT nmax;
268  C_FLOAT64 xtol;
269  C_INT i__, maxit;
270  C_FLOAT64 accrcy;
271  C_INT maxfun, msglvl;
272  C_FLOAT64 stepmx, eta;
273 
274  /* Fortran I/O blocks */
275  //remove
276  /*
277  cilist io___21 = {0, 6, 0, fmt_800, 0 };
278  cilist io___22 = {0, 6, 0, fmt_810, 0 };
279  cilist io___23 = {0, 6, 0, fmt_820, 0 };
280  cilist io___25 = {0, 6, 0, fmt_830, 0 };*/
281 
282  /* THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM */
283 
284  /* MINIMIZE F(X) */
285  /* X */
286  /* SUBJECT TO LOW <= X <= UP */
287 
288  /* WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS */
289  /* A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA */
290  /* THE LANCZOS ALGORITHM" BY S.G. NASH (TECHNICAL REPORT 378, MATH. */
291  /* THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), */
292  /* PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES
293  */
294  /* NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A */
295  /* GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW.
296  */
297  /* IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS */
298  /* ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. */
299 
300  /* SUBROUTINE PARAMETERS: */
301 
302  /* IERROR - (C_INT) ERROR CODE */
303  /* (0 => NORMAL RETURN */
304  /* (2 => MORE THAN MAXFUN EVALUATIONS */
305  /* (3 => LINE SEARCH FAILED TO FIND LOWER */
306  /* (POINT (MAY NOT BE SERIOUS) */
307  /* (-1 => ERROR IN INPUT PARAMETERS */
308  /* N - (C_INT) NUMBER OF VARIABLES */
309  /* X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL */
310  /* ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION.
311  */
312  /* G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL */
313  /* VALUE OF THE GRADIENT */
314  /* F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE */
315  /* OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE */
316  /* OF THE OBJECTIVE FUNCTION AT THE SOLUTION */
317  /* W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N */
318  /* LW - (C_INT) THE DECLARED DIMENSION OF W */
319  /* SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION */
320  /* AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE */
321  /* THE CALLING SEQUENCE */
322  /* SUBROUTINE SFUN (N, X, F, G) */
323  /* C_INT N */
324  /* DOUBLE PRECISION X(N), G(N), F */
325  /* LOW, UP - (REAL*8) VECTORS OF LENGTH AT LEAST N CONTAINING */
326  /* THE LOWER AND UPPER BOUNDS ON THE VARIABLES. IF */
327  /* THERE ARE NO BOUNDS ON A PARTICULAR VARIABLE, SET */
328  /* THE BOUNDS TO -1.D38 AND 1.D38, RESPECTIVELY. */
329  /* IPIVOT - (C_INT) WORK VECTOR OF LENGTH AT LEAST N, USED */
330  /* TO RECORD WHICH VARIABLES ARE AT THEIR BOUNDS. */
331 
332  /* THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE */
333  /* LMQNBC. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE */
334  /* OF THIS ALGORITHM SHOULD CALL LMQBC DIRECTLY. */
335 
336  /* ----------------------------------------------------------------------
337  */
338  /* THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON */
339  /* ALGORITHM. THE PARAMETERS ARE: */
340 
341  /* ETA - SEVERITY OF THE LINESEARCH */
342  /* MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS */
343  /* XTOL - DESIRED ACCURACY FOR THE SOLUTION X* */
344  /* STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH */
345  /* ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES */
346  /* MSGLVL - CONTROLS QUANTITY OF PRINTED OUTPUT */
347  /* 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. */
348  /* MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP */
349 
350  /* SET PARAMETERS FOR THE OPTIMIZATION ROUTINE */
351 
352  /* Parameter adjustments */
353  --ipivot;
354  --up;
355  --low;
356  --g;
357  --x;
358  --w;
359 
360  /* Function Body */
361  maxit = *n / 2;
362 
363  if (maxit > 50)
364  {
365  maxit = 50;
366  }
367 
368  if (maxit <= 0)
369  {
370  maxit = 1;
371  }
372 
373  msglvl = 0;
374  maxfun = *n * 150;
375  eta = .25;
376  stepmx = 10.;
377  accrcy = mchpr1_() * 100.;
378  xtol = sqrt(accrcy);
379 
380  /* MINIMIZE FUNCTION */
381 
382  CTruncatedNewton::lmqnbc_(ierror, n, &x[1], f, &g[1], &w[1], lw, sfun, &low[1], &up[1]
383  , &ipivot[1], &msglvl, &maxit, &maxfun, &eta, &stepmx, &accrcy, &
384  xtol);
385 
386  /* PRINT RESULTS */
387 
388  if (*ierror != 0)
389  {
390  // s_wsfe(&io___21);
391  // do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(C_INT));
392  // e_wsfe();
393  }
394 
395  // s_wsfe(&io___22);
396  // do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(C_FLOAT64));
397  // e_wsfe();
398  if (msglvl < 1)
399  {
400  return 0;
401  }
402 
403  // s_wsfe(&io___23);
404  // e_wsfe();
405  nmax = 10;
406 
407  if (*n < nmax)
408  {
409  nmax = *n;
410  }
411 
412  // s_wsfe(&io___25);
413  i__1 = nmax;
414 
415  for (i__ = 1; i__ <= i__1; ++i__)
416  {
417  // do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(C_INT));
418  // do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(C_FLOAT64));
419  }
420 
421  // e_wsfe();
422  return 0;
423 } /* tnbc_ */
424 
425 /* Subroutine */ int CTruncatedNewton::lmqn_(C_INT *ifail, C_INT *n, C_FLOAT64 *x,
426  C_FLOAT64 *f, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *lw, FTruncatedNewton * sfun,
427  C_INT *msglvl, C_INT *maxit, C_INT *maxfun, C_FLOAT64 *eta,
428  C_FLOAT64 *stepmx, C_FLOAT64 *accrcy, C_FLOAT64 *xtol)
429 {
430  /* Format strings */
431  /* char fmt_800[] = "(//\002 NIT NF CG\002,9x,\002F\002,21x,\002"
432  "GTG\002,//)";
433  char fmt_810[] = "(\002 \002,i3,1x,i4,1x,i4,1x,1pd22.15,2x,1pd15."
434  "8)"; */
435 
436  /* System generated locals */
437  C_INT i__1;
438 
439  /* Builtin functions */
440  // C_INT s_wsfe(cilist *), e_wsfe(void), do_fio(C_INT *, char *, ftnlen);
441 
442  /* Local variables */
443  C_FLOAT64 fold, oldf;
444  C_FLOAT64 fnew;
445  C_INT numf;
446  C_FLOAT64 peps;
447  C_INT lhyr;
448  C_FLOAT64 zero, rtol, yksk, tiny;
449  C_INT nwhy;
450  C_FLOAT64 yrsr;
451  C_INT i__;
452  C_FLOAT64 alpha, fkeep;
453  C_INT ioldg;
454  C_FLOAT64 small;
455  C_INT modet;
456  C_INT niter;
457  C_FLOAT64 gnorm, ftest, fstop, pnorm, rteps, xnorm;
458  C_INT idiagb;
459  C_FLOAT64 fm, pe, difold;
460  C_INT icycle, nlincg, nfeval;
461  C_FLOAT64 difnew;
462  C_INT nmodif;
463  C_FLOAT64 epsmch;
464  C_FLOAT64 epsred, fabstol, oldgtp;
465  C_INT ireset;
466  C_INT lreset;
467  C_FLOAT64 reltol, gtpnew;
468  C_INT nftotl;
469  C_FLOAT64 toleps;
470  C_FLOAT64 rtleps;
471  C_INT ipivot, nm1;
472  C_FLOAT64 rtolsq, tnytol, gtg, one;
473  C_INT ipk;
474  C_FLOAT64 gsk, spe;
475  C_INT isk, iyk;
476  C_INT upd1;
477 
478  /* Fortran I/O blocks */
479  //remove the blocks
480  /* cilist io___27 = {0, 6, 0, fmt_800, 0 };
481  cilist io___56 = {0, 6, 0, fmt_810, 0 };
482  cilist io___81 = {0, 6, 0, fmt_810, 0 };
483  */
484 
485  /* THIS ROUTINE IS A TRUNCATED-NEWTON METHOD. */
486  /* THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY */
487  /* QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED */
488  /* IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3).
489  */
490  /* FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TN. */
491 
492  /* THE FOLLOWING IMSL AND STANDARD FUNCTIONS ARE USED */
493 
494  /* INITIALIZE PARAMETERS AND CONSTANTS */
495 
496  /* Parameter adjustments */
497  --g;
498  --x;
499  --w;
500 
501  /* Function Body */
502  /* if (*msglvl >= -2) {
503  s_wsfe(&io___27);
504  e_wsfe();
505  } */
506  setpar_(n);
507  upd1 = TRUE_;
508  ireset = 0;
509  nfeval = 0;
510  nmodif = 0;
511  nlincg = 0;
512  fstop = *f;
513  zero = 0.;
514  one = 1.;
515  nm1 = *n - 1;
516 
517  /* WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) */
518 
519  lhyr = subscr_1.loldg;
520 
521  /* CHECK PARAMETERS AND SET CONSTANTS */
522 
523  chkucp_(&subscr_1.lwtest, maxfun, &nwhy, n, &alpha, &epsmch, eta, &peps, &
524  rteps, &rtol, &rtolsq, stepmx, &ftest, xtol, &xnorm, &x[1], lw, &
525  small, &tiny, accrcy);
526 
527  if (nwhy < 0)
528  {
529  goto L120;
530  }
531 
532  CTruncatedNewton::setucr_(&small, &nftotl, &niter, n, f, &fnew, &fm, &gtg, &oldf,
533  sfun, &g[1], &x[1]);
534  fold = fnew;
535 
536  if (*msglvl >= 1)
537  {
538  /* s_wsfe(&io___56);
539  do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(C_INT));
540  do_fio(&c__1, (char *)&nftotl, (ftnlen)sizeof(C_INT));
541  do_fio(&c__1, (char *)&nlincg, (ftnlen)sizeof(C_INT));
542  do_fio(&c__1, (char *)&fnew, (ftnlen)sizeof(C_FLOAT64));
543  do_fio(&c__1, (char *)&gtg, (ftnlen)sizeof(C_FLOAT64));
544  e_wsfe();*/
545  }
546 
547  /* CHECK FOR SMALL GRADIENT AT THE STARTING POINT. */
548 
549  ftest = one + fabs(fnew);
550 
551  if (gtg < epsmch * 1e-4 * ftest * ftest)
552  {
553  goto L90;
554  }
555 
556  /* SET INITIAL VALUES TO OTHER PARAMETERS */
557 
558  icycle = nm1;
559  toleps = rtol + rteps;
560  rtleps = rtolsq + epsmch;
561  gnorm = sqrt(gtg);
562  difnew = zero;
563  epsred = .05;
564  fkeep = fnew;
565 
566  /* SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. */
567 
568  idiagb = subscr_1.ldiagb;
569  i__1 = *n;
570 
571  for (i__ = 1; i__ <= i__1; ++i__)
572  {
573  w[idiagb] = one;
574  ++idiagb;
575  /* L10: */
576  }
577 
578  /* ..................START OF MAIN ITERATIVE LOOP.......... */
579 
580  /* COMPUTE THE NEW SEARCH DIRECTION */
581 
582  modet = *msglvl - 3;
583  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
584  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
585  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
586  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
587  c_false, &ipivot, accrcy, &gtpnew, &gnorm, &xnorm);
588 L20:
589  dcopy_(n, &g[1], &c__1, &w[subscr_1.loldg], &c__1);
590  pnorm = dnrm2_(n, &w[subscr_1.lpk], &c__1);
591  oldf = fnew;
592  oldgtp = gtpnew;
593 
594  /* PREPARE TO COMPUTE THE STEP LENGTH */
595 
596  pe = pnorm + epsmch;
597 
598  /* COMPUTE THE fabsOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH */
599 
600  reltol = rteps * (xnorm + one) / pe;
601  fabstol = -epsmch * ftest / (oldgtp - epsmch);
602 
603  /* COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN */
604  /* THE LINEAR SEARCH */
605 
606  tnytol = epsmch * (xnorm + one) / pe;
607  spe = *stepmx / pe;
608 
609  /* SET THE INITIAL STEP LENGTH. */
610 
611  alpha = step1_(&fnew, &fm, &oldgtp, &spe);
612 
613  /* PERFORM THE LINEAR SEARCH */
614 
615  CTruncatedNewton::linder_(n, sfun, &small, &epsmch, &reltol, &fabstol, &tnytol, eta, &
616  zero, &spe, &w[subscr_1.lpk], &oldgtp, &x[1], &fnew, &alpha, &g[1]
617  , &numf, &nwhy, &w[1], lw);
618 
619  fold = fnew;
620  ++niter;
621  nftotl += numf;
622  gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
623 
624  if (*msglvl >= 1)
625  {
626  //remove the printing
627  /*
628  s_wsfe(&io___81);
629  do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(C_INT));
630  do_fio(&c__1, (char *)&nftotl, (ftnlen)sizeof(C_INT));
631  do_fio(&c__1, (char *)&nlincg, (ftnlen)sizeof(C_INT));
632  do_fio(&c__1, (char *)&fnew, (ftnlen)sizeof(C_FLOAT64));
633  do_fio(&c__1, (char *)&gtg, (ftnlen)sizeof(C_FLOAT64));
634  e_wsfe();*/
635  }
636 
637  if (nwhy < 0)
638  {
639  goto L120;
640  }
641 
642  if (nwhy == 0 || nwhy == 2)
643  {
644  goto L30;
645  }
646 
647  /* THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT */
648 
649  nwhy = 3;
650  goto L100;
651 L30:
652 
653  if (nwhy <= 1)
654  {
655  goto L40;
656  }
657 
658  (*sfun)(n, &x[1], &fnew, &g[1]);
659  ++nftotl;
660 
661  /* TERMINATE IF MORE THAN MAXFUN EVALUTATIONS HAVE BEEN MADE */
662 
663 L40:
664  nwhy = 2;
665 
666  if (nftotl > *maxfun)
667  {
668  goto L110;
669  }
670 
671  nwhy = 0;
672 
673  /* SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS */
674 
675  difold = difnew;
676  difnew = oldf - fnew;
677 
678  /* IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE */
679  /* PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. */
680 
681  if (icycle != 1)
682  {
683  goto L50;
684  }
685 
686  if (difnew > difold * 2.)
687  {
688  epsred += epsred;
689  }
690 
691  if (difnew < difold * .5)
692  {
693  epsred *= .5;
694  }
695 
696 L50:
697  gnorm = sqrt(gtg);
698  ftest = one + fabs(fnew);
699  xnorm = dnrm2_(n, &x[1], &c__1);
700 
701  /* TEST FOR CONVERGENCE */
702 
703  if ((alpha * pnorm < toleps * (one + xnorm) &&
704  fabs(difnew) < rtleps * ftest &&
705  gtg < peps * ftest * ftest) ||
706  gtg < *accrcy * 1e-4 * ftest * ftest)
707  {
708  goto L90;
709  }
710 
711  /* COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE */
712  /* IN THE GRADIENTS */
713 
714  isk = subscr_1.lsk;
715  ipk = subscr_1.lpk;
716  iyk = subscr_1.lyk;
717  ioldg = subscr_1.loldg;
718  i__1 = *n;
719 
720  for (i__ = 1; i__ <= i__1; ++i__)
721  {
722  w[iyk] = g[i__] - w[ioldg];
723  w[isk] = alpha * w[ipk];
724  ++ipk;
725  ++isk;
726  ++iyk;
727  ++ioldg;
728  /* L60: */
729  }
730 
731  /* SET UP PARAMETERS USED IN UPDATING THE DIRECTION OF SEARCH. */
732 
733  yksk = ddot_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lsk], &c__1);
734  lreset = FALSE_;
735 
736  if (icycle == nm1 || difnew < epsred * (fkeep - fnew))
737  {
738  lreset = TRUE_;
739  }
740 
741  if (lreset)
742  {
743  goto L70;
744  }
745 
746  yrsr = ddot_(n, &w[subscr_1.lyr], &c__1, &w[subscr_1.lsr], &c__1);
747 
748  if (yrsr <= zero)
749  {
750  lreset = TRUE_;
751  }
752 
753 L70:
754  upd1 = FALSE_;
755 
756  /* COMPUTE THE NEW SEARCH DIRECTION */
757 
758  modet = *msglvl - 3;
759  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
760  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
761  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
762  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
763  c_false, &ipivot, accrcy, &gtpnew, &gnorm, &xnorm);
764 
765  if (lreset)
766  {
767  goto L80;
768  }
769 
770  /* STORE THE ACCUMULATED CHANGE IN THE POINT AND GRADIENT AS AN */
771  /* "AVERAGE" DIRECTION FOR PRECONDITIONING. */
772 
773  dxpy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
774  dxpy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
775  ++icycle;
776  goto L20;
777 
778  /* RESET */
779 
780 L80:
781  ++ireset;
782 
783  /* INITIALIZE THE SUM OF ALL THE CHANGES IN X. */
784 
785  dcopy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
786  dcopy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
787  fkeep = fnew;
788  icycle = 1;
789  goto L20;
790 
791  /* ...............END OF MAIN ITERATION....................... */
792 
793 L90:
794  *ifail = 0;
795  *f = fnew;
796  return 0;
797 L100:
798  oldf = fnew;
799 
800  /* LOCAL SEARCH HERE COULD BE INSTALLED HERE */
801 
802 L110:
803  *f = oldf;
804 
805  /* SET IFAIL */
806 
807 L120:
808  *ifail = nwhy;
809  return 0;
810 } /* lmqn_ */
811 
812 /* Subroutine */ int CTruncatedNewton::lmqnbc_(C_INT *ifail, C_INT *n, C_FLOAT64 *x,
813  C_FLOAT64 *f, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *lw, FTruncatedNewton *sfun,
814  C_FLOAT64 *low, C_FLOAT64 *up, C_INT *ipivot, C_INT *msglvl,
815  C_INT *maxit, C_INT *maxfun, C_FLOAT64 *eta, C_FLOAT64 *stepmx,
816  C_FLOAT64 *accrcy, C_FLOAT64 *xtol)
817 {
818  /* Format strings */
819  //remove
820  /*
821  char fmt_800[] = "(\002 THERE IS NO FEASIBLE POINT; TERMINATING A"
822  "LGORITHM\002)";
823  char fmt_810[] = "(//\002 NIT NF CG\002,9x,\002F\002,21x,"
824  "\002GTG\002,//)";
825  char fmt_820[] = "(\002 LINESEARCH RESULTS: ALPHA,PNOR"
826  "M\002,2(1pd12.4))";
827  char fmt_830[] = "(\002 UPD1 IS TRUE - TRIVIAL PRECONDITIONING"
828  "\002)";
829  char fmt_840[] = "(\002 NEWCON IS TRUE - CONSTRAINT ADDED IN LINE"
830  "SEARCH\002)";
831  */
832 
833  /* System generated locals */
834  C_INT i__1;
835  C_FLOAT64 d__1;
836 
837  /* Builtin functions */
838  // C_INT s_wsfe(cilist *), e_wsfe(void);
839  // C_INT do_fio(C_INT *, char *, ftnlen);
840 
841  /* Local variables */
842  C_FLOAT64 fold, oldf;
843  C_FLOAT64 fnew;
844  C_INT numf;
845  C_INT conv;
846  C_FLOAT64 peps;
847  C_INT lhyr;
848  C_FLOAT64 zero, rtol, yksk, tiny;
849  C_INT nwhy;
850  C_FLOAT64 yrsr;
851  C_INT i__;
852  C_FLOAT64 alpha, fkeep;
853  C_INT ioldg;
854  C_FLOAT64 small, flast;
855  C_INT modet;
856  C_INT niter;
857  C_FLOAT64 gnorm, ftest;
858  C_FLOAT64 fstop, pnorm, rteps, xnorm;
859  C_INT idiagb;
860  C_FLOAT64 fm, pe, difold;
861  C_INT icycle, nlincg, nfeval;
862  C_FLOAT64 difnew;
863  C_INT nmodif;
864  C_FLOAT64 epsmch;
865  C_FLOAT64 epsred, fabstol, oldgtp;
866  C_INT newcon;
867  C_INT ireset;
868  C_INT lreset;
869  C_FLOAT64 reltol, gtpnew;
870  C_INT nftotl;
871  C_FLOAT64 toleps;
872  C_FLOAT64 rtleps;
873  C_INT nm1;
874  C_FLOAT64 rtolsq, tnytol;
875  C_INT ier;
876  C_FLOAT64 gtg, one;
877  C_INT ipk;
878  C_FLOAT64 gsk, spe;
879  C_INT isk, iyk;
880  C_INT upd1;
881 
882  /* Fortran I/O blocks */
883  /*
884  cilist io___88 = {0, 6, 0, fmt_800, 0 };
885  cilist io___89 = {0, 6, 0, fmt_810, 0 };
886  cilist io___144 = {0, 6, 0, fmt_820, 0 };
887  cilist io___150 = {0, 6, 0, fmt_830, 0 };
888  cilist io___151 = {0, 6, 0, fmt_840, 0 }; */
889 
890  /* THIS ROUTINE IS A BOUNDS-CONSTRAINED TRUNCATED-NEWTON METHOD. */
891  /* THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY */
892  /* QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED */
893  /* IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3).
894  */
895  /* FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TNBC. */
896 
897  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE USED */
898 
899  /* CHECK THAT INITIAL X IS FEASIBLE AND THAT THE BOUNDS ARE CONSISTENT */
900 
901  /* Parameter adjustments */
902  --ipivot;
903  --up;
904  --low;
905  --g;
906  --x;
907  --w;
908 
909  /* Function Body */
910  crash_(n, &x[1], &ipivot[1], &low[1], &up[1], &ier);
911  /*if (ier != 0) {
912  s_wsfe(&io___88);
913  e_wsfe();
914  }
915  if (ier != 0) {
916  return 0;
917  }
918  if (*msglvl >= 1) {
919  s_wsfe(&io___89);
920  e_wsfe();
921  } */
922 
923  /* INITIALIZE VARIABLES */
924 
925  setpar_(n);
926  upd1 = TRUE_;
927  ireset = 0;
928  nfeval = 0;
929  nmodif = 0;
930  nlincg = 0;
931  fstop = *f;
932  conv = FALSE_;
933  zero = 0.;
934  one = 1.;
935  nm1 = *n - 1;
936 
937  /* WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) */
938 
939  lhyr = subscr_1.loldg;
940 
941  /* CHECK PARAMETERS AND SET CONSTANTS */
942 
943  chkucp_(&subscr_1.lwtest, maxfun, &nwhy, n, &alpha, &epsmch, eta, &peps, &
944  rteps, &rtol, &rtolsq, stepmx, &ftest, xtol, &xnorm, &x[1], lw, &
945  small, &tiny, accrcy);
946 
947  if (nwhy < 0)
948  {
949  goto L160;
950  }
951 
952  CTruncatedNewton::setucr_(&small, &nftotl, &niter, n, f, &fnew, &fm, &gtg, &oldf,
953  sfun, &g[1], &x[1]);
954  fold = fnew;
955  flast = fnew;
956 
957  /* TEST THE LAGRANGE MULTIPLIERS TO SEE IF THEY ARE NON-NEGATIVE. */
958  /* BECAUSE THE CONSTRAINTS ARE ONLY LOWER BOUNDS, THE COMPONENTS */
959  /* OF THE GRADIENT CORRESPONDING TO THE ACTIVE CONSTRAINTS ARE THE */
960  /* LAGRANGE MULTIPLIERS. AFTERWORDS, THE PROJECTED GRADIENT IS FORMED. */
961 
962  i__1 = *n;
963 
964  for (i__ = 1; i__ <= i__1; ++i__)
965  {
966  if (ipivot[i__] == 2)
967  {
968  goto L10;
969  }
970 
971  if (-ipivot[i__] * g[i__] >= 0.)
972  {
973  goto L10;
974  }
975 
976  ipivot[i__] = 0;
977 L10:
978  ;
979  }
980 
981  ztime_(n, &g[1], &ipivot[1]);
982  gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
983 
984  if (*msglvl >= 1)
985  {
986  monit_(n, &x[1], &fnew, &g[1], &niter, &nftotl, &nfeval, &lreset, &
987  ipivot[1]);
988  }
989 
990  /* CHECK IF THE INITIAL POINT IS A LOCAL MINIMUM. */
991 
992  ftest = one + fabs(fnew);
993 
994  if (gtg < epsmch * 1e-4 * ftest * ftest)
995  {
996  goto L130;
997  }
998 
999  /* SET INITIAL VALUES TO OTHER PARAMETERS */
1000 
1001  icycle = nm1;
1002  toleps = rtol + rteps;
1003  rtleps = rtolsq + epsmch;
1004  gnorm = sqrt(gtg);
1005  difnew = zero;
1006  epsred = .05;
1007  fkeep = fnew;
1008 
1009  /* SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. */
1010 
1011  idiagb = subscr_1.ldiagb;
1012  i__1 = *n;
1013 
1014  for (i__ = 1; i__ <= i__1; ++i__)
1015  {
1016  w[idiagb] = one;
1017  ++idiagb;
1018  /* L15: */
1019  }
1020 
1021  /* ..................START OF MAIN ITERATIVE LOOP.......... */
1022 
1023  /* COMPUTE THE NEW SEARCH DIRECTION */
1024 
1025  modet = *msglvl - 3;
1026  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
1027  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
1028  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
1029  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
1030  c_true, &ipivot[1], accrcy, &gtpnew, &gnorm, &xnorm);
1031 L20:
1032 
1033  /* added manually by Pedro Mendes 12/2/1998 */
1034  // if(callback != 0/*NULL*/) callback(fnew);
1035 
1036  dcopy_(n, &g[1], &c__1, &w[subscr_1.loldg], &c__1);
1037  pnorm = dnrm2_(n, &w[subscr_1.lpk], &c__1);
1038  oldf = fnew;
1039  oldgtp = gtpnew;
1040 
1041  /* PREPARE TO COMPUTE THE STEP LENGTH */
1042 
1043  pe = pnorm + epsmch;
1044 
1045  /* COMPUTE THE fabsOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH */
1046 
1047  reltol = rteps * (xnorm + one) / pe;
1048  fabstol = -epsmch * ftest / (oldgtp - epsmch);
1049 
1050  /* COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN */
1051  /* THE LINEAR SEARCH */
1052 
1053  tnytol = epsmch * (xnorm + one) / pe;
1054  stpmax_(stepmx, &pe, &spe, n, &x[1], &w[subscr_1.lpk], &ipivot[1], &low[1]
1055  , &up[1]);
1056 
1057  /* SET THE INITIAL STEP LENGTH. */
1058 
1059  alpha = step1_(&fnew, &fm, &oldgtp, &spe);
1060 
1061  /* PERFORM THE LINEAR SEARCH */
1062 
1063  CTruncatedNewton::linder_(n, sfun, &small, &epsmch, &reltol, &fabstol, &tnytol, eta, &
1064  zero, &spe, &w[subscr_1.lpk], &oldgtp, &x[1], &fnew, &alpha, &g[1]
1065  , &numf, &nwhy, &w[1], lw);
1066  newcon = FALSE_;
1067 
1068  if ((d__1 = alpha - spe, fabs(d__1)) > epsmch * 10.)
1069  {
1070  goto L30;
1071  }
1072 
1073  newcon = TRUE_;
1074  nwhy = 0;
1075  modz_(n, &x[1], &w[subscr_1.lpk], &ipivot[1], &epsmch, &low[1], &up[1], &
1076  flast, &fnew);
1077  flast = fnew;
1078 
1079 L30:
1080 
1081  if (*msglvl >= 3)
1082  {
1083  /* s_wsfe(&io___144);
1084  do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(C_FLOAT64));
1085  do_fio(&c__1, (char *)&pnorm, (ftnlen)sizeof(C_FLOAT64));
1086  e_wsfe();*/
1087  }
1088 
1089  fold = fnew;
1090  ++niter;
1091  nftotl += numf;
1092 
1093  /* IF REQUIRED, PRINT THE DETAILS OF THIS ITERATION */
1094 
1095  /* ############################ */
1096  if (*msglvl >= 1)
1097  {
1098  monit_(n, &x[1], &fnew, &g[1], &niter, &nftotl, &nfeval, &lreset, &
1099  ipivot[1]);
1100  }
1101 
1102  if (nwhy < 0)
1103  {
1104  goto L160;
1105  }
1106 
1107  if (nwhy == 0 || nwhy == 2)
1108  {
1109  goto L40;
1110  }
1111 
1112  /* THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT */
1113 
1114  nwhy = 3;
1115  goto L140;
1116 L40:
1117 
1118  if (nwhy <= 1)
1119  {
1120  goto L50;
1121  }
1122 
1123  (*sfun)(n, &x[1], &fnew, &g[1]);
1124  ++nftotl;
1125 
1126  /* TERMINATE IF MORE THAN MAXFUN EVALUATIONS HAVE BEEN MADE */
1127 
1128 L50:
1129  nwhy = 2;
1130 
1131  if (nftotl > *maxfun)
1132  {
1133  goto L150;
1134  }
1135 
1136  nwhy = 0;
1137 
1138  /* SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS */
1139 
1140  difold = difnew;
1141  difnew = oldf - fnew;
1142 
1143  /* IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE */
1144  /* PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. */
1145 
1146  if (icycle != 1)
1147  {
1148  goto L60;
1149  }
1150 
1151  if (difnew > difold * 2.)
1152  {
1153  epsred += epsred;
1154  }
1155 
1156  if (difnew < difold * .5)
1157  {
1158  epsred *= .5;
1159  }
1160 
1161 L60:
1162  dcopy_(n, &g[1], &c__1, &w[subscr_1.lgv], &c__1);
1163  ztime_(n, &w[subscr_1.lgv], &ipivot[1]);
1164  gtg = ddot_(n, &w[subscr_1.lgv], &c__1, &w[subscr_1.lgv], &c__1);
1165  gnorm = sqrt(gtg);
1166  ftest = one + fabs(fnew);
1167  xnorm = dnrm2_(n, &x[1], &c__1);
1168 
1169  /* TEST FOR CONVERGENCE */
1170 
1171  cnvtst_(&conv, &alpha, &pnorm, &toleps, &xnorm, &difnew, &rtleps, &ftest,
1172  &gtg, &peps, &epsmch, &gtpnew, &fnew, &flast, &g[1], &ipivot[1],
1173  n, accrcy);
1174 
1175  if (conv)
1176  {
1177  goto L130;
1178  }
1179 
1180  ztime_(n, &g[1], &ipivot[1]);
1181 
1182  /* COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE */
1183  /* IN THE GRADIENTS */
1184 
1185  if (newcon)
1186  {
1187  goto L90;
1188  }
1189 
1190  isk = subscr_1.lsk;
1191  ipk = subscr_1.lpk;
1192  iyk = subscr_1.lyk;
1193  ioldg = subscr_1.loldg;
1194  i__1 = *n;
1195 
1196  for (i__ = 1; i__ <= i__1; ++i__)
1197  {
1198  w[iyk] = g[i__] - w[ioldg];
1199  w[isk] = alpha * w[ipk];
1200  ++ipk;
1201  ++isk;
1202  ++iyk;
1203  ++ioldg;
1204  /* L70: */
1205  }
1206 
1207  /* SET UP PARAMETERS USED IN UPDATING THE PRECONDITIONING STRATEGY. */
1208 
1209  yksk = ddot_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lsk], &c__1);
1210  lreset = FALSE_;
1211 
1212  if (icycle == nm1 || difnew < epsred * (fkeep - fnew))
1213  {
1214  lreset = TRUE_;
1215  }
1216 
1217  if (lreset)
1218  {
1219  goto L80;
1220  }
1221 
1222  yrsr = ddot_(n, &w[subscr_1.lyr], &c__1, &w[subscr_1.lsr], &c__1);
1223 
1224  if (yrsr <= zero)
1225  {
1226  lreset = TRUE_;
1227  }
1228 
1229 L80:
1230  upd1 = FALSE_;
1231 
1232  /* COMPUTE THE NEW SEARCH DIRECTION */
1233 
1234 L90:
1235 
1236  if (upd1 && *msglvl >= 3)
1237  {
1238  /*
1239  s_wsfe(&io___150);
1240  e_wsfe();*/
1241  }
1242 
1243  if (newcon && *msglvl >= 3)
1244  {
1245  /*
1246  s_wsfe(&io___151);
1247  e_wsfe();*/
1248  }
1249 
1250  modet = *msglvl - 3;
1251  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
1252  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
1253  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
1254  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
1255  c_true, &ipivot[1], accrcy, &gtpnew, &gnorm, &xnorm);
1256 
1257  if (newcon)
1258  {
1259  goto L20;
1260  }
1261 
1262  if (lreset)
1263  {
1264  goto L110;
1265  }
1266 
1267  /* COMPUTE THE ACCUMULATED STEP AND ITS CORRESPONDING */
1268  /* GRADIENT DIFFERENCE. */
1269 
1270  dxpy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
1271  dxpy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
1272  ++icycle;
1273  goto L20;
1274 
1275  /* RESET */
1276 
1277 L110:
1278  ++ireset;
1279 
1280  /* INITIALIZE THE SUM OF ALL THE CHANGES IN X. */
1281 
1282  dcopy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
1283  dcopy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
1284  fkeep = fnew;
1285  icycle = 1;
1286  goto L20;
1287 
1288  /* ...............END OF MAIN ITERATION....................... */
1289 
1290 L130:
1291  *ifail = 0;
1292  *f = fnew;
1293  return 0;
1294 L140:
1295  oldf = fnew;
1296 
1297  /* LOCAL SEARCH COULD BE INSTALLED HERE */
1298 
1299 L150:
1300  *f = oldf;
1301 
1302  if (*msglvl >= 1)
1303  {
1304  monit_(n, &x[1], f, &g[1], &niter, &nftotl, &nfeval, &ireset, &ipivot[
1305  1]);
1306  }
1307 
1308  /* SET IFAIL */
1309 
1310 L160:
1311  *ifail = nwhy;
1312  return 0;
1313 } /* lmqnbc_ */
1314 
1315 /* Subroutine */ int monit_(C_INT *n, C_FLOAT64 *x, C_FLOAT64 * /* f */,
1316  C_FLOAT64 *g, C_INT * /* niter */, C_INT * /* nftotl */, C_INT * /* nfeval */,
1317  C_INT * /* ireset */, C_INT *ipivot)
1318 {
1319  /* Format strings */
1320  // char fmt_800[] = "(\002 \002,i4,1x,i4,1x,i4,1x,1pd22.15,2x,1pd15." "8)";
1321 
1322  /* System generated locals */
1323  C_INT i__1;
1324 
1325  /* Builtin functions */
1326  //C_INT s_wsfe(cilist *), do_fio(C_INT *, char *, ftnlen), e_wsfe(void);
1327 
1328  /* Local variables */
1329  C_INT i__;
1330  C_FLOAT64 gtg;
1331 
1332  /* PRINT RESULTS OF CURRENT ITERATION */
1333 
1334  /* Parameter adjustments */
1335  --ipivot;
1336  --g;
1337  --x;
1338 
1339  /* Function Body */
1340  gtg = 0.;
1341  i__1 = *n;
1342 
1343  for (i__ = 1; i__ <= i__1; ++i__)
1344  {
1345  if (ipivot[i__] != 0)
1346  {
1347  goto L10;
1348  }
1349 
1350  gtg += g[i__] * g[i__];
1351 L10:
1352  ;
1353  }
1354 
1355  //remove the printing
1356  /* s_wsfe(&io___154);
1357  do_fio(&c__1, (char *)&(*niter), (ftnlen)sizeof(C_INT));
1358  do_fio(&c__1, (char *)&(*nftotl), (ftnlen)sizeof(C_INT));
1359  do_fio(&c__1, (char *)&(*nfeval), (ftnlen)sizeof(C_INT));
1360  do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(C_FLOAT64));
1361  do_fio(&c__1, (char *)&gtg, (ftnlen)sizeof(C_FLOAT64));
1362  e_wsfe();*/
1363  return 0;
1364 } /* monit_ */
1365 
1366 /* Subroutine */ int ztime_(C_INT *n, C_FLOAT64 *x, C_INT *ipivot)
1367 {
1368  /* System generated locals */
1369  C_INT i__1;
1370 
1371  /* Local variables */
1372  C_INT i__;
1373 
1374  /* THIS ROUTINE MULTIPLIES THE VECTOR X BY THE CONSTRAINT MATRIX Z */
1375 
1376  /* Parameter adjustments */
1377  --ipivot;
1378  --x;
1379 
1380  /* Function Body */
1381  i__1 = *n;
1382 
1383  for (i__ = 1; i__ <= i__1; ++i__)
1384  {
1385  if (ipivot[i__] != 0)
1386  {
1387  x[i__] = 0.;
1388  }
1389 
1390  /* L10: */
1391  }
1392 
1393  return 0;
1394 } /* ztime_ */
1395 
1396 /* Subroutine */ int stpmax_(C_FLOAT64 *stepmx, C_FLOAT64 *pe, C_FLOAT64 *
1397  spe, C_INT *n, C_FLOAT64 *x, C_FLOAT64 *p, C_INT *ipivot,
1398  C_FLOAT64 *low, C_FLOAT64 *up)
1399 {
1400  /* System generated locals */
1401  C_INT i__1;
1402  /* Local variables */
1403  C_INT i__;
1404  C_FLOAT64 t;
1405 
1406  /* COMPUTE THE MAXIMUM ALLOWABLE STEP LENGTH */
1407 
1408  /* Parameter adjustments */
1409  --up; --low; --ipivot; --p; --x;
1410 
1411  /* Function Body */
1412  *spe = *stepmx / *pe;
1413  /* SPE IS THE STANDARD (UNCONSTRAINED) MAX STEP */
1414  i__1 = *n;
1415 
1416  for (i__ = 1; i__ <= i__1; ++i__)
1417  {
1418  if (ipivot[i__] != 0) goto L10;
1419 
1420  if (p[i__] == 0.) goto L10;
1421 
1422  if (p[i__] > 0.) goto L5;
1423 
1424  t = low[i__] - x[i__];
1425 
1426  if (t > *spe * p[i__]) *spe = t / p[i__];
1427 
1428  goto L10;
1429 L5:
1430  t = up[i__] - x[i__];
1431 
1432  if (t < *spe * p[i__]) *spe = t / p[i__];
1433 
1434 L10:;
1435  }
1436 
1437  return 0;
1438 } /* stpmax_ */
1439 
1440 /* Subroutine */ int modz_(C_INT *n, C_FLOAT64 *x, C_FLOAT64 *p, C_INT *
1441  ipivot, C_FLOAT64 *epsmch, C_FLOAT64 *low, C_FLOAT64 *up,
1442  C_FLOAT64 *flast, C_FLOAT64 *fnew)
1443 {
1444  /* System generated locals */
1445  C_INT i__1;
1446  C_FLOAT64 d__1;
1447 
1448  /* Local variables */
1449  C_INT i__;
1450  C_FLOAT64 tol;
1451 
1452  /* UPDATE THE CONSTRAINT MATRIX IF A NEW CONSTRAINT IS ENCOUNTERED */
1453 
1454  /* Parameter adjustments */
1455  --up;
1456  --low;
1457  --ipivot;
1458  --p;
1459  --x;
1460 
1461  /* Function Body */
1462  i__1 = *n;
1463 
1464  for (i__ = 1; i__ <= i__1; ++i__)
1465  {
1466  if (ipivot[i__] != 0)
1467  {
1468  goto L10;
1469  }
1470 
1471  if (p[i__] == 0.)
1472  {
1473  goto L10;
1474  }
1475 
1476  if (p[i__] > 0.)
1477  {
1478  goto L5;
1479  }
1480 
1481  tol = *epsmch * 10. * ((d__1 = low[i__], fabs(d__1)) + 1.);
1482 
1483  if (x[i__] - low[i__] > tol)
1484  {
1485  goto L10;
1486  }
1487 
1488  *flast = *fnew;
1489  ipivot[i__] = -1;
1490  x[i__] = low[i__];
1491  goto L10;
1492 L5:
1493  tol = *epsmch * 10. * ((d__1 = up[i__], fabs(d__1)) + 1.);
1494 
1495  if (up[i__] - x[i__] > tol)
1496  {
1497  goto L10;
1498  }
1499 
1500  *flast = *fnew;
1501  ipivot[i__] = 1;
1502  x[i__] = up[i__];
1503 L10:
1504  ;
1505  }
1506 
1507  return 0;
1508 } /* modz_ */
1509 
1510 /* Subroutine */ int cnvtst_(C_INT *conv, C_FLOAT64 *alpha, C_FLOAT64 *
1511  pnorm, C_FLOAT64 *toleps, C_FLOAT64 *xnorm, C_FLOAT64 *difnew,
1512  C_FLOAT64 *rtleps, C_FLOAT64 *ftest, C_FLOAT64 *gtg, C_FLOAT64 *
1513  peps, C_FLOAT64 * /* epsmch */, C_FLOAT64 *gtpnew, C_FLOAT64 *fnew,
1514  C_FLOAT64 *flast, C_FLOAT64 *g, C_INT *ipivot, C_INT *n,
1515  C_FLOAT64 *accrcy)
1516 {
1517  /* System generated locals */
1518  C_INT i__1;
1519 
1520  /* Local variables */
1521  C_FLOAT64 cmax;
1522  C_INT imax, i__;
1523  C_FLOAT64 t;
1524  C_INT ltest;
1525  C_FLOAT64 one;
1526 
1527  /* TEST FOR CONVERGENCE */
1528 
1529  /* Parameter adjustments */
1530  --ipivot;
1531  --g;
1532 
1533  /* Function Body */
1534  imax = 0;
1535  cmax = 0.;
1536  ltest = *flast - *fnew <= *gtpnew * -.5;
1537  i__1 = *n;
1538 
1539  for (i__ = 1; i__ <= i__1; ++i__)
1540  {
1541  if (ipivot[i__] == 0 || ipivot[i__] == 2)
1542  {
1543  goto L10;
1544  }
1545 
1546  t = -ipivot[i__] * g[i__];
1547 
1548  if (t >= 0.)
1549  {
1550  goto L10;
1551  }
1552 
1553  *conv = FALSE_;
1554 
1555  if (ltest)
1556  {
1557  goto L10;
1558  }
1559 
1560  if (cmax <= t)
1561  {
1562  goto L10;
1563  }
1564 
1565  cmax = t;
1566  imax = i__;
1567 L10:
1568  ;
1569  }
1570 
1571  if (imax == 0)
1572  {
1573  goto L15;
1574  }
1575 
1576  ipivot[imax] = 0;
1577  *flast = *fnew;
1578  return 0;
1579 L15:
1580  *conv = FALSE_;
1581  one = 1.;
1582 
1583  if ((*alpha * *pnorm >= *toleps * (one + *xnorm) || fabs(*difnew) >= *
1584  rtleps * *ftest || *gtg >= *peps * *ftest * *ftest) && *gtg >= *
1585  accrcy * 1e-4 * *ftest * *ftest)
1586  {
1587  return 0;
1588  }
1589 
1590  *conv = TRUE_;
1591 
1592  /* FOR DETAILS, SEE GILL, MURRAY, AND WRIGHT (1981, P. 308) AND */
1593  /* FLETCHER (1981, P. 116). THE MULTIPLIER TESTS (HERE, TESTING */
1594  /* THE SIGN OF THE COMPONENTS OF THE GRADIENT) MAY STILL NEED TO */
1595  /* MODIFIED TO INCORPORATE TOLERANCES FOR ZERO. */
1596 
1597  return 0;
1598 } /* cnvtst_ */
1599 
1600 /* Subroutine */ int crash_(C_INT *n, C_FLOAT64 *x, C_INT *ipivot,
1601  C_FLOAT64 *low, C_FLOAT64 *up, C_INT *ier)
1602 {
1603  /* System generated locals */
1604  C_INT i__1;
1605 
1606  /* Local variables */
1607  C_INT i__;
1608 
1609  /* THIS INITIALIZES THE CONSTRAINT INFORMATION, AND ENSURES THAT THE */
1610  /* INITIAL POINT SATISFIES LOW <= X <= UP. */
1611  /* THE CONSTRAINTS ARE CHECKED FOR CONSISTENCY. */
1612 
1613  /* Parameter adjustments */
1614  --up;
1615  --low;
1616  --ipivot;
1617  --x;
1618 
1619  /* Function Body */
1620  *ier = 0;
1621  i__1 = *n;
1622 
1623  for (i__ = 1; i__ <= i__1; ++i__)
1624  {
1625  if (x[i__] < low[i__])
1626  {
1627  x[i__] = low[i__];
1628  }
1629 
1630  if (x[i__] > up[i__])
1631  {
1632  x[i__] = up[i__];
1633  }
1634 
1635  ipivot[i__] = 0;
1636 
1637  if (x[i__] == low[i__])
1638  {
1639  ipivot[i__] = -1;
1640  }
1641 
1642  if (x[i__] == up[i__])
1643  {
1644  ipivot[i__] = 1;
1645  }
1646 
1647  if (up[i__] == low[i__])
1648  {
1649  ipivot[i__] = 2;
1650  }
1651 
1652  if (low[i__] > up[i__])
1653  {
1654  *ier = -i__;
1655  }
1656 
1657  /* L30: */
1658  }
1659 
1660  return 0;
1661 } /* crash_ */
1662 
1663 /* THE VECTORS SK AND YK, ALTHOUGH NOT IN THE CALL, */
1664 /* ARE USED (VIA THEIR POSITION IN W) BY THE ROUTINE MSOLVE. */
1665 
1666 /* Subroutine */ int CTruncatedNewton::modlnp_(C_INT *modet, C_FLOAT64 *zsol, C_FLOAT64 *gv,
1667  C_FLOAT64 *r__, C_FLOAT64 *v, C_FLOAT64 *diagb, C_FLOAT64 *emat,
1668  C_FLOAT64 *x, C_FLOAT64 *g, C_FLOAT64 *zk, C_INT *n, C_FLOAT64 *
1669  w, C_INT *lw, C_INT * /* niter */, C_INT *maxit, C_INT *nfeval,
1670  C_INT * /* nmodif */, C_INT *nlincg, C_INT *upd1, C_FLOAT64 *yksk,
1671  C_FLOAT64 *gsk, C_FLOAT64 *yrsr, C_INT *lreset, FTruncatedNewton *sfun,
1672  C_INT *bounds, C_INT *ipivot, C_FLOAT64 *accrcy, C_FLOAT64 *gtp,
1673  C_FLOAT64 *gnorm, C_FLOAT64 *xnorm)
1674 {
1675  /* Format strings */
1676  //remove the format
1677  /* char fmt_800[] = "(\002 \002,//,\002 ENTERING MODLNP\002)";
1678  char fmt_810[] = "(\002 \002,//,\002 ### ITERATION \002,i2,\002 #"
1679  "##\002)";
1680  char fmt_820[] = "(\002 ALPHA\002,1pd16.8)";
1681  char fmt_830[] = "(\002 G(T)Z POSITIVE AT ITERATION \002,i2,\002 "
1682  "- TRUNCATING METHOD\002,/)";
1683  char fmt_840[] = "(\002 \002,10x,\002HESSIAN NOT POSITIVE-DEFIN"
1684  "ITE\002)";
1685  char fmt_850[] = "(\002 \002,/,8x,\002MODLAN TRUNCATED AFTER \002"
1686  ",i3,\002 ITERATIONS\002,\002 RNORM = \002,1pd14.6)";
1687  char fmt_860[] = "(\002 PRECONDITIONING NOT POSITIVE-DEFINITE\002)"
1688  ;*/
1689 
1690  /* System generated locals */
1691  C_INT i__1, i__2;
1692  C_FLOAT64 d__1;
1693 
1694  /* Builtin functions */
1695  // C_INT s_wsfe(cilist *), e_wsfe(void), do_fio(C_INT *, char *, ftnlen);
1696 
1697  /* Local variables */
1698  C_FLOAT64 beta = 0.0;
1699  C_FLOAT64 qold, qnew;
1700  C_INT i__, k;
1701  C_FLOAT64 alpha, delta;
1702  C_INT first;
1703  C_FLOAT64 rzold = 0.0, qtest, pr;
1704  C_FLOAT64 rz;
1705  C_FLOAT64 rhsnrm, tol, vgv;
1706 
1707  /* THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT */
1708  /* ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH */
1709  /* DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM. WHEN THE VALUE OF THE */
1710  /* QUADRATIC MODEL IS SUFFICIENTLY REDUCED, */
1711  /* THE ITERATION IS TERMINATED. */
1712 
1713  /* PARAMETERS */
1714 
1715  /* MODET - C_INT WHICH CONTROLS AMOUNT OF OUTPUT */
1716  /* ZSOL - COMPUTED SEARCH DIRECTION */
1717  /* G - CURRENT GRADIENT */
1718  /* GV,GZ1,V - SCRATCH VECTORS */
1719  /* R - RESIDUAL */
1720  /* DIAGB,EMAT - DIAGONAL PRECONDITONING MATRIX */
1721  /* NITER - NONLINEAR ITERATION # */
1722  /* FEVAL - VALUE OF QUADRATIC FUNCTION */
1723 
1724  /* ************************************************************* */
1725  /* INITIALIZATION */
1726  /* ************************************************************* */
1727 
1728  /* GENERAL INITIALIZATION */
1729 
1730  /* Parameter adjustments */
1731  --zk;
1732  --g;
1733  --x;
1734  --emat;
1735  --diagb;
1736  --v;
1737  --r__;
1738  --gv;
1739  --zsol;
1740  --w;
1741  --ipivot;
1742 
1743  /* Function Body */
1744  if (*modet > 0)
1745  {
1746  //remove
1747  /*
1748  s_wsfe(&io___167);
1749  e_wsfe();*/
1750  }
1751 
1752  if (*maxit == 0)
1753  {
1754  return 0;
1755  }
1756 
1757  first = TRUE_;
1758  rhsnrm = *gnorm;
1759  tol = 1e-12;
1760  qold = 0.;
1761 
1762  /* INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM */
1763 
1764  CTruncatedNewton::initpc_(&diagb[1], &emat[1], n, &w[1], lw, modet, upd1, yksk, gsk, yrsr,
1765  lreset);
1766  i__1 = *n;
1767 
1768  for (i__ = 1; i__ <= i__1; ++i__)
1769  {
1770  r__[i__] = -g[i__];
1771  v[i__] = 0.;
1772  zsol[i__] = 0.;
1773  /* L10: */
1774  }
1775 
1776  /* ************************************************************ */
1777  /* MAIN ITERATION */
1778  /* ************************************************************ */
1779 
1780  i__1 = *maxit;
1781 
1782  for (k = 1; k <= i__1; ++k)
1783  {
1784  ++(*nlincg);
1785 
1786  if (*modet > 1)
1787  {
1788  /*
1789  s_wsfe(&io___174);
1790  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1791  e_wsfe();*/
1792  }
1793 
1794  /* CG ITERATION TO SOLVE SYSTEM OF EQUATIONS */
1795 
1796  if (*bounds)
1797  {
1798  ztime_(n, &r__[1], &ipivot[1]);
1799  }
1800 
1801  CTruncatedNewton::msolve_(&r__[1], &zk[1], n, &w[1], lw, upd1, yksk, gsk, yrsr, lreset,
1802  &first);
1803 
1804  if (*bounds)
1805  {
1806  ztime_(n, &zk[1], &ipivot[1]);
1807  }
1808 
1809  rz = ddot_(n, &r__[1], &c__1, &zk[1], &c__1);
1810 
1811  if (rz / rhsnrm < tol)
1812  {
1813  goto L80;
1814  }
1815 
1816  if (k == 1)
1817  {
1818  beta = 0.;
1819  }
1820 
1821  if (k > 1)
1822  {
1823  beta = rz / rzold;
1824  }
1825 
1826  i__2 = *n;
1827 
1828  for (i__ = 1; i__ <= i__2; ++i__)
1829  {
1830  v[i__] = zk[i__] + beta * v[i__];
1831  /* L20: */
1832  }
1833 
1834  if (*bounds)
1835  {
1836  ztime_(n, &v[1], &ipivot[1]);
1837  }
1838 
1839  CTruncatedNewton::gtims_(&v[1], &gv[1], n, &x[1], &g[1], &w[1], lw, sfun, &first,
1840  &delta, accrcy, xnorm);
1841 
1842  if (*bounds)
1843  {
1844  ztime_(n, &gv[1], &ipivot[1]);
1845  }
1846 
1847  ++(*nfeval);
1848  vgv = ddot_(n, &v[1], &c__1, &gv[1], &c__1);
1849 
1850  if (vgv / rhsnrm < tol)
1851  {
1852  goto L50;
1853  }
1854 
1855  ndia3_(n, &emat[1], &v[1], &gv[1], &r__[1], &vgv, modet);
1856 
1857  /* COMPUTE LINEAR STEP LENGTH */
1858 
1859  alpha = rz / vgv;
1860 
1861  if (*modet >= 1)
1862  {
1863  /*
1864  s_wsfe(&io___181);
1865  do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(C_FLOAT64));
1866  e_wsfe();*/
1867  }
1868 
1869  /* COMPUTE CURRENT SOLUTION AND RELATED VECTORS */
1870 
1871  daxpy_(n, &alpha, &v[1], &c__1, &zsol[1], &c__1);
1872  d__1 = -alpha;
1873  daxpy_(n, &d__1, &gv[1], &c__1, &r__[1], &c__1);
1874 
1875  /* TEST FOR CONVERGENCE */
1876 
1877  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1878  pr = ddot_(n, &r__[1], &c__1, &zsol[1], &c__1);
1879  qnew = (*gtp + pr) * .5;
1880  qtest = k * (1. - qold / qnew);
1881 
1882  if (qtest < 0.)
1883  {
1884  goto L70;
1885  }
1886 
1887  qold = qnew;
1888 
1889  if (qtest <= .5)
1890  {
1891  goto L70;
1892  }
1893 
1894  /* PERFORM CAUTIONARY TEST */
1895 
1896  if (*gtp > 0.)
1897  {
1898  goto L40;
1899  }
1900 
1901  rzold = rz;
1902  /* L30: */
1903  }
1904 
1905  /* TERMINATE ALGORITHM */
1906 
1907  --k;
1908  goto L70;
1909 
1910  /* TRUNCATE ALGORITHM IN CASE OF AN EMERGENCY */
1911 
1912 L40:
1913 
1914  if (*modet >= -1)
1915  {
1916  /*
1917  s_wsfe(&io___185);
1918  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1919  e_wsfe();*/
1920  }
1921 
1922  d__1 = -alpha;
1923  daxpy_(n, &d__1, &v[1], &c__1, &zsol[1], &c__1);
1924  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1925  goto L90;
1926 L50:
1927 
1928  if (*modet > -2)
1929  {
1930  /*
1931  s_wsfe(&io___186);
1932  e_wsfe();*/
1933  }
1934 
1935  /* L60: */
1936  if (k > 1)
1937  {
1938  goto L70;
1939  }
1940 
1941  CTruncatedNewton::msolve_(&g[1], &zsol[1], n, &w[1], lw, upd1, yksk, gsk, yrsr, lreset, &
1942  first);
1943  negvec_(n, &zsol[1]);
1944 
1945  if (*bounds)
1946  {
1947  ztime_(n, &zsol[1], &ipivot[1]);
1948  }
1949 
1950  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1951 L70:
1952 
1953  if (*modet >= -1)
1954  {
1955  /*
1956  s_wsfe(&io___187);
1957  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1958  do_fio(&c__1, (char *)&rnorm, (ftnlen)sizeof(C_FLOAT64));
1959  e_wsfe();*/
1960  }
1961 
1962  goto L90;
1963 L80:
1964 
1965  if (*modet >= -1)
1966  {
1967  /*
1968  s_wsfe(&io___189);
1969  e_wsfe();*/
1970  }
1971 
1972  if (k > 1)
1973  {
1974  goto L70;
1975  }
1976 
1977  dcopy_(n, &g[1], &c__1, &zsol[1], &c__1);
1978  negvec_(n, &zsol[1]);
1979 
1980  if (*bounds)
1981  {
1982  ztime_(n, &zsol[1], &ipivot[1]);
1983  }
1984 
1985  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1986  goto L70;
1987 
1988  /* STORE (OR RESTORE) DIAGONAL PRECONDITIONING */
1989 
1990 L90:
1991  dcopy_(n, &emat[1], &c__1, &diagb[1], &c__1);
1992  return 0;
1993 } /* modlnp_ */
1994 
1995 /* Subroutine */ int ndia3_(C_INT *n, C_FLOAT64 *e, C_FLOAT64 *v,
1996  C_FLOAT64 *gv, C_FLOAT64 *r__, C_FLOAT64 *vgv, C_INT *modet)
1997 {
1998  /* Format strings */
1999  // char fmt_800[] = "(\002 *** EMAT NEGATIVE: \002,1pd16.8)";
2000 
2001  /* System generated locals */
2002  C_INT i__1;
2003 
2004  /* Builtin functions */
2005  // C_INT s_wsfe(cilist *), do_fio(C_INT *, char *, ftnlen), e_wsfe(void);
2006 
2007  /* Local variables */
2008  C_INT i__;
2009  C_FLOAT64 vr;
2010 
2011  /* Fortran I/O blocks */
2012  // cilist io___192 = {0, 6, 0, fmt_800, 0 };
2013 
2014  /* UPDATE THE PRECONDITIOING MATRIX BASED ON A DIAGONAL VERSION */
2015  /* OF THE BFGS QUASI-NEWTON UPDATE. */
2016 
2017  /* Parameter adjustments */
2018  --r__;
2019  --gv;
2020  --v;
2021  --e;
2022 
2023  /* Function Body */
2024  vr = ddot_(n, &v[1], &c__1, &r__[1], &c__1);
2025  i__1 = *n;
2026 
2027  for (i__ = 1; i__ <= i__1; ++i__)
2028  {
2029  e[i__] = e[i__] - r__[i__] * r__[i__] / vr + gv[i__] * gv[i__] / *vgv;
2030 
2031  if (e[i__] > 1e-6)
2032  {
2033  goto L10;
2034  }
2035 
2036  if (*modet > -2)
2037  {
2038  /*
2039  s_wsfe(&io___192);
2040  do_fio(&c__1, (char *)&e[i__], (ftnlen)sizeof(C_FLOAT64));
2041  e_wsfe();*/
2042  }
2043 
2044  e[i__] = 1.;
2045 L10:
2046  ;
2047  }
2048 
2049  return 0;
2050 } /* ndia3_ */
2051 
2052 /* SERVICE ROUTINES FOR OPTIMIZATION */
2053 
2054 /* Subroutine */ int negvec_(C_INT *n, C_FLOAT64 *v)
2055 {
2056  /* System generated locals */
2057  C_INT i__1;
2058 
2059  /* Local variables */
2060  C_INT i__;
2061 
2062  /* NEGATIVE OF THE VECTOR V */
2063 
2064  /* Parameter adjustments */
2065  --v;
2066 
2067  /* Function Body */
2068  i__1 = *n;
2069 
2070  for (i__ = 1; i__ <= i__1; ++i__)
2071  {
2072  v[i__] = -v[i__];
2073  /* L10: */
2074  }
2075 
2076  return 0;
2077 } /* negvec_ */
2078 
2079 /* Subroutine */ int lsout_(C_INT * /* iloc */, C_INT * /* itest */, C_FLOAT64 *xmin,
2080  C_FLOAT64 * /* fmin */, C_FLOAT64 * /* gmin */, C_FLOAT64 *xw, C_FLOAT64 * /* fw */,
2081  C_FLOAT64 * /* gw */, C_FLOAT64 *u, C_FLOAT64 *a, C_FLOAT64 *b,
2082  C_FLOAT64 * /* tol */, C_FLOAT64 * /* eps */, C_FLOAT64 *scxbd, C_FLOAT64 *
2083  /* xlamda */)
2084 {
2085  /* Format strings */
2086  /* char fmt_800[] = "(///\002 OUTPUT FROM LINEAR SEARCH\002)";
2087  char fmt_810[] = "(\002 TOL AND EPS\002/2d25.14)";
2088  char fmt_820[] = "(\002 CURRENT UPPER AND LOWER BOUNDS\002/2d25."
2089  "14)";
2090  char fmt_830[] = "(\002 STRICT UPPER BOUND\002/d25.14)";
2091  char fmt_840[] = "(\002 XW, FW, GW\002/3d25.14)";
2092  char fmt_850[] = "(\002 XMIN, FMIN, GMIN\002/3d25.14)";
2093  char fmt_860[] = "(\002 NEW ESTIMATE\002/2d25.14)";
2094  char fmt_870[] = "(\002 ILOC AND ITEST\002/2i3)"; */
2095 
2096  /* Builtin functions */
2097  // C_INT s_wsfe(cilist *), e_wsfe(void), do_fio(C_INT *, char *, ftnlen);
2098 
2099  /* Local variables */
2100  C_FLOAT64 ybnd, ya, yb, yu, yw;
2101 
2102  /* ERROR PRINTOUTS FOR GETPTC */
2103 
2104  yu = *xmin + *u;
2105  ya = *a + *xmin;
2106  yb = *b + *xmin;
2107  yw = *xw + *xmin;
2108  ybnd = *scxbd + *xmin;
2109  //remove the printing
2110  /* s_wsfe(&io___199);
2111  e_wsfe();
2112  s_wsfe(&io___200);
2113  do_fio(&c__1, (char *)&(*tol), (ftnlen)sizeof(C_FLOAT64));
2114  do_fio(&c__1, (char *)&(*eps), (ftnlen)sizeof(C_FLOAT64));
2115  e_wsfe();
2116  s_wsfe(&io___201);
2117  do_fio(&c__1, (char *)&ya, (ftnlen)sizeof(C_FLOAT64));
2118  do_fio(&c__1, (char *)&yb, (ftnlen)sizeof(C_FLOAT64));
2119  e_wsfe();
2120  s_wsfe(&io___202);
2121  do_fio(&c__1, (char *)&ybnd, (ftnlen)sizeof(C_FLOAT64));
2122  e_wsfe();
2123  s_wsfe(&io___203);
2124  do_fio(&c__1, (char *)&yw, (ftnlen)sizeof(C_FLOAT64));
2125  do_fio(&c__1, (char *)&(*fw), (ftnlen)sizeof(C_FLOAT64));
2126  do_fio(&c__1, (char *)&(*gw), (ftnlen)sizeof(C_FLOAT64));
2127  e_wsfe();
2128  s_wsfe(&io___204);
2129  do_fio(&c__1, (char *)&(*xmin), (ftnlen)sizeof(C_FLOAT64));
2130  do_fio(&c__1, (char *)&(*fmin), (ftnlen)sizeof(C_FLOAT64));
2131  do_fio(&c__1, (char *)&(*gmin), (ftnlen)sizeof(C_FLOAT64));
2132  e_wsfe();
2133  s_wsfe(&io___205);
2134  do_fio(&c__1, (char *)&yu, (ftnlen)sizeof(C_FLOAT64));
2135  e_wsfe();
2136  s_wsfe(&io___206);
2137  do_fio(&c__1, (char *)&(*iloc), (ftnlen)sizeof(C_INT));
2138  do_fio(&c__1, (char *)&(*itest), (ftnlen)sizeof(C_INT));
2139  e_wsfe(); */
2140  return 0;
2141 } /* lsout_ */
2142 
2144  C_FLOAT64 *smax)
2145 {
2146  /* System generated locals */
2147  C_FLOAT64 ret_val, d__1;
2148 
2149  /* Local variables */
2150  C_FLOAT64 d__, alpha;
2151  C_FLOAT64 epsmch;
2152 
2153  /* ******************************************************** */
2154  /* STEP1 RETURNS THE LENGTH OF THE INITIAL STEP TO BE TAKEN ALONG THE */
2155  /* VECTOR P IN THE NEXT LINEAR SEARCH. */
2156  /* ******************************************************** */
2157 
2158  epsmch = mchpr1_();
2159  d__ = (d__1 = *fnew - *fm, fabs(d__1));
2160  alpha = 1.;
2161 
2162  if (d__ * 2. <= -(*gtp) && d__ >= epsmch)
2163  {
2164  alpha = d__ * -2. / *gtp;
2165  }
2166 
2167  if (alpha >= *smax)
2168  {
2169  alpha = *smax;
2170  }
2171 
2172  ret_val = alpha;
2173  return ret_val;
2174 } /* step1_ */
2175 
2177 {
2178  /* System generated locals */
2179  C_FLOAT64 ret_val;
2180 
2181  /* RETURNS THE VALUE OF EPSMCH, WHERE EPSMCH IS THE SMALLEST POSSIBLE */
2182  /* REAL NUMBER SUCH THAT 1.0 + EPSMCH .GT. 1.0 */
2183 
2184  /* FOR IEEE */
2185 
2186  ret_val = std::numeric_limits< C_FLOAT64 >::epsilon();
2187 
2188  return ret_val;
2189 } /* mchpr1_ */
2190 
2191 /* Subroutine */ int chkucp_(C_INT *lwtest, C_INT *maxfun, C_INT *nwhy,
2192  C_INT *n, C_FLOAT64 *alpha, C_FLOAT64 *epsmch, C_FLOAT64 *eta,
2193  C_FLOAT64 *peps, C_FLOAT64 *rteps, C_FLOAT64 *rtol, C_FLOAT64 *
2194  rtolsq, C_FLOAT64 *stepmx, C_FLOAT64 *test, C_FLOAT64 *xtol,
2195  C_FLOAT64 *xnorm, C_FLOAT64 *x, C_INT *lw, C_FLOAT64 *small,
2196  C_FLOAT64 *tiny, C_FLOAT64 *accrcy)
2197 {
2198  /* Local variables */
2199  /* CHECKS PARAMETERS AND SETS CONSTANTS WHICH ARE COMMON TO BOTH */
2200  /* DERIVATIVE AND NON-DERIVATIVE ALGORITHMS */
2201 
2202  /* Parameter adjustments */
2203  --x;
2204 
2205  /* Function Body */
2206  *epsmch = mchpr1_();
2207  *small = *epsmch * *epsmch;
2208  *tiny = *small;
2209  *nwhy = -1;
2210  *rteps = sqrt(*epsmch);
2211  *rtol = *xtol;
2212 
2213  if (fabs(*rtol) < *accrcy)
2214  {
2215  *rtol = *rteps * 10.;
2216  }
2217 
2218  /* CHECK FOR ERRORS IN THE INPUT PARAMETERS */
2219 
2220  if (*lw < *lwtest || *n < 1 || *rtol < 0. || *eta >= 1. || *eta < 0. || *
2221  stepmx < *rtol || *maxfun < 1)
2222  {
2223  return 0;
2224  }
2225 
2226  *nwhy = 0;
2227 
2228  /* SET CONSTANTS FOR LATER */
2229 
2230  *rtolsq = *rtol * *rtol;
2231  *peps = pow(*accrcy, c_b246);
2232  *xnorm = dnrm2_(n, &x[1], &c__1);
2233  *alpha = 0.;
2234  *test = 0.;
2235  return 0;
2236 } /* chkucp_ */
2237 
2238 /* Subroutine */ int CTruncatedNewton::setucr_(C_FLOAT64 * /* small */, C_INT *nftotl, C_INT *
2239  niter, C_INT *n, C_FLOAT64 *f, C_FLOAT64 *fnew, C_FLOAT64 *fm,
2240  C_FLOAT64 *gtg, C_FLOAT64 *oldf, FTruncatedNewton *sfun, C_FLOAT64 *g,
2241  C_FLOAT64 *x)
2242 {
2243  /* CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET */
2244  /* CONSTANTS FOR THE SUBSEQUENT MINIMIZATION */
2245 
2246  /* Parameter adjustments */
2247  --x;
2248  --g;
2249 
2250  /* Function Body */
2251  *fm = *f;
2252 
2253  /* COMPUTE THE INITIAL FUNCTION VALUE */
2254 
2255  (*sfun)(n, &x[1], fnew, &g[1]);
2256  *nftotl = 1;
2257 
2258  /* SET CONSTANTS FOR LATER */
2259 
2260  *niter = 0;
2261  *oldf = *fnew;
2262  *gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
2263  return 0;
2264 } /* setucr_ */
2265 
2266 /* Subroutine */ int CTruncatedNewton::gtims_(C_FLOAT64 *v, C_FLOAT64 *gv, C_INT *n,
2267  C_FLOAT64 *x, C_FLOAT64 *g, C_FLOAT64 *w, C_INT * /* lw */, FTruncatedNewton *sfun,
2268  C_INT *first, C_FLOAT64 *delta, C_FLOAT64 *accrcy, C_FLOAT64 *
2269  xnorm)
2270 {
2271  /* System generated locals */
2272  C_INT i__1;
2273 
2274  /* Local variables */
2275  C_FLOAT64 dinv, f;
2276  C_INT i__, ihg;
2277 
2278  /* THIS ROUTINE COMPUTES THE PRODUCT OF THE MATRIX G TIMES THE VECTOR */
2279  /* V AND STORES THE RESULT IN THE VECTOR GV (FINITE-DIFFERENCE VERSION) */
2280 
2281  /* Parameter adjustments */
2282  --g;
2283  --x;
2284  --gv;
2285  --v;
2286  --w;
2287 
2288  /* Function Body */
2289  if (!(*first))
2290  {
2291  goto L20;
2292  }
2293 
2294  *delta = sqrt(*accrcy) * (*xnorm + 1.);
2295  *first = FALSE_;
2296 L20:
2297  dinv = 1. / *delta;
2298  ihg = subscr_2.lhg;
2299  i__1 = *n;
2300 
2301  for (i__ = 1; i__ <= i__1; ++i__)
2302  {
2303  w[ihg] = x[i__] + *delta * v[i__];
2304  ++ihg;
2305  /* L30: */
2306  }
2307 
2308  (*sfun)(n, &w[subscr_2.lhg], &f, &gv[1]);
2309  i__1 = *n;
2310 
2311  for (i__ = 1; i__ <= i__1; ++i__)
2312  {
2313  gv[i__] = (gv[i__] - g[i__]) * dinv;
2314  /* L40: */
2315  }
2316 
2317  return 0;
2318 } /* gtims_ */
2319 
2320 /* Subroutine */ int CTruncatedNewton::msolve_(C_FLOAT64 *g, C_FLOAT64 *y, C_INT *n,
2321  C_FLOAT64 *w, C_INT * /* lw */, C_INT *upd1, C_FLOAT64 *yksk,
2322  C_FLOAT64 *gsk, C_FLOAT64 *yrsr, C_INT *lreset, C_INT *first)
2323 {
2324  /* THIS ROUTINE SETS UPT THE ARRAYS FOR MSLV */
2325 
2326  /* Parameter adjustments */
2327  --y;
2328  --g;
2329  --w;
2330 
2331  /* Function Body */
2332  mslv_(&g[1], &y[1], n, &w[subscr_2.lsk], &w[subscr_2.lyk], &w[
2333  subscr_2.ldiagb], &w[subscr_2.lsr], &w[subscr_2.lyr], &w[
2334  subscr_2.lhyr], &w[subscr_2.lhg], &w[subscr_2.lhyk], upd1, yksk,
2335  gsk, yrsr, lreset, first);
2336  return 0;
2337 } /* msolve_ */
2338 
2339 /* Subroutine */ int mslv_(C_FLOAT64 *g, C_FLOAT64 *y, C_INT *n,
2340  C_FLOAT64 *sk, C_FLOAT64 *yk, C_FLOAT64 *diagb, C_FLOAT64 *sr,
2341  C_FLOAT64 *yr, C_FLOAT64 *hyr, C_FLOAT64 *hg, C_FLOAT64 *hyk,
2342  C_INT *upd1, C_FLOAT64 *yksk, C_FLOAT64 *gsk, C_FLOAT64 *yrsr,
2343  C_INT *lreset, C_INT *first)
2344 {
2345  /* System generated locals */
2346  C_INT i__1;
2347 
2348  /* Local variables */
2349  C_FLOAT64 ghyk, ghyr, yksr = 0.0;
2350  C_INT i__;
2351  C_FLOAT64 ykhyk = 0.0, ykhyr = 0.0, yrhyr = 0.0, rdiagb;
2352  C_FLOAT64 one, gsr;
2353 
2354  /* THIS ROUTINE ACTS AS A PRECONDITIONING STEP FOR THE */
2355  /* LINEAR CONJUGATE-GRADIENT ROUTINE. IT IS ALSO THE */
2356  /* METHOD OF COMPUTING THE SEARCH DIRECTION FROM THE */
2357  /* GRADIENT FOR THE NON-LINEAR CONJUGATE-GRADIENT CODE. */
2358  /* IT REPRESENTS A TWO-STEP SELF-SCALED BFGS FORMULA. */
2359 
2360  /* Parameter adjustments */
2361  --hyk;
2362  --hg;
2363  --hyr;
2364  --yr;
2365  --sr;
2366  --diagb;
2367  --yk;
2368  --sk;
2369  --y;
2370  --g;
2371 
2372  /* Function Body */
2373  if (*upd1)
2374  {
2375  goto L100;
2376  }
2377 
2378  one = 1.;
2379  *gsk = ddot_(n, &g[1], &c__1, &sk[1], &c__1);
2380 
2381  if (*lreset)
2382  {
2383  goto L60;
2384  }
2385 
2386  /* COMPUTE HG AND HY WHERE H IS THE INVERSE OF THE DIAGONALS */
2387 
2388  i__1 = *n;
2389 
2390  for (i__ = 1; i__ <= i__1; ++i__)
2391  {
2392  rdiagb = 1. / diagb[i__];
2393  hg[i__] = g[i__] * rdiagb;
2394 
2395  if (*first)
2396  {
2397  hyk[i__] = yk[i__] * rdiagb;
2398  }
2399 
2400  if (*first)
2401  {
2402  hyr[i__] = yr[i__] * rdiagb;
2403  }
2404 
2405  /* L57: */
2406  }
2407 
2408  if (*first)
2409  {
2410  yksr = ddot_(n, &yk[1], &c__1, &sr[1], &c__1);
2411  }
2412 
2413  if (*first)
2414  {
2415  ykhyr = ddot_(n, &yk[1], &c__1, &hyr[1], &c__1);
2416  }
2417 
2418  gsr = ddot_(n, &g[1], &c__1, &sr[1], &c__1);
2419  ghyr = ddot_(n, &g[1], &c__1, &hyr[1], &c__1);
2420 
2421  if (*first)
2422  {
2423  yrhyr = ddot_(n, &yr[1], &c__1, &hyr[1], &c__1);
2424  }
2425 
2426  ssbfgs_(n, &one, &sr[1], &yr[1], &hg[1], &hyr[1], yrsr, &yrhyr, &gsr, &
2427  ghyr, &hg[1]);
2428 
2429  if (*first)
2430  {
2431  ssbfgs_(n, &one, &sr[1], &yr[1], &hyk[1], &hyr[1], yrsr, &yrhyr, &
2432  yksr, &ykhyr, &hyk[1]);
2433  }
2434 
2435  ykhyk = ddot_(n, &hyk[1], &c__1, &yk[1], &c__1);
2436  ghyk = ddot_(n, &hyk[1], &c__1, &g[1], &c__1);
2437  ssbfgs_(n, &one, &sk[1], &yk[1], &hg[1], &hyk[1], yksk, &ykhyk, gsk, &
2438  ghyk, &y[1]);
2439  return 0;
2440 L60:
2441 
2442  /* COMPUTE GH AND HY WHERE H IS THE INVERSE OF THE DIAGONALS */
2443 
2444  i__1 = *n;
2445 
2446  for (i__ = 1; i__ <= i__1; ++i__)
2447  {
2448  rdiagb = 1. / diagb[i__];
2449  hg[i__] = g[i__] * rdiagb;
2450 
2451  if (*first)
2452  {
2453  hyk[i__] = yk[i__] * rdiagb;
2454  }
2455 
2456  /* L65: */
2457  }
2458 
2459  if (*first)
2460  {
2461  ykhyk = ddot_(n, &yk[1], &c__1, &hyk[1], &c__1);
2462  }
2463 
2464  ghyk = ddot_(n, &g[1], &c__1, &hyk[1], &c__1);
2465  ssbfgs_(n, &one, &sk[1], &yk[1], &hg[1], &hyk[1], yksk, &ykhyk, gsk, &
2466  ghyk, &y[1]);
2467  return 0;
2468 L100:
2469  i__1 = *n;
2470 
2471  for (i__ = 1; i__ <= i__1; ++i__)
2472  {
2473  /* L110: */
2474  y[i__] = g[i__] / diagb[i__];
2475  }
2476 
2477  return 0;
2478 } /* mslv_ */
2479 
2480 /* Subroutine */ int ssbfgs_(C_INT *n, C_FLOAT64 *gamma, C_FLOAT64 *sj,
2481  C_FLOAT64 *yj, C_FLOAT64 *hjv, C_FLOAT64 *hjyj, C_FLOAT64 *yjsj,
2482  C_FLOAT64 *yjhyj, C_FLOAT64 *vsj, C_FLOAT64 *vhyj, C_FLOAT64 *
2483  hjp1v)
2484 {
2485  /* System generated locals */
2486  C_INT i__1;
2487 
2488  /* Local variables */
2489  C_FLOAT64 beta;
2490  C_INT i__;
2491  C_FLOAT64 delta;
2492 
2493  /* SELF-SCALED BFGS */
2494 
2495  /* Parameter adjustments */
2496  --hjp1v;
2497  --hjyj;
2498  --hjv;
2499  --yj;
2500  --sj;
2501 
2502  /* Function Body */
2503  delta = (*gamma * *yjhyj / *yjsj + 1.) * *vsj / *yjsj - *gamma * *vhyj / *
2504  yjsj;
2505  beta = -(*gamma) * *vsj / *yjsj;
2506  i__1 = *n;
2507 
2508  for (i__ = 1; i__ <= i__1; ++i__)
2509  {
2510  hjp1v[i__] = *gamma * hjv[i__] + delta * sj[i__] + beta * hjyj[i__];
2511  /* L10: */
2512  }
2513 
2514  return 0;
2515 } /* ssbfgs_ */
2516 
2517 /* ROUTINES TO INITIALIZE PRECONDITIONER */
2518 
2519 /* Subroutine */ int CTruncatedNewton::initpc_(C_FLOAT64 *diagb, C_FLOAT64 *emat, C_INT *n,
2520  C_FLOAT64 *w, C_INT * /* lw */, C_INT *modet, C_INT *upd1, C_FLOAT64
2521  *yksk, C_FLOAT64 * /* gsk */, C_FLOAT64 *yrsr, C_INT *lreset)
2522 {
2523  /* Parameter adjustments */
2524  --emat;
2525  --diagb;
2526  --w;
2527 
2528  /* Function Body */
2529  initp3_(&diagb[1], &emat[1], n, lreset, yksk, yrsr, &w[subscr_2.lhyk], &w[
2530  subscr_2.lsk], &w[subscr_2.lyk], &w[subscr_2.lsr], &w[
2531  subscr_2.lyr], modet, upd1);
2532  return 0;
2533 } /* initpc_ */
2534 
2535 /* Subroutine */ int initp3_(C_FLOAT64 *diagb, C_FLOAT64 *emat, C_INT *n,
2536  C_INT *lreset, C_FLOAT64 *yksk, C_FLOAT64 *yrsr, C_FLOAT64 *bsk,
2537  C_FLOAT64 *sk, C_FLOAT64 *yk, C_FLOAT64 *sr, C_FLOAT64 *yr,
2538  C_INT *modet, C_INT *upd1)
2539 {
2540  /* Format strings */
2541  // char fmt_800[] = "(\002 \002,//8x,\002DMIN =\002,1pd12.4,\002 DM"
2542  // "AX =\002,1pd12.4,\002 COND =\002,1pd12.4,/)";
2543 
2544  /* System generated locals */
2545  C_INT i__1;
2546 
2547  /* Builtin functions */
2548  // C_INT s_wsfe(cilist *), do_fio(C_INT *, char *, ftnlen), e_wsfe(void);
2549 
2550  /* Local variables */
2551  C_FLOAT64 cond;
2552  C_FLOAT64 srds, yrsk;
2553  C_INT i__;
2554  C_FLOAT64 d1, dn, td, sds;
2555 
2556  /* Fortran I/O blocks */
2557  // cilist io___235 = {0, 6, 0, fmt_800, 0 };
2558 
2559  /* Parameter adjustments */
2560  --yr;
2561  --sr;
2562  --yk;
2563  --sk;
2564  --bsk;
2565  --emat;
2566  --diagb;
2567 
2568  /* Function Body */
2569  if (*upd1)
2570  {
2571  goto L90;
2572  }
2573 
2574  if (*lreset)
2575  {
2576  goto L60;
2577  }
2578 
2579  i__1 = *n;
2580 
2581  for (i__ = 1; i__ <= i__1; ++i__)
2582  {
2583  bsk[i__] = diagb[i__] * sr[i__];
2584  /* L10: */
2585  }
2586 
2587  sds = ddot_(n, &sr[1], &c__1, &bsk[1], &c__1);
2588  srds = ddot_(n, &sk[1], &c__1, &bsk[1], &c__1);
2589  yrsk = ddot_(n, &yr[1], &c__1, &sk[1], &c__1);
2590  i__1 = *n;
2591 
2592  for (i__ = 1; i__ <= i__1; ++i__)
2593  {
2594  td = diagb[i__];
2595  bsk[i__] = td * sk[i__] - bsk[i__] * srds / sds + yr[i__] * yrsk / *
2596  yrsr;
2597  emat[i__] = td - td * td * sr[i__] * sr[i__] / sds + yr[i__] * yr[i__]
2598  / *yrsr;
2599  /* L20: */
2600  }
2601 
2602  sds = ddot_(n, &sk[1], &c__1, &bsk[1], &c__1);
2603  i__1 = *n;
2604 
2605  for (i__ = 1; i__ <= i__1; ++i__)
2606  {
2607  emat[i__] = emat[i__] - bsk[i__] * bsk[i__] / sds + yk[i__] * yk[i__]
2608  / *yksk;
2609  /* L30: */
2610  }
2611 
2612  goto L110;
2613 L60:
2614  i__1 = *n;
2615 
2616  for (i__ = 1; i__ <= i__1; ++i__)
2617  {
2618  bsk[i__] = diagb[i__] * sk[i__];
2619  /* L70: */
2620  }
2621 
2622  sds = ddot_(n, &sk[1], &c__1, &bsk[1], &c__1);
2623  i__1 = *n;
2624 
2625  for (i__ = 1; i__ <= i__1; ++i__)
2626  {
2627  td = diagb[i__];
2628  emat[i__] = td - td * td * sk[i__] * sk[i__] / sds + yk[i__] * yk[i__]
2629  / *yksk;
2630  /* L80: */
2631  }
2632 
2633  goto L110;
2634 L90:
2635  dcopy_(n, &diagb[1], &c__1, &emat[1], &c__1);
2636 L110:
2637 
2638  if (*modet < 1)
2639  {
2640  return 0;
2641  }
2642 
2643  d1 = emat[1];
2644  dn = emat[1];
2645  i__1 = *n;
2646 
2647  for (i__ = 1; i__ <= i__1; ++i__)
2648  {
2649  if (emat[i__] < d1)
2650  {
2651  d1 = emat[i__];
2652  }
2653 
2654  if (emat[i__] > dn)
2655  {
2656  dn = emat[i__];
2657  }
2658 
2659  /* L120: */
2660  }
2661 
2662  cond = dn / d1;
2663  /*s_wsfe(&io___235);
2664  do_fio(&c__1, (char *)&d1, (ftnlen)sizeof(C_FLOAT64));
2665  do_fio(&c__1, (char *)&dn, (ftnlen)sizeof(C_FLOAT64));
2666  do_fio(&c__1, (char *)&cond, (ftnlen)sizeof(C_FLOAT64));
2667  e_wsfe();*/
2668  return 0;
2669 } /* initp3_ */
2670 
2671 /* Subroutine */ int CTruncatedNewton::setpar_(C_INT *n)
2672 {
2673  C_INT i__;
2674 
2675  /* SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE */
2676 
2677  for (i__ = 1; i__ <= 14; ++i__)
2678  {
2679  subscr_3.lsub[i__ - 1] = (i__ - 1) * *n + 1;
2680  /* L10: */
2681  }
2682 
2683  subscr_3.lwtest = subscr_3.lsub[13] + *n - 1;
2684  return 0;
2685 } /* setpar_ */
2686 
2687 /* LINE SEARCH ALGORITHMS OF GILL AND MURRAY */
2688 
2689 /* Subroutine */ int CTruncatedNewton::linder_(C_INT *n, FTruncatedNewton *sfun, C_FLOAT64 *small,
2690  C_FLOAT64 *epsmch, C_FLOAT64 *reltol, C_FLOAT64 *fabstol,
2691  C_FLOAT64 *tnytol, C_FLOAT64 *eta, C_FLOAT64 * /* sftbnd */, C_FLOAT64 *
2692  xbnd, C_FLOAT64 *p, C_FLOAT64 *gtp, C_FLOAT64 *x, C_FLOAT64 *f,
2693  C_FLOAT64 *alpha, C_FLOAT64 *g, C_INT *nftotl, C_INT *iflag,
2694  C_FLOAT64 *w, C_INT * /* lw */)
2695 {
2696  /* System generated locals */
2697  C_INT i__1;
2698 
2699  /* Local variables */
2700  C_FLOAT64 oldf, fmin, gmin;
2701  C_INT numf;
2702  C_FLOAT64 step, xmin, a, b, e;
2703  C_INT i__, l;
2704  C_FLOAT64 u;
2705  C_INT itcnt;
2706  C_FLOAT64 b1;
2707  C_INT itest, nprnt;
2708  C_FLOAT64 gtest1, gtest2;
2709  C_INT lg;
2710  C_FLOAT64 fu, gu, fw, gw;
2711  C_INT lx;
2712  C_INT braktd;
2713  C_FLOAT64 ualpha, factor, scxbnd, xw;
2714  C_FLOAT64 fpresn;
2715  C_INT ientry;
2716  C_FLOAT64 rtsmll;
2717  C_INT lsprnt;
2718  C_FLOAT64 big, tol, rmu;
2719 
2720  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE */
2721  /* CALLED WITHIN LINDER */
2722 
2723  /* ALLOCATE THE ADDRESSES FOR LOCAL WORKSPACE */
2724 
2725  /* Parameter adjustments */
2726  --g;
2727  --x;
2728  --p;
2729  --w;
2730 
2731  /* Function Body */
2732  lx = 1;
2733  lg = lx + *n;
2734  lsprnt = 0;
2735  nprnt = 10000;
2736  rtsmll = sqrt(*small);
2737  big = 1. / *small;
2738  itcnt = 0;
2739 
2740  /* SET THE ESTIMATED RELATIVE PRECISION IN F(X). */
2741 
2742  fpresn = *epsmch * 10.;
2743  numf = 0;
2744  u = *alpha;
2745  fu = *f;
2746  fmin = *f;
2747  gu = *gtp;
2748  rmu = 1e-4;
2749 
2750  /* FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY. */
2751 
2752  ientry = 1;
2753 L10:
2754 
2755  /* TEST FOR TOO MANY ITERATIONS */
2756 
2757  ++itcnt;
2758  *iflag = 1;
2759 
2760  if (itcnt > 20)
2761  {
2762  goto L50;
2763  }
2764 
2765  *iflag = 0;
2766  CTruncatedNewton::getptc_(&big, small, &rtsmll, reltol, fabstol, tnytol, &fpresn, eta, &rmu,
2767  xbnd, &u, &fu, &gu, &xmin, &fmin, &gmin, &xw, &fw, &gw, &a, &b, &
2768  oldf, &b1, &scxbnd, &e, &step, &factor, &braktd, &gtest1, &gtest2,
2769  &tol, &ientry, &itest);
2770 
2771  /* LSOUT */
2772  if (lsprnt >= nprnt)
2773  {
2774  lsout_(&ientry, &itest, &xmin, &fmin, &gmin, &xw, &fw, &gw, &u, &a, &
2775  b, &tol, reltol, &scxbnd, xbnd);
2776  }
2777 
2778  /* IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE */
2779  /* CALCULATED. */
2780 
2781  if (itest != 1)
2782  {
2783  goto L30;
2784  }
2785 
2786  ualpha = xmin + u;
2787  l = lx;
2788  i__1 = *n;
2789 
2790  for (i__ = 1; i__ <= i__1; ++i__)
2791  {
2792  w[l] = x[i__] + ualpha * p[i__];
2793  ++l;
2794  /* L20: */
2795  }
2796 
2797  (*sfun)(n, &w[lx], &fu, &w[lg]);
2798  ++numf;
2799  gu = ddot_(n, &w[lg], &c__1, &p[1], &c__1);
2800 
2801  /* THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS */
2802  /* OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY */
2803  /* LOWER THAN F AT THE ORIGIN. */
2804 
2805  if (fu <= fmin && fu <= oldf - ualpha * gtest1)
2806  {
2807  dcopy_(n, &w[lg], &c__1, &g[1], &c__1);
2808  }
2809 
2810  goto L10;
2811 
2812  /* IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND */
2813 
2814 L30:
2815  *nftotl = numf;
2816  *iflag = 1;
2817 
2818  if (itest != 0)
2819  {
2820  goto L50;
2821  }
2822 
2823  /* IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE */
2824 
2825  *iflag = 0;
2826  *f = fmin;
2827  *alpha = xmin;
2828  i__1 = *n;
2829 
2830  for (i__ = 1; i__ <= i__1; ++i__)
2831  {
2832  x[i__] += *alpha * p[i__];
2833  /* L40: */
2834  }
2835 
2836 L50:
2837  return 0;
2838 } /* linder_ */
2839 
2840 /* Subroutine */
2842  rtsmll, C_FLOAT64 *reltol, C_FLOAT64 *fabstol, C_FLOAT64 *tnytol,
2843  C_FLOAT64 *fpresn, C_FLOAT64 *eta, C_FLOAT64 *rmu, C_FLOAT64 *
2844  xbnd, C_FLOAT64 *u, C_FLOAT64 *fu, C_FLOAT64 *gu, C_FLOAT64 *xmin,
2845  C_FLOAT64 *fmin, C_FLOAT64 *gmin, C_FLOAT64 *xw, C_FLOAT64 *fw,
2846  C_FLOAT64 *gw, C_FLOAT64 *a, C_FLOAT64 *b, C_FLOAT64 *oldf,
2847  C_FLOAT64 *b1, C_FLOAT64 *scxbnd, C_FLOAT64 *e, C_FLOAT64 *step,
2848  C_FLOAT64 *factor, C_INT *braktd, C_FLOAT64 *gtest1, C_FLOAT64 *
2849  gtest2, C_FLOAT64 *tol, C_INT *ientry, C_INT *itest)
2850 {
2851  /* System generated locals */
2852  C_FLOAT64 d__1, d__2;
2853 
2854  /* Local variables */
2855  C_FLOAT64 half, abgw, fabsr, five, zero, p, q, r__, s, scale,
2856  denom, three, a1, d1, d2, sumsq, point1, abgmin, chordm, eleven,
2857  chordu;
2858  C_INT convrg;
2859  C_FLOAT64 xmidpt, twotol, one;
2860 
2861  /* ************************************************************ */
2862  /* GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY */
2863  /* ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC */
2864  /* INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL */
2865  /* IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A
2866  */
2867  /* POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM. */
2868  /* THE VALUE OF THE C_INT PARAMETERS IENTRY DETERMINES THE PATH TAKEN */
2869  /* THROUGH THE CODE. */
2870  /* ************************************************************ */
2871 
2872  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE CALLED */
2873  /* WITHIN GETPTC */
2874 
2875  zero = 0.;
2876  point1 = .1;
2877  half = .5;
2878  one = 1.;
2879  three = 3.;
2880  five = 5.;
2881  eleven = 11.;
2882 
2883  /* BRANCH TO APPROPRIATE SECTION OF CODE DEPENDING ON THE */
2884  /* VALUE OF IENTRY. */
2885 
2886  switch (*ientry)
2887  {
2888  case 1: goto L10;
2889  case 2: goto L20;
2890  }
2891 
2892  /* IENTRY=1 */
2893  /* CHECK INPUT PARAMETERS */
2894 
2895 L10:
2896  *itest = 2;
2897 
2898  if (*u <= zero || *xbnd <= *tnytol || *gu > zero)
2899  {
2900  return 0;
2901  }
2902 
2903  *itest = 1;
2904 
2905  if (*xbnd < *fabstol)
2906  {
2907  *fabstol = *xbnd;
2908  }
2909 
2910  *tol = *fabstol;
2911  twotol = *tol + *tol;
2912 
2913  /* A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS */
2914  /* WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED. */
2915  /* INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF */
2916  /* FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH */
2917  /* AT VALUES FOR LATEST ESTIMATE AT MINIMUM. */
2918 
2919  *a = zero;
2920  *xw = zero;
2921  *xmin = zero;
2922  *oldf = *fu;
2923  *fmin = *fu;
2924  *fw = *fu;
2925  *gw = *gu;
2926  *gmin = *gu;
2927  *step = *u;
2928  *factor = five;
2929 
2930  /* THE MINIMUM HAS NOT YET BEEN BRACKETED. */
2931 
2932  *braktd = FALSE_;
2933 
2934  /* SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED
2935  */
2936  /* EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.) SET THE UPPER BOUND */
2937  /* ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND). */
2938 
2939  *scxbnd = *xbnd;
2940  *b = *scxbnd + *reltol * fabs(*scxbnd) + *fabstol;
2941  *e = *b + *b;
2942  *b1 = *b;
2943 
2944  /* COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA. */
2945 
2946  *gtest1 = -(*rmu) * *gu;
2947  *gtest2 = -(*eta) * *gu;
2948 
2949  /* SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION */
2950 
2951  *ientry = 2;
2952  goto L210;
2953 
2954  /* IENTRY = 2 */
2955 
2956  /* UPDATE A,B,XW, AND XMIN */
2957 
2958 L20:
2959 
2960  if (*fu > *fmin)
2961  {
2962  goto L60;
2963  }
2964 
2965  /* IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT */
2966  /* ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY. */
2967 
2968  chordu = *oldf - (*xmin + *u) * *gtest1;
2969 
2970  if (*fu <= chordu)
2971  {
2972  goto L30;
2973  }
2974 
2975  /* THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE */
2976  /* CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND */
2977  /* FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF */
2978  /* UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES */
2979  /* THE ROOT OF F(ALPHA)=CHORD(ALPHA). */
2980 
2981  chordm = *oldf - *xmin * *gtest1;
2982  *gu = -(*gmin);
2983  denom = chordm - *fmin;
2984 
2985  if (fabs(denom) >= 1e-15)
2986  {
2987  goto L25;
2988  }
2989 
2990  denom = 1e-15;
2991 
2992  if (chordm - *fmin < 0.)
2993  {
2994  denom = -denom;
2995  }
2996 
2997 L25:
2998 
2999  if (*xmin != zero)
3000  {
3001  *gu = *gmin * (chordu - *fu) / denom;
3002  }
3003 
3004  *fu = half * *u * (*gmin + *gu) + *fmin;
3005 
3006  if (*fu < *fmin)
3007  {
3008  *fu = *fmin;
3009  }
3010 
3011  goto L60;
3012 L30:
3013  *fw = *fmin;
3014  *fmin = *fu;
3015  *gw = *gmin;
3016  *gmin = *gu;
3017  *xmin += *u;
3018  *a -= *u;
3019  *b -= *u;
3020  *xw = -(*u);
3021  *scxbnd -= *u;
3022 
3023  if (*gu <= zero)
3024  {
3025  goto L40;
3026  }
3027 
3028  *b = zero;
3029  *braktd = TRUE_;
3030  goto L50;
3031 L40:
3032  *a = zero;
3033 L50:
3034  *tol = fabs(*xmin) * *reltol + *fabstol;
3035  goto L90;
3036 
3037  /* IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED */
3038  /* BUT NEW POINT MAY NOW QUALIFY AS W. */
3039 
3040 L60:
3041 
3042  if (*u < zero)
3043  {
3044  goto L70;
3045  }
3046 
3047  *b = *u;
3048  *braktd = TRUE_;
3049  goto L80;
3050 L70:
3051  *a = *u;
3052 L80:
3053  *xw = *u;
3054  *fw = *fu;
3055  *gw = *gu;
3056 L90:
3057  twotol = *tol + *tol;
3058  xmidpt = half * (*a + *b);
3059 
3060  /* CHECK TERMINATION CRITERIA */
3061 
3062  convrg = fabs(xmidpt) <= twotol - half * (*b - *a) ||
3063  (fabs(*gmin) <= *gtest2 &&
3064  *fmin < *oldf &&
3065  ((d__1 = *xmin - *xbnd, fabs(d__1)) > * tol ||
3066  !(*braktd)));
3067 
3068  if (! convrg)
3069  {
3070  goto L100;
3071  }
3072 
3073  *itest = 0;
3074 
3075  if (*xmin != zero)
3076  {
3077  return 0;
3078  }
3079 
3080  /* IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE */
3081  /* CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA- */
3082  /* UNIMODALITY CONSTANT, TOL. IF THE CHANGE IN F(X) IS LARGER THAN */
3083  /* EXPECTED, REDUCE THE VALUE OF TOL. */
3084 
3085  *itest = 3;
3086 
3087  if ((d__1 = *oldf - *fw, fabs(d__1)) <= *fpresn * (one + fabs(*oldf)))
3088  {
3089  return 0;
3090  }
3091 
3092  *tol = point1 * *tol;
3093 
3094  if (*tol < *tnytol)
3095  {
3096  return 0;
3097  }
3098 
3099  *reltol = point1 * *reltol;
3100  *fabstol = point1 * *fabstol;
3101  twotol = point1 * twotol;
3102 
3103  /* CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH */
3104 
3105 L100:
3106  r__ = zero;
3107  q = zero;
3108  s = zero;
3109 
3110  if (fabs(*e) <= *tol)
3111  {
3112  goto L150;
3113  }
3114 
3115  /* FIT CUBIC THROUGH XMIN AND XW */
3116 
3117  r__ = three * (*fmin - *fw) / *xw + *gmin + *gw;
3118  fabsr = fabs(r__);
3119  q = fabsr;
3120 
3121  if (*gw == zero || *gmin == zero)
3122  {
3123  goto L140;
3124  }
3125 
3126  /* COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY */
3127  /* WHICH AVOIDS UNDERFLOW AND OVERFLOW. */
3128 
3129  abgw = fabs(*gw);
3130  abgmin = fabs(*gmin);
3131  s = sqrt(abgmin) * sqrt(abgw);
3132 
3133  if (*gw / abgw * *gmin > zero)
3134  {
3135  goto L130;
3136  }
3137 
3138  /* COMPUTE THE SQUARE ROOT OF R*R + S*S. */
3139 
3140  sumsq = one;
3141  p = zero;
3142 
3143  if (fabsr >= s)
3144  {
3145  goto L110;
3146  }
3147 
3148  /* THERE IS A POSSIBILITY OF OVERFLOW. */
3149 
3150  if (s > *rtsmll)
3151  {
3152  p = s * *rtsmll;
3153  }
3154 
3155  if (fabsr >= p)
3156  {
3157  /* Computing 2nd power */
3158  d__1 = fabsr / s;
3159  sumsq = one + d__1 * d__1;
3160  }
3161 
3162  scale = s;
3163  goto L120;
3164 
3165  /* THERE IS A POSSIBILITY OF UNDERFLOW. */
3166 
3167 L110:
3168 
3169  if (fabsr > *rtsmll)
3170  {
3171  p = fabsr * *rtsmll;
3172  }
3173 
3174  if (s >= p)
3175  {
3176  /* Computing 2nd power */
3177  d__1 = s / fabsr;
3178  sumsq = one + d__1 * d__1;
3179  }
3180 
3181  scale = fabsr;
3182 L120:
3183  sumsq = sqrt(sumsq);
3184  q = *big;
3185 
3186  if (scale < *big / sumsq)
3187  {
3188  q = scale * sumsq;
3189  }
3190 
3191  goto L140;
3192 
3193  /* COMPUTE THE SQUARE ROOT OF R*R - S*S */
3194 
3195 L130:
3196  q = sqrt((d__1 = r__ + s, fabs(d__1))) * sqrt((d__2 = r__ - s, fabs(d__2)));
3197 
3198  if (r__ >= s || r__ <= -s)
3199  {
3200  goto L140;
3201  }
3202 
3203  r__ = zero;
3204  q = zero;
3205  goto L150;
3206 
3207  /* COMPUTE THE MINIMUM OF FITTED CUBIC */
3208 
3209 L140:
3210 
3211  if (*xw < zero)
3212  {
3213  q = -q;
3214  }
3215 
3216  s = *xw * (*gmin - r__ - q);
3217  q = *gw - *gmin + q + q;
3218 
3219  if (q > zero)
3220  {
3221  s = -s;
3222  }
3223 
3224  if (q <= zero)
3225  {
3226  q = -q;
3227  }
3228 
3229  r__ = *e;
3230 
3231  if (*b1 != *step || *braktd)
3232  {
3233  *e = *step;
3234  }
3235 
3236  /* CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH */
3237 
3238 L150:
3239  a1 = *a;
3240  *b1 = *b;
3241  *step = xmidpt;
3242 
3243  if (*braktd)
3244  {
3245  goto L160;
3246  }
3247 
3248  *step = -(*factor) * *xw;
3249 
3250  if (*step > *scxbnd)
3251  {
3252  *step = *scxbnd;
3253  }
3254 
3255  if (*step != *scxbnd)
3256  {
3257  *factor = five * *factor;
3258  }
3259 
3260  goto L170;
3261 
3262  /* IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE */
3263  /* WITHIN (A,B). */
3264 
3265 L160:
3266 
3267  if ((*a != zero || *xw >= zero) && (*b != zero || *xw <= zero))
3268  {
3269  goto L180;
3270  }
3271 
3272  /* IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE */
3273  /* WITHIN (A1,B1). */
3274 
3275  d1 = *xw;
3276  d2 = *a;
3277 
3278  if (*a == zero)
3279  {
3280  d2 = *b;
3281  }
3282 
3283  /* THIS LINE MIGHT BE */
3284  /* IF (A .EQ. ZERO) D2 = E */
3285  *u = -d1 / d2;
3286  *step = five * d2 * (point1 + one / *u) / eleven;
3287 
3288  if (*u < one)
3289  {
3290  *step = half * d2 * sqrt(*u);
3291  }
3292 
3293 L170:
3294 
3295  if (*step <= zero)
3296  {
3297  a1 = *step;
3298  }
3299 
3300  if (*step > zero)
3301  {
3302  *b1 = *step;
3303  }
3304 
3305  /* REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE */
3306  /* REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED */
3307  /* DURING THE LAST-BUT-ONE ITERATION. */
3308 
3309 L180:
3310 
3311  if (fabs(s) <= (d__1 = half * q * r__, fabs(d__1)) || s <= q * a1 || s >= q
3312  * *b1)
3313  {
3314  goto L200;
3315  }
3316 
3317  /* A CUBIC INTERPOLATION STEP */
3318 
3319  *step = s / q;
3320 
3321  /* THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B. */
3322 
3323  if (*step - *a >= twotol && *b - *step >= twotol)
3324  {
3325  goto L210;
3326  }
3327 
3328  if (xmidpt > zero)
3329  {
3330  goto L190;
3331  }
3332 
3333  *step = -(*tol);
3334  goto L210;
3335 L190:
3336  *step = *tol;
3337  goto L210;
3338 L200:
3339  *e = *b - *a;
3340 
3341  /* IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO */
3342  /* COMPUTE THE NEW POINT ON THE BOUNDARY). */
3343 
3344 L210:
3345 
3346  if (*step < *scxbnd)
3347  {
3348  goto L220;
3349  }
3350 
3351  *step = *scxbnd;
3352 
3353  /* MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND. */
3354 
3355  *scxbnd -= (*reltol * fabs(*xbnd) + *fabstol) / (one + *reltol);
3356 L220:
3357  *u = *step;
3358 
3359  if (fabs(*step) < *tol && *step < zero)
3360  {
3361  *u = -(*tol);
3362  }
3363 
3364  if (fabs(*step) < *tol && *step >= zero)
3365  {
3366  *u = *tol;
3367  }
3368 
3369  *itest = 1;
3370  return 0;
3371 } /* getptc_ */
3372 
3373 // this special BLAS is copied here without modification by Joseph
3374 
3375 /* DXPY (A VERSION OF DAXPY WITH A=1.0) */
3376 /* WRITTEN BY: STEPHEN G. NASH */
3377 /* OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. */
3378 /* GEORGE MASON UNIVERSITY */
3379 /* FAIRFAX, VA 22030 */
3380 /* ****************************************************************** */
3381 /* SPECIAL BLAS FOR Y = X+Y */
3382 /* ****************************************************************** */
3383 /* Subroutine */ int dxpy_(C_INT *n, C_FLOAT64 *dx, C_INT *incx,
3384  C_FLOAT64 *dy, C_INT *incy)
3385 {
3386  /* System generated locals */
3387  C_INT i__1;
3388 
3389  /* Local variables */
3390  C_INT i__, m, ix, iy, mp1;
3391 
3392  /* VECTOR PLUS A VECTOR. */
3393  /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */
3394  /* STEPHEN G. NASH 5/30/89. */
3395 
3396  /* Parameter adjustments */
3397  --dy;
3398  --dx;
3399 
3400  /* Function Body */
3401  if (*n <= 0)
3402  {
3403  return 0;
3404  }
3405 
3406  if (*incx == 1 && *incy == 1)
3407  {
3408  goto L20;
3409  }
3410 
3411  /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */
3412  /* NOT EQUAL TO 1 */
3413 
3414  ix = 1;
3415  iy = 1;
3416 
3417  if (*incx < 0)
3418  {
3419  ix = (-(*n) + 1) * *incx + 1;
3420  }
3421 
3422  if (*incy < 0)
3423  {
3424  iy = (-(*n) + 1) * *incy + 1;
3425  }
3426 
3427  i__1 = *n;
3428 
3429  for (i__ = 1; i__ <= i__1; ++i__)
3430  {
3431  dy[iy] += dx[ix];
3432  ix += *incx;
3433  iy += *incy;
3434  /* L10: */
3435  }
3436 
3437  return 0;
3438 
3439  /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */
3440 
3441  /* CLEAN-UP LOOP */
3442 
3443 L20:
3444  m = *n % 4;
3445 
3446  if (m == 0)
3447  {
3448  goto L40;
3449  }
3450 
3451  i__1 = m;
3452 
3453  for (i__ = 1; i__ <= i__1; ++i__)
3454  {
3455  dy[i__] += dx[i__];
3456  /* L30: */
3457  }
3458 
3459  if (*n < 4)
3460  {
3461  return 0;
3462  }
3463 
3464 L40:
3465  mp1 = m + 1;
3466  i__1 = *n;
3467 
3468  for (i__ = mp1; i__ <= i__1; i__ += 4)
3469  {
3470  dy[i__] += dx[i__];
3471  dy[i__ + 1] += dx[i__ + 1];
3472  dy[i__ + 2] += dx[i__ + 2];
3473  dy[i__ + 3] += dx[i__ + 3];
3474  /* L50: */
3475  }
3476 
3477  return 0;
3478 } /* dxpy_ */
3479 
3481 {
3482  return(pow(*ap, *bp));
3483 }
#define C_INT
Definition: copasi.h:115
C_FLOAT64 step1_(C_FLOAT64 *fnew, C_FLOAT64 *fm, C_FLOAT64 *gtp, C_FLOAT64 *smax)
int lsout_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int linder_(C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_INT *)
static C_INT c_true
int daxpy_(integer *n, doublereal *da, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int modz_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int negvec_(C_INT *, C_FLOAT64 *)
int ndia3_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
static C_FLOAT64 c_b246
int getptc_(C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
static C_INT c__1
int initpc_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
#define TRUE_
int crash_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
int ssbfgs_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int stpmax_(C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *)
doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int dxpy_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *)
#define subscr_1
int initp3_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
int mslv_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
int modlnp_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
static C_INT c_false
#define subscr_2
C_FLOAT64 mchpr1_(void)
int monit_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *)
int tn_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *)
int setucr_(C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *)
int lmqnbc_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int chkucp_(C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int gtims_(C_FLOAT64 *v, C_FLOAT64 *gv, C_INT *n, C_FLOAT64 *x, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *, FTruncatedNewton *sfun, C_INT *first, C_FLOAT64 *delta, C_FLOAT64 *accrcy, C_FLOAT64 *xnorm)
#define C_FLOAT64
Definition: copasi.h:92
#define subscr_3
int ztime_(C_INT *, C_FLOAT64 *, C_INT *)
int lmqn_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
#define FALSE_
C_FLOAT64 pow_dd(C_FLOAT64 *ap, C_FLOAT64 *bp)
int msolve_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
int cnvtst_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *)
int tnbc_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)