对称逆矩阵的计算方法

：JohnsonGuo, 时间：2000-12-28 18:17:34, ID：427389 给个思路：
1。通过传统方法，把矩阵[A I]通过行变换得到[I A-1]，具体方法其实就是高斯消元法。
2。对矩阵进行LDL'分解，于是A-1 = L'-1 D-1 L-1
其中L-1就是把L的对角线下元素取负即可，D-1就是对角线元素取倒数，L'-1 = (L-1)'
第2种方法对对称阵的效果会更好一些。
　
下面给一段程序（还没有优化）：
type
TMatrix = array of array of Extended;
ENotSquare = class(Exception);
EDimNotMatch = class(Exception);
...
procedure LDLT(m: TMatrix; var l, d: TMatrix); //LDL'分解
var
n, i, j, k: Integer;
s: Extended;
begin
n := High(m);
if High(m[0]) <> n then
Raise ENotSquare.Create('LDL''分解需要方阵！');
SetLength(l, n + 1, n + 1);
SetLength(d, n + 1, n + 1);
for i := 0 to n do
for j := 0 to n do begin
l[i, j] := Ord(i = j);
d[i, j] := 0;
end;
for i := 0 to n do begin
s := m[i, i];
for j := 0 to i - 1 do
s := s - Sqr(l[i, j]) * d[j, j];
d[i, i] := s;
for k := i + 1 to n do begin
s := m[k, i];
for j := 0 to i - 1 do
s := s - l[k, j] * l[i, j] * d[j, j];
l[k, i] := s / d[i, i];
end;
end;
end;
procedure Transpose(m: TMatrix; var Trans: TMatrix); //矩阵转置
var
i, j, n0, n1: Integer;
begin
n0 := High(m);
n1 := High(m[0]);
SetLength(Trans, n1 + 1, n0 + 1);
for i := 0 to n0 do
for j := 0 to n1 do
Trans[j, i] := m[i, j];
end;
procedure Multiply(m1, m2: TMatrix; var Res: TMatrix); //矩阵乘法
var
i, j, k, n0, n1, n2: Integer;
s: Extended;
begin
n0 := High(m1);
n1 := High(m1[0]);
if n1 <> High(m2) then
Raise EDimNotMatch.Create('矩阵维数不匹配！');
n2 := High(m2[0]);
SetLength(Res, n0 + 1, n2 + 1);
for i := 0 to n0 do
for j := 0 to n2 do begin
s := 0;
for k := 0 to n1 do
s := s + m1[i, k] * m2[k, j];
Res[i, j] := s;
end;
end;
procedure InverseSymmetry(m: TMatrix; var Inv: TMatrix); //对称阵求逆
var
l, d, lt, temp: TMatrix;
i, j, k, n: Integer;
s: Extended;
begin
n := High(m);
LDLT(m ,l, d);
SetLength(lt, n + 1, n + 1);
for i := 0 to n do
for j := 0 to n do
lt[i, j] := Ord(i = j);
for i := 0 to n do begin
d[i, i] := 1 / d[i, i];
for j := i + 1 to n do begin
s := 0;
for k := 0 to j - 1 do
s := s + l[j, k] * lt[k, i];
lt[j, i] := -s;
end;
end;
Transpose(lt, l);
Multiply(l, d, temp);
Multiply(temp, lt, Inv);
end;