PROGRAM d9r1(input,output);
(* driver for routine SCRSHO *)

(*$I MODFILE.PAS *)
(*$I BESSJ0.PAS *)

FUNCTION fx(x: real): real;
BEGIN
   fx := bessj0(x)
END;

(*$I SCRSHO.PAS *)

BEGIN
   scrsho
END.
PROGRAM d8r9(input,output,dfile);
(* driver for routine SHELL *)
CONST
   npt=100;
TYPE
   glnarray = ARRAY [1..npt] OF real;
VAR
   i,j : integer;
   a : glnarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SHELL.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to npt DO read(dfile,a[i]);
   close(dfile);
(* write original array *)
   writeln('Original array:');
   FOR i := 1 to (npt DIV 10) DO BEGIN
      FOR j := 1 to 10 DO write(a[10*(i-1)+j]:6:2);
      writeln
   END;
(* write sorted array *)
   shell(npt,a);
   writeln;
   writeln('Sorted array:');
   FOR i := 1 to (npt DIV 10) DO BEGIN
      FOR j := 1 to 10 DO write(a[10*(i-1)+j]:6:2);
      writeln
   END;
END.
PROGRAM d16r1(input,output);
(* driver for routine SHOOT *)
(* Solves for eigenvalues of spheroidal harmonics. Both
prolate and oblate case are handled simultaneously, leading
to six first-order equations. Unknown to shoot, it is
actually two independent sets of three coupled equations,
one set with c^2 positive and the other with c^2 negative.  *)
LABEL 1,2;
CONST
   nvar=6;
   n2=2;
   np=2;
   delta=1.0e-3;
   eps=1.0e-6;
   dx=1.0e-4;
TYPE
   gl2array = ARRAY [1..2] OF real;
   gln2array = gl2array;
   gl6array = ARRAY [1..6] OF real;
   glnarray = gl6array;
   glnvar = gl6array;
   glarray = gl6array;
   glindx = ARRAY [1..6] OF integer;
   glinvar = glindx;
   gln2byn2 = ARRAY [1..n2,1..n2] OF real;
   glnpbynp = gln2byn2;
VAR
   c2,factr,h1,hmin,q1,x1,x2 : real;
   i,m,n : integer;
   delv,v : gl2array;
   dv,f : glnarray;
   kmax,kount : integer;
   dxsav : real;
   xp : ARRAY [1..200] OF real;
   yp : ARRAY [1..10,1..200] OF real;

PROCEDURE load(x1: real; v: gl2array; VAR y: gl6array);
(* Programs using routine LOAD must define the types
TYPE
   gl2array = ARRAY [1..2] OF real;
   gl6array = ARRAY [1..6] OF real;
and must define the global variables
VAR
   c2,factr : real;
   m : integer;
in the main routine. *)
BEGIN
   y[3] := v[1];
   y[2] := -(y[3]-c2)*factr/2.0/(m+1.0);
   y[1] := factr+y[2]*dx;
   y[6] := v[2];
   y[5] := -(y[6]+c2)*factr/2.0/(m+1.0);
   y[4] := factr+y[5]*dx
END;

PROCEDURE score(x2: real; y: gl6array; VAR f: gl6array);
(* Programs using routine SCORE must define the types
TYPE
   gl6array = ARRAY [1..6] OF real;
and must define the global variables
VAR
   m,n : integer;
in the main routine. *)
BEGIN
   IF (((n-m) MOD 2) = 0) THEN BEGIN
      f[1] := y[2];
      f[2] := y[5]
   END ELSE BEGIN
      f[1] := y[1];
      f[2] := y[4]
   END
END;

PROCEDURE derivs(x: real; y: gl6array; VAR dydx: gl6array);
(* Programs using routine DERIVS must define the type
TYPE
   gl6array = ARRAY [1..6] OF real;
and must define the global variables
VAR
   c2 : real;
   m : integer;
in the main routine. *)
BEGIN
   dydx[1] := y[2];
   dydx[3] := 0.0;
   dydx[2] := (2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x);
   dydx[4] := y[5];
   dydx[6] := 0.0;
   dydx[5] := (2.0*x*(m+1.0)*y[5]-(y[6]+c2*x*x)*y[4])/(1.0-x*x)
END;

(*$I MODFILE.PAS *)
(*$I LUBKSB.PAS *)

(*$I LUDCMP.PAS *)

(*$I RK4.PAS *)

(*$I RKQC.PAS *)

(*$I ODEINT.PAS *)

(*$I SHOOT.PAS *)

BEGIN
1:   write('Input M,N,C-Squared:  ');
   readln(m,n,c2);
   IF ((n < m) OR (m < 0)) THEN GOTO 1;
   factr := 1.0;
   IF  (m <> 0)  THEN BEGIN
      q1 := n;
      FOR i := 1 to m DO BEGIN
         factr := -0.5*factr*(n+i)*(q1/i);
         q1 := q1-1.0
      END
   END;
   v[1] := n*(n+1)-m*(m+1)+c2/2.0;
   v[2] := n*(n+1)-m*(m+1)-c2/2.0;
   delv[1] := delta*v[1];
   delv[2] := delv[1];
   h1 := 0.1;
   hmin := 0.0;
   x1 := -1.0+dx;
   x2 := 0.0;
   writeln;
   writeln('Prolate':17,'Oblate':23);
   writeln('Mu(m,n)':11,'Error Est.':14,'Mu(m,n)':10,'Error Est.':14);
2:   shoot(nvar,v,delv,n2,x1,x2,eps,h1,hmin,f,dv);
   writeln(v[1]:12:6,dv[1]:12:6,v[2]:12:6,dv[2]:12:6);
   IF ((abs(dv[1]) > abs(eps*v[1])) OR 
      ((dv[2]) > abs(eps*v[2]))) THEN GOTO 2
END.
PROGRAM d16r2(input,output);
(* driver for routine SHOOTF *)
LABEL 1;
CONST
   nvar=3;
   n1=2;
   n2=1;
   delta=1.0e-3;
   eps=1.0e-6;
   dx=1.0e-4;
