d6_matrices.pas

ADaM is a data mining and image processing toolkit

procedure LU_Solve(A : TCompMatrix; B : TCompVector;
                   Lbound, Ubound : Integer; X : TCompVector); overload;
{ ----------------------------------------------------------------------
  LU solution for complex matrices
  ---------------------------------------------------------------------- }

function SV_Decomp(A : TMatrix; Lbound, Ubound1, Ubound2 : Integer;
                   S : TVector; V : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Singular value decomposition. Factors the matrix A (n x m, with n >= m)
  as a product U * S * V' where U is a (n x m) column-orthogonal matrix,
  S a (m x m) diagonal matrix with elements >= 0 (the singular values)
  and V a (m x m) orthogonal matrix. This routine is used in conjunction
  with SV_Solve to solve a system of equations.
  ----------------------------------------------------------------------
  Input parameters : A       = matrix
                     Lbound  = index of first matrix element
                     Ubound1 = index of last matrix element in 1st dim.
                     Ubound2 = index of last matrix element in 2nd dim.
  ----------------------------------------------------------------------
  Output parameter : A       = contains the elements of U
                     S       = vector of singular values
                     V       = orthogonal matrix
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_NON_CONV
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

procedure SV_SetZero(S : TVector; Lbound, Ubound : Integer; Tol : Float);
{ ----------------------------------------------------------------------
  Sets the singular values to zero if they are lower than a specified
  threshold.
  ----------------------------------------------------------------------
  Input parameters : S      = vector of singular values
                     Tol    = relative tolerance
                              Threshold value will be Tol * Max(S)
                     Lbound = index of first vector element
                     Ubound = index of last vector element
  ----------------------------------------------------------------------
  Output parameter : S      = modified singular values
  ---------------------------------------------------------------------- }

procedure SV_Solve(U : TMatrix; S : TVector; V : TMatrix; B : TVector;
                   Lbound, Ubound1, Ubound2 : Integer;
                   X : TVector);
{ ----------------------------------------------------------------------
  Solves a system of equations by singular value decomposition, after
  the matrix has been transformed by SV_Decomp, and the lowest singular
  values have been set to zero by SV_SetZero.
  ----------------------------------------------------------------------
  Input parameters : U, S, V = vector and matrices from SV_Decomp
                     B       = constant vector
                     Lbound,
                     Ubound1,
                     Ubound2 = as in SV_Decomp
  ----------------------------------------------------------------------
  Output parameter : X       = solution vector
                             = V * Diag(1/s(i)) * U' * B, for s(i) <> 0
  ---------------------------------------------------------------------- }

procedure SV_Approx(U : TMatrix; S : TVector; V : TMatrix;
                    Lbound, Ubound1, Ubound2 : Integer;
                    A : TMatrix);
{ ----------------------------------------------------------------------
  Approximates a matrix A by the product USV', after the lowest singular
  values have been set to zero by SV_SetZero.
  ----------------------------------------------------------------------
  Input parameters : U, S, V = vector and matrices from SV_Decomp
                     Lbound,
                     Ubound1,
                     Ubound2 = as in SV_Decomp
  ----------------------------------------------------------------------
  Output parameter : A       = approximated matrix
  ---------------------------------------------------------------------- }

function QR_Decomp(A : TMatrix; Lbound, Ubound1, Ubound2 : Integer;
                   R : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  QR decomposition. Factors the matrix A (n x m, with n >= m) as a
  product Q * R where Q is a (n x m) column-orthogonal matrix, and R
  a (m x m) upper triangular matrix. This routine is used in conjunction
  with QR_Solve to solve a system of equations.
  ----------------------------------------------------------------------
  Input parameters : A       = matrix
                     Lbound  = index of first matrix element
                     Ubound1 = index of last matrix element in 1st dim.
                     Ubound2 = index of last matrix element in 2nd dim.
  ----------------------------------------------------------------------
  Output parameter : A       = contains the elements of Q
                     R       = upper triangular matrix
  ----------------------------------------------------------------------
  Possible results : MAT_OK
                     MAT_SING
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

procedure QR_Solve(Q, R : TMatrix; B : TVector;
                   Lbound, Ubound1, Ubound2 : Integer;
                   X : TVector);
{ ----------------------------------------------------------------------
  Solves a system of equations by the QR decomposition,
  after the matrix has been transformed by QR_Decomp.
  ----------------------------------------------------------------------
  Input parameters : Q, R    = matrices from QR_Decomp
                     B       = constant vector
                     Lbound,
                     Ubound1,
                     Ubound2 = as in QR_Decomp
  ----------------------------------------------------------------------
  Output parameter : X       = solution vector
  ---------------------------------------------------------------------- }

implementation

  procedure DimVector(var V : TVector; Ubound : Integer); overload;
  var
    I : Integer;
  begin
    { Check bounds }
    if Ubound < 0 then
      begin
        V := nil;
        Exit;
      end;

    { Allocate vector }
    SetLength(V, Succ(Ubound));
    if V = nil then Exit;

    { Initialize vector }
    for I := 0 to Ubound do
      V[I] := 0.0;
  end;

  procedure DimVector(var V : TIntVector; Ubound : Integer); overload;
  var
    I : Integer;
  begin
    { Check bounds }
    if Ubound < 0 then
      begin
        V := nil;
        Exit;
      end;

    { Allocate vector }
    SetLength(V, Succ(Ubound));
    if V = nil then Exit;

    { Initialize vector }
    for I := 0 to Ubound do
      V[I] := 0;
  end;

  procedure DimVector(var V : TCompVector; Ubound : Integer); overload;
  var
    I : Integer;
  begin
    { Check bounds }
    if Ubound < 0 then
      begin
        V := nil;
        Exit;
      end;

    { Allocate vector }
    SetLength(V, Succ(Ubound));
    if V = nil then Exit;

    { Initialize vector }
    for I := 0 to Ubound do
      V[I] := C_zero;
  end;

  procedure DimVector(var V : TBoolVector; Ubound : Integer); overload;
  var
    I : Integer;
  begin
    { Check bounds }
    if Ubound < 0 then
      begin
        V := nil;
        Exit;
      end;

    { Allocate vector }
    SetLength(V, Succ(Ubound));
    if V = nil then Exit;

    { Initialize vector }
    for I := 0 to Ubound do
      V[I] := False;
  end;

  procedure DimVector(var V : TStrVector; Ubound : Integer); overload;
  var
    I : Integer;
  begin
    { Check bounds }
    if Ubound < 0 then
      begin
        V := nil;
        Exit;
      end;

    { Allocate vector }
    SetLength(V, Succ(Ubound));
    if V = nil then Exit;

    { Initialize vector }
    for I := 0 to Ubound do
      V[I] := '';
  end;

  procedure DimMatrix(var A : TMatrix; Ubound1, Ubound2 : Integer); overload;
  var
    I, J : Integer;
  begin
    if (Ubound1 < 0) or (Ubound2 < 0) then
      begin
        A := nil;
        Exit;
      end;

    { Allocate matrix }
    SetLength(A, Succ(Ubound1), Succ(Ubound2));
    if A = nil then Exit;

    { Initialize matrix }
    for I := 0 to Ubound1 do
      for J := 0 to Ubound2 do
        A[I,J] := 0.0;
  end;

  procedure DimMatrix(var A : TIntMatrix; Ubound1, Ubound2 : Integer); overload;
  var
    I, J : Integer;
  begin
    { Check bounds }
    if (Ubound1 < 0) or (Ubound2 < 0) then
      begin
        A := nil;
        Exit;
      end;

    { Allocate matrix }
    SetLength(A, Succ(Ubound1), Succ(Ubound2));
    if A = nil then Exit;

    { Initialize matrix }
    for I := 0 to Ubound1 do
      for J := 0 to Ubound2 do
        A[I,J] := 0;
  end;

  procedure DimMatrix(var A : TCompMatrix; Ubound1, Ubound2 : Integer); overload;
  var
    I, J : Integer;
  begin
    { Check bounds }
    if (Ubound1 < 0) or (Ubound2 < 0) then
         begin
           A := nil;
           Exit;
         end;

    { Allocate matrix }
    SetLength(A, Succ(Ubound1), Succ(Ubound2));
    if A = nil then Exit;

    { Initialize matrix }
    for I := 0 to Ubound1 do
      for J := 0 to Ubound2 do
        A[I,J] := C_zero;
  end;

  procedure DimMatrix(var A : TBoolMatrix; Ubound1, Ubound2 : Integer); overload;
  var
    I, J : Integer;
  begin
    { Check bounds }
    if (Ubound1 < 0) or (Ubound2 < 0) then
      begin
        A := nil;
        Exit;
      end;

    { Allocate matrix }
    SetLength(A, Succ(Ubound1), Succ(Ubound2));
    if A = nil then Exit;

    { Initialize matrix }
    for I := 0 to Ubound1 do
      for J := 0 to Ubound2 do
        A[I,J] := False;
  end;

  procedure DimMatrix(var A : TStrMatrix; Ubound1, Ubound2 : Integer); overload;
  var
    I, J : Integer;
  begin
    { Check bounds }
    if (Ubound1 < 0) or (Ubound2 < 0) then
      begin
        A := nil;
        Exit;
      end;

    { Allocate matrix }
    SetLength(A, Succ(Ubound1), Succ(Ubound2));
    if A = nil then Exit;

    { Initialize matrix }
    for I := 0 to Ubound1 do
      for J := 0 to Ubound2 do
        A[I,J] := '';
  end;

  procedure SwapRows(I, K : Integer; A : TMatrix; Lbound, Ubound : Integer);
  var
    J : Integer;
  begin
    for J := Lbound to Ubound do
      Swap(A[I,J], A[K,J]);
  end;

  procedure SwapCols(J, K : Integer; A : TMatrix; Lbound, Ubound : Integer);
  var
    I : Integer;
  begin
    for I := Lbound to Ubound do
      Swap(A[I,J], A[I,K]);
  end;

  procedure CopyVector(Dest, Source : TVector; Lbound, Ubound : Integer);
  var
    I : Integer;
  begin
    for I := Lbound to Ubound do
      Dest[I] := Source[I];
  end;

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

  procedure CopyRowFromVector(Dest : TMatrix; Source : TVector;
                              Lbound, Ubound, Row : Integer);
  var
    J : Integer;
  begin
    for J := Lbound to Ubound do
      Dest[Row,J] := Source[J];
  end;

  procedure CopyColFromVector(Dest : TMatrix; Source : TVector;
                              Lbound, Ubound, Col : Integer);
  var
    I : Integer;
  begin
    for I := Lbound to Ubound do
      Dest[I,Col] := Source[I];
  end;

  procedure CopyVectorFromRow(Dest : TVector; Source : TMatrix;
                              Lbound, Ubound, Row : Integer);
  var
    J : Integer;
  begin
    for J := Lbound to Ubound do
      Dest[J] := Source[Row,J];
  end;

  procedure CopyVectorFromCol(Dest : TVector; Source : TMatrix;
                              Lbound, Ubound, Col : Integer);
  var
    I : Integer;
  begin
    for I := Lbound to Ubound do
      Dest[I] := Source[I,Col];
  end;

  function Min(X : TVector; Lbound, Ubound : Integer) : Float; overload;
  var
    Xmin : Float;
    I    : Integer;