linalg.pas

Delphi math processing compononets and sources. Release.

  Signal := 0;

  if  N < 1  then
    begin
      Signal := -1;
      goto 3;
    end;
  if  N = 1  then
    begin
      T[1]^[1] := 1.0;
      Eigen[1] := A[1]^[1];
      goto 3;
    end;
  for i:= 1 to N  do
    begin
      for j := 1 to N  do
        T[i]^[j] := 0.0;
      T[i]^[i] := 1.0;
      Eigen[i] := A[i]^[i];
    end;

  Sigma1 := 0.0;
  OffDsq := 0.0;
  for i := 1 to N  do
    begin
      Sigma1 := Sigma1 + sqr(A[i]^[i]);
      if i < N  then
        for j := i + 1 to N  do
          OffDsq := OffDsq+sqr(A[i]^[j]);
    end;

  if  OffDsq < sqr(Eps2) then goto 2;
  S := OffDsq*2.0+Sigma1;
  Iter := 1;

  HoldIK := 1.0;
  while (Iter <= ItMax) and (HoldIK > Eps3)  do
    begin
      for i := 1 to N - 1  do
        for j:=i + 1 to N  do
          begin
            Q := abs(A[i]^[i]-A[j]^[j]);
            if Q > Eps1 then
              begin
                if abs(A[i]^[j]) <= Eps2  then  goto 1;
                P := 2.0*A[i]^[j]*(Q/(A[i]^[i]-A[j]^[j]));
                SPQ := sqrt(sqr(P)+sqr(Q));
                CSA := sqrt((1.0+Q/SPQ)/2.0);
                SNA := P/(SPQ*CSA*2.0);
              end       else
              begin
                CSA := sqrt(0.5);
                SNA := CSA;
              end;
            for k := 1 to N  do
              begin
                HoldKI := T[k]^[i];
                T[k]^[i] := HoldKI*CSA + T[k]^[j]*SNA;
                T[k]^[j] := HoldKI*SNA - T[k]^[j]*CSA;
              end;

            for k := i to N  do
              if k <= j then
                begin
                  Aik[k] := A[i]^[k];
                  A[i]^[k] := CSA*Aik[k] + SNA*A[k]^[j];
                  if k = j then
                    begin
                      A[j]^[k] := SNA*Aik[k] - CSA*A[j]^[k];
                      Aik[j]   := SNA*Aik[i] - CSA*Aik[j];
                    end;
                end     else
                begin
                  HoldIK := A[i]^[k];
                  A[i]^[k] := CSA*HoldIK + SNA*A[j]^[k];
                  A[j]^[k] := SNA*HoldIK - CSA*A[j]^[k];
                end;

            for k := 1 to j  do
              if k > i then
                A[k]^[j] := SNA*Aik[k] - CSA*A[k]^[j]
                       else
                begin
                  HoldKI := A[k]^[i];
                  A[k]^[i] := CSA*HoldKI + SNA*A[k]^[j];
                  A[k]^[j] := SNA*HoldKI - CSA*A[k]^[j];
                end;
          1:
          end; { of i, j - loop }

      Sigma2 := 0.0;
      for i := 1 to N do
        begin
          Eigen[i] := A[i]^[i];
          Sigma2 := Sigma2 + sqr(Eigen[i]);
        end;

      HoldIK := abs(1.0-Sigma1/Sigma2);
      Sigma1 := Sigma2;
      inc(Iter);
    end; { of while }

  if Iter > ItMax then Signal := ItMax;
{ sorting of eigenvalues and eigenvectors in increasing order }
2:for i := 1 to N  do
    begin
      k := i;
      for j := i to N  do
        if Eigen[j]<Eigen[k]  then  k := j;
      if k<>i  then
        begin
          P := Eigen[k];
          Eigen[k] := Eigen[i];
          Eigen[i] := P;
          for j := 1 to N  do
            begin
              P := T[j]^[k];
              T[j]^[k] := T[j]^[i];
              T[j]^[i] := P;
            end;
        end;
    end;
3:
end;{ Of Jacobi }

procedure JacobiC(N: integer; const A, B, U, V: Matrix; var E: Vector;
                    var Signal: integer);
const
  Eps1 = 6.0e-12;    Eps2  = 1.0e-20;
  Eps3 = 1.0e-10;    ItMax = 999;

