{-----------------------------------------------------------------------------
 Unit Name: ComplexLinSysSolve
 Author:    kaetzel
 Purpose:   Realizes the solution of Linear equation systems with complex
            values in both the coefficient matrix and the result vector. There
            are three algorithms included:
            Direct: LU-Decomposition as a function which gives directly the
                    Vector of the Unknowns
            Iterative: BICGSTAB(l) Bi-Conjugate-Gradients-Stabilized Method
                       which can be applied to non-symmetric, but positive
                       definite matrices
                       GMRES(m) Restarted version of the General Minimum
                       Residual algorithm which does not occupy that much
                       memory than the original GMRES, can be applied to all
                       regular matrices
-----------------------------------------------------------------------------}

unit ComplexLinSysSolve;

interface

uses JclMath,SysUtils,Math,ComplexNumbers,JvDSADialogs;

type
        //Error type specific for the unit
  E_ComplexLinearEquationSystemSolver_Error=Class(Exception);
        //Vector with complex elements and undefined length
  TComplexVector=array of TComplexCart;
        //Matrix with complex elements and undefined size
  TComplexMatrix=array of TComplexVector;

        {                      Public Functions                        }

function LU_Decomposition(CoMatrix:TComplexMatrix;ResultVector:TComplexVector;
                          ResDimension:integer):TComplexVector;
procedure BICGSTAB (CoeffMatrix:TComplexMatrix;ResultVector:TComplexVector;
                    Eps,EpsTilde:real;Nmax:integer;
                    var UnknownsVector:TComplexVector;
                    out NumberOfIterations:integer; out Restart:integer;
                    out FinalResidual:real);
procedure GMRES_m (CoeffMatrix:TComplexMatrix;ResultVector:TComplexVector;
                   Eps:real;MaxKrylov,MaxRestarts:integer;
                   var UnknownsVector:TComplexVector;
                   out Restart:integer;out InnerLoops:integer;
                   out FinalResidual:real);

implementation


{-----------------------------------------------------------------------------
  Procedure: VectorNormEuclidian
  Purpose:   Calculates the Euclidian Norm of a vector with complex elements
  Comments:  Definition -->  ||v||2=Sqrt(Sum(|vi|²))
-----------------------------------------------------------------------------}

function VectorNormEuclidian(vector:TComplexVector):extended;
var
  i:integer;
  sum:extended;
begin
  sum:=0;
  for i:=0 to High(vector) do
    sum:=sum+Power(AbsoluteValue(vector[i]),2);
  Result:=Power(sum,0.5);
end;


{-----------------------------------------------------------------------------
  Procedure: VectorAdd
  Purpose:   Calculates the addition of two vectors of the same length
  Comments:  if they are not of the same length, an error will occur
-----------------------------------------------------------------------------}

function VectorAdd(vector1,vector2:TComplexVector):TComplexVector;
var
  i:integer;
