Pisze programik, ktory dokonuje pewnych obliczen numerycznych z dziedziny chemii kwantowej. Jednym z etapow jest diagonalizacja macierzy zerojedynkowej.
Jestem na etapie w ktorym dynamiczna macierz kwadratowa NxN jest utworzona i gotowa do nastepnych obliczen. Znalazlem w sieci kilka procedur diagonalizacyjnych mniej lub bardziej zlozonych (jedna bazuje na wskaznikach i w ogole nie wiem o co w tym chodzi (jeszcze) a druga wyglada na calkiem prosta, napisana jest w pascalu.
Przed diagonalizacja macierz na przekatnej ma zera ( zanika wyznacznik).
Po diagonalizacji na przekatnej otrzymujemy wartosci wlasne (jakies liczby :))
Mam problem z implementacja tej procedury.

Czy moglby ktos sprawdzic czy dobrze zaimplementowalem ta procedure.
Z gory dzieki za wszelkie wskazowki.

Oto kod procedury pascalowskiej :

Uses Crt;

Const
 n=2;
 m=2;

Type
 Matrixnm = Array[1..n,1..m] of Double;

Var
 Txt          : Text;
 i,j,count    : Integer;
 A,Abuf,Q,R   : Matrixnm;



Procedure Qrdec(A:Matrixnm;x,y:Integer);
{Global variables : R,Q : Martices}
Var
 i,j,s : Integer;
 sum   : Double;
Begin
 For i:=1 to x do
  Begin
   {Rii = Ai dot Ai:}
   sum:=0;
   For j:=1 to y do
    sum:=sum+A[j,i]*A[j,i];
   sum:=sqrt(sum);
   R[i,i]:=sum;
   {Qi = Ai/Rii:}
   For j:=1 to y do
    Q[j,i]:=A[j,i]/R[i,i];
   {Inner loop:}
   For j:=i+1 to x do
    Begin
     {Qi dot Aj}
     sum:=0;
     For s:=1 to y do
      sum:=sum+Q[s,i]*A[s,j];
     R[i,j]:=sum;
     For s:=1 to y do
      A[s,j]:=A[s,j]-R[i,j]*Q[s,i];
    End;
  End;
End;

Procedure MatmulMat(n,m,p : Integer; A,B : Matrixnm; Var C : Matrixnm);
{Multiplies two matrices A (nxm) and B (mxp). Output is in C and the matrix
type is an Array[1..n,1..m].}
Var sum : Double; i,j,s,t : Integer;
Begin
 For i:=1 to n do
  For j:=1 to m do
   Begin
    sum:=0;
     For s:=1 to m do
      sum:=sum+A[i,s]*B[s,j];
    C[i,j]:=sum;
   End;
End;

Procedure id(var A:Matrixnm; {Inputs:} B:Matrixnm; n,m : Integer);
{Sets A=B (both nxm matrices).}
Var i,j : Integer;
Begin
 For i:=1 to n do
  For j:=1 to m do
   A[i,j]:=B[i,j];
End;

Begin
 Clrscr;
 Assign(Txt,'c:\tp\program.out');
 Rewrite(Txt);

 {Initialization:}
  Count:=0;
  A[1,1]:=1; A[1,2]:=2;
  A[2,1]:=2; A[2,2]:=5;
  id(Abuf,A,n,m);

 {Itterative loop:}
  Repeat
   QRdec(A,n,m);
   Count:=Count+1;
   MatmulMat(n,m,n,R,Q,A);
   delay(10);
   For i:=1 to 6 do
   Writeln('                       ');
   GotoXY(1,1);
   Writeln(A[1,1],' ',A[1,2]);
   Writeln(A[2,1],' ',A[2,2]);
   Writeln(Q[1,1],' ',Q[1,2]);
   Writeln(Q[2,1],' ',Q[2,2]);
   Writeln(R[1,1],' ',R[1,2]);
   Writeln(R[2,1],' ',R[2,2]);
  Until (keypressed) or (A[2,1]=0);
  Writeln('Stopped!');

 {Output to text file:}
  Writeln(Txt,'A matrix before diagonalization:');
  Writeln(Txt,Abuf[1,1],' ',Abuf[1,2]);
  Writeln(Txt,Abuf[2,1],' ',Abuf[2,2]);
  Writeln(Txt);
  Writeln(Txt,'A matrix after diagonalization:');
  Writeln(Txt,A[1,1],' ',A[1,2]);
  Writeln(Txt,A[2,1],' ',A[2,2]);
  Writeln(Txt);
  Writeln(Txt,'Number of calls to QR procedure: ',Count);
 Close(Txt);
 Readln;
