segment QMATRIX;   (* 6/20/1989 *)
  (* The quick matrix routines collection. *)
  (* minimal or no error checking is done *)
(******************************************)
(* INTERFACE -> *)
%include QMATRIX;
(******************************************)
(* IMPLEMENTATION -> *)
const
  MatMsg1 = 'Arguments out of range';
  MatMsg2 = 'Matrices not conformable';
  MatMsg3 = 'Invalid matrix dimension';
  MatMsg4 = 'EOF reached';
  MatMsg5 = 'Singular matrix';
  MatMsg6 = 'Negative definite matrix';
(************************)
procedure ZeroOut;
 (* (var m : matrix;  row1, row2, col1, col2 : integer);  *)

  (* Zero-out a sumbatrix *)

var row, col : integer;

Begin (* ZeroOut *)
  with m do
  begin
    for row := row1 to row2 do
      for col := col1 to col2 do
        el(.row, col.) := 0;
  end;
End;  (* ZeroOut *)
(******************************************)
procedure MatL2U;
 (* reflect lower half to upper half *)
  (* matrix must be square *)
const
  myname = 'MatL2U: ';
var  i, j : integer;
Begin (* MatL2U *)
  with a do
  begin
    for i := 2 to rows do
    for j := 1 to pred(i) do el(.j,i.) := el(.i,j.);
  end;
End; (* MatL2U *)
(******************************************)
procedure MatU2L; (* reflect upper half to lower half *)
  (* matrix must be square *)
const
  myname = 'MatU2L: ';
var i, j : integer;
Begin (* MatU2L *)
  with a do
  begin
    for i := 2 to rows do
    for j := 1 to pred(i) do el(.i,j.) := el(.j,i.);
  end;
End; (* MatU2L *)
(******************************************)
procedure MatDim;
const
  myname = 'MatDim: ';
Begin
  if (r > MAXROWS) or (c > MAXCOLS) or (r < 1) or (c < 1) then
  begin
    writeln (StdErr, myname, MatMsg1);
    writeln (StdErr, 'r = ', r : 4, '  c = ', c : 4);
    MatFail := TRUE;
    RETURN;
  end;
  MatFail := FALSE;
  a.rows := r;
  a.cols := c;
End;
(******************************************)
procedure MatRead;
const
  myname = 'MatRead: ';
var
  i, j : integer;
Begin
  with a do
  for i := 1 to rows do for j := 1 to cols do
  begin
    if EOF(f) then
    begin
      writeln(StdErr, myname, MatMsg4);
      MatFail := TRUE;
      RETURN;
    end;
    read (f, el(.i,j.));
  end;
  MatFail := FALSE;
End;
(******************************************)
procedure MatWrite;
const
  myname = 'MatWrite: ';
  SP = ' ';
  CPP = 12; (* max columns per page *)
var
  i, j, column : integer;
Begin
  with a do
  for i := 1 to rows do
  begin
    column := 0;
    write ('(' : 1, i : 2, ') ' : 2 );
    for j := 1 to cols do
    begin
      column := succ(column);  (* bump *)
      if (column > CPP) then   (* end of line? *)
      begin
        column := 1;           (* new line *)
        writeln (f);
        write (f, '':5);
      end;
      write (f, el(.i,j.)  : 8 : 5, '' : 2); (* at column *)
    end;
    writeln (f);
  end;
  writeln (f);
  MatFail := FALSE;
End;
(******************************************)
procedure MatMul;
   (* a := b * c *)
const
  myname = 'MatMul: ';
var
  i, j, k : integer;
  s : real;
Begin
  if (b.cols <> c.rows) then
  begin
    writeln (StdErr, myname, MatMsg2);
    MatFail := TRUE;
    RETURN;
  end;
  MatDim (a, b.rows, c.cols);
  for i := 1 to b.rows do
  for j := 1 to c.cols do
  begin
    s := 0;
    for k := 1 to b.cols do
       s := s + b.el(.i, k.) * c.el(.k, j.);
    a.el(.i, j.) :=  s;
  end;
End;
(******************************************)
procedure MatAdd;
   (* a := b + c *)
const
  myname = 'MatAdd: ';
var
  i, j : integer;
Begin
  if not ((b.rows = c.rows) and (b.cols = c.cols)) then
  begin
    writeln (StdErr, myname, MatMsg2);
    MatFail := TRUE;
    RETURN;
  end;
  MatDim (a, b.rows, b.cols);
  for i := 1 to b.rows do
  for j := 1 to b.cols do
  begin
    a.el(.i,j.) :=  b.el(.i, j.) + c.el(.i, j.);
  end;
