/*
   Author Barton Willis
   University of Nebraska at Kearney
   Copyright (C) 2005, Barton Willis

   Brief Description: Maxima code for finding the nullspace and column space 
   of a matrix, along with code that triangularizes matrices using the method
   described in ``Eigenvalues by row operations,'' by Barton Willis.

   This is free software  you can redistribute it and/or
   modify it under the terms of the GNU General Public License,
   http://www.gnu.org/copyleft/gpl.html.

   This software has NO WARRANTY, not even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
*/

put('linearalgebra,1,'version);

/* The translator generates code that doesn't work when 
   translate_fast_arrays is true. So set translate_fast_arrays
   to false.
*/

eval_when([translate, compile], translate_fast_arrays : false);

eval_when(translate,
	  declare_translated(ptriangularize_with_proviso,locate_matrix_entry,hankel,toeplitz,
			     nullspace,require_integer,columnspace,rowswap,rowop,require_symbol,
	                     mat_fullunblocker,mat_trace,mat_unblocker, column_reduce,good_pivot,
                             hipow_gzeromat_norm))$

eval_when([batch,load,translate,compile],
  put(infolevel, debug, linearalgebra),    
  load("polynomialp"),
  load("mring"),
  load("lu"),
  load("linalgcholesky"),
  load("eigens-by-jacobi"),
  load("linalg-extra"),	
  load("matrixexp"),
  load("linalg-utilities"));
  
require_integer(i, pos, fn) := 
  if integerp(i) then true else 
     error("The", pos, "argument to" ,fn, "must be an integer");

require_symbol(x, pos, fn) :=
  if symbolp(x) or subvarp(x) then true else
     error("The", pos, "argument to", fn, "must be a symbol");

request_rational_matrix(m, pos, fn) :=
  if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else
    print("Some entries in the matrix are not rational numbers. The result might be wrong.");

dotproduct(a,b) := block([scalarmatrixp : true],
  require_matrix(a,"first", "dotproduct"),
  require_matrix(b,"second", "dotproduct"),
  if matrix_size(a) # matrix_size(b) or second(matrix_size(a)) # 1 then 
     error("Arguments to dotproduct must be vectors with the same size"),
  ctranspose(a) . b);

/* I use algsys instead of linsolve because it seems that linsolve
   doesn't pay attention to tellrat information.
*/

