QReigenvalues

program qreigenvalues;
uses crt;
const maxiter = 1000;
type matrix = array of array of Double;
     vector = array of Double;

procedure QRdecomp(var A,Q:matrix;n:integer);
var v:vector;
    sum,A0,sqrsum,temp:double;
    sign:integer;
    i,j,k:integer;
begin
  for i := 0 to n - 1 do
     for j := 0 to n - 1 do
         Q[i,j] := 0;
  for i := 0 to n - 1 do
      Q[i,i] := 1;
  SetLength(v,n);
  for k := 0 to n - 1 do
  begin
    sum := 0;
    A0 := A[k,k];
    if A0 < 0 then
       sign := -1
    else
       sign := 1;
    for i := k to n - 1 do
       sum := sum + sqr(A[i,k]);
    sqrSum := sign*sqrt(sum);
    temp := sqrt(2*(sum + A0 * sqrSum));
    v[k] := (sqrSum + A0)/temp;
    for i := k + 1 to n - 1 do
       v[i] := A[i,k]/temp;
    for j := 0 to n - 1 do
    begin
      sum := 0;
      for i := k to n - 1 do
         sum := sum + v[i]*A[i,j];
      for i := k to n - 1 do
         A[i,j] := A[i,j] - 2*v[i]*sum;
    end;
    for j := 0 to n - 1 do
    begin
      sum := 0;
      for i := k to n - 1 do
         sum := sum + v[i]*Q[i,j];
      for i := k to n - 1 do
         Q[i,j] := Q[i,j] - 2*v[i]*sum;
    end;
  end;
  for i := 0 to n - 1 do
  for j := i + 1 to n - 1 do
  begin
    temp := Q[i,j];
    Q[i,j] := Q[j,i];
    Q[j,i] := temp;
  end;
end;

procedure multiplyMatrix(A,B:matrix;var C:matrix;n:integer);
var i,j,k:integer;
    sum : double;
begin
  for i := 0 to n - 1 do
      for j := 0 to n - 1 do
      begin
        sum := 0;
        for k := 0 to n - 1 do
          sum := sum + A[i,k]*B[k,j];
          C[i,j] := sum;
      end;
end;

procedure copyMatrix(A:matrix;var B:matrix;n:integer);
var i,j:integer;
begin
  for i := 0 to n - 1 do
     for j := 0 to n - 1 do
        B[i,j] := A[i,j];
end;

procedure printMatrix(var f:text;A:matrix;n:integer);
var i,j:integer;
begin
  for i := 0 to n - 1 do
    begin
      for j := 0 to n - 1 do
        write(f,A[i,j]:1:12,' ');
      writeln(f);
    end;
    writeln(f);
end;

var A,Q,R:matrix;
    i,j,n:integer;
    esc:char;
    f:text;

begin
  clrscr;
  assign(f,'eig.txt');
  {$I-}
  reset(f);
  {$I+}
  if IOResult = 0 then
    append(f)
  else
    Rewrite(f);
  repeat
    writeln('Podaj rozmiar macierzy');
    readln(n);
    SetLength(A,n,n);
    SetLength(Q,n,n);
    SetLength(R,n,n);
    for i := 0 to n- 1 do
       for j := 0 to n-1 do
       begin
         write('A[',i+1,',',j+1,']=');
         readln(A[i,j]);
       end;
    printmatrix(f,A,n);
    for i := 1 to maxiter do
    begin
      QRdecomp(A,Q,n);
      copyMatrix(A,R,n);
      multiplyMatrix(R,Q,A,n);
    end;
    printmatrix(f,A,n);
    esc := readkey;
  until esc = #27;
  close(f);
end.