Модуль реализации матричных вычислений для массивов больших размеров

Previous  Top  Next

    
 

Автор: Andrey

WEB-сайт: http://delphikingdom.com

 

В этом модуле «осели» все операции с матрицами и векторами, которые я использовал для работы. Но есть алгоритмы, которые многие, наверняка, увидят впервые: Divide – алгоритм прямого деления, MSqrt – квадратный корень, MAbs – абсолютная величина. Поскольку модуль содержит все, от элементарных операций до матричных, разобраться будет несложно:

 

Например, решение системы ЛУ (консольное приложение)

Code:

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 )

 

– присваивание значения главной диагонали.

 

 

 

Code:

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;

 

  F := 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

  FMin := HMinZero.FMin( SELF, A, B );

end;

 

function Spline.FMax( A, B : Real ) : Real;

begin

  FMax := HMinZero.FMax( SELF, A, B );

end;

 

function Spline.ZeroIn( A, B : Real ) : Real;

begin

  ZeroIn := 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

        { Рекурсия - если точность плохая }

        Int := Int( x1, x2, x3, f1, f2, f3, Eps )

           + Int( x3, x4, x5, f3, f4, f5, Eps )

     else

 

        { За решения принять значение формулы 5-го порядка }

        Int := s3;

  end;

 

begin

 

  Integr := 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 );

 

  Integr := Int( x1, x3, x5, F1, F3, F5, Eps );

 

end;

 

function Spline.Integral( A, B : Real ) : Real;

begin

  Integral := Integr( A, B );

end;

 

function Spline.Integral : Real;

begin

  Integral := 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

  PX := Coeff.Poly( xr );

end;

 

function PolyApprox.DPX( xr : Real ) : Real;

{

  Производная полинома в "xr"

}

begin

  DPX := 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

  Difference := SELF.Norm - A.Norm;

end;

 

function Vector.Difference( x : Vector ) : Real;

begin

  Difference := 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 )};

  Asimm := 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;

 

  Poly := 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 ];

 

  Max := 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 ];

 

  Min := 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;

 

  MaxIndex := 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;

 

  MinIndex := 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;

 

  Cond := 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;

 

  Norm := 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

 

           SplineIndex := 1;

 

           if X <= d1.E( 1 )^ then

              exit;

 

           SplineIndex := NPoint - 1;

 

           if X >= d1.E( NPoint )^ then

              exit;

 

           IndRecurs( 1, NPoint );

           SplineIndex := 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 );

 

        F := 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 );

     DF := 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 ) );

 

        DF := 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

  Addr := VAddr.Ptr;

  //Addr := VirtualToPtr( VAddr );

end;

 

function Matrix.Addr : MatrixVPrimPtr;

begin

  Addr := 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

     E := 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

     E := 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 ];

 

  Summ := 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 ] );

 

  SummSqr := 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

  Addr := 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 ) ) );

  Alpha := 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.

 

 

Этот модуль используется почти во всех реализованных мной численных алгоритмах и методах. Те части, которые писал не я – приводятся без изменений(по возможности) стиля и комментариев.

Code:

unit HNVString;

 

interface

 

{*********************************************************}

{*                  OPSTRING.PAS 1.20                    *}

{*     Copyright (c) TurboPower Software 1987, 1992.     *}

{*                 All rights reserved.                  *}

{*********************************************************}

 

Type

  CharSet =

     Set of Char;

 

function WordCount(S : string; WordDelims : CharSet) : Byte;

   {-Given a set of word delimiters, return number of words in S}

 

function WordPosition(N : Byte; S : string; WordDelims : CharSet) : Byte;

   {-Given a set of word delimiters, return start position of N'th word in S}

 

function ExtractWord(N : Byte; S : string; WordDelims : CharSet) : string;

   {-Given a set of word delimiters, return the N'th word in S}

 

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.

©Drkb::03965

http://delphiworld.narod.ru/

DelphiWorld 6.0