PROGRAM d9r13(input,output);
(* driver for routine MNEWT *)

(*$I MODFILE.PAS *)

CONST
   ntrial = 5;
   tolx = 1.0e-6;
   n = 4;
   np = n;
   tolf = 1.0e-6;
TYPE
   RealArrayNP = ARRAY [1..n] OF real;
   RealArrayNPbyNP = ARRAY [1..n,1..n] OF real;
   IntegerArrayNP = ARRAY [1..n] OF integer;
VAR
   i,j,k,kk,l: integer;
   xx: real;
   x,beta: RealArrayNP;
   alpha: RealArrayNPbyNP;

PROCEDURE usrfun(VAR x: RealArrayNP; n: integer; VAR alpha: RealArrayNPbyNP;
       VAR beta: RealArrayNP);
(* Programs using routine USRFUN must define the types
TYPE
   RealArrayNP = ARRAY [1..n] OF real;
   RealArrayNPbyNP = ARRAY [1..n,1..n] OF real;
in the main routine. *)
BEGIN
   alpha[1,1] := 2.0*x[1];
   alpha[1,2] := 2.0*x[2];
   alpha[1,3] := 2.0*x[3];
   alpha[1,4] := -3.0;
   alpha[2,1] := 2.0*x[1];
   alpha[2,2] := 2.0*x[2];
   alpha[2,3] := 2.0*x[3];
   alpha[2,4] := -2.0*x[4];
   alpha[3,1] := 1.0;
   alpha[3,2] := -1.0;
   alpha[3,3] := 0.0;
   alpha[3,4] := 0.0;
   alpha[4,1] := 0.0;
   alpha[4,2] := 1.0;
   alpha[4,3] := -1.0;
   alpha[4,4] := 0.0;
   beta[1] := -sqr(x[1])-sqr(x[2])-sqr(x[3])+3.0*x[4];
   beta[2] := -sqr(x[1])-sqr(x[2])-sqr(x[3])+sqr(x[4])+1.0;
   beta[3] := -x[1]+x[2];
   beta[4] := -x[2]+x[3]
END;

(*$I LUBKSB.PAS *)

(*$I LUDCMP.PAS *)

(*$I MNEWT.PAS *)

BEGIN
   FOR l := 0 TO 1 DO BEGIN
      kk := 2*l-1;
      FOR k := 1 TO 3 DO BEGIN
         xx := 0.2*k*kk;
         writeln('Starting vector number',k:2);
         FOR i := 1 TO 4 DO
            x[i] := xx+0.2*i;
         FOR i := 1 TO 4 DO
            writeln('x[':7,i:1,']  :=  ',x[i]:5:2);
         writeln;
         FOR j := 1 TO ntrial DO BEGIN
            mnewt(1,x,n,tolx,tolf);
            usrfun(x,n,alpha,beta);
            writeln('i':5,'x[i]':13,'f':13);
            FOR i := 1 TO n DO 
               writeln(i:5,'   ',x[i]:12,'   ',-beta[i]:12);
            writeln;
            writeln('press RETURN to continue...');
            readln
         END
      END
   END
END.