var
  l,       m,
  i,       j,       k,       Iter: integer;
  Hii,     Hjj,     Aij,     Aji,
  Bij,     Bji,     Hii_Hjj, HijxHji,
  Norm,    OffDsq,  C,       CC,
  SS,      S_R,     S_I,     T2,
  Bki,     Aki,     Bik,     Aik,
  Uki,     Vki,     Ujk,     Vjk,
  signum,  Ek             : RealType;

label
  1, 2, 3;

begin{of JacobiC }

  Signal := 0;

  if N < 1  then
    begin
      Signal := -1;
      goto 3;
    end;

  if N = 1  then
    begin
      U[1]^[1] := 1.0;
      V[1]^[1] := 0.0;
      E[1]     := A[1]^[1];
      goto 3;
    end;

  for i:= 1 to N  do
  begin
    for j := 1 to N  do
    begin
      U[i]^[j] := 0.0;
      V[i]^[j] := 0.0;
    end;
    U[i]^[i] := 1.0;
    V[i]^[i] := 0.0;
    E[i]     := A[i]^[i];
  end;

  Norm   := 0.0;
  OffDsq := 0.0;
{
  Norm   is the L2 norm of the matrix A + i*B
  OffDsq is the L2 norm of the off-diagonal elements
}
  for i := 1 to N  do
  begin
    Norm := Norm + sqr(A[i]^[i]);
    for j := 1 to N  do
    if j <> i then OffDsq := OffDsq + sqr(A[i]^[j]) + sqr(B[i]^[j]);
  end;
  Norm := Norm + OffDsq;

  if OffDsq < Eps3*Norm then goto 2;

  Iter := 1;

  while (Iter <= ItMax) and (OffDsq > Eps3*Norm)  do
  begin{of iteration loop}
    for i := 1 to N - 1  do
      for j:=i + 1 to N  do
      begin
        Aji     := A[j]^[i];
        Bji     := B[j]^[i];
        HijxHji := sqr(Aji) + sqr(Bji);
        if HijxHji < Eps2  then  goto 1;
          Hii     := A[i]^[i];
          Hjj     := A[j]^[j];
          Hii_Hjj := A[i]^[i]-A[j]^[j];
          Aij     := A[i]^[j];
          Bij     := B[i]^[j];
          if abs(Hii_Hjj) > Eps1 then
          begin
            T2  := 4.0*HijxHji/sqr(Hii_Hjj);
            S_R := 2.0*Aji/(Hii_Hjj)/
                            sqrt(2.0*(1.0 + T2 + sqrt(1.0 + T2)));
            S_I := 2.0*Bji/(Hii_Hjj)/
                            sqrt(2.0*(1.0 + T2 + sqrt(1.0 + T2)));
            SS  := sqr(S_R) + sqr(S_I);
            CC  := 1.0 - SS;
            C   := sqrt(abs(CC));
          end
          else{ degenerate case }
          begin
            CC  := 0.5;
            SS  := CC;
            C   := sqrt(CC);
            if Hii_Hjj < 0.0 then signum := -1.0 else
                                  signum :=  1.0;
            S_R := signum*Aji/sqrt(2.0*HijxHji);
            S_I := signum*Bji/sqrt(2.0*HijxHji);
          end;
          for k := 1 to N  do
          begin
            Uki := U[k]^[i];
            Vki := V[k]^[i];
            U[k]^[i] := Uki*C      + U[k]^[j]*S_R - V[k]^[j]*S_I;
            V[k]^[i] := Vki*C      + U[k]^[j]*S_I + V[k]^[j]*S_R;
            U[k]^[j] := U[k]^[j]*C - Uki*S_R      - Vki*S_I;
            V[k]^[j] := V[k]^[j]*C + Uki*S_I      - Vki*S_R;
          end;
          for k := 1 to N  do
          if (k <> i) and (k <> j) then
          begin
            Aki := A[k]^[i];
            Bki := B[k]^[i];
            A[k]^[i] := Aki*C      + A[k]^[j]*S_R - B[k]^[j]*S_I;
            B[k]^[i] := Bki*C      + A[k]^[j]*S_I + B[k]^[j]*S_R;
            A[k]^[j] := A[k]^[j]*C - Aki*S_R      - Bki*S_I;
            B[k]^[j] := B[k]^[j]*C + Aki*S_I      - Bki*S_R;
          end;
          for k := 1 to N  do
          if (k <> i) and (k <> j) then
          begin
            Aik := A[i]^[k];
            Bik := B[i]^[k];
            A[i]^[k] := Aik*C      + A[j]^[k]*S_R + B[j]^[k]*S_I;
            B[i]^[k] := Bik*C      - A[j]^[k]*S_I + B[j]^[k]*S_R;
            A[j]^[k] := A[j]^[k]*C - Aik*S_R      + Bik*S_I;
            B[j]^[k] := B[j]^[k]*C - Aik*S_I      - Bik*S_R;
          end;
          A[i]^[i] := Hii*CC + Hjj*SS + 2.0*C*(S_R*Aji + S_I*Bji);
          A[j]^[j] := Hjj*CC + Hii*SS - 2.0*C*(S_R*Aji + S_I*Bji);
{             _ _
 SS -> S*S or S*S
          A[j]^[i] := Aji*CC - Aij*SS + C*S_R*(Hjj - Hii);
          B[j]^[i] := Bji*CC - Bij*SS + C*S_I*(Hjj - Hii);
          A[i]^[j] := Aij*CC - Aji*SS + C*S_R*(Hjj - Hii);
          B[i]^[j] := Bij*CC - Bji*SS - C*S_I*(Hjj - Hii);
}
          A[j]^[i] := 0.0;
          B[j]^[i] := 0.0;
          A[i]^[j] := 0.0;
          B[i]^[j] := 0.0;

        1:
        end; { of i, j - loop }
  OffDsq := 0.0;
  for i := 1 to N  do
  begin
    E[i] := A[i]^[i];
    for j := 1 to N  do
    if j <> i then OffDsq := OffDsq + sqr(A[i]^[j]) + sqr(B[i]^[j]);
  end;
  inc(Iter);
  end; { of while iteration loop}

  if Iter > ItMax then Signal := ItMax;
{
 sorting of eigenvalues and eigenvectors in increasing order
}
2:for i := 1 to N  do
  begin
    k := i;
    for j := i to N  do
      if E[j] < E[k]  then  k := j;
    if k<>i  then
    begin
      Ek   := E[k];
      E[k] := E[i];
      E[i] := Ek;
      for j := 1 to N  do
      begin
        Ujk      := U[j]^[k];
        U[j]^[k] := U[j]^[i];
        U[j]^[i] := Ujk;
        Vjk      := V[j]^[k];
        V[j]^[k] := V[j]^[i];
        V[j]^[i] := Vjk;
      end;
    end;
  end;
