QR алгоритм

[Rocky09]

Мне нужен QR алгоритм разложения матрицы(можно с поиском собственных значений),причем не теорию а готовый исходный код. Везде искал, в интернете этого нет. -_-
Скиньте если не жалко)

[Aritaborian]

Хм... Для начала, скажите, на каком именно языке. А впрочем...
http://rosettacode.org/wiki/QR_decomposition
http://www.alglib.net/matrixops/qr.php
http://www.iaa.ncku.edu.tw/~dychiang/la ... ition.html
http://glowingpython.blogspot.com/2011/ ... numpy.html
Ну, и на финал...
http://reference.wolfram.com/mathematic ... ition.html
Или я тут написал чушь, или вы не умеете искать информацию в Сети.

[Rocky09]

на Pascal
мне уже помогли немного,но какието ошибки в обратной матрице,и q матрице в самом алгоритме,вот

type
  TMatrix=array[1..4,1..4] of real;
  
  
  var
  i, j: integer;
   a: TMatrix;                            {матрица коэффициентов}
    n,d: integer;                            {размерность системы}
   data_file: text;                       {файл данных}
    b:TMatrix;
    c:tMatrix;
   k:real;
   q,r:tmatrix;


function MultMatrix(var a,b:tmatrix):tmatrix;//функция умножения матриц
var c:tmatrix;i,j:integer;
begin
for i:=1 to 4 do
for j:=1 to 4 do
c[i,j]:=a[i,j]*b[i,j];
multmatrix:=c;
end;

procedure GetMatr(a:tmatrix; var b:tmatrix; n,i,j:integer); 
{ Вычеркивание из матрицы строки и столбца } 
var ki,kj,di,dj:integer; 
 begin 
    di:=0; 
    for ki:=1 to (n-1) do 
    begin 
       if (ki=i) then di:=1; 
       dj:=0; 
       for kj:=1 to (n-1) do 
       begin 
          if (kj=j) then dj:=1; 
          b[ki,kj]:=a[ki+di,kj+dj]; 
       end; 
    end; 
 end;

 Function Determinant(a:TMatrix;n:integer):real;
 { Вычисление определителя матрицы }
 var i,j,k:longint;d:real;
     b:TMatrix;
   begin
   d:=0; k:=1;
   if (n<1) then
     begin
     writeln('Determinant: Cann''t run. N=',n); halt;
     end;
   if (n=1)
     then d:=a[1,1]
   else if (n=2)
     then d:=a[1,1]*a[2,2]-a[2,1]*a[1,2]
   else { n>2 }
     for i:=1 to n do
       begin
       GetMatr(a,b,n,i,1);
       
       d:=d+k*a[i,1]*Determinant(b,n-1);
       k:=-k;
       end;
   Determinant:=d;
   end;

function inverse(var a:tmatrix):tmatrix;//функция обратной матрицы
   var c:tmatrix;i,j:integer;
   begin
   for i:=1 to 4 do
   for j:=1 to 4 do //ищем матрицу алгебраических дополнений
   begin
   getmatr(a,b,n,i,j); 
   if ((i+j mod 2)<>0) then c[i,j]:=-1*(determinant(b,n-1))
   else c[i,j]:=determinant(b,n-1);
   end;
   
   for i:=1 to 4 do
   for j:=1 to 4 do //транспонируем матрицу алгебраичесских дополнений
   begin
   k:=c[i,j];
   c[i,j]:=c[j,i];
   c[j,i]:=k;
   end;
   
   k:=-1;
   for i:=1 to 4 do 
   begin
   for j:=1 to 4 do 
   begin
   c[j,i]:=k*c[j,i]/d;
   k:=k*-1;
   end;
   end;
   inverse:=c;
  end;
   
procedure QR(const A:TMatrix;var q,r:tmatrix); //функцию qr переделал в процедуру,на выход эти две q и r
var
  i,j,k:integer;
  E,T:TMatrix;
begin
  for i:=1 to N do E[i,i]:=1;
  R:=A;
  for i:=1 to N do
  begin
    for j:=1 to N do Q[i,j]:=1;
  end;
  for k:=1 to N-1 do
  begin
    for i:=k+1 to N do
    begin
      T:=E;
      if R[i,k]<>0 then
        begin
          T[k,k]:=R[k,k]/(sqr(R[k,k])+sqr(A[i,k]));
          T[i,i]:=T[k,k];
          T[k,i]:=R[i,k]/(sqr(R[k,k])+sqr(A[i,k]));
          T[i,k]:=T[k,i];
        end;
      R:=MultMatrix(T,R);
      Q:=MultMatrix(T,Q);
    end;
  end;
  
  
  q:=Inverse(Q);
  
  
end;
   

   begin
   assign(data_file, 'a.txt');
   reset(data_file);
   
   read(data_file, n);
   for i:=1 to n do begin
      for j:=1 to n do begin
         read(data_file, a[i,j]);
      end;  
   end;

   close(data_file);
   
   writeln('matrica'); 
   for i:=1 to n do 
     begin
     for j:=1 to n do   
     write(a[i,j],' ');     
     writeln;
     end;
     
  qr(a,q,r);
  writeln('matrica q'); 
   for i:=1 to 4 do 
     begin
     for j:=1 to 4 do   
     write(q[i,j]:5:3,' ');     
     writeln;
     end;
     
     writeln('matrica r'); 
   for i:=1 to 4 do 
     begin
     for j:=1 to 4 do   
     write(r[i,j]:5:3,' ');     
     writeln;
     end;
      
end.