TYPE
   gl1array = ARRAY [1..1] OF real;
   gln2array = gl1array;
   gl2array = ARRAY [1..2] OF real;
   gln1array = gl2array;
   gl3array = ARRAY [1..3] OF real;
   glnarray = gl3array;
   glarray = gl3array;
   glnvar = gl3array;
   glnpbynp = ARRAY [1..3,1..3] OF real;
   glnvarbynvar = glnpbynp;
   glindx = ARRAY [1..3] OF integer;
VAR
   c2,factr,h1,hmin : real;
   q1,x1,x2,xf : real;
   i,m,n : integer;
   v1,delv1,dv1 : gln2array;
   v2,delv2,dv2 : gln1array;
   f : glnvar;
   kmax,kount : integer;
   dxsav : real;
   xp : ARRAY [1..200] OF real;
   yp : ARRAY [1..10,1..200] OF real;

PROCEDURE load1(x1: real; v1: gl1array; VAR y: gl3array);
(* Programs using routine LOAD1 must define the types
TYPE
   gl1array = ARRAY [1..1] OF real;
   gl3array = ARRAY [1..3] OF real;
and must declare the variables
VAR
   c2,factr : real;
   m : integer;
in the main routine. *)
BEGIN
   y[3] := v1[1];
   y[2] := -(y[3]-c2)*factr/2.0/(m+1.0);
   y[1] := factr+y[2]*dx
END;

PROCEDURE load2(x2: real; v2: gl2array; VAR y: gl3array);
(* Programs using routine LOAD2 must define the types
TYPE
   gl2array = ARRAY [1..2] OF real;
   gl3array = ARRAY [1..3] OF real;
and must declare the variables
   c2 : real;
   m : integer;
in the main routine. *)
BEGIN
   y[3] := v2[2];
   y[1] := v2[1];
   y[2] := (y[3]-c2)*y[1]/2.0/(m+1.0)
END;

PROCEDURE score(xf: real; y: gl3array; VAR f: gl3array);
(* Programs using routine SCORE must define the type
TYPE
   gl3array = ARRAY [1..3] OF real;
in the main routine. *)
VAR
   i : integer;
BEGIN
   FOR i := 1 to 3 DO BEGIN
      f[i] := y[i]
   END
END;

PROCEDURE derivs(x: real; y: gl3array; VAR dydx: gl3array);
(* Programs using routine DERIVS must define the type
TYPE
   gl3array = ARRAY [1..3] OF real;
and must declare the variables
   c2 : real;
   m : integer;
in the main routine. *)
BEGIN
   dydx[1] := y[2];
   dydx[3] := 0.0;
   dydx[2] := (2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x)
END;

(*$I MODFILE.PAS *)
(*$I LUBKSB.PAS *)

(*$I LUDCMP.PAS *)

(*$I RK4.PAS *)

(*$I RKQC.PAS *)

(*$I ODEINT.PAS *)

(*$I SHOOTF.PAS *)

BEGIN
1:   write('Input M,N,C-SQUARED: ');
   readln(m,n,c2);
   IF ((n < m) OR (m < 0)) THEN BEGIN
      writeln('Improper arguments');
      GOTO 1
   END;
   factr := 1.0;
   IF (m <> 0) THEN BEGIN
      q1 := n;
      FOR i := 1 to m DO BEGIN
         factr := -0.5*factr*(n+i)*(q1/i);
         q1 := q1-1.0
      END
   END;
   v1[1] := n*(n+1)-m*(m+1)+c2/2.0;
   IF (((n-m) MOD 2) = 0) THEN v2[1] := factr
   ELSE v2[1] := -factr;
   v2[2] := v1[1];
   delv1[1] := delta*v1[1];
   delv2[1] := delta*factr;
   delv2[2] := delv1[1];
   h1 := 0.1;
   hmin := 0.0;
   x1 := -1.0+dx;
   x2 := 1.0-dx;
   xf := 0.0;
   writeln;
   writeln('mu(-1)':26,'y(1-dx)':20,'mu(+1)':19);
   REPEAT
      shootf(nvar,v1,v2,delv1,delv2,n1,n2,x1,x2,
         xf,eps,h1,hmin,f,dv1,dv2);
      writeln;
      writeln('v ':6,v1[1]:20:6,v2[1]:20:6,v2[2]:20:6);
      writeln('dv':6,dv1[1]:20:6,dv2[1]:20:6,dv2[2]:20:6);
   UNTIL (abs(dv1[1]) <= abs(eps*v1[1]))
END.
PROGRAM d10r12(input,output);
(* driver for routine SIMPLX *)
(* incorporates examples discussed in text *)
CONST
   n=4;
   m=4;
   np=5;   (* np >= n+1 *)
   mp=6;   (* mp >= m+2 *)
   m1=2;   (* m1+m2+m3=m *)
   m2=1;
   m3=1;
   nm1m2=7;   (* nm1m2=n+m1+m2 *)
TYPE
   twochar = PACKED ARRAY [1..2] OF char;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnarray = ARRAY [1..n] OF integer;
   glmarray = ARRAY [1..m] OF integer;
   glnparray = ARRAY [1..np] OF integer;
   glmparray = ARRAY [1..mp] OF integer;
VAR
   i,icase,j : integer;
   izrov : glnarray;
   iposv : glmarray;
   a : glmpbynp;
   txt : ARRAY [1..nm1m2] OF twochar;

(*$I MODFILE.PAS *)
(*$I SIMP1.PAS *)

(*$I SIMP2.PAS *)

(*$I SIMP3.PAS *)

(*$I SIMPLX.PAS *)

