home *** CD-ROM | disk | FTP | other *** search
/ Turbo Toolbox / Turbo_Toolbox.iso / turbo5 / hilb.pas < prev    next >
Pascal/Delphi Source File  |  1988-10-09  |  6KB  |  246 lines
  1. {$N-}
  2. program Hilb;
  3. {
  4.  
  5.   The program performs simultaneous solution by Gauss-Jordan
  6.   elimination.
  7.  
  8.   --------------------------------------------------
  9.   From: Pascal Programs for Scientists and Engineers
  10.  
  11.   Alan R. Miller, Sybex
  12.   n x n inverse hilbert matrix
  13.   solution is 1 1 1 1 1
  14.   double precision version
  15.   --------------------------------------------------
  16.  
  17.   INSTRUCTIONS
  18.   1.  Compile and run the program using the $N- (Numeric Processing :
  19.       Software) compiler directive.
  20.   2.  if you have a math coprocessor in your computer, compile and run the
  21.       program using the $N+ (Numeric Processing : Hardware) compiler
  22.       directive.  Compare the speed and precision of the results to those
  23.       of example 1.
  24. }
  25.  
  26. const
  27.   maxr = 10;
  28.   maxc = 10;
  29.  
  30. type
  31. {$IFOPT N+}                        { use extended type if using 80x87 }
  32.   real  = extended;
  33. {$ENDIF}
  34.   ary   = array[1..maxr] of real;
  35.   arys  = array[1..maxc] of real;
  36.   ary2s = array[1..maxr, 1..maxc] of real;
  37.  
  38. var
  39.   y          : arys;
  40.   coef       : arys;
  41.   a, b       : ary2s;
  42.   n, m, i, j : integer;
  43.   error      : boolean;
  44.  
  45. procedure gaussj
  46.   (var b     : ary2s;  (* square matrix of coefficients *)
  47.     y        : arys;  (* constant vector *)
  48.     var coef : arys;  (* solution vector *)
  49.     ncol     : integer;  (* order of matrix *)
  50.     var error: boolean); (* true if matrix singular *)
  51.  
  52. (*  Gauss Jordan matrix inversion and solution *)
  53. (*  Adapted from McCormick  *)
  54. (*  Feb  8, 81 *)
  55. (*   B(N,N) coefficient matrix, becomes inverse *)
  56. (*   Y(N)   original constant vector *)
  57. (*   W(N,M) constant vector(s) become solution vector *)
  58. (*   DETERM is the determinant *)
  59. (*   ERROR = 1 if singular *)
  60. (*   INDEX(N,3) *)
  61. (*   NV is number of constant vectors *)
  62.  
  63. var
  64.   w    : array[1..maxc, 1..maxc] of real;
  65.   index: array[1..maxc, 1..3] of integer;
  66.   i, j, k, l, nv, irow, icol, n, l1   : integer;
  67.   determ, pivot, hold, sum, t, ab, big: real;
  68.  
  69. procedure swap(var a, b: real);
  70.  
  71. var
  72.   hold: real;
  73.  
  74. begin  (* swap *)
  75.   hold := a;
  76.   a := b;
  77.   b := hold
  78. end  (* procedure swap *);
  79.  
  80.  
  81. begin     (* Gauss-Jordan main program *)
  82.   error := false;
  83.   nv := 1 (* single constant vector *);
  84.   n := ncol;
  85.   for i := 1 to n do
  86.     begin
  87.       w[i, 1] := y[i] (* copy constant vector *);
  88.       index[i, 3] := 0
  89.     end;
  90.   determ := 1.0;
  91.   for i := 1 to n do
  92.     begin
  93.       (* search for largest element *)
  94.       big := 0.0;
  95.       for j := 1 to n do
  96.         begin
  97.           if index[j, 3] <> 1 then
  98.             begin
  99.               for k := 1 to n do
  100.                 begin
  101.                   if index[k, 3] > 1 then
  102.                     begin
  103.                       writeln(' ERROR: matrix singular');
  104.                       error := true;
  105.                       exit;         (* abort *)
  106.                     end;
  107.                   if index[k, 3] < 1 then
  108.                     if abs(b[j, k]) > big then
  109.                       begin
  110.                         irow := j;
  111.                         icol := k;
  112.                         big := abs(b[j, k])
  113.                       end
  114.                 end (* k loop *)
  115.             end
  116.         end (* j loop *);
  117.       index[icol, 3] := index[icol, 3] + 1;
  118.       index[i, 1] := irow;
  119.       index[i, 2] := icol;
  120.  
  121.   (* interchange rows to put pivot on diagonal *)
  122.   if irow <> icol then
  123.     begin
  124.       determ := - determ;
  125.       for l := 1 to n do
  126.         swap(b[irow, l], b[icol, l]);
  127.       if nv > 0 then
  128.         for l := 1 to nv do
  129.           swap(w[irow, l], w[icol, l])
  130.     end; (* if irow <> icol *)
  131.  
  132.       (* divide pivot row by pivot column *)
  133.       pivot := b[icol, icol];
  134.       determ := determ * pivot;
  135.       b[icol, icol] := 1.0;
  136.       for l := 1 to n do
  137.         b[icol, l] := b[icol, l] / pivot;
  138.       if nv > 0 then
  139.         for l := 1 to nv do
  140.           w[icol, l] := w[icol, l] / pivot;
  141.       (*  reduce nonpivot rows *)
  142.       for l1 := 1 to n do
  143.         begin
  144.           if l1 <> icol then
  145.             begin
  146.               t := b[l1, icol];
  147.               b[l1, icol] := 0.0;
  148.               for l := 1 to n do
  149.                 b[l1, l] := b[l1, l] - b[icol, l] * t;
  150.               if nv > 0 then
  151.                 for l := 1 to nv do
  152.                   w[l1, l] := w[l1, l] - w[icol, l] * t;
  153.             end   (* if l1 <> icol *)
  154.         end
  155.     end (* i loop *);
  156.  
  157.   if error then exit;
  158.   (* interchange columns *)
  159.   for i := 1 to n do
  160.     begin
  161.       l := n - i + 1;
  162.       if index[l, 1] <> index[l, 2] then
  163.         begin
  164.           irow := index[l, 1];
  165.           icol := index[l, 2];
  166.           for k := 1 to n do
  167.             swap(b[k, irow], b[k, icol])
  168.         end (* if index *)
  169.     end  (* i loop *);
  170.   for k := 1 to n do
  171.     if index[k, 3] <> 1 then
  172.       begin
  173.         writeln(' ERROR: matrix singular');
  174.         error := true;
  175.         exit;   (* abort *)
  176.       end;
  177.   for i := 1 to n do
  178.     coef[i] := w[i, 1];
  179. end (* procedure gaussj *);
  180.  
  181.  
  182. procedure get_data(var a : ary2s;
  183.                    var y : arys;
  184.                    var n, m : integer);
  185.  
  186. (* setup n-by-n hilbert matrix *)
  187.  
  188. var
  189.   i, j : integer;
  190.  
  191. begin
  192.   for i := 1 to n do
  193.     begin
  194.       a[n,i] := 1.0/(n + i - 1);
  195.       a[i,n] := a[n,i]
  196.     end;
  197.   a[n,n] := 1.0/(2*n -1);
  198.   for i := 1 to n do
  199.     begin
  200.       y[i] := 0.0;
  201.       for j := 1 to n do
  202.         y[i] := y[i] + a[i,j]
  203.     end;
  204.   writeln;
  205.   if n < 7 then
  206.     begin
  207.       for i:= 1 to n  do
  208.         begin
  209.           for j:= 1 to m do
  210.             write( a[i,j] :7:5, '  ');
  211.           writeln( ' : ', y[i] :7:5)
  212.         end;
  213.       writeln
  214.     end  (* if n<7 *)
  215. end (* procedure get_data *);
  216.  
  217. procedure write_data;
  218.  
  219. (* print out the answers *)
  220.  
  221. var
  222.   i : integer;
  223.  
  224. begin
  225.   for i := 1 to m do
  226.     write( coef[i] :13:9);
  227.   writeln;
  228. end (* write_data *);
  229.  
  230.  
  231. begin  (* main program *)
  232.   a[1,1] := 1.0;
  233.   n := 2;
  234.   m := n;
  235.   repeat
  236.     get_data (a, y, n, m);
  237.     for i := 1 to n do
  238.       for j := 1 to n do
  239.         b[i,j] := a[i,j] (* setup work array *);
  240.     gaussj (b, y, coef, n, error);
  241.     if not error then write_data;
  242.     n := n+1;
  243.     m := n
  244.   until n > maxr;
  245. end.
  246.