急求：矩阵QR分解的pascal源程序！！！！(200分)

越快越好，12号前给答案，再另给300分！

程序可以贴出来，或者email：OnlyNow@china.com 注明你的ID，以便给分。

C语言的源程序你要吗?

我已经自己解决了，100行而已～

这分收不回来吧？？？？
算了，那就送人好了：）

{ **********************************************************************
* Unit MATRICES.PAS *
* Version 2.0 *
* (c) J. Debord, May 2000 *
**********************************************************************
This unit implements dynamic allocation of vectors and matrices in
Pascal, together with various matrix operations.
Dynamic allocation is allowed by declaring arrays as pointers. There
are 8 types available :
PVector, PMatrix for floating point arrays
PIntVector, PIntMatrix for integer arrays
PBoolVector, PBoolMatrix for boolean arrays
PStrVector, PStrMatrix for string arrays (255 char.)
To use these arrays in your programs, you must :
(1) Declare variables of the appropriate type, e.g.
var
V : PVector;
A : PMatrix;
(2) Allocate each array BEFORE using it :
DimVector(V, N);
creates vector V[0..N]
DimMatrix(A, N, M);
creates matrix A[0..N, 0..M]
where N, M are two integer variables
If the allocation succeeds, all array elements are initialized
to zero (for numeric arrays), False (for boolean arrays), or
the null string (for string arrays). Otherwise, the pointer is
initialized to NIL.
(3) Use arrays as in standard Turbo Pascal, with the following
exceptions :
(a) You must use the indirection operator (^) to reference any
array element, i.e. write V^ and A^^[J] instead of
V and A[I,J].
(b) You cannot use the assignment operator =) to copy the
contents of an array into another array. Writing B := A
simply makes B point to the same memory block than A. You
must use one of the provided Copy... procedures (see their
do
cumentation in the interface part of the unit).
In addition, note that :
(a) All arrays begin
at index 0, so that the 0-indexed element
is always present, even if youdo
n't use it.
(b) A matrix is declared as an array of vectors, so that A^
denotes the I-th vector of matrix A and may be used as any
vector.
(4) Deallocate arrays when you no longer need them. This will free
the corresponding memory :
DelVector(V, N);
DelMatrix(A, N, M);
For more information, read the comments of each routine in the
interface part of the unit, and check the demo programs.
**********************************************************************
References :
1) 'Basic Programs for Scientists and Engineers' by A.R. Miller :
GaussJordan, InvMat
2) Borland's Numerical Methods Toolbox : Det
3) 'Numerical Recipes' by Press et al. : Cholesky, LU, SVD
4) 'Matrix Computations' by Golub &
Van Loan : QR_Decomp &
QR_Solve
(Pascal implementation contributed by Mark Vaughan)
********************************************************************** }
unit utMatrices;
interface
uses
FMath;
{ **********************************************************************
This section defines some error codes.
********************************************************************** }
const
MAT_OK = 0;
{ No error }
MAT_SINGUL = -1;
{ Singular matrix }
MAT_NON_CONV = -2;
{ Non convergence of iterative procedure }
MAT_NOT_PD = -3;
{ Matrix not positive definite }
{ **********************************************************************
This section defines the vector and matrix types. Maximal sizes are
given for a 16-bit compiler (TP / BP / Delphi 1). Higher values may
be used with the 32-bit compilers (Delphi 2-4, FPK, GPC).
********************************************************************** }
const
{$IFDEF EXTENDEDREAL}
MAX_FLT = 6552;
{ Max size of real vector }
{$else
}
{$IFDEF SINGLEREAL}
MAX_FLT = 16382;
{$else
}
{$IFDEF PASCALREAL}
MAX_FLT = 10921;
{$else
}
{$DEFINEdo
UBLEREAL}
MAX_FLT = 8190;
{$ENDIF}
{$ENDIF}
{$ENDIF}
MAX_INT = 16382;
{ Max size of integer vector }
MAX_BOOL = 32766;
{ Max size of boolean vector }
MAX_STR = 254;
{ Max size of string vector }
MAX_VEC = 16382;
{ Max number of vectors in a matrix }
TVector = array[0..MAX_FLT] of Float;
TIntVector = array[0..MAX_INT] of Integer;
TBoolVector = array[0..MAX_BOOL] of Boolean;
TStrVector = array[0..MAX_STR] of string;
PVector = ^TVector;
PIntVector = ^TIntVector;
PBoolVector = ^TBoolVector;
PStrVector = ^TStrVector;
TMatrix = array[0..MAX_VEC] of PVector;
TIntMatrix = array[0..MAX_VEC] of PIntVector;
TBoolMatrix = array[0..MAX_VEC] of PBoolVector;
TStrMatrix = array[0..MAX_VEC] of PStrVector;
PMatrix = ^TMatrix;
PIntMatrix = ^TIntMatrix;
PBoolMatrix = ^TBoolMatrix;
PStrMatrix = ^TStrMatrix;
{ **********************************************************************
Memory allocation routines
********************************************************************** }
procedure DimVector(var V: PVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Creates floating point vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DimIntVector(var V: PIntVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Creates integer vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DimBoolVector(var V: PBoolVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Creates boolean vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DimStrVector(var V: PStrVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Creates string vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DimMatrix(var A: PMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Creates floating point matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DimIntMatrix(var A: PIntMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Creates integer matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DimBoolMatrix(var A: PBoolMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Creates boolean matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DimStrMatrix(var A: PStrMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Creates string matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
{ **********************************************************************
Memory deallocation routines
********************************************************************** }
procedure DelVector(var V: PVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Deletes floating point vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DelIntVector(var V: PIntVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Deletes integer vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DelBoolVector(var V: PBoolVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Deletes boolean vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DelStrVector(var V: PStrVector;
Ubound: Integer);
{ ----------------------------------------------------------------------
Deletes string vector V[0..Ubound]
---------------------------------------------------------------------- }
procedure DelMatrix(var A: PMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Deletes floating point matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DelIntMatrix(var A: PIntMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Deletes integer matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DelBoolMatrix(var A: PBoolMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Deletes boolean matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
procedure DelStrMatrix(var A: PStrMatrix;
Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Deletes string matrix A[0..Ubound1, 0..Ubound2]
---------------------------------------------------------------------- }
{ **********************************************************************
Routines for copying vectors and matrices
----------------------------------------------------------------------
Lbound, Ubound : indices of first and last vector elements
Lbound1, Lbound2 : indices of first matrix element in each dimension
Ubound1, Ubound2 : indices of last matrix element in each dimension
********************************************************************** }
procedure SwapRows(I, K: Integer;
A: PMatrix;
Lbound, Ubound: Integer);
{ ----------------------------------------------------------------------
Exchanges rows I and K of matrix A
---------------------------------------------------------------------- }
procedure SwapCols(J, K: Integer;
A: PMatrix;
Lbound, Ubound: Integer);
{ ----------------------------------------------------------------------
Exchanges columns J and K of matrix A
---------------------------------------------------------------------- }
procedure CopyVector(Dest, Source: PVector;
Lbound, Ubound: Integer);
{ ----------------------------------------------------------------------
Copies vector Source into vector Dest
---------------------------------------------------------------------- }
procedure CopyMatrix(Dest, Source: PMatrix;
Lbound1, Lbound2, Ubound1, Ubound2: Integer);
{ ----------------------------------------------------------------------
Copies matrix Source into matrix Dest
---------------------------------------------------------------------- }
procedure CopyRowFromVector(Dest: PMatrix;
Source: PVector;
Lbound, Ubound, Row: Integer);
{ ----------------------------------------------------------------------
Copies vector Source into line Row of matrix Dest
---------------------------------------------------------------------- }
procedure CopyColFromVector(Dest: PMatrix;
Source: PVector;
Lbound, Ubound, Col: Integer);
{ ----------------------------------------------------------------------
Copies vector Source into column Col of matrix Dest
---------------------------------------------------------------------- }
procedure CopyVectorFromRow(Dest: PVector;
Source: PMatrix;
Lbound, Ubound, Row: Integer);
{ ----------------------------------------------------------------------
Copies line Row of matrix Source into vector Dest
---------------------------------------------------------------------- }
procedure CopyVectorFromCol(Dest: PVector;
Source: PMatrix;
Lbound, Ubound, Col: Integer);
{ ----------------------------------------------------------------------
Copies column Col of matrix Source into vector Dest
---------------------------------------------------------------------- }
{ **********************************************************************
Vector and matrix functions
********************************************************************** }
function Min(X: PVector;
Lbound, Ubound: Integer): Float;
{ ----------------------------------------------------------------------
Returns the lowest value of vector X
---------------------------------------------------------------------- }
function Max(X: PVector;
Lbound, Ubound: Integer): Float;
{ ----------------------------------------------------------------------
Returns the highest value of vector X
---------------------------------------------------------------------- }
function IntMin(X: PIntVector;
Lbound, Ubound: Integer): Integer;
{ ----------------------------------------------------------------------
Returns the lowest value of integer vector X
---------------------------------------------------------------------- }
function IntMax(X: PIntVector;
Lbound, Ubound: Integer): Integer;
{ ----------------------------------------------------------------------
Returns the highest value of integer vector X
---------------------------------------------------------------------- }
procedure Transpose(A: PMatrix;
Lbound1, Lbound2, Ubound1, Ubound2: Integer;
A_t: PMatrix);
{ ----------------------------------------------------------------------
Transposes a matrix
----------------------------------------------------------------------
Input parameters : A = original matrix
Lbound1,
Lbound2 = indices of 1st matrix elem. in each dim.
Ubound1,
Ubound2 = indices of last matrix elem. in each dim.
----------------------------------------------------------------------
Output parameter : A_t = transposed matrix
---------------------------------------------------------------------- }
function GaussJordan(A: PMatrix;
B: PVector;
Lbound, Ubound: Integer;
A_inv: PMatrix;
X: PVector): Integer;
{ ----------------------------------------------------------------------
Solves a system of linear equations by the Gauss-Jordan method
----------------------------------------------------------------------
Input parameters : A = system matrix
B = constant vector
Lbound = index of first matrix element
Ubound = index of last matrix element
----------------------------------------------------------------------
Output parameters : A_inv = inverse matrix
X = solution vector
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_SINGUL
---------------------------------------------------------------------- }
function InvMat(A: PMatrix;
Lbound, Ubound: Integer;
A_inv: PMatrix): Integer;
{ ----------------------------------------------------------------------
Computes the inverse of a square matrix by the Gauss-Jordan method
----------------------------------------------------------------------
Parameters : as in Gauss-Jordan
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_SINGUL
---------------------------------------------------------------------- }
function Det(A: PMatrix;
Lbound, Ubound: Integer): Float;
{ ----------------------------------------------------------------------
Computes the determinant of a square matrix
----------------------------------------------------------------------
Parameters : as in Gauss-Jordan
----------------------------------------------------------------------
NB : This procedure destroys the original matrix A
---------------------------------------------------------------------- }
function Cholesky(A: PMatrix;
Lbound, Ubound: Integer;
L: PMatrix): Integer;
{ ----------------------------------------------------------------------
Cholesky decomposition. Factors the symmetric positive definite matrix
A as a product L * L', where L is a lower triangular matrix. This
procedure may be used as a test of positive definiteness.
----------------------------------------------------------------------
Input parameters : A = matrix
Lbound = index of first matrix element
Ubound = index of last matrix element
----------------------------------------------------------------------
Output parameter : L = Cholesky factor of matrix A
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_NOT_PD
---------------------------------------------------------------------- }
function LU_Decomp(A: PMatrix;
Lbound, Ubound: Integer): Integer;
{ ----------------------------------------------------------------------
LU decomposition. Factors the square matrix A as a product L * U,
where L is a lower triangular matrix (with unit diagonal terms) and U
is an upper triangular matrix. This routine is used in conjunction
with LU_Solve to solve a system of equations.
----------------------------------------------------------------------
Input parameters : A = matrix
Lbound = index of first matrix element
Ubound = index of last matrix element
----------------------------------------------------------------------
Output parameter : A = contains the elements of L and U
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_SINGUL
----------------------------------------------------------------------
NB : This procedure destroys the original matrix A
---------------------------------------------------------------------- }
procedure LU_Solve(A: PMatrix;
B: PVector;
Lbound, Ubound: Integer;
X: PVector);
{ ----------------------------------------------------------------------
Solves a system of equations whose matrix has been transformed by
LU_Decomp
----------------------------------------------------------------------
Input parameters : A = result from LU_Decomp
B = constant vector
Lbound,
Ubound = as in LU_Decomp
----------------------------------------------------------------------
Output parameter : X = solution vector
---------------------------------------------------------------------- }
function SV_Decomp(A: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
S: PVector;
V: PMatrix): Integer;
{ ----------------------------------------------------------------------
Singular value decomposition. Factors the matrix A (n x m, with n >= m)
as a product U * S * V' where U is a (n x m) column-orthogonal matrix,
S a (m x m) diagonal matrix with elements >= 0 (the singular values)
and V a (m x m) orthogonal matrix. This routine is used in conjunction
with SV_Solve to solve a system of equations.
----------------------------------------------------------------------
Input parameters : A = matrix
Lbound = index of first matrix element
Ubound1 = index of last matrix element in 1st dim.
Ubound2 = index of last matrix element in 2nd dim.
----------------------------------------------------------------------
Output parameter : A = contains the elements of U
S = vector of singular values
V = orthogonal matrix
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_NON_CONV
----------------------------------------------------------------------
NB : This procedure destroys the original matrix A
---------------------------------------------------------------------- }
procedure SV_SetZero(S: PVector;
Lbound, Ubound: Integer;
Tol: Float);
{ ----------------------------------------------------------------------
Sets the singular values to zero if they are lower than a specified
threshold.
----------------------------------------------------------------------
Input parameters : S = vector of singular values
Tol = relative tolerance
Threshold value will be Tol * Max(S)
Lbound = index of first vector element
Ubound = index of last vector element
----------------------------------------------------------------------
Output parameter : S = modified singular values
---------------------------------------------------------------------- }
procedure SV_Solve(U: PMatrix;
S: PVector;
V: PMatrix;
B: PVector;
Lbound, Ubound1, Ubound2: Integer;
X: PVector);
{ ----------------------------------------------------------------------
Solves a system of equations by singular value decomposition, after
the matrix has been transformed by SV_Decomp, and the lowest singular
values have been set to zero by SV_SetZero.
----------------------------------------------------------------------
Input parameters : U, S, V = vector and matrices from SV_Decomp
B = constant vector
Lbound,
Ubound1,
Ubound2 = as in SV_Decomp
----------------------------------------------------------------------
Output parameter : X = solution vector
= V * Diag(1/s(i)) * U' * B, for s(i) <> 0
---------------------------------------------------------------------- }
procedure SV_Approx(U: PMatrix;
S: PVector;
V: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
A: PMatrix);
{ ----------------------------------------------------------------------
Approximates a matrix A by the product USV', after the lowest singular
values have been set to zero by SV_SetZero.
----------------------------------------------------------------------
Input parameters : U, S, V = vector and matrices from SV_Decomp
Lbound,
Ubound1,
Ubound2 = as in SV_Decomp
----------------------------------------------------------------------
Output parameter : A = approximated matrix
---------------------------------------------------------------------- }
function QR_Decomp(A: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
R: PMatrix): Integer;
{ ----------------------------------------------------------------------
QR decomposition. Factors the matrix A (n x m, with n >= m) as a
product Q * R where Q is a (n x m) column-orthogonal matrix, and R
a (m x m) upper triangular matrix. This routine is used in conjunction
with QR_Solve to solve a system of equations.
----------------------------------------------------------------------
Input parameters : A = matrix
Lbound = index of first matrix element
Ubound1 = index of last matrix element in 1st dim.
Ubound2 = index of last matrix element in 2nd dim.
----------------------------------------------------------------------
Output parameter : A = contains the elements of Q
R = upper triangular matrix
----------------------------------------------------------------------
Possible results : MAT_OK
MAT_SING
----------------------------------------------------------------------
NB : This procedure destroys the original matrix A
---------------------------------------------------------------------- }
procedure QR_Solve(Q, R: PMatrix;
B: PVector;
Lbound, Ubound1, Ubound2: Integer;
X: PVector);
{ ----------------------------------------------------------------------
Solves a system of equations by the QR decomposition,
after the matrix has been transformed by QR_Decomp.
----------------------------------------------------------------------
Input parameters : Q, R = matrices from QR_Decomp
B = constant vector
Lbound,
Ubound1,
Ubound2 = as in QR_Decomp
----------------------------------------------------------------------
Output parameter : X = solution vector
---------------------------------------------------------------------- }
implementation
const
{ Used by LU procedures }
LastDim: Integer = 1;
{ Dimension of the last system solved }
Index: PIntVector = nil;
{ Records the row permutations }
procedure DimVector(var V: PVector;
Ubound: Integer);
var
I: Integer;
begin
{ Check bounds }
if (Ubound < 0) or (Ubound > MAX_FLT) then
begin
V := nil;
Exit;
end;

{ Allocate vector }
GetMem(V, Succ(Ubound) * SizeOf(Float));
if V = nil then
Exit;
{ Initialize vector }
for I := 0 to Ubounddo
V^ := 0.0;
end;

procedure DimIntVector(var V: PIntVector;
Ubound: Integer);
var
I: Integer;
begin
{ Check bounds }
if (Ubound < 0) or (Ubound > MAX_INT) then
begin
V := nil;
Exit;
end;

{ Allocate vector }
GetMem(V, Succ(Ubound) * SizeOf(Integer));
if V <> nil then
for I := 0 to Ubounddo
V^ := 0;
end;

procedure DimBoolVector(var V: PBoolVector;
Ubound: Integer);
var
I: Integer;
begin
{ Check bounds }
if (Ubound < 0) or (Ubound > MAX_BOOL) then
begin
V := nil;
Exit;
end;

{ Allocate vector }
GetMem(V, Succ(Ubound) * SizeOf(Boolean));
if V <> nil then
Exit;
{ Initialize vector }
for I := 0 to Ubounddo
V^ := False;
end;

procedure DimStrVector(var V: PStrVector;
Ubound: Integer);
var
I: Integer;
begin
{ Check bounds }
if (Ubound < 0) or (Ubound > MAX_STR) then
begin
V := nil;
Exit;
end;

{ Allocate vector }
GetMem(V, Succ(Ubound) * 256);
if V <> nil then
Exit;
{ Initialize vector }
for I := 0 to Ubounddo
V^ := '';
end;

procedure DimMatrix(var A: PMatrix;
Ubound1, Ubound2: Integer);
var
I, J: Integer;
RowSize: Word;
begin
if (Ubound1 < 0) or (Ubound1 > MAX_VEC) or
(Ubound2 < 0) or (Ubound2 > MAX_FLT) then
begin
A := nil;
Exit;
end;

{ Allocate matrix }
GetMem(A, Succ(Ubound1) * SizeOf(PVector));
if A = nil then
Exit;
{ Size of a row }
RowSize := Succ(Ubound2) * SizeOf(Float);
{ Allocate each row }
for I := 0 to Ubound1do
begin
GetMem(A^, RowSize);
if A^ = nil then
begin
A := nil;
Exit;
end;
end;

{ Initialize matrix }
for I := 0 to Ubound1do
for J := 0 to Ubound2do
A^^[J] := 0.0;
end;

procedure DimIntMatrix(var A: PIntMatrix;
Ubound1, Ubound2: Integer);
var
I, J: Integer;
RowSize: Word;
begin
{ Check bounds }
if (Ubound1 < 0) or (Ubound1 > MAX_VEC) or
(Ubound2 < 0) or (Ubound2 > MAX_INT) then
begin
A := nil;
Exit;
end;

{ Allocate matrix }
GetMem(A, Succ(Ubound1) * SizeOf(PIntVector));
if A = nil then
Exit;
{ Size of a row }
RowSize := Succ(Ubound2) * SizeOf(Integer);
{ Allocate each row }
for I := 0 to Ubound1do
begin
GetMem(A^, RowSize);
if A^ = nil then
begin
A := nil;
Exit;
end;
end;

{ Initialize matrix }
for I := 0 to Ubound1do
for J := 0 to Ubound2do
A^^[J] := 0;
end;

procedure DimBoolMatrix(var A: PBoolMatrix;
Ubound1, Ubound2: Integer);
var
I, J: Integer;
RowSize: Word;
begin
{ Check bounds }
if (Ubound1 < 0) or (Ubound1 > MAX_VEC) or
(Ubound2 < 0) or (Ubound2 > MAX_BOOL) then
begin
A := nil;
Exit;
end;

{ Allocate matrix }
GetMem(A, Succ(Ubound1) * SizeOf(PBoolVector));
if A = nil then
Exit;
{ Size of a row }
RowSize := Succ(Ubound2) * SizeOf(Boolean);
{ Allocate each row }
for I := 0 to Ubound1do
begin
GetMem(A^, RowSize);
if A^ = nil then
begin
A := nil;
Exit;
end;
end;

{ Initialize matrix }
for I := 0 to Ubound1do
for J := 0 to Ubound2do
A^^[J] := False;
end;

procedure DimStrMatrix(var A: PStrMatrix;
Ubound1, Ubound2: Integer);
var
I, J: Integer;
RowSize: Word;
begin
{ Check bounds }
if (Ubound1 < 0) or (Ubound1 > MAX_VEC) or
(Ubound2 < 0) or (Ubound2 > MAX_STR) then
begin
A := nil;
Exit;
end;

{ Allocate matrix }
GetMem(A, Succ(Ubound1) * SizeOf(PStrVector));
if A = nil then
Exit;
{ Size of a row }
RowSize := Succ(Ubound2) * 256;
{ Allocate each row }
for I := 0 to Ubound1do
begin
GetMem(A^, RowSize);
if A^ = nil then
begin
A := nil;
Exit;
end;
end;

{ Initialize matrix }
for I := 0 to Ubound1do
for J := 0 to Ubound2do
A^^[J] := '';
end;

procedure DelVector(var V: PVector;
Ubound: Integer);
begin
if V <> nil then
begin
FreeMem(V, Succ(Ubound) * SizeOf(Float));
V := nil;
end;
end;

procedure DelIntVector(var V: PIntVector;
Ubound: Integer);
begin
if V <> nil then
begin
FreeMem(V, Succ(Ubound) * SizeOf(Integer));
V := nil;
end;
end;

procedure DelBoolVector(var V: PBoolVector;
Ubound: Integer);
begin
if V <> nil then
begin
FreeMem(V, Succ(Ubound) * SizeOf(Boolean));
V := nil;
end;
end;

procedure DelStrVector(var V: PStrVector;
Ubound: Integer);
begin
if V <> nil then
begin
FreeMem(V, Succ(Ubound) * 256);
V := nil;
end;
end;

procedure DelMatrix(var A: PMatrix;
Ubound1, Ubound2: Integer);
var
I: Integer;
RowSize: Word;
begin
if A <> nil then
begin
RowSize := Succ(Ubound2) * SizeOf(Float);
for I := Ubound1do
wnto 0do
FreeMem(A^, RowSize);
FreeMem(A, Succ(Ubound1) * SizeOf(PVector));
A := nil;
end;
end;

procedure DelIntMatrix(var A: PIntMatrix;
Ubound1, Ubound2: Integer);
var
I: Integer;
RowSize: Word;
begin
if A <> nil then
begin
RowSize := Succ(Ubound2) * SizeOf(Integer);
for I := Ubound1do
wnto 0do
FreeMem(A^, RowSize);
FreeMem(A, Succ(Ubound1) * SizeOf(PIntVector));
A := nil;
end;
end;

procedure DelBoolMatrix(var A: PBoolMatrix;
Ubound1, Ubound2: Integer);
var
I: Integer;
RowSize: Word;
begin
if A <> nil then
begin
RowSize := Succ(Ubound2) * SizeOf(Boolean);
for I := Ubound1do
wnto 0do
FreeMem(A^, RowSize);
FreeMem(A, Succ(Ubound1) * SizeOf(PBoolVector));
A := nil;
end;
end;

procedure DelStrMatrix(var A: PStrMatrix;
Ubound1, Ubound2: Integer);
var
I: Integer;
RowSize: Word;
begin
if A <> nil then
begin
RowSize := Succ(Ubound2) * 256;
for I := Ubound1do
wnto 0do
FreeMem(A^, RowSize);
FreeMem(A, Succ(Ubound1) * SizeOf(PStrVector));
A := nil;
end;
end;

procedure SwapRows(I, K: Integer;
A: PMatrix;
Lbound, Ubound: Integer);
var
J: Integer;
begin
for J := Lbound to Ubounddo
FSwap(A^^[J], A^[K]^[J]);
end;

procedure SwapCols(J, K: Integer;
A: PMatrix;
Lbound, Ubound: Integer);
var
I: Integer;
begin
for I := Lbound to Ubounddo
FSwap(A^^[J], A^^[K]);
end;

procedure CopyVector(Dest, Source: PVector;
Lbound, Ubound: Integer);
var
I: Integer;
begin
for I := Lbound to Ubounddo
Dest^ := Source^;
end;

procedure CopyMatrix(Dest, Source: PMatrix;
Lbound1, Lbound2, Ubound1, Ubound2: Integer);
var
I, J: Integer;
begin
for I := Lbound1 to Ubound1do
for J := Lbound2 to Ubound2do
Dest^^[J] := Source^^[J];
end;

procedure CopyRowFromVector(Dest: PMatrix;
Source: PVector;
Lbound, Ubound, Row: Integer);
var
J: Integer;
begin
for J := Lbound to Ubounddo
Dest^[Row]^[J] := Source^[J];
end;

procedure CopyColFromVector(Dest: PMatrix;
Source: PVector;
Lbound, Ubound, Col: Integer);
var
I: Integer;
begin
for I := Lbound to Ubounddo
Dest^^[Col] := Source^;
end;

procedure CopyVectorFromRow(Dest: PVector;
Source: PMatrix;
Lbound, Ubound, Row: Integer);
var
J: Integer;
begin
for J := Lbound to Ubounddo
Dest^[J] := Source^[Row]^[J];
end;

procedure CopyVectorFromCol(Dest: PVector;
Source: PMatrix;
Lbound, Ubound, Col: Integer);
var
I: Integer;
begin
for I := Lbound to Ubounddo
Dest^ := Source^^[Col];
end;

function Min(X: PVector;
Lbound, Ubound: Integer): Float;
var
Xmin: Float;
I: Integer;
begin
Xmin := X^[Lbound];
for I := Succ(Lbound) to Ubounddo
if X^ < Xmin then
Xmin := X^;
Min := Xmin;
end;

function Max(X: PVector;
Lbound, Ubound: Integer): Float;
var
Xmax: Float;
I: Integer;
begin
Xmax := X^[Lbound];
for I := Succ(Lbound) to Ubounddo
if X^ > Xmax then
Xmax := X^;
Max := Xmax;
end;

function IntMin(X: PIntVector;
Lbound, Ubound: Integer): Integer;
var
I, Xmin: Integer;
begin
Xmin := X^[Lbound];
for I := Succ(Lbound) to Ubounddo
if X^ < Xmin then
Xmin := X^;
IntMin := Xmin;
end;

function IntMax(X: PIntVector;
Lbound, Ubound: Integer): Integer;
var
I, Xmax: Integer;
begin
Xmax := X^[Lbound];
for I := Succ(Lbound) to Ubounddo
if X^ > Xmax then
Xmax := X^;
IntMax := Xmax;
end;

procedure Transpose(A: PMatrix;
Lbound1, Lbound2, Ubound1, Ubound2: Integer;
A_t: PMatrix);
var
I, J: Integer;
begin
for I := Lbound1 to Ubound1do
for J := Lbound2 to Ubound2do
A_t^[J]^ := A^^[J];
end;

function GaussJordan(A: PMatrix;
B: PVector;
Lbound, Ubound: Integer;
A_inv: PMatrix;
X: PVector): Integer;
var
I, J, K: Integer;
Pvt, T: Float;
PRow, PCol: PIntVector;
{ Store line and column of pivot }
begin
DimIntVector(PRow, Ubound);
DimIntVector(PCol, Ubound);
{ Copy A into A_inv and B into X }
CopyMatrix(A_inv, A, Lbound, Lbound, Ubound, Ubound);
CopyVector(X, B, Lbound, Ubound);
K := Lbound;
while K <= Ubounddo
begin
{ Search for largest pivot in submatrix A_inv[K..Ubound, K..Ubound] }
Pvt := A_inv^[K]^[K];
PRow^[K] := K;
PCol^[K] := K;
for I := K to Ubounddo
for J := K to Ubounddo
if Abs(A_inv^^[J]) > Abs(Pvt) then
begin
Pvt := A_inv^^[J];
PRow^[K] := I;
PCol^[K] := J;
end;

{ Pivot too weak ==> quasi-singular matrix }
if Abs(Pvt) < MACHEP then
begin
DelIntVector(PRow, Ubound);
DelIntVector(PCol, Ubound);
GaussJordan := MAT_SINGUL;
Exit;
end;

{ Exchange current row (K) with pivot row }
if PRow^[K] <> K then
begin
SwapRows(PRow^[K], K, A_inv, Lbound, Ubound);
FSwap(X^[PRow^[K]], X^[K]);
end;

{ Exchange current column (K) with pivot column }
if PCol^[K] <> K then
SwapCols(PCol^[K], K, A_inv, Lbound, Ubound);
{ Transform pivot row }
A_inv^[K]^[K] := 1.0;
for J := Lbound to Ubounddo
A_inv^[K]^[J] := A_inv^[K]^[J] / Pvt;
X^[K] := X^[K] / Pvt;
{ Transform other rows }
for I := Lbound to Ubounddo
if I <> K then
begin
T := A_inv^^[K];
A_inv^^[K] := 0.0;
for J := Lbound to Ubounddo
A_inv^^[J] := A_inv^^[J] - T * A_inv^[K]^[J];
X^ := X^ - T * X^[K];
end;
Inc(K);
end;

{ Rearrange inverse matrix }
for I := Ubounddo
wnto Lbounddo
if PCol^ <> I then
begin
SwapRows(PCol^, I, A_inv, Lbound, Ubound);
FSwap(X^[PCol^], X^);
end;
for J := Ubounddo
wnto Lbounddo
if PRow^[J] <> J then
SwapCols(PRow^[J], J, A_inv, Lbound, Ubound);
DelIntVector(PRow, Ubound);
DelIntVector(PCol, Ubound);
GaussJordan := MAT_OK;
end;

function InvMat(A: PMatrix;
Lbound, Ubound: Integer;
A_inv: PMatrix): Integer;
var
I, J, K: Integer;
Pvt, T: Float;
PRow, PCol: PIntVector;
{ Store line and column of pivot }
begin
DimIntVector(PRow, Ubound);
DimIntVector(PCol, Ubound);
{ Copy A into A_inv }
CopyMatrix(A_inv, A, Lbound, Lbound, Ubound, Ubound);
K := Lbound;
while K <= Ubounddo
begin
{ Search for largest pivot in submatrix A_inv[K..Ubound, K..Ubound] }
Pvt := A_inv^[K]^[K];
PRow^[K] := K;
PCol^[K] := K;
for I := K to Ubounddo
for J := K to Ubounddo
if Abs(A_inv^^[J]) > Abs(Pvt) then
begin
Pvt := A_inv^^[J];
PRow^[K] := I;
PCol^[K] := J;
end;

{ Pivot too weak ==> quasi-singular matrix }
if Abs(Pvt) < MACHEP then
begin
DelIntVector(PRow, Ubound);
DelIntVector(PCol, Ubound);
InvMat := MAT_SINGUL;
Exit;
end;

{ Exchange current row (K) with pivot row }
if PRow^[K] <> K then
SwapRows(PRow^[K], K, A_inv, Lbound, Ubound);
{ Exchange current column (K) with pivot column }
if PCol^[K] <> K then
SwapCols(PCol^[K], K, A_inv, Lbound, Ubound);
{ Transform pivot row }
A_inv^[K]^[K] := 1.0;
for J := Lbound to Ubounddo
A_inv^[K]^[J] := A_inv^[K]^[J] / Pvt;
{ Transform other rows }
for I := Lbound to Ubounddo
if I <> K then
begin
T := A_inv^^[K];
A_inv^^[K] := 0.0;
for J := Lbound to Ubounddo
A_inv^^[J] := A_inv^^[J] - T * A_inv^[K]^[J];
end;
Inc(K);
end;

{ Rearrange inverse matrix }
for I := Ubounddo
wnto Lbounddo
if PCol^ <> I then
SwapRows(PCol^, I, A_inv, Lbound, Ubound);
for J := Ubounddo
wnto Lbounddo
if PRow^[J] <> J then
SwapCols(PRow^[J], J, A_inv, Lbound, Ubound);
DelIntVector(PRow, Ubound);
DelIntVector(PCol, Ubound);
InvMat := MAT_OK;
end;

function Det(A: PMatrix;
Lbound, Ubound: Integer): Float;
var
D, T: Float;
{ Partial determinant &amp;
multiplier }
I, J, K: Integer;
{ Loop variables }
ZeroDet: Boolean;
{ Flags a null determinant }
begin
ZeroDet := False;
D := 1.0;
K := Lbound;
{ Make the matrix upper triangular }
while not (ZeroDet) and (K < Ubound)do
begin
{ If diagonal element is zero then
switch rows }
if Abs(A^[K]^[K]) < MACHEP then
begin
ZeroDet := True;
I := K;
{ Try to find a row with a non-zero element in this column }
while ZeroDet and (I < Ubound)do
begin
I := Succ(I);
if Abs(A^^[K]) > MACHEP then
begin
{ Switch these two rows }
SwapRows(I, K, A, Lbound, Ubound);
ZeroDet := False;
{ Switching rows changes the sign of the determinant }
D := -D;
end;
end;
end;

if not (ZeroDet) then
for I := Succ(K) to Ubounddo
if Abs(A^^[K]) > MACHEP then
begin
{ Make the K element of this row zero }
T := -A^^[K] / A^[K]^[K];
for J := 1 to Ubounddo
A^^[J] := A^^[J] + T * A^[K]^[J];
end;

D := D * A^[K]^[K];
{ Multiply the diagonal term into D }
Inc(K);
end;

if ZeroDet then
Det := 0.0
else
Det := D * A^[Ubound]^[Ubound];
end;

function Cholesky(A: PMatrix;
Lbound, Ubound: Integer;
L: PMatrix): Integer;
var
I, J, K: Integer;
Sum: Float;
begin
for K := Lbound to Ubounddo
begin
Sum := A^[K]^[K];
for J := Lbound to K - 1do
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 := K + 1 to Ubounddo
begin
Sum := A^^[K];
for J := Lbound to K - 1do
Sum := Sum - L^^[J] * L^[K]^[J];
L^^[K] := Sum / L^[K]^[K];
end;
end;
Cholesky := MAT_OK;
end;

function LU_Decomp(A: PMatrix;
Lbound, Ubound: Integer): Integer;
const
TINY = 1.0E-20;
var
I, Imax, J, K: Integer;
Pvt, T, Sum: Float;
V: PVector;
begin
DimVector(V, Ubound);
{ Reallocate Index }
if Index <> nil then
DelIntVector(Index, LastDim);
DimIntVector(Index, Ubound);
LastDim := Ubound;
for I := Lbound to Ubounddo
begin
Pvt := 0.0;
for J := Lbound to Ubounddo
if Abs(A^^[J]) > Pvt then
Pvt := Abs(A^^[J]);
if Pvt < MACHEP then
begin
DelVector(V, Ubound);
LU_Decomp := MAT_SINGUL;
Exit;
end;
V^ := 1.0 / Pvt;
end;
for J := Lbound to Ubounddo
begin
for I := Lbound to Pred(J)do
begin
Sum := A^^[J];
for K := Lbound to Pred(I)do
Sum := Sum - A^^[K] * A^[K]^[J];
A^^[J] := Sum;
end;
Pvt := 0.0;
for I := J to Ubounddo
begin
Sum := A^^[J];
for K := Lbound to Pred(J)do
Sum := Sum - A^^[K] * A^[K]^[J];
A^^[J] := Sum;
T := V^ * 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] := TINY;
if J <> Ubound then
begin
T := 1.0 / A^[J]^[J];
for I := Succ(J) to Ubounddo
A^^[J] := A^^[J] * T;
end;
end;
DelVector(V, Ubound);
LU_Decomp := MAT_OK;
end;

procedure LU_Solve(A: PMatrix;
B: PVector;
Lbound, Ubound: Integer;
X: PVector);
var
I, Ip, J, K: Integer;
Sum: Float;
begin
K := Pred(Lbound);
CopyVector(X, B, Lbound, Ubound);
for I := Lbound to Ubounddo
begin
Ip := Index^;
Sum := X^[Ip];
X^[Ip] := X^;
if K >= Lbound then
for J := K to Pred(I)do
Sum := Sum - A^^[J] * X^[J]
else
if Sum <> 0.0 then
K := I;
X^ := Sum;
end;
for I := Ubounddo
wnto Lbounddo
begin
Sum := X^;
if I < Ubound then
for J := Succ(I) to Ubounddo
Sum := Sum - A^^[J] * X^[J];
X^ := Sum / A^^;
end;
end;

function SV_Decomp(A: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
S: PVector;
V: PMatrix): Integer;
label
1, 2, 3;
var
I, Its, J, JJ, K, L, N: Integer;
Anorm, C, F, G, H, Sum, Scale, T, X, Y, Z: Float;
R: PVector;
begin
G := 0.0;
Scale := 0.0;
Anorm := 0.0;
DimVector(R, Ubound2);
for I := Lbound to Ubound2do
begin
L := I + 1;
R^ := Scale * G;
G := 0.0;
Sum := 0.0;
Scale := 0.0;
if I <= Ubound1 then
begin
for K := I to Ubound1do
Scale := Scale + Abs(A^[K]^);
if Scale <> 0.0 then
begin
for K := I to Ubound1do
begin
A^[K]^ := A^[K]^ / Scale;
Sum := Sum + A^[K]^ * A^[K]^;
end;
F := A^^;
G := -Sgn(F) * Sqrt(Sum);
H := F * G - Sum;
A^^ := F - G;
if I <> Ubound2 then
begin
for J := L to Ubound2do
begin
Sum := 0.0;
for K := I to Ubound1do
Sum := Sum + A^[K]^ * A^[K]^[J];
F := Sum / H;
for K := I to Ubound1do
A^[K]^[J] := A^[K]^[J] + F * A^[K]^;
end;
end;
for K := I to Ubound1do
A^[K]^ := Scale * A^[K]^;
end;
end;
S^ := Scale * G;
G := 0.0;
Sum := 0.0;
Scale := 0.0;
if (I <= Ubound1) and (I <> Ubound2) then
begin
for K := L to Ubound2do
Scale := Scale + Abs(A^^[K]);
if Scale <> 0.0 then
begin
for K := L to Ubound2do
begin
A^^[K] := A^^[K] / Scale;
Sum := Sum + A^^[K] * A^^[K];
end;
F := A^^[L];
G := -Sgn(F) * Sqrt(Sum);
H := F * G - Sum;
A^^[L] := F - G;
for K := L to Ubound2do
R^[K] := A^^[K] / H;
if I <> Ubound1 then
for J := L to Ubound1do
begin
Sum := 0.0;
for K := L to Ubound2do
Sum := Sum + A^[J]^[K] * A^^[K];
for K := L to Ubound2do
A^[J]^[K] := A^[J]^[K] + Sum * R^[K];
end;
for K := L to Ubound2do
A^^[K] := Scale * A^^[K];
end;
end;
Anorm := FMax(Anorm, Abs(S^) + Abs(R^));
end;
for I := Ubound2do
wnto Lbounddo
begin
if I < Ubound2 then
begin
if G <> 0.0 then
begin
for J := L to Ubound2do
V^[J]^ := (A^^[J] / A^^[L]) / G;
for J := L to Ubound2do
begin
Sum := 0.0;
for K := L to Ubound2do
Sum := Sum + A^^[K] * V^[K]^[J];
for K := L to Ubound2do
V^[K]^[J] := V^[K]^[J] + Sum * V^[K]^;
end;
end;
for J := L to Ubound2do
begin
V^^[J] := 0.0;
V^[J]^ := 0.0;
end;
end;
V^^ := 1.0;
G := R^;
L := I;
end;
for I := Ubound2do
wnto Lbounddo
begin
L := I + 1;
G := S^;
if I < Ubound2 then
for J := L to Ubound2do
A^^[J] := 0.0;
if G <> 0.0 then
begin
G := 1.0 / G;
if I <> Ubound2 then
for J := L to Ubound2do
begin
Sum := 0.0;
for K := L to Ubound1do
Sum := Sum + A^[K]^ * A^[K]^[J];
F := (Sum / A^^) * G;
for K := I to Ubound1do
A^[K]^[J] := A^[K]^[J] + F * A^[K]^;
end;
for J := I to Ubound1do
A^[J]^ := A^[J]^ * G;
end
else
for J := I to Ubound1do
A^[J]^ := 0.0;
A^^ := A^^ + 1.0;
end;
for K := Ubound2do
wnto Lbounddo
begin
for Its := 1 to 30do
begin
for L := Kdo
wnto Lbounddo
begin
N := L - 1;
if (Abs(R^[L]) + Anorm) = Anorm then
goto 2;
if (Abs(S^[N]) + Anorm) = Anorm then
goto 1;
end;
1: T := 1.0;
for I := L to Kdo
begin
F := T * R^;
if (Abs(F) + Anorm) <> Anorm then
begin
G := S^;
H := Pythag(F, G);
S^ := H;
H := 1.0 / H;
C := G * H;
T := -(F * H);
for J := Lbound to Ubound1do
begin
Y := A^[J]^[N];
Z := A^[J]^;
A^[J]^[N] := (Y * C) + (Z * T);
A^[J]^ := -(Y * T) + (Z * C);
end;
end;
end;
2: Z := S^[K];
if L = K then
begin
if Z < 0.0 then
begin
S^[K] := -Z;
for J := Lbound to Ubound2do
V^[J]^[K] := -V^[J]^[K];
end;
goto 3
end;
if Its = 30 then
begin
DelVector(R, Ubound2);
SV_Decomp := MAT_NON_CONV;
Exit;
end;
X := S^[L];
N := K - 1;
Y := S^[N];
G := R^[N];
H := R^[K];
F := ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2.0 * H * Y);
G := Pythag(F, 1.0);
F := ((X - Z) * (X + Z) + H * ((Y / (F + Sgn(F) * Abs(G))) - H)) / X;
C := 1.0;
T := 1.0;
for J := L to Ndo
begin
I := J + 1;
G := R^;
Y := S^;
H := T * G;
G := C * G;
Z := Pythag(F, H);
R^[J] := Z;
C := F / Z;
T := H / Z;
F := (X * C) + (G * T);
G := -(X * T) + (G * C);
H := Y * T;
Y := Y * C;
for JJ := Lbound to Ubound2do
begin
X := V^[JJ]^[J];
Z := V^[JJ]^;
V^[JJ]^[J] := (X * C) + (Z * T);
V^[JJ]^ := -(X * T) + (Z * C);
end;
Z := Pythag(F, H);
S^[J] := Z;
if Z <> 0.0 then
begin
Z := 1.0 / Z;
C := F * Z;
T := H * Z;
end;
F := (C * G) + (T * Y);
X := -(T * G) + (C * Y);
for JJ := Lbound to Ubound1do
begin
Y := A^[JJ]^[J];
Z := A^[JJ]^;
A^[JJ]^[J] := (Y * C) + (Z * T);
A^[JJ]^ := -(Y * T) + (Z * C);
end
end;
R^[L] := 0.0;
R^[K] := F;
S^[K] := X;
end;
3:
end;
DelVector(R, Ubound2);
SV_Decomp := MAT_OK;
end;

procedure SV_SetZero(S: PVector;
Lbound, Ubound: Integer;
Tol: Float);
var
Threshold: Float;
I: Integer;
begin
Threshold := Tol * Max(S, Lbound, Ubound);
for I := Lbound to Ubounddo
if S^ < Threshold then
S^ := 0.0;
end;

procedure SV_Solve(U: PMatrix;
S: PVector;
V: PMatrix;
B: PVector;
Lbound, Ubound1, Ubound2: Integer;
X: PVector);
var
I, J, JJ: Integer;
Sum: Float;
Tmp: PVector;
begin
DimVector(Tmp, Ubound2);
for J := Lbound to Ubound2do
begin
Sum := 0.0;
if S^[J] > 0.0 then
begin
for I := Lbound to Ubound1do
Sum := Sum + U^^[J] * B^;
Sum := Sum / S^[J];
end;
Tmp^[J] := Sum;
end;
for J := Lbound to Ubound2do
begin
Sum := 0.0;
for JJ := Lbound to Ubound2do
Sum := Sum + V^[J]^[JJ] * Tmp^[JJ];
X^[J] := Sum;
end;
DelVector(Tmp, Ubound2);
end;

procedure SV_Approx(U: PMatrix;
S: PVector;
V: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
A: PMatrix);
var
I, J, K: Integer;
begin
for I := Lbound to Ubound1do
for J := Lbound to Ubound2do
begin
A^^[J] := 0.0;
for K := Lbound to Ubound2do
if S^[K] > 0.0 then
A^^[J] := A^^[J] + U^^[K] * V^[J]^[K];
end;
end;

function QR_Decomp(A: PMatrix;
Lbound, Ubound1, Ubound2: Integer;
R: PMatrix): Integer;
var
I, J, K: Integer;
Sum: Float;
begin
for K := Lbound to Ubound2do
begin
{ Compute the "k"th diagonal entry in R }
Sum := 0.0;
for I := Lbound to Ubound1do
Sum := Sum + Sqr(A^^[K]);
if Sum = 0.0 then
begin
QR_Decomp := MAT_SINGUL;
Exit;
end;

R^[K]^[K] := Sqrt(Sum);
{ Divide the entries in the "k"th column of A by the computed "k"th }
{ diagonal element of R. this begin
s the process of overwriting A }
{ with Q . . . }
for I := Lbound to Ubound1do
A^^[K] := A^^[K] / R^[K]^[K];
for J := (K + 1) to Ubound2do
begin
{ Complete the remainder of the row entries in R }
Sum := 0.0;
for I := Lbound to Ubound1do
Sum := Sum + A^^[K] * A^^[J];
R^[K]^[J] := Sum;
{ Update the column entries of the Q/A matrix }
for I := Lbound to Ubound1do
A^^[J] := A^^[J] - A^^[K] * R^[K]^[J];
end;
end;

QR_Decomp := MAT_OK;
end;

procedure QR_Solve(Q, R: PMatrix;
B: PVector;
Lbound, Ubound1, Ubound2: Integer;
X: PVector);
var
I, J: Integer;
Sum: Float;
begin
{ Form Q'B and store the result in X }
for J := Lbound to Ubound2do
begin
X^[J] := 0.0;
for I := Lbound to Ubound1do
X^[J] := X^[J] + Q^^[J] * B^;
end;

{ Update X with the solution vector }
X^[Ubound2] := X^[Ubound2] / R^[Ubound2]^[Ubound2];
for I := (Ubound2 - 1)do
wnto Lbounddo
begin
Sum := 0.0;
for J := (I + 1) to Ubound2do
Sum := Sum + R^^[J] * X^[J];
X^ := (X^ - Sum) / R^^;
end;
end;

end.

多人接受答案了。