End;
(******************************************)
procedure MatSub;
   (* a := b - c *)
const
  myname = 'MatSub: ';
var
  i, j : integer;
Begin
  if not ((b.rows = c.rows) and (b.cols = c.cols)) then
  begin
    writeln (StdErr, myname, MatMsg2);
    MatFail := TRUE;
    RETURN;
  end;
  MatDim (a, b.rows, b.cols);
  for i := 1 to b.rows do
  for j := 1 to b.cols do
  begin
    a.el(.i,j.) :=  b.el(.i, j.) - c.el(.i, j.);
  end;
End;
(******************************************)
procedure MatScale;
   (* a := b * c *)
const
  myname = 'MatScale: ';
var
  i, j : integer;
Begin
  MatDim (a, b.rows, b.cols);
  for i := 1 to b.rows do
  for j := 1 to b.cols do
  begin
    a.el(.i, j.) := b.el(.i, j.) * c;
  end;
End;
(******************************************)
procedure MatTran;
  (* a := b' *)
const
  myname = 'MatTran: ';
var
  i, j : integer;
Begin
  MatDim (a, b.cols, b.rows);
  with a do
  for i := 1 to b.rows do
  for j := 1 to b.cols do el(.j, i.) := b.el(.i, j.);
End;
(******************************************)
procedure MatFill;
var
  i, j : integer;
Begin
  with a do
  for i := 1 to rows do
  for j := 1 to cols do
  begin
    el(.i, j.) :=  b;
  end;
  MatFail := FALSE;
End;

(* ============================================================= *)
(*
  Invert a symmetric positive definite matrix in situ
  using Cholesky  decomposition
  Translated for the matrix procedures. 6/20/1989
*)
(* ============================================================= *)
(* By  Bahjat Qaqish               2/8 /1988                     *)
(* procs: MatSInv, MatChol, MatInvL, MatSInv, MatSHInv   *)
(* all what is needed to invert a symmetric positive definite    *)
(* matrix in situ                                                *)
(* ============================================================= *)
procedure MatChol ;
(*
       Cholesky decomposition of a symmetric positive definite
       matrix A = L L' . L is returned in the
       lower triangle of A, upper A not accessed.
*)
const
  epsilon = 1.0E-9;
  myname = 'MatChol: ';
var
    k, i, j, n :integer ;
    s, w       : real ;

begin  (* MatChol *)
   n := a.rows;
   with a do
   for k := 1 to n   do
   begin ;   (* k-loop *)
     for i := 1 to pred(k) do
     begin ;     (* i-loop *)
       s := 0   ;
       for j := 1 to pred(i)  do
         s := s + el(.i, j.) * el(.k, j.);
         w := el(.i,i.);
         if (w < 0) then
         begin
           MatFail := TRUE;
           writeln (StdErr, myname, MatMsg6);
           RETURN;
         end;
         if (w < epsilon) then
         begin
           MatFail := TRUE;
           writeln (StdErr, myname, MatMsg5);
           RETURN;
        end;
       el(.k, i.) := (el(.k, i.)-s) / w;
     end ;       (* i-loop *)
     s := 0;
     for j := 1 to k-1 do s := s + sqr (el(.k,j.));
     w := el(.k,k.) - s;
     if (w < 0) then
     begin
       MatFail := TRUE;
       writeln (StdErr, myname, MatMsg6);
       RETURN;
     end;
     el(.k, k.) :=  sqrt (w);
   end ; (* k-loop *)
   MatFail := FALSE;
end ; (*  MatChol  *)
(* ============================================================= *)
procedure MatInvL ;
(*
  by Bahjat Qaqish 1/21/1987
  invert a lower triangular matrix L in situ
  assumes L is invertible. Upper L not changed
*)
const
  epsilon = 1.0E-9;
  myname = 'MatInvL: ';
var
  i,j,p, n : integer ;
  s,w   : real ;

begin (* MatInvL  *)
  n := L.rows;
  with L do
  for i:= 1 to n do
  begin
    for j:=1 to pred(i) do
    begin
      s:=0 ;
      for p:=j to pred(i) do
        s := s + el(.i,p.)*el(.p,j.) ;
      w := el(.i,i.);
      if (abs(w) < epsilon) then
      begin
        MatFail := TRUE;
        writeln (StdErr, myname, MatMsg5);
        RETURN;
      end;
      el(.i,j.) := -s / w;
    end ;
    w := el(.i,i.);
    if (abs(w) < epsilon) then
    begin
      MatFail := TRUE;
      writeln (StdErr, myname, MatMsg5);
      RETURN;
    end;
    el(.i,i.) := 1 / w;
  end ;
  MatFail := FALSE;
end ; (* MatInvL  *)
(* ============================================================= *)
procedure MatLpL ;
(*
  by Bahjat Qaqish 1/21/1987
  compute L'L , where L is lower triangular
  result is returned in L
*)
var
  i,j,k: integer ;
  s    : real ;

begin (* MatLpL *)
  with L do
  (* first compute lower triangle of LPL *)
  for i:= 1 to rows do
  for j:= 1 to i do
  begin
    s:=0 ;
    for k := i to rows do s:=s+el(.k,i.)*el(.k,j.);
    el(.i,j.) := s;
  end ;
  (* next miror image into upper triangle *)
  MatL2U (L);
end ; (* MatLpL *)
(* ============================================================= *)
procedure MatSInv;
(*
  by Bahjat Qaqish    1/21/1987
  inverts a symmetric positive definite matrix A in situ
*)
begin  (* MatSInv *)
  MatChol    (A);  (* A =LL' (L in A ) *)
  if MatFail then RETURN;
  MatInvL  (A);  (* A <- inverse of A *)
  if MatFail then RETURN;
  MatLpL     (A);  (* A <- A'A *)
end ;  (* MatSInv *)
(* ============================================================= *)
procedure MatSHInv ;
(*
  by Bahjat Qaqish    1/21/1987
  returns the inverse of the Cholesky root of A  into lower A
*)
begin  (* MatSHInv *)
  MatChol   (A) ;  (* A    <- lower Cholesky root A *)
  MatInvL (A) ;  (* A    <- inverse of A *)
end ;  (* MatSHInv *)
(* ============================================================= *)
(************************)
procedure blockout;
var
  row, col : integer;
Begin
  with a do
  begin
    for row := r1 to r2 do
    for col := c1 to c2 do
    el(.row, col.) := 0
  end
End;
(************************)
procedure VecWrite;
var
  i : integer;

Begin (* VecWrite *)
  with v do
  for i := 1 to rows do
  writeln (f, i:8, '':2, ve(.i.) : 8 : 5 );
End;  (* VecWrite *)
(************************)
procedure MatxVec;
        (* b <- A * x *)
        (* must: A.cols == x.rows *)
var
  i, j : integer;
  w : real;
Begin (* MatxVec *)
  if (A.cols <> x.rows) then
  begin
    writeln(StdErr, 'MatxVec: nonconformable dimensions.');
    HALT;
  end;
  b.rows := A.rows;
  with A do
  for i := 1 to rows do
  begin
    w := 0;
    for j := 1 to cols do w := w + el(.i,j.) * x.ve(.j.);
    b.ve(.i.) := w;
  end;
End; (* MatxVec *)
(************************)
procedure VecFill;
var
  i : integer;
Begin (* VecFill *)
  with x do
    for i := 1 to rows do ve(.i.) := r;
End; (* VecFill *)
(************************)
function  VecSumAbs;
var
  i : integer;
  w : real;
Begin (* VecSumAbs *)
  w := 0;
  with x do
    for i := 1 to rows do w := w + abs(ve(.i.));
  VecSumAbs  := w;
End; (* VecSumAbs *)
(************************)
procedure AddOP;
      (* m  <- m  + outer product of v *)
var
  i, j : integer;  c : real;

Begin (*  AddOP *)
  if (m.rows <> m.cols) or (m.rows <> v.rows) then
  begin
    writeln(StdErr, 'AddOP: Nonconformable');
    HALT;
  end;
  with m  do
  begin
    for i := 1 to rows do el(.i,i.) := el(.i,i.) + sqr(v.ve(.i.));
    for i := 1 to pred(rows) do
    for j := succ(i) to rows do
    begin
      c := v.ve(.i.)*v.ve(.j.);
      el(.i,j.) := el(.i,j.) + c;
      el(.j,i.) := el(.j,i.) + c;
    end;
  end;
End;  (*  AddOP *)
(************************)
procedure  VecAdd;
var
  i : integer;

Begin (*   VecAdd *)
  if (v2.rows <> v3.rows) then
  begin
    writeln(StdErr, 'VecAdd: Nonconformable');
    HALT;
  end;
  v1.rows := v2.rows;
  with v1 do
  for i := 1 to rows do ve(.i.) := v2.ve(.i.) + v3.ve(.i.);
End;  (*   VecAdd *)
(************************)

(*======================*)
procedure ludcmp(var a    : marray;
                     n    : integer;
                 var indx : ivarray;
                 var d    : real;
                 var fail : BOOLEAN);
(* LU decomposition. Source: Numerical Recipes *)
const
   tiny=1.0e-20;
var
   k,j,imax,i: integer;
   sum,dum,big: real;
   vv: varray;
Begin  (* ludcmp *)
   fail := FALSE;
   d := 1.0;
   for i := 1 to n do begin
      big := 0.0;
      for j := 1 to n do
        if (abs(a(.i,j.)) > big) then big := abs(a(.i,j.));
      if (big = 0.0) then begin
         fail := TRUE; RETURN;
      end;
      vv(.i.) := 1.0/big
   end;
   for j := 1 to n do begin
      if (j > 1) then begin
         for i := 1 to j-1 do begin
            sum := a(.i,j.);
            if (i > 1) then begin
               for k := 1 to i-1 do begin
                  sum := sum-a(.i,k.)*a(.k,j.)
               end;
               a(.i,j.) := sum
            end
         end
      end;
      big := 0.0;
      for i := j to n do begin
         sum := a(.i,j.);
         if (j > 1) then begin
            for k := 1 to j-1 do begin
               sum := sum-a(.i,k.)*a(.k,j.)
            end;
            a(.i,j.) := sum
         end;
         dum := vv(.i.)*abs(sum);
         if (dum > big) then begin
            big := dum;
            imax := i
         end
      end;
      if (j <> imax) then begin
         for k := 1 to n do begin
            dum := a(.imax,k.);
            a(.imax,k.) := a(.j,k.);
            a(.j,k.) := dum
         end;
         d := -d;
         vv(.imax.) := vv(.j.)
      end;
      indx(.j.) := imax;
      if (j <> n) then begin
         if (a(.j,j.) = 0.0) then a(.j,j.) := tiny;
         dum := 1.0/a(.j,j.);
         for i := j+1 to n do begin
            a(.i,j.) := a(.i,j.)*dum
         end
      end
   end;
   if (a(.n,n.) = 0.0) then a(.n,n.) := tiny
End;   (* ludcmp *)
(**********************)
procedure lubksb(var a     : marray;
                 n         : integer;
                 var indx  : ivarray;
                 var b     : varray);
(* LU back-substitution. Source: Numerical Recipes *)
var
   j,ip,ii,i: integer;
   sum: real;
Begin  (* lubksb *)
   ii := 0;
   for i := 1 to n do begin
      ip := indx(.i.);
      sum := b(.ip.);
      b(.ip.) := b(.i.);
      if  (ii <> 0) then begin
         for j := ii to i-1 do begin
            sum := sum-a(.i,j.)*b(.j.)
         end
      end else if (sum <> 0.0) then begin
         ii := i
      end;
      b(.i.) := sum
   end;
   for i := n downto 1 do begin
      sum := b(.i.);
      if (i < n) then begin
         for j := i+1 to n do begin
            sum := sum-a(.i,j.)*b(.j.)
         end
      end;
      b(.i.) := sum/a(.i,i.)
   end
End;   (* lubksb *)
(*======================*)
procedure Invert;  /* (var ainv, a1 : matrix);    */
  (* ainv <- inverse(a) *)
const
  msg1   = 'Invert: Matrix is not square';

var fail : BOOLEAN;
    i, j, n : integer;
    x, b : vector;
    indx : ivarray;
    a : matrix;
    sign : real;

Begin (* Invert *)
  n := a1.rows;
  if (n <> a1.cols) then
  begin
    writeln(StdErr, msg1);  HALT;
  end;
  ainv.rows := n;
  ainv.cols := n;

  a := a1;  (* copy a1 *)
  ludcmp (a.el, n, indx, sign, fail);
  if fail then begin
    writeln ('Invert: Singular matrix.');HALT  end;

  b.rows := n;
  i := 1;    (* first column *)
  b.ve(.i.) := 1;
  for j := 2 to n do b.ve(.j.) := 0;
  x := b;
  lubksb (a.el, n, indx, x.ve);
  for j := 1 to n do ainv.el(.j,i.) :=  x.ve(.j.);

  for i := 2 to n do     (* the rest *)
  begin
    b.ve(.pred(i).) := 0;
    b.ve(.i.) := 1;
    x := b;
    lubksb (a.el, n, indx, x.ve);
    for j := 1 to n do ainv.el(.j,i.) :=  x.ve(.j.);
  end;

End;  (* Invert *)
(*======================*)
.   (* end of segment *)
(*----- QMATRIX ends   here -----------*)