1/*
2   Author Barton Willis
3   University of Nebraska at Kearney
4   Copyright (C) 2005, Barton Willis
5
6   Brief Description: Maxima code for finding the nullspace and column space
7   of a matrix, along with code that triangularizes matrices using the method
8   described in ``Eigenvalues by row operations,'' by Barton Willis.
9
10   This is free software  you can redistribute it and/or
11   modify it under the terms of the GNU General Public License,
12   http://www.gnu.org/copyleft/gpl.html.
13
14   This software has NO WARRANTY, not even the implied warranty of
15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
16*/
17
18put('linearalgebra,1,'version);
19
20/* The translator generates code that doesn't work when
21   translate_fast_arrays is true. So set translate_fast_arrays
22   to false.
23*/
24
25eval_when([translate, compile], translate_fast_arrays : false);
26
27eval_when(translate,
28	  declare_translated(ptriangularize_with_proviso,locate_matrix_entry,hankel,toeplitz,
29			     nullspace,require_integer,columnspace,rowswap,rowop,require_symbol,
30	                     mat_fullunblocker,mat_trace,mat_unblocker, column_reduce,good_pivot,
31                             hipow_gzero, mat_norm))$
32
33eval_when([batch,load,translate,compile],
34  put('infolevel, 'silent, 'linearalgebra),
35  load("load-linearalgebra-lisp-files"));
36
37require_integer(i, pos, fn) :=
38  if integerp(i) then true else
39     error("The", pos, "argument to" ,fn, "must be an integer");
40
41require_symbol(x, pos, fn) :=
42  if symbolp(x) or subvarp(x) then true else
43     error("The", pos, "argument to", fn, "must be a symbol");
44
45request_rational_matrix(m, pos, fn) :=
46  if every('identity, map(lambda([s], every('ratnump,s)), args(m))) then true else
47    print("Some entries in the matrix are not rational numbers. The result might be wrong.");
48
49dotproduct(a,b) := block([scalarmatrixp : true],
50  require_matrix(a,"first", "dotproduct"),
51  require_matrix(b,"second", "dotproduct"),
52  if matrix_size(a) # matrix_size(b) or second(matrix_size(a)) # 1 then
53     error("Arguments to dotproduct must be vectors with the same size"),
54  ctranspose(a) . b);
55
56nullspace(m) := block([nr, nc, acc : set(), proviso : true, pv, prederror : false],
57
58  require_unblockedmatrix(m, "first", "nullspace"),
59  /* nc and nr are the sizes of the transpose of m */
60
61   [nc, nr] : matrix_size(m),
62   m : triangularize(addcol(transpose(m), ident(nr))),
63   for row : 1 thru nr do (
64     pv : locate_matrix_entry(m, row, 1, row, nc, lambda([s], compare(s,0) # "="), 'bool),
65
66     if not(listp(pv)) then (
67       acc : adjoin(transpose(genmatrix(lambda([ii,j], m[row,j+nc]),1,nr)),acc))
68     else (
69       proviso : proviso and notequal(m[row, second(pv)],0))),
70
71   if proviso # true then (
72     print("Proviso: ",proviso),
73     put(concat(outchar, linenum), proviso, 'proviso)),
74
75   subst('span, 'set, subset(acc, lambda([s], some(lambda([x], compare(x,0) # "="), s)))));
76
77
78nullity(m) := length(nullspace(m));
79
80orthogonal_complement([v]) :=  block([sz],
81  require_unblockedmatrix(first(v)," ","orthogonal_complement"),
82  sz : matrix_size(first(v)),
83  if second(sz) # 1 then error("Each argument must be a column vector"),
84  for vi in v do (
85    require_matrix(vi," ","orthogonal_complement"),
86    if matrix_size(vi) # sz then error("Each vector must have the same dimension")),
87  nullspace(transpose(rreduce('addcol, v))));
88
89locate_matrix_entry(m, r1, c1, r2, c2, fn, rel) := block([im, cm, mf, e,
90     nr, nc, ok : true, frel],
91  require_unblockedmatrix(m, "first", "locate_matrix_entry"),
92  require_integer(r1, "second", "locate_matrix_entry"),
93  require_integer(r2, "second", "locate_matrix_entry"),
94  require_integer(c1, "third", "locate_matrix_entry"),
95  require_integer(r2, "fourth", "locate_matrix_entry"),
96  nr : length(m),
97  nc : length(first(m)),
98
99  if (c1 > c2) or (r1 > r2) or (r1 < 1) or (c1 < 1) or (r2 < 1) or (c2 < 1) or
100     (c1 > nc) or (c2 > nc) or (r1 > nr) or (r2 > nr) then error("Bogus submatrix"),
101
102  mf : if rel = 'min then 'inf else if rel = 'max then 'minf else 0,
103  frel : if rel = 'max then ">" else if  rel = 'min then "<" else if rel = 'bool
104    then lambda([x,y],x) else error("Last argument must be 'max' or 'min'"),
105  im : 0,
106  cm : 0,
107  for i : r1 thru r2 while ok do (
108    for j : c1 thru c2 while ok do (
109       e : apply(fn, [m[i,j]]),
110       if is(apply(frel, [e, mf])) then (
111         im : i,
112         cm : j,
113         if rel = 'bool then ok : false,
114         mf : e))),
115  if im > 0 and cm > 0 then [im, cm] else false);
116
117columnspace(a) := block([m, nc, nr, cs : set(), pos, x, proviso : true, algebraic : true, c_min],
118  require_unblockedmatrix(a, "first","columnspace"),
119  [nr, nc] : matrix_size(a),
120  if nc = 0 or nr = 0 then (
121    error("The argument to 'columnspace' must be a matrix with one or more rows and columns")),
122
123  m : triangularize(a),
124  c_min : 1,
125  for i : 1 thru nr while c_min <= nc do (
126   if listp(pos : locate_matrix_entry(m, i, c_min, i, nc, lambda([x], x # 0), 'bool)) then (
127     c_min : second(pos) + 1,
128     if constantp(x : part(part(m, first(pos)),second(pos))) = false then (
129       proviso : proviso and notequal(ratsimp(x), 0)),
130     cs : adjoin(col(a,second(pos)) ,cs))),
131
132  if proviso # true then (
133    print("Proviso: ",proviso),
134    put(concat(outchar, linenum), proviso, 'proviso)),
135  funmake('span, listify(cs)));
136
137/* Maxima's rank function doesn't complain with symbolic matrices.
138I don't approve of rank(matrix([a,b],[c,d])) --> 2. This version will
139print the proviso warnings.
140*/
141
142linalg_rank(m) := length(columnspace(m));
143
144rowswap(m,i,j) := block([n, p, r],
145  local(p),
146  require_matrix(m, "first", "rowswap"),
147  require_integer(i, "second", "rowswap"),
148  require_integer(j, "third", "rowswap"),
149  n : length(m),
150  if (i < 1) or (i > n) or (j > n) then error("Array index out of bounds"),
151  p : copymatrix(m),
152  r : p[i],
153  p[i] : p[j],
154  p[j] : r,
155  p);
156
157columnswap(m,i,j) := transpose(rowswap(transpose(m),i,j));
158
159/*  row(m,i) <-- row(m,i) - theta * row(m,i) */
160
161rowop(m,i,j,theta) := block([p : copymatrix(m), listarith : true],
162  local(p),
163  p[i] : p[i] - theta * p[j],
164  p);
165
166/*  col(m,i) <-- col(m,i) - theta * col(m,i) */
167
168columnop(m,i,j, theta) := transpose(rowop(transpose(m),i, j, theta));
169
170hipow_gzero(e,x) := block([n : hipow(e,x)], if n > 0 then n else 'inf);
171good_pivot(e,x) := freeof(x,e) and e # 0;
172
173ptriangularize(m,v) := block([p : ptriangularize_with_proviso(m,v)],
174   if not(emptyp(p[2])) then print("Proviso: ",p[2]),
175   p[1]);
176
177ptriangularize_with_proviso(m,v) := block([nr, nc, proviso : [], mp],
178  require_unblockedmatrix(m, "first", "ptriangularize"),
179  require_symbol(v, "second", "ptriangularize"),
180  nr : length(m),
181  nc : length(first(m)),
182  for i : 1 thru min(nr, nc) do (
183     mp : column_reduce(m,i,v),
184     m : first(mp),
185     proviso : append(proviso, second(mp))),
186  proviso : setify(proviso),
187  [expand(m), ratsimp(proviso)]);
188
189column_reduce(m,i,x) := block([nc, nr, pos, q, proviso : []],
190   /* require_matrix(m, "first", "column_reduce"), */
191   nc : length(first(m)),
192   nr : length(m),
193   /* require_integer(i, "second", "column_reduce"),
194   if (i < 1) or (i > (length(first(m)))) then error("Bogus matrix column"), */
195
196   while not(good_pivot(m[i,i],x)) and (i+1 <= nr) and
197        listp(pos : locate_matrix_entry(m,i+1,i,nr,i, lambda([e],e # 0), 'bool)) do (
198     pos : locate_matrix_entry(m,i,i,nr,i,lambda([e], good_pivot(e,x)), 'bool),
199     if listp(pos) then (
200       m : rowswap(m,i,pos[1]))
201     else (
202       m : expand(m),
203       pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], hipow(e,x)), 'max),
204       m : rowswap(m,i,pos[1]),
205       pos : locate_matrix_entry(m, i+1, i, nr, i, lambda([e], hipow_gzero(e,x)), 'min),
206       if listp(pos) then (
207         q : divide(m[i,i], m[pos[1],i],x),
208         m : rowop(m,i,pos[1],first(q))))),
209
210     pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], good_pivot(e,x)), 'bool),
211
212     if listp(pos) then (
213        m : rowswap(m,i,first(pos)),
214        for j : i+1 thru nr do (
215          if not(constantp(ratnumer(m[i,i]))) then (
216            proviso : cons(ratnumer(m[i,i]) # 0, proviso)),
217           if not(constantp(ratdenom(m[i,i]))) then (
218            proviso : cons(ratdenom(m[i,i]) # 0, proviso)),
219          m : rowop(m,j,i,m[j,i]/m[i,i])))
220     else (
221       pos : locate_matrix_entry(m, i, i, nr, i, lambda([e], e # 0), 'bool),
222       if listp(pos) then m : rowswap(m,i, first(pos))),
223    [m, proviso]);
224
225mat_norm(m, p) := block([listarith : true],
226  require_matrix(m,"first","mat_norm"),
227  if blockmatrixp(m) then m : mat_fullunblocker(m),
228  if p = 1 then lmax(xreduce("+", map(lambda([s], map('cabs, s)), if args(m) = [] then [[]] else args(m))))
229  else if p = 'inf then lmax(map(lambda([s], xreduce("+", map('cabs, s))), args(m)))
230  else if p = 'frobenius then sqrt(xreduce("+", lreduce('append, args(cabs(m)^2))))
231  else error("Not able to compute the ",p," norm"));
232
233mat_fullunblocker(m) := block(
234  require_matrix(m, "first", "mat_unblocker"),
235  if blockmatrixp(m) then (
236     mat_fullunblocker(lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m)))))
237   else m);
238
239mat_unblocker(m):=block(
240  require_matrix(m, "first", "mat_unblocker"),
241  if blockmatrixp(m) then (
242     lreduce('addrow, args(map(lambda([x], lreduce('addcol, x)), m))))
243   else m);
244
245mat_trace(m) := block([n, acc : 0],
246  if not matrixp(m) then funmake('mat_trace, [m]) else (
247    n : matrix_size(m),
248    if first(n) # second(n) then error("The first argument to 'mat_trace' must be a square matrix"),
249    n : first(n),
250    for i from 1 thru n do (
251      acc : acc + if matrixp(m[i,i]) then mat_trace(m[i,i]) else m[i,i]),
252    acc));
253
254kronecker_product(a,b) := (
255  require_matrix(a,"first", "kronecker_product"),
256  require_matrix(b,"second", "kronecker_product"),
257  mat_unblocker(outermap('matrix_element_mult, a,b)));
258
259diag_matrix([d]) := block([f, n, prederror : false],
260
261  if every(lambda([s], matrixp(s) and not(blockmatrixp(s))), d) then (
262    f : lambda([i,j], if i = j then inpart(d,i)
263	       else zeromatrix(first(matrix_size(inpart(d,i))), second(matrix_size(inpart(d,j))))))
264
265  else if some(lambda([s], matrixp(s) or blockmatrixp(s)), d) then
266     error("All arguments to 'diag_matrix' must either be unblocked matrices or non-matrices")
267
268  else f : lambda([i,j], if i = j then inpart(d,i) else 0),
269
270  n : length(d),
271  genmatrix(f, n, n));
272
273hilbert_matrix(n) := block([row, mat : []],
274  mode_declare(n,integer),
275  require_posinteger(n,"first","hilbert_matrix"),
276  row : makelist(1/i,i,1,n),
277  for i : 1 thru n do (
278    mat : cons(row, mat),
279    row : endcons(1/(n + i), rest(row))),
280  funmake('matrix, reverse(mat)));
281
282hankel([q]) := block([col,row,m,n,partswitch : false],
283  if length(q) > 2 or length(q) < 1 then error("The function 'hankel' requires one or two arguments"),
284  col : inpart(q,1),
285  row : if length(q) = 2 then inpart(q,2) else map(lambda([x],0), col),
286
287  require_list(row,"first","hankel"),
288  require_list(col,"second","hankel"),
289  m : length(row),
290  n : length(col),
291  genmatrix(lambda([i,j],if i+j-1 <= n then inpart(col,i+j-1) else inpart(row,i+j-n)),n,m));
292
293toeplitz([q]) := block([col,row,m,n,partswitch : false],
294  if length(q) > 2 or length(q) < 1 then error("The function 'toeplitz' requires one or two arguments"),
295  col : inpart(q,1),
296  row : if length(q) = 2 then inpart(q,2) else map('conjugate, col),
297
298  require_list(row,"first","toeplitz"),
299  require_list(col,"second","toeplitz"),
300  m : length(row),
301  n : length(col),
302  genmatrix(lambda([i,j], if i -j >= 0 then inpart(col, i-j+1) else inpart(row,j-i+1)),n,m));
303
304polytocompanion(p,x) := block([n],
305  if not polynomialp(p,[x], lambda([e], freeof(x,e))) then
306     error("First argument to 'polytocompanion' must be a polynomial"),
307  p : expand(p),
308  n : hipow(p,x),
309  p : multthru(p / coeff(p,x,n)),
310  genmatrix(lambda([i,j], if i-j=1 then 1 else if j = n then -coeff(p,x,i-1) else 0),n));
311
312moore_penrose_pseudoinverse(m) := block([mm, z : ?gensym(), scalarmatrixp : false],
313  require_unblockedmatrix(m, "first", "moore_penrose_pseudoinverse"),
314  mm : m . ctranspose(m),
315  limit(mm : ctranspose(m) . (mm - z * identfor(mm))^^-1,z,0));
316