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.