BEGIN
   txt[1] := 'x1'; txt[2] := 'x2'; txt[3] := 'x3';
   txt[4] := 'x4'; txt[5] := 'y1'; txt[6] := 'y2';
   txt[7] := 'y3';
   a[1,1] := 0.0; a[1,2] := 1.0; a[1,3] := 1.0;
   a[1,4] := 3.0; a[1,5] := -0.5;
   a[2,1] := 740.0; a[2,2] := -1.0; a[2,3] := 0.0;
   a[2,4] := -2.0; a[2,5] := 0.0;
   a[3,1] := 0.0; a[3,2] := 0.0; a[3,3] := -2.0;
   a[3,4] := 0.0; a[3,5] := 7.0;
   a[4,1] := 0.5; a[4,2] := 0.0; a[4,3] := -1.0;
   a[4,4] := 1.0; a[4,5] := -2.0;
   a[5,1] := 9.0; a[5,2] := -1.0; a[5,3] := -1.0;
   a[5,4] := -1.0; a[5,5] := -1.0;
   simplx(a,m,n,mp,np,m1,m2,m3,icase,izrov,iposv);
   writeln;
   IF (icase = 1) THEN BEGIN
      writeln('unbounded objective function')
   END ELSE IF (icase = -1) THEN BEGIN
      writeln('no solutions satisfy constraints given')
   END ELSE BEGIN
      write(' ':11);
      FOR i := 1 to n DO BEGIN
         IF (izrov[i] <= nm1m2) THEN BEGIN
            write(txt[izrov[i]]:10)
         END
      END;
      writeln;
      FOR i := 1 to m+1 DO BEGIN
         IF (i > 1) THEN BEGIN
            write(txt[iposv[i-1]])
         END ELSE BEGIN
            write('  ')
         END;
         FOR j := 1 to (n+1) DO BEGIN
            IF (j=1) THEN write(a[i,j]:10:2);
            IF (j>1) THEN BEGIN
               IF (izrov[j-1] <= nm1m2) THEN
                  write(a[i,j]:10:2)
            END
         END;
         writeln
      END
   END
END.
PROGRAM d12r4(input,output);
(* driver for routine SINFT *)
LABEL 1,99;
CONST
   eps=1.0e-3;
   np=16;
   np2=18;      (* np2=np+2 *)
   width=30.0;
   pi=3.1415926;
TYPE
   gldarray=ARRAY [1..np2] OF real;
   glyarray=gldarray;
VAR
   big,per,scal,small : real;
   i,j,nlim : integer;
   data,size : gldarray;

(*$I MODFILE.PAS *)
(*$I FOUR1.PAS *)

(*$I REALFT.PAS *)

(*$I SINFT.PAS *)

BEGIN
1:   writeln('period of sinusoid in channels (3-',np:2,')');
   readln(per);
   IF (per <= 0.0) THEN GOTO 99;
   FOR i := 1 to np DO BEGIN
      data[i] := sin(2.0*pi*(i-1)/per)
   END;
   sinft(data,np);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF  (data[i] < small) THEN  small := data[i];
      IF  (data[i] > big) THEN  big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   writeln('press RETURN to continue ...');
   readln;
   sinft(data,np);
   big := -1.0e10;
   small := 1.0e10;
   FOR i := 1 to np DO BEGIN
      IF (data[i] < small) THEN small := data[i];
      IF (data[i] > big) THEN big := data[i]
   END;
   scal := width/(big-small);
   FOR i := 1 to np DO BEGIN
      nlim := round(scal*(data[i]-small)+eps);
      write(i:4,' ');
      FOR j := 1 to nlim+1 DO write('*');
      writeln
   END;
   GOTO 1;
99:
END.
PROGRAM d13r20(input,output);
(* driver for routine SMOOFT *)
CONST
   n=100;
   hash=0.05;
   scale=100.0;
   pts=10.0;
   m=258;  (* 2 + first integral power of 2 that is greater or
         equal to (n+2*pts) *)
TYPE
   glyarray = ARRAY [1..m] OF real;
   gldarray = glyarray;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   bar : real;
   i,ii,idum,j,k,nstp : integer;
   y : glyarray;
   txt : PACKED ARRAY [1..50] OF char;

(*$I MODFILE.PAS *)
(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I FOUR1.PAS *)

(*$I REALFT.PAS *)

(*$I SMOOFT.PAS *)

BEGIN
   gliset := 0;
   idum := -7;
   FOR i := 1 to n DO BEGIN
      y[i] := 3.0*i/n*exp(-3.0*i/n);
      y[i] := y[i]+hash*gasdev(idum)
   END;
   FOR k := 1 to 3 DO BEGIN
      nstp := n DIV 20;
      writeln('data':8,'graph':11);
      FOR i := 1 to 20 DO BEGIN
         ii := nstp*(i-1)+1;
         FOR j := 1 to 50 DO BEGIN
            txt[j] := ' '
         END;
         bar := round(scale*y[ii]);
         FOR j := 1 to 64 DO BEGIN
            IF (j <= bar) THEN txt[j] := '*'
         END;
         writeln (y[ii]:10:6,' ':4,txt)
      END;
      writeln(' press return to smooth ...');
      readln;
      smooft(y,n,pts)
   END
END.
PROGRAM d6r30(input,output,dfile);
(* driver for routine SNCNDN *)
VAR
   cn,dn,em,emmc,sn,uu,val : real;
   i,nval : integer;
   txt : string[26];
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SNCNDN.PAS *)

BEGIN
   glopen(dfile,'fncval.dat');
   REPEAT readln(dfile,txt) UNTIL (txt = 'Jacobian Elliptic Function');
   readln(dfile,nval);
   writeln(txt);
   writeln('mc':4,'u':8,'actual':16,'sn':13,
      'sn^2+cn^2':15,'(mc)*(sn^2)+dn^2':18);
   FOR i := 1 to nval DO BEGIN
      readln(dfile,em,uu,val);
      emmc := 1.0-em;
      sncndn(uu,emmc,sn,cn,dn);
      writeln(emmc:5:2,uu:8:2,val:15:5,sn:15:5,
         (sn*sn+cn*cn):12:5,(em*sn*sn+dn*dn):14:5)
   END;
   close(dfile)
