2015-07-02
520
Подпрограмма SVDDecomposition осуществляет SVD-разложение прямоугольной матрицы размером MxN. На выходе подпрограмма возвращает массив сингулярных значений, упорядоченных по убыванию, и, по требованию, матрицы U и V T, причем возможно как возвращение только левых и правых сингулярных векторов, так и полных матриц размером MxM и NxN (в зависимости от параметров UNeeded и VTNeeded). Обратите внимание, что в переменной VT алгоритм возвращает не матрицу V, а матрицу V T.
Более подробно параметры подпрограммы описаны в комментариях к ней.
(***********************************************************************
Ýòîò êîä ñãåíåðèðîâàí òðàíñëÿòîðîì AlgoPascal â ðàìêàõ ïðîåêòà ALGLIB
Óñëîâèÿ, íà êîòîðûõ ðàñïðîñòðàíÿåòñÿ ýòîò êîä, ìîæíî ïîëó÷èòü ïî àäðåñó
|
|
|
http://alglib.sources.ru/copyrules.php
***********************************************************************)
uses Math, Ap;
(*
Ýòè ïîäïðîãðàììû äîëæåí îïðåäåëèòü ïðîãðàììèñò:
procedure ToBidiagonal(var A: TReal2DArray;
M: Integer;
N: Integer;
var TauQ: TReal1DArray;
var TauP: TReal1DArray);
procedure UnpackQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TAUQ: TReal1DArray;
QColumns: Integer;
var Q: TReal2DArray);
procedure MultiplyByQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauQ: TReal1DArray;
var Z: TReal2DArray;
ZRows: Integer;
ZColumns: Integer;
FromTheRight: Boolean;
DoTranspose: Boolean);
procedure UnpackPTFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TAUP: TReal1DArray;
PTRows: Integer;
var PT: TReal2DArray);
procedure MultiplyByPFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauP: TReal1DArray;
var Z: TReal2DArray;
ZRows: Integer;
ZColumns: Integer;
FromTheRight: Boolean;
DoTranspose: Boolean);
procedure UnpackDiagonalsFromBidiagonal(const B: TReal2DArray;
M: Integer;
N: Integer;
var IsUpper: Boolean;
var D: TReal1DArray;
var E: TReal1DArray);
function BidiagonalSVDDecomposition(var D: TReal1DArray;
E: TReal1DArray;
N: Integer;
IsUpper: Boolean;
IsFractionalAccuracyRequired: Boolean;
var U: TReal2DArray;
NRU: Integer;
var C: TReal2DArray;
NCC: Integer;
var VT: TReal2DArray;
NCVT: Integer):Boolean;
procedure QRDecomposition(var A: TReal2DArray;
M: Integer;
N: Integer;
var Tau: TReal1DArray);
procedure UnpackQFromQR(const QR: TReal2DArray;
M: Integer;
N: Integer;
const Tau: TReal1DArray;
QColumns: Integer;
var Q: TReal2DArray);
procedure LQDecomposition(var A: TReal2DArray;
M: Integer;
N: Integer;
var Tau: TReal1DArray);
procedure UnpackQFromLQ(const LQ: TReal2DArray;
M: Integer;
N: Integer;
const Tau: TReal1DArray;
QRows: Integer;
var Q: TReal2DArray);
procedure CopyAndTranspose(const A: TReal2DArray;
IS1: Integer;
IS2: Integer;
JS1: Integer;
JS2: Integer;
var B: TReal2DArray;
ID1: Integer;
ID2: Integer;
JD1: Integer;
JD2: Integer);
procedure InplaceTranspose(var A: TReal2DArray;
I1: Integer;
I2: Integer;
J1: Integer;
J2: Integer;
var WORK: TReal1DArray);
procedure MatrixMatrixMultiply(const A: TReal2DArray;
AI1: Integer;
AI2: Integer;
AJ1: Integer;
AJ2: Integer;
TransA: Boolean;
const B: TReal2DArray;
BI1: Integer;
BI2: Integer;
BJ1: Integer;
BJ2: Integer;
TransB: Boolean;
|
|
|
Alpha: Double;
var C: TReal2DArray;
CI1: Integer;
CI2: Integer;
CJ1: Integer;
CJ2: Integer;
Beta: Double;
var WORK: TReal1DArray);
procedure CopyMatrix(const A: TReal2DArray;
IS1: Integer;
IS2: Integer;
JS1: Integer;
JS2: Integer;
var B: TReal2DArray;
ID1: Integer;
ID2: Integer;
JD1: Integer;
JD2: Integer);
*)
function SVDDecomposition(A: TReal2DArray;
M: Integer;
N: Integer;
UNeeded: Integer;
VTNeeded: Integer;
AdditionalMemory: Integer;
var W: TReal1DArray;
var U: TReal2DArray;
var VT: TReal2DArray):Boolean;forward;
(*************************************************************************
Ñèíãóëÿðíîå ðàçëîæåíèå ïðÿìîóãîëüíîé ìàòðèöû.
Àëãîðèòì âû÷èñëÿåò ñèíãóëÿðíîå ðàçëîæåíèå ìàòðèöû ðàçìåðîì MxN:
A = U * S * V^T
Àëãîðèòì íàõîäèò ñèíãóëÿðíûå çíà÷åíèÿ è, îïöèîíàëüíî, ìàòðèöû U è V^T.
Ïðè ýòîì âîçìîæíî êàê íàõîæäåíèå ïåðâûõ min(M,N) ñòîëáöîâ ìàòðèöû U è
ñòðîê ìàòðèöû V^T (ñèíãóëÿðíûõ âåêòîðîâ), òàê è ïîëíûõ ìàòðèö U è V^T
(ðàçìåðîì MxM è NxN).
Îáðàòèòå âíèìàíèå, ÷òî ïîäïðîãðàììà âîçâðàùàåò ìàòðèöó V^T, à íå V.
Âõîäíûå ïàðàìåòðû:
A - ðàçëàãàåìàÿ ìàòðèöà.
Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..M, 1..N]
M - ÷èñëî ñòðîê â ìàòðèöå A
N - ÷èñëî ñòîëáöîâ â ìàòðèöå A
UNeeded - 0, 1 èëè 2. Ïîäðîáíåå ñì. îïèñàíèå ïàðàìåòðà U
VTNeeded - 0, 1 èëè 2. Ïîäðîáíåå ñì. îïèñàíèå ïàðàìåòðà VT
AdditionalMemory-
åñëè ïàðàìåòð:
* ðàâåí 0, òî àëãîðèòì íå èñïîëüçóåò äîïîëíèòåëüíóþ
ïàìÿòü (ìåíüøå òðåáîâàíèÿ ê ðåñóðñàì, íèæå ñêîðîñòü).
* ðàâåí 1, òî àëãîðèòì ìîæåò èñïîëüçîâàòü äîïîëíèòåëüíóþ
ïàìÿòü â ðàçìåðå min(M,N)*min(M,N) âåùåñòâåííûõ ÷èñåë.
 ðÿäå ñëó÷àåâ óâåëè÷èâàåò ñêîðîñòü àëãîðèòìà.
