Category : Pascal Source Code
Archive   : PAS-SCI.ZIP
Filename : SOLVIT.PAS

PROGRAM solvit(output);
(* Feb 8.0, 81 *)
(* from Pascal Programs for Scientists and Engineers *)
(* Alan R. Miller, Sybex *)
(* n x n inverse hilbert matrix *)
(* solution is 1 1 1 1 1 *)
(* double precision version *)

(* Pascal program to perform simultaneous solution *)
(* by Gauss-Jordan elimination *)

CONST
maxr = 10;
maxc = 10;

TYPE
ary = ARRAY[1..maxr] OF real;
arys = ARRAY[1..maxc] OF real;
ary2s = ARRAY[1..maxr, 1..maxc] OF real;

VAR
y : arys;
coef : arys;
a, b : ary2s;
n, m, i, j : integer;
error : boolean;

PROCEDURE gaussj
(VAR b : ary2s; (* square matrix of coefficients *)
y : arys; (* constant vector *)
VAR coef : arys; (* solution vector *)
ncol : integer; (* order of matrix *)
VAR error: boolean); (* true if matrix singular *)

(* Gauss Jordan matrix inversion and solution *)
(* Adapted from McCormick *)
(* Feb 8, 81 *)
(* B(N,N) coefficient matrix, becomes inverse *)
(* Y(N) original constant vector *)
(* W(N,M) constant vector(s) become solution vector *)
(* DETERM is the determinant *)
(* ERROR = 1 if singular *)
(* INDEX(N,3) *)
(* NV is number of constant vectors *)

LABEL
99,98;

VAR
w : ARRAY[1..maxc, 1..maxc] OF real;
index: ARRAY[1..maxc, 1..3] OF integer;
i, j, k, l, nv, irow, icol, n, l1 : integer;
determ, pivot, hold, sum, t, ab, big: real;

PROCEDURE swap(VAR a, b: real);

VAR
hold: real;

BEGIN (* swap *)
hold := a;
a := b;
b := hold
END (* procedure swap *);


BEGIN (* Gauss-Jordan main program *)
(* put gausj2 here *)
(* PROCEDURE gausj2;

LABEL 98;

VAR
i, j, k, l, l1: integer;

BEGIN *) (* procedure gausj2 *)
(* actual start of gaussj *)
error := false;
nv := 1 (* single constant vector *);
n := ncol;
FOR i := 1 TO n DO
BEGIN
w[i, 1] := y[i] (* copy constant vector *);
index[i, 3] := 0
END;
determ := 1.0;
FOR i := 1 TO n DO
BEGIN
(* search for largest element *)
big := 0.0;
FOR j := 1 TO n DO
BEGIN
IF index[j, 3] <> 1 THEN
BEGIN
FOR k := 1 TO n DO
BEGIN
IF index[k, 3] > 1 THEN
BEGIN
writeln(' ERROR: matrix singular');
error := true;
GOTO 98 (* abort *)
END;
IF index[k, 3] < 1 THEN
IF abs(b[j, k]) > big THEN
BEGIN
irow := j;
icol := k;
big := abs(b[j, k])
END
END (* k loop *)
END
END (* j loop *);
index[icol, 3] := index[icol, 3] + 1;
index[i, 1] := irow;
index[i, 2] := icol;

(* gausj3 *) (* further subdivision of gaussj *);
(*
PROCEDURE gausj3;

VAR
l: integer;

BEGIN *) (* procedure gausj3 *)

(* interchange rows to put pivot on diagonal *)
IF irow <> icol THEN
BEGIN
determ := - determ;
FOR l := 1 TO n DO
swap(b[irow, l], b[icol, l]);
IF nv > 0 THEN
FOR l := 1 TO nv DO
swap(w[irow, l], w[icol, l])
END (* if irow <> icol *)
(* END *) (* gausj3 *);

(* divide pivot row by pivot column *)
pivot := b[icol, icol];
determ := determ * pivot;
b[icol, icol] := 1.0;
FOR l := 1 TO n DO
b[icol, l] := b[icol, l] / pivot;
IF nv > 0 THEN
FOR l := 1 TO nv DO
w[icol, l] := w[icol, l] / pivot;
(* reduce nonpivot rows *)
FOR l1 := 1 TO n DO
BEGIN
IF l1 <> icol THEN
BEGIN
t := b[l1, icol];
b[l1, icol] := 0.0;
FOR l := 1 TO n DO
b[l1, l] := b[l1, l] - b[icol, l] * t;
IF nv > 0 THEN
FOR l := 1 TO nv DO
w[l1, l] := w[l1, l] - w[icol, l] * t;
END (* IF l1 <> icol *)
END
END (* i loop *);
98:
(* END *) (* gausj2 *);

(* gausj2 *) (* first half of gaussj *);
IF error THEN GOTO 99;
(* interchange columns *)
FOR i := 1 TO n DO
BEGIN
l := n - i + 1;
IF index[l, 1] <> index[l, 2] THEN
BEGIN
irow := index[l, 1];
icol := index[l, 2];
FOR k := 1 TO n DO
swap(b[k, irow], b[k, icol])
END (* if index *)
END (* i loop *);
FOR k := 1 TO n DO
IF index[k, 3] <> 1 THEN
BEGIN
writeln(' ERROR: matrix singular');
error := true;
GOTO 99 (* abort *)
END;
FOR i := 1 TO n DO
coef[i] := w[i, 1];
99:
END (* procedure gaussj *);


PROCEDURE get_data(VAR a : ary2s;
VAR y : arys;
VAR n, m : integer);

(* setup n-by-n hilbert matrix *)

VAR
i, j : integer;

BEGIN
FOR i := 1 TO n DO
BEGIN
a[n,i] := 1.0/(n + i - 1);
a[i,n] := a[n,i]
END;
a[n,n] := 1.0/(2*n -1);
FOR i := 1 TO n DO
BEGIN
y[i] := 0.0;
FOR j := 1 TO n DO
y[i] := y[i] + a[i,j]
END;
writeln;
IF n < 7 THEN
BEGIN
FOR i:= 1 TO n DO
BEGIN
FOR j:= 1 TO m DO
write( a[i,j] :7:5, ' ');
writeln( ' : ', y[i] :7:5)
END;
writeln
END (* if n<7 *)
END (* procedure get_data *);

PROCEDURE write_data;

(* print out the answers *)

VAR
i : integer;

BEGIN
FOR i := 1 TO m DO
write( coef[i] :13:9);
writeln;
END (* write_data *);

(* PROCEDURE gaussj
(VAR a : ary2s;
y : arys;
VAR coef : arys;
ncol : integer;
VAR error : boolean);
extern; *)
(* F GAUSSJ.PAS *)

BEGIN (* main program *)
a[1,1] := 1.0;
n := 2;
m := n;
REPEAT
get_data (a, y, n, m);
FOR i := 1 TO n DO
FOR j := 1 TO n DO
b[i,j] := a[i,j] (* setup work array *);
gaussj (b, y, coef, n, error);
IF not error THEN write_data;
n := n+1;
m := n
UNTIL n > maxr
END.