COPASI API  4.16.103
CPraxis.cpp
Go to the documentation of this file.
1 // Begin CVS Header
2 // $Source: /Volumes/Home/Users/shoops/cvs/copasi_dev/copasi/optimization/CPraxis.cpp,v $
3 // $Revision: 1.15 $
4 // $Name: $
5 // $Author: shoops $
6 // $Date: 2012/04/23 21:11:20 $
7 // End CVS Header
8 
9 // Copyright (C) 2012 - 2010 by Pedro Mendes, Virginia Tech Intellectual
10 // Properties, Inc., University of Heidelberg, and The University
11 // of Manchester.
12 // All rights reserved.
13 
14 // Copyright (C) 2008 by Pedro Mendes, Virginia Tech Intellectual
15 // Properties, Inc., EML Research, gGmbH, University of Heidelberg,
16 // and The University of Manchester.
17 // All rights reserved.
18 
19 // Copyright (C) 2001 - 2007 by Pedro Mendes, Virginia Tech Intellectual
20 // Properties, Inc. and EML Research, gGmbH.
21 // All rights reserved.
22 
23 /* praxis.f -- translated by f2c (version 19970805).
24  You must link the resulting object file with the libraries:
25  -lf2c -lm (in that order)
26  */
27 
28 #include <cmath>
29 #include <stdio.h>
30 
31 #include "copasi.h"
32 
33 #include "CPraxis.h"
34 
36 
37 int maprnt_(C_INT *, C_FLOAT64 *, C_INT *, C_INT *);
38 int vcprnt_(C_INT *, C_FLOAT64 *, C_INT *);
39 int sort_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *);
40 
41 /* Common Block Declarations */
42 
43 #define TRUE_ (1)
44 #define FALSE_ (0)
45 #define dmax(a1,a2) ((a1 > a2) ? a1 : a2)
46 #define dmin(a1,a2) ((a1 < a2) ? a1 : a2)
47 
48 /* Table of constant values */
49 
50 static C_INT c__1 = 1;
51 static C_INT c__2 = 2;
52 static bool c_false = FALSE_;
53 static C_INT c__4 = 4;
54 static bool c_true = TRUE_;
55 static C_INT c__3 = 3;
56 static C_INT c__0 = 0;
57 
59  mpRandom(NULL)
60 {
62 }
63 
65 {
67 }
68 
70  C_INT *n, C_INT *prin, C_FLOAT64 *x, FPraxis *f, C_FLOAT64 *fmin)
71 {
72  /* System generated locals */
73  C_INT i__1, i__2, i__3;
74  C_FLOAT64 ret_val, d__1;
75 
76  /* Local variables */
77  static C_FLOAT64 scbd;
78  static C_INT idim;
79  static bool illc;
80  static C_INT klmk;
81  static C_FLOAT64 d__[100], h__, ldfac;
82  static C_INT i__, j, k;
83  static C_FLOAT64 s, t, y[100], large, z__[100], small, value, f1;
84  static C_INT k2;
85  static C_FLOAT64 m2, m4, t2, df, dn;
86  static C_INT kl, ii;
87  static C_FLOAT64 sf;
88  static C_INT kt;
89  static C_FLOAT64 sl, vlarge;
90  static C_FLOAT64 vsmall;
91  static C_INT km1, im1;
92  static C_FLOAT64 dni, lds;
93  static C_INT ktm;
94 
95  C_FLOAT64 lastValue = std::numeric_limits<C_FLOAT64>::infinity();
96 
97  /* LAST MODIFIED 3/1/73 */
98 
99  /* PRAXIS RETURNS THE MINIMUM OF THE FUNCTION F(X,N) OF N VARIABLES */
100  /* USING THE PRINCIPAL AXIS METHOD. THE GRADIENT OF THE FUNCTION IS */
101  /* NOT REQUIRED. */
102 
103  /* FOR A DESCRIPTION OF THE ALGORITHM, SEE CHAPTER SEVEN OF */
104  /* "ALGORITHMS FOR FINDING ZEROS AND EXTREMA OF FUNCTIONS WITHOUT */
105  /* CALCULATING DERIVATIVES" BY RICHARD P BRENT. */
106 
107  /* THE PARAMETERS ARE: */
108  /* T0 IS A TOLERANCE. PRAXIS ATTEMPTS TO RETURN PRAXIS=F(X) */
109  /* SUCH THAT IF X0 IS THE TRUE LOCAL MINIMUM NEAR X, THEN */
110  /* NORM(X-X0) < T0 + SQUAREROOT(MACHEP)*NORM(X). */
111  /* MACHEP IS THE MACHINE PRECISION, THE SMALLEST NUMBER SUCH THAT */
112  /* 1 + MACHEP > 1. MACHEP SHOULD BE 16.**-13 (ABOUT */
113  /* 2.22D-16) FOR REAL*8 ARITHMETIC ON THE IBM 360. */
114  /* H0 IS THE MAXIMUM STEP SIZE. H0 SHOULD BE SET TO ABOUT THE */
115  /* MAXIMUM DISTANCE FROM THE INITIAL GUESS TO THE MINIMUM. */
116  /* (IF H0 IS SET TOO LARGE OR TOO SMALL, THE INITIAL RATE OF */
117  /* CONVERGENCE MAY BE SLOW.) */
118  /* N (AT LEAST TWO) IS THE NUMBER OF VARIABLES UPON WHICH */
119  /* THE FUNCTION DEPENDS. */
120  /* PRIN CONTROLS THE PRINTING OF INTERMEDIATE RESULTS. */
121  /* IF PRIN=0, NOTHING IS PRINTED. */
122  /* IF PRIN=1, F IS PRINTED AFTER EVERY N+1 OR N+2 LINEAR */
123  /* MINIMIZATIONS. FINAL X IS PRINTED, BUT INTERMEDIATE X IS */
124  /* PRINTED ONLY IF N IS AT MOST 4. */
125  /* IF PRIN=2, THE SCALE FACTORS AND THE PRINCIPAL VALUES OF */
126  /* THE APPROXIMATING QUADRATIC FORM ARE ALSO PRINTED. */
127  /* IF PRIN=3, X IS ALSO PRINTED AFTER EVERY FEW LINEAR */
128  /* MINIMIZATIONS. */
129  /* IF PRIN=4, THE PRINCIPAL VECTORS OF THE APPROXIMATING */
130  /* QUADRATIC FORM ARE ALSO PRINTED. */
131  /* X IS AN ARRAY CONTAINING ON ENTRY A GUESS OF THE POINT OF */
132  /* MINIMUM, ON RETURN THE ESTIMATED POINT OF MINIMUM. */
133  /* F(X,N) IS THE FUNCTION TO BE MINIMIZED. F SHOULD BE A REAL*8 */
134  /* FUNCTION DECLARED EXTERNAL IN THE CALLING PROGRAM. */
135  /* FMIN IS AN ESTIMATE OF THE MINIMUM, USED ONLY IN PRINTING */
136  /* INTERMEDIATE RESULTS. */
137  /* THE APPROXIMATING QUADRATIC FORM IS */
138  /* Q(X') = F(X,N) + (1/2) * (X'-X)-TRANSPOSE * A * (X'-X) */
139  /* WHERE X IS THE BEST ESTIMATE OF THE MINIMUM AND A IS */
140  /* INVERSE(V-TRANSPOSE) * D * INVERSE(V) */
141  /* (V(*,*) IS THE MATRIX OF SEARCH DIRECTIONS; D(*) IS THE ARRAY */
142  /* OF SECOND DIFFERENCES). IF F HAS CONTINUOUS SECOND DERIVATIVES */
143  /* NEAR X0, A WILL TEND TO THE HESSIAN OF F AT X0 AS X APPROACHES X0. */
144 
145  /* IT IS ASSUMED THAT ON FLOATING-POINT UNDERFLOW THE RESULT IS SET */
146  /* TO ZERO. */
147  /* THE USER SHOULD OBSERVE THE COMMENT ON HEURISTIC NUMBERS AFTER */
148  /* THE INITIALIZATION OF MACHINE DEPENDENT NUMBERS. */
149 
150  /* .....IF N>20 OR IF N<20 AND YOU NEED MORE SPACE, CHANGE '20' TO THE */
151  /* LARGEST VALUE OF N IN THE NEXT CARD, IN THE CARD 'IDIM=20', AND */
152  /* IN THE DIMENSION STATEMENTS IN SUBROUTINES MINFIT,MIN,FLIN,QUAD. */
153 
154  /* .....INITIALIZATION..... */
155  /* MACHINE DEPENDENT NUMBERS: */
156 
157  /* Parameter adjustments */
158  --x;
159 
160  /* Function Body */
161  small = *machep * *machep;
162  vsmall = small * small;
163  large = 1. / small;
164  vlarge = 1. / vsmall;
165  m2 = sqrt(*machep);
166  m4 = sqrt(m2);
167 
168  /* HEURISTIC NUMBERS: */
169  /* IF THE AXES MAY BE BADLY SCALED (WHICH IS TO BE AVOIDED IF */
170  /* POSSIBLE), THEN SET SCBD=10. OTHERWISE SET SCBD=1. */
171  /* IF THE PROBLEM IS KNOWN TO BE ILL-CONDITIONED, SET ILLC=TRUE. */
172  /* OTHERWISE SET ILLC=FALSE. */
173  /* KTM IS THE NUMBER OF ITERATIONS WITHOUT IMPROVEMENT BEFORE THE */
174  /* ALGORITHM TERMINATES. KTM=4 IS VERY CAUTIOUS; USUALLY KTM=1 */
175  /* IS SATISFACTORY. */
176 
177  scbd = 1.;
178  illc = FALSE_;
179  ktm = 1;
180 
181  ldfac = .01;
182 
183  if (illc)
184  {
185  ldfac = .1;
186  }
187 
188  kt = 0;
189  global_1.nl = 0;
190  global_1.nf = 1;
191  global_1.fx = (*f)(&x[1], n);
192  q_1.qf1 = global_1.fx;
193  t = small + fabs(*t0);
194  t2 = t;
195  global_1.dmin__ = small;
196  h__ = *h0;
197 
198  if (h__ < t * 100)
199  {
200  h__ = t * 100;
201  }
202 
203  global_1.ldt = h__;
204  /* .....THE FIRST SET OF SEARCH DIRECTIONS V IS THE IDENTITY MATRIX.....
205  */
206  i__1 = *n;
207 
208  for (i__ = 1; i__ <= i__1; ++i__)
209  {
210  i__2 = *n;
211 
212  for (j = 1; j <= i__2; ++j)
213  {
214  /* L10: */
215  q_1.v[i__ + j * 100 - 101] = 0.;
216  }
217 
218  /* L20: */
219  q_1.v[i__ + i__ * 100 - 101] = 1.;
220  }
221 
222  d__[0] = 0.;
223  q_1.qd0 = 0.;
224  i__1 = *n;
225 
226  for (i__ = 1; i__ <= i__1; ++i__)
227  {
228  q_1.q0[i__ - 1] = x[i__];
229  /* L30: */
230  q_1.q1[i__ - 1] = x[i__];
231  }
232 
233  if (*prin > 0)
234  {
235  print_(n, &x[1], prin, fmin);
236  }
237 
238  /* .....THE MAIN LOOP STARTS HERE..... */
239 L40:
240  sf = d__[0];
241  d__[0] = 0.;
242  s = 0.;
243 
244  /* .....MINIMIZE ALONG THE FIRST DIRECTION V(*,1). */
245  /* FX MUST BE PASSED TO MIN BY VALUE. */
246  value = global_1.fx;
247  min_(n, &c__1, &c__2, d__, &s, &value, &c_false, f, &x[1], &t,
248  machep, &h__);
249 
250  if (s > 0.)
251  {
252  goto L50;
253  }
254 
255  i__1 = *n;
256 
257  for (i__ = 1; i__ <= i__1; ++i__)
258  {
259  /* L45: */
260  q_1.v[i__ - 1] = -q_1.v[i__ - 1];
261  }
262 
263 L50:
264 
265  if (sf > d__[0] * .9 && sf * .9 < d__[0])
266  {
267  goto L70;
268  }
269 
270  i__1 = *n;
271 
272  for (i__ = 2; i__ <= i__1; ++i__)
273  {
274  /* L60: */
275  d__[i__ - 1] = 0.;
276  }
277 
278  /* .....THE INNER LOOP STARTS HERE..... */
279 L70:
280  i__1 = *n;
281 
282  for (k = 2; k <= i__1; ++k)
283  {
284  i__2 = *n;
285 
286  for (i__ = 1; i__ <= i__2; ++i__)
287  {
288  /* L75: */
289  y[i__ - 1] = x[i__];
290  }
291 
292  sf = global_1.fx;
293 
294  if (kt > 0)
295  {
296  illc = TRUE_;
297  }
298 
299 L80:
300  kl = k;
301  df = 0.;
302 
303  /* .....A RANDOM STEP FOLLOWS (TO AVOID RESOLUTION VALLEYS). */
304  /* PRAXIS ASSUMES THAT RANDOM RETURNS A RANDOM NUMBER UNIFORMLY */
305  /* DISTRIBUTED IN (0,1). */
306 
307  if (! illc)
308  {
309  goto L95;
310  }
311 
312  i__2 = *n;
313 
314  for (i__ = 1; i__ <= i__2; ++i__)
315  {
316  s = (global_1.ldt * .1 + t2 * pow(10.0, (C_FLOAT64)kt)) * (mpRandom->getRandomCC() - 0.5);
317  z__[i__ - 1] = s;
318  i__3 = *n;
319 
320  for (j = 1; j <= i__3; ++j)
321  {
322  /* L85: */
323  x[j] += s * q_1.v[j + i__ * 100 - 101];
324  }
325 
326  /* L90: */
327  }
328 
329  global_1.fx = (*f)(&x[1], n);
330  ++global_1.nf;
331 
332  /* .....MINIMIZE ALONG THE "NON-CONJUGATE" DIRECTIONS V(*,K),...,V(*,N
333  ) */
334 
335 L95:
336  i__2 = *n;
337 
338  for (k2 = k; k2 <= i__2; ++k2)
339  {
340  sl = global_1.fx;
341  s = 0.;
342  value = global_1.fx;
343  min_(n, &k2, &c__2, &d__[k2 - 1], &s, &value, &c_false, f, &
344  x[1], &t, machep, &h__);
345 
346  if (illc)
347  {
348  goto L97;
349  }
350 
351  s = sl - global_1.fx;
352  goto L99;
353 L97:
354  /* Computing 2nd power */
355  d__1 = s + z__[k2 - 1];
356  s = d__[k2 - 1] * (d__1 * d__1);
357 L99:
358 
359  if (df > s)
360  {
361  goto L105;
362  }
363 
364  df = s;
365  kl = k2;
366 L105:
367  ;
368  }
369 
370  if (illc || df >= (d__1 = *machep * 100 * global_1.fx, fabs(d__1)))
371  {
372  goto L110;
373  }
374 
375  /* .....IF THERE WAS NOT MUCH IMPROVEMENT ON THE FIRST TRY, SET */
376  /* ILLC=TRUE AND START THE INNER LOOP AGAIN..... */
377 
378  illc = TRUE_;
379  goto L80;
380 L110:
381 
382  if (k == 2 && *prin > 1)
383  {
384  vcprnt_(&c__1, d__, n);
385  }
386 
387  /* .....MINIMIZE ALONG THE "CONJUGATE" DIRECTIONS V(*,1),...,V(*,K-1)
388  */
389 
390  km1 = k - 1;
391  i__2 = km1;
392 
393  for (k2 = 1; k2 <= i__2; ++k2)
394  {
395  s = 0.;
396  value = global_1.fx;
397  min_(n, &k2, &c__2, &d__[k2 - 1], &s, &value, &c_false, f, &
398  x[1], &t, machep, &h__);
399  /* L120: */
400  }
401 
402  f1 = global_1.fx;
403  global_1.fx = sf;
404  lds = 0.;
405  i__2 = *n;
406 
407  for (i__ = 1; i__ <= i__2; ++i__)
408  {
409  sl = x[i__];
410  x[i__] = y[i__ - 1];
411  sl -= y[i__ - 1];
412  y[i__ - 1] = sl;
413  /* L130: */
414  lds += sl * sl;
415  }
416 
417  lds = sqrt(lds);
418 
419  if (lds <= small)
420  {
421  goto L160;
422  }
423 
424  /* .....DISCARD DIRECTION V(*,KL). */
425  /* IF NO RANDOM STEP WAS TAKEN, V(*,KL) IS THE "NON-CONJUGATE" */
426  /* DIRECTION ALONG WHICH THE GREATEST IMPROVEMENT WAS MADE..... */
427 
428  klmk = kl - k;
429 
430  if (klmk < 1)
431  {
432  goto L141;
433  }
434 
435  i__2 = klmk;
436 
437  for (ii = 1; ii <= i__2; ++ii)
438  {
439  i__ = kl - ii;
440  i__3 = *n;
441 
442  for (j = 1; j <= i__3; ++j)
443  {
444  /* L135: */
445  q_1.v[j + (i__ + 1) * 100 - 101] = q_1.v[j + i__ * 100 - 101];
446  }
447 
448  /* L140: */
449  d__[i__] = d__[i__ - 1];
450  }
451 
452 L141:
453  d__[k - 1] = 0.;
454  i__2 = *n;
455 
456  for (i__ = 1; i__ <= i__2; ++i__)
457  {
458  /* L145: */
459  q_1.v[i__ + k * 100 - 101] = y[i__ - 1] / lds;
460  }
461 
462  /* .....MINIMIZE ALONG THE NEW "CONJUGATE" DIRECTION V(*,K), WHICH IS */
463  /* THE NORMALIZED VECTOR: (NEW X) - (0LD X)..... */
464 
465  value = f1;
466  min_(n, &k, &c__4, &d__[k - 1], &lds, &value, &c_true, f, &x[1],
467  &t, machep, &h__);
468 
469  if (lds > 0.)
470  {
471  goto L160;
472  }
473 
474  lds = -lds;
475  i__2 = *n;
476 
477  for (i__ = 1; i__ <= i__2; ++i__)
478  {
479  /* L150: */
480  q_1.v[i__ + k * 100 - 101] = -q_1.v[i__ + k * 100 - 101];
481  }
482 
483 L160:
484  global_1.ldt = ldfac * global_1.ldt;
485 
486  if (global_1.ldt < lds)
487  {
488  global_1.ldt = lds;
489  }
490 
491  if (*prin > 0)
492  {
493  print_(n, &x[1], prin, fmin);
494  }
495 
496  t2 = 0.;
497  i__2 = *n;
498 
499  for (i__ = 1; i__ <= i__2; ++i__)
500  {
501  /* L165: */
502  /* Computing 2nd power */
503  d__1 = x[i__];
504  t2 += d__1 * d__1;
505  }
506 
507  t2 = m2 * sqrt(t2) + t;
508 
509  /* .....SEE WHETHER THE LENGTH OF THE STEP TAKEN SINCE STARTING THE */
510  /* INNER LOOP EXCEEDS HALF THE TOLERANCE..... */
511 
512  if (global_1.ldt > t2 * .5f)
513  {
514  kt = -1;
515  }
516 
517  ++kt;
518 
519  if (kt > ktm)
520  {
521  goto L400;
522  }
523 
524  /* L170: */
525  }
526 
527  /* added manually by Pedro Mendes 11/11/1998 */
528  // if(callback != 0/*NULL*/) callback(global_1.fx);
529  /* .....THE INNER LOOP ENDS HERE. */
530 
531  /* TRY QUADRATIC EXTRAPOLATION IN CASE WE ARE IN A CURVED VALLEY. */
532 
533  /* L171: */
534  quad_(n, f, &x[1], &t, machep, &h__);
535  dn = 0.;
536  i__1 = *n;
537 
538  for (i__ = 1; i__ <= i__1; ++i__)
539  {
540  d__[i__ - 1] = 1. / sqrt(d__[i__ - 1]);
541 
542  if (dn < d__[i__ - 1])
543  {
544  dn = d__[i__ - 1];
545  }
546 
547  /* L175: */
548  }
549 
550  if (*prin > 3)
551  {
552  maprnt_(&c__1, q_1.v, &idim, n);
553  }
554 
555  i__1 = *n;
556 
557  for (j = 1; j <= i__1; ++j)
558  {
559  s = d__[j - 1] / dn;
560  i__2 = *n;
561 
562  for (i__ = 1; i__ <= i__2; ++i__)
563  {
564  /* L180: */
565  q_1.v[i__ + j * 100 - 101] = s * q_1.v[i__ + j * 100 - 101];
566  }
567  }
568 
569  /* .....SCALE THE AXES TO TRY TO REDUCE THE CONDITION NUMBER..... */
570 
571  if (scbd <= 1.)
572  {
573  goto L200;
574  }
575 
576  s = vlarge;
577  i__2 = *n;
578 
579  for (i__ = 1; i__ <= i__2; ++i__)
580  {
581  sl = 0.;
582  i__1 = *n;
583 
584  for (j = 1; j <= i__1; ++j)
585  {
586  /* L182: */
587  sl += q_1.v[i__ + j * 100 - 101] * q_1.v[i__ + j * 100 - 101];
588  }
589 
590  z__[i__ - 1] = sqrt(sl);
591 
592  if (z__[i__ - 1] < m4)
593  {
594  z__[i__ - 1] = m4;
595  }
596 
597  if (s > z__[i__ - 1])
598  {
599  s = z__[i__ - 1];
600  }
601 
602  /* L185: */
603  }
604 
605  i__2 = *n;
606 
607  for (i__ = 1; i__ <= i__2; ++i__)
608  {
609  sl = s / z__[i__ - 1];
610  z__[i__ - 1] = 1. / sl;
611 
612  if (z__[i__ - 1] <= scbd)
613  {
614  goto L189;
615  }
616 
617  sl = 1. / scbd;
618  z__[i__ - 1] = scbd;
619 L189:
620  i__1 = *n;
621 
622  for (j = 1; j <= i__1; ++j)
623  {
624  /* L190: */
625  q_1.v[i__ + j * 100 - 101] = sl * q_1.v[i__ + j * 100 - 101];
626  }
627 
628  /* L195: */
629  }
630 
631  /* .....CALCULATE A NEW SET OF ORTHOGONAL DIRECTIONS BEFORE REPEATING */
632  /* THE MAIN LOOP. */
633  /* FIRST TRANSPOSE V FOR MINFIT: */
634 
635 L200:
636  i__2 = *n;
637 
638  for (i__ = 2; i__ <= i__2; ++i__)
639  {
640  im1 = i__ - 1;
641  i__1 = im1;
642 
643  for (j = 1; j <= i__1; ++j)
644  {
645  s = q_1.v[i__ + j * 100 - 101];
646  q_1.v[i__ + j * 100 - 101] = q_1.v[j + i__ * 100 - 101];
647  /* L210: */
648  q_1.v[j + i__ * 100 - 101] = s;
649  }
650 
651  /* L220: */
652  }
653 
654  /* .....CALL MINFIT TO FIND THE SINGULAR VALUE DECOMPOSITION OF V. */
655  /* THIS GIVES THE PRINCIPAL VALUES AND PRINCIPAL DIRECTIONS OF THE */
656  /* APPROXIMATING QUADRATIC FORM WITHOUT SQUARING THE CONDITION */
657  /* NUMBER..... */
658 
659  minfit_(&idim, n, machep, &vsmall, q_1.v, d__);
660 
661  /* .....UNSCALE THE AXES..... */
662 
663  if (scbd <= 1.)
664  {
665  goto L250;
666  }
667 
668  i__2 = *n;
669 
670  for (i__ = 1; i__ <= i__2; ++i__)
671  {
672  s = z__[i__ - 1];
673  i__1 = *n;
674 
675  for (j = 1; j <= i__1; ++j)
676  {
677  /* L225: */
678  q_1.v[i__ + j * 100 - 101] = s * q_1.v[i__ + j * 100 - 101];
679  }
680 
681  /* L230: */
682  }
683 
684  i__2 = *n;
685 
686  for (i__ = 1; i__ <= i__2; ++i__)
687  {
688  s = 0.;
689  i__1 = *n;
690 
691  for (j = 1; j <= i__1; ++j)
692  {
693  /* L235: */
694  /* Computing 2nd power */
695  d__1 = q_1.v[j + i__ * 100 - 101];
696  s += d__1 * d__1;
697  }
698 
699  s = sqrt(s);
700  d__[i__ - 1] = s * d__[i__ - 1];
701  s = 1 / s;
702  i__1 = *n;
703 
704  for (j = 1; j <= i__1; ++j)
705  {
706  /* L240: */
707  q_1.v[j + i__ * 100 - 101] = s * q_1.v[j + i__ * 100 - 101];
708  }
709 
710  /* L245: */
711  }
712 
713 L250:
714  i__2 = *n;
715 
716  for (i__ = 1; i__ <= i__2; ++i__)
717  {
718  dni = dn * d__[i__ - 1];
719 
720  if (dni > large)
721  {
722  goto L265;
723  }
724 
725  if (dni < small)
726  {
727  goto L260;
728  }
729 
730  d__[i__ - 1] = 1 / (dni * dni);
731  goto L270;
732 L260:
733  d__[i__ - 1] = vlarge;
734  goto L270;
735 L265:
736  d__[i__ - 1] = vsmall;
737 L270:
738  ;
739  }
740 
741  /* .....SORT THE EIGENVALUES AND EIGENVECTORS..... */
742 
743  sort_(&idim, n, d__, q_1.v);
744  global_1.dmin__ = d__[*n - 1];
745 
746  if (global_1.dmin__ < small)
747  {
748  global_1.dmin__ = small;
749  }
750 
751  illc = FALSE_;
752 
753  if (m2 * d__[0] > global_1.dmin__)
754  {
755  illc = TRUE_;
756  }
757 
758  if (*prin > 1 && scbd > 1.)
759  {
760  vcprnt_(&c__2, z__, n);
761  }
762 
763  if (*prin > 1)
764  {
765  vcprnt_(&c__3, d__, n);
766  }
767 
768  if (*prin > 3)
769  {
770  maprnt_(&c__2, q_1.v, &idim, n);
771  }
772 
773  /* .....THE MAIN LOOP ENDS HERE..... */
774 
775  /* added manually by Pedro Mendes 29/1/1998 */
776  /* if(callback != 0) callback(global_1.fx);*/
777 
778  // We need a stopping condition for 1 dimensional problems.
779  // We compare the values between the last and current steps.
780  if (*n > 1 ||
781  global_1.fx < lastValue)
782  {
783  lastValue = global_1.fx;
784  goto L40;
785  }
786 
787  /* .....RETURN..... */
788 
789 L400:
790 
791  if (*prin > 0)
792  {
793  vcprnt_(&c__4, &x[1], n);
794  }
795 
796  ret_val = global_1.fx;
797  return ret_val;
798 } /* praxis_ */
799 
800 /* Subroutine */ int CPraxis::minfit_(C_INT *m, C_INT *n, C_FLOAT64 *machep,
801  C_FLOAT64 *tol, C_FLOAT64 *ab, C_FLOAT64 *q)
802 {
803  /* System generated locals */
804  C_INT ab_dim1, ab_offset, i__1, i__2, i__3;
805  C_FLOAT64 d__1, d__2;
806 
807  /* Local variables */
808  static C_FLOAT64 temp, c__, e[100], f, g, h__;
809  static C_INT i__, j, k, l;
810  static C_FLOAT64 s, x, y, z__;
811  static C_INT l2, ii, kk, kt, ll2, lp1;
812  static C_FLOAT64 eps;
813 
814  /* ...AN IMPROVED VERSION OF MINFIT (SEE GOLUB AND REINSCH, 1969) */
815  /* RESTRICTED TO M=N,P=0. */
816  /* THE SINGULAR VALUES OF THE ARRAY AB ARE RETURNED IN Q AND AB IS */
817  /* OVERWRITTEN WITH THE ORTHOGONAL MATRIX V SUCH THAT U.DIAG(Q) = AB.V,
818  */
819  /* WHERE U IS ANOTHER ORTHOGONAL MATRIX. */
820  /* ...HOUSEHOLDER'S REDUCTION TO BIDIAGONAL FORM... */
821  /* Parameter adjustments */
822  --q;
823  ab_dim1 = *m;
824  ab_offset = ab_dim1 + 1;
825  ab -= ab_offset;
826 
827  /* Function Body */
828  if (*n == 1)
829  {
830  goto L200;
831  }
832 
833  eps = *machep;
834  g = 0.;
835  x = 0.;
836  i__1 = *n;
837 
838  for (i__ = 1; i__ <= i__1; ++i__)
839  {
840  e[i__ - 1] = g;
841  s = 0.;
842  l = i__ + 1;
843  i__2 = *n;
844 
845  for (j = i__; j <= i__2; ++j)
846  {
847  /* L1: */
848  /* Computing 2nd power */
849  d__1 = ab[j + i__ * ab_dim1];
850  s += d__1 * d__1;
851  }
852 
853  g = 0.;
854 
855  if (s < *tol)
856  {
857  goto L4;
858  }
859 
860  f = ab[i__ + i__ * ab_dim1];
861  g = sqrt(s);
862 
863  if (f >= 0.)
864  {
865  g = -g;
866  }
867 
868  h__ = f * g - s;
869  ab[i__ + i__ * ab_dim1] = f - g;
870 
871  if (l > *n)
872  {
873  goto L4;
874  }
875 
876  i__2 = *n;
877 
878  for (j = l; j <= i__2; ++j)
879  {
880  f = 0.;
881  i__3 = *n;
882 
883  for (k = i__; k <= i__3; ++k)
884  {
885  /* L2: */
886  f += ab[k + i__ * ab_dim1] * ab[k + j * ab_dim1];
887  }
888 
889  f /= h__;
890  i__3 = *n;
891 
892  for (k = i__; k <= i__3; ++k)
893  {
894  /* L3: */
895  ab[k + j * ab_dim1] += f * ab[k + i__ * ab_dim1];
896  }
897  }
898 
899 L4:
900  q[i__] = g;
901  s = 0.;
902 
903  if (i__ == *n)
904  {
905  goto L6;
906  }
907 
908  i__3 = *n;
909 
910  for (j = l; j <= i__3; ++j)
911  {
912  /* L5: */
913  s += ab[i__ + j * ab_dim1] * ab[i__ + j * ab_dim1];
914  }
915 
916 L6:
917  g = 0.;
918 
919  if (s < *tol)
920  {
921  goto L10;
922  }
923 
924  if (i__ == *n)
925  {
926  goto L16;
927  }
928 
929  f = ab[i__ + (i__ + 1) * ab_dim1];
930 L16:
931  g = sqrt(s);
932 
933  if (f >= 0.)
934  {
935  g = -g;
936  }
937 
938  h__ = f * g - s;
939 
940  if (i__ == *n)
941  {
942  goto L10;
943  }
944 
945  ab[i__ + (i__ + 1) * ab_dim1] = f - g;
946  i__3 = *n;
947 
948  for (j = l; j <= i__3; ++j)
949  {
950  /* L7: */
951  e[j - 1] = ab[i__ + j * ab_dim1] / h__;
952  }
953 
954  i__3 = *n;
955 
956  for (j = l; j <= i__3; ++j)
957  {
958  s = 0.;
959  i__2 = *n;
960 
961  for (k = l; k <= i__2; ++k)
962  {
963  /* L8: */
964  s += ab[j + k * ab_dim1] * ab[i__ + k * ab_dim1];
965  }
966 
967  i__2 = *n;
968 
969  for (k = l; k <= i__2; ++k)
970  {
971  /* L9: */
972  ab[j + k * ab_dim1] += s * e[k - 1];
973  }
974  }
975 
976 L10:
977  y = (d__1 = q[i__], fabs(d__1)) + (d__2 = e[i__ - 1], fabs(d__2));
978 
979  /* L11: */
980  if (y > x)
981  {
982  x = y;
983  }
984  }
985 
986  /* ...ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS... */
987  ab[*n + *n * ab_dim1] = 1.;
988  g = e[*n - 1];
989  l = *n;
990  i__1 = *n;
991 
992  for (ii = 2; ii <= i__1; ++ii)
993  {
994  i__ = *n - ii + 1;
995 
996  if (g == 0.)
997  {
998  goto L23;
999  }
1000 
1001  h__ = ab[i__ + (i__ + 1) * ab_dim1] * g;
1002  i__2 = *n;
1003 
1004  for (j = l; j <= i__2; ++j)
1005  {
1006  /* L20: */
1007  ab[j + i__ * ab_dim1] = ab[i__ + j * ab_dim1] / h__;
1008  }
1009 
1010  i__2 = *n;
1011 
1012  for (j = l; j <= i__2; ++j)
1013  {
1014  s = 0.;
1015  i__3 = *n;
1016 
1017  for (k = l; k <= i__3; ++k)
1018  {
1019  /* L21: */
1020  s += ab[i__ + k * ab_dim1] * ab[k + j * ab_dim1];
1021  }
1022 
1023  i__3 = *n;
1024 
1025  for (k = l; k <= i__3; ++k)
1026  {
1027  /* L22: */
1028  ab[k + j * ab_dim1] += s * ab[k + i__ * ab_dim1];
1029  }
1030  }
1031 
1032 L23:
1033  i__3 = *n;
1034 
1035  for (j = l; j <= i__3; ++j)
1036  {
1037  ab[i__ + j * ab_dim1] = 0.;
1038  /* L24: */
1039  ab[j + i__ * ab_dim1] = 0.;
1040  }
1041 
1042  ab[i__ + i__ * ab_dim1] = 1.;
1043  g = e[i__ - 1];
1044  /* L25: */
1045  l = i__;
1046  }
1047 
1048  /* ...DIAGONALIZATION OF THE BIDIAGONAL FORM... */
1049  /* L100: */
1050  eps *= x;
1051  i__1 = *n;
1052 
1053  for (kk = 1; kk <= i__1; ++kk)
1054  {
1055  k = *n - kk + 1;
1056  kt = 0;
1057 L101:
1058  ++kt;
1059 
1060  if (kt <= 30)
1061  {
1062  goto L102;
1063  }
1064 
1065  e[k - 1] = 0.;
1066  printf("QR FAILED\n");
1067 L102:
1068  i__3 = k;
1069 
1070  for (ll2 = 1; ll2 <= i__3; ++ll2)
1071  {
1072  l2 = k - ll2 + 1;
1073  l = l2;
1074 
1075  if ((d__1 = e[l - 1], fabs(d__1)) <= eps)
1076  {
1077  goto L120;
1078  }
1079 
1080  if (l == 1)
1081  {
1082  goto L103;
1083  }
1084 
1085  if ((d__1 = q[l - 1], fabs(d__1)) <= eps)
1086  {
1087  goto L110;
1088  }
1089 
1090 L103:
1091  ;
1092  }
1093 
1094  /* ...CANCELLATION OF E(L) IF L>1... */
1095 L110:
1096  c__ = 0.;
1097  s = 1.;
1098  i__3 = k;
1099 
1100  for (i__ = l; i__ <= i__3; ++i__)
1101  {
1102  f = s * e[i__ - 1];
1103  e[i__ - 1] = c__ * e[i__ - 1];
1104 
1105  if (fabs(f) <= eps)
1106  {
1107  goto L120;
1108  }
1109 
1110  g = q[i__];
1111 
1112  /* ...Q(I) = H = DSQRT(G*G + F*F)... */
1113  if (fabs(f) < fabs(g))
1114  {
1115  goto L113;
1116  }
1117 
1118  if (f != 0.)
1119  {
1120  goto L112;
1121  }
1122  else
1123  {
1124  goto L111;
1125  }
1126 
1127 L111:
1128  h__ = 0.;
1129  goto L114;
1130 L112:
1131  /* Computing 2nd power */
1132  d__1 = g / f;
1133  h__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1134  goto L114;
1135 L113:
1136  /* Computing 2nd power */
1137  d__1 = f / g;
1138  h__ = fabs(g) * sqrt(d__1 * d__1 + 1);
1139 L114:
1140  q[i__] = h__;
1141 
1142  if (h__ != 0.)
1143  {
1144  goto L115;
1145  }
1146 
1147  g = 1.;
1148  h__ = 1.;
1149 L115:
1150  c__ = g / h__;
1151  /* L116: */
1152  s = -f / h__;
1153  }
1154 
1155  /* ...TEST FOR CONVERGENCE... */
1156 L120:
1157  z__ = q[k];
1158 
1159  if (l == k)
1160  {
1161  goto L140;
1162  }
1163 
1164  /* ...SHIFT FROM BOTTOM 2*2 MINOR... */
1165  x = q[l];
1166  y = q[k - 1];
1167  g = e[k - 2];
1168  h__ = e[k - 1];
1169  f = ((y - z__) * (y + z__) + (g - h__) * (g + h__)) / (h__ * 2 * y);
1170  g = sqrt(f * f + 1.);
1171  temp = f - g;
1172 
1173  if (f >= 0.)
1174  {
1175  temp = f + g;
1176  }
1177 
1178  f = ((x - z__) * (x + z__) + h__ * (y / temp - h__)) / x;
1179  /* ...NEXT QR TRANSFORMATION... */
1180  c__ = 1.;
1181  s = 1.;
1182  lp1 = l + 1;
1183 
1184  if (lp1 > k)
1185  {
1186  goto L133;
1187  }
1188 
1189  i__3 = k;
1190 
1191  for (i__ = lp1; i__ <= i__3; ++i__)
1192  {
1193  g = e[i__ - 1];
1194  y = q[i__];
1195  h__ = s * g;
1196  g *= c__;
1197 
1198  if (fabs(f) < fabs(h__))
1199  {
1200  goto L123;
1201  }
1202 
1203  if (f != 0.)
1204  {
1205  goto L122;
1206  }
1207  else
1208  {
1209  goto L121;
1210  }
1211 
1212 L121:
1213  z__ = 0.;
1214  goto L124;
1215 L122:
1216  /* Computing 2nd power */
1217  d__1 = h__ / f;
1218  z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1219  goto L124;
1220 L123:
1221  /* Computing 2nd power */
1222  d__1 = f / h__;
1223  z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
1224 L124:
1225  e[i__ - 2] = z__;
1226 
1227  if (z__ != 0.)
1228  {
1229  goto L125;
1230  }
1231 
1232  f = 1.;
1233  z__ = 1.;
1234 L125:
1235  c__ = f / z__;
1236  s = h__ / z__;
1237  f = x * c__ + g * s;
1238  g = -x * s + g * c__;
1239  h__ = y * s;
1240  y *= c__;
1241  i__2 = *n;
1242 
1243  for (j = 1; j <= i__2; ++j)
1244  {
1245  x = ab[j + (i__ - 1) * ab_dim1];
1246  z__ = ab[j + i__ * ab_dim1];
1247  ab[j + (i__ - 1) * ab_dim1] = x * c__ + z__ * s;
1248  /* L126: */
1249  ab[j + i__ * ab_dim1] = -x * s + z__ * c__;
1250  }
1251 
1252  if (fabs(f) < fabs(h__))
1253  {
1254  goto L129;
1255  }
1256 
1257  if (f != 0.)
1258  {
1259  goto L128;
1260  }
1261  else
1262  {
1263  goto L127;
1264  }
1265 
1266 L127:
1267  z__ = 0.;
1268  goto L130;
1269 L128:
1270  /* Computing 2nd power */
1271  d__1 = h__ / f;
1272  z__ = fabs(f) * sqrt(d__1 * d__1 + 1);
1273  goto L130;
1274 L129:
1275  /* Computing 2nd power */
1276  d__1 = f / h__;
1277  z__ = fabs(h__) * sqrt(d__1 * d__1 + 1);
1278 L130:
1279  q[i__ - 1] = z__;
1280 
1281  if (z__ != 0.)
1282  {
1283  goto L131;
1284  }
1285 
1286  f = 1.;
1287  z__ = 1.;
1288 L131:
1289  c__ = f / z__;
1290  s = h__ / z__;
1291  f = c__ * g + s * y;
1292  /* L132: */
1293  x = -s * g + c__ * y;
1294  }
1295 
1296 L133:
1297  e[l - 1] = 0.;
1298  e[k - 1] = f;
1299  q[k] = x;
1300  goto L101;
1301  /* ...CONVERGENCE: Q(K) IS MADE NON-NEGATIVE... */
1302 L140:
1303 
1304  if (z__ >= 0.)
1305  {
1306  goto L150;
1307  }
1308 
1309  q[k] = -z__;
1310  i__3 = *n;
1311 
1312  for (j = 1; j <= i__3; ++j)
1313  {
1314  /* L141: */
1315  ab[j + k * ab_dim1] = -ab[j + k * ab_dim1];
1316  }
1317 
1318 L150:
1319  ;
1320  }
1321 
1322  return 0;
1323 L200:
1324  q[1] = ab[ab_dim1 + 1];
1325  ab[ab_dim1 + 1] = 1.;
1326  return 0;
1327 } /* minfit_ */
1328 
1329 /* Subroutine */ int CPraxis::min_(C_INT *n, C_INT *j, C_INT *nits, C_FLOAT64 *
1330  d2, C_FLOAT64 *x1, C_FLOAT64 *f1, bool *fk, FPraxis *f, C_FLOAT64 *
1331  x, C_FLOAT64 *t, C_FLOAT64 *machep, C_FLOAT64 *h__)
1332 {
1333  /* System generated locals */
1334  C_INT i__1;
1335  C_FLOAT64 d__1, d__2;
1336 
1337  /* Local variables */
1338  static C_FLOAT64 temp;
1339  static C_INT i__, k;
1340  static C_FLOAT64 s, small, d1, f0, f2, m2, m4, t2, x2, fm;
1341  static bool dz;
1342  static C_FLOAT64 xm, sf1, sx1;
1343 
1344  /* ...THE SUBROUTINE MIN MINIMIZES F FROM X IN THE DIRECTION V(*,J) UNLESS */
1345  /* J IS LESS THAN 1, WHEN A QUADRATIC SEARCH IS MADE IN THE PLANE */
1346  /* DEFINED BY Q0,Q1,X. */
1347  /* D2 IS EITHER ZERO OR AN APPROXIMATION TO HALF F". */
1348  /* ON ENTRY, X1 IS AN ESTIMATE OF THE DISTANCE FROM X TO THE MINIMUM */
1349  /* ALONG V(*,J) (OR, IF J=0, A CURVE). ON RETURN, X1 IS THE DISTANCE */
1350  /* FOUND. */
1351  /* IF FK=.TRUE., THEN F1 IS FLIN(X1). OTHERWISE X1 AND F1 ARE IGNORED */
1352  /* ON ENTRY UNLESS FINAL FX IS GREATER THAN F1. */
1353  /* NITS CONTROLS THE NUMBER OF TIMES AN ATTEMPT WILL BE MADE TO HALVE */
1354  /* THE INTERVAL. */
1355  /* Parameter adjustments */
1356  --x;
1357 
1358  /* Function Body */
1359  /* Computing 2nd power */
1360  d__1 = *machep;
1361  small = d__1 * d__1;
1362  m2 = sqrt(*machep);
1363  m4 = sqrt(m2);
1364  sf1 = *f1;
1365  sx1 = *x1;
1366  k = 0;
1367  xm = 0.;
1368  fm = global_1.fx;
1369  f0 = global_1.fx;
1370  dz = *d2 < *machep;
1371  /* ...FIND THE STEP SIZE... */
1372  s = 0.;
1373  i__1 = *n;
1374 
1375  for (i__ = 1; i__ <= i__1; ++i__)
1376  {
1377  /* L1: */
1378  /* Computing 2nd power */
1379  d__1 = x[i__];
1380  s += d__1 * d__1;
1381  }
1382 
1383  s = sqrt(s);
1384  temp = *d2;
1385 
1386  if (dz)
1387  {
1388  temp = global_1.dmin__;
1389  }
1390 
1391  t2 = m4 * sqrt(fabs(global_1.fx) / temp + s * global_1.ldt) + m2 *
1392  global_1.ldt;
1393  s = m4 * s + *t;
1394 
1395  if (dz && t2 > s)
1396  {
1397  t2 = s;
1398  }
1399 
1400  t2 = dmax(t2, small);
1401  /* Computing MIN */
1402  d__1 = t2, d__2 = *h__ * .01;
1403  t2 = dmin(d__1, d__2);
1404 
1405  if (!(*fk) || *f1 > fm)
1406  {
1407  goto L2;
1408  }
1409 
1410  xm = *x1;
1411  fm = *f1;
1412 L2:
1413 
1414  if (*fk && fabs(*x1) >= t2)
1415  {
1416  goto L3;
1417  }
1418 
1419  temp = 1.;
1420 
1421  if (*x1 < 0.)
1422  {
1423  temp = -1.;
1424  }
1425 
1426  *x1 = temp * t2;
1427  *f1 = flin_(n, j, x1, f, &x[1], &global_1.nf);
1428 L3:
1429 
1430  if (*f1 > fm)
1431  {
1432  goto L4;
1433  }
1434 
1435  xm = *x1;
1436  fm = *f1;
1437 L4:
1438 
1439  if (! dz)
1440  {
1441  goto L6;
1442  }
1443 
1444  /* ...EVALUATE FLIN AT ANOTHER POINT AND ESTIMATE THE SECOND DERIVATIVE...
1445  */
1446  x2 = -(*x1);
1447 
1448  if (f0 >= *f1)
1449  {
1450  x2 = *x1 * 2.;
1451  }
1452 
1453  f2 = flin_(n, j, &x2, f, &x[1], &global_1.nf);
1454 
1455  if (f2 > fm)
1456  {
1457  goto L5;
1458  }
1459 
1460  xm = x2;
1461  fm = f2;
1462 L5:
1463  *d2 = (x2 * (*f1 - f0) - *x1 * (f2 - f0)) / (*x1 * x2 * (*x1 - x2));
1464  /* ...ESTIMATE THE FIRST DERIVATIVE AT 0... */
1465 L6:
1466  d1 = (*f1 - f0) / *x1 - *x1 * *d2;
1467  dz = TRUE_;
1468 
1469  /* ...PREDICT THE MINIMUM... */
1470  if (*d2 > small)
1471  {
1472  goto L7;
1473  }
1474 
1475  x2 = *h__;
1476 
1477  if (d1 >= 0.)
1478  {
1479  x2 = -x2;
1480  }
1481 
1482  goto L8;
1483 L7:
1484  x2 = d1 * -.5 / *d2;
1485 L8:
1486 
1487  if (fabs(x2) <= *h__)
1488  {
1489  goto L11;
1490  }
1491 
1492  if (x2 <= 0.)
1493  {
1494  goto L9;
1495  }
1496  else
1497  {
1498  goto L10;
1499  }
1500 
1501 L9:
1502  x2 = -(*h__);
1503  goto L11;
1504 L10:
1505  x2 = *h__;
1506  /* ...EVALUATE F AT THE PREDICTED MINIMUM... */
1507 L11:
1508  f2 = flin_(n, j, &x2, f, &x[1], &global_1.nf);
1509 
1510  if (k >= *nits || f2 <= f0)
1511  {
1512  goto L12;
1513  }
1514 
1515  /* ...NO SUCCESS, SO TRY AGAIN... */
1516  ++k;
1517 
1518  if (f0 < *f1 && *x1 * x2 > 0.)
1519  {
1520  goto L4;
1521  }
1522 
1523  x2 *= .5;
1524  goto L11;
1525  /* ...INCREMENT THE ONE-DIMENSIONAL SEARCH COUNTER... */
1526 L12:
1527  ++global_1.nl;
1528 
1529  if (f2 <= fm)
1530  {
1531  goto L13;
1532  }
1533 
1534  x2 = xm;
1535  goto L14;
1536 L13:
1537  fm = f2;
1538  /* ...GET A NEW ESTIMATE OF THE SECOND DERIVATIVE... */
1539 L14:
1540 
1541  if ((d__1 = x2 * (x2 - *x1), fabs(d__1)) <= small)
1542  {
1543  goto L15;
1544  }
1545 
1546  *d2 = (x2 * (*f1 - f0) - *x1 * (fm - f0)) / (*x1 * x2 * (*x1 - x2));
1547  goto L16;
1548 L15:
1549 
1550  if (k > 0)
1551  {
1552  *d2 = 0.;
1553  }
1554 
1555 L16:
1556 
1557  if (*d2 <= small)
1558  {
1559  *d2 = small;
1560  }
1561 
1562  *x1 = x2;
1563  global_1.fx = fm;
1564 
1565  if (sf1 >= global_1.fx)
1566  {
1567  goto L17;
1568  }
1569 
1570  global_1.fx = sf1;
1571  *x1 = sx1;
1572  /* ...UPDATE X FOR LINEAR BUT NOT PARABOLIC SEARCH... */
1573 L17:
1574 
1575  if (*j == 0)
1576  {
1577  return 0;
1578  }
1579 
1580  i__1 = *n;
1581 
1582  for (i__ = 1; i__ <= i__1; ++i__)
1583  {
1584  /* L18: */
1585  x[i__] += *x1 * q_1.v[i__ + *j * 100 - 101];
1586  }
1587 
1588  return 0;
1589 } /* min_ */
1590 
1592  C_INT *nf)
1593 {
1594  /* System generated locals */
1595  C_INT i__1;
1596  C_FLOAT64 ret_val;
1597 
1598  /* Local variables */
1599  static C_INT i__;
1600  static C_FLOAT64 t[100];
1601 
1602  /* ...FLIN IS THE FUNCTION OF ONE REAL VARIABLE L THAT IS MINIMIZED */
1603  /* BY THE SUBROUTINE MIN... */
1604  /* Parameter adjustments */
1605  --x;
1606 
1607  /* Function Body */
1608  if (*j == 0)
1609  {
1610  goto L2;
1611  }
1612 
1613  /* ...THE SEARCH IS LINEAR... */
1614  i__1 = *n;
1615 
1616  for (i__ = 1; i__ <= i__1; ++i__)
1617  {
1618  /* L1: */
1619  t[i__ - 1] = x[i__] + *l * q_1.v[i__ + *j * 100 - 101];
1620  }
1621 
1622  goto L4;
1623  /* ...THE SEARCH IS ALONG A PARABOLIC SPACE CURVE... */
1624 L2:
1625  q_1.qa = *l * (*l - q_1.qd1) / (q_1.qd0 * (q_1.qd0 + q_1.qd1));
1626  q_1.qb = (*l + q_1.qd0) * (q_1.qd1 - *l) / (q_1.qd0 * q_1.qd1);
1627  q_1.qc = *l * (*l + q_1.qd0) / (q_1.qd1 * (q_1.qd0 + q_1.qd1));
1628  i__1 = *n;
1629 
1630  for (i__ = 1; i__ <= i__1; ++i__)
1631  {
1632  /* L3: */
1633  t[i__ -1] = q_1.qa * q_1.q0[i__ - 1] + q_1.qb * x[i__] + q_1.qc *
1634  q_1.q1[i__ - 1];
1635  }
1636 
1637  /* ...THE FUNCTION EVALUATION COUNTER NF IS INCREMENTED... */
1638 L4:
1639  ++(*nf);
1640  ret_val = (*f)(t, n);
1641  return ret_val;
1642 } /* flin_ */
1643 
1644 /* Subroutine */ int sort_(C_INT *m, C_INT *n, C_FLOAT64 *d__,
1645  C_FLOAT64 *v)
1646 {
1647  /* System generated locals */
1648  C_INT v_dim1, v_offset, i__1, i__2;
1649 
1650  /* Local variables */
1651  static C_INT i__, j, k;
1652  static C_FLOAT64 s;
1653  static C_INT ip1, nm1;
1654 
1655  /* ...SORTS THE ELEMENTS OF D(N) INTO DESCENDING ORDER AND MOVES THE */
1656  /* CORRESPONDING COLUMNS OF V(N,N). */
1657  /* M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM. */
1658  /* Parameter adjustments */
1659  v_dim1 = *m;
1660  v_offset = v_dim1 + 1;
1661  v -= v_offset;
1662  --d__;
1663 
1664  /* Function Body */
1665  if (*n == 1)
1666  {
1667  return 0;
1668  }
1669 
1670  nm1 = *n - 1;
1671  i__1 = nm1;
1672 
1673  for (i__ = 1; i__ <= i__1; ++i__)
1674  {
1675  k = i__;
1676  s = d__[i__];
1677  ip1 = i__ + 1;
1678  i__2 = *n;
1679 
1680  for (j = ip1; j <= i__2; ++j)
1681  {
1682  if (d__[j] <= s)
1683  {
1684  goto L1;
1685  }
1686 
1687  k = j;
1688  s = d__[j];
1689 L1:
1690  ;
1691  }
1692 
1693  if (k <= i__)
1694  {
1695  goto L3;
1696  }
1697 
1698  d__[k] = d__[i__];
1699  d__[i__] = s;
1700  i__2 = *n;
1701 
1702  for (j = 1; j <= i__2; ++j)
1703  {
1704  s = v[j + i__ * v_dim1];
1705  v[j + i__ * v_dim1] = v[j + k * v_dim1];
1706  /* L2: */
1707  v[j + k * v_dim1] = s;
1708  }
1709 
1710 L3:
1711  ;
1712  }
1713 
1714  return 0;
1715 } /* sort_ */
1716 
1717 /* Subroutine */ int CPraxis::quad_(C_INT *n, FPraxis *f, C_FLOAT64 *x, C_FLOAT64 *t,
1718  C_FLOAT64 *machep, C_FLOAT64 *h__)
1719 {
1720  /* System generated locals */
1721  C_INT i__1;
1722  C_FLOAT64 d__1;
1723 
1724  /* Local variables */
1725  static C_INT i__;
1726  static C_FLOAT64 l, s, value;
1727  /* ...QUAD LOOKS FOR THE MINIMUM OF F ALONG A CURVE DEFINED BY Q0,Q1,X...
1728  */
1729  /* Parameter adjustments */
1730  --x;
1731 
1732  /* Function Body */
1733  s = global_1.fx;
1734  global_1.fx = q_1.qf1;
1735  q_1.qf1 = s;
1736  q_1.qd1 = 0.;
1737  i__1 = *n;
1738 
1739  for (i__ = 1; i__ <= i__1; ++i__)
1740  {
1741  s = x[i__];
1742  l = q_1.q1[i__ - 1];
1743  x[i__] = l;
1744  q_1.q1[i__ - 1] = s;
1745  /* L1: */
1746  /* Computing 2nd power */
1747  d__1 = s - l;
1748  q_1.qd1 += d__1 * d__1;
1749  }
1750 
1751  q_1.qd1 = sqrt(q_1.qd1);
1752  l = q_1.qd1;
1753  s = 0.;
1754 
1755  if (q_1.qd0 <= 0. || q_1.qd1 <= 0. || global_1.nl < *n * 3 * *n)
1756  {
1757  goto L2;
1758  }
1759 
1760  value = q_1.qf1;
1761  min_(n, &c__0, &c__2, &s, &l, &value, &c_true, f, &x[1], t, machep,
1762  h__);
1763  q_1.qa = l * (l - q_1.qd1) / (q_1.qd0 * (q_1.qd0 + q_1.qd1));
1764  q_1.qb = (l + q_1.qd0) * (q_1.qd1 - l) / (q_1.qd0 * q_1.qd1);
1765  q_1.qc = l * (l + q_1.qd0) / (q_1.qd1 * (q_1.qd0 + q_1.qd1));
1766  goto L3;
1767 L2:
1768  global_1.fx = q_1.qf1;
1769  q_1.qa = 0.;
1770  q_1.qb = q_1.qa;
1771  q_1.qc = 1.;
1772 L3:
1773  q_1.qd0 = q_1.qd1;
1774  i__1 = *n;
1775 
1776  for (i__ = 1; i__ <= i__1; ++i__)
1777  {
1778  s = q_1.q0[i__ - 1];
1779  q_1.q0[i__ - 1] = x[i__];
1780  /* L4: */
1781  x[i__] = q_1.qa * s + q_1.qb * x[i__] + q_1.qc * q_1.q1[i__ - 1];
1782  }
1783 
1784  return 0;
1785 } /* quad_ */
1786 
1787 /* Subroutine */ int vcprnt_(C_INT *option, C_FLOAT64 *v, C_INT *n)
1788 {
1789  /* System generated locals */
1790  C_INT i__1;
1791 
1792  /* Local variables */
1793  C_INT i;
1794 
1795  /* Parameter adjustments */
1796  --v;
1797 
1798  /* Function Body */
1799  switch ((int)*option)
1800  {
1801  case 1: goto L1;
1802  case 2: goto L2;
1803  case 3: goto L3;
1804  case 4: goto L4;
1805  }
1806 
1807 L1:
1808  printf("THE SECOND DIFFERENCE ARRAY D[*] IS :\n");
1809  i__1 = *n;
1810 
1811  for (i = 1; i <= i__1; ++i)
1812  {
1813  printf("%g\n", v[i]);
1814  }
1815 
1816  return 0;
1817 L2:
1818  printf("THE SCALE FACTORS ARE:\n");
1819  i__1 = *n;
1820 
1821  for (i = 1; i <= i__1; ++i)
1822  {
1823  printf("%g\n", v[i]);
1824  }
1825 
1826  return 0;
1827 L3:
1828  printf("THE APPROXIMATING QUADRATIC FORM HAS THE PRINCEPAL VALUES:\n");
1829  i__1 = *n;
1830 
1831  for (i = 1; i <= i__1; ++i)
1832  {
1833  printf("%g\n", v[i]);
1834  }
1835 
1836  return 0;
1837 L4:
1838  printf("x is:\n");
1839  i__1 = *n;
1840 
1841  for (i = 1; i <= i__1; ++i)
1842  {
1843  printf("%g\n", v[i]);
1844  }
1845 
1846  return 0;
1847 } /* vcprnt_ */
1848 
1849 /* Subroutine */ int CPraxis::print_(C_INT * /* n */ , C_FLOAT64 * /* x */ , C_INT * /* prin */,
1850  C_FLOAT64 * /* fmin */)
1851 {
1852 #ifdef XXXX
1853  /* System generated locals */
1854  C_INT i__1;
1855  C_FLOAT64 d__1;
1856 
1857  /* Local variables */
1858  C_INT i;
1859  C_FLOAT64 ln = 0.0;
1860 
1861  /* Builtin functions */
1862  C_FLOAT64 d_lg10(C_FLOAT64 *);
1863  /* Parameter adjustments */
1864  --x;
1865 
1866  /* Function Body */
1867  printf("AFTER ");
1868  printf("%d", global_1.nl);
1869  printf(" LINEAR SEARCHES, THE FUNCTION HAS BEEN EVALUATED ");
1870  printf("%d TIMES.\n", global_1.nf);
1871  printf("THE SMALLEST VALUE FOUND IS F(X) = ");
1872  printf("%g\n", global_1.fx);
1873 
1874  if (global_1.fx <= *fmin)
1875  {
1876  goto L1;
1877  }
1878 
1879  d__1 = global_1.fx - *fmin;
1880  // ln = d_lg10(&d__1);
1881  printf("log (f(x)) - ");
1882  printf("%g", *fmin);
1883  printf(" = ");
1884  printf("%g\n", ln);
1885  goto L2;
1886 L1:
1887  printf("LOG (F(x)) -- ");
1888  printf("%g", *fmin);
1889  printf(" IS UNDERFINED\n");
1890 L2:
1891 
1892  if (*n > 4 && *prin <= 2)
1893  {
1894  return 0;
1895  }
1896 
1897  i__1 = *n;
1898 
1899  for (i = 1; i <= i__1; ++i)
1900  {
1901  printf("x is:");
1902  printf("%g\n", x[i]);
1903  }
1904 
1905 #endif // XXXX
1906 
1907  return 0;
1908 } /* print_ */
1909 
1910 /* Subroutine */ int maprnt_(C_INT * /* option */ , C_FLOAT64 * /* v */, C_INT * /* m */,
1911  C_INT * /* n */)
1912 {
1913 #ifdef XXXX
1914  /* System generated locals */
1915  C_INT v_dim1, v_offset, i__1, i__2;
1916 
1917  /* Local variables */
1918  C_INT i, j, low, upp;
1919 
1920  /* ...THE SUBROUTINE MAPRNT PRINTS THE COLUMNS OF THE NXN MATRIX V */
1921  /* WITH A HEADING AS SPECIFIED BY OPTION. */
1922  /* M IS THE ROW DIMENSION OF V AS DECLARED IN THE CALLING PROGRAM... */
1923  /* Parameter adjustments */
1924  v_dim1 = *m;
1925  v_offset = v_dim1 + 1;
1926  v -= v_offset;
1927 
1928  /* Function Body */
1929  low = 1;
1930  upp = 5;
1931 
1932  switch ((int)*option)
1933  {
1934  case 1: goto L1;
1935  case 2: goto L2;
1936  }
1937 
1938 L1:
1939  printf("HE NEW DIRECTIONS ARE:\n");
1940  goto L3;
1941 L2:
1942  printf("AND THE PRINCIPAL AXES:\n");
1943 L3:
1944 
1945  if (*n < upp)
1946  {
1947  upp = *n;
1948  }
1949 
1950  i__1 = *n;
1951 
1952  for (i = 1; i <= i__1; ++i)
1953  {
1954  /* L4: */
1955  printf("%3d", i);
1956  i__2 = upp;
1957 
1958  for (j = low; j <= i__2; ++j)
1959  {
1960  printf(" %12g", v[i*v_dim1 + j]);
1961  }
1962 
1963  printf("\n");
1964  }
1965 
1966  low += 5;
1967 
1968  if (*n < low)
1969  {
1970  return 0;
1971  }
1972 
1973  upp += 5;
1974  goto L3;
1975 
1976 #endif // XXXX
1977 
1978  return 0;
1979 } /* maprnt_ */
C_FLOAT64 praxis_(C_FLOAT64 *t0, C_FLOAT64 *machep, C_FLOAT64 *h0, C_INT *n, C_INT *prin, C_FLOAT64 *x, FPraxis *f, C_FLOAT64 *fmin)
Definition: CPraxis.cpp:69
int maprnt_(C_INT *, C_FLOAT64 *, C_INT *, C_INT *)
Definition: CPraxis.cpp:1910
#define C_INT
Definition: copasi.h:115
static bool c_false
Definition: CPraxis.cpp:52
#define pdelete(p)
Definition: copasi.h:215
C_FLOAT64 qd1
Definition: CPraxis.h:73
#define FALSE_
Definition: CPraxis.cpp:44
~CPraxis()
Definition: CPraxis.cpp:64
#define dmin(a1, a2)
Definition: CPraxis.cpp:46
#define dmax(a1, a2)
Definition: CPraxis.cpp:45
static bool c_true
Definition: CPraxis.cpp:54
int sort_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *)
Definition: CPraxis.cpp:1644
CRandom * mpRandom
Definition: CPraxis.h:81
int minfit_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
Definition: CPraxis.cpp:800
C_FLOAT64 fx
Definition: CPraxis.h:67
int vcprnt_(C_INT *, C_FLOAT64 *, C_INT *)
Definition: CPraxis.cpp:1787
Q q_1
Definition: CPraxis.h:84
C_FLOAT64 q1[100]
Definition: CPraxis.h:73
static CRandom * createGenerator(CRandom::Type type=CRandom::mt19937, unsigned C_INT32 seed=0)
Definition: CRandom.cpp:49
virtual C_FLOAT64 getRandomCC()
Definition: CRandom.cpp:235
C_FLOAT64 v[10000]
Definition: CPraxis.h:73
int quad_(C_INT *, FPraxis *f, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
Definition: CPraxis.cpp:1717
int min_(C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, bool *, FPraxis *f, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
Definition: CPraxis.cpp:1329
C_FLOAT64 qc
Definition: CPraxis.h:73
C_FLOAT64 dmin__
Definition: CPraxis.h:67
static C_INT c__3
Definition: CPraxis.cpp:55
int print_(C_INT *n, C_FLOAT64 *x, C_INT *prin, C_FLOAT64 *fmin)
Definition: CPraxis.cpp:1849
#define TRUE_
Definition: CPraxis.cpp:43
C_FLOAT64 qb
Definition: CPraxis.h:73
#define C_FLOAT64
Definition: copasi.h:92
C_FLOAT64 qd0
Definition: CPraxis.h:73
static C_INT c__4
Definition: CPraxis.cpp:53
static C_INT c__0
Definition: CPraxis.cpp:56
C_FLOAT64 flin_(C_INT *, C_INT *, C_FLOAT64 *, FPraxis *, C_FLOAT64 *, C_INT *)
Definition: CPraxis.cpp:1591
C_FLOAT64 q0[100]
Definition: CPraxis.h:73
C_FLOAT64 ldt
Definition: CPraxis.h:67
static C_INT c__2
Definition: CPraxis.cpp:51
Global global_1
Definition: CPraxis.h:83
CPraxis()
Definition: CPraxis.cpp:58
static C_INT c__1
Definition: CPraxis.cpp:50
C_FLOAT64 qf1
Definition: CPraxis.h:73
C_FLOAT64 qa
Definition: CPraxis.h:73