begin
  if length(Vector1)<>length(vector2) then
    Raise E_ComplexLinearEquationSystemSolver_Error.Create('Wrong defintion'+
          ' in function VectorAdd! Dimensions must agree!')
  else
    begin
      SetLength(Result,length(vector1));
      for i:=0 to High(vector1) do
        Result[i]:=ComplexAdd(vector1[i],vector2[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: VectorSub
  Purpose:   Calculates the subtraction of two vectors of the same length
  Comments:  if they are not of the same length, an error will occur
-----------------------------------------------------------------------------}

function VectorSub(vector1,vector2:TComplexVector):TComplexVector;
var
  i:integer;
begin
  if length(Vector1)<>length(vector2) then
    Raise E_ComplexLinearEquationSystemSolver_Error.Create('Wrong defintion'+
          ' in function VectorSub! Dimensions must agree!')
  else
    begin
      SetLength(Result,length(vector1));
      for i:=0 to High(vector1) do
        Result[i]:=ComplexSub(vector1[i],vector2[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: VectorScalarProduct
  Purpose:   Calculates the scalar product of two vectors with complex elements
  Comments:  Definition: e*v=Sum(êi*vi) with ê is the conjugate complex value
             ATTENTION! WITH COMPLEX VALUES e*v IS NOT EQUAL v*e AS WITH REAL
             VALUES!
-----------------------------------------------------------------------------}

function VectorScalarProduct(vector1,vector2:TComplexVector):TComplexCart;
var
  i:integer;
  sum:TComplexCart;
begin
  if length(Vector1)<>length(vector2) then
    Raise E_ComplexLinearEquationSystemSolver_Error.Create('Wrong defintion'+
          ' in function VectorScalarProduct! Dimensions must agree!')
  else
    begin
      sum:=ComplexZeroCart;
      for i:=0 to High(vector1) do
        sum:=ComplexAdd(sum,ComplexMul(Conjugate(vector1[i]),vector2[i]));
      Result:=sum;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: VectorProduct
  Purpose:   Calculates the product of a vector with one single complex scalar
  Comments:
-----------------------------------------------------------------------------}

function VectorProduct(vector1:TComplexVector;
                       factor:TComplexCart):TComplexVector;
var
  i:integer;
begin
  SetLength(Result,length(vector1));
  for i:=0 to High(vector1) do
    Result[i]:=ComplexMul(factor,vector1[i]);
end;


{-----------------------------------------------------------------------------
  Procedure: VectorRealProduct
  Purpose:   Calculates the product of a vector with one single real scalar
  Comments:
-----------------------------------------------------------------------------}

function VectorRealProduct(vector1:TComplexVector;
                           factor:extended):TComplexVector;
var
  i:integer;
begin
  SetLength(Result,length(vector1));
  for i:=0 to High(vector1) do
    Result[i]:=ComplexRealMul(vector1[i],factor);
end;


{-----------------------------------------------------------------------------
  Procedure: MatrixVectorProduct
  Purpose:   Computes the Product of a Matrix with a vector with the same
             length A*b.
  Comments:  The result is a vector with complex elements.
-----------------------------------------------------------------------------}

function MatrixVectorProduct(matrix:TComplexMatrix;
                             vector:TComplexVector):TComplexVector;
var
  i,j:integer;
  sum:TComplexCart;
begin
  if length(vector)<>length(Matrix[0]) then
    Raise E_ComplexLinearEquationSystemSolver_Error.Create('Wrong defintion'+
          ' in function MatrixVectorProduct! Dimensions must agree!')
  else
    begin
      if length(Matrix)<>length(Matrix[0]) then
        Raise E_ComplexLinearEquationSystemSolver_Error.Create('Wrong '+
              'definition in function MatrixVectorProduct! Matrix is not '+
              'quadratic!')
      else
        begin
          SetLength(Result,length(vector));
          for i:=0 to High(vector) do
            begin
              sum:=ComplexZeroCart;
              for j:=0 to High(Matrix) do
                sum:=ComplexAdd(sum,ComplexMul(vector[j],Matrix[i,j]));
              Result[i]:=sum;
            end;
        end;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: ScalingMatrix
  Purpose:   Calculates the scaled matrix with a certain scaling vector to
             improve condition number of the matrix
  Comments:
-----------------------------------------------------------------------------}

procedure ScalingMatrix (ScaleVector:array of extended;
                         var matrix:TComplexMatrix);
var
  i,j:integer;
begin
  for i:=0 to High(ScaleVector) do
    begin
      for j:=0 to High(matrix) do
        matrix[i,j]:=ComplexRealMul(matrix[i,j],ScaleVector[i]);
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: ScalingVector
  Purpose:   Calculates the scaled vector with a certain scaling vector
  Comments:
-----------------------------------------------------------------------------}

procedure ScalingVector (ScaleVector:array of extended;
                         var vector:TComplexVector);
var
  i:integer;
begin
  for i:=0 to High(ScaleVector) do
    vector[i]:=ComplexRealMul(vector[i],ScaleVector[i]);
end;


{-----------------------------------------------------------------------------
  Procedure: ZeroFill
  Purpose:   Overloaded version to fill either a vector or a matrix with
             complex zeros
  Comments:
-----------------------------------------------------------------------------}

procedure ZeroFill (var vector:TComplexVector); overload;
var
  i:integer;
begin
  for i:=0 to high(vector) do
    vector[i]:=ComplexZeroCart;
end;

procedure ZeroFill (var matrix:TComplexMatrix); overload;
var
  i,j:integer;
