Rename specvar integer-info to *integer-info*
[maxima.git] / doc / info / linearalgebra.texi
bloba065786a337f48444fb97f73fb71a5e00394b987
1 @menu
2 * Introduction to linearalgebra::
3 * Functions and Variables for linearalgebra::
4 @end menu
6 @c -----------------------------------------------------------------------------
7 @node Introduction to linearalgebra, Functions and Variables for linearalgebra, Package linearalgebra, Package linearalgebra
8 @section Introduction to linearalgebra
9 @c -----------------------------------------------------------------------------
11 @code{linearalgebra} is a collection of functions for linear algebra.
13 Example:
15 @c ===beg===
16 @c M : matrix ([1, 2], [1, 2]);
17 @c nullspace (M);
18 @c columnspace (M);
19 @c ptriangularize (M - z*ident(2), z);
20 @c M : matrix ([1, 2, 3], [4, 5, 6], [7, 8, 9]) - z*ident(3);
21 @c MM : ptriangularize (M, z);
22 @c algebraic : true;
23 @c tellrat (MM [3, 3]);
24 @c MM : ratsimp (MM);
25 @c nullspace (MM);
26 @c M : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12], 
27 @c                    [13, 14, 15, 16]);
28 @c columnspace (M);
29 @c apply ('orthogonal_complement, args (nullspace (transpose (M))));
30 @c ===end===
31 @example
32 (%i1) M : matrix ([1, 2], [1, 2]);
33                             [ 1  2 ]
34 (%o1)                       [      ]
35                             [ 1  2 ]
36 (%i2) nullspace (M);
37                                [  1  ]
38                                [     ]
39 (%o2)                     span([   1 ])
40                                [ - - ]
41                                [   2 ]
42 (%i3) columnspace (M);
43                                 [ 1 ]
44 (%o3)                      span([   ])
45                                 [ 1 ]
46 (%i4) ptriangularize (M - z*ident(2), z);
47                          [ 1   2 - z   ]
48 (%o4)                    [             ]
49                          [           2 ]
50                          [ 0  3 z - z  ]
51 (%i5) M : matrix ([1, 2, 3], [4, 5, 6], [7, 8, 9]) - z*ident(3);
52                      [ 1 - z    2      3   ]
53                      [                     ]
54 (%o5)                [   4    5 - z    6   ]
55                      [                     ]
56                      [   7      8    9 - z ]
57 (%i6) MM : ptriangularize (M, z);
58               [ 4  5 - z            6            ]
59               [                                  ]
60               [                2                 ]
61               [     66        z    102 z   132   ]
62               [ 0   --      - -- + ----- + ---   ]
63 (%o6)         [     49        7     49     49    ]
64               [                                  ]
65               [               3        2         ]
66               [           49 z    245 z    147 z ]
67               [ 0    0    ----- - ------ - ----- ]
68               [            264      88      44   ]
69 (%i7) algebraic : true;
70 (%o7)                         true
71 (%i8) tellrat (MM [3, 3]);
72                          3       2
73 (%o8)                  [z  - 15 z  - 18 z]
74 (%i9) MM : ratsimp (MM);
75                [ 4  5 - z           6           ]
76                [                                ]
77                [                2               ]
78 (%o9)          [     66      7 z  - 102 z - 132 ]
79                [ 0   --    - ------------------ ]
80                [     49              49         ]
81                [                                ]
82                [ 0    0             0           ]
83 (%i10) nullspace (MM);
84                         [        1         ]
85                         [                  ]
86                         [   2              ]
87                         [  z  - 14 z - 16  ]
88                         [  --------------  ]
89 (%o10)             span([        8         ])
90                         [                  ]
91                         [    2             ]
92                         [   z  - 18 z - 12 ]
93                         [ - -------------- ]
94                         [         12       ]
95 (%i11) M : matrix ([1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12],
96                    [13, 14, 15, 16]);
97                        [ 1   2   3   4  ]
98                        [                ]
99                        [ 5   6   7   8  ]
100 (%o11)                 [                ]
101                        [ 9   10  11  12 ]
102                        [                ]
103                        [ 13  14  15  16 ]
104 (%i12) columnspace (M);
105                            [ 1  ]  [ 2  ]
106                            [    ]  [    ]
107                            [ 5  ]  [ 6  ]
108 (%o12)                span([    ], [    ])
109                            [ 9  ]  [ 10 ]
110                            [    ]  [    ]
111                            [ 13 ]  [ 14 ]
112 (%i13) apply ('orthogonal_complement, args (nullspace (transpose (M))));
113                            [ 0 ]  [  1  ]
114                            [   ]  [     ]
115                            [ 1 ]  [  0  ]
116 (%o13)                span([   ], [     ])
117                            [ 2 ]  [ - 1 ]
118                            [   ]  [     ]
119                            [ 3 ]  [ - 2 ]
120 @end example
122 @opencatbox{Categories:}
123 @category{Linear algebra}
124 @category{Share packages}
125 @category{Package linearalgebra}
126 @closecatbox
128 @c -----------------------------------------------------------------------------
129 @need 800
130 @node Functions and Variables for linearalgebra,  , Introduction to linearalgebra, Package linearalgebra
131 @section Functions and Variables for linearalgebra
132 @c -----------------------------------------------------------------------------
134 @c -----------------------------------------------------------------------------
135 @anchor{addmatrices}
136 @deffn {Function} addmatrices (@var{f}, @var{M_1}, @dots{}, @var{M_n})
138 @c REWORD -- THE RESULT IS NOT GENERALLY THE SUM OF M_1, ..., M_N
139 Using the function @var{f} as the addition function, return the sum of the
140 matrices @var{M_1}, @dots{}, @var{M_n}.  The function @var{f} must accept any
141 number of arguments (a Maxima nary function).
143 Examples:
145 @c ===beg===
146 @c m1 : matrix([1,2],[3,4])$
147 @c m2 : matrix([7,8],[9,10])$
148 @c addmatrices('max,m1,m2);
149 @c addmatrices('max,m1,m2,5*m1);
150 @c ===end===
151 @example
152 (%i1) m1 : matrix([1,2],[3,4])$
153 (%i2) m2 : matrix([7,8],[9,10])$
154 (%i3) addmatrices('max,m1,m2);
155 (%o3) matrix([7,8],[9,10])
156 (%i4) addmatrices('max,m1,m2,5*m1);
157 (%o4) matrix([7,10],[15,20])
158 @end example
160 @opencatbox{Categories:}
161 @category{Package linearalgebra}
162 @closecatbox
163 @end deffn
165 @c -----------------------------------------------------------------------------
166 @anchor{blockmatrixp}
167 @deffn {Function} blockmatrixp (@var{M})
169 Return true if and only if @var{M} is a matrix and every entry of 
170 @var{M} is a matrix.
172 @opencatbox{Categories:}
173 @category{Package linearalgebra}
174 @category{Predicate functions}
175 @closecatbox
176 @end deffn
178 @c -----------------------------------------------------------------------------
179 @anchor{columnop}
180 @deffn {Function} columnop (@var{M}, @var{i}, @var{j}, @var{theta})
182 If @var{M} is a matrix, return the matrix that results from doing the column
183 operation @code{C_i <- C_i - @var{theta} * C_j}. If @var{M} doesn't have a row
184 @var{i} or @var{j}, signal an error.
186 @opencatbox{Categories:}
187 @category{Package linearalgebra}
188 @closecatbox
189 @end deffn
191 @c -----------------------------------------------------------------------------
192 @anchor{columnswap}
193 @deffn {Function} columnswap (@var{M}, @var{i}, @var{j})
195 If @var{M} is a matrix, swap columns @var{i} and @var{j}.  If @var{M} doesn't
196 have a column @var{i} or @var{j}, signal an error.
198 @opencatbox{Categories:}
199 @category{Package linearalgebra}
200 @closecatbox
201 @end deffn
203 @c -----------------------------------------------------------------------------
204 @anchor{columnspace}
205 @deffn {Function} columnspace (@var{M})
207 If @var{M} is a matrix, return @code{span (v_1, ..., v_n)}, where the set
208 @code{@{v_1, ..., v_n@}} is a basis for the column space of @var{M}.  The span 
209 of the empty set is @code{@{0@}}.  Thus, when the column space has only
210 one member, return @code{span ()}.
212 @opencatbox{Categories:}
213 @category{Package linearalgebra}
214 @closecatbox
215 @end deffn
217 @c -----------------------------------------------------------------------------
218 @anchor{cholesky}
219 @deffn  {Function} cholesky @
220 @fname{cholesky} (@var{M}) @
221 @fname{cholesky} (@var{M}, @var{field})
223 Return the Cholesky factorization of the matrix selfadjoint (or hermitian)
224 matrix @var{M}.  The second argument defaults to 'generalring.' For a
225 description of the possible values for @var{field}, see @code{lu_factor}.
227 @opencatbox{Categories:}
228 @category{Matrix decompositions}
229 @category{Package linearalgebra}
230 @closecatbox
231 @end deffn
233 @c -----------------------------------------------------------------------------
234 @anchor{ctranspose}
235 @deffn {Function} ctranspose (@var{M})
237 Return the complex conjugate transpose of the matrix @var{M}.  The function
238 @code{ctranspose} uses @code{matrix_element_transpose} to transpose each matrix
239 element.
241 @opencatbox{Categories:}
242 @category{Package linearalgebra}
243 @closecatbox
244 @end deffn
246 @c -----------------------------------------------------------------------------
247 @anchor{diag_matrix}
248 @deffn {Function} diag_matrix (@var{d_1}, @var{d_2}, @dots{}, @var{d_n})
250 Return a diagonal matrix with diagonal entries @var{d_1}, @var{d_2}, @dots{},
251 @var{d_n}.  When the diagonal entries are matrices, the zero entries of the
252 returned matrix are zero matrices of the appropriate size; for example:
254 @c ===beg===
255 @c diag_matrix(diag_matrix(1,2),diag_matrix(3,4));
256 @c diag_matrix(p,q);
257 @c ===end===
258 @example
259 (%i1) diag_matrix(diag_matrix(1,2),diag_matrix(3,4));
261                             [ [ 1  0 ]  [ 0  0 ] ]
262                             [ [      ]  [      ] ]
263                             [ [ 0  2 ]  [ 0  0 ] ]
264 (%o1)                       [                    ]
265                             [ [ 0  0 ]  [ 3  0 ] ]
266                             [ [      ]  [      ] ]
267                             [ [ 0  0 ]  [ 0  4 ] ]
268 (%i2) diag_matrix(p,q);
270                                    [ p  0 ]
271 (%o2)                              [      ]
272                                    [ 0  q ]
273 @end example
275 @opencatbox{Categories:}
276 @category{Package linearalgebra}
277 @closecatbox
278 @end deffn
280 @c -----------------------------------------------------------------------------
281 @anchor{dotproduct}
282 @deffn {Function} dotproduct (@var{u}, @var{v})
284 Return the dotproduct of vectors @var{u} and @var{v}.  This is the same as
285 @code{conjugate (transpose (@var{u})) .  @var{v}}.  The arguments @var{u} and
286 @var{v} must be column vectors.
288 @opencatbox{Categories:}
289 @category{Package linearalgebra}
290 @closecatbox
291 @end deffn
293 @c -----------------------------------------------------------------------------
294 @anchor{eigens_by_jacobi}
295 @deffn  {Function} eigens_by_jacobi @
296 @fname{eigens_by_jacobi} (@var{A}) @
297 @fname{eigens_by_jacobi} (@var{A}, @var{field_type})
299 Computes the eigenvalues and eigenvectors of @var{A} by the method of Jacobi
300 rotations.  @var{A} must be a symmetric matrix (but it need not be positive
301 definite nor positive semidefinite).  @var{field_type} indicates the
302 computational field, either @code{floatfield} or @code{bigfloatfield}.
303 If @var{field_type} is not specified, it defaults to @code{floatfield}.
305 The elements of @var{A} must be numbers or expressions which evaluate to numbers
306 via @code{float} or @code{bfloat} (depending on @var{field_type}).
308 Examples:
310 @c ===beg===
311 @c S : matrix ([1/sqrt(2), 1/sqrt(2)], [- 1/sqrt(2), 1/sqrt(2)]);
312 @c L : matrix ([sqrt(3), 0], [0, sqrt(5)]);
313 @c M : S . L . transpose (S);
314 @c eigens_by_jacobi (M);
315 @c float ([[sqrt(3), sqrt(5)], S]);
316 @c eigens_by_jacobi (M, bigfloatfield);
317 @c ===end===
318 @example
319 (%i1) S: matrix([1/sqrt(2), 1/sqrt(2)],[-1/sqrt(2), 1/sqrt(2)]);
320                      [     1         1    ]
321                      [  -------   ------- ]
322                      [  sqrt(2)   sqrt(2) ]
323 (%o1)                [                    ]
324                      [      1        1    ]
325                      [ - -------  ------- ]
326                      [   sqrt(2)  sqrt(2) ]
327 (%i2) L : matrix ([sqrt(3), 0], [0, sqrt(5)]);
328                       [ sqrt(3)     0    ]
329 (%o2)                 [                  ]
330                       [    0     sqrt(5) ]
331 (%i3) M : S . L . transpose (S);
332             [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
333             [ ------- + -------  ------- - ------- ]
334             [    2         2        2         2    ]
335 (%o3)       [                                      ]
336             [ sqrt(5)   sqrt(3)  sqrt(5)   sqrt(3) ]
337             [ ------- - -------  ------- + ------- ]
338             [    2         2        2         2    ]
339 (%i4) eigens_by_jacobi (M);
340 The largest percent change was 0.1454972243679
341 The largest percent change was 0.0
342 number of sweeps: 2
343 number of rotations: 1
344 (%o4) [[1.732050807568877, 2.23606797749979], 
345                         [  0.70710678118655   0.70710678118655 ]
346                         [                                      ]]
347                         [ - 0.70710678118655  0.70710678118655 ]
348 (%i5) float ([[sqrt(3), sqrt(5)], S]);
349 (%o5) [[1.732050807568877, 2.23606797749979], 
350                         [  0.70710678118655   0.70710678118655 ]
351                         [                                      ]]
352                         [ - 0.70710678118655  0.70710678118655 ]
353 (%i6) eigens_by_jacobi (M, bigfloatfield);
354 The largest percent change was 1.454972243679028b-1
355 The largest percent change was 0.0b0
356 number of sweeps: 2
357 number of rotations: 1
358 (%o6) [[1.732050807568877b0, 2.23606797749979b0], 
359                 [  7.071067811865475b-1   7.071067811865475b-1 ]
360                 [                                              ]]
361                 [ - 7.071067811865475b-1  7.071067811865475b-1 ]
362 @end example
364 @opencatbox{Categories:}
365 @category{Matrix decompositions}
366 @category{Package linearalgebra}
367 @closecatbox
368 @end deffn
370 @c -----------------------------------------------------------------------------
371 @anchor{get_lu_factors}
372 @deffn {Function} get_lu_factors (@var{x}) 
374 When @code{@var{x} = lu_factor (@var{A})}, then @code{get_lu_factors} returns a
375 list of the form @code{[P, L, U]}, where @var{P} is a permutation matrix,
376 @var{L} is lower triangular with ones on the diagonal, and @var{U} is upper
377 triangular, and @code{@var{A} = @var{P} @var{L} @var{U}}.
379 @opencatbox{Categories:}
380 @category{Package linearalgebra}
381 @closecatbox
382 @end deffn
384 @c -----------------------------------------------------------------------------
385 @anchor{hankel}
386 @deffn  {Function} hankel @
387 @fname{hankel} (@var{col}) @
388 @fname{hankel} (@var{col}, @var{row})
390 Return a Hankel matrix @var{H}.  The first column of @var{H} is @var{col};
391 except for the first entry, the last row of @var{H} is @var{row}.  The
392 default for @var{row} is the zero vector with the same length as @var{col}.
394 @opencatbox{Categories:}
395 @category{Package linearalgebra}
396 @closecatbox
397 @end deffn
399 @c -----------------------------------------------------------------------------
400 @anchor{hessian}
401 @deffn {Function} hessian (@var{f}, @var{x})
403 Returns the Hessian matrix of @var{f} with respect to the list of variables
404 @var{x}.  The @code{(i, j)}-th element of the Hessian matrix is
405 @code{diff(@var{f}, @var{x}[i], 1, @var{x}[j], 1)}.
407 Examples:
409 @c ===beg===
410 @c hessian (x * sin (y), [x, y]);
411 @c depends (F, [a, b]);
412 @c hessian (F, [a, b]);
413 @c ===end===
414 @example
415 (%i1) hessian (x * sin (y), [x, y]);
416                      [   0       cos(y)   ]
417 (%o1)                [                    ]
418                      [ cos(y)  - x sin(y) ]
419 (%i2) depends (F, [a, b]);
420 (%o2)                       [F(a, b)]
421 (%i3) hessian (F, [a, b]);
422                         [   2      2   ]
423                         [  d F    d F  ]
424                         [  ---   ----- ]
425                         [    2   da db ]
426                         [  da          ]
427 (%o3)                   [              ]
428                         [   2      2   ]
429                         [  d F    d F  ]
430                         [ -----   ---  ]
431                         [ da db     2  ]
432                         [         db   ]
433 @end example
435 @opencatbox{Categories:}
436 @category{Differential calculus}
437 @category{Package linearalgebra}
438 @closecatbox
439 @end deffn
441 @c -----------------------------------------------------------------------------
442 @anchor{hilbert_matrix}
443 @deffn {Function} hilbert_matrix (@var{n})
445 Return the @var{n} by @var{n} Hilbert matrix.  When @var{n} isn't a positive
446 integer, signal an error.
448 @opencatbox{Categories:}
449 @category{Package linearalgebra}
450 @closecatbox
451 @end deffn
453 @c -----------------------------------------------------------------------------
454 @anchor{identfor}
455 @deffn  {Function} identfor @
456 @fname{identfor} (@var{M}) @
457 @fname{identfor} (@var{M}, @var{fld})
459 Return an identity matrix that has the same shape as the matrix
460 @var{M}.  The diagonal entries of the identity matrix are the 
461 multiplicative identity of the field @var{fld}; the default for
462 @var{fld} is @var{generalring}.
464 The first argument @var{M} should be a square matrix or a non-matrix.  When
465 @var{M} is a matrix, each entry of @var{M} can be a square matrix -- thus
466 @var{M} can be a blocked Maxima matrix.  The matrix can be blocked to any
467 (finite) depth.
469 See also @mref{zerofor}
471 @opencatbox{Categories:}
472 @category{Package linearalgebra}
473 @closecatbox
474 @end deffn
476 @c -----------------------------------------------------------------------------
477 @anchor{invert_by_lu}
478 @deffn {Function} invert_by_lu (@var{M}, @var{(rng generalring)})
480 Invert a matrix @var{M} by using the LU factorization.  The LU factorization
481 is done using the ring @var{rng}.
483 @opencatbox{Categories:}
484 @category{Package linearalgebra}
485 @closecatbox
486 @end deffn
488 @c -----------------------------------------------------------------------------
489 @anchor{jacobian}
490 @deffn {Function} jacobian (@var{f}, @var{x})
492 Returns the Jacobian matrix of the list of functions @var{f} with respect to
493 the list of variables @var{x}.  The @code{(i, j)}-th element of the Jacobian
494 matrix is @code{diff(@var{f}[i], @var{x}[j])}.
496 Examples:
498 @c ===beg===
499 @c jacobian ([sin (u - v), sin (u * v)], [u, v]);
500 @c depends ([F, G], [y, z]);
501 @c jacobian ([F, G], [y, z]);
502 @c ===end===
503 @example
504 (%i1) jacobian ([sin (u - v), sin (u * v)], [u, v]);
505                   [ cos(v - u)  - cos(v - u) ]
506 (%o1)             [                          ]
507                   [ v cos(u v)   u cos(u v)  ]
508 (%i2) depends ([F, G], [y, z]);
509 (%o2)                  [F(y, z), G(y, z)]
510 (%i3) jacobian ([F, G], [y, z]);
511                            [ dF  dF ]
512                            [ --  -- ]
513                            [ dy  dz ]
514 (%o3)                      [        ]
515                            [ dG  dG ]
516                            [ --  -- ]
517                            [ dy  dz ]
518 @end example
520 @opencatbox{Categories:}
521 @category{Differential calculus}
522 @category{Package linearalgebra}
523 @closecatbox
524 @end deffn
526 @c -----------------------------------------------------------------------------
527 @anchor{kronecker_product}
528 @deffn {Function} kronecker_product (@var{A}, @var{B})
530 Return the Kronecker product of the matrices @var{A} and @var{B}.
532 @opencatbox{Categories:}
533 @category{Package linearalgebra}
534 @closecatbox
535 @end deffn
537 @c -----------------------------------------------------------------------------
538 @anchor{locate_matrix_entry}
539 @deffn {Function} locate_matrix_entry (@var{M}, @var{r_1}, @var{c_1}, @var{r_2}, @var{c_2}, @var{f}, @var{rel})
541 The first argument must be a matrix; the arguments
542 @var{r_1} through @var{c_2} determine a sub-matrix of @var{M} that consists of
543 rows @var{r_1} through @var{r_2} and columns @var{c_1} through @var{c_2}.
545 Find an entry in the sub-matrix @var{M} that satisfies some property.
546 Three cases:
548 (1) @code{@var{rel} = 'bool} and @var{f} a predicate: 
550 Scan the sub-matrix from left to right then top to bottom,
551 and return the index of the first entry that satisfies the 
552 predicate @var{f}.  If no matrix entry satisfies @var{f}, return @code{false}.
554 (2) @code{@var{rel} = 'max} and @var{f} real-valued:
556 Scan the sub-matrix looking for an entry that maximizes @var{f}.
557 Return the index of a maximizing entry.
559 (3) @code{@var{rel} = 'min} and @var{f} real-valued:
561 Scan the sub-matrix looking for an entry that minimizes @var{f}.
562 Return the index of a minimizing entry.
564 @opencatbox{Categories:}
565 @category{Package linearalgebra}
566 @closecatbox
567 @end deffn
569 @c -----------------------------------------------------------------------------
570 @anchor{lu_backsub}
571 @deffn {Function} lu_backsub (@var{M}, @var{b})
573 When @code{@var{M} = lu_factor (@var{A}, @var{field})},
574 then @code{lu_backsub (@var{M}, @var{b})} solves the linear
575 system @code{@var{A} @var{x} = @var{b}}.
577 The @var{n} by @var{m} matrix @code{@var{b}}, with @var{n} the number of 
578 rows of the matrix @code{@var{A}}, contains one right hand side per column. If 
579 there is only one right hand side then @code{@var{b}} must be a @var{n} by 1
580 matrix.
582 Each column of the matrix @code{@var{x}=lu_backsub (@var{M}, @var{b})} is the 
583 solution corresponding to the respective column of @code{@var{b}}.
585 Examples:
587 @c ===beg===
588 @c A : matrix ([1 - z, 3], [3, 8 - z]);
589 @c M : lu_factor (A,generalring);
590 @c b : matrix([a],[c]);
591 @c x : lu_backsub(M,b);
592 @c ratsimp(A . x - b);
593 @c B : matrix([a,d],[c,f]);
594 @c x : lu_backsub(M,B);
595 @c ratsimp(A . x - B);
596 @c ===end===
597 @example
598 (%i1) A : matrix ([1 - z, 3], [3, 8 - z]);
599                                [ 1 - z    3   ]
600 (%o1)                          [              ]
601                                [   3    8 - z ]
602 (%i2) M : lu_factor (A,generalring);
603                [ 1 - z          3         ]
604                [                          ]
605 (%o2)         [[   3              9       ], [1, 2], generalring]
606                [ -----  (- z) - ----- + 8 ]
607                [ 1 - z          1 - z     ]
608 (%i3) b : matrix([a],[c]);
609                                      [ a ]
610 (%o3)                                [   ]
611                                      [ c ]
612 (%i4) x : lu_backsub(M,b);
613                            [               3 a     ]
614                            [       3 (c - -----)   ]
615                            [              1 - z    ]
616                            [ a - ----------------- ]
617                            [               9       ]
618                            [     (- z) - ----- + 8 ]
619                            [             1 - z     ]
620                            [ --------------------- ]
621 (%o4)                      [         1 - z         ]
622                            [                       ]
623                            [            3 a        ]
624                            [       c - -----       ]
625                            [           1 - z       ]
626                            [   -----------------   ]
627                            [             9         ]
628                            [   (- z) - ----- + 8   ]
629                            [           1 - z       ]
630 (%i5) ratsimp(A . x - b);
631                                      [ 0 ]
632 (%o5)                                [   ]
633                                      [ 0 ]
634 (%i6) B : matrix([a,d],[c,f]);
635                                    [ a  d ]
636 (%o6)                              [      ]
637                                    [ c  f ]
638 (%i7) x : lu_backsub(M,B);
639                [               3 a                    3 d     ]
640                [       3 (c - -----)          3 (f - -----)   ]
641                [              1 - z                  1 - z    ]
642                [ a - -----------------  d - ----------------- ]
643                [               9                      9       ]
644                [     (- z) - ----- + 8      (- z) - ----- + 8 ]
645                [             1 - z                  1 - z     ]
646                [ ---------------------  --------------------- ]
647 (%o7)          [         1 - z                  1 - z         ]
648                [                                              ]
649                [            3 a                    3 d        ]
650                [       c - -----              f - -----       ]
651                [           1 - z                  1 - z       ]
652                [   -----------------      -----------------   ]
653                [             9                      9         ]
654                [   (- z) - ----- + 8      (- z) - ----- + 8   ]
655                [           1 - z                  1 - z       ]
656 (%i8) ratsimp(A . x - B);
657                                    [ 0  0 ]
658 (%o8)                              [      ]
659                                    [ 0  0 ]
660 @end example
662 @opencatbox{Categories:}
663 @category{Package linearalgebra}
664 @closecatbox
665 @end deffn
667 @c -----------------------------------------------------------------------------
668 @anchor{lu_factor}
669 @deffn {Function} lu_factor (@var{M}, @var{field})
671 Return a list of the form @code{[@var{LU}, @var{perm}, @var{fld}]}, or
672 @code{[@var{LU}, @var{perm}, @var{fld}, @var{lower-cnd} @var{upper-cnd}]}, where
674 (1) The matrix @var{LU} contains the factorization of @var{M} in a packed form.
675     Packed form means three things: First, the rows of @var{LU} are permuted
676     according to the list @var{perm}.  If, for example, @var{perm} is the list
677     @code{[3,2,1]}, the actual first row of the @var{LU} factorization is the
678     third row of the matrix @var{LU}.  Second, the lower triangular factor of
679     m is the lower triangular part of @var{LU} with the diagonal entries
680     replaced by all ones.  Third, the upper triangular factor of @var{M} is the
681     upper triangular part of @var{LU}.
683 (2) When the field is either @code{floatfield} or @code{complexfield}, the
684     numbers @var{lower-cnd} and @var{upper-cnd} are lower and upper bounds for
685     the infinity norm condition number of @var{M}.  For all fields, the
686     condition number might not be estimated; for such fields, @code{lu_factor}
687     returns a two item list.  Both the lower and upper bounds can differ from
688     their true values by arbitrarily large factors.  (See also @mref{mat_cond}.)
689    
690   The argument @var{M} must be a square matrix.
692   The optional argument @var{fld} must be a symbol that determines a ring or
693   field.  The pre-defined fields and rings are:
695   (a) @code{generalring}      -- the ring of Maxima expressions,
697   (b) @code{floatfield}       -- the field of floating point numbers of the
698                                  type double,
700   (c) @code{complexfield}     -- the field of complex floating point numbers of
701                                  the type double,
703   (d) @code{crering}          -- the ring of Maxima CRE expressions,
705   (e) @code{rationalfield}    -- the field of rational numbers,
707   (f) @code{runningerror}     -- track the all floating point rounding errors,
709   (g) @code{noncommutingring} -- the ring of Maxima expressions where
710                                  multiplication is the non-commutative dot
711                                  operator.
713 When the field is @code{floatfield}, @code{complexfield}, or
714 @code{runningerror}, the algorithm uses partial pivoting; for all
715 other fields, rows are switched only when needed to avoid a zero
716 pivot.
718 Floating point addition arithmetic isn't associative, so the meaning
719 of 'field' differs from the mathematical definition.
721 A member of the field @code{runningerror} is a two member Maxima list
722 of the form @code{[x,n]},where @var{x} is a floating point number and
723 @code{n} is an integer.  The relative difference between the 'true'
724 value of @code{x} and @code{x} is approximately bounded by the machine
725 epsilon times @code{n}.  The running error bound drops some terms that
726 of the order the square of the machine epsilon.
728 There is no user-interface for defining a new field.  A user that is
729 familiar with Common Lisp should be able to define a new field.  To do
730 this, a user must define functions for the arithmetic operations and
731 functions for converting from the field representation to Maxima and
732 back.  Additionally, for ordered fields (where partial pivoting will be
733 used), a user must define functions for the magnitude and for
734 comparing field members.  After that all that remains is to define a
735 Common Lisp structure @code{mring}.  The file @code{mring} has many
736 examples.
738 To compute the factorization, the first task is to convert each matrix
739 entry to a member of the indicated field.  When conversion isn't
740 possible, the factorization halts with an error message.  Members of
741 the field needn't be Maxima expressions.  Members of the
742 @code{complexfield}, for example, are Common Lisp complex numbers.  Thus
743 after computing the factorization, the matrix entries must be
744 converted to Maxima expressions.
746 See also  @mref{get_lu_factors}.
748 Examples:
750 @c ===beg===
751 @c w[i,j] := random (1.0) + %i * random (1.0);
752 @c showtime : true$
753 @c M : genmatrix (w, 100, 100)$
754 @c lu_factor (M, complexfield)$
755 @c lu_factor (M, generalring)$
756 @c showtime : false$
757 @c M : matrix ([1 - z, 3], [3, 8 - z]);
758 @c lu_factor (M, generalring);
759 @c get_lu_factors (%);
760 @c %[1] . %[2] . %[3];
761 @c ===end===
762 @example
763 (%i1) w[i,j] := random (1.0) + %i * random (1.0);
764 (%o1)          w     := random(1.) + %i random(1.)
765                 i, j
766 (%i2) showtime : true$
767 Evaluation took 0.00 seconds (0.00 elapsed)
768 (%i3) M : genmatrix (w, 100, 100)$
769 Evaluation took 7.40 seconds (8.23 elapsed)
770 (%i4) lu_factor (M, complexfield)$
771 Evaluation took 28.71 seconds (35.00 elapsed)
772 (%i5) lu_factor (M, generalring)$
773 Evaluation took 109.24 seconds (152.10 elapsed)
774 (%i6) showtime : false$
776 (%i7) M : matrix ([1 - z, 3], [3, 8 - z]); 
777                         [ 1 - z    3   ]
778 (%o7)                   [              ]
779                         [   3    8 - z ]
780 (%i8) lu_factor (M, generalring);
781           [ 1 - z         3        ]
782           [                        ]
783 (%o8)    [[   3            9       ], [1, 2], generalring]
784           [ -----  - z - ----- + 8 ]
785           [ 1 - z        1 - z     ]
786 (%i9) get_lu_factors (%);
787                   [   1    0 ]  [ 1 - z         3        ]
788         [ 1  0 ]  [          ]  [                        ]
789 (%o9)  [[      ], [   3      ], [                9       ]]
790         [ 0  1 ]  [ -----  1 ]  [   0    - z - ----- + 8 ]
791                   [ 1 - z    ]  [              1 - z     ]
792 (%i10) %[1] . %[2] . %[3];
793                         [ 1 - z    3   ]
794 (%o10)                  [              ]
795                         [   3    8 - z ]
796 @end example
798 @opencatbox{Categories:}
799 @category{Matrix decompositions}
800 @category{Package linearalgebra}
801 @closecatbox
802 @end deffn
804 @c -----------------------------------------------------------------------------
805 @anchor{mat_cond}
806 @deffn  {Function} mat_cond @
807 @fname{mat_cond} (@var{M}, 1) @
808 @fname{mat_cond} (@var{M}, inf)
810 Return the @var{p}-norm matrix condition number of the matrix
811 @var{m}.  The allowed values for @var{p} are 1 and @var{inf}.  This
812 function uses the LU factorization to invert the matrix @var{m}.  Thus
813 the running time for @code{mat_cond} is proportional to the cube of
814 the matrix size; @code{lu_factor} determines lower and upper bounds
815 for the infinity norm condition number in time proportional to the
816 square of the matrix size.
818 @opencatbox{Categories:}
819 @category{Package linearalgebra}
820 @closecatbox
821 @end deffn
823 @c -----------------------------------------------------------------------------
824 @anchor{mat_norm}
825 @deffn  {Function} mat_norm @
826 @fname{mat_norm} (@var{M}, 1) @
827 @fname{mat_norm} (@var{M}, inf) @
828 @fname{mat_norm} (@var{M}, frobenius)
830 Return the matrix @var{p}-norm of the matrix @var{M}.  The allowed values for
831 @var{p} are 1, @code{inf}, and @code{frobenius} (the Frobenius matrix norm).
832 The matrix @var{M} should be an unblocked matrix.
834 @opencatbox{Categories:}
835 @category{Package linearalgebra}
836 @closecatbox
837 @end deffn
839 @c -----------------------------------------------------------------------------
840 @anchor{linearalgebra_matrixp}
841 @deffn  {Function} matrixp @
842 @fname{matrixp} (@var{e}, @var{p}) @
843 @fname{matrixp} (@var{e})
845 Given an optional argument @var{p}, return @code{true} if @var{e} is 
846 a matrix and @var{p} evaluates to @code{true} for every matrix element.
847 When @code{matrixp} is not given an optional argument, return @code{true} 
848 if @code{e} is a matrix.  In all other cases, return @code{false}.
850 See also @mref{blockmatrixp}
852 @opencatbox{Categories:}
853 @category{Package linearalgebra}
854 @category{Predicate functions}
855 @closecatbox
856 @end deffn
858 @c -----------------------------------------------------------------------------
859 @anchor{matrix_size}
860 @deffn {Function} matrix_size (@var{M})
862 Return a two member list that gives the number of rows and columns, respectively
863 of the matrix @var{M}.
865 @opencatbox{Categories:}
866 @category{Package linearalgebra}
867 @closecatbox
868 @end deffn
870 @c -----------------------------------------------------------------------------
871 @anchor{mat_fullunblocker}
872 @deffn {Function} mat_fullunblocker (@var{M})
874 If @var{M} is a block matrix, unblock the matrix to all levels.  If @var{M} is
875 a matrix, return @var{M}; otherwise, signal an error.
877 @opencatbox{Categories:}
878 @category{Package linearalgebra}
879 @closecatbox
880 @end deffn
882 @c -----------------------------------------------------------------------------
883 @anchor{mat_trace}
884 @deffn {Function} mat_trace (@var{M})
886 Return the trace of the matrix @var{M}.  If @var{M} isn't a matrix, return a
887 noun form.  When @var{M} is a block matrix, @code{mat_trace(M)} returns
888 the same value as does @code{mat_trace(mat_unblocker(m))}.
890 @opencatbox{Categories:}
891 @category{Package linearalgebra}
892 @closecatbox
893 @end deffn
895 @c -----------------------------------------------------------------------------
896 @anchor{mat_unblocker}
897 @deffn {Function} mat_unblocker (@var{M})
899 If @var{M} is a block matrix, unblock @var{M} one level.  If @var{M} is a
900 matrix, @code{mat_unblocker (M)} returns @var{M}; otherwise, signal an error.
902 Thus if each entry of @var{M} is matrix, @code{mat_unblocker (M)} returns an 
903 unblocked matrix, but if each entry of @var{M} is a block matrix,
904 @code{mat_unblocker (M)} returns a block matrix with one less level of blocking.
906 If you use block matrices, most likely you'll want to set
907 @code{matrix_element_mult} to @code{"."} and @code{matrix_element_transpose} to
908 @code{'transpose}.  See also @mref{mat_fullunblocker}.
910 Example:
912 @c ===beg===
913 @c A : matrix ([1, 2], [3, 4]);
914 @c B : matrix ([7, 8], [9, 10]);
915 @c matrix ([A, B]);
916 @c mat_unblocker (%);
917 @c ===end===
918 @example
919 (%i1) A : matrix ([1, 2], [3, 4]);
920                             [ 1  2 ]
921 (%o1)                       [      ]
922                             [ 3  4 ]
923 (%i2) B : matrix ([7, 8], [9, 10]);
924                             [ 7  8  ]
925 (%o2)                       [       ]
926                             [ 9  10 ]
927 (%i3) matrix ([A, B]);
928 @group
929                      [ [ 1  2 ]  [ 7  8  ] ]
930 (%o3)                [ [      ]  [       ] ]
931                      [ [ 3  4 ]  [ 9  10 ] ]
932 @end group
933 (%i4) mat_unblocker (%);
934                          [ 1  2  7  8  ]
935 (%o4)                    [             ]
936                          [ 3  4  9  10 ]
937 @end example
939 @opencatbox{Categories:}
940 @category{Package linearalgebra}
941 @closecatbox
942 @end deffn
944 @c -----------------------------------------------------------------------------
945 @anchor{nullspace}
946 @deffn {Function} nullspace (@var{M})
948 If @var{M} is a matrix, return @code{span (v_1, ..., v_n)}, where the set
949 @code{@{v_1, ..., v_n@}} is a basis for the nullspace of @var{M}.  The span of
950 the empty set is @code{@{0@}}.  Thus, when the nullspace has only one member,
951 return @code{span ()}.
953 @opencatbox{Categories:}
954 @category{Package linearalgebra}
955 @closecatbox
956 @end deffn
958 @c -----------------------------------------------------------------------------
959 @anchor{nullity}
960 @deffn {Function} nullity (@var{M})
962 If @var{M} is a matrix, return the dimension of the nullspace of @var{M}.
964 @opencatbox{Categories:}
965 @category{Package linearalgebra}
966 @closecatbox
967 @end deffn
969 @c -----------------------------------------------------------------------------
970 @anchor{orthogonal_complement}
971 @deffn {Function} orthogonal_complement (@var{v_1}, @dots{}, @var{v_n})
973 Return @code{span (u_1, ..., u_m)}, where the set @code{@{u_1, ..., u_m@}} is a 
974 basis for the orthogonal complement of the set @code{(v_1, ..., v_n)}.
976 Each vector @var{v_1} through @var{v_n} must be a column vector.
978 @opencatbox{Categories:}
979 @category{Package linearalgebra}
980 @closecatbox
981 @end deffn
983 @c -----------------------------------------------------------------------------
984 @anchor{polytocompanion}
985 @deffn {Function} polytocompanion (@var{p}, @var{x})
987 If @var{p} is a polynomial in @var{x}, return the companion matrix of @var{p}.
988 For a monic polynomial @var{p} of degree @var{n}, we have
989 @code{@var{p} = (-1)^@var{n} charpoly (polytocompanion (@var{p}, @var{x}))}.
991 When @var{p} isn't a polynomial in @var{x}, signal an error.
993 @opencatbox{Categories:}
994 @category{Package linearalgebra}
995 @closecatbox
996 @end deffn
998 @c -----------------------------------------------------------------------------
999 @anchor{ptringularize}
1000 @deffn {Function} ptriangularize (@var{M}, @var{v})
1002 If @var{M} is a matrix with each entry a polynomial in @var{v}, return 
1003 a matrix @var{M2} such that
1005 (1) @var{M2} is upper triangular,
1007 (2) @code{@var{M2} = @var{E_n} ... @var{E_1} @var{M}},
1008 where @var{E_1} through @var{E_n} are elementary matrices 
1009 whose entries are polynomials in @var{v},
1011 (3) @code{|det (@var{M})| = |det (@var{M2})|},
1013 Note: This function doesn't check that every entry is a polynomial in @var{v}.
1015 @opencatbox{Categories:}
1016 @category{Package linearalgebra}
1017 @closecatbox
1018 @end deffn
1020 @c -----------------------------------------------------------------------------
1021 @anchor{rowp}
1022 @deffn {Function} rowop (@var{M}, @var{i}, @var{j}, @var{theta})
1024 If @var{M} is a matrix, return the matrix that results from doing the
1025 row operation @code{R_i <- R_i - theta * R_j}.  If @var{M} doesn't have a row
1026 @var{i} or @var{j}, signal an error.
1028 @opencatbox{Categories:}
1029 @category{Package linearalgebra}
1030 @closecatbox
1032 @end deffn
1034 @anchor{linalg_rank}
1035 @deffn {Function} linalg_rank (@var{M})
1037 Return the rank of the matrix @var{M}. This function is equivalent to
1038 function @mref{rank}, but it uses a different algorithm: it finds the
1039 @mref{columnspace} of the matrix and counts its elements, since the rank
1040 of a matrix is the dimension of its column space. 
1042 @c ===beg===
1043 @c linalg_rank(matrix([1,2],[2,4]));
1044 @c linalg_rank(matrix([1,b],[c,d]));
1045 @c ===end===
1046 @example
1047 @group
1048 (%i1) linalg_rank(matrix([1,2],[2,4]));
1049 (%o1)                           1
1050 @end group
1051 @group
1052 (%i2) linalg_rank(matrix([1,b],[c,d]));
1053 (%o2)                           2
1054 @end group
1055 @end example
1057 @opencatbox{Categories:}
1058 @category{Package linearalgebra}
1059 @closecatbox
1060 @end deffn
1062 @c -----------------------------------------------------------------------------
1063 @anchor{rowswap}
1064 @deffn {Function} rowswap (@var{M}, @var{i}, @var{j})
1066 If @var{M} is a matrix, swap rows @var{i} and @var{j}.  If @var{M} doesn't
1067 have a row @var{i} or @var{j}, signal an error.
1069 @opencatbox{Categories:}
1070 @category{Package linearalgebra}
1071 @closecatbox
1072 @end deffn
1074 @c -----------------------------------------------------------------------------
1075 @anchor{toeplitz}
1076 @deffn  {Function} toeplitz @
1077 @fname{toeplitz} (@var{col}) @
1078 @fname{toeplitz} (@var{col}, @var{row})
1080 Return a Toeplitz matrix @var{T}.  The first first column of @var{T} is
1081 @var{col}; except for the first entry, the first row of @var{T} is @var{row}.
1082 The default for @var{row} is complex conjugate of @var{col}.  Example:
1084 @c ===beg===
1085 @c toeplitz([1,2,3],[x,y,z]);
1086 @c toeplitz([1,1+%i]);
1087 @c ===end===
1088 @example
1089 (%i1)  toeplitz([1,2,3],[x,y,z]);
1090 @group
1091                                   [ 1  y  z ]
1092                                   [         ]
1093 (%o1)                             [ 2  1  y ]
1094                                   [         ]
1095                                   [ 3  2  1 ]
1096 @end group
1097 (%i2)  toeplitz([1,1+%i]);
1099                               [   1     1 - %I ]
1100 (%o2)                         [                ]
1101                               [ %I + 1    1    ]
1102 @end example
1104 @opencatbox{Categories:}
1105 @category{Package linearalgebra}
1106 @closecatbox
1107 @end deffn
1109 @c -----------------------------------------------------------------------------
1110 @anchor{vandermonde_matrix}
1111 @deffn {Function} vandermonde_matrix ([@var{x_1}, ..., @var{x_n}])
1113 Return a @var{n} by @var{n} matrix whose @var{i}-th row is 
1114 @code{[1, @var{x_i}, @var{x_i}^2, ... @var{x_i}^(@var{n}-1)]}.
1116 @opencatbox{Categories:}
1117 @category{Package linearalgebra}
1118 @closecatbox
1119 @end deffn
1121 @c -----------------------------------------------------------------------------
1122 @anchor{zerofor}
1123 @deffn {Function} zerofor @
1124 @fname{zerofor} (@var{M}) @
1125 @fname{zerofor} (@var{M}, @var{fld})
1127 Return a zero  matrix that has the same shape as the matrix
1128 @var{M}.  Every entry of the zero matrix is the
1129 additive identity of the field @var{fld}; the default for
1130 @var{fld} is @var{generalring}.
1132 The first argument @var{M} should be a square matrix or a
1133 non-matrix.  When @var{M} is a matrix, each entry of @var{M} can be a
1134 square matrix -- thus @var{M} can be a blocked Maxima matrix.  The
1135 matrix can be blocked to any (finite) depth.
1137 See also @mref{identfor}
1139 @opencatbox{Categories:}
1140 @category{Package linearalgebra}
1141 @closecatbox
1142 @end deffn
1144 @c -----------------------------------------------------------------------------
1145 @anchor{zeromatrixp}
1146 @deffn {Function} zeromatrixp (@var{M})
1148 If @var{M} is not a block matrix, return @code{true} if
1149 @code{is (equal (@var{e}, 0))} is true for each element @var{e} of the matrix
1150 @var{M}.  If @var{M} is a block matrix, return @code{true} if @code{zeromatrixp}
1151 evaluates to @code{true} for each element of @var{e}.
1153 @opencatbox{Categories:}
1154 @category{Package linearalgebra}
1155 @category{Predicate functions}
1156 @closecatbox
1157 @end deffn