home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / msdos / math / praxis.arc / PRAXIS.INC < prev    next >
Text File  |  1987-07-15  |  24KB  |  565 lines
  1. (****************************************************************************)
  2. (*              P R O C E D U R E      P R A X I S                          *)
  3. (*                                                                          *)
  4. (*   P R A X I S  is for the minimization of a function in several          *)
  5. (*   variables. The algorithm used is a modification of a conjugate         *)
  6. (*   gradient method developed by Powell. Changes are due to Brent,         *)
  7. (*   who gives an ALGOL W program.                                          *)
  8. (*   Users who are interested in more of the details should read            *)
  9. (*          - Powell, M. J. D., 1962. An efficient method for finding       *)
  10. (*            the minimum of a function of several variables without        *)
  11. (*            calculating derivatives, Computer Journal 7, 155-162.         *)
  12. (*          - Brent, R. P., 1973. Algorithms for minimization without       *)
  13. (*            derivatives. Prentice Hall, Englewood Cliffs.                 *)
  14. (*   If you have any further comments, questions, etc. please feel free     *)
  15. (*   and contact me.                                                        *)
  16. (*                                                                          *)
  17. (*                           Karl Gegenfurtner     02/17/86                 *)
  18. (*                                                 TURBO - P Version 1.0    *)
  19. (****************************************************************************)
  20. (*  The use of PRAXIS is fairly simple. There are only two parameters:      *)
  21. (*      - X is of type vector and contains on input an initial guess        *)
  22. (*        of the solution. on output X holds the final solution of the      *)
  23. (*        system of equations.                                              *)
  24. (*      - N is of type integer and gives the number of unknown parameters   *)
  25. (*  The result of PRAXIS is the least calculated value of the function F    *)
  26. (*  F is one of the global parameters used by PRAXIS:                       *)
  27. (*    - F(X,N) is declared forward and is the function to be minimized      *)
  28. (*  All other globals are optional, i.e. you can change them or else        *)
  29. (*  PRAXIS will use some default values, which are adequate for most        *)
  30. (*  problems under consideration.                                           *)
  31. (*    - Prin controls the printout from PRAXIS                              *)
  32. (*           0:  no printout at all                                         *)
  33. (*           1:  only initial and final values                              *)
  34. (*           2:  detailed map of the minimization process                   *)
  35. (*           3:  also prints eigenvalues and eigenvectors of the            *)
  36. (*               direction matices in use (for insiders only).              *)
  37. (*    - T is the tolerance for the precision of the solution                *)
  38. (*      PRAXIS returns if the criterion                                     *)
  39. (*             2 * ||x(k)-x(k-1)|| <= Sqrt(MachEps) * ||x(k)|| + T          *)
  40. (*             is fulfilled more than KTM times, where                      *)
  41. (*    - KTM is another global parameter. It's default value is 1 and        *)
  42. (*      a value of 4 leads to a very(!) cautious stopping criterion         *)
  43. (*    - MachEps is the relative machine precision and is                    *)
  44. (*      1.0 E-15 using an 8087 and TURBO-87 and                             *)
  45. (*      1.0 E-06 without 8087.                                              *)
  46. (*    - H is a steplength parameter and shoul be set to the expected        *)
  47. (*      distance to the solution. An exceptional large or small value       *)
  48. (*      of H leads to slower convergence on the first few iterations.       *)
  49. (*    - Scbd is a scaling parameter and should be set to about 10.          *)
  50. (*      The default is 1 and with that value no scaling is done at all      *)
  51. (*      It's only necessary when the parameters are scaled very different   *)
  52. (*    - IllC is a Boolean variable and should be set to TRUE if the         *)
  53. (*      problem is known to be illconditioned.                              *)
  54. (****************************************************************************)
  55.  
  56.  
  57. CONST MaxPar = 20;
  58.  
  59. TYPE  Vector = ARRAY[1..MaxPar] OF REAL;
  60.       Matrix = ARRAY[1..MaxPar] OF Vector;
  61.       PrString = STRING[80];
  62.  
  63. CONST Macheps: REAL = 1.0E-08;  (* set to 1.0E-13 for TURBO-87 *)
  64.       H:       REAL = 1.0E+00;
  65.       T:       REAL = 1.0E-05;
  66.       Prin:    INTEGER = 2;
  67.       IllC:    BOOLEAN = FALSE;
  68.       Scbd:    REAL = 1;
  69.       Ktm:     INTEGER = 2;
  70.  
  71.  
  72. FUNCTION F(VAR x:Vector;n:INTEGER):REAL; FORWARD;
  73.  
  74. FUNCTION Power(a,b: REAL):REAL;
  75. BEGIN    Power := Exp(b*Ln(a));
  76. END;     (* Power *)
  77.  
  78.  
  79. PROCEDURE MinFit(N:INTEGER;Eps,Tol:REAL;VAR ab:Matrix;VAR q:Vector);
  80. LABEL     TestFsplitting,
  81.           TestFconvergence,
  82.           Convergence,
  83.           Cancellation;
  84. VAR       l, kt,
  85.           l2,
  86.           i, j, k: INTEGER;
  87.           c, f, g,
  88.           h, s, x,
  89.           y, z:    REAL;
  90.           e:       Vector;
  91. BEGIN     (* Householders reduction to bidiagonal form *)
  92.           x:= 0; g:= 0;
  93.           FOR i:= 1 TO N DO
  94.           BEGIN e[i]:= g; s:= 0; l:= i+1;
  95.                 FOR j:= i TO N DO
  96.                     s:= s + Sqr(ab[j,i]);
  97.                 IF s<Tol THEN g:= 0 ELSE
  98.                 BEGIN f:= ab[i,i];
  99.                       IF f<0
  100.                          THEN g:= Sqrt(s)
  101.                          ELSE g:= - Sqrt(s);
  102.                       h:= f*g-s;  ab[i,i]:= f - g;
  103.                       FOR j:= l TO N DO
  104.                       BEGIN f:= 0;
  105.                             FOR k:= i TO N DO
  106.                                 f:= f + ab[k,i]*ab[k,j];
  107.                             f:= f/h;
  108.                             FOR k:= i TO N DO
  109.                                 ab[k,j]:= ab[k,j] + f*ab[k,i];
  110.                       END; (* j *)
  111.                 END; (* IF *)
  112.                 q[i]:= g; s:= 0;
  113.                 IF i<=N
  114.                    THEN FOR j:= l TO N DO
  115.                             s:= s + Sqr(ab[i,j]);
  116.                 IF s<Tol THEN g:= 0 ELSE
  117.                 BEGIN f:= ab[i,i+1];
  118.                       IF f<0
  119.                          THEN g:= Sqrt(s)
  120.                          ELSE g:= - Sqrt(s);
  121.                       h:= f*g-s;  ab[i,i+1]:= f-g;
  122.                       FOR j:= l TO N DO e[j]:= ab[i,j]/h;
  123.                       FOR j:= l TO N DO
  124.                       BEGIN s:= 0;
  125.                             FOR k:= l TO N DO s:= s + ab[j,k]*ab[i,k];
  126.                             FOR k:= l TO N DO ab[j,k]:= ab[j,k] + s*e[k];
  127.                       END; (* J *)
  128.                 END; (* IF *)
  129.                 y:= Abs(q[i])+Abs(e[i]);
  130.                 IF y > x THEN x:= y;
  131.           END; (* I *)
  132.           (* accumulation of right hand transformations *)
  133.           FOR i:= N DOWNTO 1 DO
  134.           BEGIN IF g<>0.0 THEN
  135.                 BEGIN h:= ab[i,i+1]*g;
  136.                       FOR j:= l TO N DO ab[j,i]:= ab[i,j]/h;
  137.                       FOR j:= l TO N DO
  138.                       BEGIN s:= 0;
  139.                             FOR k:= l TO N DO s:= s + ab[i,k]*ab[k,j];
  140.                             FOR k:= l TO N DO ab[k,j]:= ab[k,j] + s*ab[k,i];
  141.                       END; (* J *)
  142.                 END; (* IF *)
  143.                 FOR j:= l TO N DO
  144.                 BEGIN ab[j,i]:= 0;
  145.                       ab[i,j]:= 0;
  146.                 END;
  147.                 ab[i,i]:= 1; g:= e[i]; l:= i;
  148.           END; (* I *)
  149.           (* diagonalization to bidiagonal form *)
  150.           eps:= eps*x;
  151.           FOR k:= N DOWNTO 1 DO
  152.           BEGIN kt:= 0;
  153. TestFsplitting: kt:= kt + 1;
  154.                 IF kt > 30 THEN
  155.                 BEGIN e[k]:= 0;
  156.                       WriteLn('+++ QR failed');
  157.                 END;
  158.                 FOR l2:= k DOWNTO 1 DO
  159.                 BEGIN l:= l2;
  160.                       IF Abs(e[l])<=Eps THEN Goto TestFconvergence;
  161.                       IF Abs(q[l-1])<=Eps THEN Goto Cancellation;
  162.                 END; (* L2 *)
  163. Cancellation:   c:= 0; s:= 1;
  164.                 FOR i:= l TO k DO
  165.                 BEGIN f:= s*e[i]; e[i]:= c*e[i];
  166.                       IF Abs(f)<=Eps THEN Goto TestFconvergence;
  167.                       g:= q[i];
  168.                       IF Abs(f) < Abs(g)
  169.                          THEN h:= Abs(g)*Sqrt(1+sqr(f/g))
  170.                          ELSE IF f <> 0
  171.                                  THEN h:= Abs(f)*Sqrt(1+Sqr(g/f))
  172.                                  ELSE h:= 0;
  173.                       q[i]:= h;
  174.                       IF h = 0 THEN
  175.                       BEGIN h:= 1;
  176.                             g:= 1;
  177.                       END;
  178.                       c:= g/h; s:= -f/h;
  179.                 END; (* I *)
  180. TestFconvergence:
  181.                 z:= q[k];
  182.                 IF l=k THEN Goto convergence;
  183.                 (* shift from bottom 2*2 minor *)
  184.                 x:= q[l]; y:= q[k-1]; g:= e[k-1];  h:= e[k];
  185.                 f:= ((y-z)*(y+z) + (g-h)*(g+h))/(2*h*y);
  186.                 g:= Sqrt(Sqr(f)+1);
  187.                 IF f<=0
  188.                    THEN f:= ((x-z)*(x+z)+h*(y/(f-g)-h))/x
  189.                    ELSE f:= ((x-z)*(x+z)+h*(y/(f+g)-h))/x;
  190.                 (* next QR transformation *)
  191.                 s:= 1; c:= 1;
  192.                 FOR i:= l+1 TO k DO
  193.                 BEGIN g:= e[i]; y:= q[i]; h:= s*g; g:= g*c;
  194.                       IF Abs(f)<Abs(h)
  195.                          THEN z:= Abs(h)*Sqrt(1+Sqr(f/h))
  196.                          ELSE IF f<>0
  197.                                  THEN z:= Abs(f)*Sqrt(1+Sqr(h/f))
  198.                                  ELSE z:= 0;
  199.                       e[i-1]:= z;
  200.                       IF z=0 THEN
  201.                       BEGIN f:= 1;
  202.                             z:= 1;
  203.                       END;
  204.                       c:= f/z;  s:= h/z;
  205.                       f:= x*c + g*s;  g:= - x*s + g*c; h:= y*s;
  206.                       y:= y*c;
  207.                       FOR j:= 1 TO N DO
  208.                       BEGIN x:= ab[j,i-1]; z:= ab[j,i];
  209.                             ab[j,i-1]:= x*c + z*s;
  210.                             ab[j,i]:= - x*s + z*c;
  211.                       END;
  212.                       IF Abs(f)<Abs(h)
  213.                          THEN z:= Abs(h)*Sqrt(1+Sqr(f/h))
  214.                          ELSE IF f<>0
  215.                                  THEN z:= Abs(f)*Sqrt(1+Sqr(h/f))
  216.                                  ELSE z:= 0;
  217.                       q[i-1]:= z;
  218.                       IF z=0 THEN
  219.                       BEGIN z:= 1;
  220.                             f:= 1;
  221.                       END;
  222.                       c:= f/z; s:= h/z;
  223.                       f:= c*g + s*y; x:= -s*g + c*y;
  224.                 END; (* I *)
  225.                 e[l]:= 0; e[k]:= f; q[k]:= x;
  226.                 Goto TestFsplitting;
  227. Convergence:    IF z<0 THEN
  228.                 BEGIN q[k]:= -z;
  229.                       FOR j:= 1 TO N DO ab[j,k]:= - ab[j,k];
  230.                 END;
  231.           END; (* K *)
  232. END;      (* Minfit *)
  233.  
  234. FUNCTION Praxis(VAR X:Vector; N:INTEGER):REAL;
  235. LABEL    l0, l1, l2;
  236. VAR      i, j,
  237.          k, k2,
  238.          nl, nf, kl, kt: INTEGER;
  239.          s, sl, dn, dmin,
  240.          fx, f1,
  241.          lds, ldt, sf, df,
  242.          qf1, qd0, qd1,
  243.          qa, qb, qc,
  244.          m2, m4,
  245.          small, vsmall,
  246.          large, vlarge,
  247.          ldfac, t2:      REAL;
  248.          d, y, z,
  249.          q0, q1:         Vector;
  250.          v:              Matrix;
  251.  
  252. PROCEDURE SORT;                         (* sort d and v in descending order *)
  253. VAR       k, i, j:       INTEGER;
  254.           s:             REAL;
  255. BEGIN     FOR i:= 1 TO N-1 DO
  256.           BEGIN k:= i; s:= d[i];
  257.                 FOR j:= i+1 TO N DO
  258.                     IF d[j] > s THEN
  259.                     BEGIN k:= j; s:= d[j];
  260.                     END;
  261.                 IF k>i THEN
  262.                 BEGIN d[k]:= d[i]; d[i]:= s;
  263.                       FOR j:= 1 TO N DO
  264.                       BEGIN s:= v[j,i];
  265.                             v[j,i]:= v[j,k];
  266.                             v[j,k]:= s;
  267.                       END;
  268.                 END; (* IF *)
  269.           END; (* FOR *)
  270. END;      (* sort *)
  271.  
  272. PROCEDURE Print;                                   (* print a line of traces *)
  273. VAR       i:             INTEGER;
  274. BEGIN     WriteLn('------------------------------------------------------');
  275.           WriteLn('Chi Square reduced to ... ', fx);
  276.           WriteLn(' ... after ',nf,' function calls ...');
  277.           WriteLn(' ... including ',nl,' linear searches.');
  278.           WriteLn('Current Values of X ...');
  279.           FOR i:= 1 TO N DO
  280.               Write(X[i]);
  281.           WriteLn;
  282. END;      (* print *)
  283.  
  284. PROCEDURE MatPrint(s:PrString;v:Matrix;n,m:INTEGER);
  285. VAR       k, i:          INTEGER;
  286. BEGIN     WriteLn;
  287.           WriteLn(s);
  288.           FOR k:= 1 TO N DO
  289.           BEGIN FOR i:= 1 TO m DO
  290.                     Write(v[k,i]);
  291.                 WriteLn;
  292.           END;
  293.           WriteLn;
  294. END;      (* MatPrint *)
  295.  
  296. PROCEDURE VecPrint(s:PrString;v:Vector;n:INTEGER);
  297. VAR       i:             INTEGER;
  298. BEGIN     WriteLn;
  299.           WriteLn(s);
  300.           FOR i:= 1 TO n DO
  301.               Write(v[i]);
  302.           WriteLn;
  303. END;      (* VecPrint *)
  304.  
  305. PROCEDURE Min(j, Nits:INTEGER; VAR d2, x1:REAL; f1:REAL;fk:BOOLEAN);
  306. LABEL     l0, l1;
  307. VAR       k, i:          INTEGER;
  308.           dz:            BOOLEAN;
  309.           x2, xm, f0,
  310.           f2, fm, d1, t2,
  311.           s, sf1, sx1:   REAL;
  312.  FUNCTION Flin(l:REAL):REAL;
  313.  VAR      i:             INTEGER;
  314.           t:             Vector;
  315.  BEGIN    IF j>0 THEN                     (* linear search *)
  316.              FOR i:= 1 TO N DO
  317.                  t[i]:= x[i]+l*v[i,j]
  318.              ELSE BEGIN              (* search along parabolic space curve *)
  319.                   qa:= l*(l-qd1)/(qd0*(qd0+qd1));
  320.                   qb:= (l+qd0)*(qd1-l)/(qd0*qd1);
  321.                   qc:= l*(l+qd0)/(qd1*(qd0+qd1));
  322.                   FOR i:= 1 TO N DO
  323.                       t[i]:= qa*q0[i]+qb*x[i]+qc*q1[i];
  324.              END; (* ELSE *)
  325.           nf:= nf+1;
  326.           Flin:= F(t, N);
  327.  END;     (* Flin *)
  328. BEGIN     (* Min *)
  329.           sf1:= f1; sx1:= x1;
  330.           k:= 0; xm:= 0; fm:= fx; f0:= fx; dz:= (d2<MachEps);
  331.           (* Find step size *)
  332.           s:= 0; FOR i:= 1 TO N DO s:= s + Sqr(X[i]);
  333.           s:= Sqrt(s);
  334.           IF dz
  335.              THEN t2:= m4*Sqrt(Abs(fx)/Dmin + s*ldt) + m2*ldt
  336.              ELSE t2:= m4*Sqrt(Abs(fx)/D2 + s*ldt) + m2*ldt;
  337.           s:= m4*s + t;
  338.           IF dz AND (t2>s) THEN t2:= s;
  339.           IF (t2<small) THEN t2:= small;
  340.           IF (t2>(0.01*h)) THEN t2:= 0.01*h;
  341.           IF fk AND (f1<=fm) THEN BEGIN xm:= x1; fm:= f1; END;
  342.           IF NOT fk OR (Abs(x1)<t2) THEN
  343.           BEGIN IF x1>0 THEN x1:= t2 ELSE x1:= - t2;
  344.                 f1:= Flin(x1);
  345.           END;
  346.           IF f1<=fm THEN BEGIN xm:= x1; fm:= f1; END;
  347. l0:       IF dz THEN
  348.           BEGIN IF f0<f1 THEN x2:= - x1 ELSE x2:= 2*x1;
  349.                 f2:= Flin(x2);
  350.                 IF f2 <= fm THEN BEGIN xm:= x2; fm:= f2; END;
  351.                 d2:= (x2*(f1-f0) - x1*(f2-f0))/(x1*x2*(x1-x2));
  352.           END;
  353.           d1:= (f1-f0)/x1 - x1*d2; dz:= TRUE;
  354.           IF d2 <= small
  355.              THEN IF d1<0
  356.                      THEN x2:= h
  357.                      ELSE x2:= -h
  358.              ELSE x2:= -0.5*d1/d2;
  359.           IF Abs(x2)>h
  360.              THEN IF x2>0
  361.                      THEN x2:= h
  362.                      ELSE x2:= -h;
  363. l1:       f2:= Flin(x2);
  364.           IF (k<Nits) AND (f2>f0) THEN
  365.           BEGIN k:= k + 1;
  366.                 IF (f0<f1) AND ((x1*x2)>0) THEN GOTO l0;
  367.                 x2:= 0.5*x2; GOTO l1;
  368.           END;
  369.           nl:= nl + 1;
  370.           IF f2>fm THEN x2:= xm ELSE fm:= f2;
  371.           IF Abs(x2*(x2-x1))>small
  372.              THEN d2:= (x2*(f1-f0) - x1*(fm-f0))/(x1*x2*(x1-x2))
  373.              ELSE IF k>0
  374.                      THEN d2:= 0;
  375.           IF d2<=small THEN d2:= small;
  376.           x1:= x2; fx:= fm;
  377.           IF sf1<fx THEN BEGIN fx:= sf1; x1:= sx1; END;
  378.           IF j>0
  379.              THEN FOR i:= 1 TO N DO
  380.                       x[i]:= x[i] + x1*v[i,j];
  381. END;      (* Min *)
  382.  
  383. PROCEDURE Quad;        (* look for a minimum along the curve q0, q1, q2 *)
  384. VAR       i:             INTEGER;
  385.           l, s:          REAL;
  386. BEGIN     s:= fx; fx:= qf1; qf1:= s; qd1:= 0;
  387.           FOR i:= 1 TO N DO
  388.           BEGIN s:= x[i];  l:= q1[i];  x[i]:= l;  q1[i]:= s;
  389.                 qd1:= qd1 + Sqr(s - l);
  390.           END;
  391.           s:= 0; qd1:= Sqrt(qd1); l:= qd1;
  392.           IF (qd0>0) AND (qd1>0) AND (nl >= (3*N*N)) THEN
  393.           BEGIN Min(0, 2, s, l, qf1, TRUE);
  394.                 qa:= l*(l-qd1)/(qd0*(qd0+qd1));
  395.                 qb:= (l+qd0)*(qd1-l)/(qd0*qd1);
  396.                 qc:= l*(l+qd0)/(qd1*(qd0+qd1));
  397.           END ELSE
  398.           BEGIN fx:= qf1; qa:= 0; qb:= 0; qc:= 1;
  399.           END;
  400.           qd0:= qd1;
  401.           FOR i:= 1 TO N DO
  402.           BEGIN s:= q0[i];  q0[i]:= x[i];
  403.                 x[i]:= qa*s + qb*x[i] + qc*q1[i];
  404.           END;
  405. END;      (* quad *)
  406.  
  407. BEGIN     (****  P R A X I S  ****)
  408.           (* initialization *)
  409.           small:= Sqr(MachEps); vsmall:= Sqr(small);
  410.           large:= 1.0/small;    vlarge:= 1.0/vsmall;
  411.           m2:= Sqrt(MachEps);   m4:= Sqrt(m2);
  412.           IF IllC THEN ldfac:= 0.1 ELSE ldfac:= 0.01;
  413.           nl:= 0; kt:= 0; nf:= 1; fx:= f(X, N); qf1:= fx;
  414.           t2:= small + Abs(t); t:= t2; dmin:= small;
  415.           IF h<(100*t) THEN h:= 100*t; ldt:= h;
  416.           FOR i:= 1 TO N DO FOR j:= 1 TO N DO
  417.               IF i=j THEN v[i,j]:= 1 ELSE v[i,j]:= 0;
  418.           d[1]:= 0; qd0:= 0;
  419.           FOR i:= 1 TO N DO q1[i]:= x[i];
  420.           IF Prin > 1 THEN
  421.           BEGIN WriteLn;WriteLn('---------- enter function praxis ------------');
  422.                 WriteLn('Current paramter settings are:');
  423.                 Writeln('... scaling ... ',scbd);
  424.                 WriteLn('... MachEps ... ',MachEps);
  425.                 WriteLn('...   Tol   ... ',t);
  426.                 WriteLn('... MaxStep ... ',h);
  427.                 WriteLn('...   IllC  ... ',IllC);
  428.                 WriteLn('...   Ktm   ... ',ktm);
  429.           END;
  430.           IF Prin > 0 THEN Print;
  431.  
  432. l0:       (* main loop *)
  433.           sf:= d[1]; s:= 0; d[1]:= 0;
  434.           (* minimize along first direction *)
  435.           Min(1, 2, d[1], s, fx, FALSE);
  436.           IF s<= 0
  437.              THEN FOR i:= 1 TO N DO
  438.                       v[i,1]:= - v[i,1];
  439.           IF (sf<= (0.9*d[1])) OR ((0.9*sf)>=d[1])
  440.              THEN FOR i:= 2 TO N DO
  441.                       d[i]:= 0;
  442.           FOR k:= 2 TO N DO
  443.           BEGIN FOR i:= 1 TO N DO y[i]:= x[i];
  444.                 sf:= fx;
  445.                 IllC:= IllC OR (kt>0);
  446. l1:             kl:= k; df:= 0;
  447.                 IF IllC THEN   (* random step to get off resolution valley *)
  448.                 BEGIN FOR i:= 1 TO N DO
  449.                       BEGIN z[i]:= (0.1*ldt + t2*Power(10,kt))*(Random-0.5);
  450.                             s:= z[i];
  451.                             FOR j:= 1 TO N DO x[j]:= x[j]+s*v[j,i];
  452.                       END; (* I *)
  453.                       fx:= F(x, N);   nf:= nf + 1;
  454.                 END; (* IF *)
  455.                 FOR k2:= k TO N DO  (* minimize along non-conjugate directions *)
  456.                 BEGIN sl:= fx; s:= 0;
  457.                       Min(k2, 2, d[k2], s, fx, FALSE);
  458.                       IF IllC
  459.                          THEN s:= d[k2]*Sqr(s+z[k2])
  460.                          ELSE s:= sl - fx;
  461.                       IF df<s THEN BEGIN df:= s;  kl:= k2; END;
  462.                 END; (* K2 *)
  463.                 IF NOT IllC AND (df < Abs(100*MachEps*fx)) THEN
  464.                 BEGIN IllC:= TRUE; GOTO l1;
  465.                 END;
  466.                 IF (k=2) AND (Prin>1) THEN VecPrint('New Direction ...', d, N);
  467.                 FOR k2:= 1 TO k-1 DO (* minimize along conjugate directions *)
  468.                 BEGIN s:= 0;
  469.                       Min(k2, 2, d[k2], s, fx, FALSE);
  470.                 END; (* K2 *)
  471.                 f1:= fx; fx:= sf; lds:= 0;
  472.                 FOR i:= 1 TO N DO
  473.                 BEGIN sl:= x[i]; x[i]:= y[i]; y[i]:= sl - y[i]; sl:= y[i];
  474.                       lds:= lds + Sqr(sl);
  475.                 END;
  476.                 lds:= Sqrt(lds);
  477.                 IF lds>small THEN
  478.                 BEGIN FOR i:= kl-1 DOWNTO k DO
  479.                       BEGIN FOR j:= 1 TO N DO v[j,i+1]:= v[j,i];
  480.                             d[i+1]:= d[i];
  481.                       END;
  482.                       d[k]:= 0;
  483.                       FOR i:= 1 TO N DO v[i,k]:= y[i]/lds;
  484.                       Min(k, 4, d[k], lds, f1, TRUE);
  485.                       IF lds<=0 THEN
  486.                       BEGIN lds:= -lds;
  487.                             FOR i:= 1 TO N DO v[i,k]:= -v[i,k];
  488.                       END;
  489.                 END; (* IF *)
  490.                 ldt:= ldfac*ldt; IF ldt<lds THEN ldt:= lds;
  491.                 IF Prin > 1 THEN Print;
  492.                 t2:= 0; FOR i:= 1 TO N DO t2:= t2 + Sqr(x[i]);
  493.                 t2:= m2*Sqrt(t2) + t;
  494.                 IF ldt > (0.5*t2) THEN kt:= 0 ELSE kt:= kt+1;
  495.                 IF kt>ktm THEN GOTO l2;
  496.           END; (* K *)
  497.           (*  try quadratic extrapolation in case    *)
  498.           (*  we are stuck in a curved valley        *)
  499.           Quad;
  500.           dn:= 0;
  501.           FOR i:= 1 TO N DO
  502.           BEGIN d[i]:= 1.0/Sqrt(d[i]);
  503.                 IF dn<d[i] THEN dn:= d[i];
  504.           END;
  505.           IF Prin>2 THEN MatPrint('New Matrix of Directions ...', v, n, n);
  506.           FOR j:= 1 TO N DO
  507.           BEGIN s:= d[j]/dn;
  508.                 FOR i:= 1 TO N DO v[i,j]:= s*v[i,j];
  509.           END;
  510.           IF scbd > 1 THEN  (* scale axis to reduce condition number *)
  511.           BEGIN s:= vlarge;
  512.                 FOR i:= 1 TO N DO
  513.                 BEGIN sl:= 0;
  514.                       FOR j:= 1 TO N DO sl:= sl + Sqr(v[i,j]);
  515.                       z[i]:= Sqrt(sl);
  516.                       IF z[i]<m4 THEN z[i]:= m4;
  517.                       IF s>z[i] THEN s:= z[i];
  518.                 END; (* I *)
  519.                 FOR i:= 1 TO N DO
  520.                 BEGIN sl:= s/z[i];  z[i]:= 1.0/sl;
  521.                       IF z[i] > scbd THEN
  522.                       BEGIN sl:= 1/scbd;
  523.                             z[i]:= scbd;
  524.                       END;
  525.                 END;
  526.           END; (* IF *)
  527.           FOR i:= 2 TO N DO  FOR j:= 1 TO i-1 DO
  528.           BEGIN s:= v[i,j]; v[i,j]:= v[j,i]; v[j,i]:= s;
  529.           END;
  530.           MinFit(N, MachEps, Vsmall, v, d);
  531.           IF scbd>1 THEN
  532.           BEGIN FOR i:= 1 TO N DO
  533.                 BEGIN s:= z[i];
  534.                       FOR j:= 1 TO n DO v[i,j]:= v[i,j]*s;
  535.                 END;
  536.                 FOR i:= 1 TO N DO
  537.                 BEGIN s:= 0;
  538.                       FOR j:= 1 TO N DO s:= s + Sqr(v[j,i]);
  539.                       s:= Sqrt(s);  d[i]:= s*d[i];  s:= 1.0/s;
  540.                       FOR j:= 1 TO N DO v[j,i]:= s*v[j,i];
  541.                 END;
  542.           END;
  543.           FOR i:= 1 TO N DO
  544.           BEGIN IF (dn*d[i])>Large
  545.                    THEN d[i]:= vsmall
  546.                    ELSE IF (dn*d[i])<small
  547.                            THEN d[i]:= Vlarge
  548.                            ELSE d[i]:= Power(dn*d[i],-2);
  549.           END;
  550.           Sort;   (* the new eigenvalues and eigenvectors *)
  551.           Dmin:= d[n];
  552.           IF Dmin < small THEN Dmin:= small;
  553.           IllC:= (m2*d[1]) > Dmin;
  554.           IF (Prin > 2) AND (scbd > 1)
  555.              THEN VecPrint('Scale Factors ...', z, n);
  556.           IF Prin > 2 THEN VecPrint('Eigenvalues of A ...', d, N);
  557.           IF Prin > 2 THEN MatPrint('Eigenvectors of A ...', v, n, n);
  558.           GOTO l0;  (* back to main loop *)
  559. l2:       IF Prin > 0 THEN
  560.           BEGIN VecPrint('Final solution is ...', x, n);
  561.                 WriteLn; WriteLn('ChiSq reduced to ',fx,' after ',nf, ' function calls.');
  562.           END;
  563.           Praxis:= fx;
  564. END;      (****  Praxis  ****)
  565.