begin
  for i:=0 to high(matrix) do
    begin
      for j:=0 to high(matrix[0]) do
        matrix[i,j]:=ComplexZeroCart;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: LU_Decomposition
  Purpose:   Solves the linear equation system A*x=b by direct LU-Decomposition
             and forward and backward substitution
  Comments:  Refers to: "Numerical Recipes in C", 1993, p.43-50
-----------------------------------------------------------------------------}

function LU_Decomposition(CoMatrix:TComplexMatrix;ResultVector:TComplexVector;
                          ResDimension:integer):TComplexVector;
const
   tiny=1.0e-20;
var
   k,j,imax,i,ii,ip,n: integer;
   dum,big: extended;
   indx:array of integer;
   vv: array of extended;
   sum,change,dum2:TComplexCart;
   Vector:TComplexVector;
   Matrix:TComplexMatrix;
begin
        //Initialization

  n:=ResDimension;
  SetLength(vv,n);
  SetLength(indx,n);
  SetLength(Matrix,n,n);
  SetLength(Vector,n);

        //Creation of the Sub-Matrix and the Sub-Vector
        //This is necessary as dynamic arrays are not a copy as parameter of
        //the function, so the changed values will be given back
        //If one needs to further work with the arrays, they should not be
        //changed

  Vector:=Copy(ResultVector,0,n);
  for i:=0 to n-1 do
    Matrix[i]:=Copy(CoMatrix[i],0,n);

        //Scaling with maximal element of each row, scaling factors are stored
        //in vector vv

  for i:=1 to n do
    begin
      big:=0;
      for j:=1 to n do
        begin
          if (AbsoluteValue(Matrix[i-1,j-1]) > big) then
            big:=AbsoluteValue(Matrix[i-1,j-1]);
        end;
      if (big=0) then
        Raise E_ComplexLinearEquationSystemSolver_Error.Create('Error in LU-'+
              'Decomposition! Matrix is singular!');
      vv[i-1] := 1.0/big;
    end;

        //Crout's Algorithm to decompose the Matrix A in two Matrices L and U
        //where L is lower triangular and U is upper triangular, the elements
        //of both are stored in the matrix A, so no additional memory is
        //required

  for j:=1 to n do
    begin
      for i:=1 to j-1 do
        begin
          sum:=Matrix[i-1,j-1];
          for k:=1 to i-1 do
            begin
              sum:=ComplexSub(sum,ComplexMul(Matrix[i-1,k-1],Matrix[k-1,j-1]));
            end;
          Matrix[i-1,j-1]:=sum;
        end;
                //search for largest pivot element in the column
      big:=0.0;
      for i:=j to n do
        begin
          sum:=Matrix[i-1,j-1];
          for k:=1 to j-1 do
            begin
              sum:=ComplexSub(sum,ComplexMul(Matrix[i-1,k-1],Matrix[k-1,j-1]));
            end;
          Matrix[i-1,j-1]:=sum;
          dum:=vv[i-1]*AbsoluteValue(sum);
          if (dum>big) then
            begin
              big:=dum;
              imax:=i;
            end;
        end;
                //Change of rows if necessary
      if (j<>imax) then
        begin
          for k:=1 to n do
            begin
              change:=Matrix[imax-1,k-1];
              Matrix[imax-1,k-1]:=Matrix[j-1,k-1];
              Matrix[j-1,k-1]:=change;
            end;
                //Change of scaling factor
          vv[imax-1]:=vv[j-1];
        end;
      indx[j-1]:=imax;

                //If Pivot element is equal to zero it can be replaced by 1e-20
                //but be aware, this is a bad attempt to handle singular
                //matrices
      if ((Matrix[j-1,j-1].Re=0)and(Matrix[j-1,j-1].Im=0)) then
          Matrix[j-1,j-1].Re:=tiny;
                //Division by Pivot element
      if (j<>n) then
        begin
          dum2:=InverseCart(Matrix[j-1,j-1]);
          for i:=j+1 to n do
            begin
              Matrix[i-1,j-1]:=ComplexMul(Matrix[i-1,j-1],dum2);
            end;
        end;
    end;
  SetLength(vv,0);
                //end of LU-Decomposition

        //Starting forward and backward substitution

  ii:=0;
                //Forward substitution
  for i:=1 to n do
    begin
      ip:=indx[i-1];
      sum:=Vector[ip-1];
      Vector[ip-1]:=Vector[i-1];
      if (ii<>0) then
        begin
          for j:=ii to i-1 do
            begin
              sum:=ComplexSub(sum,ComplexMul(Matrix[i-1,j-1],
                              Vector[j-1]));
            end;
        end
      else
        begin
          if ((sum.Re<>0.0)or(sum.Im<>0)) then
            begin
              ii:=i;
            end;
        end;
      Vector[i-1]:=sum;
    end;
                //Backward substitution
  for i:=n downto 1 do
    begin
      sum:=Vector[i-1];
      if (i<n) then
        begin
          for j:=i+1 to n do
            begin
              sum:=ComplexSub(sum,ComplexMul(Matrix[i-1,j-1],
                              Vector[j-1]));
            end;
        end;
                //Results are first stored in the Result Vector
      Vector[i-1]:=ComplexDiv(sum,Matrix[i-1,i-1]);
    end;
  SetLength(Result,ResDimension);
                //Now the results are passed to the unknowns vector
  Result:=Vector;
