linalg.pas

Delphi math processing compononets and sources. Release.

      for i:=mn downto 1  do
        begin
          l := i+1;
          G := W[i];
          if i<>N  then
            for j:=l to N  do
              U[i]^[j] := 0.0;
          if G<>0.0  then
            begin
              if i<>mn  then
                for j:=l to N  do
                  begin
                    S := 0.0;
                    for k:=l to M  do
                      S := S+U[k]^[i]*U[k]^[j];
                    F := (S/U[i]^[i]) / G;
                    for k:=i to M  do
                      U[k]^[j] := U[k]^[j]+F*U[k]^[i];
                  end;
              for j:=i to M  do
                U[j]^[i] := U[j]^[i]/G;
            end      else
            for j:=i to M  do
              U[j]^[i] := 0.0;

          U[i]^[i] := U[i]^[i]+1.0;

        end;
    end;

  { Diagonalization of the two-diagonal form }
  for k:=N downto 1  do
    begin
      k1  := k-1;
      its := 0;

      repeat
        ExitKey  := 0;
        l        := k+1;
        while (ExitKey=0) and (l>1)  do
          begin
            dec(l);
            l1 := l-1;
            if  abs(RV1[l])+ANorm = ANorm   then   ExitKey:=1
            else  if  l1>0  then
              if abs(W[l1])+ANorm = ANorm   then   ExitKey:=2;
          end;

        if ExitKey<>1  then
          begin
            C       := 0.0;
            S       := 1.0;
            ExitKey := 0;
            i       := l;
            while (ExitKey=0) and (i<=k)  do
              begin
                F       :=  S*RV1[i];
                RV1[i]  :=  C*RV1[i];
                if abs(F)+ANorm = ANorm  then   ExitKey := 1
                                         else
                  begin
                    G := W[i];
                    H := SrX2Y2(F,G);
                    W[i] := H;
                    C := G/H;
                    S := -F/H;
                    if  MatU  then
                      for j:=1 to M  do
                        begin
                          Y          :=  U[j]^[l1];
                          Z          :=  U[j]^[i];
                          U[j]^[l1]  :=  Y*C+Z*S;
                          U[j]^[i]   :=  -Y*S+Z*C
                        end;
                    inc(i);
                  end;
              end;
          end;

        { Convergency  Check }
        Z := W[k];
        if l<>k  then
          begin
            if its=Itnlimit  then
              begin
                RetCode := k;    Exit;
              end;
            inc(its);
            X  :=  W[l];
            Y  :=  W[k1];
            G  :=  RV1[k1];
            H  :=  RV1[k];
            F  :=  ((Y-Z)*(Y+Z) + (G-H)*(G+H)) /(2.0*H*Y);
            if  abs(F)<=1.0  then   GG := Sign2(sqrt(F*F+1.0),F)
                             else   GG := F*sqrt(1.0+1.0/F/F);
            F  :=  ((X-Z)*(X+Z) + H*(Y/(F+GG)-H)) / X;

       { Next  QR - Transformation  }
            C  :=  1.0;
            S  :=  1.0;
            for i1:=l to k1  do
              begin
                i := i1+1;
                G := RV1[i];
                Y := W[i];
                H := S*G;
                G := C*G;
                Z := SrX2Y2(F,H);
                RV1[i1] := Z;
                C := F/Z;
                S := H/Z;
                F := X*C+G*S;
                G := -X*S+G*C;
                H := Y*S;
                Y := Y*C;
                if  MatV  then
                  for j:=1 to N  do
                    begin
                      X          :=  V[j]^[i1];
                      Z          :=  V[j]^[i];
                      V[j]^[i1]  :=  X*C+Z*S;
                      V[j]^[i]   :=  -X*S+Z*C;
                    end;

                Z := SrX2Y2(F,H);
                W[i1] := Z;
                if Z<>0.0  then
                  begin
                    C := F/Z;
                    S := H/Z;
                  end;
                F := C*G+S*Y;
                X := -S*G+C*Y;
                if  MatU  then
                  for J:=1 to M  do
                    begin
                      Y          :=  U[j]^[i1];
                      Z          :=  U[j]^[i];
                      U[j]^[i1]  :=  Y*C+Z*S;
                      U[j]^[i]   :=  -Y*S+Z*C;
                    end;

              end;

            RV1[l] := 0.0;
            RV1[k] := F;
            W[k]   := X;
          end

          else if Z<0.0  then
            begin
              W[k] := -Z;
              if  MatV  then
                for j:=1 to N  do
                  V[j]^[k] := -V[j]^[k];
            end;
        until  l=k;
    end;{ Of the kk-loop }

  if NA = N  then
    for i:=1 to M  do
      for j:=1 to N  do
      begin
        G := V[i]^[j];
        V[i]^[j] := U[i]^[j];
        U[i]^[j] := G;
      end;
