type
  opr_t = procedure(var b:double; a:double);

// row/col step... hence apparent duplication along row or col.
// only dim=1 can be expanded and the other dim must be 1 (a is
// a scalar) or the same as b i.e. a must be either a scalar or a
// row vector of the same size as a row of b or a column vector of
// the same size as a column of b.
// NOTE if the argument is a left value in add() mul()...
// then this implement the summation (product) operator over rows or lines
{}function xpd(a:mtx; b:mtx):mtx;
begin
  xpd.nam := a.nam;
  xpd.dc  := a.dc;
  xpd.dl  := a.dl;
  xpd.e   := a.e;
  xpd.col := b.col;
  xpd.lgn := b.lgn;
  if (lgn(a)=1) and (col(a)=col(b)) then begin
    xpd.dc := 0
  end else if (col(a)=1) and (lgn(a)=lgn(b)) then begin
    xpd.dl := 0
  end else if (lgn(a)=1) and (col(a)=1) then begin
    xpd.dc := 0;
    xpd.dl := 0
  end else assert((lgn(a)=col(a)) and (lgn(b)=col(b)));
  xpd.smx := true
end;

// generic operator used to implement arithmetic operations: + - .* ./
{}procedure opr(a:mtx; b:mtx; o:opr_t);
var
  i, j, nl, nc, ia, ib, ia1, ib1: longint;
begin
  if (a.lgn<b.lgn)or(a.col<b.col) then begin
    opr(xpd(a,b),b,o)
  end else if (a.lgn>b.lgn)or(a.col>b.col) then begin
    opr(a,xpd(b,a),o)
    //opr(xpd(b,a),a,o)  // swap a and b
  end else begin
    nl := a.lgn;
    nc := a.col;
    assert(nl=b.lgn);
    assert(nc=b.col);
    ia1 := 1;
    ib1 := 1;
    for i := 1 to nl do begin
      ia := ia1;
      ib := ib1;
      for j := 1 to nc do begin
        o(a.e[ia],b.e[ib]);
        ia := ia + a.dl;
        ib := ib + b.dl
      end;
      ia1 := ia1 + a.dc;
      ib1 := ib1 + b.dc
    end
  end
end;

// C style operators b += a for mtx with soft dimensions
{}procedure opradd(var b:double; a:double);
begin
  b := b + a
end;

// C style operators b -= a for mtx with soft dimensions
{}procedure oprsub(var b:double; a:double);
begin
  b := b - a
end;

// C style operators b *= a for mtx with soft dimensions
{}procedure oprmul(var b:double; a:double);
begin
  b := b * a
end;

// C style operators b /= a for mtx with soft dimensions
{}procedure oprdvd(var b:double; a:double);
begin
  b := b / a
end;