3:
end; { Of JacobiC }

{ Fast Inversion of the matrix A by the Gauss-Jordan elimination algorithm }
procedure FastInverse(N: integer; var A: Matrix; var J0: IVector;
                         var Det: RealType; var Signal: integer);
const
  Eps = 1.0e-16;
var
  i,j,k,i0     : integer;
  A0,B         : RealType;
  Ai,Ai0       : VcPtr;

begin
  Det := 1.0;
  Signal := 0;
  if  N<=1  then
    begin
      if abs(A[1]^[1])<Eps  then
        begin
          Signal := 1;
          Exit;
        end;
      A[1]^[1] := 1.0/A[1]^[1];
      Det      := A[1]^[1];
      Exit;
    end;

  for i:=1 to N  do
    begin
      {  1. Finding of the pivot element(in A0)  ;
         i0  is the number of the pivot row    }
      A0 := 0.0;
      i0 := 0;
      for j:=i to N  do
        if abs(A[j]^[i])>A0  then
          begin
            A0 := abs(A[j]^[i]);
            i0 := j;
          end;
      if A0<Eps  then
        begin
          Signal := i;  { Degeneracy is found here }
          Exit;
        end;

      {  2. Swapping rows   }
      J0[i] := i0;
      B     := 1.0/A[i0]^[i];
      Ai    := A[i0];
      Ai0   := A[i];
      A[i]  := Ai;
      A[i0] := Ai0;
      for j:=1 to N  do
        Ai^[j] := Ai^[j]*B;
      Ai^[i] := B;

      {  3. Subtracting the rows   }
      for j:=1 to N  do
        if i<>j  then
          begin
            Ai0 := A[j];
            B   := Ai0^[i];
            Ai0^[i] := 0.0;
            for k:=1 to N  do
              Ai0^[k] := Ai0^[k] - B*Ai^[k];
          end;

      Det := Det/Ai^[i];

    end;

    {  4.  Back Swapping the columns    }
    for i:=N downto 1  do
      begin
        j := J0[i];
        if j<>i  then
          begin
            Det := -Det;
            for k:=1 to N  do
              begin
                B        := A[k]^[i];
                A[k]^[i] := A[k]^[j];
                A[k]^[j] := B
              end;
          end;
      end;
