d6_matrices.pas

ADaM is a data mining and image processing toolkit

  begin
    Xmin := X[Lbound];
    for I := Succ(Lbound) to Ubound do
      if X[I] < Xmin then Xmin := X[I];
    Min := Xmin;
  end;

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

  function Min(X : TIntVector; Lbound, Ubound : Integer) : Integer; overload;
  var
    I, Xmin : Integer;
  begin
    Xmin := X[Lbound];
    for I := Succ(Lbound) to Ubound do
      if X[I] < Xmin then Xmin := X[I];
    Min := Xmin;
  end;

  function Max(X : TIntVector; Lbound, Ubound : Integer) : Integer; overload;
  var
    I, Xmax : Integer;
  begin
    Xmax := X[Lbound];
    for I := Succ(Lbound) to Ubound do
      if X[I] > Xmax then Xmax := X[I];
    Max := Xmax;
  end;

  procedure Transpose(A : TIntMatrix;
                      Lbound1, Lbound2, Ubound1, Ubound2 : Integer;
                      A_t : TIntMatrix); overload;
  var
    I, J : Integer;
  begin
    for I := Lbound1 to Ubound1 do
      for J := Lbound2 to Ubound2 do
        A_t[J,I] := A[I,J];
  end;

  procedure Transpose(A : TMatrix;
                      Lbound1, Lbound2, Ubound1, Ubound2 : Integer;
                      A_t : TMatrix); overload;
  var
    I, J : Integer;
  begin
    for I := Lbound1 to Ubound1 do
      for J := Lbound2 to Ubound2 do
        A_t[J,I] := A[I,J];
  end;

  function GaussJordan(A              : TMatrix;
                       B              : TVector;
                       Lbound, Ubound : Integer;
                       A_inv          : TMatrix;
                       X              : TVector;
                       var Det        : Float) : Integer; overload;
  var
    I, J, K    : Integer;     { Loop variables }
    Ik, Jk     : Integer;     { Pivot coordinates }
    Pvt        : Float;       { Pivot }
    T          : Float;       { Auxiliary variable }
    PRow, PCol : TIntVector;  { Store line and column of pivot }
    MCol       : TVector;     { Stores a column of the matrix }
  begin
    DimVector(PRow, Ubound);
    DimVector(PCol, Ubound);
    DimVector(MCol, Ubound);

    { Copy A into A_inv and B into X }
    CopyMatrix(A_inv, A, Lbound, Lbound, Ubound, Ubound);
    CopyVector(X, B, Lbound, Ubound);

    Det := 1.0;
    K := Lbound;
    while K <= Ubound do
      begin
        { Search for largest pivot in submatrix A_inv[K..Ubound, K..Ubound] }
        Pvt := A_inv[K,K];
        Ik := K;
        Jk := K;
        for I := K to Ubound do
          for J := K to Ubound do
            if Abs(A_inv[I,J]) > Abs(Pvt) then
              begin
                Pvt := A_inv[I,J];
                Ik := I;
                Jk := J;
              end;

        { Pivot too small ==> quasi-singular matrix }
        if Abs(Pvt) < MACHEP then
          begin
            GaussJordan := MAT_SINGUL;
            Exit;
          end;

        { Save pivot position }
        PRow[K] := Ik;
        PCol[K] := Jk;

        { Update determinant }
        Det := Det * Pvt;
        if Ik <> K then Det := - Det;
        if Jk <> K then Det := - Det;

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

        { Exchange current column (K) with pivot column (Jk) }
        if Jk <> K then
          SwapCols(Jk, K, A_inv, Lbound, Ubound);

        { Store col. K of A_inv into MCol and set this col. to 0 }
        for I := Lbound to Ubound do
          if I <> K then
            begin
              MCol[I] := A_inv[I, K];
              A_inv[I, K] := 0.0;
            end
          else
            begin
              MCol[I] := 0.0;
              A_inv[I, K] := 1.0;
            end;

        { Transform pivot row }
        for J := Lbound to Ubound do
          A_inv[K,J] := A_inv[K,J] / Pvt;
        X[K] := X[K] / Pvt;

        { Transform other rows }
        for I := Lbound to Ubound do
          if I <> K then
            begin
              T := MCol[I];
              for J := Lbound to Ubound do
                A_inv[I,J] := A_inv[I,J] - T * A_inv[K,J];
              X[I] := X[I] - T * X[K];
            end;
        Inc(K);
      end;

    { Exchange rows of inverse matrix and solution vector }
    for I := Ubound downto Lbound do
      begin
        Ik := PCol[I];
        if Ik <> I then
          begin
            SwapRows(Ik, I, A_inv, Lbound, Ubound);
            Swap(X[Ik], X[I]);
          end;
      end;

    { Exchange columns of inverse matrix }
    for J := Ubound downto Lbound do
      begin
        Jk := PRow[J];
        if Jk <> J then
          SwapCols(Jk, J, A_inv, Lbound, Ubound);
      end;

    GaussJordan := MAT_OK;
  end;

  function GaussJordan(A, B                     : TMatrix;
                       Lbound, Ubound1, Ubound2 : Integer;
                       A_inv, X                 : TMatrix;
                       var Det                  : Float) : Integer; overload;
  var
    I, J, K    : Integer;     { Loop variables }
    Ik, Jk     : Integer;     { Pivot coordinates }
    Pvt        : Float;       { Pivot }
    T          : Float;       { Auxiliary variable }
    PRow, PCol : TIntVector;  { Store line and column of pivot }
    MCol       : TVector;     { Stores a column of the matrix }
  begin
    DimVector(PRow, Ubound1);
    DimVector(PCol, Ubound1);
    DimVector(MCol, Ubound1);

    { Copy A into A_inv and B into X }
    CopyMatrix(A_inv, A, Lbound, Lbound, Ubound1, Ubound1);
    CopyMatrix(X, B, Lbound, Lbound, Ubound1, Ubound2);

    Det := 1.0;
    K := Lbound;
    while K <= Ubound1 do
      begin
        { Search for largest pivot in submatrix A_inv[K..Ubound1, K..Ubound1] }
        Pvt := A_inv[K,K];
        Ik := K;
        Jk := K;
        for I := K to Ubound1 do
          for J := K to Ubound1 do
            if Abs(A_inv[I,J]) > Abs(Pvt) then
              begin
                Pvt := A_inv[I,J];
                Ik := I;
                Jk := J;
              end;

        { Pivot too small ==> quasi-singular matrix }
        if Abs(Pvt) < MACHEP then
          begin
            GaussJordan := MAT_SINGUL;
            Exit;
          end;

        { Save pivot position }
        PRow[K] := Ik;
        PCol[K] := Jk;

        { Update determinant }
        Det := Det * Pvt;
        if Ik <> K then Det := - Det;
        if Jk <> K then Det := - Det;

        { Exchange current row (K) with pivot row (Ik) }
        if Ik <> K then
          begin
            SwapRows(Ik, K, A_inv, Lbound, Ubound1);
            SwapRows(Ik, K, X, Lbound, Ubound2);
          end;

        { Exchange current column (K) with pivot column (Jk) }
        if Jk <> K then
          SwapCols(Jk, K, A_inv, Lbound, Ubound1);

        { Store col. K of A_inv into MCol and set this col. to 0 }
        for I := Lbound to Ubound1 do
          if I <> K then
            begin
              MCol[I] := A_inv[I, K];
              A_inv[I, K] := 0.0;
            end
          else
            begin
              MCol[I] := 0.0;
              A_inv[I, K] := 1.0;
            end;

        { Transform pivot row }
        for J := Lbound to Ubound1 do
          A_inv[K,J] := A_inv[K,J] / Pvt;
        for J := Lbound to Ubound2 do
          X[K,J] := X[K,J] / Pvt;

        { Transform other rows }
        for I := Lbound to Ubound1 do
          if I <> K then
            begin
              T := MCol[I];
              for J := Lbound to Ubound1 do
                A_inv[I,J] := A_inv[I,J] - T * A_inv[K,J];
              for J := Lbound to Ubound2 do
                X[I,J] := X[I,J] - T * X[K,J];
            end;

        Inc(K);
      end;

    { Exchange lines of inverse and solution matrices }
    for I := Ubound1 downto Lbound do
      begin
        Ik := PCol[I];
        if Ik <> I then
          begin
            SwapRows(Ik, I, A_inv, Lbound, Ubound1);
            SwapRows(Ik, I, X, Lbound, Ubound2);
          end;
      end;

    { Exchange columns of inverse matrix }
    for J := Ubound1 downto Lbound do
      begin
        Jk := PRow[J];
        if Jk <> J then
          SwapCols(Jk, J, A_inv, Lbound, Ubound1);
      end;

    GaussJordan := MAT_OK;
  end;

  function InvMat(A              : TMatrix;
                  Lbound, Ubound : Integer;
                  A_inv          : TMatrix) : Integer;
  var
    I, J, K    : Integer;     { Loop variables }
    Ik, Jk     : Integer;     { Pivot coordinates }
    Pvt        : Float;       { Pivot }
    T          : Float;       { Auxiliary variable }
    PRow, PCol : TIntVector;  { Store line and column of pivot }
    MCol       : TVector;     { Stores a column of the matrix }
  begin
    DimVector(PRow, Ubound);
    DimVector(PCol, Ubound);
    DimVector(MCol, Ubound);

    { Copy A into A_inv }
    CopyMatrix(A_inv, A, Lbound, Lbound, Ubound, Ubound);

    K := Lbound;
    while K <= Ubound do
      begin
        { Search for largest pivot in submatrix A_inv[K..Ubound, K..Ubound] }
        Pvt := A_inv[K,K];
        Ik := K;
        Jk := K;
        for I := K to Ubound do
          for J := K to Ubound do
            if Abs(A_inv[I,J]) > Abs(Pvt) then
              begin
                Pvt := A_inv[I,J];
                Ik := I;
                Jk := J;
              end;

        { Pivot too small ==> quasi-singular matrix }
        if Abs(Pvt) < MACHEP then
          begin
            InvMat := MAT_SINGUL;
            Exit;
          end;

        { Save pivot position }
        PRow[K] := Ik;
        PCol[K] := Jk;

        { Exchange current row (K) with pivot row (Ik) }
        if Ik <> K then
          SwapRows(Ik, K, A_inv, Lbound, Ubound);

        { Exchange current column (K) with pivot column (Jk) }
        if Jk <> K then
          SwapCols(Jk, K, A_inv, Lbound, Ubound);

        { Store col. K of A_inv into MCol and set this col. to 0 }
        for I := Lbound to Ubound do
          if I <> K then
            begin
              MCol[I] := A_inv[I, K];
              A_inv[I, K] := 0.0;
            end
          else
            begin
              MCol[I] := 0.0;
              A_inv[I, K] := 1.0;
            end;

        { Transform pivot row }
        for J := Lbound to Ubound do
          A_inv[K,J] := A_inv[K,J] / Pvt;

        { Transform other rows }
        for I := Lbound to Ubound do
          if I <> K then
            begin
              T := MCol[I];
              for J := Lbound to Ubound do
                A_inv[I,J] := A_inv[I,J] - T * A_inv[K,J];
            end;
        Inc(K);
      end;

    { Exchange lines of inverse matrix }
    for I := Ubound downto Lbound do
      begin
        Ik := PCol[I];
        if Ik <> I then
          SwapRows(Ik, I, A_inv, Lbound, Ubound);
      end;

    { Exchange columns of inverse matrix }
    for J := Ubound downto Lbound do
      begin
        Jk := PRow[J];
        if Jk <> J then
          SwapCols(Jk, J, A_inv, Lbound, Ubound);
      end;

    InvMat := MAT_OK;
  end;