?? matrices.pas
字號:
Lbound, Ubound : Integer;
A_inv : TMatrix;
var Det : Float) : Integer;
var
I, J, K : Integer; { Loop variables }
Ik, Jk : Integer; { Pivot coordinates }
Pvt : Float; { Pivot }
T : Float; { Auxiliary variable }
PRow, PCol : TIntVector; { Store line and column of pivot }
MCol : TVector; { Stores a column of the matrix }
begin
DimVector(PRow, Ubound);
DimVector(PCol, Ubound);
DimVector(MCol, Ubound);
{ Copy A into A_inv }
CopyMatrix(A_inv, A, Lbound, Lbound, Ubound, Ubound);
Det := 1.0;
K := Lbound;
while K <= Ubound do
begin
{ Search for largest pivot in submatrix A_inv[K..Ubound, K..Ubound] }
Pvt := A_inv[K,K];
Ik := K;
Jk := K;
for I := K to Ubound do
for J := K to Ubound do
if Abs(A_inv[I,J]) > Abs(Pvt) then
begin
Pvt := A_inv[I,J];
Ik := I;
Jk := J;
end;
{ Pivot too small ==> quasi-singular matrix }
if Abs(Pvt) < MACHEP then
begin
InvDet := MAT_SINGUL;
Exit;
end;
{ Save pivot position }
PRow[K] := Ik;
PCol[K] := Jk;
{ Update determinant }
Det := Det * Pvt;
if Ik <> K then Det := - Det;
if Jk <> K then Det := - Det;
{ Exchange current row (K) with pivot row (Ik) }
if Ik <> K then
SwapRows(Ik, K, A_inv, Lbound, Ubound);
{ Exchange current column (K) with pivot column (Jk) }
if Jk <> K then
SwapCols(Jk, K, A_inv, Lbound, Ubound);
{ Store col. K of A_inv into MCol and set this col. to 0 }
for I := Lbound to Ubound do
if I <> K then
begin
MCol[I] := A_inv[I, K];
A_inv[I, K] := 0.0;
end
else
begin
MCol[I] := 0.0;
A_inv[I, K] := 1.0;
end;
{ Transform pivot row }
for J := Lbound to Ubound do
A_inv[K,J] := A_inv[K,J] / Pvt;
{ Transform other rows }
for I := Lbound to Ubound do
if I <> K then
begin
T := MCol[I];
for J := Lbound to Ubound do
A_inv[I,J] := A_inv[I,J] - T * A_inv[K,J];
end;
Inc(K);
end;
{ Exchange lines of inverse matrix }
for I := Ubound downto Lbound do
begin
Ik := PCol[I];
if Ik <> I then
SwapRows(Ik, I, A_inv, Lbound, Ubound);
end;
{ Exchange columns of inverse matrix }
for J := Ubound downto Lbound do
begin
Jk := PRow[J];
if Jk <> J then
SwapCols(Jk, J, A_inv, Lbound, Ubound);
end;
InvDet := MAT_OK;
end;
function Det(A : TMatrix; Lbound, Ubound : Integer) : Float;
var
I, J, K : Integer; { Loop variables }
Ik, Jk : Integer; { Pivot coordinates }
Pvt : Float; { Pivot }
T : Float; { Auxiliary variable }
PRow, PCol : TIntVector; { Store line and column of pivot }
MCol : TVector; { Stores a column of the matrix }
D : Float; { Determinant }
begin
DimVector(PRow, Ubound);
DimVector(PCol, Ubound);
DimVector(MCol, Ubound);
D := 1.0;
K := Lbound;
while K <= Ubound do
begin
{ Search for largest pivot in submatrix A[K..Ubound, K..Ubound] }
Pvt := A[K,K];
Ik := K;
Jk := K;
for I := K to Ubound do
for J := K to Ubound do
if Abs(A[I,J]) > Abs(Pvt) then
begin
Pvt := A[I,J];
Ik := I;
Jk := J;
end;
{ Pivot too small ==> quasi-singular matrix }
if Abs(Pvt) < MACHEP then
begin
Det := 0.0;
Exit;
end;
{ Save pivot position }
PRow[K] := Ik;
PCol[K] := Jk;
{ Update determinant }
D := D * Pvt;
if Ik <> K then D := - D;
if Jk <> K then D := - D;
{ Exchange current row (K) with pivot row (Ik) }
if Ik <> K then
SwapRows(Ik, K, A, Lbound, Ubound);
{ Exchange current column (K) with pivot column (Jk) }
if Jk <> K then
SwapCols(Jk, K, A, Lbound, Ubound);
{ Store col. K of A into MCol and set this col. to 0 }
for I := Lbound to Ubound do
if I <> K then
begin
MCol[I] := A[I, K];
A[I, K] := 0.0;
end
else
begin
MCol[I] := 0.0;
A[I, K] := 1.0;
end;
{ Transform pivot row }
for J := Lbound to Ubound do
A[K,J] := A[K,J] / Pvt;
{ Transform other rows }
for I := Lbound to Ubound do
if I <> K then
begin
T := MCol[I];
for J := Lbound to Ubound do
A[I,J] := A[I,J] - T * A[K,J];
end;
Inc(K);
end;
Det := D;
end;
function Cholesky(A : TMatrix; Lbound, Ubound : Integer;
L : TMatrix) : Integer;
var
I, J, K : Integer;
Sum : Float;
begin
for K := Lbound to Ubound do
begin
Sum := A[K,K];
for J := Lbound to Pred(K) do
Sum := Sum - Sqr(L[K,J]);
if Sum <= 0.0 then
begin
Cholesky := MAT_NOT_PD;
Exit;
end;
L[K,K] := Sqrt(Sum);
for I := Succ(K) to Ubound do
begin
Sum := A[I,K];
for J := Lbound to Pred(K) do
Sum := Sum - L[I,J] * L[K,J];
L[I,K] := Sum / L[K,K];
end;
end;
Cholesky := MAT_OK;
end;
var { Used by LU procedures }
Index : TIntVector; { Records the row permutations }
function LU_Decomp(A : TMatrix; Lbound, Ubound : Integer) : Integer; overload;
var
I, Imax, J, K : Integer;
Pvt, T, Sum : Float;
V : TVector;
begin
DimVector(V, Ubound);
DimVector(Index, Ubound);
for I := Lbound to Ubound do
begin
Pvt := 0.0;
for J := Lbound to Ubound do
if Abs(A[I,J]) > Pvt then
Pvt := Abs(A[I,J]);
if Pvt < MACHEP then
begin
LU_Decomp := MAT_SINGUL;
Exit;
end;
V[I] := 1.0 / Pvt;
end;
for J := Lbound to Ubound do
begin
for I := Lbound to Pred(J) do
begin
Sum := A[I,J];
for K := Lbound to Pred(I) do
Sum := Sum - A[I,K] * A[K,J];
A[I,J] := Sum;
end;
Imax := 0;
Pvt := 0.0;
for I := J to Ubound do
begin
Sum := A[I,J];
for K := Lbound to Pred(J) do
Sum := Sum - A[I,K] * A[K,J];
A[I,J] := Sum;
T := V[I] * Abs(Sum);
if T > Pvt then
begin
Pvt := T;
Imax := I;
end;
end;
if J <> Imax then
begin
SwapRows(Imax, J, A, Lbound, Ubound);
V[Imax] := V[J];
end;
Index[J] := Imax;
if A[J,J] = 0.0 then
A[J,J] := MACHEP;
if J <> Ubound then
begin
T := 1.0 / A[J,J];
for I := Succ(J) to Ubound do
A[I,J] := A[I,J] * T;
end;
end;
LU_Decomp := MAT_OK;
end;
procedure LU_Solve(A : TMatrix; B : TVector; Lbound, Ubound : Integer;
X : TVector); overload;
var
I, Ip, J, K : Integer;
Sum : Float;
begin
K := Pred(Lbound);
CopyVector(X, B, Lbound, Ubound);
for I := Lbound to Ubound do
begin
Ip := Index[I];
Sum := X[Ip];
X[Ip] := X[I];
if K >= Lbound then
for J := K to Pred(I) do
Sum := Sum - A[I,J] * X[J]
else if Sum <> 0.0 then
K := I;
X[I] := Sum;
end;
for I := Ubound downto Lbound do
begin
Sum := X[I];
if I < Ubound then
for J := Succ(I) to Ubound do
Sum := Sum - A[I,J] * X[J];
X[I] := Sum / A[I,I];
end;
end;
function LU_Decomp(A : TCompMatrix; Lbound, Ubound : Integer) : Integer; overload;
var
I, Imax, J, K : Integer;
C, Pvt, T : Float;
Sum : Complex;
V : TVector;
begin
DimVector(V, Ubound);
DimVector(Index, Ubound);
for I := Lbound to Ubound do
begin
Pvt := 0.0;
for J := Lbound to Ubound do
begin
C := CAbs(A[I,J]);
if C > Pvt then Pvt := C;
end;
if Pvt < MACHEP then
begin
LU_Decomp := MAT_SINGUL;
Exit;
end;
V[I] := 1.0 / Pvt;
end;
for J := Lbound to Ubound do
begin
for I := Lbound to Pred(J) do
begin
Sum := A[I,J];
for K := Lbound to Pred(I) do
{ Sum := Sum - A[I,K] * A[K,J]; }
Sum := CSub(Sum, CMult(A[I,K], A[K,J]));
A[I,J] := Sum;
end;
Imax := 0;
Pvt := 0.0;
for I := J to Ubound do
begin
Sum := A[I,J];
for K := Lbound to Pred(J) do
{ Sum := Sum - A[I,K] * A[K,J]; }
Sum := CSub(Sum, CMult(A[I,K], A[K,J]));
A[I,J] := Sum;
T := V[I] * CAbs(Sum);
if T > Pvt then
begin
Pvt := T;
Imax := I;
end;
end;
if J <> Imax then
begin
{ SwapRows(Imax, J, A, Lbound, Ubound); }
for K := Lbound to Ubound do
Swap(A[Imax,K], A[J,K]);
V[Imax] := V[J];
end;
Index[J] := Imax;
if CAbs(A[J,J]) = 0.0 then
A[J,J] := Cmplx(MACHEP, MACHEP);
if J <> Ubound then
for I := Succ(J) to Ubound do
{ A[I,J] := A[I,J] / A[J,J]; }
A[I,J] := CDiv(A[I,J], A[J,J]);
end;
LU_Decomp := MAT_OK;
end;
procedure LU_Solve(A : TCompMatrix; B : TCompVector;
Lbound, Ubound : Integer; X : TCompVector); overload;
var
I, Ip, J, K : Integer;
Sum : Complex;
begin
K := Pred(Lbound);
{ CopyVector(X, B, Lbound, Ubound); }
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -