d6_matrices.pas

ADaM is a data mining and image processing toolkit

  function InvDet(A              : TMatrix;
                  Lbound, Ubound : Integer;
                  A_inv          : TMatrix;
                  var Det        : Float) : 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);

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

    InvDet := MAT_OK;
  end;

  function Det(A : TMatrix; Lbound, Ubound : Integer) : Float;
  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 }
    D          : Float;       { Determinant }
  begin
    DimVector(PRow, Ubound);
    DimVector(PCol, Ubound);
    DimVector(MCol, Ubound);

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

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

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

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

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

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

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

        { Transform pivot row }
        for J := Lbound to Ubound do
          A[K,J] := A[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[I,J] := A[I,J] - T * A[K,J];
            end;
        Inc(K);
      end;

    Det := D;
  end;

  function Cholesky(A : TMatrix; Lbound, Ubound : Integer;
                    L : TMatrix) : Integer;
  var
    I, J, K : Integer;
    Sum     : Float;
  begin
    for K := Lbound to Ubound do
      begin
        Sum := A[K,K];
        for J := Lbound to Pred(K) do
          Sum := Sum - Sqr(L[K,J]);

        if Sum <= 0.0 then
          begin
            Cholesky := MAT_NOT_PD;
            Exit;
          end;

        L[K,K] := Sqrt(Sum);
        for I := Succ(K) to Ubound do
          begin
            Sum := A[I,K];
            for J := Lbound to Pred(K) do
              Sum := Sum - L[I,J] * L[K,J];
            L[I,K] := Sum / L[K,K];
          end;
      end;
    Cholesky := MAT_OK;
  end;

var                    { Used by LU procedures }
  Index : TIntVector;  { Records the row permutations }

  function LU_Decomp(A : TMatrix; Lbound, Ubound : Integer) : Integer; overload;
  var
    I, Imax, J, K : Integer;
    Pvt, T, Sum   : Float;
    V             : TVector;
  begin
    DimVector(V, Ubound);
    DimVector(Index, Ubound);

    for I := Lbound to Ubound do
      begin
        Pvt := 0.0;
        for J := Lbound to Ubound do
          if Abs(A[I,J]) > Pvt then
            Pvt := Abs(A[I,J]);
        if Pvt < MACHEP then
          begin
            LU_Decomp := MAT_SINGUL;
            Exit;
          end;
        V[I] := 1.0 / Pvt;
      end;
    for J := Lbound to Ubound do
      begin
        for I := Lbound to Pred(J) do
          begin
            Sum := A[I,J];
            for K := Lbound to Pred(I) do
              Sum := Sum - A[I,K] * A[K,J];
            A[I,J] := Sum;
          end;
        Imax := 0;
        Pvt := 0.0;
        for I := J to Ubound do
          begin
            Sum := A[I,J];
            for K := Lbound to Pred(J) do
              Sum := Sum - A[I,K] * A[K,J];
            A[I,J] := Sum;
            T := V[I] * Abs(Sum);
            if T > Pvt then
              begin
                Pvt := T;
                Imax := I;
              end;
          end;
        if J <> Imax then
          begin
            SwapRows(Imax, J, A, Lbound, Ubound);
            V[Imax] := V[J];
          end;
        Index[J] := Imax;
        if A[J,J] = 0.0 then
          A[J,J] := MACHEP;
        if J <> Ubound then
          begin
            T := 1.0 / A[J,J];
            for I := Succ(J) to Ubound do
              A[I,J] := A[I,J] * T;
          end;
      end;
    LU_Decomp := MAT_OK;
  end;

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

  function LU_Decomp(A : TCompMatrix; Lbound, Ubound : Integer) : Integer; overload;
  var
    I, Imax, J, K : Integer;
    C, Pvt, T     : Float;
    Sum           : Complex;
    V             : TVector;
  begin
    DimVector(V, Ubound);
    DimVector(Index, Ubound);

    for I := Lbound to Ubound do
      begin
        Pvt := 0.0;
        for J := Lbound to Ubound do
          begin
            C := CAbs(A[I,J]);
            if C > Pvt then Pvt := C;
          end;
        if Pvt < MACHEP then
          begin
            LU_Decomp := MAT_SINGUL;
            Exit;
          end;
        V[I] := 1.0 / Pvt;
      end;
    for J := Lbound to Ubound do
      begin
        for I := Lbound to Pred(J) do
          begin
            Sum := A[I,J];
            for K := Lbound to Pred(I) do
              { Sum := Sum - A[I,K] * A[K,J]; }
              Sum := CSub(Sum, CMult(A[I,K], A[K,J]));
            A[I,J] := Sum;
          end;
        Imax := 0;
        Pvt := 0.0;
        for I := J to Ubound do
          begin
            Sum := A[I,J];
            for K := Lbound to Pred(J) do
              { Sum := Sum - A[I,K] * A[K,J]; }
              Sum := CSub(Sum, CMult(A[I,K], A[K,J]));
            A[I,J] := Sum;
            T := V[I] * CAbs(Sum);
            if T > Pvt then
              begin
                Pvt := T;
                Imax := I;
              end;
          end;
        if J <> Imax then
          begin
            { SwapRows(Imax, J, A, Lbound, Ubound); }
            for K := Lbound to Ubound do
              Swap(A[Imax,K], A[J,K]);
            V[Imax] := V[J];
          end;
        Index[J] := Imax;
        if CAbs(A[J,J]) = 0.0 then
          A[J,J] := Cmplx(MACHEP, MACHEP);
        if J <> Ubound then
          for I := Succ(J) to Ubound do
            { A[I,J] := A[I,J] / A[J,J]; }
            A[I,J] := CDiv(A[I,J], A[J,J]);
      end;
    LU_Decomp := MAT_OK;
  end;

  procedure LU_Solve(A : TCompMatrix; B : TCompVector;
                     Lbound, Ubound : Integer; X : TCompVector); overload;
  var
    I, Ip, J, K : Integer;
    Sum         : Complex;
  begin