home *** CD-ROM | disk | FTP | other *** search
/ ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng) / ProfitPress-MegaCDROM2.B6I / PROG / PASCAL / NRPAS13.ZIP / JACOBI.DEM < prev    next >
Encoding:
Text File  |  1991-04-29  |  4.4 KB  |  130 lines
  1. PROGRAM d11r1(input,output);
  2. (* driver for routine JACOBI *)
  3. CONST
  4.    np=10;
  5.    nmat=3;
  6. TYPE
  7.    glnpnp = ARRAY [1..np,1..np] OF real;
  8.    glnp = ARRAY [1..np] OF real;
  9. VAR
  10.    i,j,k,kk,l,ll,nrot : integer;
  11.    a,b,c,v : glnpnp;
  12.    d,r : glnp;
  13.    num : ARRAY [1..3] OF integer;
  14.  
  15. (*$I MODFILE.PAS *)
  16. (*$I JACOBI.PAS *)
  17.  
  18. BEGIN
  19.    num[1] := 3; num[2] := 5; num[3] := 10;
  20.    a[1,1] := 1.0; a[1,2] := 2.0; a[1,3] := 3.0;
  21.    a[2,1] := 2.0; a[2,2] := 2.0; a[2,3] := 3.0;
  22.    a[3,1] := 3.0; a[3,2] := 3.0; a[3,3] := 3.0;
  23.    b[1,1] := -2.0; b[1,2] := -1.0; b[1,3] := 0.0;
  24.    b[1,4] := 1.0; b[1,5] := 2.0;
  25.    b[2,1] := -1.0; b[2,2] := -1.0; b[2,3] := 0.0;
  26.    b[2,4] := 1.0; b[2,5] := 2.0;
  27.    b[3,1] := 0.0; b[3,2] := 0.0; b[3,3] := 0.0;
  28.    b[3,4] := 1.0; b[3,5] := 2.0;
  29.    b[4,1] := 1.0; b[4,2] := 1.0; b[4,3] := 1.0;
  30.    b[4,4] := 1.0; b[4,5] := 2.0;
  31.    b[5,1] := 2.0; b[5,2] := 2.0; b[5,3] := 2.0;
  32.    b[5,4] := 2.0; b[5,5] := 2.0;
  33.    c[1,1] := 5.0; c[1,2] := 4.0; c[1,3] := 3.0;
  34.    c[1,4] := 2.0; c[1,5] := 1.0; c[1,6] := 0.0;
  35.    c[1,7] := -1.0; c[1,8] := -2.0; c[1,9] := -3.0;
  36.    c[1,10] := -4.0;
  37.    c[2,1] := 4.0; c[2,2] := 5.0; c[2,3] := 4.0;
  38.    c[2,4] := 3.0; c[2,5] := 2.0; c[2,6] := 1.0;
  39.    c[2,7] := 0.0; c[2,8] := -1.0; c[2,9] := -2.0;
  40.    c[2,10] := -3.0;
  41.    c[3,1] := 3.0; c[3,2] := 4.0; c[3,3] := 5.0;
  42.    c[3,4] := 4.0; c[3,5] := 3.0; c[3,6] := 2.0;
  43.    c[3,7] := 1.0; c[3,8] := 0.0; c[3,9] := -1.0;
  44.    c[3,10] := -2.0;
  45.    c[4,1] := 2.0; c[4,2] := 3.0; c[4,3] := 4.0;
  46.    c[4,4] := 5.0; c[4,5] := 4.0; c[4,6] := 3.0;
  47.    c[4,7] := 2.0; c[4,8] := 1.0; c[4,9] := 0.0;
  48.    c[4,10] := -1.0;
  49.    c[5,1] := 1.0; c[5,2] := 2.0; c[5,3] := 3.0;
  50.    c[5,4] := 4.0; c[5,5] := 5.0; c[5,6] := 4.0;
  51.    c[5,7] := 3.0; c[5,8] := 2.0; c[5,9] := 1.0;
  52.    c[5,10] := 0.0;
  53.    c[6,1] := 0.0; c[6,2] := 1.0; c[6,3] := 2.0;
  54.    c[6,4] := 3.0; c[6,5] := 4.0; c[6,6] := 5.0;
  55.    c[6,7] := 4.0; c[6,8] := 3.0; c[6,9] := 2.0;
  56.    c[6,10] := 1.0;
  57.    c[7,1] := -1.0; c[7,2] := 0.0; c[7,3] := 1.0;
  58.    c[7,4] := 2.0; c[7,5] := 3.0; c[7,6] := 4.0;
  59.    c[7,7] := 5.0; c[7,8] := 4.0; c[7,9] := 3.0;
  60.    c[7,10] := 2.0;
  61.    c[8,1] := -2.0; c[8,2] := -1.0; c[8,3] := 0.0;
  62.    c[8,4] := 1.0; c[8,5] := 2.0; c[8,6] := 3.0;
  63.    c[8,7] := 4.0; c[8,8] := 5.0; c[8,9] := 4.0;
  64.    c[8,10] := 3.0;
  65.    c[9,1] := -3.0; c[9,2] := -2.0; c[9,3] := -1.0;
  66.    c[9,4] := 0.0; c[9,5] := 1.0; c[9,6] := 2.0;
  67.    c[9,7] := 3.0; c[9,8] := 4.0; c[9,9] := 5.0;
  68.    c[9,10] := 4.0;
  69.    c[10,1] := -4.0; c[10,2] := -3.0; c[10,3] := -2.0;
  70.    c[10,4] := -1.0; c[10,5] := 0.0; c[10,6] := 1.0;
  71.    c[10,7] := 2.0; c[10,8] := 3.0; c[10,9] := 4.0;
  72.    c[10,10] := 5.0;
  73.    FOR i := 1 to nmat DO BEGIN
  74.       IF (i = 1) THEN BEGIN
  75.          jacobi(a,num[i],d,v,nrot)
  76.       END ELSE IF (i = 2) THEN BEGIN
  77.          jacobi(b,num[i],d,v,nrot)
  78.       END ELSE IF (i = 3) THEN BEGIN
  79.          jacobi(c,num[i],d,v,nrot)
  80.       END;
  81.       writeln('matrix number',i:2);
  82.       writeln('number of jacobi rotations:',nrot:3);
  83.       writeln('eigenvalues:');
  84.       FOR j := 1 to num[i] DO BEGIN
  85.          write(d[j]:12:6);
  86.          IF ((j MOD 5) = 0) THEN writeln
  87.       END;
  88.       writeln;
  89.       writeln('eigenvectors:');
  90.       FOR j := 1 to num[i] DO BEGIN
  91.          writeln('number':9,j:3);
  92.          FOR k := 1 to num[i] DO BEGIN
  93.             write(v[k,j]:12:6);         
  94.             IF ((k MOD 5) = 0) THEN writeln
  95.          END;
  96.          writeln
  97.       END;
  98. (* eigenvector test *)
  99.       writeln('eigenvector test');
  100.       FOR j := 1 to num[i] DO BEGIN
  101.          FOR l := 1 to num[i] DO BEGIN
  102.             r[l] := 0.0;
  103.             FOR k := 1 to num[i] DO BEGIN
  104.                IF (k > l) THEN BEGIN
  105.                   kk := l;
  106.                   ll := k
  107.                END ELSE BEGIN
  108.                   kk := k;
  109.                   ll := l
  110.                END;
  111.                IF  (i = 1)  THEN BEGIN
  112.                   r[l] := r[l]+a[ll,kk]*v[k,j]
  113.                END ELSE IF (i = 2) THEN BEGIN
  114.                   r[l] := r[l]+b[ll,kk]*v[k,j]
  115.                END ELSE IF (i = 3) THEN BEGIN
  116.                   r[l] := r[l]+c[ll,kk]*v[k,j]
  117.                END
  118.             END
  119.          END;
  120.          writeln('vector number',j:3);
  121.          writeln('vector':11,'mtrx*vec.':14,'ratio':10);
  122.          FOR l := 1 to num[i] DO BEGIN
  123.             writeln(v[l,j]:12:6,r[l]:12:6,r[l]/v[l,j]:12:6)
  124.          END
  125.       END;
  126.       writeln('press return to continue...');
  127.       readln
  128.    END
  129. END.
  130.