end;


{-----------------------------------------------------------------------------
  Procedure: BICGSTAB
  Purpose:   Realizes the iterative method Bi-Conjugate-Gradient-Stabilized to
             solve linear equation system A*x=b
  Comments:  Refers to: Meister, A. "Numerik linearer Gleichungssysteme", (1999)
             and Sleijpen, Fokkema El.Trans.Num.Anal. 1 (1993), p.11-32
             Contains all the improvements and restart tests which are proposed
             there!
-----------------------------------------------------------------------------}

procedure BICGSTAB (CoeffMatrix:TComplexMatrix;ResultVector:TComplexVector;
                    Eps,EpsTilde:real;Nmax:integer;
                    var UnknownsVector:TComplexVector;
                    out NumberOfIterations:integer; out Restart:integer;
                    out FinalResidual:real);
var
  p,residual,residual0,s,t,v:TComplexVector;
  alpha,beta,sigma,omega,aid:TComplexCart;
  j,i:integer;
  sum,NormOfB,RelTol,RelTol2:extended;
  d:array of extended;
begin

        {             Initializing the aid vectors                    }

        //First Residual
  SetLength(residual0,length(ResultVector));
        //Search Direction
  SetLength(p,length(ResultVector));
        //Actual Residual
  SetLength(residual,length(ResultVector));
        //Preconditioning vector
  SetLength(d,length(ResultVector));
        //Orthogonal basis of Krylov subspace
  SetLength(v,Length(ResultVector));
  SetLength(s,Length(ResultVector));
  SetLength(t,Length(ResultVector));

        {             Preconditioning of the system                   }

        //Scaling of each row by means of the Norm of the sum of the absolute
        //values of each element
  for i:=0 to High(d) do
    begin
      sum:=0;
      for j:=0 to High(CoeffMatrix) do
        sum:=sum+AbsoluteValue(CoeffMatrix[i,j]);
      d[i]:=1/sum;
    end;
  ScalingMatrix(d,CoeffMatrix);
  ScalingVector(d,ResultVector);
  SetLength(d,0);

        {        Starting values for the vector of the Unkowns        }

        //First try is zero for all values
  for j:=0 to High(UnknownsVector) do
    UnknownsVector[j]:=ComplexZeroCart;


  NormOfB:=VectorNormEuclidian(ResultVector);
  if NormOfB=0 then
   begin
     for i:=0 to high(UnknownsVector) do
       UnknownsVector[i]:=ComplexZeroCart;
     Restart:=0;
     FinalResidual:=0;
     ShowMessage(Format('WARNING!'+
        '%s Result Vector is equal to zero! %s '+
        'Solution will be set to zero!',
        [#10#13#10#13,#10#13#10#13] ),dckActiveForm);
     exit;
   end;
  RelTol:=Eps*NormOfB;
  RelTol2:=Epstilde*NormOfB;

        //Calculation of the first residual
  residual0:=VectorSub(ResultVector,MatrixVectorProduct(CoeffMatrix,
                       UnknownsVector));
        //First direction of the search
  p:=residual0;
        //Actual residual
  residual:=residual0;
  j:=0;
        //This variable counts the restarts of the solution process, it has been
        //involved to control the efficiency of the algorithm
  Restart:=0;

        {                      Iteration Loop                         }

  while VectorNormEuclidian(residual) > RelTol do
    begin
                //Second stop criterion, since it does not involve the real
                //solution, an error message will occur
      if j>Nmax then
        begin
          Raise E_ComplexLinearEquationSystemSolver_Error.Create('Maximum num'+
                'ber of iterations '+IntToStr(Nmax)+' exceeded in procedure'+
                ' BICGSTAB! Try higher Nmax or change starting values!');
          exit;
        end;
      v:=MatrixVectorProduct(CoeffMatrix,p);
      sigma:=VectorScalarProduct(residual0,v);

        {                      Convergence Test                       }

                //Test whether the starting values lead to convergence
                //of the algorithm otherwise the actual parameters are used as
                //starting values and a restart will be made, this is just
                //because the residual of the initial values
                //occurs in every iteration step
      if AbsoluteValue(sigma) > RelTol2*VectorNormEuclidian(v)*
         VectorNormEuclidian(residual0) then
        begin
          alpha:=ComplexDiv(VectorScalarProduct(residual0,residual),sigma);
          s:=VectorSub(residual,VectorProduct(v,alpha));

                //if this value is small enough nearly the exact solution is
                //obtained the equations get much simpler here, the loop will
                //also be stopped at the next run
          if VectorNormEuclidian(s)> RelTol then
            begin
              t:=MatrixVectorProduct(CoeffMatrix,s);
              omega:=ComplexDiv(VectorScalarProduct(t,s),
                                VectorScalarProduct(t,t));
                        //New values for the unknowns are calculated
              UnknownsVector:=VectorAdd(UnknownsVector,VectorAdd(VectorProduct(
                              p,alpha),VectorProduct(s,omega)));
                        //This aid variable has been introduced to save memory
                        //requirement by not needing two residual vectors
                        //for j and j+1
              aid:=VectorScalarProduct(residual0,residual);
                        //New residual
              residual:=VectorSub(s,VectorProduct(t,omega));
              beta:=ComplexMul(ComplexDiv(alpha,omega),ComplexDiv(
                     VectorScalarProduct(residual0,residual),aid));
                        //New search direction
              p:=VectorAdd(residual,VectorProduct(VectorSub(p,VectorProduct(v,
                           omega)),beta));
              j:=j+1
            end
          else
            begin
                        //This is the more simple way of calculation when s
                        //gets small enough, the exact solution is obtained
              UnknownsVector:=VectorAdd(UnknownsVector,VectorProduct(p,alpha));
              residual:=s;
              j:=j+1;
            end;
        end
      else
        begin
                //This causes a restart with the actual parameters of Unknowns
                //as starting values, if j is still zero, other starting values
                //will be tried
          if j=0 then
            begin
              case Restart of
                0: begin
                     for i:=0 to high(UnknownsVector) do
                       UnknownsVector[i]:=ImaginaryUnitCart;
                     residual0:=VectorSub(ResultVector,MatrixVectorProduct(
                                CoeffMatrix,UnknownsVector));
                     p:=residual0;
                     Restart:=Restart+1;
                     continue;
                   end;
                1: begin
                     for i:=0 to high(UnknownsVector) do
                       UnknownsVector[i]:=RealOne;
                     residual0:=VectorSub(ResultVector,MatrixVectorProduct(
                                CoeffMatrix,UnknownsVector));
                     p:=residual0;
                     Restart:=Restart+1;
                     continue;
                   end;
                2: begin
                     for i:=0 to high(UnknownsVector) do
                       UnknownsVector[i]:=RealNegOne;
                     residual0:=VectorSub(ResultVector,MatrixVectorProduct(
                                CoeffMatrix,UnknownsVector));
                     p:=residual0;
                     Restart:=Restart+1;
                     continue;
                   end;
                3: begin
                     for i:=0 to high(UnknownsVector) do
                       UnknownsVector[i]:=NegImaginaryUnit;
                     residual0:=VectorSub(ResultVector,MatrixVectorProduct(
                                CoeffMatrix,UnknownsVector));
                     p:=residual0;
                     Restart:=Restart+1;
                     continue;
                   end;
                else
                  Raise E_ComplexLinearEquationSystemSolver_Error.Create('Err'+
                        'or while trying to find iterative solution for equat'+
                        'ion system! Should try other starting values!');
                end;
            end
          else
            begin
              Restart:=Restart+1;
              j:=j+1;
              residual0:=VectorSub(ResultVector,MatrixVectorProduct(CoeffMatrix,
                         UnknownsVector));
              p:=residual0;
              continue;
            end;
        end;
                //This condition checks whether the relative change in the
                //residual is below a certain boundary value, so the program
                //will need too much iterations a restart is produced
      if AbsoluteValue(alpha)*VectorNormEuclidian(v)/
         VectorNormEuclidian(residual) < 1e-6 then
        begin
          Restart:=Restart+1;
          residual0:=VectorSub(ResultVector,MatrixVectorProduct(CoeffMatrix,
                    UnknownsVector));
          p:=residual0;
          continue;
        end;
                //Regular control of round-up errors by comparing the
                //actual residual with the real residual b-A*x, if the
                //difference is great a restart is caused
      if (j mod 50) =0 then
        if VectorNormEuclidian(VectorSub(residual,VectorSub(ResultVector,
           MatrixVectorProduct(CoeffMatrix,UnknownsVector)))) > 0.000000001 then
        begin
          Restart:=Restart+1;
          residual0:=VectorSub(ResultVector,MatrixVectorProduct(CoeffMatrix,
                    UnknownsVector));
          p:=residual0;
          continue;
        end;
    end;
  NumberOfIterations:=j;
  FinalResidual:=VectorNormEuclidian(residual)/NormOfB;
end;


{-----------------------------------------------------------------------------
  Procedure: GMRES_m
  Purpose:   Calculates the solution to the system A*x=b with the restarted
             General Minimum Residual
  Comments:  Refers to Meister, A. "Numerik linearer Gleichungssysteme"
             and Matlab Help Version 6.1
-----------------------------------------------------------------------------}

procedure GMRES_m (CoeffMatrix:TComplexMatrix;ResultVector:TComplexVector;
                   Eps:real;MaxKrylov,MaxRestarts:integer;
                   var UnknownsVector:TComplexVector;
                   out Restart:integer;out InnerLoops:integer;
                   out FinalResidual:real);
var
  Stop:                                     Boolean;
  i,j,k,n,count:                            integer;
  NormOfB,NormOfRes,RelTol,sum,MinRes,rho:  extended;
  d:                                        array of extended;
  phibar,sum2,rt,c,s,aid,qv:                TComplexCart;
  Hessenberg,Residual,f,u,aidv,x0,vrf:      TComplexVector;
  Q,R,W,V:                                  TComplexMatrix;

begin
        //Test if maximum amount of Krylov subspaces is needed
 if MaxKrylov>length(CoeffMatrix) then
   MaxKrylov:=length(CoeffMatrix);

        {                 Preconditioning of the system                }

        //Scaling of each row by means of the Norm of the sum of the absolute
        //values of each element
  SetLength(d,length(ResultVector));
  for i:=0 to High(d) do
    begin
      sum:=0;
      for j:=0 to High(CoeffMatrix) do
        sum:=sum+AbsoluteValue(CoeffMatrix[i,j]);
      d[i]:=1/sum;
    end;
  ScalingMatrix(d,CoeffMatrix);
  ScalingVector(d,ResultVector);
  SetLength(d,0);

        //Check if the result vector is completely zero, then solution is set to
        //zero
 NormOfB:=VectorNormEuclidian(ResultVector);
 if NormOfB=0 then
   begin
     for i:=0 to high(UnknownsVector) do
       UnknownsVector[i]:=ComplexZeroCart;
     Restart:=0;
     FinalResidual:=0;
     ShowMessage(Format('WARNING!'+
        '%s Result Vector is equal to zero! %s '+
        'Solution will be set to zero!',
        [#10#13#10#13,#10#13#10#13] ),dckActiveForm);
     exit;
   end;

        {              Initialization of the algorithm                 }

  for i:=0 to high(UnknownsVector) do
    UnknownsVector[i]:=ComplexZeroCart;
  RelTol:=eps*NormOfB;
        //First Residual
  Residual:=ResultVector;
  NormOfRes:=VectorNormEuclidian(Residual);

        //Check if start values are a good enough solution
  if NormOfRes<=RelTol then
    begin
      FinalResidual:=NormOfRes/NormOfB;
      Restart:=0;
      ShowMessage(Format('WARNING!'+
        '%s Initial values (Zero) are already a good enough approximation! %s '+
        'Maybe eps should be changed!',
        [#10#13#10#13,#10#13#10#13] ),dckActiveForm);
      exit;
    end;

  SetLength(u,length(CoeffMatrix));
        //Orthonormal basis vectors of the Krylov Subspaces
  SetLength(V,MaxKrylov+1,length(CoeffMatrix));
        //Vector of the upper Hessenberg Matrix
  SetLength(Hessenberg,MaxKrylov+1);
        //Matrix containing the Given's Rotations
  SetLength(Q,MaxKrylov+1,MaxKrylov+1);
        //Upper triangular Matrix for Resolution
  SetLength(R,MaxKrylov,MaxKrylov);
  SetLength(f,MaxKrylov);
  SetLength(W,length(CoeffMatrix),MaxKrylov);
        //Meantime UnknownsVector which is actualized by every restart
  SetLength(x0,length(UnknownsVector));
  x0:=Copy(UnknownsVector,0,length(UnknownsVector));
  Stop:=False;

        {   Outer Loop until Maximum Number of Restarts are exceeded   }

  for n:=1 to MaxRestarts do
    begin
                //Standard Values for all Matrices
      ZeroFill(V);
      ZeroFill(Hessenberg);
      ZeroFill(Q);
      ZeroFill(R);
      ZeroFill(f);
      ZeroFill(W);
      Hessenberg[0].Re:=NormOfRes;
      V[0]:=VectorRealProduct(Residual,(1/(Hessenberg[0].Re)));
      Q[0,0].Re:=1;
      phibar:=Hessenberg[0];
      Count:=0;

                {Inner Loop for the maximum number of Krylov subspaces}

      for j:=1 to MaxKrylov do
        begin
          Count:=Count+1;
          u:=MatrixVectorProduct(CoeffMatrix,v[j-1]);
          for k:=1 to j do
            begin
              Hessenberg[k-1]:=VectorScalarProduct(V[k-1],u);
              u:=VectorSub(u,VectorProduct(V[k-1],Hessenberg[k-1]));
            end;
          Hessenberg[j].Re:=VectorNormEuclidian(u);
          Hessenberg[j].Im:=0;
          V[j]:=VectorProduct(u,InverseCart(Hessenberg[j]));
          for k:=1 to j do
            begin
              sum2:=ComplexZeroCart;
              for i:=1 to j do
                sum2:=ComplexAdd(sum2,ComplexMul(Q[k-1,i-1],Hessenberg[i-1]));
              R[k-1,j-1]:=sum2;
            end;
          rt:=R[j-1,j-1];

                //Givens Rotation
          if Hessenberg[j].Re<1e-20 then
            begin
              c:=RealOne;
              s:=ComplexZeroCart;
            end
          else
            begin
              if AbsoluteValue(Hessenberg[j])>AbsoluteValue(rt) then
                begin
                  aid:=ComplexDiv(rt,Hessenberg[j]);
                  s.Re:=1/(sqrt(1+sqr(AbsoluteValue(aid))));
                  s.Im:=0;
                  c:=ComplexRealMul(aid,-1*s.Re);
                end
              else
                begin
                  aid:=ComplexDiv(Hessenberg[j],rt);
                  c.Re:=1/(sqrt(1+sqr(AbsoluteValue(aid))));
                  c.Im:=0;
                  s:=ComplexRealMul(aid,-1*c.Re);
                end;
            end;
          R[j-1,j-1]:=ComplexSub(ComplexMul(Conjugate(c),rt),
                                 ComplexMul(Conjugate(s),Hessenberg[j]));

                //Calculate the Rotation Matrix
          for k:=1 to j do
            begin
              qv:=Q[j-1,k-1];
              Q[j-1,k-1]:=ComplexMul(Conjugate(c),qv);
              Q[j,k-1]:=ComplexMul(s,qv);
            end;
          Q[j-1,j].Re:=-s.Re;
          Q[j-1,j].Im:=s.Im;
          Q[j,j]:=c;
          f[j-1]:=ComplexMul(Conjugate(c),phibar);
          phibar:=ComplexMul(s,phibar);

          if j<MaxKrylov then
            begin
              if AbsoluteValue(f[j-1])=0 then
                begin
                  Raise E_ComplexLinearEquationSystemSolver_Error.Create('Erro'+
                       'r while evaluating GMRES-Solution, Process stagnates!');
                  exit;
                end;
              for i:=0 to high(W) do
                begin
                  sum2:=ComplexZeroCart;
                  for k:=1 to j-1 do
                    sum2:=ComplexAdd(sum2,ComplexMul(W[i,k-1],R[k-1,j-1]));
                  W[i,j-1]:=ComplexDiv(ComplexSub(v[j-1,i],sum2),R[j-1,j-1]);
                end;
                        //New iterate for inner cycle
              for i:=0 to high(UnknownsVector) do
                UnknownsVector[i]:=ComplexAdd(UnknownsVector[i],
                                              ComplexMul(f[j-1],W[i,j-1]));
            end
                        //Maximum Number of Krylov Subspaces reached
          else
            begin
              SetLength(aidv,j);
              SetLength(vrf,length(x0));
                        //Solution of equation system R*aidv=f
                        //R is upper triangular, therefore only backsubstitution
                        //is needed
              for i:=j downto 1 do
                begin
                  sum2:=f[i-1];
                  if (i<j) then
                    begin
                      for k:=i+1 to j do
                        sum2:=ComplexSub(sum2,ComplexMul(R[i-1,k-1],
                             aidv[k-1]));
                    end;
                  aidv[i-1]:=ComplexDiv(sum2,R[i-1,i-1]);
                 end;
                        //Calculate the new iterate for the outer loop
              for i:=1 to length(V[0]) do
                begin
                  sum2:=ComplexZeroCart;
                  for k:=1 to j do
                    sum2:=ComplexAdd(sum2,ComplexMul(V[k-1,i-1],aidv[k-1]));
                  vrf[i-1]:=sum2;
                end;
              SetLength(aidv,0);
              UnknownsVector:=VectorAdd(x0,vrf);
              SetLength(vrf,0);
            end;

                //New Residual
          Residual:=VectorSub(ResultVector,
                              MatrixVectorProduct(CoeffMatrix,UnknownsVector));
          NormOfRes:=VectorNormEuclidian(Residual);

          if NormOfRes<RelTol then
            begin
              if Count<MaxKrylov then
                begin
                        //Accurate calculation of the UnknownsVector
                  SetLength(aidv,j);
                  SetLength(vrf,length(x0));

                        //Solution of equation system R*aidv=f
                        //R is upper triangular, therefore only backsubstitution
                        //is needed
                  for i:=j downto 1 do
                    begin
                      sum2:=f[i-1];
                        if (i<j) then
                          begin
                            for k:=i+1 to j do
                              sum2:=ComplexSub(sum2,ComplexMul(R[i-1,k-1],
                                               aidv[k-1]));
                          end;
                        aidv[i-1]:=ComplexDiv(sum2,R[i-1,i-1]);
                    end;
                        //Calculate the result
                  for i:=1 to length(V[0]) do
                    begin
                      sum2:=ComplexZeroCart;
                      for k:=1 to j do
                        sum2:=ComplexAdd(sum2,ComplexMul(V[k-1,i-1],aidv[k-1]));
                      vrf[i-1]:=sum2;
                    end;
                  SetLength(aidv,0);
                  UnknownsVector:=VectorAdd(x0,vrf);
                  SetLength(vrf,0);
                  NormOfRes:=VectorNormEuclidian(VectorSub(ResultVector,
                             MatrixVectorProduct(CoeffMatrix,UnknownsVector)));
                end;
              if NormOfRes<RelTol then
                begin
                        //Stop of the inner loop
                  Restart:=n-1;
                  Stop:=True;
                  break;
                end;
            end;
        end;
                //Stop of the outer loop
      if Stop then break;
      x0:=Copy(UnknownsVector,0,length(UnknownsVector));
      Residual:=VectorSub(ResultVector,MatrixVectorProduct(CoeffMatrix,
                                                           UnknownsVector));
    end;
        //Calculate result values
  FinalResidual:=NormOfRes/NormOfB;
  InnerLoops:=Count;
end;

end.
