home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / sigm / vol292 / scilib.lbr / SOLVIT.PQS / SOLVIT.PAS
Encoding:
Pascal/Delphi Source File  |  1986-12-17  |  6.3 KB  |  264 lines
  1. PROGRAM solvit(output);
  2.  (* Feb  8.0, 81 *)
  3.  (* from Pascal Programs for Scientists and Engineers *)
  4.  (* Alan R. Miller, Sybex *)
  5.  (* n x n inverse hilbert matrix *)
  6.  (* solution is 1 1 1 1 1 *)
  7.  (* double precision version *)
  8.  
  9. (* Pascal program to perform simultaneous solution *)
  10. (* by Gauss-Jordan elimination *)
  11.  
  12. CONST
  13.   maxr = 10;
  14.   maxc = 10;
  15.  
  16. TYPE
  17.   ary   = ARRAY[1..maxr] OF real;
  18.   arys  = ARRAY[1..maxc] OF real;
  19.   ary2s = ARRAY[1..maxr, 1..maxc] OF real;
  20.  
  21. VAR
  22.   y          : arys;
  23.   coef       : arys;
  24.   a, b       : ary2s;
  25.   n, m, i, j : integer;
  26.   error      : boolean;
  27.  
  28. PROCEDURE gaussj
  29.   (VAR b     : ary2s; (* square matrix of coefficients *)
  30.     y        : arys;  (* constant vector *)
  31.     VAR coef : arys;  (* solution vector *)
  32.     ncol     : integer;  (* order of matrix *)
  33.     VAR error: boolean); (* true if matrix singular *)
  34.  
  35. (*  Gauss Jordan matrix inversion and solution *)
  36. (*  Adapted from McCormick  *)
  37. (*  Feb  8, 81 *)
  38. (*   B(N,N) coefficient matrix, becomes inverse *)
  39. (*   Y(N)   original constant vector *)
  40. (*   W(N,M) constant vector(s) become solution vector *)
  41. (*   DETERM is the determinant *)
  42. (*   ERROR = 1 if singular *)
  43. (*   INDEX(N,3) *)
  44. (*   NV is number of constant vectors *)
  45.  
  46. LABEL 
  47.    99,98;
  48.  
  49. VAR
  50.   w    : ARRAY[1..maxc, 1..maxc] OF real;
  51.   index: ARRAY[1..maxc, 1..3] OF integer;
  52.   i, j, k, l, nv, irow, icol, n, l1   : integer;
  53.   determ, pivot, hold, sum, t, ab, big: real;
  54.  
  55. PROCEDURE swap(VAR a, b: real);
  56.  
  57. VAR
  58.   hold: real;
  59.  
  60. BEGIN  (* swap *)
  61.   hold := a;
  62.   a := b;
  63.   b := hold
  64. END  (* procedure swap *);
  65.  
  66.  
  67. BEGIN     (* Gauss-Jordan main program *)
  68. (* put gausj2 here *)
  69. (* PROCEDURE gausj2;
  70.  
  71. LABEL 98;
  72.  
  73. VAR
  74.   i, j, k, l, l1: integer;
  75.  
  76. BEGIN *)  (* procedure gausj2 *)
  77.   (* actual start of gaussj *)
  78.   error := false;
  79.   nv := 1 (* single constant vector *);
  80.   n := ncol;
  81.   FOR i := 1 TO n DO
  82.     BEGIN
  83.       w[i, 1] := y[i] (* copy constant vector *);
  84.       index[i, 3] := 0
  85.     END;
  86.   determ := 1.0;
  87.   FOR i := 1 TO n DO
  88.     BEGIN
  89.       (* search for largest element *)
  90.       big := 0.0;
  91.       FOR j := 1 TO n DO
  92.         BEGIN
  93.           IF index[j, 3] <> 1 THEN
  94.             BEGIN
  95.               FOR k := 1 TO n DO
  96.                 BEGIN
  97.                   IF index[k, 3] > 1 THEN
  98.                     BEGIN
  99.                       writeln(' ERROR: matrix singular');
  100.                       error := true;
  101.                       GOTO 98        (* abort *)
  102.                     END;
  103.                   IF index[k, 3] < 1 THEN
  104.                     IF abs(b[j, k]) > big THEN
  105.                       BEGIN
  106.                         irow := j;
  107.                         icol := k;
  108.                         big := abs(b[j, k])
  109.                       END
  110.                 END (* k loop *)
  111.             END
  112.         END (* j loop *);
  113.       index[icol, 3] := index[icol, 3] + 1;
  114.       index[i, 1] := irow;
  115.       index[i, 2] := icol;
  116.  
  117. (*    gausj3 *) (* further subdivision of gaussj *);
  118. (*
  119. PROCEDURE gausj3;
  120.  
  121. VAR
  122.   l: integer;
  123.  
  124. BEGIN *) (* procedure gausj3 *)
  125.  
  126.   (* interchange rows to put pivot on diagonal *)
  127.   IF irow <> icol THEN
  128.     BEGIN
  129.       determ := - determ;
  130.       FOR l := 1 TO n DO
  131.         swap(b[irow, l], b[icol, l]);
  132.       IF nv > 0 THEN
  133.         FOR l := 1 TO nv DO
  134.           swap(w[irow, l], w[icol, l])
  135.     END (* if irow <> icol *)
  136. (* END *) (* gausj3 *);
  137.  
  138.       (* divide pivot row by pivot column *)
  139.       pivot := b[icol, icol];
  140.       determ := determ * pivot;
  141.       b[icol, icol] := 1.0;
  142.       FOR l := 1 TO n DO
  143.         b[icol, l] := b[icol, l] / pivot;
  144.       IF nv > 0 THEN
  145.         FOR l := 1 TO nv DO
  146.           w[icol, l] := w[icol, l] / pivot;
  147.       (*  reduce nonpivot rows *)
  148.       FOR l1 := 1 TO n DO
  149.         BEGIN
  150.           IF l1 <> icol THEN
  151.             BEGIN
  152.               t := b[l1, icol];
  153.               b[l1, icol] := 0.0;
  154.               FOR l := 1 TO n DO
  155.                 b[l1, l] := b[l1, l] - b[icol, l] * t;
  156.               IF nv > 0 THEN
  157.                 FOR l := 1 TO nv DO
  158.                   w[l1, l] := w[l1, l] - w[icol, l] * t;
  159.             END   (* IF l1 <> icol *)
  160.         END
  161.     END (* i loop *);
  162.   98:
  163. (* END *) (* gausj2 *);
  164.  
  165. (*  gausj2 *)  (* first half of gaussj *);
  166.   IF error THEN GOTO 99;
  167.   (* interchange columns *)
  168.   FOR i := 1 TO n DO
  169.     BEGIN
  170.       l := n - i + 1;
  171.       IF index[l, 1] <> index[l, 2] THEN
  172.         BEGIN
  173.           irow := index[l, 1];
  174.           icol := index[l, 2];
  175.           FOR k := 1 TO n DO
  176.             swap(b[k, irow], b[k, icol])
  177.         END (* if index *)
  178.     END  (* i loop *);
  179.   FOR k := 1 TO n DO
  180.     IF index[k, 3] <> 1 THEN
  181.       BEGIN
  182.         writeln(' ERROR: matrix singular');
  183.         error := true;
  184.         GOTO 99   (* abort *)
  185.       END;
  186.   FOR i := 1 TO n DO
  187.     coef[i] := w[i, 1];
  188.   99:
  189. END (* procedure gaussj *);
  190.  
  191.  
  192. PROCEDURE get_data(VAR a : ary2s;
  193.                    VAR y : arys;
  194.                    VAR n, m : integer);
  195.  
  196. (* setup n-by-n hilbert matrix *)
  197.  
  198. VAR
  199.   i, j : integer;
  200.  
  201. BEGIN
  202.   FOR i := 1 TO n DO
  203.     BEGIN
  204.       a[n,i] := 1.0/(n + i - 1);
  205.       a[i,n] := a[n,i]
  206.     END;
  207.   a[n,n] := 1.0/(2*n -1);
  208.   FOR i := 1 TO n DO
  209.     BEGIN
  210.       y[i] := 0.0;
  211.       FOR j := 1 TO n DO
  212.         y[i] := y[i] + a[i,j]
  213.     END;
  214.   writeln;
  215.   IF n < 7 THEN
  216.     BEGIN
  217.       FOR i:= 1 TO n  DO
  218.         BEGIN
  219.           FOR j:= 1 TO m DO
  220.             write( a[i,j] :7:5, '  ');
  221.           writeln( ' : ', y[i] :7:5)
  222.         END;
  223.       writeln
  224.     END  (* if n<7 *)
  225. END (* procedure get_data *);
  226.  
  227. PROCEDURE write_data;
  228.  
  229. (* print out the answers *)
  230.  
  231. VAR
  232.   i : integer;
  233.  
  234. BEGIN
  235.   FOR i := 1 TO m DO
  236.     write( coef[i] :13:9);
  237.   writeln;
  238. END (* write_data *);
  239.  
  240. (* PROCEDURE gaussj
  241.    (VAR a       : ary2s; 
  242.     y           : arys; 
  243.     VAR coef    : arys; 
  244.     ncol        : integer; 
  245.     VAR error   : boolean);
  246. extern; *)
  247. (* F GAUSSJ.PAS *)
  248.  
  249. BEGIN  (* main program *)
  250.   a[1,1] := 1.0;
  251.   n := 2;
  252.   m := n;
  253.   REPEAT
  254.     get_data (a, y, n, m);
  255.     FOR i := 1 TO n DO
  256.       FOR j := 1 TO n DO
  257.         b[i,j] := a[i,j] (* setup work array *);
  258.     gaussj (b, y, coef, n, error);
  259.     IF not error THEN write_data;
  260.     n := n+1;
  261.     m := n
  262.   UNTIL n > maxr
  263. END.
  264. NLEAST4  PAS !ùLEAST5  PAS #φ