hej , mam mały problem z przetłumaczeniem jednej z procedur z Fortrana na Delphi jest to program do diagonalizacji macierzy, nie znam Fortrana i gdzieś robię jakiś błąd , oto kod Fortrana :
SUBROUTINE jacobi(a,n,np,d,v,nrot)
INTEGER n,np,nrot,NMAX
REAL a(np,np),d(np),v(np,np)
PARAMETER (NMAX=500)
INTEGER i,ip,iq,j
REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX)
do 12 ip=1,n
do 11 iq=1,n
v(ip,iq)=0.
enddo 11
v(ip,ip)=1.
enddo 12
do 13 ip=1,n
b(ip)=a(ip,ip)
d(ip)=b(ip)
z(ip)=0.
enddo 13
nrot=0
do 24 i=1,50
sm=0.
do 15 ip=1,n-1
do 14 iq=ip+1,n
sm=sm+abs(a(ip,iq))
enddo 14
enddo 15
if(sm.eq.0.)return
if(i.lt.4)then
tresh=0.2*sm/n**2
else
tresh=0.
endif
do 22 ip=1,n-1
do 21 iq=ip+1,n
g=100.*abs(a(ip,iq))
if((i.gt.4).and.(abs(d(ip))+g.eq.abs(d(ip)))
* .and.(abs(d(iq))+g.eq.abs(d(iq))))then
a(ip,iq)=0.
else if(abs(a(ip,iq)).gt.tresh)then
h=d(iq)-d(ip)
if(abs(h)+g.eq.abs(h))then
t=a(ip,iq)/h t = 1/(2θ)
else
theta=0.5*h/a(ip,iq) Equation (11.1.10).
t=1./(abs(theta)+sqrt(1.+theta**2))
if(theta.lt.0.)t=-t
endif
c=1./sqrt(1+t**2)
s=t*c
tau=s/(1.+c)
h=t*a(ip,iq)
z(ip)=z(ip)-h
z(iq)=z(iq)+h
d(ip)=d(ip)-h
d(iq)=d(iq)+h
a(ip,iq)=0.
do 16 j=1,ip-1
g=a(j,ip)
h=a(j,iq)
a(j,ip)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
enddo 16
do 17 j=ip+1,iq-1
g=a(ip,j)
h=a(j,iq)
a(ip,j)=g-s*(h+g*tau)
a(j,iq)=h+s*(g-h*tau)
enddo 17
do 18 j=iq+1,n
g=a(ip,j)
h=a(iq,j)
a(ip,j)=g-s*(h+g*tau)
a(iq,j)=h+s*(g-h*tau)
enddo 18
do 19 j=1,n
g=v(j,ip)
h=v(j,iq)
v(j,ip)=g-s*(h+g*tau)
v(j,iq)=h+s*(g-h*tau)
enddo 19
nrot=nrot+1
endif
enddo 21
enddo 22
do 23 ip=1,n
b(ip)=b(ip)+z(ip)
d(ip)=b(ip)
z(ip)=0.
enddo 23
enddo 24
END
a to jest mój w Delphi :
procedure jakobi1(var A,V:t1;var B,D,Z:t2;var N:integer);
var
i,j,k,m,nrot:integer;
sm,Tresh,g,h,teta,t,c,s,tau:real;
BEGIN
// zerowanie elementow macierzy
for i :=1 to N do
begin
for j :=1 to N do
begin
V[i,j] := 0;
end;
V[i,i] := 1;
end;
// wektory i macierze
for i :=1 to N do
begin
b[i] := A[i,i];
d[i] := b[i];
Z[i] := 0;
end;
//-----
nRot := 0;
// glowna petla
for k :=1 to 50 do
Begin
Sm := 0;
// sumowanie elementow
for i :=1 to N-1 do
for j := i + 1 to N do
begin
Sm := Sm + Abs(A[i,j]);
end;
// koniec gdy :
if Sm = 0 then
break;
// przeskalowanie
if k < 4 then
Tresh := 0.2*Sm/(N*N)
else
Tresh := 0;
//petla 2
for i :=1 to N-1 do
begin
for j :=i +1 to N do
Begin
g := 100*Abs(A[i,j]);
// gdy element pozadiagonalny jest maly wowczas
if (k<4)and( Abs(D[i])+g=Abs(D[i]) )and( Abs(D[j])+g=Abs(D[j])) then
A[i,j] := 0
else
BEGIN
if Abs(A[i,j])>tresh then
h := D[i] - D[j];
if Abs(h) + g = Abs(h) then
t := A[i,j]/h
else
begin
teta := 0.5*h/A[i,j];
t := 1/( Abs(teta) + SQRT(1+teta*teta) );
if teta<0 then
t := -t;
end;
c := 1/SQRT(1+t*t);
s := t*c;
tau := s/(1+c);
h := t*A[i,j];
Z[i] := Z[i] - h;
Z[j] := Z[j] + h;
D[i] := D[i] - h;
D[j] := D[j] + h;
A[i,j] := 0;
for m :=1 to i-1 do
begin
g := A[m,i];
h := A[m,j];
A[m,i] := g - s*(h + g*tau);
A[m,j] := g + s*(h - g*tau);
end;
for m :=i+1 to j-1 do
begin
g := A[i,m];
h := A[m,j];
A[i,m] := g - s*(h + g*tau);
A[m,j] := g + s*(h - g*tau);
end;
for m:=i+1 to N do
begin
g := A[i,m];
h := A[j,m];
A[i,m] := g - s*(h + g*tau);
A[j,m] := g + s*(h - g*tau);
end;
for m:=1 to N do
begin
g := V[m,i];
h := V[m,j];
v[m,i] := g - s*(h + g*tau);
V[m,j] := g + s*(h - g*tau);
end;
Nrot := Nrot + 1;
end;
end;
END;
for i :=1 to N do
begin
B[i] := B[i] + Z[i];
D[i] := B[i];
Z[i] := 0;
end;
End; // po k iteracjach
END; // po pr.
Dzieki :)
</delphi>