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;