Title: How to invert a matrix
type
RCOMat = array of array of Extended;
var
DimMat: integer;
procedure InvertMatrix(var aa: RCOMat);
var
numb, nula1, ipiv, indxr, indxc: array of Integer;
i, j, l, kod, jmax, k, ll, icol, irow: Integer;
amax, d, c, pomos, big, dum, pivinv: Double;
ind: Boolean;
begin
for j := 0 to Pred(DimMat) do ipiv[j] := 0;
irow := 1;
icol := 1;
for i := 0 to Pred(DimMat) do
begin
big := 0;
for j := 0 to Pred(DimMat) do
begin
if (ipiv[j] 1) then
begin
for k := 0 to Pred(DimMat) do
begin
if (ipiv[k] = 0) then
if (Abs(aa[j, k]) = big) then
begin
big := Abs(aa[j, k]);
irow := j;
icol := k;
end
else;
ipiv[icol] := ipiv[icol] 1;
if (irow icol) then
begin
for l := 0 to Pred(DimMat) do
begin
dum := aa[irow, l];
aa[irow, l] := aa[icol, l];
aa[icol, l] := dum;
end;
for l := 0 to Pred(DimMat) do
begin
dum := aa[irow DimMat 1, l];
aa[irow DimMat 1, l] := aa[icol DimMat 1, l];
aa[icol DimMat 1, l] := dum;
end;
indxr[i] := irow;
indxc[i] := icol;
if (aa[icol, icol] = 0) then;
pivinv := 1.0 / aa[icol, icol];
aa[icol, icol] := 1.0;
for l := 0 to Pred(DimMat) do aa[icol, l] := aa[icol, l] * pivinv;
for l := 0 to Pred(DimMat) do aa[icol DimMat 1, l] :=
aa[icol DimMat 1, l] * pivinv;
for ll := 0 to Pred(DimMat) do
begin
if (ll icol) then
begin
dum := aa[ll, icol];
aa[ll, icol] := 0.0;
for l := 0 to Pred(DimMat) do aa[ll, l] := aa[ll, l] - aa[icol, l] * dum;
for l := 0 to Pred(DimMat) do aa[ll DimMat 1, l] :=
aa[ll DimMat 1, l] - aa[icol DimMat 1, l] * dum;
end;
end;
end;
for l := Pred(DimMat) downto 0 do
begin
if (indxr[l] indxc[l]) then
begin
for k := 0 to Pred(DimMat) do
begin
dum := aa[k, indxr[l]];
aa[k, indxr[l]] := aa[k, indxc[l]];
aa[k, indxc[l]] := dum;
end;