246 lines
6.1 KiB
Plaintext
246 lines
6.1 KiB
Plaintext
{$N-}
|
||
program Hilb;
|
||
{
|
||
|
||
The program performs simultaneous solution by Gauss-Jordan
|
||
elimination.
|
||
|
||
--------------------------------------------------
|
||
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
|
||
--------------------------------------------------
|
||
|
||
INSTRUCTIONS
|
||
1. Compile and run the program using the $N- (Numeric Processing :
|
||
Software) compiler directive.
|
||
2. if you have a math coprocessor in your computer, compile and run the
|
||
program using the $N+ (Numeric Processing : Hardware) compiler
|
||
directive. Compare the speed and precision of the results to those
|
||
of example 1.
|
||
}
|
||
|
||
const
|
||
maxr = 10;
|
||
maxc = 10;
|
||
|
||
type
|
||
{$IFOPT N+} { use extended type if using 80x87 }
|
||
real = extended;
|
||
{$ENDIF}
|
||
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 *)
|
||
|
||
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 *)
|
||
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;
|
||
exit; (* 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;
|
||
|
||
(* 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 *)
|
||
|
||
(* 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 *);
|
||
|
||
if error then exit;
|
||
(* 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;
|
||
exit; (* abort *)
|
||
end;
|
||
for i := 1 to n do
|
||
coef[i] := w[i, 1];
|
||
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 *);
|
||
|
||
|
||
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.
|
||
|