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.I < prev    next >
Text File  |  1987-07-15  |  23KB  |  548 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 3 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. function power(a,b: real):real;
  58. begin    power := exp(b*ln(a));
  59. end;     (* power *)
  60.  
  61.  
  62. procedure minfit(n:integer;eps,tol:real;var ab:matrix;var q:vector);
  63. label     1,
  64.           2,
  65.           3,
  66.           4;
  67. var       l, kt,
  68.           l2,
  69.           i, j, k: integer;
  70.           c, f, g,
  71.           h, s, x,
  72.           y, z:    real;
  73.           e:       vector;
  74. begin     (* householders reduction to bidiagonal form *)
  75.           x:= 0; g:= 0;
  76.           for i:= 1 to n do
  77.           begin e[i]:= g; s:= 0; l:= i+1;
  78.                 for j:= i to n do
  79.                     s:= s + sqr(ab[j,i]);
  80.                 if s<tol then g:= 0 else
  81.                 begin f:= ab[i,i];
  82.                       if f<0
  83.                          then g:= sqrt(s)
  84.                          else g:= - sqrt(s);
  85.                       h:= f*g-s;  ab[i,i]:= f - g;
  86.                       for j:= l to n do
  87.                       begin f:= 0;
  88.                             for k:= i to n do
  89.                                 f:= f + ab[k,i]*ab[k,j];
  90.                             f:= f/h;
  91.                             for k:= i to n do
  92.                                 ab[k,j]:= ab[k,j] + f*ab[k,i];
  93.                       end; (* j *)
  94.                 end; (* if *)
  95.                 q[i]:= g; s:= 0;
  96.                 if i<=n
  97.                    then for j:= l to n do
  98.                             s:= s + sqr(ab[i,j]);
  99.                 if s<tol then g:= 0 else
  100.                 begin f:= ab[i,i+1];
  101.                       if f<0
  102.                          then g:= sqrt(s)
  103.                          else g:= - sqrt(s);
  104.                       h:= f*g-s;  ab[i,i+1]:= f-g;
  105.                       for j:= l to n do e[j]:= ab[i,j]/h;
  106.                       for j:= l to n do
  107.                       begin s:= 0;
  108.                             for k:= l to n do s:= s + ab[j,k]*ab[i,k];
  109.                             for k:= l to n do ab[j,k]:= ab[j,k] + s*e[k];
  110.                       end; (* j *)
  111.                 end; (* if *)
  112.                 y:= abs(q[i])+abs(e[i]);
  113.                 if y > x then x:= y;
  114.           end; (* i *)
  115.           (* accumulation of right hand transformations *)
  116.           for i:= n downto 1 do
  117.           begin if g<>0.0 then
  118.                 begin h:= ab[i,i+1]*g;
  119.                       for j:= l to n do ab[j,i]:= ab[i,j]/h;
  120.                       for j:= l to n do
  121.                       begin s:= 0;
  122.                             for k:= l to n do s:= s + ab[i,k]*ab[k,j];
  123.                             for k:= l to n do ab[k,j]:= ab[k,j] + s*ab[k,i];
  124.                       end; (* j *)
  125.                 end; (* if *)
  126.                 for j:= l to n do
  127.                 begin ab[j,i]:= 0;
  128.                       ab[i,j]:= 0;
  129.                 end;
  130.                 ab[i,i]:= 1; g:= e[i]; l:= i;
  131.           end; (* i *)
  132.           (* diagonalization to bidiagonal form *)
  133.           eps:= eps*x;
  134.           for k:= n downto 1 do
  135.           begin kt:= 0;
  136. 1: kt:= kt + 1;
  137.                 if kt > 30 then
  138.                 begin e[k]:= 0;
  139.                       writeln('+++ qr failed');
  140.                 end;
  141.                 for l2:= k downto 1 do
  142.                 begin l:= l2;
  143.                       if abs(e[l])<=eps then goto 2;
  144.                       if abs(q[l-1])<=eps then goto 4;
  145.                 end; (* l2 *)
  146. 4:   c:= 0; s:= 1;
  147.                 for i:= l to k do
  148.                 begin f:= s*e[i]; e[i]:= c*e[i];
  149.                       if abs(f)<=eps then goto 2;
  150.                       g:= q[i];
  151.                       if abs(f) < abs(g)
  152.                          then h:= abs(g)*sqrt(1+sqr(f/g))
  153.                          else if f <> 0
  154.                                  then h:= abs(f)*sqrt(1+sqr(g/f))
  155.                                  else h:= 0;
  156.                       q[i]:= h;
  157.                       if h = 0 then
  158.                       begin h:= 1;
  159.                             g:= 1;
  160.                       end;
  161.                       c:= g/h; s:= -f/h;
  162.                 end; (* i *)
  163. 2:
  164.                 z:= q[k];
  165.                 if l=k then goto 3;
  166.                 (* shift from bottom 2*2 minor *)
  167.                 x:= q[l]; y:= q[k-1]; g:= e[k-1];  h:= e[k];
  168.                 f:= ((y-z)*(y+z) + (g-h)*(g+h))/(2*h*y);
  169.                 g:= sqrt(sqr(f)+1);
  170.                 if f<=0
  171.                    then f:= ((x-z)*(x+z)+h*(y/(f-g)-h))/x
  172.                    else f:= ((x-z)*(x+z)+h*(y/(f+g)-h))/x;
  173.                 (* next qr transformation *)
  174.                 s:= 1; c:= 1;
  175.                 for i:= l+1 to k do
  176.                 begin g:= e[i]; y:= q[i]; h:= s*g; g:= g*c;
  177.                       if abs(f)<abs(h)
  178.                          then z:= abs(h)*sqrt(1+sqr(f/h))
  179.                          else if f<>0
  180.                                  then z:= abs(f)*sqrt(1+sqr(h/f))
  181.                                  else z:= 0;
  182.                       e[i-1]:= z;
  183.                       if z=0 then
  184.                       begin f:= 1;
  185.                             z:= 1;
  186.                       end;
  187.                       c:= f/z;  s:= h/z;
  188.                       f:= x*c + g*s;  g:= - x*s + g*c; h:= y*s;
  189.                       y:= y*c;
  190.                       for j:= 1 to n do
  191.                       begin x:= ab[j,i-1]; z:= ab[j,i];
  192.                             ab[j,i-1]:= x*c + z*s;
  193.                             ab[j,i]:= - x*s + z*c;
  194.                       end;
  195.                       if abs(f)<abs(h)
  196.                          then z:= abs(h)*sqrt(1+sqr(f/h))
  197.                          else if f<>0
  198.                                  then z:= abs(f)*sqrt(1+sqr(h/f))
  199.                                  else z:= 0;
  200.                       q[i-1]:= z;
  201.                       if z=0 then
  202.                       begin z:= 1;
  203.                             f:= 1;
  204.                       end;
  205.                       c:= f/z; s:= h/z;
  206.                       f:= c*g + s*y; x:= -s*g + c*y;
  207.                 end; (* i *)
  208.                 e[l]:= 0; e[k]:= f; q[k]:= x;
  209.                 goto 1;
  210. 3:    if z<0 then
  211.                 begin q[k]:= -z;
  212.                       for j:= 1 to n do ab[j,k]:= - ab[j,k];
  213.                 end;
  214.           end; (* k *)
  215. end;      (* minfit *)
  216.  
  217. function praxis(var x:vector; n:integer):real;
  218. label    5, 6, 7;
  219. var      i, j,
  220.          k, k2,
  221.          nl, nf, kl, kt: integer;
  222.          s, sl, dn, dmin,
  223.          fx, f1,
  224.          lds, ldt, sf, df,
  225.          qf1, qd0, qd1,
  226.          qa, qb, qc,
  227.          m2, m4,
  228.          small, vsmall,
  229.          large, vlarge,
  230.          ldfac, t2:      real;
  231.          d, y, z,
  232.          q0, q1:         vector;
  233.          v:              matrix;
  234.  
  235. procedure sort;                         (* sort d and v in descending order *)
  236. var       k, i, j:       integer;
  237.           s:             real;
  238. begin     for i:= 1 to n-1 do
  239.           begin k:= i; s:= d[i];
  240.                 for j:= i+1 to n do
  241.                     if d[j] > s then
  242.                     begin k:= j; s:= d[j];
  243.                     end;
  244.                 if k>i then
  245.                 begin d[k]:= d[i]; d[i]:= s;
  246.                       for j:= 1 to n do
  247.                       begin s:= v[j,i];
  248.                             v[j,i]:= v[j,k];
  249.                             v[j,k]:= s;
  250.                       end;
  251.                 end; (* if *)
  252.           end; (* for *)
  253. end;      (* sort *)
  254.  
  255. procedure print;                                   (* print a line of traces *)
  256. var       i:             integer;
  257. begin     writeln('------------------------------------------------------');
  258.           writeln('chi square reduced to ... ', fx);
  259.           writeln(' ... after ',nf,' function calls ...');
  260.           writeln(' ... including ',nl,' linear searches.');
  261.           writeln('current values of x ...');
  262.           for i:= 1 to n do
  263.               write(x[i]);
  264.           writeln;
  265. end;      (* print *)
  266.  
  267. procedure matprint(s:prstring;v:matrix;n,m:integer);
  268. var       k, i:          integer;
  269. begin     writeln;
  270.           writeln(s);
  271.           for k:= 1 to n do
  272.           begin for i:= 1 to m do
  273.                     write(v[k,i]);
  274.                 writeln;
  275.           end;
  276.           writeln;
  277. end;      (* matprint *)
  278.  
  279. procedure vecprint(s:prstring;v:vector;n:integer);
  280. var       i:             integer;
  281. begin     writeln;
  282.           writeln(s);
  283.           for i:= 1 to n do
  284.               write(v[i]);
  285.           writeln;
  286. end;      (* vecprint *)
  287.  
  288. procedure min(j, nits:integer; var d2, x1:real; f1:real;fk:boolean);
  289. label     5, 6;
  290. var       k, i:          integer;
  291.           dz:            boolean;
  292.           x2, xm, f0,
  293.           f2, fm, d1, t2,
  294.           s, sf1, sx1:   real;
  295.  function flin(l:real):real;
  296.  var      i:             integer;
  297.           t:             vector;
  298.  begin    if j>0 then                     (* linear search *)
  299.              for i:= 1 to n do
  300.                  t[i]:= x[i]+l*v[i,j]
  301.              else begin              (* search along parabolic space curve *)
  302.                   qa:= l*(l-qd1)/(qd0*(qd0+qd1));
  303.                   qb:= (l+qd0)*(qd1-l)/(qd0*qd1);
  304.                   qc:= l*(l+qd0)/(qd1*(qd0+qd1));
  305.                   for i:= 1 to n do
  306.                       t[i]:= qa*q0[i]+qb*x[i]+qc*q1[i];
  307.              end; (* else *)
  308.           nf:= nf+1;
  309.           flin:= f(t, n);
  310.  end;     (* flin *)
  311. begin     (* min *)
  312.           sf1:= f1; sx1:= x1;
  313.           k:= 0; xm:= 0; fm:= fx; f0:= fx; dz:= (d2<macheps);
  314.           (* find step size *)
  315.           s:= 0; for i:= 1 to n do s:= s + sqr(x[i]);
  316.           s:= sqrt(s);
  317.           if dz
  318.              then t2:= m4*sqrt(abs(fx)/dmin + s*ldt) + m2*ldt
  319.              else t2:= m4*sqrt(abs(fx)/d2 + s*ldt) + m2*ldt;
  320.           s:= m4*s + t;
  321.           if dz and (t2>s) then t2:= s;
  322.           if (t2<small) then t2:= small;
  323.           if (t2>(0.01*h)) then t2:= 0.01*h;
  324.           if fk and (f1<=fm) then begin xm:= x1; fm:= f1; end;
  325.           if not fk or (abs(x1)<t2) then
  326.           begin if x1>0 then x1:= t2 else x1:= - t2;
  327.                 f1:= flin(x1);
  328.           end;
  329.           if f1<=fm then begin xm:= x1; fm:= f1; end;
  330. 5:       if dz then
  331.           begin if f0<f1 then x2:= - x1 else x2:= 2*x1;
  332.                 f2:= flin(x2);
  333.                 if f2 <= fm then begin xm:= x2; fm:= f2; end;
  334.                 d2:= (x2*(f1-f0) - x1*(f2-f0))/(x1*x2*(x1-x2));
  335.           end;
  336.           d1:= (f1-f0)/x1 - x1*d2; dz:= true;
  337.           if d2 <= small
  338.              then if d1<0
  339.                      then x2:= h
  340.                      else x2:= -h
  341.              else x2:= -0.5*d1/d2;
  342.           if abs(x2)>h
  343.              then if x2>0
  344.                      then x2:= h
  345.                      else x2:= -h;
  346. 6:       f2:= flin(x2);
  347.           if (k<nits) and (f2>f0) then
  348.           begin k:= k + 1;
  349.                 if (f0<f1) and ((x1*x2)>0) then goto 5;
  350.                 x2:= 0.5*x2; goto 6;
  351.           end;
  352.           nl:= nl + 1;
  353.           if f2>fm then x2:= xm else fm:= f2;
  354.           if abs(x2*(x2-x1))>small
  355.              then d2:= (x2*(f1-f0) - x1*(fm-f0))/(x1*x2*(x1-x2))
  356.              else if k>0
  357.                      then d2:= 0;
  358.           if d2<=small then d2:= small;
  359.           x1:= x2; fx:= fm;
  360.           if sf1<fx then begin fx:= sf1; x1:= sx1; end;
  361.           if j>0
  362.              then for i:= 1 to n do
  363.                       x[i]:= x[i] + x1*v[i,j];
  364. end;      (* min *)
  365.  
  366. procedure quad;        (* look for a minimum along the curve q0, q1, q2 *)
  367. var       i:             integer;
  368.           l, s:          real;
  369. begin     s:= fx; fx:= qf1; qf1:= s; qd1:= 0;
  370.           for i:= 1 to n do
  371.           begin s:= x[i];  l:= q1[i];  x[i]:= l;  q1[i]:= s;
  372.                 qd1:= qd1 + sqr(s - l);
  373.           end;
  374.           s:= 0; qd1:= sqrt(qd1); l:= qd1;
  375.           if (qd0>0) and (qd1>0) and (nl >= (3*n*n)) then
  376.           begin min(0, 2, s, l, qf1, true);
  377.                 qa:= l*(l-qd1)/(qd0*(qd0+qd1));
  378.                 qb:= (l+qd0)*(qd1-l)/(qd0*qd1);
  379.                 qc:= l*(l+qd0)/(qd1*(qd0+qd1));
  380.           end else
  381.           begin fx:= qf1; qa:= 0; qb:= 0; qc:= 1;
  382.           end;
  383.           qd0:= qd1;
  384.           for i:= 1 to n do
  385.           begin s:= q0[i];  q0[i]:= x[i];
  386.                 x[i]:= qa*s + qb*x[i] + qc*q1[i];
  387.           end;
  388. end;      (* quad *)
  389.  
  390. begin     (****  p r a x i s  ****)
  391.           (* initialization *)
  392.           small:= sqr(macheps); vsmall:= sqr(small);
  393.           large:= 1.0/small;    vlarge:= 1.0/vsmall;
  394.           m2:= sqrt(macheps);   m4:= sqrt(m2);
  395.           if illc then ldfac:= 0.1 else ldfac:= 0.01;
  396.           nl:= 0; kt:= 0; nf:= 1; fx:= f(x, n); qf1:= fx;
  397.           t2:= small + abs(t); t:= t2; dmin:= small;
  398.           if h<(100*t) then h:= 100*t; ldt:= h;
  399.           for i:= 1 to n do for j:= 1 to n do
  400.               if i=j then v[i,j]:= 1 else v[i,j]:= 0;
  401.           d[1]:= 0; qd0:= 0;
  402.           for i:= 1 to n do q1[i]:= x[i];
  403.           if prin > 1 then
  404.           begin writeln;writeln('---------- enter function praxis ------------');
  405.                 writeln('current paramter settings are:');
  406.                 writeln('... scaling ... ',scbd);
  407.                 writeln('... macheps ... ',macheps);
  408.                 writeln('...   tol   ... ',t);
  409.                 writeln('... maxstep ... ',h);
  410.                 writeln('...   illc  ... ',illc);
  411.                 writeln('...   ktm   ... ',ktm);
  412.           end;
  413.           if prin > 0 then print;
  414.  
  415. 5:       (* main loop *)
  416.           sf:= d[1]; s:= 0; d[1]:= 0;
  417.           (* minimize along first direction *)
  418.           min(1, 2, d[1], s, fx, false);
  419.           if s<= 0
  420.              then for i:= 1 to n do
  421.                       v[i,1]:= - v[i,1];
  422.           if (sf<= (0.9*d[1])) or ((0.9*sf)>=d[1])
  423.              then for i:= 2 to n do
  424.                       d[i]:= 0;
  425.           for k:= 2 to n do
  426.           begin for i:= 1 to n do y[i]:= x[i];
  427.                 sf:= fx;
  428.                 illc:= illc or (kt>0);
  429. 6:             kl:= k; df:= 0;
  430.                 if illc then   (* random step to get off resolution valley *)
  431.                 begin for i:= 1 to n do
  432.                       begin z[i]:= (0.1*ldt + t2*power(10,kt))*(random(1)-0.5);
  433.                             s:= z[i];
  434.                             for j:= 1 to n do x[j]:= x[j]+s*v[j,i];
  435.                       end; (* i *)
  436.                       fx:= f(x, n);   nf:= nf + 1;
  437.                 end; (* if *)
  438.                 for k2:= k to n do  (* minimize along non-conjugate directions *)
  439.                 begin sl:= fx; s:= 0;
  440.                       min(k2, 2, d[k2], s, fx, false);
  441.                       if illc
  442.                          then s:= d[k2]*sqr(s+z[k2])
  443.                          else s:= sl - fx;
  444.                       if df<s then begin df:= s;  kl:= k2; end;
  445.                 end; (* k2 *)
  446.                 if not illc and (df < abs(100*macheps*fx)) then
  447.                 begin illc:= true; goto 6;
  448.                 end;
  449.                 if (k=2) and (prin>1) then vecprint('new direction ...', d, n);
  450.                 for k2:= 1 to k-1 do (* minimize along conjugate directions *)
  451.                 begin s:= 0;
  452.                       min(k2, 2, d[k2], s, fx, false);
  453.                 end; (* k2 *)
  454.                 f1:= fx; fx:= sf; lds:= 0;
  455.                 for i:= 1 to n do
  456.                 begin sl:= x[i]; x[i]:= y[i]; y[i]:= sl - y[i]; sl:= y[i];
  457.                       lds:= lds + sqr(sl);
  458.                 end;
  459.                 lds:= sqrt(lds);
  460.                 if lds>small then
  461.                 begin for i:= kl-1 downto k do
  462.                       begin for j:= 1 to n do v[j,i+1]:= v[j,i];
  463.                             d[i+1]:= d[i];
  464.                       end;
  465.                       d[k]:= 0;
  466.                       for i:= 1 to n do v[i,k]:= y[i]/lds;
  467.                       min(k, 4, d[k], lds, f1, true);
  468.                       if lds<=0 then
  469.                       begin lds:= -lds;
  470.                             for i:= 1 to n do v[i,k]:= -v[i,k];
  471.                       end;
  472.                 end; (* if *)
  473.                 ldt:= ldfac*ldt; if ldt<lds then ldt:= lds;
  474.                 if prin > 1 then print;
  475.                 t2:= 0; for i:= 1 to n do t2:= t2 + sqr(x[i]);
  476.                 t2:= m2*sqrt(t2) + t;
  477.                 if ldt > (0.5*t2) then kt:= 0 else kt:= kt+1;
  478.                 if kt>ktm then goto 7;
  479.           end; (* k *)
  480.           (*  try quadratic extrapolation in case    *)
  481.           (*  we are stuck in a curved valley        *)
  482.           quad;
  483.           dn:= 0;
  484.           for i:= 1 to n do
  485.           begin d[i]:= 1.0/sqrt(d[i]);
  486.                 if dn<d[i] then dn:= d[i];
  487.           end;
  488.           if prin>2 then matprint('new matrix of directions ...', v, n, n);
  489.           for j:= 1 to n do
  490.           begin s:= d[j]/dn;
  491.                 for i:= 1 to n do v[i,j]:= s*v[i,j];
  492.           end;
  493.           if scbd > 1 then  (* scale axis to reduce condition number *)
  494.           begin s:= vlarge;
  495.                 for i:= 1 to n do
  496.                 begin sl:= 0;
  497.                       for j:= 1 to n do sl:= sl + sqr(v[i,j]);
  498.                       z[i]:= sqrt(sl);
  499.                       if z[i]<m4 then z[i]:= m4;
  500.                       if s>z[i] then s:= z[i];
  501.                 end; (* i *)
  502.                 for i:= 1 to n do
  503.                 begin sl:= s/z[i];  z[i]:= 1.0/sl;
  504.                       if z[i] > scbd then
  505.                       begin sl:= 1/scbd;
  506.                             z[i]:= scbd;
  507.                       end;
  508.                 end;
  509.           end; (* if *)
  510.           for i:= 2 to n do  for j:= 1 to i-1 do
  511.           begin s:= v[i,j]; v[i,j]:= v[j,i]; v[j,i]:= s;
  512.           end;
  513.           minfit(n, macheps, vsmall, v, d);
  514.           if scbd>1 then
  515.           begin for i:= 1 to n do
  516.                 begin s:= z[i];
  517.                       for j:= 1 to n do v[i,j]:= v[i,j]*s;
  518.                 end;
  519.                 for i:= 1 to n do
  520.                 begin s:= 0;
  521.                       for j:= 1 to n do s:= s + sqr(v[j,i]);
  522.                       s:= sqrt(s);  d[i]:= s*d[i];  s:= 1.0/s;
  523.                       for j:= 1 to n do v[j,i]:= s*v[j,i];
  524.                 end;
  525.           end;
  526.           for i:= 1 to n do
  527.           begin if (dn*d[i])>large
  528.                    then d[i]:= vsmall
  529.                    else if (dn*d[i])<small
  530.                            then d[i]:= vlarge
  531.                            else d[i]:= power(dn*d[i],-2);
  532.           end;
  533.           sort;   (* the new eigenvalues and eigenvectors *)
  534.           dmin:= d[n];
  535.           if dmin < small then dmin:= small;
  536.           illc:= (m2*d[1]) > dmin;
  537.           if (prin > 2) and (scbd > 1)
  538.              then vecprint('scale factors ...', z, n);
  539.           if prin > 2 then vecprint('eigenvalues of a ...', d, n);
  540.           if prin > 2 then matprint('eigenvectors of a ...', v, n, n);
  541.           goto 5;  (* back to main loop *)
  542. 7:        if prin > 0 then
  543.           begin vecprint('final solution is ...', x, n);
  544.                 writeln; writeln('chisq reduced to ',fx,' after ',nf, ' function calls.');
  545.           end;
  546.           praxis:= fx;
  547. end;      (****  praxis  ****)
  548.