END.
PROGRAM d17r1(input,output);
(* driver for routine SOR *)
CONST
   jmax=11;
   pi=3.1415926;
TYPE
   double = real;
   gljmax = ARRAY [1..jmax,1..jmax] OF double;
VAR
   i,j,midl : integer;
   rjac : double;
   a,b,c,d,e,f,u : gljmax;

(*$I SOR.PAS *)

BEGIN
   FOR i := 1 to jmax DO BEGIN
      FOR j := 1 to jmax DO BEGIN
         a[i,j] := 1.0;
         b[i,j] := 1.0;
         c[i,j] := 1.0;
         d[i,j] := 1.0;
         e[i,j] := -4.0;
         f[i,j] := 0.0;
         u[i,j] := 0.0
      END
   END;
   midl := (jmax DIV 2)+1;
   f[midl,midl] := 2.0;
   rjac := cos(pi/jmax);
   sor(a,b,c,d,e,f,u,jmax,rjac);
   writeln('SOR Solution:');
   FOR i := 1 to jmax DO BEGIN
      FOR j := 1 to jmax DO write(u[i,j]:7:2);
      writeln
   END;
   writeln;
   writeln('Test that solution satisfies difference eqns:');
   FOR i := 2 to (jmax-1) DO BEGIN
      FOR j := 2 to (jmax-1) DO BEGIN
         f[i,j] := u[i+1,j]+u[i-1,j]+u[i,j+1]
            +u[i,j-1]-4.0*u[i,j]
      END;
      write(' ':7);
      FOR j := 2 to (jmax-1) DO write(f[i,j]:7:2);
      writeln
   END
END.
PROGRAM d8r3(input,output,dfile);
(* driver for routine SORT *)
CONST
   np=100;
TYPE
   glsarray= ARRAY [1..np] OF real;
VAR
   i,j : integer;
   a : glsarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SORT.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to 100 DO BEGIN
      read(dfile,a[i]);
   END;
   close(dfile);
(* writeln original array *)
   writeln ('original array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
(* writeln sorted array *)
   sort(np,a);
   writeln ('sorted array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END
END.

PROGRAM d8r4(input,output,dfile);
(* driver for routine SORT2 *)
CONST
   np=100;
TYPE
   glsarray=ARRAY [1..np] OF real;
VAR
   i,j : integer;
   a,b : glsarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SORT2.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to 100 DO BEGIN
      read(dfile,a[i])
   END;
   close(dfile);
(* generate b-array *)
   FOR i := 1 to 100 DO BEGIN
      b[i] := i-1
   END;
(* sort a and mix b *)
   sort2(100,a,b);
   writeln('after sorting a and mixing b, array a is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('... and array b is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(b[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('press return to continue...');
   readln;
(* sort b and mix a *)
   sort2(100,b,a);
   writeln('after sorting b and mixing a, array a is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln ('... and array b is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(b[10*(i-1)+j]:6:2)
      END;
      writeln
   END
END.

PROGRAM d8r6(input,output,dfile);
(* driver for routine SORT3 *)
CONST
   nlen=64;
TYPE
   glsarray = ARRAY [1..nlen] OF real;
   gliarray = ARRAY [1..nlen] OF integer;
VAR
   i,j : integer;
   a,b,c,wksp : glsarray;
   indx : gliarray;
   amsg1 : PACKED ARRAY [1..40] OF char;
   amsg2 : PACKED ARRAY [1..24] OF char;
   n1,n2 : integer;
   amsg,bmsg,cmsg : PACKED ARRAY [1..nlen] OF char;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I INDEXX.PAS *)

(*$I SORT3.PAS *)

BEGIN
   amsg1 := 'i''d rather have a bottle in front of me ';
   n1 := 40;
   amsg2 := 'than a frontal lobotomy.';
   n2 := 24;
   FOR i := 1 to n1 DO amsg[i] := amsg1[i];
   FOR i := 1 to n2 DO amsg[n1+i] := amsg2[i];
   writeln;
   writeln ('original message:');
   writeln (amsg);
(* read array of random numbers *)
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to nlen DO read(dfile,a[i]);
   close(dfile);
(* create array b and array c *)
   FOR i := 1 to nlen DO BEGIN
      b[i] := i;
      c[i] := nlen+1-i
   END;
(* sort array a while mixing ib and ic *)
   sort3(nlen,a,b,c,wksp,indx);
(* scramble message according to array b *)
   FOR i := 1 to nlen DO BEGIN
      j := round(b[i]);
      bmsg[i] := amsg[j]
   END;
   writeln;
   writeln ('scrambled message:');
   writeln (bmsg);
(* unscramble according to array c *)
   FOR i := 1 to nlen DO BEGIN
      j := round(c[i]);
      cmsg[j] := bmsg[i]
   END;
   writeln;
   writeln ('mirrored message:');
   writeln (cmsg)
END.

PROGRAM d8r1(input,output,dfile);
(* driver for routine PIKSRT *)
CONST
   np=100;
TYPE
   glsarray = ARRAY [1..np] OF real;
VAR
   i,j : integer;
   a : glsarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I PIKSRT.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to 100 DO BEGIN
      read(dfile,a[i])
   END;
   close(dfile);
(* write original array *)
   writeln ('original array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
(* write sorted array *)
   piksrt(np,a);
   writeln('sorted array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END
END.

PROGRAM d8r2(input,output,dfile);
(* driver for routine PIKSR2 *)
CONST
   np=100;
TYPE
   glsarray = ARRAY [1..np] OF real;
VAR
   i,j : integer;
   a,b : glsarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I PIKSR2.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to np DO BEGIN
      read(dfile,a[i])
   END;
   close(dfile);
(* generate b-array *)
   FOR i := 1 to np DO BEGIN
      b[i] := i-1
   END;
(* sort a and mix b *)
   piksr2(np,a,b);
   writeln('after sorting a and mixing b, array a is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('... and array b is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(b[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('press return to continue ...');
   readln;
(* sort b and mix a *)
   piksr2(np,b,a);
   writeln('after sorting b and mixing a, array a is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
   writeln('... and array b is:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(b[10*(i-1)+j]:6:2)
      END;
      writeln
   END
END.

PROGRAM d8r8(input,output,dfile);
(* driver for routine QCKSRT *)
CONST
   np = 100;
TYPE
   glarray = ARRAY [1..np] OF real;
VAR
   i,j : integer;
   a : glarray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I QCKSRT.PAS *)

BEGIN
   glopen(dfile,'tarray.dat');
   readln(dfile);
   FOR i := 1 to 100 DO BEGIN
      read(dfile,a[i])
   END;
   close(dfile);
(* print original array *)
   writeln('original array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END;
(* writeln sorted array *)
   qcksrt(np,a);
   writeln;
   writeln ('sorted array:');
   FOR i := 1 to 10 DO BEGIN
      FOR j := 1 to 10 DO BEGIN
         write(a[10*(i-1)+j]:6:2)
      END;
      writeln
   END
END.
PROGRAM d2r10(input,output);
(* driver for routine SPARSE *)
CONST
   n=20;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   i,ii : integer;
   rsq : real;
   b,bcmp,x : glnarray;

PROCEDURE asub(xin: glnarray; VAR xout: glnarray; n: integer);
(* Programs using routine ASUB must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   i: integer;
BEGIN
   xout[1] := xin[1]+2.0*xin[2];
   xout[n] := -2.0*xin[n-1]+xin[n];
   FOR i := 2 to n-1 DO BEGIN
      xout[i] := -2.0*xin[i-1]+xin[i]+2.0*xin[i+1]
   END
END;

PROCEDURE atsub(xin: glnarray; VAR xout: glnarray; n: integer);
(* Programs using routine ATSUB must define the type
TYPE
   glnarray = ARRAY [1..n] OF real;
in the main routine. *)
VAR
   i: integer;
BEGIN
   xout[1] := xin[1]-2.0*xin[2];
   xout[n] := 2.0*xin[n-1]+xin[n];
   FOR i := 2 to n-1 DO BEGIN
      xout[i] := 2.0*xin[i-1]+xin[i]-2.0*xin[i+1]
   END
END;

(*$I MODFILE.PAS *)
(*$I SPARSE.PAS *)

BEGIN
   FOR i := 1 to n DO BEGIN
      x[i] := 0.0;
      b[i] := 1.0
   END;
   b[1] := 3.0;
   b[n] := -1.0;
   sparse(b,n,x,rsq);
   writeln('sum-squared residual:',rsq:15);
   writeln;
   writeln('solution vector:');
   FOR ii := 1 to (n DIV 5) DO BEGIN
      FOR i := (5*(ii-1)+1) to (5*ii) DO write(x[i]:12:6);
      writeln
   END;
   IF ((n MOD 5) > 0) THEN BEGIN
      FOR i := 1 to (n MOD 5) DO write(x[5*(n DIV 5)+i]:12:6)
   END;
   writeln;
   asub(x,bcmp,n);
   writeln;
   writeln('press RETURN to continue...');
   readln;
   writeln('test of solution vector:');
   writeln('a*x':9,'b':12);
   FOR i := 1 to n DO BEGIN
      writeln(bcmp[i]:12:6,b[i]:12:6)
   END
END.
PROGRAM d12r8(input,output,dfile);
(* driver for routine SPCTRM *)
CONST
   m=16;
   m4=64;   (* m4=4*m *)
TYPE
   glmarray = ARRAY [1..m] OF real;
   gl4marray = ARRAY [1..m4] OF real;
   gldarray = gl4marray;
VAR
   j,k : integer;
   ovrlap : boolean;
   p,q,w2 : glmarray;
   w1 : gl4marray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I FOUR1.PAS *)

(*$I SPCTRM.PAS *)

BEGIN
   glopen(dfile,'spctrl.dat');
   k := 8;
   ovrlap := true;
   spctrm(p,m,k,ovrlap,w1,w2);
   close(dfile);
   glopen(dfile,'spctrl.dat');
   k := 16;
   ovrlap := false;
   spctrm(q,m,k,ovrlap,w1,w2);
   close(dfile);
   writeln('Spectrum of data in file SPCTRL.DAT');
   writeln(' ':13,'overlapped ',' ':5,'non-overlapped');
   FOR j := 1 to m DO BEGIN
      writeln (j:6,' ':5,p[j]:13,' ':5,q[j]:13)
   END
END.
PROGRAM d13r16(input,output,dfile);
(* driver for routine SPEAR *)
CONST
   ndat=20;
   nmon=12;
TYPE
   cityname = string[15];
   monthname = string[4];
   narray = ARRAY [1..ndat] OF real;
   glsarray = narray;
VAR
   d,probd,probrs,rs,zd : real;
   i,j : integer;
   ave,data1,data2 : narray;
   wksp1,wksp2,zlat : narray;
   rays : ARRAY [1..ndat,1..nmon] OF real;
   city : ARRAY [1..ndat] OF cityname;
   mon : ARRAY [1..nmon] OF monthname;
   txt : string[64];
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SORT2.PAS *)

(*$I CRANK.PAS *)

(*$I ERFCC.PAS *)

(*$I GAMMLN.PAS *)

(*$I BETACF.PAS *)

(*$I BETAI.PAS *)

(*$I SPEAR.PAS *)

BEGIN
   glopen(dfile,'table2.dat');
   readln(dfile);
   readln(dfile,txt);
   read(dfile,city[1]);
   FOR i := 1 to nmon DO read(dfile,mon[i]);
   readln(dfile);
   readln(dfile);
   FOR i := 1 to ndat DO BEGIN
      read(dfile,city[i]);
      FOR j := 1 to nmon DO read(dfile,rays[i,j]);
      read(dfile,ave[i]);
      read(dfile,zlat[i]);
      readln(dfile)
   END;
   close(dfile);
   writeln(txt);
   write(' ':15);
   FOR i := 1 to 12 DO write(mon[i]:4);
   writeln;
   FOR i := 1 to ndat DO BEGIN
      write(city[i]:15);
      FOR j := 1 to 12 DO write(round(rays[i,j]):4);
      writeln
   END;
(* check temperature correlations between different months *)
   writeln;
   writeln('Are sunny summer places also sunny winter places?');
   writeln('Check correlation of sampled U.S. solar radiation');
   writeln('(july with other months)');
   writeln;
   writeln('month','d':9,'st. dev.':14,'probd':11,
         'spearman-r':15,'probrs':10);
   FOR i := 1 to ndat DO data1[i] := rays[i,1];
   FOR j := 1 to 12 DO BEGIN
      FOR i := 1 to ndat DO data2[i] := rays[i,j];
      spear(data1,data2,ndat,wksp1,wksp2,d,zd,probd,rs,probrs);
      writeln(mon[j],d:12:2,zd:12:6,probd:12:6,
            rs:13:6,probrs:12:6)
   END
END.
PROGRAM d3r12(input,output);
(* driver for routine SPLIE2 *)
CONST
   m=10;
   n=10;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glmbyn = ARRAY [1..m,1..n] OF real;
VAR
   i,j : integer;
   x1x2 : real;
   x1,x2 : glnarray;
   y,y2 : glmbyn;

(*$I MODFILE.PAS *)
(*$I SPLINE.PAS *)

(*$I SPLIE2.PAS *)

BEGIN
   FOR i := 1 to m DO BEGIN
      x1[i] := 0.2*i
   END;
   FOR i := 1 to n DO BEGIN
      x2[i] := 0.2*i
   END;
   FOR i := 1 to m DO BEGIN
      FOR j := 1 to n DO BEGIN
         x1x2 := x1[i]*x2[j];
         y[i,j] := sqr(x1x2)
      END
   END;
   splie2(x1,x2,y,m,n,y2);
   writeln;
   writeln('second derivatives from SPLIE2');
   writeln('natural spline assumed');
   FOR i := 1 to 5 DO BEGIN
      FOR j := 1 to 5 DO write(y2[i,j]:12:6);
      writeln
   END;
   writeln;
   writeln('actual second derivatives');
   FOR i := 1 to 5 DO BEGIN
      FOR j := 1 to 5 DO BEGIN
         x1x2 := x1[i]*x2[i];
         y2[i,j] := 2.0*sqr(x1[i])
      END;
      FOR j := 1 to 5 DO write(y2[i,j]:12:6);
      writeln
   END
END.
PROGRAM d3r13(input,output);
(* driver for routine SPLIN2 *)
CONST
   m=10;
   n=10;
TYPE
   glnarray = ARRAY [1..n] OF real;
   glmbyn = ARRAY [1..m,1..n] OF real;
VAR
   f,ff,x1x2,xx1,xx2 : real;
   i,j : integer;
   x1,x2 : glnarray;
   y,y2 : glmbyn;

(*$I MODFILE.PAS *)
(*$I SPLINT.PAS *)

(*$I SPLINE.PAS *)

(*$I SPLIE2.PAS *)

(*$I SPLIN2.PAS *)

BEGIN
   FOR i := 1 to m DO BEGIN
      x1[i] := 0.2*i
   END;
   FOR i := 1 to n DO BEGIN
      x2[i] := 0.2*i
   END;
   FOR i := 1 to m DO BEGIN
      FOR j := 1 to n DO BEGIN
         x1x2 := x1[i]*x2[j];
         y[i,j] := x1x2*exp(-x1x2)
      END
   END;
   splie2(x1,x2,y,m,n,y2);
   writeln('x1':9,'x2':12,'splin2':14,'actual':12);
   FOR i := 1 to 10 DO BEGIN
      xx1 := 0.1*i;
      xx2 := sqr(xx1);
      splin2(x1,x2,y,y2,m,n,xx1,xx2,f);
      x1x2 := xx1*xx2;
      ff := x1x2*exp(-x1x2);
      writeln(xx1:12:6,xx2:12:6,f:12:6,ff:12:6)
   END
END.
PROGRAM d3r3 (input,output);
(* driver for routine SPLINE *)
CONST
   n=20;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..n] OF real;
VAR
   i : integer;
   yp1,ypn : real;
   x,y,y2 : glnarray;

(*$I MODFILE.PAS *)
(*$I SPLINE.PAS *)

BEGIN
   writeln ('second-derivatives for sin(x) from 0 to pi');
(* generate array for interpolation *)
   FOR i := 1 to 20 DO BEGIN
      x[i] := i*pi/n;
      y[i] := sin(x[i])
   END;
(* calculate 2nd derivative with spline *)
   yp1 := cos(x[1]);
   ypn := cos(x[n]);
   spline(x,y,n,yp1,ypn,y2);
(* test result *)
   writeln ('spline':23,'actual':16);
   writeln ('number':11,'2nd deriv':14,'2nd deriv':16);
   FOR i := 1 to n DO BEGIN
      writeln (i:8,y2[i]:16:6,-sin(x[i]):16:6)
   END
END.
PROGRAM d3r4 (input,output);
(* driver for routine SPLINT *)
LABEL 1;
CONST
   np=10;
   pi=3.1415926;
TYPE
   glnarray = ARRAY [1..np] OF real;
VAR
   i,nfunc : integer;
   f,x,y,yp1,yp2,ypn : real;
   xa,ya,y2 : glnarray;

(*$I MODFILE.PAS *)
(*$I SPLINE.PAS *)

(*$I SPLINT.PAS *)

BEGIN
   FOR nfunc := 1 to 2 DO BEGIN
      IF (nfunc = 1) THEN BEGIN
         writeln;
         writeln ('sine function from 0 to pi');
         FOR i := 1 to np DO BEGIN
            xa[i] := i*pi/np;
            ya[i] := sin(xa[i])
         END;
         yp1 := cos(xa[1]);
         ypn := cos(xa[np])
      END ELSE IF (nfunc = 2) THEN BEGIN
         writeln;
         writeln ('exponential function from 0 to 1');
         FOR i := 1 to np DO BEGIN
            xa[i] := 1.0*i/np;
            ya[i] := exp(xa[i])
         END;
         yp1 := exp(xa[1]);
         ypn := exp(xa[np])
      END ELSE BEGIN
         GOTO 1
      END;
(* call spline to get second derivatives *)
      spline(xa,ya,np,yp1,yp2,y2);
(* call splint FOR interpolations *)
      writeln;
      writeln ('x':9,'f(x)':13,'interpolation':17);
      FOR i := 1 to 10 DO BEGIN
         IF (nfunc = 1) THEN BEGIN
            x := (-0.05+i/10.0)*pi;
            f := sin(x)
         END ELSE IF (nfunc = 2) THEN BEGIN
            x := -0.05+i/10.0;
            f := exp(x)
         END;
         splint(xa,ya,y2,np,x,y);
         writeln (x:12:6,f:12:6,y:12:6)
      END;
      writeln;
      writeln ('***********************************');
      writeln ('press RETURN');
      readln
   END;
1:
END.
PROGRAM d2r8(input,output,dfile);
(* driver for routine SVBKSB *)
LABEL 10,99;
CONST
   np=20;
   mp=20;
TYPE
   glnparray = ARRAY [1..np] OF real;
   glmparray = ARRAY [1..mp] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
VAR
   j,k,l,m,n : integer;
   wmax,wmin : real;
   a,b,u : glmpbynp;
   v : glnpbynp;
   w,x : glnparray;
   c : glmparray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SVDCMP.PAS *)

(*$I SVBKSB.PAS *)

BEGIN
   glopen(dfile,'matrx1.dat');
10:   readln(dfile);
   readln(dfile);
   readln(dfile,n,m);
   readln(dfile);
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n-1 DO read(dfile,a[k,l]);
      readln(dfile,a[k,n])
   END;
   readln(dfile);
   FOR l := 1 to m DO BEGIN
      FOR k := 1 to n-1 DO read(dfile,b[k,l]);
      readln(dfile,b[n,l])
   END;
(* copy a into u *)
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n DO BEGIN
         u[k,l] := a[k,l]
      END
   END;
(* decompose matrix a *)
   svdcmp(u,n,n,np,np,w,v);
(* find maximum singular value *)
   wmax := 0.0;
   FOR k := 1 to n DO BEGIN
      IF  (w[k] > wmax) THEN  wmax := w[k]
   END;
(* define "small" *)
   wmin := wmax*(1.0e-6);
(* zero the "small" singular values *)
   FOR k := 1 to n DO BEGIN
      IF  (w[k] < wmin) THEN  w[k] := 0.0
   END;
(* backsubstitute for each right-hand side vector *)
   FOR l := 1 to m DO BEGIN
      writeln;
      writeln('Vector number ',l:2);
      FOR k := 1 to n DO BEGIN
         c[k] := b[k,l]
      END;
      svbksb(u,w,v,n,n,np,np,c,x);
      writeln ('    solution vector is:');
      FOR k := 1 to n-1 DO write(x[k]:12:6);
      writeln(x[n]:12:6);
      writeln ('    original right-hand side vector:');
      FOR k := 1 to n-1 DO write(c[k]:12:6);
      writeln(c[n]:12:6);
      writeln ('    result of (matrix)*(sol''n vector):');
      FOR k := 1 to n DO BEGIN
         c[k] := 0.0;
         FOR j := 1 to n DO BEGIN
            c[k] := c[k]+a[k,j]*x[j]
         END
      END;
      FOR k := 1 to n-1 DO write(c[k]:12:6);
      writeln(c[n]:12:6)
   END;
   writeln ('***********************************');
   IF eof(dfile) THEN GOTO 99;
   writeln ('press RETURN for next problem');
   readln;
   GOTO 10;
99:   close(dfile)
END.
PROGRAM d2r9 (input,output,dfile);
(* driver for routine SVDCMP *)
LABEL 10,99;
CONST
   np = 20;
   mp = 20;
TYPE
   glnparray = ARRAY [1..np] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
VAR
   j,k,l,m,n : integer;
   a,u : glmpbynp;
   v : glnpbynp;
   w : glnparray;
   dfile : text;

(*$I MODFILE.PAS *)
(*$I SVDCMP.PAS *)

BEGIN
(* read input matrices *)
   glopen(dfile,'matrx3.dat');
10:   readln(dfile);
   readln(dfile);
   readln(dfile,m,n);
   readln(dfile);
(* copy original matrix into u *)
   FOR k := 1 to m DO BEGIN
      FOR l := 1 to n DO BEGIN
         read(dfile,a[k,l]);
         u[k,l] := a[k,l]
      END;
      readln(dfile)
   END;
   IF (n > m) THEN BEGIN
      FOR k := m+1 to n DO BEGIN
         FOR l := 1 to n DO BEGIN
            a[k,l] := 0.0;
            u[k,l] := 0.0
         END
      END;
      m := n
   END;
(* perform decomposition *)
   svdcmp(u,m,n,np,np,w,v);
(* write results *)
   writeln ('Decomposition matrices:');
   writeln ('Matrix u');
   FOR k := 1 to m DO BEGIN
      FOR l := 1 to n DO BEGIN
         write (u[k,l]:12:6);
      END;
      writeln
   END;
   writeln ('Diagonal of matrix w');
   FOR k := 1 to n DO BEGIN
      write(w[k]:12:6)
   END;
   writeln;
   writeln ('Matrix v-transpose');
   FOR k := 1 to n DO BEGIN
      FOR l := 1 to n DO BEGIN
         write(v[l,k]:12:6)
      END;
      writeln
   END;
   writeln;
   writeln ('Check product against original matrix:');
   writeln ('Original matrix:');
   FOR k := 1 to m DO BEGIN
      FOR l := 1 to n DO BEGIN
         write (a[k,l]:12:6)
      END;
      writeln
   END;
   writeln ('Product u*w*(v-transpose):');
   FOR k := 1 to m DO BEGIN
      FOR l := 1 to n DO BEGIN
         a[k,l] := 0.0;
         FOR j := 1 to n DO BEGIN
            a[k,l] := a[k,l]+u[k,j]*w[j]*v[l,j]
         END
      END;
      FOR l := 1 to n-1 DO write(a[k,l]:12:6);
      writeln(a[k,n]:12:6);
   END;
   writeln ('***********************************');
   IF eof(dfile) THEN GOTO 99;
   writeln ('press RETURN for next problem');
   readln;
   GOTO 10;
99:   close(dfile)
END.
PROGRAM d14r4(input,output);
(* driver for routine SVDFIT *)
(* polynomial fit *)
CONST
   npt=100;
   spread=0.02;
   npol=5;
   mp=npt;
   np=npol;
   ncvm=npol;
TYPE
   glndata = ARRAY [1..npt] OF real;
   glmma = ARRAY [1..npol] OF real;
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glmpbynp = ARRAY [1..mp,1..np] OF real;
   glnparray = glmma;
   glmparray = glndata;
   (* The previous 2 declarations are required because of the absence *)
   (* of conformant arrays in most Pascal implementations. *)
   glcvm = ARRAY [1..ncvm,1..ncvm] OF real;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   chisq : real;
   i,idum : integer;
   x,y,sig : glndata;
   a : glmma;
   cvm : glcvm;
   u : glmpbynp;
   v : glnpbynp;
   w : glnparray;

(*$I MODFILE.PAS *)
(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I SVDCMP.PAS *)

(*$I SVBKSB.PAS *)

(*$I SVDVAR.PAS *)

PROCEDURE func(x: real; VAR p: glmma; mma: integer);
(* This is essentially FPOLY renamed. *)
VAR
   j: integer;
BEGIN
   p[1] := 1.0;
   FOR j := 2 to mma DO p[j] := p[j-1]*x
END;

(*$I SVDFIT.PAS *)

BEGIN
   gliset := 0;
   idum := -911;
   FOR i := 1 to npt DO BEGIN
      x[i] := 0.02*i;
      y[i] := 1.0+x[i]*(2.0+x[i]*(3.0+x[i]*(4.0+x[i]*5.0)));
      y[i] := y[i]*(1.0+spread*gasdev(idum));
      sig[i] := y[i]*spread
   END;
   svdfit(x,y,sig,npt,a,npol,u,v,w,mp,np,chisq);
   svdvar(v,npol,np,w,cvm,npol);
   writeln;
   writeln('Polynomial fit:');
   FOR i := 1 to npol DO BEGIN
      writeln(a[i]:12:6,'  +-',sqrt(cvm[i,i]):10:6)
   END;
   writeln('Chi-squared',chisq:12:6);
END.
PROGRAM d14r5(input,output);
(* driver for routine SVDVAR *)
CONST
   np=6;
   ma=3;
   ncvm=ma;
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF real;
   glnparray = ARRAY [1..np] OF real;
   glcvm = ARRAY [1..ncvm,1..ncvm] OF real;
   mabyma = ARRAY [1..ma,1..ma] OF real;
VAR
   i,j : integer;
   v : glnpbynp;
   w : glnparray;
   cvm : glcvm;
   tru : mabyma;

(*$I MODFILE.PAS *)
(*$I SVDVAR.PAS *)

BEGIN
   w[1] := 0.0; w[2] := 1.0; w[3] := 2.0;
   w[4] := 3.0; w[5] := 4.0; w[6] := 5.0;
   v[1,1] := 1.0; v[1,2] := 1.0; v[1,3] := 1.0;
   v[1,4] := 1.0; v[1,5] := 1.0; v[1,6] := 1.0;
   v[2,1] := 2.0; v[2,2] := 2.0; v[2,3] := 2.0;
   v[2,4] := 2.0; v[2,5] := 2.0; v[2,6] := 2.0;
   v[3,1] := 3.0; v[3,2] := 3.0; v[3,3] := 3.0;
   v[3,4] := 3.0; v[3,5] := 3.0; v[3,6] := 3.0;
   v[4,1] := 4.0; v[4,2] := 4.0; v[4,3] := 4.0;
   v[4,4] := 4.0; v[4,5] := 4.0; v[4,6] := 4.0;
   v[5,1] := 5.0; v[5,2] := 5.0; v[5,3] := 5.0;
   v[5,4] := 5.0; v[5,5] := 5.0; v[5,6] := 5.0;
   v[6,1] := 6.0; v[6,2] := 6.0; v[6,3] := 6.0;
   v[6,4] := 6.0; v[6,5] := 6.0; v[6,6] := 6.0;
   tru[1,1] := 1.25; tru[1,2] := 2.5; tru[1,3] := 3.75;
   tru[2,1] := 2.5; tru[2,2] := 5.0; tru[2,3] := 7.5;
   tru[3,1] := 3.75; tru[3,2] := 7.5; tru[3,3] := 11.25;
   writeln;
   writeln('matrix v');
   FOR i := 1 to np DO BEGIN
      FOR j := 1 to np DO write(v[i,j]:12:6);
      writeln
   END;
   writeln;
   writeln('vector w');
   FOR i := 1 to np DO write(w[i]:12:6);
   writeln;
   svdvar(v,ma,np,w,cvm,ncvm);
   writeln;
   writeln('covariance matrix from svdvar');
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(cvm[i,j]:12:6);
      writeln
   END;
   writeln;
   writeln('expected covariance matrix');
   FOR i := 1 to ma DO BEGIN
      FOR j := 1 to ma DO write(tru[i,j]:12:6);
      writeln
   END
END.