End.

a tu po implementacji tej procedury do mojego programu://kod czesciowy//
Matrix1 : --> jest to moja wyjsciowa macierz kwadratowa ktora ma zostac poddana procesowi diagonalizacji.

..........
type
  MatrixNN = Array of array of extended;
..........
private
   
    Procedure Qrdec(var A:Matrixnn; x,y:Integer);
    Procedure MatmulMat(n,m,p : Integer; A,B : MatrixNN; Var C : MatrixNN);
    Procedure id(var Matrix1:MatrixNN;var B:MatrixNN; n,m : Integer);

..........
var
  Form1              : TForm1;
  Matrix1,A,B,Abuf,Q,R,C : MatrixNN;

implementation
..........
procedure TForm1.OblRunClick(Sender: TObject);
begin

   id(Matrix1,B,NrAtomu,NrAtomu);
   Qrdec(Matrix1,High(Matrix1),High(Matrix1));
  end;

Procedure TForm1.Qrdec(var A:MatrixNN;x,y:Integer);
{Global variables : R,Q : Martices}
Var
 i,j,s : Integer;
 sum   : Double;
Begin
SetLength(A,NrAtomu,NrAtomu);
SetLength(Q,NrAtomu,NrAtomu);
SetLength(R,NrAtomu,NrAtomu);
 For i:=0 to x do
  Begin
   {Rii = Ai dot Ai:}
   sum:=0;
   For j:=0 to y do
    sum:=sum+A[j,i]*A[j,i];
   sum:=sqrt(sum);

   R[i,i]:=sum;
   {Qi = Ai/Rii:}
    
   For j:=0 to y do
    Q[j,i]:=A[j,i]/R[i,i];
   {Inner loop:}
   For j:=i+1 to x do
    Begin
     {Qi dot Aj}
     sum:=0;
     For s:=0 to y do
      sum:=sum+Q[s,i]*A[s,j];
     R[i,j]:=sum;
     For s:=0 to y do
      A[s,j]:=A[s,j]-R[i,j]*Q[s,i];
    End;
  End;
  MatmulMat(NrAtomu,NrAtomu,NrAtomu,A,B, C);
End;

Procedure TForm1.id(var Matrix1:MatrixNN; var B:MatrixNN; n,m : Integer);
{Sets A=B (both nxm matrices).}
Var i,j : Integer;
Begin

SetLength(B,NrAtomu, NrAtomu);
 For i:=0 to n-1 do
  For j:=0 to m-1 do
   begin
    B[i,j]:=Matrix1[i,j];
   end;
End;

Procedure TForm1.MatmulMat(n,m,p : Integer; A,B : MatrixNN; Var C : MatrixNN);
{Multiplies two matrices A (nxm) and B (mxp). Output is in C and the matrix
type is an Array[1..n,1..m].}
Var sum : Double; i,j,s : Integer;
Begin

SetLength(C,NrAtomu, NrAtomu);
 For i:=0 to n-1 do
  For j:=0 to m-1 do
   Begin
    sum:=0;
     For s:=0 to m-1 do
      sum:=sum+A[i,s]*B[s,j];
    C[i,j]:=sum;
   End;
End;