* ðàâåí 2, òî àëãîðèòì ìîæåò èñïîëüçîâàòü äîïîëíèòåëüíóþ
ïàìÿòü â ðàçìåðå M*min(M,N) âåùåñòâåííûõ ÷èñåë.  ðÿäå
ñëó÷àåâ ýòî ïîçâîëÿåò äîáèâàòüñÿ ìàêñèìàëüíîé ñêîðîñòè.
|
|
|
Ðåêîìåíäóåìîå çíà÷åíèå ïàðàìåòðà: 2.
Âûõîäíûå ïàðàìåòðû:
W - ñîäåðæèò ñèíãóëÿðíûå çíà÷åíèÿ, óïîðÿäî÷åííûå ïî óáûâàíèþ.
U - åñëè UNeeded=0, íå èçìåíÿåòñÿ. Ëåâûå ñèíãóëÿðíûå âåêòîðû
íå âû÷èñëÿþòñÿ.
åñëè UNeeded=1, ñîäåðæèò ëåâûå ñèíãóëÿðíûå âåêòîðû (ïåðâûå
min(M,N) ñòîëáöîâ ìàòðèöû U). Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ
[1..M, 1..Min(M,N)].
åñëè UNeeded=2, ñîäåðæèò ïîëíóþ ìàòðèöó U. Ìàññèâ ñ íóìåðà-
öèåé ýëåìåíòîâ [1..M, 1..M].
VT - åñëè VTNeeded=0, íå èçìåíÿåòñÿ. Ïðàâûå ñèíãóëÿðíûå âåêòîðû
íå âû÷èñëÿþòñÿ.
åñëè VTNeeded=1, ñîäåðæèò ïðàâûå ñèíãóëÿðíûå âåêòîðû
(ïåðâûå min(M,N) ñòðîê ìàòðèöû V^T). Ìàññèâ ñ íóìåðàöèåé
ýëåìåíòîâ [1..min(M,N), 1..N].
åñëè VTNeeded=2, ñîäåðæèò ïîëíóþ ìàòðèöó V^T. Ìàññèâ ñ
|
|
|
íóìåðàöèåé ýëåìåíòîâ [1..N, 1..N].
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************)
function SVDDecomposition(A: TReal2DArray;
M: Integer;
N: Integer;
UNeeded: Integer;
VTNeeded: Integer;
AdditionalMemory: Integer;
var W: TReal1DArray;
var U: TReal2DArray;
var VT: TReal2DArray):Boolean;
var
TauQ: TReal1DArray;
TauP: TReal1DArray;
Tau: TReal1DArray;
E: TReal1DArray;
WORK: TReal1DArray;
T2: TReal2DArray;
IsUpper: Boolean;
MinMN: Integer;
NCU: Integer;
NRVT: Integer;
NRU: Integer;
NCVT: Integer;
I: Integer;
J: Integer;
IM1: Integer;
SM: Double;
begin
A:= DynamicArrayCopy(A);
Result:= True;
if (M=0) or (N=0) then
begin
Exit;
end;
Assert((UNeeded>=0) and (UNeeded<=2), 'SVDDecomposition: wrong parameters!');
Assert((VTNeeded>=0) and (VTNeeded<=2), 'SVDDecomposition: wrong parameters!');
Assert((AdditionalMemory>=0) and (AdditionalMemory<=2), 'SVDDecomposition: wrong parameters!');
//
// initialize
//
MinMN:= Min(M, N);
SetLength(W, MinMN+1);
NCU:= 0;
NRU:= 0;
if UNeeded=1 then
begin
NRU:= M;
NCU:= MinMN;
SetLength(U, NRU+1, NCU+1);
end;
if UNeeded=2 then
begin
NRU:= M;
NCU:= M;
SetLength(U, NRU+1, NCU+1);
end;
NRVT:= 0;
NCVT:= 0;
if VTNeeded=1 then
begin
NRVT:= MinMN;
NCVT:= N;
SetLength(VT, NRVT+1, NCVT+1);
end;
if VTNeeded=2 then
begin
NRVT:= N;
NCVT:= N;
SetLength(VT, NRVT+1, NCVT+1);
end;
//
// M much larger than N
// Use bidiagonal reduction with QR-decomposition
//
if M>1.6*N then
begin
if UNeeded=0 then
begin
//
// No left singular vectors to be computed
//
QRDecomposition(A, M, N, Tau);
I:=2;
while I<=N do
begin
J:=1;
while J<=I-1 do
begin
A[I,J]:= 0;
Inc(J);
end;
Inc(I);
end;
ToBidiagonal(A, N, N, TauQ, TauP);
UnpackPTFromBidiagonal(A, N, N, TauP, NRVT, VT);
UnpackDiagonalsFromBidiagonal(A, N, N, IsUpper, W, E);
Result:= BidiagonalSVDDecomposition(W, E, N, IsUpper, False, U, 0, A, 0, VT, NCVT);
Exit;
end
else
begin
//
// Left singular vectors (may be full matrix U) to be computed
//
QRDecomposition(A, M, N, Tau);
UnpackQFromQR(A, M, N, Tau, NCU, U);
I:=2;
while I<=N do
begin
J:=1;
while J<=I-1 do
begin
A[I,J]:= 0;
Inc(J);
end;
Inc(I);
end;
ToBidiagonal(A, N, N, TauQ, TauP);
UnpackPTFromBidiagonal(A, N, N, TauP, NRVT, VT);
UnpackDiagonalsFromBidiagonal(A, N, N, IsUpper, W, E);
if AdditionalMemory<1 then
begin
//
// No additional memory can be used
//
MultiplyByQFromBidiagonal(A, N, N, TauQ, U, M, N, True, False);
Result:= BidiagonalSVDDecomposition(W, E, N, IsUpper, False, U, M, A, 0, VT, NCVT);
end
else
begin
//
// Large U. Transforming intermediate matrix T2
//
SetLength(WORK, Max(M, N)+1);
UnpackQFromBidiagonal(A, N, N, TauQ, N, T2);
CopyMatrix(U, 1, M, 1, N, A, 1, M, 1, N);
InplaceTranspose(T2, 1, N, 1, N, WORK);
Result:= BidiagonalSVDDecomposition(W, E, N, IsUpper, False, U, 0, T2, N, VT, NCVT);
MatrixMatrixMultiply(A, 1, M, 1, N, False, T2, 1, N, 1, N, True, 1.0, U, 1, M, 1, N, 0.0, WORK);
end;
Exit;
end;
end;
//
// N much larger than M
// Use bidiagonal reduction with LQ-decomposition
//
if N>1.6*M then
begin
if VTNeeded=0 then
begin
//
// No right singular vectors to be computed
//
LQDecomposition(A, M, N, Tau);
I:=1;
while I<=M-1 do
begin
J:=I+1;
while J<=M do
begin
A[I,J]:= 0;
Inc(J);
end;
Inc(I);
end;
ToBidiagonal(A, M, M, TauQ, TauP);
UnpackQFromBidiagonal(A, M, M, TauQ, NCU, U);
UnpackDiagonalsFromBidiagonal(A, M, M, IsUpper, W, E);
SetLength(WORK, M+1);
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
Result:= BidiagonalSVDDecomposition(W, E, M, IsUpper, False, A, 0, U, NRU, VT, 0);
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
Exit;
end
else
begin
//
// Right singular vectors (may be full matrix VT) to be computed
//
LQDecomposition(A, M, N, Tau);
UnpackQFromLQ(A, M, N, Tau, NRVT, VT);
I:=1;
while I<=M-1 do
begin
J:=I+1;
while J<=M do
begin
A[I,J]:= 0;
Inc(J);
end;
Inc(I);
end;
ToBidiagonal(A, M, M, TauQ, TauP);
UnpackQFromBidiagonal(A, M, M, TauQ, NCU, U);
UnpackDiagonalsFromBidiagonal(A, M, M, IsUpper, W, E);
SetLength(WORK, Max(M, N)+1);
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
if AdditionalMemory<1 then
begin
//
// No additional memory available
//
MultiplyByPFromBidiagonal(A, M, M, TauP, VT, M, N, False, True);
Result:= BidiagonalSVDDecomposition(W, E, M, IsUpper, False, A, 0, U, NRU, VT, N);
end
else
begin
//
// Large VT. Transforming intermediate matrix T2
//
UnpackPTFromBidiagonal(A, M, M, TauP, M, T2);
Result:= BidiagonalSVDDecomposition(W, E, M, IsUpper, False, A, 0, U, NRU, T2, M);
CopyMatrix(VT, 1, M, 1, N, A, 1, M, 1, N);
MatrixMatrixMultiply(T2, 1, M, 1, M, False, A, 1, M, 1, N, False, 1.0, VT, 1, M, 1, N, 0.0, WORK);
end;
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
Exit;
end;
end;
//
// M<=N
// We can use inplace transposition of U to get rid of columnwise operations
//
if M<=N then
begin
ToBidiagonal(A, M, N, TauQ, TauP);
UnpackQFromBidiagonal(A, M, N, TauQ, NCU, U);
UnpackPTFromBidiagonal(A, M, N, TauP, NRVT, VT);
UnpackDiagonalsFromBidiagonal(A, M, N, IsUpper, W, E);
SetLength(WORK, M+1);
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
Result:= BidiagonalSVDDecomposition(W, E, MinMN, IsUpper, False, A, 0, U, NRU, VT, NCVT);
InplaceTranspose(U, 1, NRU, 1, NCU, WORK);
Exit;
end;
//
// Simple bidiagonal reduction
//
ToBidiagonal(A, M, N, TauQ, TauP);
UnpackQFromBidiagonal(A, M, N, TauQ, NCU, U);
UnpackPTFromBidiagonal(A, M, N, TauP, NRVT, VT);
UnpackDiagonalsFromBidiagonal(A, M, N, IsUpper, W, E);
if (AdditionalMemory<2) or (UNeeded=0) then
begin
//
// We cant use additional memory or there is no need in such operations
//
Result:= BidiagonalSVDDecomposition(W, E, MinMN, IsUpper, False, U, NRU, A, 0, VT, NCVT);
end
else
begin
//
// We can use additional memory
//
SetLength(T2, MinMN+1, M+1);
CopyAndTranspose(U, 1, M, 1, MinMN, T2, 1, MinMN, 1, M);
Result:= BidiagonalSVDDecomposition(W, E, MinMN, IsUpper, False, U, 0, T2, M, VT, NCVT);
CopyAndTranspose(T2, 1, MinMN, 1, M, U, 1, M, 1, MinMN);
end;
end;
Приведение прямоугольной матрицы к двухдиагональной форме
(***********************************************************************
Ýòîò êîä ñãåíåðèðîâàí òðàíñëÿòîðîì AlgoPascal â ðàìêàõ ïðîåêòà ALGLIB
Óñëîâèÿ, íà êîòîðûõ ðàñïðîñòðàíÿåòñÿ ýòîò êîä, ìîæíî ïîëó÷èòü ïî àäðåñó
http://alglib.sources.ru/copyrules.php
***********************************************************************)
uses Math, Ap;
(*
Ýòè ïîäïðîãðàììû äîëæåí îïðåäåëèòü ïðîãðàììèñò:
procedure GenerateReflection(var X: TReal1DArray;
N: Integer;
var Tau: Double);
procedure ApplyReflectionFromTheLeft(var C: TReal2DArray;
Tau: Double;
const V: TReal1DArray;
M1: Integer;
M2: Integer;
N1: Integer;
N2: Integer;
var WORK: TReal1DArray);
procedure ApplyReflectionFromTheRight(var C: TReal2DArray;
Tau: Double;
const V: TReal1DArray;
M1: Integer;
M2: Integer;
N1: Integer;
N2: Integer;
var WORK: TReal1DArray);
*)
procedure ToBidiagonal(var A: TReal2DArray;
M: Integer;
N: Integer;
var TauQ: TReal1DArray;
var TauP: TReal1DArray);forward;
procedure UnpackQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauQ: TReal1DArray;
QColumns: Integer;
var Q: TReal2DArray);forward;
procedure MultiplyByQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauQ: TReal1DArray;
var Z: TReal2DArray;
ZRows: Integer;
ZColumns: Integer;
FromTheRight: Boolean;
DoTranspose: Boolean);forward;
procedure UnpackPTFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauP: TReal1DArray;
PTRows: Integer;
var PT: TReal2DArray);forward;
procedure MultiplyByPFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauP: TReal1DArray;
var Z: TReal2DArray;
ZRows: Integer;
ZColumns: Integer;
FromTheRight: Boolean;
DoTranspose: Boolean);forward;
procedure UnpackDiagonalsFromBidiagonal(const B: TReal2DArray;
M: Integer;
N: Integer;
var IsUpper: Boolean;
var D: TReal1DArray;
var E: TReal1DArray);forward;
(*************************************************************************
Ïðèâåäåíèå ïðÿìîóãîëüíîé ìàòðèöû ê äâóõäèàãîíàëüíîìó âèäó
Àëãîðèòì ïðèâîäèò ïðÿìîóãîëüíóþ ìàòðèöó A ê äâóõäèàãîíàëüíîìó âèäó
îðòîãîíàëüíûìè ïðåîáðàçîâàíèÿìè P è Q, òàêèìè, ÷òî A = Q*B*P'
Âõîäíûå ïàðàìåòðû:
A - èñõîäíàÿ ìàòðèöà A. Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..M, 1..N]
M - ÷èñëî ñòðîê â ìàòðèöå
N - ÷èñëî ñòîëáöîâ â ìàòðèöå
Âûõîäíûå ïàðàìåòðû:
A - ìàòðèöû Q, B, P â óïàêîâàííîé ôîðìå (ñì. íèæå).
TauQ- ñêàëÿðíûå ìíîæèòåëè, ó÷àñòâóþùèå â ôîðìèðîâàíèè ìàòðèöû Q
TauP- ñêàëÿðíûå ìíîæèòåëè, ó÷àñòâóþùèå â ôîðìèðîâàíèè ìàòðèöû Q
 ðåçóëüòàòå ðàáîòû àëãîðèòìà ãëàâíàÿ äèàãîíàëü ìàòðèöû A è îäíà èç ïîáî÷íûõ
