⟦ff3ea2c50⟧

TextFile

program HilbDemo87;
æ
          TURBO-87 DEMONSTRATION PROGRAM  Version 1.00A

      This program demonstrates the increased speed and precision
      of the TURBO-87 compiler:

      --------------------------------------------------
      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
      --------------------------------------------------

      The program performs simultaneous solution by Gauss-Jordan
      elimination.

      INSTRUCTIONS
      1.  Compile the program using the TURBO-87.COM compiler.
      2.  Type Ctrl-C to interrupt the program.
å

CONST
  maxr = 10;
  maxc = 10;

TYPE
  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 *)

LABEL 
  99,98;

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;
                      GOTO 98        (* 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 *);
  98:

  IF error THEN GOTO 99;
  (* 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;
        GOTO 99   (* abort *)
      END;
  FOR i := 1 TO n DO
    coefÆiÅ := wÆi, 1Å;
  99:
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

⟦ff3ea2c50⟧ TextFile

Derivation

TextFile