nullspace(m) := block([n, x, e, vo, %r, proviso, nc, nr, pos,
                       scalarmatrixp : false,
                       domxmxops : true],

  require_unblockedmatrix(m, "first", "nullspace"),
  nc : length(first(m)),
  nr : length(m),
  if nc = 0 or nr = 0 then 
    error("The argument to 'nullspace' must be a matrix with one or more rows and columns"),  
  
  m : ptriangularize_with_proviso(m, ?gensym()),
  proviso : second(m),
  m : first(m),

  for i : 1 thru nr do (
   if listp(pos : locate_matrix_entry(m, i, 1, i, nc, lambda([x],x # 0), 'bool)) then (
      if not(constantp(x : part(part(m, first(pos)),second(pos)))) then 
          proviso : adjoin(x # 0,proviso))),
 
  if not(emptyp(proviso)) then print("Proviso: ",proviso),
     
  x : funmake('matrix, makelist([?gensym()],i,1,nc)),
  e : lreduce('append, args(m.x)),
  e : map('rat,e),
  %r : setify(%rnum_list),
  e : first(algsys(e,(listofvars(x)))),
  %r : listify(setdifference(setify(%rnum_list), %r)),
  vo : if emptyp(%r) then [ ] else map("=",%r, makelist(0,i,1,length(%r))),
  x : subst(e, x),
  x : sort(makelist(subst(vo, subst(vi = 1, x)), vi, %r)),
  funmake('span, x));

nullity(m) := length(nullspace(m));

orthogonal_complement([v]) :=  block([sz],
  require_unblockedmatrix(first(v)," ","orthogonal_complement"),
  sz : matrix_size(first(v)),
  if second(sz) # 1 then error("Each argument must be a column vector"),
  for vi in v do (
    require_matrix(vi," ","orthogonal_complement"),
    if matrix_size(vi) # sz then error("Each vector must have the same dimension")),
  nullspace(funmake('matrix, map('first, map('transpose, v)))));
 
locate_matrix_entry(m, r1, c1, r2, c2, fn, rel) := block([im, cm, mf, e, 
     nr, nc, ok : true, frel],
  require_unblockedmatrix(m, "first", "locate_matrix_entry"),
  require_integer(r1, "second", "locate_matrix_entry"),
  require_integer(r2, "second", "locate_matrix_entry"),
  require_integer(c1, "third", "locate_matrix_entry"),
  require_integer(r2, "fourth", "locate_matrix_entry"),
  nr : length(m),
  nc : length(first(m)),
 
  if (c1 > c2) or (r1 > r2) or (r1 < 1) or (c1 < 1) or (r2 < 1) or (c2 < 1) or 
     (c1 > nc) or (c2 > nc) or (r1 > nr) or (r2 > nr) then error("Bogus submatrix"),

  mf : if rel = 'min then 'inf else if rel = 'max then 'minf else 0,	
  frel : if rel = 'max then ">" else if  rel = 'min then "<" else if rel = 'bool 
    then lambda([x,y],x) else error("Last argument must be 'max' or 'min'"),
  im : 0,
  cm : 0,
  for i : r1 thru r2 while ok do (
    for j : c1 thru c2 while ok do (
       e : apply(fn, [m[i,j]]),
       if is(apply(frel, [e, mf])) then (
         im : i,
         cm : j,
         if rel = 'bool then ok : false,
         mf : e))),
  if im > 0 and cm > 0 then [im, cm] else false);
         
columnspace(a) := block([m, nc, nr, cs : [], pos, x, proviso : [], algebraic : true],
  require_unblockedmatrix(a, "first","columnspace"),
  /* request_rational_matrix(a, "first", "columnspace"), */
  nc : length(first(a)),
  nr : length(a),
  if nc = 0 or nr = 0 then 
    error("The argument to 'columnspace' must be a matrix with one or more rows and columns"),  

  m : ptriangularize_with_proviso(a, ?gensym()),
  proviso : second(m),
  m : first(m),
  for i : 1 thru nr do (
   if listp(pos : locate_matrix_entry(m, i, 1, i, nc, lambda([x],x # 0), 'bool)) then (
      if not(constantp(x : part(part(m, first(pos)),second(pos)))) then 
          proviso : adjoin(x # 0,proviso),
       cs : cons(col(a,second(pos)),cs))),
  cs : sort(cs),
  if not(emptyp(proviso)) then print("Proviso: ",proviso),
  funmake('span,cs));

/* Maxima's rank function doesn't complain with symbolic matrices. 
I don't approve of rank(matrix([a,b],[c,d])) --> 2.  
*/

rank(m) := length(columnspace(m));

rowswap(m,i,j) := block([n, p, r],
  require_matrix(m, "first", "rowswap"),
  require_integer(i, "second", "rowswap"),
  require_integer(j, "third", "rowswap"),
  n : length(m),
  if (i < 1) or (i > n) or (j > n) then error("Array index out of bounds"),
  p : copymatrix(m),
  r : p[i],
  p[i] : p[j],
  p[j] : r,
  p);

columnswap(m,i,j) := transpose(rowswap(transpose(m),i,j));

/*  row(m,i) <-- row(m,i) - theta * row(m,i) */

rowop(m,i,j,theta) := block([p : copymatrix(m), listarith : true],
  p[i] : p[i] - theta * p[j],
  p);

/*  col(m,i) <-- col(m,i) - theta * col(m,i) */       

columnop(m,i,j,theta) := transpose(rowop(transpose(m),i,j,theta));

hipow_gzero(e,x) := block([n : hipow(e,x)], if n > 0 then n else 'inf);
good_pivot(e,x) := freeof(x,e) and e # 0;

ptriangularize(m,v) := block([p : ptriangularize_with_proviso(m,v)],
   if not(emptyp(p[2])) then print("Proviso: ",p[2]),
   p[1]); 
  
ptriangularize_with_proviso(m,v) := block([nr, nc, proviso : [], mp],
  require_unblockedmatrix(m, "first", "ptriangularize"),
  require_symbol(v, "second", "ptriangularize"),
  nr : length(m),
  nc : length(first(m)),
  for i : 1 thru min(nr, nc) do (
     mp : column_reduce(m,i,v),
     m : first(mp),
     proviso : append(proviso, second(mp))),
  proviso : setify(proviso),
  [expand(m), proviso]);

column_reduce(m,i,x) := block([nc, nr, pos, q, proviso : []],
   /* require_matrix(m, "first", "column_reduce"), */
   nc : length(first(m)),
   nr : length(m),
   /* require_integer(i, "second", "column_reduce"), 
   if (i < 1) or (i > (length(first(m)))) then error("Bogus matrix column"), */
 
   while not(good_pivot(m[i,i],x)) and (i+1 <= nr) and
        listp(pos : locate_matrix_entry(m,i+1,i,nr,i, lambda([e],e # 0), 'bool)) do (
     pos : locate_matrix_entry(m,i,i,nr,i,lambda([e], good_pivot(e,x)), 'bool),
     if listp(pos) then (
       m : rowswap(m,i,pos[1]))
     else (
       m : expand(m),
       pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], hipow(e,x)), 'max),
       m : rowswap(m,i,pos[1]),
       pos : locate_matrix_entry(m, i+1, i, nr, i, lambda([e], hipow_gzero(e,x)), 'min),
       if listp(pos) then (
         q : divide(m[i,i], m[pos[1],i],x),
         m : rowop(m,i,pos[1],first(q))))),
                
     pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], good_pivot(e,x)), 'bool),

     if listp(pos) then (
        m : rowswap(m,i,first(pos)),
        for j : i+1 thru nr do (
          if not(constantp(ratnumer(m[i,i]))) then (
            proviso : cons(ratnumer(m[i,i]) # 0, proviso)),
           if not(constantp(ratdenom(m[i,i]))) then (
            proviso : cons(ratdenom(m[i,i]) # 0, proviso)),
          m : rowop(m,j,i,m[j,i]/m[i,i])))
     else (
       pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], e # 0), 'bool),
       if listp(pos) then m : rowswap(m,i, first(pos))),
    [m, proviso]);

mat_norm(m, p) := block([listarith : true, d],
  require_matrix(m,"first","mat_norm"),
  if blockmatrixp(m) then m : mat_fullunblocker(m),
  if p = 1 then mat_norm(transpose(m), 'inf)
  else if p = 'inf then (
    d : 0,
    for ri in m do (
      d : max(d, tree_reduce("+", map('cabs, ri)))),
    d)
  else if p = 'frobenius then tree_reduce("+", lreduce('append, args(cabs(m)^2)))
  else error("Not able to compute the ",p," norm"));

mat_fullunblocker(m) := block(
  require_matrix(m, "first", "mat_unblocker"),
  if blockmatrixp(m) then (
     mat_fullunblocker(lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m)))))
   else m);

mat_unblocker(m):=block(
  require_matrix(m, "first", "mat_unblocker"),
  if blockmatrixp(m) then (
     lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m))))     
   else m);

mat_trace(m) := block([n, acc : 0],
  if not matrixp(m) then funmake('mat_trace, [m]) else (
    n : matrix_size(m),
    if first(n) # second(n) then error("The first argument to 'mat_trace' must be a square matrix"),
    n : first(n),
    for i from 1 thru n do (
      acc : acc + if matrixp(m[i,i]) then mat_trace(m[i,i]) else m[i,i]),
    acc));

kronecker_product(a,b) := (
  require_matrix(a,"first", "kronecker_product"),
  require_matrix(b,"second", "kronecker_product"),
  mat_unblocker(outermap('matrix_element_mult, a,b)));
  
diag_matrix([d]) := genmatrix(lambda([i,j], if i = j then inpart(d,i) else zerofor(inpart(d,i))),
                                      length(d),length(d));

hilbert_matrix(n) := (
  require_posinteger(n,"first","hilbert_matrix"),
  funmake('matrix, makelist(makelist(1/(i + j - 1),i,1, n),j, 1, n)));

/* To avoid a "too many arguments" error, use use funmake('matrix, ....) instead of
   apply('matrix, ...).
*/

vandermonde_matrix(vars) := block([n],
  require_list(vars,"first","vandermonde_matrix"),
  n : length(vars),
  funmake('matrix, makelist(makelist(vars[i]^j,j,0,n-1),i,1,n)));

hankel([q]) := block([col,row,m,n,partswitch : false],
  if length(q) > 2 or length(q) < 1 then error("The function 'hankel' requires one or two arguments"),
  col : inpart(q,1),
  row : if length(q) = 2 then inpart(q,2) else map(lambda([x],0), col),

  require_list(row,"first","hankel"),
  require_list(col,"second","hankel"),  
  m : length(row),
  n : length(col),
  genmatrix(lambda([i,j],if i+j-1 <= n then inpart(col,i+j-1) else inpart(row,i+j-n)),n,m));

toeplitz([q]) := block([col,row,m,n,partswitch : false],
  if length(q) > 2 or length(q) < 1 then error("The function 'toeplitz' requires one or two arguments"),
  col : inpart(q,1),
  row : if length(q) = 2 then inpart(q,2) else map('conjugate, col),

  require_list(row,"first","toeplitz"),
  require_list(col,"second","toeplitz"),  
  m : length(row),
  n : length(col),
  genmatrix(lambda([i,j], if i -j >= 0 then inpart(col, i-j+1) else inpart(row,j-i+1)),n,m));

polytocompanion(p,x) := block([n],
  if not polynomialp(p,[x], lambda([e], freeof(x,e))) then 
     error("First argument to 'polytocompanion' must be a polynomial"),
  p : expand(p),
  n : hipow(p,x),
  p : multthru(p / coeff(p,x,n)),
  genmatrix(lambda([i,j], if i-j=1 then 1 else if j = n then -coeff(p,x,i-1) else 0),n));