home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / JACOBI.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1991-04-29  |  3.1 KB  |  100 lines
  1. PROCEDURE jacobi(VAR a: glnpnp; n: integer;VAR d: glnp;
  2.        VAR v: glnpnp; VAR nrot: integer);
  3. (* Programs using routine JACOBI must define the types
  4. TYPE
  5.    glnpnp = ARRAY [1..np,1..np] OF real;
  6.    glnp = ARRAY [1..np] OF real;
  7. where 'np by np' is the physical dimension of the array
  8. a into which all arrays are loaded for analysis. *)
  9. LABEL 99;
  10. CONST
  11.    nmax=100;
  12. VAR
  13.    j,iq,ip,i: integer;
  14.    tresh,theta,tau,t,sm,s,h,g,c: real;
  15.    b,z: ARRAY [1..nmax] OF real;
  16. BEGIN
  17.    FOR ip := 1 TO n DO BEGIN
  18.       FOR iq := 1 TO n DO BEGIN
  19.          v[ip,iq] := 0.0
  20.       END;
  21.       v[ip,ip] := 1.0
  22.    END;
  23.    FOR ip := 1 TO n DO BEGIN
  24.       b[ip] := a[ip,ip];
  25.       d[ip] := b[ip];
  26.       z[ip] := 0.0
  27.    END;
  28.    nrot := 0;
  29.    FOR i := 1 TO 50 DO BEGIN
  30.       sm := 0.0;
  31.       FOR ip := 1 TO n-1 DO BEGIN
  32.          FOR iq := ip+1 TO n DO BEGIN
  33.             sm := sm+abs(a[ip,iq])
  34.          END
  35.       END;
  36.       IF (sm = 0.0) THEN GOTO 99;
  37.       IF (i < 4) THEN tresh := 0.2*sm/sqr(n)
  38.       ELSE tresh := 0.0;
  39.       FOR ip := 1 TO n-1 DO BEGIN
  40.          FOR iq := ip+1 TO n DO BEGIN
  41.             g := 100.0*abs(a[ip,iq]);
  42.             IF ((i > 4) AND ((abs(d[ip])+g) = abs(d[ip]))
  43.                AND ((abs(d[iq])+g) = abs(d[iq]))) THEN
  44.                a[ip,iq] := 0.0
  45.             ELSE IF (abs(a[ip,iq]) > tresh) THEN BEGIN
  46.                h := d[iq]-d[ip];
  47.                IF ((abs(h)+g) = abs(h)) THEN BEGIN
  48.                   t := a[ip,iq]/h
  49.                END ELSE BEGIN
  50.                   theta := 0.5*h/a[ip,iq];
  51.                   t := 1.0/(abs(theta)+sqrt(1.0+sqr(theta)));
  52.                   IF (theta < 0.0) THEN t := -t
  53.                END;
  54.                c := 1.0/sqrt(1+sqr(t));
  55.                s := t*c;
  56.                tau := s/(1.0+c);
  57.                h := t*a[ip,iq];
  58.                z[ip] := z[ip]-h;
  59.                z[iq] := z[iq]+h;
  60.                d[ip] := d[ip]-h;
  61.                d[iq] := d[iq]+h;
  62.                a[ip,iq] := 0.0;
  63.                FOR j := 1 TO ip-1 DO BEGIN
  64.                   g := a[j,ip];
  65.                   h := a[j,iq];
  66.                   a[j,ip] := g-s*(h+g*tau);
  67.                   a[j,iq] := h+s*(g-h*tau)
  68.                END;
  69.                FOR j := ip+1 TO iq-1 DO BEGIN
  70.                   g := a[ip,j];
  71.                   h := a[j,iq];
  72.                   a[ip,j] := g-s*(h+g*tau);
  73.                   a[j,iq] := h+s*(g-h*tau)
  74.                END;
  75.                FOR j := iq+1 TO n DO BEGIN
  76.                   g := a[ip,j];
  77.                   h := a[iq,j];
  78.                   a[ip,j] := g-s*(h+g*tau);
  79.                   a[iq,j] := h+s*(g-h*tau)
  80.                END;
  81.                FOR j := 1 TO n DO BEGIN
  82.                   g := v[j,ip];
  83.                   h := v[j,iq];
  84.                   v[j,ip] := g-s*(h+g*tau);
  85.                   v[j,iq] := h+s*(g-h*tau)
  86.                END;
  87.                nrot := nrot+1
  88.             END
  89.          END
  90.       END
  91.    END;
  92.    FOR ip := 1 TO n DO BEGIN
  93.       b[ip] := b[ip]+z[ip];
  94.       d[ip] := b[ip];
  95.       z[ip] := 0.0
  96.    END;
  97.    writeln('pause in routine JACOBI');
  98.    writeln('50 iterations should not happen'); readln;
  99. 99:   END;
  100.