##
## "solving" inconsistent linear systems of equations
## in the infinity-norm
##
## Input: real mxn matrix A, where m >= n,
##        and a m-dimensional vector b of real constants.
##
## Output: n-dimensional vector x, such that ||Ax - b|| = min
##         where ||.|| is the infinity-norm
##

incsolve := proc (A,b)

  local m, n, i, j, c, d, w, x, s, ss;

#
# validity checks
#

  if type (A, `matrix`) = false then
    ERROR (`first parameter must be of type matrix`);
  fi;

  if type (b, `vector`) = false then
    ERROR (`second parameter must be of type vector`);
  fi;

  m := linalg[rowdim](A);
  n := linalg[coldim](A);

  if m < n then
    ERROR (`the number of rows must be > or = the number of columns`);
  fi;

  if m <> linalg[vectdim](b) then
    ERROR (`incompatible dimensions`);
  fi;

#
# set up local vectors
#

  x := linalg[vector](n);
  d := linalg[vector](m);

  c := {};   # set of constraints

#
# generate linear equations, and constraints
#

  for i from 1 to m do

    d[i] := -b[i];

    for j from 1 to n do

      d[i] := d[i] + A[i,j]*x[j];

    od;

    c := c union { w>=d[i], w>=-d[i] };

  od;

#
# solve linear program
#

  ss := simplex[minimize](w,c);

#
# assign solution to output vector
#

  for i from 1 to n+1 do

    s := op (i,ss);

    if op (1,s) = 'w' then
      print (s);
    else
      assign (s);
    fi;

  od;

  RETURN (evalm(x));

end;

incsolve_weighted := proc (A,b)

  local m, n, i, j, c, d, w, x, s, ss;

#
# validity checks
#

  if type (A, `matrix`) = false then
    ERROR (`first parameter must be of type matrix`);
  fi;

  if type (b, `vector`) = false then
    ERROR (`second parameter must be of type vector`);
  fi;

  m := linalg[rowdim](A);
  n := linalg[coldim](A);

  if m < n then
    ERROR (`the number of rows must be > or = the number of columns`);
  fi;

  if m <> linalg[vectdim](b) then
    ERROR (`incompatible dimensions`);
  fi;

#
# set up local vectors
#

  x := linalg[vector](n);
  d := linalg[vector](m);

  c := {};   # set of constraints

#
# generate linear equations, and constraints
#

  for i from 1 to m do

    d[i] := -b[i];

    for j from 1 to n do

      d[i] := d[i] + A[i,j]*x[j];

    od;

    c := c union { abs(b[i])*w>=d[i], abs(b[i])*w>=-d[i] };

  od;

#
# solve linear program
#

  ss := simplex[minimize](w,c);

#
# assign solution to output vector
#

  for i from 1 to n+1 do

    s := op (i,ss);

    if op (1,s) = 'w' then
      print (s);
    else
      assign (s);
    fi;

  od;

  RETURN (evalm(x));

end;