奇异值分解
procedure SVDCMP(A:matrx2; M, N:integer;var W:array of real; V:matrx2);
label 1,2,3;
var
RV1:array [0..100] of real;
I,J,L,ITS,JJ,NM,K:integer;
G,S,SCALE1,C,F,H,X,Y,Z,Sgn,ANORM,AAAAA:real;
begin
If M < N Then
ShowMessage('You must augment A with extra zero rows.');
G:=0;
SCALE1:=0;
ANORM:=0;
For I:=1 To N do
begin
L:=I + 1;
RV1[I]:=SCALE1 * G;
G:=0;
S:=0;
SCALE1:=0;
If I <= M Then
begin
For K:=I To M do
SCALE1:=SCALE1 + Abs(A[K, I]);
If SCALE1 <> 0 Then
begin
For K:=I To M do
begin
A[K, I]:=A[K, I] / SCALE1;
S:=S + A[K, I] * A[K, I];
end;
F:=A[I, I];
If F > 0 Then
SGN:=1
else
SGN:=-1;
G:=-Sqrt(S) * Sgn;
H:=F * G - S;
A[I, I]:=F - G;
If I <> N Then
begin
For J:=L To N DO
begin
S:=0;
For K:=I To M do
S:=S + A[K, I] * A[K, J];
F:=S / H;
For K:=I To M do
A[K, J]:=A[K, J] + F * A[K, I];
end;
end;
For K:=I To M do
A[K, I]:=SCALE1 * A[K, I];
end;
end;
W[I]:=SCALE1 * G;
G:=0;
S:=0;
SCALE1:=0;
If (I <= M) And (I <> N) Then
begin
For K:=L To N do
SCALE1:=SCALE1 + Abs(A[I, K]);
If SCALE1 <> 0 Then
begin
For K:=L To N do
begin
A[I, K]:=A[I, K] / SCALE1;
S:=S + A[I, K] * A[I, K];
end;
F:=A[I, L];
If F > 0 Then
SGN:=1
else
SGN:=-1;
G:=-Sqrt(S) * Sgn;
H:=F * G - S;
A[I, L]:=F - G;
For K:=L To N do
RV1[K]:=A[I, K] / H;
If I <> M Then
begin
For J:=L To M do
begin
S:=0;
For K:=L To N do
S:=S + A[J, K] * A[I, K];
For K:=L To N do
A[J, K]:=A[J, K] + S * RV1[K];
end;
end;
For K:=L To N do
A[I, K]:=SCALE1 * A[I, K];
end;
end;
If ANORM > Abs(W[I]) + Abs(RV1[I]) Then
ANORM:=ANORM
Else
ANORM:=Abs(W[I]) + Abs(RV1[I]);
end;
For I:=N DownTo 1 do
begin
If I < N Then
begin
If G <> 0 Then
begin
For J:=L To N do
V[J, I]:=(A[I, J] / A[I, L]) / G;
For J:=L To N do
begin
S:=0;
For K:=L To N do
S:=S + A[I, K] * V[K, J];
For K:=L To N do
V[K, J]:=V[K, J] + S * V[K, I];
end;
end;
For J:=L To N do
begin
V[I, J]:=0;
V[J, I]:=0;
end;
end;
V[I, I]:=1;
G:=RV1[I];
L:=I;
end;
For I:=N DownTo 1 do
begin
L:=I + 1;
G:=W[I];
If I < N Then
begin
For J:=L To N do
A[I, J]:=0;
end;
If G <> 0 Then
begin
G:=1 / G;
If I <> N Then
begin
For J:=L To N do
begin
S:=0;
For K:=L To M do
S:=S + A[K, I] * A[K, J];
F:=(S / A[I, I]) * G;
For K:=I To M do
A[K, J]:=A[K, J] + F * A[K, I];
end;
end;
For J:=I To M do
begin
A[J, I]:=A[J, I] * G;
end;
end
Else
For J:=I To M do
A[J, I]:=0;
A[I, I]:=A[I, I] + 1;
end;
For K:=N DownTo 1 do
begin
For ITS:=1 To 30 do
begin
For L:=K DownTo 1 do
begin
NM:=L - 1;
If Abs(RV1[L]) + ANORM = ANORM Then GoTo 2;
If Abs(W[NM]) + ANORM = ANORM Then GoTo 1;
end;
1: C:=0;
S:=1;
For I:=L To K do
begin
F:=S * RV1[I];
If Abs(F) + ANORM <> ANORM Then
begin
G:=W[I];
H:=Sqrt(F * F + G * G);
W[I]:=H;
H:=1 / H;
C:=(G * H);
S:=-(F * H);
For J:=1 To M do
begin
Y:=A[J, NM];
Z:=A[J, I];
A[J, NM]:=(Y * C) + (Z * S);
A[J, I]:=-(Y * S) + (Z * C);
end;
end;
end;
2: Z:=W[K];
If L = K Then
begin
If Z < 0 Then
begin
W[K]:=-Z;
For J:=1 To N do
V[J, K]:=-V[J, K];
end;
GoTo 3;
end;
If ITS = 30 Then
ShowMessage('No convergence in 30 iterations');
X:=W[L];
NM:=K - 1;
Y:=W[NM];
G:=RV1[NM];
H:=RV1[K];
F:=((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2 * H * Y);
G:=Sqrt(F * F + 1);
If F > 0 Then
Sgn:=1
else
Sgn:=-1;
F:=((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn)) - H)) / X;
C:=1;
S:=1;
For J:=L To NM do
begin
I:=J + 1;
G:=RV1[I];
Y:=W[I];
H:=S * G;
G:=G * C;
Z:=Sqrt(F * F + H * H);
RV1[J]:=Z;
C:=F / Z;
S:=H / Z;
F:=(X * C) + (G * S);
G:=-(X * S) + (G * C);
H:=Y * S;
Y:=Y * C;
For NM:=1 To N do
begin
X:=V[NM, J];
Z:=V[NM, I];
V[NM, J]:=(X * C) + (Z * S);
V[NM, I]:=-(X * S) + (Z * C);
end;
Z:=Sqrt(F * F + H * H);
W[J]:=Z;
If Z <> 0 Then
begin
Z:=1 / Z;
C:=F * Z;
S:=H * Z;
end;
F:=(C * G) + (S * Y);
X:=-(S * G) + (C * Y);
For NM:=1 To M do
begin
Y:=A[NM, J];
Z:=A[NM, I];
A[NM, J]:=(Y * C) + (Z * S);
A[NM, I]:=-(Y * S) + (Z * C);
end;
end;
RV1[L]:=0;
RV1[K]:=F;
W[K]:=X;
end;
3: AAAAA:=1;
end;
end;
procedure SVBKSB(U:matrx2; W:array of real; V:matrx2;
M, N:integer; B:array of real;var X:array of reaL);
var
TMP:array[0..100] of real;
I,J,JJ:integer;
S:real;
begin
For J:=1 To N do
begin
S:=0;
If W[J] <> 0 Then
begin
For I:=1 To M do
S:=S + U[I, J] * B[I];
S:=S / W[J];
end;
TMP[J]:=S;
end;
For J:=1 To N do
begin
S:=0;
For JJ:=1 To N do
S:=S + V[J, JJ] * TMP[JJ];
X[J]:=S;
end;
end;