PROGRAM d17r2(input,output);
(* driver for routine ADI *)
LABEL 1;
CONST
   jmax=11;
   pi=3.1415926;
TYPE
   double = real;
   gljmax = ARRAY [1..jmax,1..jmax] OF double;
VAR
   alim,alpha,beta,eps : double;
   i,j,k,mid,twotok : integer;
   a,b,c,d,e,f,g,u : gljmax;

(*$I ADI.PAS *)

BEGIN
   FOR i := 1 to jmax DO BEGIN
      FOR j := 1 to jmax DO BEGIN
         a[i,j] := -1.0;
         b[i,j] := 2.0;
         c[i,j] := -1.0;
         d[i,j] := -1.0;
         e[i,j] := 2.0;
         f[i,j] := -1.0;
         g[i,j] := 0.0;
         u[i,j] := 0.0
      END
   END;
   mid := (jmax DIV 2)+1;
   g[mid,mid] := 2.0;
   alpha := 2.0*(1.0-cos(pi/jmax));
   beta := 2.0*(1.0-cos((jmax-1)*pi/jmax));
   alim := ln(4.0*jmax/pi);
   k := 0;
   twotok := 1;
   REPEAT
      k := k+1;
      twotok := 2*twotok;
   UNTIL twotok >= alim;
   eps := 1.0e-4;
   adi(a,b,c,d,e,f,g,u,jmax,k,alpha,beta,eps);
   writeln('ADI 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
         g[i,j] := -4.0*u[i,j]+u[i+1,j]
            +u[i-1,j]+u[i,j-1]+u[i,j+1]
      END;
      write(' ':7);
      FOR j := 2 to (jmax-1) DO write(g[i,j]:7:2);
      writeln
   END
END.