Модуль реализации матричных вычислений для массивов больших размеров
Материал из DRKB
В этом модуле «осели» все операции с матрицами и векторами, которые я использовал для работы. Но есть алгоритмы, которые многие, наверняка, увидят впервые: Divide - алгоритм прямого деления, MSqrt - квадратный корень, MAbs - абсолютная величина. Поскольку модуль содержит все, от элементарных операций до матричных, разобраться будет несложно:
Например, решение системы ЛУ (консольное приложение)
var
N: Integer;
A: Matrix;
b, x: Vector;
begin
N := {..};
A.Init(N, N);
b.Init(N);
x.Init(N); // или x.Init(B); или x.InitRow(A);
{...}
{ формирование A и b }
{...}
x.Divide( b, A );
x.Print;
{...}
end.
Некоторые алгоритмы требуют пояснения, например:
Matrix.E(i, j: LongWord)
или
Vector.E(i: Integer): RealPtr
(RealPtr = ^Real) функция для вычисления адреса элемента матрицы/вектора. Перешла из ДОС когда в модуле использовался алгоритм управления виртуальной памятью для больших размерностей.
Matrix.Multiple(X, Y: Vector)
Результатом, которого является произведение вектора X на транспонированный вектор Y - матрица ранга 1.
Matrix.Invert(A: Matrix)
- если A[N,M], и N <> M то результат - матрица размера [M,N] - псевдообратная = A+.
Matrix.Addition(A: Matrix; B: Real)
- добавление числа в главную диагональ.
Matrix.Diag(r: Real)
- присваивание значения главной диагонали.
unit HMatrixW;
{
Матричная алгебра 1996-1998,
HNV 1996 .
}
interface
uses
SysUtils,
HNVString,
Math,
// VHeapLow,
Windows,
Dialogs;
const
MaxRow = 214748364;
{ Максимальное к-во элементов в векторе }
NumbSpline: longword = 1;
Indicator: word = $118F;
Condition: Boolean = True;
InitCount: Integer = 0;
InitSize: Integer = 0;
var
NIteration: longword; { Число итераций процесса }
Temp: Real;
MathEps: Real; { Относительная машинная точность }
SqrtEps: Real; { и корень из нее }
Eps: Real;
MinEps: Real;
type
VirtualPtr = record
Ptr: Pointer;
Size: longword;
end;
RealPtr = ^Real;
ArrayType =
array[1..1] of Real;
ArrayPtr = ^ArrayType;
VectorPrim = { Вектор } object
VAddr: VirtualPtr; { Виртуальный адрес (tail of DOS) }
FInit: word; { Флаг инициализации (tail of DOS) }
n: longword; { Размерность вектора }
function Addr: ArrayPtr; { Адрес начала вектора }
function E(I: longword): RealPtr; { Адрес элемента }
{ _ }
procedure Copy(X: VectorPrim); { = X }
{ _ _ }
procedure Substrac(X, Y: VectorPrim); { X - Y }
{ _ _ }
procedure Addition(X, Y: VectorPrim); { X + Y }
procedure Init(Size: longword); overload;
procedure Init(X: VectorPrim); overload;
procedure Free;
function Summ: Real; { Сумма элементов вектора }
function SummSqr: Real; { Сумма квадратов элементов }
procedure Sort(a1: VectorPrim);
{ Сортировка вектора по возрастанию }
procedure SortAbs(a1: VectorPrim); { -<>- по модулю }
procedure SortBy(a1: VectorPrim);
{ Сортировка вектора по возрастанию }
procedure Clear; { Очистить нулями }
procedure Print;
procedure PrintF(M: Byte);
procedure ScMultiple(X: Real; y: VectorPrim);
{ Умножение на скаляр }
procedure ScDivide(X: Real; y: VectorPrim);
{ Деление на скаляр }
procedure SetAll(X: Real); { Все эл-ты = X }
procedure Multiple(x, y: VectorPrim); overload;
{ поэлементное умножение }
end;
MatrixPrimArray =
array[1..1] of VectorPrim;
MatrixVPrimPtr = ^MatrixPrimArray;
Matrix = { Матрица } object
VAddr: VirtualPtr; { Виртуальный адрес массива векторов
матрица хранится по строкам. Максималь-
ная размерность MaxRow * MaxRow (tail of DOS)}
FInit: word; { Флаг инициализации (tail of DOS) }
n, m: longword; { Размерности матрицы }
function Addr: MatrixVPrimPtr;
{ Адрес массива векторов-строк }
function E(i, j: longword): RealPtr; { Адрес элемента }
procedure Copy(B: Matrix); { = B }
procedure Multiple(A, B: Matrix); overload; { A * B }
procedure Multiple(A, B: VectorPrim); overload;
procedure Substrac(A, B: Matrix); { A - B }
procedure Addition(A, B: Matrix); overload; { A + B }
procedure Addition(A: Matrix; B: VectorPrim); overload;
procedure Addition(A: Matrix; B: Real); overload;
procedure Divide(Bx, Ax: Matrix); { B / A }
procedure Diag(B: VectorPrim); overload; { Diag = b }
procedure Diag(r: Real); overload; { Diag(*) = r }
procedure Col(i: longword; B: VectorPrim); overload; { Столбец = b }
procedure Row(i: longword; B: VectorPrim); overload; { Строка = b }
procedure Col(i: longword; r: Real); overload; { Столбец = r }
procedure Row(i: longword; r: Real); overload; { Строка = r }
procedure Clear; { Обнуление эл-тов }
procedure Invert(A: Matrix);
{ Вычисление обратной/псевдообратной матрицы }
procedure Print;
procedure PrintF(M2: Byte); { M2 - к-во знач. цифр }
procedure ScMultiple(X: Real; Y: Matrix);
{ умножение на скаляр }
procedure Trans(A: Matrix); { Транспонирование }
procedure MSqrt(AX: Matrix);
{ Вычисление корня для }
{ положительно-определенной матрицы }
procedure MAbs(AX: Matrix); { Матричная функция ABS }
function Cond: Real; { числo обусловленности (1-Ok)}
function Norm: Real; { Норма матрицы }
function Alpha: Real; { Параметр регуляризации }
procedure Init(Size1, Size2: longword); overload;
{ Размещение матрицы }
procedure Init(A: Matrix); overload;
procedure InitTrans(A: Matrix);
procedure InitSqr(A: Matrix);
procedure ReadFromFile(FN: string); { Читать из файла TXT,
Init не требуется }
procedure MSqr(A: Matrix); { = Trans( A ) * A }
procedure Simmetry(A: Matrix); { = ( Trans( A ) + A ) / 2 }
function Asimm: Real; { Ассимметрия матрицы }
function Difference(A: Matrix): Real;
procedure Free; { Освобождение памяти }
end;
Vector = object(VectorPrim)
procedure Divide(x: Vector; AB: Matrix); overload; { =x / AB }
procedure Divide(x, y: Vector); overload;
procedure Multiple(A: Matrix; x: Vector); overload; { =A * x }
procedure Diag(A: Matrix); { =Диаголаль }
procedure Col(i: longword; A: Matrix); { =Столбец }
procedure Row(i: longword; A: Matrix); { =Строка }
function Max: Real; { Максимальный элемент }
function Min: Real; { Минимальный элемент }
procedure Abs(x: Vector); { Абсолютная величина }
function MinIndex: longword; { Индекс мин.эл-та }
function MaxIndex: longword; { Индекс макс. эл-та }
function Poly(x: Real): Real; { Полином x }
procedure IntPower(x: Vector; P: Integer);
{ Целая степенть(поэлементно) }
procedure SubVector(x: Vector; From, Count: longword);
procedure InitRow(A: Matrix); overload;
procedure InitRow(i: Integer; A: Matrix); overload;
procedure InitCol(A: Matrix); overload;
procedure InitCol(i: Integer; A: Matrix); overload;
procedure FillBodyTo(XMax: Real);
procedure FillBody;
function Difference(x: Vector): Real;
procedure TableDiff(x, y: Vector); // Две программы из
procedure Smoothe(y: Vector); // Фортран-пакета(старые)
end;
SplineCoeffType = record
NPoint: longword;
h, d, bb, d1, dg, d2: Vector;
end;
SplineCoeffPtr = ^SplineCoeffType;
Spline = { Кубические сплайны } object
VAddr: VirtualPtr;
function Addr: SplineCoeffPtr;
procedure Init(X, Y: Vector);
function F(X: Real): Real; { Значение сплайна в точке X }
function Integral(A, B: Real): Real;
overload; { Интеграл от А до В }
function Integral: Real; overload; { Интеграл }
function FMin(A, B: Real): Real;
function FMax(A, B: Real): Real;
function ZeroIn(A, B: Real): Real;
function DF(X: Real): Real;
{ Значение первой производной }
procedure Free;
private
function Integr(A, B: Real): Real; { Интеграл от А до В }
end;
PolyApprox = { Аппроксимация полиномом } object
Coeff: Vector;
dCoeff: Vector;
PlanM: Matrix;
procedure Init(xt, yt: Vector; NP: Integer);
function PX(xr: Real): Real;
function DPX(xr: Real): Real;
procedure Free;
private
XMin, XMax: Real;
end;
Eigen = { Собственные значения ( для симметричных матриц ) }
object
private
function F(Alpha: Real): Real;
public
Values: Vector;
Vectors: Matrix;
AbsValues: Vector;
Cond: Real;
D, R, P: Matrix; { D, R и P - факторы : D = R * R; A = D + R + P; }
Gamma: Matrix; { Матрица ошибок извлеч. из A }
procedure Init(A: Matrix);
procedure Free;
end;
{*--------------------------------------------------------------------*}
implementation
uses
HMinZero; // одномерный поиск
procedure Error(E: word);
begin
ShowMessage('Internal Erorr : ' + IntToStr(E));
Halt(E);
end;
procedure GetMem(var P: Pointer; Size: Integer); { (tail of DOS) }
begin
Inc(InitSize, Size);
Inc(InitCount);
System.GetMem(P, Size);
end;
procedure FreeMem(P: Pointer; Size: Integer); { (tail of DOS) }
begin
Dec(InitSize, Size);
Dec(InitCount);
System.FreeMem(P, Size);
end;
procedure Vector.TableDiff(x, y: Vector);
{
C SUBROUTINE DDGT3
C
C PURPOSE
C TO COMPUTE A VECTOR OF DERIVATIVE VALUES GIVEN VECTORS OF
C ARGUMENT VALUES AND CORRESPONDING FUNCTION VALUES.
C
C USAGE
C CALL DDGT3(X,Y,Z,NDIM,IER)
C
C DESCRIPTION OF PARAMETERS
C X - GIVEN VECTOR OF DOUBLE PRECISION ARGUMENT VALUES
C (DIMENSION NDIM)
C Y - GIVEN VECTOR OF DOUBLE PRECISION FUNCTION VALUES
C CORRESPONDING TO X (DIMENSION NDIM)
C Z - RESULTING VECTOR OF DOUBLE PRECISION DERIVATIVE
C VALUES (DIMENSION NDIM)
C NDIM - DIMENSION OF VECTORS X,Y AND Z
C IER - RESULTING ERROR PARAMETER
C IER = -1 - NDIM IS LESS THAN 3
C IER = 0 - NO ERROR
C IER POSITIVE - X(IER) = X(IER-1) OR X(IER) =
C X(IER-2)
C
C REMARKS
C (1) IF IER = -1,2,3, THEN THERE IS NO COMPUTATION.
C (2) IF IER = 4,...,N, THEN THE DERIVATIVE VALUES Z(1)
C ,..., Z(IER-1) HAVE BEEN COMPUTED.
C (3) Z CAN HAVE THE SAME STORAGE ALLOCATION AS X OR Y. IF
C X OR Y IS DISTINCT FROM Z, THEN IT IS NOT DESTROYED.
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C EXCEPT AT THE ENDPOINTS X(1) AND X(NDIM), Z(I) IS THE
C DERIVATIVE AT X(I) OF THE LAGRANGIAN INTERPOLATION
C POLYNOMIAL OF DEGREE 2 RELEVANT TO THE 3 SUCCESSIVE POINTS
C (X(I+K),Y(I+K)) K = -1,0,1. (SEE HILDEBRAND, F.B.,
C INTRODUCTION TO NUMERICAL ANALYSIS, MC GRAW-HILL, NEW YORK/
C TORONTO/LONDON, 1956, PP. 64-68.)
C
C ..................................................................
C
SUBROUTINE DDGT3(X,Y,Z,NDIM,IER)}
var
DY1, DY2, DY3, A, B: Real;
i: Integer;
begin
{ TEST OF DIMENSION AND ERROR EXIT IN CASE NDIM IS LESS THAN 3}
if x.n <> y.n then
Error(8001);
{ PREPARE DIFFERENTIATION LOOP}
A := X.e(1)^;
B := Y.e(1)^;
I := 2;
DY2 := X.e(2)^ - A;
if DY2 = 0 then
Error(8002);
DY2 := (Y.e(2)^ - B) / DY2;
{ START DIFFERENTIATION LOOP}
for I := 3 to x.n do
begin
A := X.e(I)^ - A;
if A = 0 then
Error(8002);
A := (Y.e(I)^ - B) / A;
B := X.e(I)^ - X.e(I - 1)^;
if B = 0 then
Error(8002);
DY1 := DY2;
DY2 := (Y.e(I)^ - Y.e(I - 1)^) / B;
DY3 := A;
A := X.e(I - 1)^;
B := Y.e(I - 1)^;
if (I - 3) <= 0 then
e(1)^ := DY1 + DY3 - DY2;
e(I - 1)^ := DY1 + DY2 - DY3;
end;
e(n)^ := DY2 + DY3 - DY1;
end;
procedure Vector.Smoothe(y: Vector);
{
C SUBROUTINE DSE15
C
C PURPOSE
C TO COMPUTE A VECTOR OF SMOOTHED FUNCTION VALUES GIVEN A
C VECTOR OF FUNCTION VALUES WHOSE ENTRIES CORRESPOND TO
C EQUIDISTANTLY SPACED ARGUMENT VALUES.
C
C USAGE
C CALL DSE15(Y,Z,NDIM,IER)
C
C DESCRIPTION OF PARAMETERS
C Y - GIVEN VECTOR OF DOUBLE PRECISION FUNCTION VALUES
C (DIMENSION NDIM)
C Z - RESULTING VECTOR OF DOUBLE PRECISION SMOOTHED
C FUNCTION VALUES (DIMENSION NDIM)
C NDIM - DIMENSION OF VECTORS Y AND Z
C IER - RESULTING ERROR PARAMETER
C IER = -1 - NDIM IS LESS THAN 5
C IER = 0 - NO ERROR
C
C REMARKS
C (1) IF IER=-1 THERE HAS BEEN NO COMPUTATION.
C (2) Z CAN HAVE THE SAME STORAGE ALLOCATION AS Y. IF Y IS
C DISTINCT FROM Z, THEN IT IS NOT DESTROYED.
C
C SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED
C NONE
C
C METHOD
C IF X IS THE (SUPPRESSED) VECTOR OF ARGUMENT VALUES, THEN
C EXCEPT AT THE POINTS X(1),X(2),X(NDIM-1) AND X(NDIM), EACH
C SMOOTHED VALUE Z(I) IS OBTAINED BY EVALUATING AT X(I) THE
C LEAST-SQUARES POLYNOMIAL OF DEGREE 1 RELEVANT TO THE 5
C SUCCESSIVE POINTS (X(I+K),Y(I+K)) K = -2,-1,...,2. (SEE
C HILDEBRAND, F.B., INTRODUCTION TO NUMERICAL ANALYSIS,
C MC GRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP. 295-302.)
C
C ..................................................................
C
SUBROUTINE DSE15(Y,Z,NDIM,IER)
}
var
A, B, C: Real;
I, NDIM: Integer;
begin
{ TEST OF DIMENSION }
NDIM := y.n;
if (NDIM < 5) then
Error(8004);
{ PREPARE LOOP }
A := Y.e(1)^ + Y.e(1)^;
C := Y.e(2)^ + Y.e(2)^;
B := 0.2 * (A + Y.e(1)^ + C + Y.e(3)^ - Y.e(5)^);
C := 0.1 * (A + A + C + Y.e(2)^ + Y.e(3)^ + Y.e(3)^ + Y.e(4)^);
{ START LOOP }
for I := 5 to NDIM do
begin
A := B;
B := C;
C := 0.2E0 * (Y.e(I - 4)^ + Y.e(I - 3)^ + Y.e(I - 2)^
+ Y.e(I - 1)^ + Y.e(I)^);
e(I - 4)^ := A;
end;
{ UPDATE LAST FOUR COMPONENTS}
A := Y.e(NDIM)^ + Y.e(NDIM)^;
A := 0.10 * (A + A + Y.e(NDIM - 1)^ + Y.e(NDIM - 1)^ +
Y.e(NDIM - 1)^ + Y.e(NDIM - 2)^ + Y.e(NDIM - 2)^ + Y.e(NDIM - 3)^);
e(NDIM - 3)^ := B;
e(NDIM - 2)^ := C;
e(NDIM - 1)^ := A;
e(NDIM)^ := A + A - C;
end;
procedure Jacobi(AT, S: Matrix; L: Vector);
var
i, j, k, p, q: word;
SinT, CosT, Sin2T, Cos2T, Lambda, Mu, Omega, AMax, App,
Aqq, Apq, Temp: Real;
T: Vector;
A, B, R: Matrix;
begin
A.Init(AT);
B.Init(A);
S.Clear;
for i := 1 to S.n do
S.E(i, i)^ := 1;
R.Init(S);
T.InitRow(R);
NIteration := 0;
repeat
Inc(NIteration);
B.Copy(A);
R.Copy(S);
p := 1;
q := 2;
AMax := A.E(p, q)^;
T.Diag(A);
Temp := T.Max;
for i := 1 to A.n do
for j := i + 1 to A.n do
if Abs(A.E(i, j)^) > AMax then
begin
p := i;
q := j;
AMax := Abs(A.E(p, q)^);
end;
if (AMax + Temp) = Temp then
Break;
App := A.E(p, p)^;
Aqq := A.E(q, q)^;
Apq := A.E(p, q)^;
Lambda := -Apq;
Mu := (App - Aqq) / 2.0;
Omega := Lambda / Sqrt(Sqr(Lambda) + Sqr(Mu));
if Mu < 0 then
Omega := -Omega;
SinT := Omega / Sqrt(2 * (1 + Sqrt(1 - Sqr(Omega))));
Sin2T := Sqr(SinT);
CosT := Sqrt(1 - Sin2T);
Cos2T := Sqr(CosT);
for k := 1 to A.n do
begin
B.E(p, k)^ := A.E(p, k)^ * CosT - A.E(q, k)^ * SinT;
B.E(q, k)^ := A.E(p, k)^ * SinT + A.E(q, k)^ * CosT;
end;
for i := 1 to A.n do
begin
B.E(i, p)^ := A.E(i, p)^ * CosT - A.E(i, q)^ * SinT;
B.E(i, q)^ := A.E(i, p)^ * SinT + A.E(i, q)^ * CosT;
R.E(i, p)^ := S.E(i, p)^ * CosT - S.E(i, q)^ * SinT;
R.E(i, q)^ := S.E(i, p)^ * SinT + S.E(i, q)^ * CosT;
end;
B.E(p, p)^ := App * Cos2T + Aqq * Sin2T - 2 * Apq * SinT * CosT;
B.E(q, q)^ := App * Sin2T + Aqq * Cos2T + 2 * Apq * SinT * CosT;
B.E(p, q)^ := 0;
B.E(q, p)^ := 0;
A.Copy(B);
S.Copy(R);
until False;
L.Diag(A);
B.Free;
R.Free;
A.Free;
T.Free;
end;
function Eigen.F(Alpha: Real): Real;
var
Q, T, X: Vector;
St, S: Matrix;
Z: Real;
begin
Inc(NIteration);
S.Init(Gamma);
St.InitTrans(Vectors);
Q.InitRow(Gamma);
T.Init(Q);
X.Init(Q);
S.Clear;
S.Diag(Exp(Alpha));
S.Multiple(Vectors, S);
S.Multiple(S, St);
T.Diag(Gamma);
T.SetAll(1);
Q.Multiple(Gamma, T);
S.Addition(Gamma, S);
X.Divide(Q, S);
T.Substrac(X, T);
Z := T.SummSqr;
Result := Ln(Z);
Q.Free;
T.Free;
S.Free;
ST.Free;
X.Free;
end;
procedure Eigen.Init(A: Matrix);
var
Iter: longword;
Q, T, X: Vector;
St, S: Matrix;
xZ: Real;
begin
if A.n <> A.m then { Матрица квадратная ? }
Error(4001);
if A.n < 2 then
Error(4003);
Values.InitRow(A);
Vectors.Init(A);
Jacobi(A, Vectors, Values);
AbsValues.Init(Values);
AbsValues.Abs(Values);
if AbsValues.Min <> 0 then
Cond := AbsValues.Max / AbsValues.Min
else
Cond := 1E+38;
D.Init(A);
R.Init(D);
P.Init(R);
Gamma.Init(A);
Q.InitRow(A);
X.Init(Q);
T.Init(Q);
St.Init(A);
S.Init(St);
Iter := NIteration;
D.Clear;
D.Diag(1);
R.ScMultiple(4, R);
R.Addition(R, D);
R.MSqrt(R);
R.Substrac(R, D);
R.ScMultiple(0.5, R);
D.MSqr(R);
P.Substrac(A, R);
P.Substrac(P, D);
NIteration := Iter + NIteration;
xZ := Exp(Ln(MathEps) * (Ln(1 / MathEps) /
(Ln(Cond) + Ln(1 / MathEps))));
Q.SetAll(Exp(xZ));
S.Clear;
S.Diag(Q);
St.Trans(Vectors);
S.Multiple(Vectors, S);
S.Multiple(S, St);
Gamma.Copy(S);
Q.Free;
T.Free;
X.Free;
St.Free;
S.Free;
end;
procedure Vector.Divide(x, y: Vector);
var
i: longword;
begin
if x.n <> y.n then
Error(4010);
for i := 1 to x.n do
E(i)^ := x.E(i)^ / y.E(i)^;
end;
procedure Eigen.Free;
begin
Values.Free;
Vectors.Free;
AbsValues.Free;
D.Free;
R.Free;
P.Free;
Gamma.Free;
end;
procedure Matrix.Simmetry(A: Matrix);
var
T: Matrix;
begin
if A.n <> A.m then
Error(5001);
if n <> A.n then
Error(5002);
T.InitTrans(A);
T.Addition(T, A);
ScMultiple(0.5, T);
T.Free;
end;
procedure Matrix.MAbs(AX: Matrix);
var
T: Matrix;
begin
T.InitSqr(AX);
MSqrt(T);
T.Free;
end;
function Spline.FMin(A, B: Real): Real;
begin
Result := HMinZero.FMin(SELF, A, B);
end;
function Spline.FMax(A, B: Real): Real;
begin
Result := HMinZero.FMax(SELF, A, B);
end;
function Spline.ZeroIn(A, B: Real): Real;
begin
Result := FZeroIn(SELF, A, B);
end;
function Spline.Integr(A, B: Real): Real;
var
x1, x3, x5, F1, F3, F5: Real;
function Int(x1, x3, x5, F1, F3, F5, Eps: Real): Real;
var
h1, h2, f2, f4, x2, x4, s1, s2, s3: Real;
q1, q2, st: Real;
begin
h1 := X3 - X1;
h2 := h1 / 2.;
x2 := (X3 + X1) / 2.;
x4 := (X5 + X3) / 2.;
f2 := F(x2);
f4 := F(x4);
{ Формула Симпсона }
s1 := h1 * (f1 + 4. * f2 + f3) / 6.;
s2 := h1 * (f3 + 4. * f4 + f5) / 6.;
{ Формула Ньютона-Коттеса 5-го порядка }
s3 := h2 * 2. * (7. * f1 + 32. * f2 + 12. * f3 + 32. * f4 + 7. * f5) / 45.;
{ Условие окончания }
st := h2 / Abs(B - A);
q1 := Abs(s1 + s2 - s3);
q2 := 0.5 * Eps * Abs(s1 + s2 + s3) * st + Eps;
if (q1 > q2) and (h2 > Eps + Abs(B - A) * Eps) then
{ Рекурсия - если точность плохая }
Result := Int(x1, x2, x3, f1, f2, f3, Eps) + Int(
x3, x4, x5, f3, f4, f5, Eps)
else
{ За решения принять значение формулы 5-го порядка }
Result := s3;
end;
begin
Result := 0;
if A = B then
exit;
x1 := A;
x3 := B - (B - A) * 0.5;
x5 := B;
F1 := F(x1);
F3 := F(x3);
F5 := F(x5);
Result := Int(x1, x3, x5, F1, F3, F5, Eps);
end;
function Spline.Integral(A, B: Real): Real;
begin
Result := Integr(A, B);
end;
function Spline.Integral: Real;
begin
Result := Integr(Addr^.d1.e(1)^, Addr^.d1.e(Addr^.d1.n)^);
end;
procedure Vector.FillBodyTo(XMax: Real);
var
x, S: Real;
i: Integer;
begin
x := Addr^[1];
S := (XMax - x) / (n - 1);
for i := 1 to n do
begin
Addr^[i] := x;
x := x + S;
end;
end;
procedure Vector.FillBody;
begin
FillBodyTo(Addr^[n]);
end;
procedure Vector.SubVector(x: Vector; From, Count: longword);
var
i, j: longword;
begin
if (n < Count) or ((Count + From - 1) > x.n) then
Error(2004);
if (x.n < From) then
Error(2005);
j := 0;
for i := From to From + Count - 1 do
begin
Inc(j);
Addr^[j] := x.e(i)^;
end;
end;
procedure PolyApprox.Init(xt, yt: Vector; NP: Integer);
{
Аппроксимация зависимости Y( x )
полиномом степени NP - 1
Вход :
X, Y;
}
var
A: Matrix; { Матрица плана }
AT: Matrix; { --//-- Транспонированная }
B: Matrix;
BY: Vector;
t: Vector;
i, j: longword;
Z: Real;
begin
Coeff.Init(NP);
A.Init(xt.n, Coeff.n);
{
Формирование матрицы плана
}
t.Init(xt.n);
XMin := xt.Min;
XMax := xt.Max;
for i := 1 to Coeff.n do
begin
for j := 1 to xt.n do
t.E(j)^ := xt.e(j)^;
t.IntPower(t, i - 1);
A.Col(i, t);
end;
{ Обычный метод НК }
PlanM.Init(A); // Запомним матрицу плана
AT.InitTrans(A);
B.InitSqr(A); // B = AT * A
Coeff.Multiple(AT, yt); // Coeff = AT * y
Coeff.Divide(Coeff, B); { Coeff = Coeff / B }
{ но, можно и так : }
{ вычисление псевдообратной матрицы }
// AT.Invert( A );
// Coeff.Multiple( AT, yt ); // Coeff = A+ * yt
if NP > 1 then
begin
dCoeff.Init(Coeff.n - 1);
for i := 2 to Coeff.n do
dCoeff.e(i - 1)^ := Coeff.e(i)^ * (i - 1);
end;
t.Free;
A.Free;
AT.Free;
end;
function PolyApprox.PX(xr: Real): Real;
begin
Result := Coeff.Poly(xr);
end;
function PolyApprox.DPX(xr: Real): Real;
{
Производная полинома в "xr"
}
begin
Result := dCoeff.Poly(xr);
end;
procedure PolyApprox.Free;
begin
if Coeff.n > 1 then
dCoeff.Free;
Coeff.Free;
PlanM.Free;
end;
function Matrix.Difference(A: Matrix): Real;
begin
Result := SELF.Norm - A.Norm;
end;
function Vector.Difference(x: Vector): Real;
begin
Result := Sqrt(SELF.SummSqr) - Sqrt(x.SummSqr);
end;
function Matrix.Asimm: Real;
{
Ассимметрия матрицы
}
var
i, j: longword;
S, T: Real;
begin
if N <> M then
Error(2002);
S := 0;
T := S;
for i := 1 to N do
for j := 1 to N do
if i <> j then
begin
S := S + Math.IntPower(E(i, j)^ - E(j, i)^, 2);
{ T := T + ( E(i,j)^ + E(j,i)^ ) / 2; }
end;
S := Sqrt(S) {/ ( T + 1 )};
Result := S;
end;
procedure Vector.IntPower(x: Vector; P: Integer);
{
Поэлементное возведение в целую степень
}
var
i: longword;
t: Vector;
begin
t.Init(x.N);
for i := 1 to x.n do
t.E(i)^ := Math.IntPower(x.E(i)^, P);
Copy(t);
t.Free;
end;
procedure Matrix.MSqr(A: Matrix);
{ T
= A * A ( квадрат матрицы )
}
var
T: Matrix;
begin
if (N <> M) or (M <> A.M) then
Error(2001);
T.InitTrans(A);
Multiple(T, A);
T.Free;
end;
procedure Matrix.InitSqr(A: Matrix);
{ T
= A * A ( квадрат матрицы )
}
begin
Init(A.m, A.m);
MSqr(A);
end;
function Vector.Poly(x: Real): Real; { Полином x }
{
Poly = .e(1) + x * ( .e(2) + x * ( .e(3) + x * ( .e(4) + ...x * .e(n))...
}
var
i: Integer;
S: Real;
begin
S := e(n)^;
for i := n - 1 downto 1 do
S := e(i)^ + S * x;
Result := S;
end;
procedure Matrix.ReadFromFile(FN: string);
{
Чтение матрицы из текстового файла
( по строкам )
}
var
f: Text;
S: string;
i, j: longword;
Y: Matrix;
const
Delim: CharSet =
[' ', ';', #9];
begin
if FileExists(FN) then
begin
Assign(f, FN);
Reset(f);
n := 0;
j := n;
while not EOF(f) do
begin
ReadLn(f, S);
if (Trim(S) = '') then
Continue;
if (S[1] = '''') then
Continue;
Inc(n);
m := WordCount(S, Delim);
if j = 0 then
j := m
else
if j <> m then
begin
ShowMessage('Error in input file.(1)');
Halt;
end;
end;
Y.Init(n, m);
Y.Clear;
Reset(f);
n := 0;
while not EOF(f) do
begin
ReadLn(f, S);
if (Trim(S) = '') then
Continue;
if (S[1] = '''') then
Continue;
Inc(n);
for i := 1 to WordCount(S, Delim) do
Y.E(n, i)^ := StrToFloat(ExtractWord(i, S, Delim));
end;
Init(n, m);
Copy(Y);
Y.Free;
Close(f);
end
else
begin
ShowMessage('File : "' + FN + '" - not found.');
Halt;
end;
end;
function Vector.Max: Real;
{
Максимальный элемент вектора
}
var
i: longword;
P: ArrayPtr;
X: Real;
begin
P := Addr;
X := P^[1];
for i := 1 to n do
if P^[i] > X then
X := P^[i];
Result := X;
end;
function Vector.Min: Real;
{
Минимальный элемент вектора
}
var
i: longword;
P: ArrayPtr;
X: Real;
begin
P := Addr;
X := P^[1];
for i := 1 to n do
if P^[i] < X then
X := P^[i];
Result := X;
end;
procedure Vector.Abs(x: Vector);
{
.e(i)^ = Abs( x.e(i)^ )
}
var
i: longword;
P1: ArrayPtr;
P2: ArrayPtr;
begin
P1 := Addr;
P2 := x.Addr;
for i := 1 to n do
P1^[i] := System.Abs(P2^[i]);
end;
function Vector.MaxIndex: longword;
{
Индекс максимального элемента
}
var
i: longword;
j: longword;
P: ArrayPtr;
X: Real;
begin
P := Addr;
X := P^[1];
j := 1;
for i := 1 to n do
if P^[i] > X then
begin
X := P^[i];
j := i;
end;
Result := j;
end;
function Vector.MinIndex: longword;
{
Индекс минимального элемента
}
var
i: longword;
j: longword;
P: ArrayPtr;
X: Real;
begin
P := Addr;
X := P^[1];
j := 1;
for i := 1 to n do
if P^[i] < X then
begin
X := P^[i];
j := i;
end;
Result := j;
end;
function Matrix.Cond: Real;
{
Число обусловленности матрицы( произведение нормы обратной и нормы исходной )
-1
Cond = ||.|| * ||.||.
}
var
B, A: Matrix;
n0, n1: Real;
begin
n0 := SELF.Norm;
if n = m then
begin
B.Init(SELF);
A.Init(B);
B.Clear;
B.Diag(1);
B.Divide(B, A);
end
else
Error(1011);
n1 := B.Norm;
Result := n1 * n0;
B.Free;
A.Free;
end;
function Matrix.Norm: Real;
{
Норма Фробениуса :
= Sqrt( Summ( Sqr( .e( *, * ) ) ) )
}
var
i: Integer;
S: Real;
begin
S := 0;
for i := 1 to n do
S := S + Addr^[i].SummSqr;
Result := Sqrt(S);
end;
procedure Matrix.MSqrt(AX: Matrix);
{
HNV 1997
Квадратный корень матрицы A.
( предполагается, что матрица A -
положительно определенная. )
R такая, что T
R * R = A
Симметричное/асимметричное разложение матрицы.
Ассиметрия характерна для плохо обусловленных.
}
var
Y, X, R, T: Matrix;
Diff, DiffOld: Real;
Alphax, EAlpha: Real;
Reject: Boolean;
begin
Y.Init(n, n);
T.Init(Y);
X.Init(T);
R.Init(X);
X.ScMultiple(0.5, AX);
X.Addition(X, 1); { B = AX / 2 + 1 }
NIteration := 0;
Y.MSqr(X);
Y.Substrac(Y, AX);
DiffOld := Y.Norm;
Alphax := 0.5;
{ Обычная реализация метода Ньютона }
repeat
Inc(NIteration);
Y.MSqr(X);
Y.Substrac(Y, AX);
T.Trans(X);
Y.Divide(Y, T);
repeat
T.ScMultiple(Alphax, Y);
R.Substrac(X, T);
T.MSqr(R);
T.Substrac(T, AX);
Diff := T.Norm;
Reject := Diff > DiffOld;
if Reject then
Alphax := -Alphax / 17
else
Break;
EAlpha := Alphax * Alphax + 1;
until EAlpha = 1;
if (not Reject) and (Abs(Alphax) * (DiffOld - Diff) >
(DiffOld + 1) * Eps) then
begin
X.Copy(R);
DiffOld := Diff;
if Abs(Alphax) < 0.5 then
Alphax := Alphax * 2;
end
else
Break;
until NIteration > 1023;
// Обычно сходится за 20-30 итераций...
Trans(X); // - результат
Y.Free;
X.Free;
R.Free;
T.Free;
end;
procedure Matrix.Init(Size1, Size2: longword);
var
i: longword;
begin
if Size1 = 0 then
Error(2018);
m := Size2;
n := Size1;
VAddr.Size := n * SizeOf(VectorPrim);
GetMem(VAddr.Ptr, VAddr.Size);
// VAddr := GetVMem( n * SizeOf( VectorPrim ) );
{ Разпределить массив векторов-строк }
for i := 1 to n do
Addr^[i].Init(m);
FInit := Indicator; { Флаг инициализации }
end;
procedure Matrix.Init(A: Matrix);
begin
Init(A.n, A.m);
Copy(A);
end;
procedure Matrix.InitTrans(A: Matrix);
begin
Init(A.m, A.n);
Trans(A);
end;
procedure Matrix.Trans(A: Matrix);
{
Транспонирование матрицы A
}
var
i: longword;
x: Vector;
B: Matrix;
begin
if (n <> A.m) or (m <> A.n) then
Error(116);
x.Init(A.n);
B.Init(A.m, A.n);
for i := 1 to A.m do
begin
x.Col(i, A);
B.Row(i, x);
end;
Copy(B);
B.Free;
x.Free;
end;
procedure VectorPrim.Init(Size: longword);
begin
if Size = 0 then
Error(2017);
n := Size;
{ Распределить память для вектора }
VAddr.Size := n * SizeOf(real);
GetMem(VAddr.Ptr, VAddr.Size);
// VAddr := GetVMem( n * SizeOf( Real ) );
FInit := Indicator;
end;
procedure Vector.InitRow(A: Matrix);
begin
Init(A.m);
end;
procedure Vector.InitCol(A: Matrix);
begin
Init(A.n);
end;
procedure Vector.InitRow(i: Integer; A: Matrix);
begin
InitRow(A);
Row(i, A);
end;
procedure Vector.InitCol(i: Integer; A: Matrix);
begin
InitCol(A);
Col(i, A);
end;
procedure VectorPrim.Init(X: VectorPrim);
begin
n := X.n;
{ Распределить память для вектора }
VAddr.Size := n * SizeOf(real);
GetMem(VAddr.Ptr, VAddr.Size);
// VAddr := GetVMem( n * SizeOf( Real ) );
FInit := Indicator;
Copy(X);
end;
procedure VectorPrim.SetAll(X: Real);
var
i: longword;
t: ArrayPtr;
begin
t := Addr;
for i := 1 to n do
t^[i] := X;
end;
procedure Spline.Init(X, Y: Vector);
procedure SplineXY;
{
Вычисление коэффициентов кубического сплайна функции
заданной таблицей значений Y( x ).
}
var
i, j: longword;
f: Real;
begin
NumbSpline := 1;
with Addr^ do
begin
for i := 1 to NPoint - 1 do
begin
h.E(i)^ := x.E(i + 1)^ - x.E(i)^;
d.E(i)^ := (y.E(i + 1)^ - y.E(i)^) / h.E(i)^;
end;
dg.E(1)^ := 2.;
dg.E(NPoint)^ := 2.;
d2.E(1)^ := 1.;
d1.E(NPoint - 1)^ := 1.;
bb.E(1)^ := d.E(1)^ * 3.;
bb.E(NPoint)^ := d.E(NPoint - 1)^ * 3.;
for i := 2 to NPoint - 1 do
begin
bb.E(i)^ :=
(d.E(i)^ * h.E(i - 1)^ + d.E(i - 1)^ * h.E(i)^) * 3.;
d1.E(i - 1)^ := h.E(i)^;
d2.E(i)^ := h.E(i - 1)^;
dg.E(i)^ := 2. * (h.E(i)^ + h.E(i - 1)^);
end;
{
Решение системы линейных уравнений с 3-х диагональной
матрицей коэффициентов
}
for i := 1 to NPoint - 1 do
begin
f := d1.E(i)^ / dg.E(i)^;
dg.E(i + 1)^ := dg.E(i + 1)^ - d2.E(i)^ * f;
bb.E(i + 1)^ := bb.E(i + 1)^ - bb.E(i)^ * f;
end;
bb.E(NPoint)^ := bb.E(NPoint)^ / dg.E(NPoint)^;
for i := 1 to NPoint - 1 do
begin
j := NPoint - i;
bb.E(j)^ :=
(bb.E(j)^ - d2.E(j)^ * bb.E(j + 1)^) / dg.E(j)^;
end;
for i := 1 to NPoint do
begin
d1.Copy(x); { d1 - исходные значения X }
d2.Copy(y); { d2 - ---:---- Y }
end;
end;
end;
begin
VAddr.Size := SizeOf(SplineCoeffType);
GetMem(VAddr.Ptr, SizeOf(SplineCoeffType));
// VAddr := GetVMem( SizeOf( SplineCoeffType ) );
with Addr^ do
begin
NPoint := X.n;
H.Init(X.n);
D.Init(X.n);
bb.Init(X.n);
d1.Init(X.n);
dg.Init(X.n);
d2.Init(X.n);
end;
SplineXY;
end;
function Spline.F(X: Real): Real;
{
Вычисление значения сплайн-функции в точке X.
}
function SplineIndex: longword;
{
Вычисляет номер сплайна для точки 'X'.
}
var
NS: longword;
procedure IndRecurs(K, L: longword);
var
M: longword;
begin
with Addr^ do
begin
NS := K;
if (L - K) <= 1 then
Exit;
M := ((K + L) div 2);
if ((X - D1.E(K)^) * (X - D1.E(M)^)) <= 0.0 then
IndRecurs(K, M)
else
IndRecurs(M, L);
end;
end;
begin
with Addr^ do
begin
Result := 1;
if X <= d1.E(1)^ then
exit;
Result := NPoint - 1;
if X >= d1.E(NPoint)^ then
exit;
IndRecurs(1, NPoint);
Result := NS;
end;
end;
var
j: longword;
f1, t, s: Real;
begin
with Addr^ do
begin
j := SplineIndex;
t := (X - d1.E(j)^) / h.E(j)^;
f1 := 1. - t;
s := t * d2.E(j + 1)^ + f1 * d2.E(j)^;
s := s + h.E(j)^ * f1 * t * ((bb.E(j)^ - d.E(j)^) *
f1 - (bb.E(j + 1)^ - d.E(j)^) * t);
Result := s;
end;
end;
function Spline.DF(X: Real): Real;
{
Вычисление первой производной от сплайн-функции в точке 'X'.
}
var
Delta: Real;
Temp: Real;
PA: PolyApprox;
xS, yS: Vector;
procedure PolyCalc(n: Integer);
begin
xS.Init(5);
yS.Init(5);
xS.SubVector(Addr^.D1, n, 5);
yS.SubVector(Addr^.D2, n, 5);
PA.Init(xS, yS, 5);
Result := PA.DPX(X);
PA.Free;
yS.Free;
xS.Free;
end;
begin
Delta := Abs(X) * SqrtEps + SqrtEps;
with Addr^ do
begin
if X - Delta <= D1.E(1)^ then
Temp := -3 * F(X) + 4 * F(X + Delta) - F(X + 2 * Delta)
else
if X + Delta >= D1.E(Npoint)^ then
Temp := 3 * F(X) - 4 * F(X - Delta) + F(X - 2 * Delta)
else
Temp := (F(X + Delta) - F(X - Delta));
Result := Temp / Delta / 2;
end;
end;
procedure VectorPrim.Free;
begin
FInit := 0;
FreeMem(VAddr.Ptr, VAddr.Size);
VAddr.Ptr := nil;
// FreeVMem( VAddr );
// VAddr := 0;
end;
procedure Matrix.Free;
var
i: longword;
begin
FInit := 0;
for i := 1 to n do
Addr^[i].Free;
FreeMem(VAddr.Ptr, VAddr.Size);
VAddr.Ptr := nil;
//FreeVMem( VAddr );
// VAddr := 0;
end;
procedure Spline.Free;
begin
with Addr^ do
begin
H.Free;
D.Free;
bb.Free;
d1.Free;
dg.Free;
d2.Free;
end;
FreeMem(VAddr.Ptr, VAddr.Size);
VAddr.Ptr := nil;
//FreeVMem( VAddr );
// VAddr := 0;
end;
function VectorPrim.Addr: ArrayPtr;
begin
Result := VAddr.Ptr;
//Addr := VirtualToPtr( VAddr );
end;
function Matrix.Addr: MatrixVPrimPtr;
begin
Result := VAddr.Ptr;
//Addr := VirtualToPtr( VAddr );
end;
function Matrix.E(I, J: longword): RealPtr;
begin
if (i >= 1) and (i <= n) and (j >= 1) and (j <= m) then
Result := System.Addr(Addr^[i].E(j)^) { Адрес элемента I,J }
else
Error(4);
end;
function VectorPrim.E(I: longword): RealPtr;
begin
if (i >= 1) and (i <= n) then
Result := System.Addr(Addr^[i]) { Адрес элемента I }
else
Error(5);
end;
procedure VectorPrim.Copy(X: VectorPrim);
begin
if FInit <> Indicator then
Error(100);
Move(X.Addr^, Addr^, X.n * SizeOf(real)); { Быстрая пересылка }
end;
procedure Matrix.Copy(B: Matrix);
var
i: longword;
begin
if FInit <> Indicator then
Error(101);
for i := 1 to n do
Addr^[i].Copy(B.Addr^[i]); { Векторное копирование }
end;
procedure Matrix.Multiple(A, B: Matrix);
{
Вычисление произведения двух матриц
= A * B
}
var
i, j: longword;
TX: Matrix;
t1, t2: Vector;
begin
if FInit <> Indicator then
Error(102);
TX.Init(A.n, B.m);
t1.Init(A.m);
t2.Init(B.n);
for i := 1 to A.n do
begin
t1.Row(i, A);
for j := 1 to B.m do
begin
t2.Col(j, B);
t2.Multiple(t2, t1);
TX.E(i, j)^ := t2.Summ;
end;
end;
Copy(TX);
TX.Free;
t1.Free;
t2.Free;
end;
procedure Matrix.Multiple(A, B: VectorPrim);
{
T
Вычисление произведения двух векторов A * B
}
var
i, j: longword;
begin
if FInit <> Indicator then
Error(102);
if A.n <> B.n then
Error(4001);
if n <> m then
Error(4002);
if n <> A.n then
Error(4003);
for i := 1 to n do
for j := 1 to n do
E(i, j)^ := A.E(i)^ * B.E(j)^;
end;
procedure Vector.Multiple(A: Matrix; x: Vector);
{
Вычисление произведения матрицы на вектор
= A * x
}
var
i: longword;
T1, T2: Vector;
begin
if x.n <> A.m then
Error(6);
if FInit <> Indicator then
Error(103);
T1.Init(A.n);
T2.Init(A.m);
for i := 1 to A.n do
begin
T2.Row(i, A); { Быстрая пересылка }
T2.Multiple(T2, x); { Умножение }
T1.Addr^[i] := T2.Summ; { Суммирование }
end;
Copy(T1);
T1.Free;
T2.Free;
end;
procedure VectorPrim.Multiple(x, y: VectorPrim);
var
i: longword;
yp, xp, Sp: ArrayPtr;
begin
if x.n <> y.n then
Error(121);
Sp := Addr;
xp := x.Addr;
yp := y.Addr;
for i := 1 to x.n do
Sp^[i] := xp^[i] * yp^[i];
end;
procedure Vector.Divide(x: Vector; AB: Matrix);
{
Решение системы линейных алгебраических уравнений
вида:
A * y = b
методом прямого деления.
}
var
b: Matrix;
begin
if (x.n <> AB.m) or (AB.m <> AB.n) then
Error(7);
if FInit <> Indicator then { Новый вектор }
Error(104);
b.Init(AB.n, 1);
b.Col(1, x);
b.Divide(b, AB);
Col(1, b);
b.Free;
end;
procedure Matrix.Divide(Bx, Ax: Matrix);
{
Решение матричной системы вида :
A * Y = B, где
A[n,n], Y[n,m], B[n,m] m<=n
может быть применена для обращения матрицы
если В = Е.
}
var
i, j, k: longword;
d, t: Real;
a, b: Matrix;
FBreak: Boolean;
begin
if FInit <> Indicator then { Новая матрица }
Error(105);
a.Init(Ax.n, Ax.m);
b.Init(Bx.n, Bx.m);
a.Copy(Ax);
b.Copy(Bx);
Condition := True;
MinEps := 1.0;
FBreak := False;
for i := 1 to a.n do
begin
if System.Abs(a.E(i, i)^) = 0.0 then
begin
Condition := Condition and (i = 1);
for j := 1 to a.n do
begin
if System.Abs(a.E(j, i)^) <> 0.0 then
begin
a.Addr^[i].Addition(a.Addr^[i], a.Addr^[j]);
b.Addr^[i].Addition(b.Addr^[i], b.Addr^[j]);
Break;
end;
end;
FBreak := True;
end;
if FBreak then
Continue;
for k := 1 to b.m do
b.E(i, k)^ := b.E(i, k)^ / a.E(i, i)^;
for k := a.n downto i do
a.E(i, k)^ := a.E(i, k)^ / a.E(i, i)^;
for j := 1 to a.n do
begin
d := a.E(j, i)^;
if (j <> i) and (System.Abs(d) <> 0.0) then
begin
for k := i + 1 to a.n do
begin
t := a.E(j, k)^ / d - a.E(i, k)^;
a.E(j, k)^ := t;
end;
for k := 1 to b.m do
b.E(j, k)^ := b.E(j, k)^ / d - b.E(i, k)^;
if j < i then
a.E(j, j)^ := a.E(j, j)^ / d;
end;
end;
end;
for i := 1 to a.n do
for k := 1 to b.m do
if a.E(i, i)^ <> 0 then
E(i, k)^ := b.E(i, k)^ / a.E(i, i)^
else
E(i, k)^ := 0;
a.Free;
b.Free;
end;
procedure Vector.Diag(A: Matrix);
{
= Diag( A );
}
var
i: longword;
P: ArrayPtr;
begin
if FInit <> Indicator then { Новый вектор }
Error(106);
P := Addr;
for i := 1 to n do
P^[i] := A.E(i, i)^;
end;
procedure Vector.Col(i: longword; A: Matrix);
{
= A( *, i )
}
var
j: longword;
P: ArrayPtr;
begin
if FInit <> Indicator then { Проверк инициализации }
Error(107);
P := Addr;
for j := 1 to A.n do
P^[j] := A.E(j, i)^;
end;
procedure Vector.Row(i: longword; A: Matrix);
{
= A( i, * );
}
begin
Copy(A.Addr^[i]);
end;
procedure VectorPrim.Addition(x, y: VectorPrim);
{
= x.e(*)^ + y.e(*)^;
}
var
j: longword;
xP, yP, sP: ArrayPtr;
begin
if (x.n <> y.n) then
Error(9);
if FInit <> Indicator then
Error(108);
xP := x.Addr;
yP := y.Addr;
sP := Addr;
for j := 1 to n do
sP^[j] := xP^[j] + yP^[j];
end;
procedure VectorPrim.Substrac(x, y: VectorPrim);
{
= x.e(*)^ - y.e(*)^;
}
var
j: longword;
xP, yP, sP: ArrayPtr;
begin
if (x.n <> y.n) then
Error(10);
if FInit <> Indicator then
Error(109);
xP := x.Addr;
yP := y.Addr;
sP := Addr;
for j := 1 to n do
sP^[j] := xP^[j] - yP^[j];
end;
function VectorPrim.Summ: Real;
{
Сумма эл-тов вектора
}
var
i: longword;
S: Real;
T: VectorPrim;
AP: ArrayPtr;
begin
S := 0.0;
T.Init(n);
AP := T.Addr;
Move(Addr^, AP^, n * SizeOf(real));
T.SortAbs(T); { Сортировка по возрастанию }
for i := 1 to n do { Суммируем начиная с меньших эл-тов }
S := S + AP^[i];
Result := S;
T.Free;
end;
function VectorPrim.SummSqr: Real;
{
Сумма квадратов эл-тов вектора
}
var
i: longword;
S: Real;
T: VectorPrim;
TP: ArrayPtr;
begin
S := 0.0;
T.Init(n);
TP := T.Addr;
Move(Addr^, TP^, SizeOf(real) * n);
T.SortAbs(T); { сортировка по возрастанию }
for i := 1 to n do
S := S + Sqr(TP^[i]);
Result := S;
T.Free;
end;
procedure Matrix.Addition(A, B: Matrix);
{
= A( *, * ) + B( *, * );
}
var
i: longword;
begin
if (A.n <> B.n) or (A.m <> B.m) then
Error(11);
if FInit <> Indicator then
Error(110);
for i := 1 to n do
Addr^[i].Addition(A.Addr^[i], B.Addr^[i]);
end;
procedure Matrix.Addition(A: Matrix; B: VectorPrim);
{
= A( *, * ) + B( *, * );
}
var
i: longword;
begin
if (A.n <> n) or (n <> m) then
Error(11);
if FInit <> Indicator then
Error(110);
for i := 1 to n do
E(i, i)^ := A.E(i, i)^ + B.E(i)^;
end;
procedure Matrix.Addition(A: Matrix; B: Real);
{
= A( *, * ) + B( *, * );
}
var
i: longword;
begin
if (A.n <> A.m) then
Error(11);
if FInit <> Indicator then
Error(110);
Copy(A);
for i := 1 to n do
E(i, i)^ := A.E(i, i)^ + B;
end;
procedure Matrix.Substrac(A, B: Matrix);
{
= A( *, * ) - B( *, * );
}
var
i: longword;
begin
if (A.n <> B.n) or (A.m <> B.m) then
Error(11);
if FInit <> Indicator then
Error(111);
for i := 1 to n do
Addr^[i].Substrac(A.Addr^[i], B.Addr^[i]);
end;
procedure Matrix.Diag(b: VectorPrim);
{
Diag(*) = b(*);
}
var
i: longword;
begin
if m <> b.n then
Error(12);
for i := 1 to n do
E(i, i)^ := b.E(i)^;
end;
procedure Matrix.Diag(r: Real);
{
Diag(*) = r;
}
var
i: longword;
begin
for i := 1 to n do
E(i, i)^ := r;
end;
procedure Matrix.Col(i: longword; b: VectorPrim);
{
.e( *, i ) = b( * );
}
var
j: longword;
begin
if n <> b.n then
Error(13);
for j := 1 to n do
E(j, i)^ := b.E(j)^;
end;
procedure Matrix.Row(i: longword; b: VectorPrim);
{
.e( i, * ) = b( * );
}
begin
if m <> b.n then
Error(13);
Addr^[i].Copy(b);
end;
procedure Matrix.Col(i: longword; r: Real);
{
.e( *, i ) = r;
}
var
j: longword;
begin
for j := 1 to n do
E(j, i)^ := r;
end;
procedure Matrix.Row(i: longword; r: Real);
{
.e( i, * ) = r;
}
var
j: longword;
begin
for j := 1 to m do
E(i, j)^ := r;
end;
function Spline.Addr: SplineCoeffPtr;
begin
Result := VAddr.Ptr;
//Addr := VirtualToPtr( VAddr );
end;
procedure VectorPrim.Sort(a1: VectorPrim);
{
Quick-сортировка по возрастанию.
}
var
a: ArrayPtr;
t: Vector;
procedure SortPrim(l, r: longword);
var
i, j: longword;
x, y: Real;
begin
i := l;
j := r;
x := a^[(l + r) div 2];
repeat
while a^[i] < x do
i := i + 1;
while x < a^[j] do
j := j - 1;
if i <= j then
begin
y := a^[i];
a^[i] := a^[j];
a^[j] := y;
i := i + 1;
j := j - 1;
end;
until i > j;
if l < j then
SortPrim(l, j);
if i < r then
SortPrim(i, r);
end;
begin {quicksort}
;
t.Init(a1.n);
t.Copy(a1);
a := t.Addr;
SortPrim(1, N);
Copy(t);
t.Free;
end;
procedure VectorPrim.SortBy(a1: VectorPrim);
{
Quick-сортировка по возрастанию.
}
var
a, b: ArrayPtr;
t: Vector;
procedure SortPrim(l, r: longword);
var
i, j: longword;
x, y: Real;
begin
i := l;
j := r;
x := a^[(l + r) div 2];
repeat
while a^[i] < x do
i := i + 1;
while x < a^[j] do
j := j - 1;
if i <= j then
begin
y := a^[i];
a^[i] := a^[j];
a^[j] := y;
y := b^[i];
b^[i] := b^[j];
b^[j] := y;
i := i + 1;
j := j - 1;
end;
until i > j;
if l < j then
SortPrim(l, j);
if i < r then
SortPrim(i, r);
end;
begin {quicksort}
t.Init(a1.n);
t.Copy(a1);
a := t.Addr;
b := Addr;
SortPrim(1, N);
{ Copy( t );}
t.Free;
end;
procedure VectorPrim.SortAbs(a1: VectorPrim);
{
Quick-сортировка по возрастанию ABS.
}
var
a: ArrayPtr;
t: Vector;
procedure SortPrim(l, r: longword);
var
i, j: longword;
x, y: Real;
begin
i := l;
j := r;
x := Abs(a^[(l + r) div 2]);
repeat
while Abs(a^[i]) < x do
i := i + 1;
while x < Abs(a^[j]) do
j := j - 1;
if i <= j then
begin
y := a^[i];
a^[i] := a^[j];
a^[j] := y;
i := i + 1;
j := j - 1;
end;
until i > j;
if l < j then
SortPrim(l, j);
if i < r then
SortPrim(i, r);
end;
begin {quicksort}
t.Init(a1.n);
t.Copy(a1);
a := t.Addr;
SortPrim(1, N);
Copy(t);
t.Free;
end;
procedure VectorPrim.Clear;
begin
SetAll(0);
end;
procedure Matrix.Clear;
var
i: longword;
begin
for i := 1 to n do
Addr^[i].Clear;
end;
function Matrix.Alpha: Real;
var
QCond, xZ, HNorm: Real;
Temp: Vector;
begin
Temp.Init(n);
QCond := SELF.Cond;
Temp.Diag(SELF);
HNorm := Sqrt(Temp.SummSqr);
xZ := Ln(MathEps) * (Ln(1 / MathEps) / (Ln(QCond) + Ln(1 / MathEps)));
Result := Exp(xZ) * HNorm;
Temp.Free;
end;
procedure Matrix.Invert(A: Matrix);
{
Метод деления для вычисления обратной/псевдообратной матрицы 'A'.
}
var
B: Matrix;
Q: Matrix;
R: Matrix;
QCond, RCond: Real;
begin
Q.InitSqr(A);
// Комбинированный метод корня и понижениея числа обусловленности
Q.Addition(Q, Q.Alpha);
R.Init(Q);
B.InitTrans(A);
QCond := Q.Cond;
// WriteLn('Cond(Q)=', QCond );
R.MSqrt(Q); // квадратный корень Q
RCond := R.Cond;
// WriteLn('Cond(R)=', RCond );
// Вычисление псевдообратной матрицы
B.Divide(B, R);
R.Trans(R);
Divide(B, R);
Q.Free;
R.Free;
B.Free;
end;
procedure VectorPrim.Print;
{ (tail of DOS) }
var
i: longword;
l: Integer;
begin
l := Trunc((79 - n) / n);
if l < 9 then
begin
WriteLn('Warning, length of output digit < 9. ');
for i := 1 to n do
WriteLn(i: 5, ' ', E(i)^);
exit;
end;
for i := 1 to n do
Write(E(i)^: l, ' ');
WriteLn;
end;
procedure VectorPrim.PrintF(M: Byte);
{
(tail of DOS)
}
var
i: longword;
l: Byte;
begin
l := Trunc((79 - n) / n);
if l - M < 2 then
begin
WriteLn('Warning, invalid format');
exit;
end;
if l < 7 then
begin
WriteLn('Warning, length of fixed digit < 7. ');
exit;
end;
if l > 18 then
l := 18;
for i := 1 to n do
Write(E(i)^: l: M, ' ');
WriteLn;
end;
procedure Matrix.Print;
{
(tail of DOS)
}
var
i: longword;
begin
for i := 1 to n do
Addr^[i].Print;
WriteLn;
end;
procedure Matrix.PrintF(M2: Byte);
{
(tail of DOS)
}
var
i: longword;
begin
for i := 1 to n do
Addr^[i].PrintF(M2);
WriteLn;
end;
procedure VectorPrim.ScMultiple(X: Real; y: VectorPrim);
var
i: longword;
t: ArrayPtr;
yp: ArrayPtr;
begin
t := Addr;
yp := y.Addr;
for i := 1 to n do
t^[i] := X * yp^[i]; { Быстрое умножение }
end;
procedure VectorPrim.ScDivide(X: Real; y: VectorPrim);
var
i: longword;
t: ArrayPtr;
yp: ArrayPtr;
begin
t := Addr;
yp := y.Addr;
for i := 1 to n do
t^[i] := yp^[i] / X; { Быстрое умножение }
end;
procedure Matrix.ScMultiple(X: Real; Y: Matrix);
var
i: longword;
begin
for i := 1 to n do
Addr^[i].ScMultiple(X, Y.Addr^[i]);
end;
begin
{
Вычисление точности:
MathEps = относительная машинная точность.
( -log( MathEps ) ~= к-во знаков в мантиссе. )
SqrtEps = MathEps ^ 1/2;
( - точность половина знаков мантиссы. )
Eps = MathEps ^ 3/4;
( --//-- 3/4 знаков мантиссы. )
}
MathEps := 1;
repeat
MathEps := MathEps / 2;
Temp := 1 + MathEps;
until Temp = 1;
SqrtEps := Sqrt(MathEps);
Eps := Sqrt(SqrtEps) * SqrtEps;
MinEps := 1;
Randomize;
Indicator := Random($7FFF);
//InitVHeap( '$$Matr.vhp', true );
end.
Этот модуль используется почти во всех реализованных мной численных алгоритмах и методах. Те части, которые писал не я - приводятся без изменений(по возможности) стиля и комментариев.
unit HNVString;
interface
{*********************************************************}
{* OPSTRING.PAS 1.20 *}
{* Copyright (c) TurboPower Software 1987, 1992. *}
{* All rights reserved. *}
{*********************************************************}
type
CharSet =
set of char;
{ Given a set of word delimiters, return number of words in S}
function WordCount(S: string; WordDelims: CharSet): Byte;
{ Given a set of word delimiters, return start position of N'th word in S}
function WordPosition(N: Byte; S: string; WordDelims: CharSet): Byte;
{ Given a set of word delimiters, return the N'th word in S}
function ExtractWord(N: Byte; S: string; WordDelims: CharSet): string;
function RightJust(S: string; N: word): string;
function HexB(B: Byte): string;
function HexText(var Buff; N: Integer): string;
implementation
const
Dig: array[0..$F] of char = '0123456789ABCDEF';
{ Convert BYTE to Hex String }
function HexB(B: Byte): string;
var
T: string;
begin
T := '';
T := T + Dig[B shr 4];
T := T + Dig[B and $0F];
HexB := T;
end;
{ Convert BUFFER to HEX String }
function HexText(var Buff; N: Integer): string;
type
ByteArray =
array of Byte;
var
i: Integer;
S: string;
B: ByteArray Absolute Buff;
begin
S := '';
for i := 1 to N do
S := S + HexB(B[i]);
HexText := S;
end;
function RightJust(S: string; N: word): string;
var
T: string;
i: Integer;
j: Integer;
begin
T := '';
j := N;
if Length(S) > j then
T := S
else
begin
j := j - Length(S);
for i := 1 to j do
T := T + ' ';
T := T + S;
end;
RightJust := T;
end;
function WordCount(S: string; WordDelims: CharSet): Byte;
{-Given a set of word delimiters, return number of words in S}
var
Count: Byte; {!!.02}
I: word; {!!.02}
SLen: Integer;
begin
SLen := Length(S);
Count := 0;
I := 1;
while I <= SLen do
begin
{skip over delimiters}
while (I <= SLen) and (S[I] in WordDelims) do
Inc(I);
{if we're not beyond end of S, we're at the start of a word}
if I <= SLen then
Inc(Count);
{find the end of the current word}
while (I <= SLen) and not (S[I] in WordDelims) do
Inc(I);
end;
WordCount := Count;
end;
function WordPosition(N: Byte; S: string; WordDelims: CharSet): Byte;
{-Given a set of word delimiters, return start position of N'th word in S}
var
Count: Byte; {!!.02}
I: word; {!!.02}
SLen: word;
begin
Count := 0;
I := 1;
SLen := Length(S);
WordPosition := 0;
while (I <= SLen) and (Count <> N) do
begin
{skip over delimiters}
while (I <= SLen) and (S[I] in WordDelims) do
Inc(I);
{if we're not beyond end of S, we're at the start of a word}
if I <= SLen then
Inc(Count);
{if not finished, find the end of the current word}
if Count <> N then
while (I <= SLen) and not (S[I] in WordDelims) do
Inc(I)
else
WordPosition := I;
end;
end;
function ExtractWord(N: Byte; S: string; WordDelims: CharSet): string;
{-Given a set of word delimiters, return the N'th word in S}
var
I: word; {!!.12}
SLen: Integer;
Temp: string;
begin
SLen := Length(S);
Temp := '';
I := WordPosition(N, S, WordDelims);
if I <> 0 then
{find the end of the current word}
while (I <= SLen) and not (S[I] in WordDelims) do
begin
{add the I'th character to result}
Temp := Temp + S[i];
Inc(I);
end;
ExtractWord := Temp;
end;
end.
Author: Andrey
Source: http://delphikingdom.com
ID: 03965