program dr24le(input, output);
{dr24le.pas == linear equations by conjugate gradients

          Copyright 1988 J.C.Nash
}

const
     banner = 'dr24le.pas -- linear equations by conjugate gradients';

#include "common.inc.pas"  // constants, types & procedures
#include "calceps.pas"   
#include "matrixin.pas"  
#include "vectorin.pas"  
#include "matmul.pas"    
#include "alg24.pas"     
#include "resids.pas"    

var
  A : rmatrix;
  Y : rvector; {RHS}
  Bvec : rvector; {solution}
  avec : smatvec; {for matrixin only}
  sym : boolean; {to tell if matrix symmetric}
  ch : char;
  i, j, n, itcount : integer;
  ssmin, t, s : real;

begin
  startup;
  write('Order of problem = ');
  readln(infile,n);
  if not isConsole then writeln(n);
  writeln(confile,'order of matrix = ',n);
  writeln('Coefficient matrix');
  matrixin(n, n, A, avec, sym);
  if not sym then halt;
  writeln('RHS vector');
  writeln(confile,'RHS vector');
  vectorin(n, Y);
  writeln('Initial guess for solution');
  writeln(confile,'Initial guess for solution');
  vectorin(n, Bvec);
  itcount:=10*n; {safety setting}
  lecg( n, A, Y, Bvec, itcount, ssmin);
  writeln('Solution after ',itcount,' iterations. Est. sumsquares ',ssmin);
  writeln(confile,
        'Solution after ',itcount,' iterations. Est. sumsquares ',ssmin);
  for i:=1 to n do
  begin
    write(Bvec[i]:10:5,' ');
    write(confile,Bvec[i]:10:5,' ');
    if (7 * (i div 7) = i) and (i<n) then
    begin
    writeln;
    writeln(confile);
    end;
  end;
  writeln;
  writeln(confile);
  s:=resids(n, n, A, Y, Bvec, true);
    {local files auto closed when procedure ends}
end. {dr24le.pas}
