⟦807b1c8a1⟧

TextFile

æ$N-å
program Hilb;
æ

  The program performs simultaneous solution by Gauss-Jordan
  elimination.

  --------------------------------------------------
  From: Pascal Programs for Scientists and Engineers

  Alan R. Miller, Sybex
  n x n inverse hilbert matrix
  solution is 1 1 1 1 1
  double precision version
  --------------------------------------------------

  INSTRUCTIONS
  1.  Compile and run the program using the $N- (Numeric Processing :
      Software) compiler directive.
  2.  if you have a math coprocessor in your computer, compile and run the
      program using the $N+ (Numeric Processing : Hardware) compiler
      directive.  Compare the speed and precision of the results to those
      of example 1.
å

const
  maxr = 10;
  maxc = 10;

type
æ$IFOPT N+å                        æ use extended type if using 80x87 å
  real  = extended;
æ$ENDIFå
  ary   = arrayÆ1..maxrÅ of real;
  arys  = arrayÆ1..maxcÅ of real;
  ary2s = arrayÆ1..maxr, 1..maxcÅ of real;

var
  y          : arys;
  coef       : arys;
  a, b       : ary2s;
  n, m, i, j : integer;
  error      : boolean;

procedure gaussj
  (var b     : ary2s;  (* square matrix of coefficients *)
    y        : arys;  (* constant vector *)
    var coef : arys;  (* solution vector *)
    ncol     : integer;  (* order of matrix *)
    var error: boolean); (* true if matrix singular *)

(*  Gauss Jordan matrix inversion and solution *)
(*  Adapted from McCormick  *)
(*  Feb  8, 81 *)
(*   B(N,N) coefficient matrix, becomes inverse *)
(*   Y(N)   original constant vector *)
(*   W(N,M) constant vector(s) become solution vector *)
(*   DETERM is the determinant *)
(*   ERROR = 1 if singular *)
(*   INDEX(N,3) *)
(*   NV is number of constant vectors *)

var
  w    : arrayÆ1..maxc, 1..maxcÅ of real;
  index: arrayÆ1..maxc, 1..3Å of integer;
  i, j, k, l, nv, irow, icol, n, l1   : integer;
  determ, pivot, hold, sum, t, ab, big: real;

procedure swap(var a, b: real);

var
  hold: real;

begin  (* swap *)
  hold := a;
  a := b;
  b := hold
end  (* procedure swap *);


begin     (* Gauss-Jordan main program *)
  error := false;
  nv := 1 (* single constant vector *);
  n := ncol;
  for i := 1 to n do
    begin
      wÆi, 1Å := yÆiÅ (* copy constant vector *);
      indexÆi, 3Å := 0
    end;
  determ := 1.0;
  for i := 1 to n do
    begin
      (* search for largest element *)
      big := 0.0;
      for j := 1 to n do
        begin
          if indexÆj, 3Å <> 1 then
            begin
              for k := 1 to n do
                begin
                  if indexÆk, 3Å > 1 then
                    begin
                      writeln(' ERROR: matrix singular');
                      error := true;
                      exit;         (* abort *)
                    end;
                  if indexÆk, 3Å < 1 then
                    if abs(bÆj, kÅ) > big then
                      begin
                        irow := j;
                        icol := k;
                        big := abs(bÆj, kÅ)
                      end
                end (* k loop *)
            end
        end (* j loop *);
      indexÆicol, 3Å := indexÆicol, 3Å + 1;
      indexÆi, 1Å := irow;
      indexÆi, 2Å := icol;

  (* interchange rows to put pivot on diagonal *)
  if irow <> icol then
    begin
      determ := - determ;
      for l := 1 to n do
        swap(bÆirow, lÅ, bÆicol, lÅ);
      if nv > 0 then
        for l := 1 to nv do
          swap(wÆirow, lÅ, wÆicol, lÅ)
    end; (* if irow <> icol *)

      (* divide pivot row by pivot column *)
      pivot := bÆicol, icolÅ;
      determ := determ * pivot;
      bÆicol, icolÅ := 1.0;
      for l := 1 to n do
        bÆicol, lÅ := bÆicol, lÅ / pivot;
      if nv > 0 then
        for l := 1 to nv do
          wÆicol, lÅ := wÆicol, lÅ / pivot;
      (*  reduce nonpivot rows *)
      for l1 := 1 to n do
        begin
          if l1 <> icol then
            begin
              t := bÆl1, icolÅ;
              bÆl1, icolÅ := 0.0;
              for l := 1 to n do
                bÆl1, lÅ := bÆl1, lÅ - bÆicol, lÅ * t;
              if nv > 0 then
                for l := 1 to nv do
                  wÆl1, lÅ := wÆl1, lÅ - wÆicol, lÅ * t;
            end   (* if l1 <> icol *)
        end
    end (* i loop *);

  if error then exit;
  (* interchange columns *)
  for i := 1 to n do
    begin
      l := n - i + 1;
      if indexÆl, 1Å <> indexÆl, 2Å then
        begin
          irow := indexÆl, 1Å;
          icol := indexÆl, 2Å;
          for k := 1 to n do
            swap(bÆk, irowÅ, bÆk, icolÅ)
        end (* if index *)
    end  (* i loop *);
  for k := 1 to n do
    if indexÆk, 3Å <> 1 then
      begin
        writeln(' ERROR: matrix singular');
        error := true;
        exit;   (* abort *)
      end;
  for i := 1 to n do
    coefÆiÅ := wÆi, 1Å;
end (* procedure gaussj *);


procedure get_data(var a : ary2s;
                   var y : arys;
                   var n, m : integer);

(* setup n-by-n hilbert matrix *)

var
  i, j : integer;

begin
  for i := 1 to n do
    begin
      aÆn,iÅ := 1.0/(n + i - 1);
      aÆi,nÅ := aÆn,iÅ
    end;
  aÆn,nÅ := 1.0/(2*n -1);
  for i := 1 to n do
    begin
      yÆiÅ := 0.0;
      for j := 1 to n do
        yÆiÅ := yÆiÅ + aÆi,jÅ
    end;
  writeln;
  if n < 7 then
    begin
      for i:= 1 to n  do
        begin
          for j:= 1 to m do
            write( aÆi,jÅ :7:5, '  ');
          writeln( ' : ', yÆiÅ :7:5)
        end;
      writeln
    end  (* if n<7 *)
end (* procedure get_data *);

procedure write_data;

(* print out the answers *)

var
  i : integer;

begin
  for i := 1 to m do
    write( coefÆiÅ :13:9);
  writeln;
end (* write_data *);


begin  (* main program *)
  aÆ1,1Å := 1.0;
  n := 2;
  m := n;
  repeat
    get_data (a, y, n, m);
    for i := 1 to n do
      for j := 1 to n do
        bÆi,jÅ := aÆi,jÅ (* setup work array *);
    gaussj (b, y, coef, n, error);
    if not error then write_data;
    n := n+1;
    m := n
  until n > maxr;
end.
«eof»
DataMuseum.dk

CP/M

⟦807b1c8a1⟧ TextFile

Derivation

TextFile
