Delphi/fortran - tłumaczenie kodu

0

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>
0

numerical recipies tez dostepne jest w pascalu

http://chaos.fiz.uni-lj.si/~znidaricm/modelska/pascal/
plik jacobi.pas

1 użytkowników online, w tym zalogowanych: 0, gości: 1