{   prog4a.pas -- Solve Least Linear System by Gauss Seidel Method
    Written Feb 20, 2001 by Eric Olson for Mathematics 483

    For more information see Step 13 in First Steps in Numerical
    Analysis, 2nd Edition, by Hosking, Joe, Joyce and Turner. }

program main;

const MAXn=40;

type matrix=array [0..MAXn-1] of array [0..MAXn-1] of real;
    vector=array [0..MAXn-1] of real;

var a : matrix;
    b : vector;
    n : integer;

function gauss_seidel(var x:vector):integer;
    var y:vector;
        s:real;
        iter:integer;
        i,j:integer;
    label gauss_seidel_exit;
begin
    iter:=0;
    repeat
        iter:=iter+1;
        s:=0.0;
        for i:=0 to n-1 do begin
            y[i]:=x[i]
        end;
        for i:=0 to n-1 do begin
            x[i]:=b[i];
            for j:=0 to i-1 do begin
                x[i]:=x[i]-a[i][j]*x[j]
            end;
            for j:=i+1 to n-1 do begin
                x[i]:=x[i]-a[i][j]*y[j]
            end;
            x[i]:=x[i]/a[i][i];
            s:=s+abs(y[i]-x[i])
        end;
        if (s>1.0e+32) or (iter>10000) then begin
            gauss_seidel:=0;
            goto gauss_seidel_exit
        end;
    until s<=1.0e-5;
    gauss_seidel:=iter;
gauss_seidel_exit:
end;

procedure outmatrix(var a:matrix);
    var i,j:integer;
begin
    write('[');
    for i:=0 to n-1 do begin
        writeln;
        for j:=0 to n-1 do begin
            write(a[i][j],' ')
        end
    end;
    writeln(']')
end;

procedure outvector(var x:vector);
    var j:integer;
begin
    write('[');
    for j:=0 to n-1 do begin
        write(x[j],' ')
    end;
    writeln(']');
end;

var x:vector;
    i,j,k,m:integer;

begin
    repeat
        read(n);
        if n<=0 then halt(0);
        if n>MAXn then begin
            writeln('Capacity exceeded!');
            halt(1)
        end;
        for i:=0 to n-1 do begin
            for j:=0 to n-1 do begin
                read(a[i][j])
            end
        end;
        write('A = ');
        outmatrix(a);
        writeln;
        read(m);
        for i:=0 to m-1 do begin
            for j:=0 to n-1 do begin
                read(b[j])
            end;
            write('b = '); outvector(b);
            for j:=0 to n-1 do x[j]:=b[j];
            k:=gauss_seidel(x);
            if k<>0 then begin
                write('x = '); outvector(x);
                writeln('Solution converged in',k:5,' iterations.');
                writeln;
            end else begin
                writeln('Failed convergence.');
                writeln;
            end
        end
    until FALSE
end.