CONST
   np = 20;
TYPE
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
   RealArrayNP = ARRAY [1..np] OF real;
{Pocita vlastni cisla a vlastni vektory symetricke matice A typu N x N
D vraci vlastni cisla ve svych prvnich N elementech a sloupce matice
V jsou na vystupu rovny normalizovanym vlastnim vektorum matice A.
NROT je pocet provedenych Jacobiho rotaci.}
PROCEDURE jacobi(VAR a: RealArrayNPbyNP;
                     n: integer;
                 VAR d: RealArrayNP;
                 VAR v: RealArrayNPbyNP;
              VAR nrot: integer);
LABEL 99;
VAR
   j,iq,ip,i: integer;
   tresh,theta,tau,t,sm,s,h,g,c: real;
   b,z: ^RealArrayNP;
BEGIN
   new(b);
   new(z);
   FOR ip := 1 TO n DO BEGIN            {inicializuje jednotkovou matici}
      FOR iq := 1 TO n DO v[ip,iq] := 0.0;
      v[ip,ip] := 1.0
   END;
   FOR ip := 1 TO n DO BEGIN          {incializuje b a d = diag(a)}
      b^[ip] := a[ip,ip];
      d[ip] := b^[ip];
      z^[ip] := 0.0                   {z obsahuje cleny t*a_pq }
   END;
   nrot := 0;
   FOR i := 1 TO 50 DO BEGIN  {max 50 pruchodu - rotace vsech nediag.elem.}
      Writeln('ipruch= ',i, ' nrot =', nrot);
      sm := 0.0;
      FOR ip := 1 TO n-1 DO           {suma mimodiagonalnich prvku}
         FOR iq := ip+1 TO n DO
            sm := sm+abs(a[ip,iq]);
      IF sm = 0.0 THEN GOTO 99;    {normalni konec pri automatickem
                             polozeni podteceni --> = 0 - kvadrat. konverg.}
      IF i < 4 THEN tresh := 0.2*sm/sqr(n) {prvni 3 pruchody-min.pro rotaci}
      ELSE tresh := 0.0;                    {... potom}
      FOR ip := 1 TO n-1 DO BEGIN
         FOR iq := ip+1 TO n DO BEGIN
            g := 100.0*abs(a[ip,iq]);
            IF (i > 4) AND (abs(d[ip])+g = abs(d[ip]))  {po prvnich 4
                                 pruchodech preskakujeme rotace pro male
                                 mimodiagonalni elementy}
               AND (abs(d[iq])+g = abs(d[iq])) THEN a[ip,iq] := 0.0
            ELSE IF abs(a[ip,iq]) > tresh THEN BEGIN
               h := d[iq]-d[ip];
               IF abs(h)+g = abs(h) THEN
                  t := a[ip,iq]/h         {t=1/2*theta - t pro velka theta}
               ELSE BEGIN
                  theta := 0.5*h/a[ip,iq];     {t pro ostatni theta}
                  t := 1.0/(abs(theta)+sqrt(1.0+sqr(theta)));
                  IF theta < 0.0 THEN t := -t
               END;
               c := 1.0/sqrt(1+sqr(t));
               s := t*c;
               tau := s/(1.0+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.0;
               FOR j := 1 TO ip-1 DO BEGIN    {elementy 1 <= j < p}
                  g := a[j,ip];
                  h := a[j,iq];
                  a[j,ip] := g-s*(h+g*tau);
                  a[j,iq] := h+s*(g-h*tau)
               END;
               FOR j := ip+1 TO iq-1 DO BEGIN   {elementy p < j < q}
                  g := a[ip,j];
                  h := a[j,iq];
                  a[ip,j] := g-s*(h+g*tau);
                  a[j,iq] := h+s*(g-h*tau)
               END;
               FOR j := iq+1 TO n DO BEGIN      {elementy q < j <= 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)
               END;
               FOR j := 1 TO n DO BEGIN   {vynasobeni matice v rotaci}
                  g := v[j,ip];
                  h := v[j,iq];
                  v[j,ip] := g-s*(h+g*tau);
                  v[j,iq] := h+s*(g-h*tau)
               END;
               nrot := nrot+1
            END
         END
      END;
      FOR ip := 1 TO n DO BEGIN
         b^[ip] := b^[ip]+z^[ip];
         d[ip] := b^[ip];              {dava do d sumu t*a_pq}
         z^[ip] := 0.0                 {reinicializuje z}
      END
   END;
   writeln('pause in routine JACOBI');
   writeln('50 iterations should not happen');
   readln;
99:
   dispose(z);
   dispose(b)
END;
VAR a,apom,vmat: RealArrayNPbyNP;
    n,i,j,k, nrot: integer;
    d, vnas: RealArrayNP;
    fvst:Text;
    vstnam:String;
Begin
    Writeln('Zadej jmeno vstupniho souboru');
    Readln(vstnam);
    ASSIGN(fvst,vstnam);
    RESET(fvst);
    Readln(fvst,n);
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Read(fvst,a[i,j]); Readln(fvst)
        END;
    FOR i:=1 TO n DO
        FOR j:=1 TO n DO
                   apom[i,j] := a[i,j];
   Jacobi(a,n,d,vmat,nrot);
   Writeln(' K uprave zapotrebi ',nrot,'  rotaci');
   Readln;
   FOR i:=1 TO n DO BEGIN
        Writeln(I:2,'-te vlastni cislo je ',d[i]);
        Writeln(' Vlastni vektor je ');
        FOR j:=1 TO n DO Write('  ',VMAT[j,i],' ');
        Writeln;
        Writeln(' A * vl.vektor je ');
        FOR j:=1 TO n DO BEGIN
              vnas[j] := 0.;
              FOR K:=1 to n DO
                  vnas[j] := vnas[j] + apom[j,k]*vmat[k,i];
              Write('  ',VNAS[j],'  ');
              END;
        Writeln; Writeln; END;
   Readln;

end.