end;{ SVD }

{ sorting }
procedure OrderSVD(M,N: integer; const U,V: Matrix; var W: Vector;
                      MatU, MatV: boolean );
var
  i,k,j: integer;
  P  : RealType;

begin

  {  External loop of the re-ordering  }
  for i:=1 to N-1  do
    begin
      k := i;
      P := W[i];

      { Internal loop: finding of the index of the largest
       singular value among the remaining }
      for j:=i+1 to N  do
        if  W[j]>P  then
          begin
            k := j;   P := W[j];
          end;

      if  k<>i  then
        begin
          { swapping  of the singular values }
          W[k] := W[i];
          W[i] := P;
          for j:=1 to M  do
            begin
              { swapping of the U's columns (if needed) }
              if  MatU  then
                begin
                  P        := U[j]^[i];
                  U[j]^[i] := U[j]^[k];
                  U[j]^[k] := P;
                end;
              { swapping of the V's columns (if needed)}
              if  MatV  then
                begin
                  P        := V[j]^[i];
                  V[j]^[i] := V[j]^[k];
                  V[j]^[k] := P;
                end;
            end;
        end;
    end;
end;

{ Perturbated Cholesky Decomposition }
procedure CholDecomp(N : integer; var HDiag: Vector; MaxOff,MachEps: RealType;
                        const L: Matrix; var MaxAdd: RealType);
var
  i,j,k       : integer;
  MinL, MinL2, S, MinLjj: RealType;
  MaxOffl     : RealType;

begin
  MaxOffl := MaxOff;
  MinL   := sqrt(sqrt(MachEps))*MaxOffl;
  if MaxOffl=0.0  then
    begin
      for i:=1 to N  do
        if abs(HDiag[i])>MaxOffl  then
                    MaxOffl := abs(HDiag[i]);
      MaxOffl := sqrt(MaxOffl);
    end;

  MinL2  := sqrt(MachEps)*MaxOffl;
  MaxAdd := 0.0;
  for j:=1 to N  do
    begin
      S := 0.0;
      if j>1  then
        for i:=1 to j-1  do
          S := S+sqr(L[j]^[i]);
      L[j]^[j] := HDiag[j]-S;
      MinLjj := 0.0;
      if j<N  then
        for i:=j+1 to N  do
          begin
            S := 0.0;
            if j>1  then
              for k:=1 to j-1  do
                S := S+L[i]^[k]*L[j]^[k];
            L[i]^[j] := L[j]^[i]-S;
            if abs(L[i]^[j])>MinLjj  then
                   MinLjj := abs(L[i]^[j]);
          end;
      if MinLjj/MaxOffl>MinL  then
                   MinLjj := MinLjj/MaxOffl
            else   MinLjj := MinL;
      if L[j]^[j]>sqr(MinLjj)  then
                 L[j]^[j] := sqrt(L[j]^[j])
                               else
        begin
          if MinL2>MinLjj  then  MinLjj := MinL2;
          if sqr(MinLjj)-L[j]^[j]>MaxAdd  then
                         MaxAdd := sqr(MinLjj)-L[j]^[j];
          L[j]^[j] := MinLjj;
        end;
      if j<N  then
        for i:=j+1 to N  do
          L[i]^[j] := L[i]^[j]/L[j]^[j];
    end;
end;

{ Cholesky   L - Solution  of L*Y  =  B (given  B)  }
procedure LSolve(N: integer; const L: Matrix;  var B,Y: Vector );
var
  i,j: integer;
begin
  Y[1] := B[1]/L[1]^[1];
  if  N>1  then
    for i:=2 to N  do
      begin
        Y[i] := B[i];
        for j:=1 to i-1  do
          Y[i] := Y[i]-L[i]^[j]*Y[j];
        Y[i] := Y[i]/L[i]^[i];
      end;
end;

{ Cholesky   LT - Solution  of LT*X  =  Y (given Y )   }
procedure LTSolve(N: integer; const L: Matrix;  var Y,X: Vector );
var
  i,j: integer;
begin
  X[N] := Y[N]/L[N]^[N];
  if  N>1  then
    for i:=N-1 downto 1  do
      begin
        X[i] := Y[i];
        for j:=i+1 to N  do
          X[i] := X[i]-L[j]^[i]*X[j];
        X[i] := X[i]/L[i]^[i];
      end;
end;

{ Solution of the equation    L*LT*S = G by the  Cholesky  method    }
procedure ChSolve(N: integer; const L: Matrix; var G,S: Vector );
var
  i: integer;
begin
  LSolve (N,L,G,S);
  LTSolve(N,L,S,S);
  for i:=1 to N  do
    S[i] := -S[i];
end;

{ @lastmod(10.10.2002) }
procedure CholBandDec(N,p: integer; const H,L: Matrix);
var
  i,j,k,m: integer;
  p1,i1,j1: integer;
  S: RealType;
begin
  p1 := p+1;
  for j:=1 to N  do
    begin
      S  := 0.0;
      j1 := p+1-j;
      if j>1  then
        for i:=Max(j-p,1) to j-1  do
          S := S+sqr(L[i+j1]^[i]);
      L[p1]^[j] := sqrt(H[p1]^[j]-S);
      if j<N  then
        for i:=j+1 to Min(j+p,N)  do
          begin
            S  := 0.0;
            i1 := p+1-i;
            m  := Max(i-p,1);
            if j>m  then
              for k:=m to j-1  do
                S := S+L[k+i1]^[k]*L[k+j1]^[k];
            L[j+i1]^[j] := (H[j+i1]^[j]-S)/L[p1]^[j];
          end;
    end;
end;

{ Cholesky   L - Solution  of L*Y  =  B (given  B)  }
procedure LBandSolve(N,p: integer; const L: Matrix; var B,Y: Vector );
var
  i,j,p1,i1: integer;
  S: RealType;
begin
  p1 := p+1;
  Y[1] := B[1]/L[p1]^[1];
  if  N>1  then
    for i:=2 to N  do
      begin
        S  := B[i];
        i1 := p+1-i;
        for j:=Max(i-p,1) to i-1  do
          S := S - L[j+i1]^[j]*Y[j];
        Y[i] := S/L[p1]^[i];
      end;
end;

{ Cholesky   LT - Solution  of LT*X  =  Y (given  Y )   }
procedure LTBandSolve(N,p: integer;  const L: Matrix; var Y,X: Vector );
var
  i,j,p1,i1: integer;
begin
  p1 := p+1;
  X[N] := Y[N]/L[p1]^[N];
  if  N>1  then
    for i:=N-1 downto 1  do
      begin
        X[i] := Y[i];
        i1   := p+i+1;
        for j:=i+1 to Min(i+p,N)  do
          X[i] := X[i]-L[i1-j]^[i]*X[j];
        X[i] := X[i]/L[p1]^[i];
      end;
end;

{  Solution of the equation    L*LT*S = G   by the  Cholesky  method    }
procedure ChBandSolve(N,p: integer; const L: Matrix; var G,S: Vector );
begin
  LBandSolve (N,p,L,G,S);
  LTBandSolve(N,p,L,S,S);
end;

{ QR - Decomposition of the  M  }
procedure QRDecomp(N: integer; const M: Matrix; var M1,M2: Vector;
                    var Sing: boolean);
var
  i,j,k: integer;
  Etha,Norm,Tau: RealType;
begin

  Sing := false;  {  set true if the degeneracy of M is found }
  if N>1  then
    for k:=1 to N-1  do
      begin
        Etha := 0.0;
        for i:=k to N  do
          if abs(M[i]^[k])>Etha  then
            Etha := abs(M[i]^[k]);
        if Etha=0.0  then
          begin
            {  degeneracy is found  }
            M1[k] := 0.0;
            M2[k] := 0.0;
            Sing  := true;
          end        else
          begin
            Norm := 0.0;
            for i:=k to N  do
              begin
                M[i]^[k] := M[i]^[k]/Etha;
                Norm := Norm + sqr(M[i]^[k]);
              end;
            Norm := sqrt(Norm);
            if M[k]^[k]<0.0  then  Norm := -Norm;
            M[k]^[k] := M[k]^[k] + Norm;
            M1[k]    := Norm*M[k]^[k];
            M2[k]    := -Etha*Norm;
            for j:=k+1 to N  do
              begin
                Tau := 0.0;
                for i:=k to N  do
                  Tau := Tau + M[i]^[k]*M[i]^[j];
                Tau := Tau/M1[k];
                for i:=k to N  do
                  M[i]^[j] := M[i]^[j] - Tau*M[i]^[k];
              end;
          end;
      end;

  if M[N]^[N]=0.0  then  Sing := true;
  M2[N] := M[N]^[N];

end;

{ This routine solves  R*X = B  for  X , where  R  is placed in the  M  by the
   QRDecomp  routine ;   X  is returned in the  B }
procedure RSolve(N: integer; const M: Matrix; var M2,B: Vector );
var
  i,j: integer;
begin
  B[N] := B[N]/M2[N];
  if  N>1  then
    for i:=N-1 downto 1  do
      begin
        for j:=i+1 to N  do
          B[i] := B[i]-M[i]^[j]*B[j];
        B[i] := B[i]/M2[i];
      end;
end;

{ This routine solves the equation  (QR)*X = B,
  where Q and R are placed in the  M  by the QRDecomp  routine;
  the computed  X  is returned   in the  B .}
procedure QRSolve(N: integer; const M: Matrix; var M1,M2,B: Vector );
var
  i,j: integer;
  Tau: RealType;
begin
  {  1.  B := QT*B  }
  if  N>1  then
    for j:=1 to N-1  do
      begin
        {  B := Qj*B  }
        Tau := 0.0;
        for i:=j to N  do
          Tau := Tau + M[i]^[j]*B[i];
        Tau := Tau/M1[j];
        for i:=j to N  do
          B[i] := B[i] - Tau*M[i]^[j];
      end;
  {  2.  B := R**(-1)*B  }
  RSolve(N,M,M2,B);
end;

{ Gram-Schmidt Orthogonalization                      }
{ On input:                                           }
{ VT -  row vectors;  Vi[k] = VT[i,k]                 }
{ Nvec - number of vectors to be normalized           }
{ On output:                                          }
{ VT: <Vi|Vi> = 1, <Vi|Vj> = 0, i <> j                }
{ Nvec is the actual number of orthogonalized vectors }
procedure Gram_Schmidt(var Nvec: IntType; N: IntType; var VT: Matrix);
var
  m, k, i: IntType;
  vk, vm: VcPtr;
  x, nn: RealType;
begin
  m := 1;
  while m <= Nvec do
  begin
    // find the largest vector
    i := m;
    nn := 0.0;
    for k := m to Nvec do
    begin
      x := Norm2V(N, VT[k]^); // <Vk|Vk> = |Vk|**2
      if x > nn then
      begin
        nn := x;
        i := k;
      end;
    end;
    // move it in the beginning
    vm := VT[i];
    VT[i] := VT[m];
    VT[m] := vm;
    if nn > MinReal then
    begin // |Vm| > 0
      x := 1.0/sqrt(nn);    // 1/|Vm|
      VxR(N, vm^,vm^,x);   // m-th vector (row) is normalized
      for k := m+1 to Nvec do
      begin
        vk := VT[k];
        x := VTxV(N, vk^,vm^);          // x := <Vk|Vm>
        CombineV(N,vk^,vk^,vm^,1.0,-x); // Vk := Vk - x*Vm
      end;
      inc(m);
    end else // Vm = 0
    begin
      Nvec := m-1;
    end;
  end;
end;

end.