çàìåùàþòñÿ áèäèàãîíàëüíîé ìàòðèöåé B, à â îñòàëüíûõ ýëåìåíòàõ õðàíÿòñÿ
ýëåìåíòàðíûå îòðàæåíèÿ, ôîðìèðóþùèå ìàòðèöû Q è P ðàçìåðîì MxM è NxN
ñîîòâåòñòâåííî.
Åñëè M>=N, òî ìàòðèöà B - âåðõíÿÿ äâóõäèàãîíàëüíàÿ ðàçìåðîì MxN è õðàíèòñÿ
â ñîîòâåòñòâóþùèõ ýëåìåíòàõ ìàòðèöû A. Ìàòðèöà Q ïðåäñòàâëÿåòñÿ â âèäå
ïðîèçâåäåíèÿ ýëåìåíòàðíûõ îòðàæåíèé Q = H(1)*H(2)*...*H(n), ãäå H(i) =
= 1 - tau*v*v'. Çäåñü tau - ñêàëÿðíàÿ âåëè÷èíà, õðàíÿùàÿ â TauQ[i], à
âåêòîð v èìååò ñëåäóþùóþ ñòðóêòóðó: v(1:i-1)=0, v(i)=1, v(i+1:m) õðàíèòñÿ
â ýëåìåíòàõ A(i+1:m,i). Ìàòðèöà P èìååò âèä P = G(1)*G(2)*...*G(n-1), ãäå
G(i) = 1 - tau*u*u'. Tau õðàíèòñÿ â TauP[i], u(1:i)=0, u(i+1)=1, u(i+2:n)
õðàíèòñÿ â ýëåìåíòàõ A(i,i+2:n).
Åñëè M<N, òî ìàòðèöà B - íèæíÿÿ äâóõäèàãîíàëüíàÿ ðàçìåðîì MxN, è õðàíèòñÿ
â ñîîòâåòñòâóþùèõ ýëåìåíòàõ ìàòðèöû A. Q = H(1)*H(2)*...*H(m-1), ãäå
H(i) = 1 - tau*v*v', tau õðàíèòñÿ â TauQ, v(1:i)=0, v(i+1)=1, v(i+2:m)
õðàíèòñÿ â A(i+2:m,i). P = G(1)*G(2)*...*G(m), G(i) = 1 - tau*u*u', tau
õðàíèòñÿ â TauP, u(1:i-1)=0, u(i)=1, u(i+1:n) õðàíèòñÿ â A(i,i+1:n).
ÏÐÈÌÅÐ:
m=6, n=5 (m > n): m=5, n=6 (m < n):
(d e u1 u1 u1) (d u1 u1 u1 u1 u1)
(v1 d e u2 u2) (e d u2 u2 u2 u2)
(v1 v2 d e u3) (v1 e d u3 u3 u3)
(v1 v2 v3 d e) (v1 v2 e d u4 u4)
(v1 v2 v3 v4 d) (v1 v2 v3 e d u5)
(v1 v2 v3 v4 v5)
çäåñü vi è ui îáîçíà÷àåò âåêòîðû, ôîðìèðóþùèå H(i) è G(i), à d è e -
äèàãîíàëüíûå è âíåäèàãîíàëüíûå ýëåìåíòû ìàòðèöû B.
*************************************************************************)
procedure ToBidiagonal(var A: TReal2DArray;
M: Integer;
N: Integer;
var TauQ: TReal1DArray;
var TauP: TReal1DArray);
var
WORK: TReal1DArray;
T: TReal1DArray;
MinMN: Integer;
MaxMN: Integer;
I: Integer;
LTau: Double;
MMIP1: Integer;
NMI: Integer;
IP1: Integer;
NMIP1: Integer;
MMI: Integer;
i_: Integer;
i1_: Integer;
begin
MinMN:= Min(M, N);
MaxMN:= Max(M, N);
SetLength(Work, MaxMN+1);
SetLength(T, MaxMN+1);
SetLength(TauP, MinMN+1);
SetLength(TauQ, MinMN+1);
if M>=N then
begin
//
// Reduce to upper bidiagonal form
//
I:=1;
while I<=N do
begin
//
// Generate elementary reflector H(i) to annihilate A(i+1:m,i)
//
MMIP1:= M-I+1;
i1_:= (I) - (1);
for i_:= 1 to MMIP1 do
begin
T[i_]:= A[i_+i1_,I];
end;
GenerateReflection(T, MMIP1, LTau);
TauQ[I]:= LTau;
i1_:= (1) - (I);
for i_:= I to M do
begin
A[i_,I]:= T[i_+i1_];
end;
T[1]:= 1;
//
// Apply H(i) to A(i:m,i+1:n) from the left
//
ApplyReflectionFromTheLeft(A, LTau, T, I, M, I+1, N, WORK);
if I<N then
begin
//
// Generate elementary reflector G(i) to annihilate
// A(i,i+2:n)
//
NMI:= N-I;
IP1:= I+1;
i1_:= (IP1) - (1);
for i_:= 1 to NMI do
begin
T[i_]:= A[I,i_+i1_];
end;
GenerateReflection(T, NMI, LTau);
TauP[I]:= LTau;
i1_:= (1) - (IP1);
for i_:= IP1 to N do
begin
A[I,i_]:= T[i_+i1_];
end;
T[1]:= 1;
//
// Apply G(i) to A(i+1:m,i+1:n) from the right
//
ApplyReflectionFromTheRight(A, LTau, T, I+1, M, I+1, N, WORK);
end
else
begin
TAUP[I]:= 0;
end;
Inc(I);
end;
end
else
begin
//
// Reduce to lower bidiagonal form
//
I:=1;
while I<=M do
begin
//
// Generate elementary reflector G(i) to annihilate A(i,i+1:n)
//
NMIP1:= N-I+1;
i1_:= (I) - (1);
for i_:= 1 to NMIP1 do
begin
T[i_]:= A[I,i_+i1_];
end;
GenerateReflection(T, NMIP1, LTau);
TauP[I]:= LTau;
i1_:= (1) - (I);
for i_:= I to N do
begin
A[I,i_]:= T[i_+i1_];
end;
T[1]:= 1;
//
// Apply G(i) to A(i+1:m,i:n) from the right
//
ApplyReflectionFromTheRight(A, LTau, T, I+1, M, I, N, WORK);
if I<M then
begin
//
// Generate elementary reflector H(i) to annihilate
// A(i+2:m,i)
//
MMI:= M-I;
IP1:= I+1;
i1_:= (IP1) - (1);
for i_:= 1 to MMI do
begin
T[i_]:= A[i_+i1_,I];
end;
GenerateReflection(T, MMI, LTau);
TauQ[I]:= LTau;
i1_:= (1) - (IP1);
for i_:= IP1 to M do
begin
A[i_,I]:= T[i_+i1_];
end;
T[1]:= 1;
//
// Apply H(i) to A(i+1:m,i+1:n) from the left
//
ApplyReflectionFromTheLeft(A, LTau, T, I+1, M, I+1, N, WORK);
end
else
begin
TAUQ[I]:= 0;
end;
Inc(I);
end;
end;
end;
(*************************************************************************
×àñòè÷íàÿ "ðàñïàêîâêà" ìàòðèöû Q, ïðèâîäÿùåé ìàòðèöó A ê äâóõäèàãîíàëüíîé
ôîðìå.
Âõîäíûå ïàðàìåòðû:
QP - ìàòðèöû Q è P â óïàêîâàííîé ôîðìå.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
M - ÷èñëî ñòðîê â ìàòðèöå A
N - ÷èñëî ñòîëáöîâ â ìàòðèöå A
TAUQ - ñêàëÿðíûå ìíîæèòåëè, ôîðìèðóþùèå ìàòðèöó Q.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
QColumns- òðåáóåìîå ÷èñëî ñòîëáöîâ ìàòðèöû Q. M >= QColumns >= 0
Âûõîäíûå ïàðàìåòðû:
Q - QColumns ïåðâûõ ñòîëáöîâ ìàòðèöû Q.
Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..M, 1..QColumns]
Åñëè QColumns=0, òî ìàññèâ íå èçìåíÿåòñÿ.
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************)
procedure UnpackQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauQ: TReal1DArray;
QColumns: Integer;
var Q: TReal2DArray);
var
I: Integer;
J: Integer;
IP1: Integer;
V: TReal1DArray;
WORK: TReal1DArray;
VM: Integer;
i_: Integer;
i1_: Integer;
begin
Assert(QColumns<=M, 'UnpackQFromBidiagonal: QColumns>M!');
if (M=0) or (N=0) or (QColumns=0) then
begin
Exit;
end;
//
// init
//
SetLength(Q, M+1, QColumns+1);
SetLength(V, M+1);
SetLength(WORK, QColumns+1);
//
// prepare Q
//
I:=1;
while I<=M do
begin
J:=1;
while J<=QColumns do
begin
if I=J then
begin
Q[I,J]:= 1;
end
else
begin
Q[I,J]:= 0;
end;
Inc(J);
end;
Inc(I);
end;
if M>=N then
begin
I:=Min(N, QColumns);
while I>=1 do
begin
VM:= M-I+1;
i1_:= (I) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[i_+i1_,I];
end;
V[1]:= 1;
ApplyReflectionFromTheLeft(Q, TAUQ[I], V, I, M, 1, QColumns, WORK);
Dec(I);
end;
end
else
begin
I:=Min(M-1, QColumns-1);
while I>=1 do
begin
VM:= M-I;
IP1:= I+1;
i1_:= (IP1) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[i_+i1_,I];
end;
V[1]:= 1;
ApplyReflectionFromTheLeft(Q, TAUQ[I], V, I+1, M, 1, QColumns, WORK);
Dec(I);
end;
end;
end;
(*************************************************************************
Óìíîæåíèå íà ìàòðèöó Q, ïðèâîäÿùóþ ìàòðèöó A ê äâóõäèàãîíàëüíîé ôîðìå.
Àëãîðèòì ïîçâîëÿåò îñóùåñòâèòü ïî âûáîðó óìíîæåíèå ñïðàâà èëè ñëåâà, íà
ìàòðèöó Q èëè íà ìàòðèöó Q'.
Âõîäíûå ïàðàìåòðû:
QP - ìàòðèöû Q è P â óïàêîâàííîé ôîðìå.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
M - ÷èñëî ñòðîê â ìàòðèöå A
N - ÷èñëî ñòîëáöîâ â ìàòðèöå A
TAUQ - ñêàëÿðíûå ìíîæèòåëè, ôîðìèðóþùèå ìàòðèöó Q.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
Z - óìíîæàåìàÿ ìàòðèöà
ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..ZRows,1..ZColumns]
ZRows - ÷èñëî ñòðîê â ìàòðèöå Z. Åñëè FromTheRight=False, òî
ZRows=M, èíà÷å ìîæåò áûòü ëþáûì.
ZColumns - ÷èñëî ñòîëáöîâ â ìàòðèöå Z. Åñëè FromTheRight=True, òî
ZColumns=M, èíà÷å ìîæåò áûòü ëþáûì.
FromTheRight- óìíîæåíèå ñïðàâà èëè ñëåâà
DoTranspose - óìíîæåíèå íà Q èëè íà Q'
Âûõîäíûå ïàðàìåòðû:
Z - Ïðîèçâåäåíèå Z è Q
Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..ZRows,1..ZColumns]
Åñëè ZRows=0 èëè ZColumns=0, òî ìàññèâ íå èçìåíÿåòñÿ.
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************)
procedure MultiplyByQFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauQ: TReal1DArray;
var Z: TReal2DArray;
ZRows: Integer;
ZColumns: Integer;
FromTheRight: Boolean;
DoTranspose: Boolean);
var
I: Integer;
J: Integer;
IP1: Integer;
I1: Integer;
I2: Integer;
IStep: Integer;
V: TReal1DArray;
WORK: TReal1DArray;
VM: Integer;
Mx: Integer;
i_: Integer;
i1_: Integer;
begin
if (M<=0) or (N<=0) or (ZRows<=0) or (ZColumns<=0) then
begin
Exit;
end;
Assert(FromTheRight and (ZColumns=M) or not FromTheRight and (ZRows=M), 'MultiplyByQFromBidiagonal: incorrect Z size!');
//
// init
//
Mx:= Max(M, N);
Mx:= Max(Mx, ZRows);
Mx:= Max(Mx, ZColumns);
SetLength(V, Mx+1);
SetLength(WORK, Mx+1);
if M>=N then
begin
//
// setup
//
if FromTheRight then
begin
I1:= 1;
I2:= N;
IStep:= +1;
end
else
begin
I1:= N;
I2:= 1;
IStep:= -1;
end;
if DoTranspose then
begin
I:= I1;
I1:= I2;
I2:= I;
IStep:= -IStep;
end;
//
// Process
//
I:= I1;
repeat
VM:= M-I+1;
i1_:= (I) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[i_+i1_,I];
end;
V[1]:= 1;
if FromTheRight then
begin
ApplyReflectionFromTheRight(Z, TAUQ[I], V, 1, ZRows, I, M, WORK);
end
else
begin
ApplyReflectionFromTheLeft(Z, TAUQ[I], V, I, M, 1, ZColumns, WORK);
end;
I:= I+IStep;
until I=I2+IStep;
end
else
begin
//
// setup
//
if FromTheRight then
begin
I1:= 1;
I2:= M-1;
IStep:= +1;
end
else
begin
I1:= M-1;
I2:= 1;
IStep:= -1;
end;
if DoTranspose then
begin
I:= I1;
I1:= I2;
I2:= I;
IStep:= -IStep;
end;
//
// Process
//
if M-1>0 then
begin
I:= I1;
repeat
VM:= M-I;
IP1:= I+1;
i1_:= (IP1) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[i_+i1_,I];
end;
V[1]:= 1;
if FromTheRight then
begin
ApplyReflectionFromTheRight(Z, TAUQ[I], V, 1, ZRows, I+1, M, WORK);
end
else
begin
ApplyReflectionFromTheLeft(Z, TAUQ[I], V, I+1, M, 1, ZColumns, WORK);
end;
I:= I+IStep;
until I=I2+IStep;
end;
end;
end;
(*************************************************************************
×àñòè÷íàÿ "ðàñïàêîâêà" ìàòðèöû P, ïðèâîäÿùåé ìàòðèöó A ê äâóõäèàãîíàëüíîé
ôîðìå. Ïîäïðîãðàììà âûâîäèò òðàíñïîíèðîâàííóþ ìàòðèöó P.
Âõîäíûå ïàðàìåòðû:
QP - ìàòðèöû Q è P â óïàêîâàííîé ôîðìå.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
M - ÷èñëî ñòðîê â ìàòðèöå A
N - ÷èñëî ñòîëáöîâ â ìàòðèöå A
TAUP - ñêàëÿðíûå ìíîæèòåëè, ôîðìèðóþùèå ìàòðèöó P.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
PTRows - òðåáóåìîå ÷èñëî ñòðîê ìàòðèöû P^T. N >= PTRows >= 0
Âûõîäíûå ïàðàìåòðû:
PT - PTRows ïåðâûõ ñòîëáöîâ ìàòðèöû PT.
Ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..PTRows, 1..N]
Åñëè PTRows=0, òî ìàññèâ íå èçìåíÿåòñÿ.
-- ALGLIB --
Copyright 2005 by Bochkanov Sergey
*************************************************************************)
procedure UnpackPTFromBidiagonal(const QP: TReal2DArray;
M: Integer;
N: Integer;
const TauP: TReal1DArray;
PTRows: Integer;
var PT: TReal2DArray);
var
I: Integer;
J: Integer;
IP1: Integer;
V: TReal1DArray;
WORK: TReal1DArray;
VM: Integer;
i_: Integer;
i1_: Integer;
begin
Assert(PTRows<=N, 'UnpackPTFromBidiagonal: PTRows>N!');
if (M=0) or (N=0) or (PTRows=0) then
begin
Exit;
end;
//
// init
//
SetLength(PT, PTRows+1, N+1);
SetLength(V, N+1);
SetLength(WORK, PTRows+1);
//
// prepare PT
//
I:=1;
while I<=PTRows do
begin
J:=1;
while J<=N do
begin
if I=J then
begin
PT[I,J]:= 1;
end
else
begin
PT[I,J]:= 0;
end;
Inc(J);
end;
Inc(I);
end;
if M>=N then
begin
I:=Min(N-1, PTRows-1);
while I>=1 do
begin
VM:= N-I;
IP1:= I+1;
i1_:= (IP1) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[I,i_+i1_];
end;
V[1]:= 1;
ApplyReflectionFromTheRight(PT, TAUP[I], V, 1, PTRows, I+1, N, WORK);
Dec(I);
end;
end
else
begin
I:=Min(M, PTRows);
while I>=1 do
begin
VM:= N-I+1;
i1_:= (I) - (1);
for i_:= 1 to VM do
begin
V[i_]:= QP[I,i_+i1_];
end;
V[1]:= 1;
ApplyReflectionFromTheRight(PT, TAUP[I], V, 1, PTRows, I, N, WORK);
Dec(I);
end;
end;
end;
(*************************************************************************
Óìíîæåíèå íà ìàòðèöó P, ïðèâîäÿùóþ ìàòðèöó A ê äâóõäèàãîíàëüíîé ôîðìå.
Àëãîðèòì ïîçâîëÿåò îñóùåñòâèòü ïî âûáîðó óìíîæåíèå ñïðàâà èëè ñëåâà, íà
ìàòðèöó P èëè íà ìàòðèöó P'
Âõîäíûå ïàðàìåòðû:
QP - ìàòðèöû Q è P â óïàêîâàííîé ôîðìå.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
M - ÷èñëî ñòðîê â ìàòðèöå A
N - ÷èñëî ñòîëáöîâ â ìàòðèöå A
TAUP - ñêàëÿðíûå ìíîæèòåëè, ôîðìèðóþùèå ìàòðèöó P.
Ðåçóëüòàò ðàáîòû ïîäïðîãðàììû ToBidiagonal.
Z - óìíîæàåìàÿ ìàòðèöà
ìàññèâ ñ íóìåðàöèåé ýëåìåíòîâ [1..ZRows,1..ZColumns]
ZRows - ÷èñëî ñòðîê â ìàòðèöå Z. Åñëè FromTheRight=False, òî
ZRows=N, èíà÷å ìîæåò áûòü ëþáûì.
ZColumns - ÷èñëî ñòîëáöîâ â ìàòðèöå Z. Åñëè FromTheRight=True, òî
ZColumns=N, èíà÷å ìîæåò áû
