```PROGRAM Solvgv;
{ simultaneous solution by Gauss-Jordan
elimination with multiple constant vectors }
{ From Borland Pascal Programs for Scientists and Engineers }
{ by Alan R. Miller, Copyright C 1993, SYBEX Inc }
USES WinCrt; { Crt for non-windows version}

CONST
Maxr = 7;
Maxc = 7;

TYPE
FltPt = Real;
Ary2s = ARRAY[1..Maxr, 1..Maxc] OF FltPt;

VAR
A, Y: Ary2s;
N, Nvec: Integer;
Error: Boolean;
Determ: FltPt;

PROCEDURE Get_Data(VAR A: Ary2s;
VAR Y: Ary2s;
VAR N, Nvec: Integer);
{ Get values for N, Nvec, and arrays A, Y }

VAR
I, J: Integer;

BEGIN
WriteLn;
REPEAT
Write(' How many equations? ');
UNTIL N < Maxr;
IF N > 1 THEN
BEGIN
Write(' How many constant vectors? ');
FOR I := 1 TO N DO
BEGIN
FOR J := 1 TO N DO
BEGIN
Write(J:3, ': ');
END;
IF Nvec > 0 THEN
BEGIN
FOR J := 1 TO Nvec DO
BEGIN
Write(', C: ');
END;
WriteLn;
END
ELSE {matrix determinant only}
WriteLn;
END; { I loop }
WriteLn;
Write('           Matrix');
IF Nvec > 0 THEN
Write('          Constants');
WriteLn;
FOR I:= 1 TO N  DO
BEGIN
FOR J:= 1 TO N DO
Write(A[I,J]:7:4, '  ');
FOR J := 1 TO Nvec DO
Write(' :', Y[I,J]:7:4);
WriteLn
END; { I loop }
WriteLn
END  { if N>1 }
END; { procedure Get_Data }

PROCEDURE Write_Data;
{ print out the answers }
VAR
I, J: Integer;

BEGIN
IF Nvec > 0 THEN
BEGIN
WriteLn(' Solution');
FOR I := 1 TO N DO
BEGIN
FOR J := 1 TO Nvec DO
Write(Y[I,J]:9:5);
WriteLn
END
END { if }
ELSE
BEGIN
WriteLn('     Inverse');
FOR I := 1 TO N DO
BEGIN
FOR J := 1 TO N DO
Write('  ',A[I,J]:10);
WriteLn
END;
WriteLn;
Write(' Determinant is ', Determ:12)
END; { ELSE }
WriteLn
END; { Write_Data }

PROCEDURE Gaussjv
(VAR B : Ary2s; { square matrix of coefficients }
VAR W : Ary2s;  { constant vector matrix }
VAR Determ: FltPt;  { the determinant }
Ncol      : Integer;  { order of matrix }
Nv        : Integer;  { number of constants }
VAR Error: Boolean); { True of matrix singular }

{      Gauss Jordan matrix inversion and solution
B(N,N) coefficient matrix, becomes inverse
W(N,M) constant vector(s) become solution vector
Determ is determinant
Error = 1 if singular
Index(N,3)
Nv is number of constant vectors }

VAR
Index: ARRAY[1..Maxc, 1..3] OF Integer;
I, J, K, L, Irow, Icol, N, L1: Integer;
Pivot, T, Big: FltPt;

PROCEDURE Swap(VAR A, B: FltPt);

VAR
Hold: FltPt;

BEGIN  { Swap }
Hold := A;
A := B;
B := Hold
END; { procedure Swap }

BEGIN     { Gauss-Jordan main program }
Error := False;
N := Ncol;
FOR I := 1 TO N DO
Index[I, 3] := 0;
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  {procedure Gaussj }
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 }
{ 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  {procedure Gaussj }
END;
END; { procedure Gaussjv }

BEGIN  { main program }
WriteLn;
WriteLn(' Simultaneous solution by Gauss-Jordan');
WriteLn
(' Multiple constant vectors, or matrix inverse');
REPEAT
Get_Data(A, Y, N, Nvec);
IF N > 1 THEN
BEGIN
Gaussjv(A, Y, Determ, N, Nvec, Error);
IF NOT Error THEN  Write_Data
END
UNTIL N < 2;
DoneWinCrt  { for Windows version only }
END.

```