end;{   Of the procedure FastInverse  }

procedure SVD(NA,M,N: integer; const A,U,V:  Matrix; var W,RV1: Vector;
                 MatU, MatV:  boolean;  RetCode:  integer );
const
  ItnLimit = 30;
var
  i,j,k,l,i1,k1,l1,its,mn,ExitKey: integer;
  C,G,F,X,S,H,Y,Z,Scale,ANorm,GG: RealType;
begin

  RetCode := 0;

  if NA=M  then
    for i:=1 to M  do
      for j:=1 to N  do
        U[i]^[j] := A[i]^[j]
  else
    for i:=1 to M  do
      for j:=1 to N  do
        U[i]^[j] := A[j]^[i];

  G     :=  0.0;
  Scale :=  0.0;
  ANorm :=  0.0;

  for i:=1 to N  do
    begin
      l      := i+1;
      RV1[i] := Scale*G;
      G      := 0.0;
      S      := 0.0;
      Scale  := 0.0;
      if i<=M  then
        begin
          for k:=1 to M  do
            Scale := Scale+abs(U[k]^[i]);
          if Scale<>0.0  then
            begin
              for k:=i to M  do
                begin
                  U[k]^[i] := U[k]^[i]/Scale;
                  S        := S+sqr(U[k]^[i]);
                end;
              F := U[i]^[i];
              G := -Sign2(sqrt(S),F);
              H := F*G-S;
              U[i]^[i] := F-G;
              if i<>N  then
                for j:=l to N  do
                  begin
                    S := 0.0;
                    for k:=i to M  do
                      S := S+U[k]^[i]*U[k]^[j];
                    F := S/H;
                    for k:=i to M  do
                      U[k]^[j] := U[k]^[j]+F*U[k]^[i];
                  end;
              for k:=i to M  do
                U[k]^[i] := Scale*U[k]^[i];
            end;
        end;

      W[i]  := Scale*G;
      G     := 0.0;
      S     := 0.0;
      Scale := 0.0;

      if (i<=M) and (i<>N)  then
        begin
          for k:=l to N  do
            Scale := Scale+abs(U[i]^[k]);
          if Scale<>0.0  then
            begin
              for k:=l to N  do
                begin
                  U[i]^[k] := U[i]^[k]/Scale;
                  S        := S+sqr(U[i]^[k]);
                end;
              F := U[i]^[l];
              G := -Sign2(sqrt(S),F);
              H := F*G-S;
              U[i]^[l] := F-G;
              for k:=l to N  do
                RV1[k] := U[i]^[k]/H;
              if i<>M  then
                for j:=l to M  do
                  begin
                    S := 0.0;
                    for k:=l to N  do
                      S := S+U[j]^[k]*U[i]^[k];
                    for k:=l to N  do
                      U[j]^[k] := U[j]^[k]+S*RV1[k];
                  end;
              for k:=l to N  do
                U[i]^[k] := Scale*U[i]^[k];
            end;
        end;

      ANorm := AMax1(ANorm,abs(W[i])+abs(RV1[i]));

    end;{ of the i-loop}

  { Accumulation of the right-hand transformations   }
  if  MatV  then
    for i:=N downto 1  do
      begin
        if i<>N  then
          begin
            if G<>0.0  then
              begin
                for j:=l to N  do
                  V[j]^[i] := (U[i]^[j]/U[i]^[l]) / G;
                for j:=l to N  do
                  begin
                    S := 0.0;
                    for k:=l to N  do
                      S := S+U[i]^[k]*V[k]^[j];
                    for k:=l to N  do
                      V[k]^[j] := V[k]^[j]+S*V[k]^[i];
                  end;
              end;
            for j:=l to N  do
              begin
                V[i]^[j] := 0.0;
                V[j]^[i] := 0.0;
              end;
          end;

        V[i]^[i] := 1.0;
        G        := RV1[i];
        l        := i;

      end;

  { Accumulation of the left-hand transformations }
  if  MatU  then
    begin
      mn := N;
      if M<N  then  mn := M;