Introduction to Matrices and Linear Algebra
===========================================

Dot
---

The operator . represents noncommutative multiplication and scalar product. When the operands are 1-column or 1-row matrices a and b, the expression a.b is equivalent to sum (a[i]*b[i], i, 1, length(a)). If a and b are not complex, this is the scalar product, also called the inner product or dot product, of a and b. The scalar product is defined as conjugate(a).b when a and b are complex; innerproduct in the eigen package provides the complex scalar product.

When the operands are more general matrices, the product is the matrix product a and b. The number of rows of b must equal the number of columns of a, and the result has number of rows equal to the number of rows of a and number of columns equal to the number of columns of b.

To distinguish . as an arithmetic operator from the decimal point in a floating point number, it may be necessary to leave spaces on either side. For example, 5.e3 is 5000.0 but 5 . e3 is 5 times e3.

There are several flags which govern the simplification of expressions involving ., namely dot, dot0nscsimp, dot0simp, dot1simp, dotassoc, dotconstrules, dotdistrib, dotexptsimp, dotident, and dotscrules.

Vectors
-------

vect is a package of functions for vector analysis. load ("vect") loads this package, and demo ("vect") displays a demonstration.

The vector analysis package can combine and simplify symbolic expressions including dot products and cross products, together with the gradient, divergence, curl, and Laplacian operators.  The distribution of these operators over sums or products is governed by several flags, as are various other expansions, including expansion into components in any specific orthogonal coordinate systems. There are also functions for deriving the scalar or vector potential of a field.

The vect package contains these functions: vectorsimp, scalefactors, express, potential, and vectorpotential.

Warning: the vect package declares the dot operator . to be a commutative operator.

eigen
-----

The package eigen contains several functions devoted to the symbolic computation of eigenvalues and eigenvectors. Maxima loads the package automatically if one of the functions eigenvalues or eigenvectors is invoked. The package may be loaded explicitly as load ("eigen").

demo ("eigen") displays a demonstration of the capabilities of this package. batch ("eigen") executes the same demonstration, but without the user prompt between successive computations.

The functions in the eigen package are innerproduct, unitvector, columnvector, gramschmidt, eigenvalues, eigenvectors, uniteigenvectors, and similaritytransform.

Functions and Variables for Matrices and Linear Algebra
=======================================================

 -- Function: addcol (M, LIST_1, ..., LIST_N)
     Appends the column(s) given by the one or more lists (or matrices) onto the matrix M.

 -- Function: addrow (M, LIST_1, ..., LIST_N)
     Appends the row(s) given by the one or more lists (or matrices) onto the matrix M.

 -- Function: adjoint (M)
     Returns the adjoint of the matrix M. The adjoint matrix is the transpose of the matrix of cofactors of M.

 -- Function: augcoefmatrix ([EQN_1, ..., EQN_M], [X_1, ..., X_N])
     Returns the augmented coefficient matrix for the variables X_1, ..., X_N of the system of linear equations EQN_1, ..., EQN_M.  This is the coefficient matrix with a column adjoined for the constant terms in each equation (i.e., those terms not dependent upon X_1, ..., X_N).

          (%i1) m: [2*x - (a - 1)*y = 5*b, c + b*y + a*x = 0]$
          (%i2) augcoefmatrix (m, [x, y]);
                                 [ 2  1 - a  - 5 b ]
          (%o2)                  [                 ]
                                 [ a    b      c   ]

 -- Function: charpoly (M, X)
     Returns the characteristic polynomial for the matrix M with respect to variable X.  That is, determinant (M - diagmatrix (length (M), X)).

          (%i1) a: matrix ([3, 1], [2, 4]);
                                      [ 3  1 ]
          (%o1)                       [      ]
                                      [ 2  4 ]
          (%i2) expand (charpoly (a, lambda));
                                     2
          (%o2)                lambda  - 7 lambda + 10
          (%i3) (programmode: true, solve (%));
          (%o3)               [lambda = 5, lambda = 2]
          (%i4) matrix ([x1], [x2]);
                                       [ x1 ]
          (%o4)                        [    ]
                                       [ x2 ]
          (%i5) ev (a . % - lambda*%, %th(2)[1]);
                                    [ x2 - 2 x1 ]
          (%o5)                     [           ]
                                    [ 2 x1 - x2 ]
          (%i6) %[1, 1] = 0;
          (%o6)                     x2 - 2 x1 = 0
          (%i7) x2^2 + x1^2 = 1;
                                      2     2
          (%o7)                     x2  + x1  = 1
          (%i8) solve ([%th(2), %], [x1, x2]);
                            1               2
          (%o8) [[x1 = - -------, x2 = - -------],
                         sqrt(5)         sqrt(5)

                                                       1             2
                                              [x1 = -------, x2 = -------]]
                                                    sqrt(5)       sqrt(5)

 -- Function: coefmatrix ([EQN_1, ..., EQN_M], [X_1, ..., X_N])
     Returns the coefficient matrix for the variables X_1, ..., X_N of the system of linear equations EQN_1, ..., EQN_M.

          (%i1) coefmatrix([2*x-(a-1)*y+5*b = 0, b*y+a*x = 3], [x,y]);
                                           [ 2  1 - a ]
          (%o1)                            [          ]
                                           [ a    b   ]

 -- Function: col (M, I)
     Returns the I'th column of the matrix M. The return value is a matrix.

 -- Function: columnvector (L)
 -- Function: covect (L)
     Returns a matrix of one column and length (L) rows, containing the elements of the list L.

     covect is a synonym for columnvector.

     load ("eigen") loads this function.

     This is useful if you want to use parts of the outputs of the functions in this package in matrix calculations.

     Example:

          (%i1) load ("eigen")$
          Warning - you are redefining the Macsyma function eigenvalues
          Warning - you are redefining the Macsyma function eigenvectors
          (%i2) columnvector ([aa, bb, cc, dd]);
                                       [ aa ]
                                       [    ]
                                       [ bb ]
          (%o2)                        [    ]
                                       [ cc ]
                                       [    ]
                                       [ dd ]

 -- Function: conjugate (X)
     Returns the complex conjugate of X.

          (%i1) declare ([aa, bb], real, cc, complex, ii, imaginary);

          (%o1)                         done
          (%i2) conjugate (aa + bb*%i);

          (%o2)                      aa - %i bb
          (%i3) conjugate (cc);

          (%o3)                     conjugate(cc)
          (%i4) conjugate (ii);

          (%o4)                         - ii
          (%i5) conjugate (xx + yy);

          (%o5)             conjugate(yy) + conjugate(xx)

 -- Function: copymatrix (M)
     Returns a copy of the matrix M.  This is the only way to make a copy aside from copying M element by element.

     Note that an assignment of one matrix to another, as in m2: m1, does not copy m1. An assignment m2 [i,j]: x or setelmx (x, i, j, m2 also modifies m1 [i,j]. Creating a copy with copymatrix and then using assignment creates a separate, modified copy.

 -- Function: determinant (M)
     Computes the determinant of M by a method similar to Gaussian elimination.

     The form of the result depends upon the setting of the switch ratmx.

     There is a special routine for computing sparse determinants which is called when the switches ratmx and sparse are both true.

 -- Option variable: detout
     Default value: false

     When detout is true, the determinant of a matrix whose inverse is computed is factored out of the inverse.

     For this switch to have an effect doallmxops and doscmxops should be false (see their descriptions).  Alternatively this switch can be given to ev which causes the other two to be set correctly.

     Example:

          (%i1) m: matrix ([a, b], [c, d]);
                                      [ a  b ]
          (%o1)                       [      ]
                                      [ c  d ]
          (%i2) detout: true$
          (%i3) doallmxops: false$
          (%i4) doscmxops: false$
          (%i5) invert (m);
                                    [  d   - b ]
                                    [          ]
                                    [ - c   a  ]
          (%o5)                     ------------
                                     a d - b c

 -- Function: diagmatrix (N, X)
     Returns a diagonal matrix of size N by N with the diagonal elements all equal to X. diagmatrix (N, 1) returns an identity matrix (same as ident (N)).

     N must evaluate to an integer, otherwise diagmatrix complains with an error message.

     X can be any kind of expression, including another matrix. If X is a matrix, it is not copied; all diagonal elements refer to the same instance, X.

 -- Option variable: doallmxops
     Default value: true

     When doallmxops is true, all operations relating to matrices are carried out. When it is false then the setting of the individual dot switches govern which operations are performed.

 -- Option variable: domxexpt
     Default value: true

     When domxexpt is true, a matrix exponential, exp (M) where M is a matrix, is interpreted as a matrix with element [i,j equal to exp (m[i,j]). Otherwise exp (M) evaluates to exp (EV(M).

     domxexpt affects all expressions of the form BASE^POWER where BASE is an expression assumed scalar or constant, and POWER is a list or matrix.

     Example:

          (%i1) m: matrix ([1, %i], [a+b, %pi]);
                                   [   1    %i  ]
          (%o1)                    [            ]
                                   [ b + a  %pi ]
          (%i2) domxexpt: false$
          (%i3) (1 - c)^m;
                                       [   1    %i  ]
                                       [            ]
                                       [ b + a  %pi ]
          (%o3)                 (1 - c)
          (%i4) domxexpt: true$
          (%i5) (1 - c)^m;
                            [                      %i  ]
                            [    1 - c      (1 - c)    ]
          (%o5)             [                          ]
                            [        b + a         %pi ]
                            [ (1 - c)       (1 - c)    ]

 -- Option variable: domxmxops
     Default value: true

     When domxmxops is true, all matrix-matrix or matrix-list operations are carried out (but not scalar-matrix operations); if this switch is false such operations are not carried out.

 -- Option variable: domxnctimes
     Default value: false

     When domxnctimes is true, non-commutative products of matrices are carried out.

 -- Option variable: dontfactor
     Default value: []

     dontfactor may be set to a list of variables with respect to which factoring is not to occur.  (The list is initially empty.) Factoring also will not take place with respect to any variables which are less important, according the variable ordering assumed for canonical rational expression (CRE) form, than those on the dontfactor list.

 -- Option variable: doscmxops
     Default value: false

     When doscmxops is true, scalar-matrix operations are carried out.

 -- Option variable: doscmxplus
     Default value: false

     When doscmxplus is true, scalar-matrix operations yield a matrix result.  This switch is not subsumed under doallmxops.

 -- Option variable: dot0nscsimp
     Default value: true

     When dot0nscsimp is true, a non-commutative product of zero and a nonscalar term is simplified to a commutative product.

 -- Option variable: dot0simp
     Default value: true

     When dot0simp is true, a non-commutative product of zero and a scalar term is simplified to a commutative product.

 -- Option variable: dot1simp
     Default value: true

     When dot1simp is true, a non-commutative product of one and another term is simplified to a commutative product.

 -- Option variable: dotassoc
     Default value: true

     When dotassoc is true, an expression (A.B).C simplifies to A.(B.C).

 -- Option variable: dotconstrules
     Default value: true

     When dotconstrules is true, a non-commutative product of a constant and another term is simplified to a commutative product. Turning on this flag effectively turns on dot0simp, dot0nscsimp, and dot1simp as well.

 -- Option variable: dotdistrib
     Default value: false

     When dotdistrib is true, an expression A.(B + C) simplifies to A.B + A.C.

 -- Option variable: dotexptsimp
     Default value: true

     When dotexptsimp is true, an expression A.A simplifies to A^^2.

 -- Option variable: dotident
     Default value: 1

     dotident is the value returned by X^^0.

 -- Option variable: dotscrules
     Default value: false

     When dotscrules is true, an expression A.SC or SC.A simplifies to SC*A and A.(SC*B) simplifies to SC*(A.B).

 -- Function: echelon (M)
     Returns the echelon form of the matrix M, as produced by Gaussian elimination. The echelon form is computed from M by elementary row operations such that the first non-zero element in each row in the resulting matrix is one and the column elements under the first one in each row are all zero.

     triangularize also carries out Gaussian elimination, but it does not normalize the leading non-zero element in each row.

     lu_factor and cholesky are other functions which yield triangularized matrices.

          (%i1) M: matrix ([3, 7, aa, bb], [-1, 8, 5, 2], [9, 2, 11, 4]);
                                 [  3   7  aa  bb ]
                                 [                ]
          (%o1)                  [ - 1  8  5   2  ]
                                 [                ]
                                 [  9   2  11  4  ]
          (%i2) echelon (M);
                            [ 1  - 8  - 5      - 2     ]
                            [                          ]
                            [         28       11      ]
                            [ 0   1   --       --      ]
          (%o2)             [         37       37      ]
                            [                          ]
                            [              37 bb - 119 ]
                            [ 0   0    1   ----------- ]
                            [              37 aa - 313 ]

 -- Function: eigenvalues (M)
 -- Function: eivals (M)
     Returns a list of two lists containing the eigenvalues of the matrix M. The first sublist of the return value is the list of eigenvalues of the matrix, and the second sublist is the list of the multiplicities of the eigenvalues in the corresponding order.

     eivals is a synonym for eigenvalues.

     eigenvalues calls the function solve to find the roots of the characteristic polynomial of the matrix. Sometimes solve may not be able to find the roots of the polynomial; in that case some other functions in this package (except innerproduct, unitvector, columnvector and gramschmidt) will not work.

     In some cases the eigenvalues found by solve may be complicated expressions. (This may happen when solve returns a not-so-obviously real expression for an eigenvalue which is known to be real.) It may be possible to simplify the eigenvalues using some other functions.

     The package eigen.mac is loaded automatically when eigenvalues or eigenvectors is referenced. If eigen.mac is not already loaded, load ("eigen") loads it. After loading, all functions and variables in the package are available.

 -- Function: eigenvectors (M)
 -- Function: eivects (M)
     takes a matrix M as its argument and returns a list of lists the first sublist of which is the output of eigenvalues and the other sublists of which are the eigenvectors of the matrix corresponding to those eigenvalues respectively. The calculated eigenvectors and the unit eigenvectors of the matrix are the right eigenvectors and the right unit eigenvectors respectively.

     eivects is a synonym for eigenvectors.

     The package eigen.mac is loaded automatically when eigenvalues or eigenvectors is referenced. If eigen.mac is not already loaded, load ("eigen") loads it. After loading, all functions and variables in the package are available.

     The flags that affect this function are:

     nondiagonalizable is set to true or false depending on whether the matrix is nondiagonalizable or diagonalizable after eigenvectors returns.

     hermitianmatrix when true, causes the degenerate eigenvectors of the Hermitian matrix to be orthogonalized using the Gram-Schmidt algorithm.

     knowneigvals when true causes the eigen package to assume the eigenvalues of the matrix are known to the user and stored under the global name listeigvals.  listeigvals should be set to a list similar to the output eigenvalues.

     The function algsys is used here to solve for the eigenvectors. Sometimes if the eigenvalues are messy, algsys may not be able to find a solution. In some cases, it may be possible to simplify the eigenvalues by first finding them using eigenvalues command and then using other functions to reduce them to something simpler. Following simplification, eigenvectors can be called again with the knowneigvals flag set to true.

 -- Function: ematrix (M, N, X, I, J)
     Returns an M by N matrix, all elements of which are zero except for the [I, J] element which is X.

 -- Function: entermatrix (M, N)
     Returns an M by N matrix, reading the elements interactively.

     If N is equal to M, Maxima prompts for the type of the matrix (diagonal, symmetric, antisymmetric, or general) and for each element. Each response is terminated by a semicolon ; or dollar sign $.

     If N is not equal to M, Maxima prompts for each element.

     The elements may be any expressions, which are evaluated. entermatrix evaluates its arguments.

          (%i1) n: 3$
          (%i2) m: entermatrix (n, n)$

          Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric
          4. General
          Answer 1, 2, 3 or 4 :
          1$
          Row 1 Column 1:
          (a+b)^n$
          Row 2 Column 2:
          (a+b)^(n+1)$
          Row 3 Column 3:
          (a+b)^(n+2)$

          Matrix entered.
          (%i3) m;
                          [        3                     ]
                          [ (b + a)      0         0     ]
                          [                              ]
          (%o3)           [                  4           ]
                          [    0      (b + a)      0     ]
                          [                              ]
                          [                            5 ]
                          [    0         0      (b + a)  ]

 -- Function: genmatrix (A, I_2, J_2, I_1, J_1)
 -- Function: genmatrix (A, I_2, J_2, I_1)
 -- Function: genmatrix (A, I_2, J_2)
     Returns a matrix generated from A, taking element A[I_1,J_1] as the upper-left element and A[I_2,J_2] as the lower-right element of the matrix. Here A is a declared array (created by array but not by make_array) or an undeclared array, or an array function, or a lambda expression of two arguments. (An array function is created like other functions with := or define, but arguments are enclosed in square brackets instead of parentheses.)

     If J_1 is omitted, it is assumed equal to I_1. If both J_1 and I_1 are omitted, both are assumed equal to 1.

     If a selected element i,j of the array is undefined, the matrix will contain a symbolic element A[i,j].

     Examples:

          (%i1) h [i, j] := 1 / (i + j - 1);
                                              1
          (%o1)                  h     := ---------
                                  i, j    i + j - 1
          (%i2) genmatrix (h, 3, 3);
                                     [    1  1 ]
                                     [ 1  -  - ]
                                     [    2  3 ]
                                     [         ]
                                     [ 1  1  1 ]
          (%o2)                      [ -  -  - ]
                                     [ 2  3  4 ]
                                     [         ]
                                     [ 1  1  1 ]
                                     [ -  -  - ]
                                     [ 3  4  5 ]
          (%i3) array (a, fixnum, 2, 2);
          (%o3)                           a
          (%i4) a [1, 1] : %e;
          (%o4)                          %e
          (%i5) a [2, 2] : %pi;
          (%o5)                          %pi
          (%i6) genmatrix (a, 2, 2);
                                     [ %e   0  ]
          (%o6)                      [         ]
                                     [ 0   %pi ]
          (%i7) genmatrix (lambda ([i, j], j - i), 3, 3);
                                   [  0    1   2 ]
                                   [             ]
          (%o7)                    [ - 1   0   1 ]
                                   [             ]
                                   [ - 2  - 1  0 ]
          (%i8) genmatrix (B, 2, 2);
                                  [ B      B     ]
                                  [  1, 1   1, 2 ]
          (%o8)                   [              ]
                                  [ B      B     ]
                                  [  2, 1   2, 2 ]

 -- Function: gramschmidt (X)
 -- Function: gschmit (X)
     Carries out the Gram-Schmidt orthogonalization algorithm on X, which is either a matrix or a list of lists. X is not modified by gramschmidt.

     If X is a matrix, the algorithm is applied to the rows of X. If X is a list of lists, the algorithm is applied to the sublists, which must have equal numbers of elements. In either case, the return value is a list of lists, the sublists of which are orthogonal and span the same space as X. If the dimension of the span of X is less than the number of rows or sublists, some sublists of the return value are zero.

     factor is called at each stage of the algorithm to simplify intermediate results. As a consequence, the return value may contain factored integers.

     gschmit (note spelling) is a synonym for gramschmidt.

     load ("eigen") loads this function.

     Example:

          (%i1) load ("eigen")$
          Warning - you are redefining the Macsyma function eigenvalues
          Warning - you are redefining the Macsyma function eigenvectors
          (%i2) x: matrix ([1, 2, 3], [9, 18, 30], [12, 48, 60]);
                                   [ 1   2   3  ]
                                   [            ]
          (%o2)                    [ 9   18  30 ]
                                   [            ]
                                   [ 12  48  60 ]
          (%i3) y: gramschmidt (x);
                                 2      2            4     3
                                3      3   3 5      2  3  2  3
          (%o3)  [[1, 2, 3], [- ---, - --, ---], [- ----, ----, 0]]
                                2 7    7   2 7       5     5
          (%i4) i: innerproduct$
          (%i5) [i (y[1], y[2]), i (y[2], y[3]), i (y[3], y[1])];
          (%o5)                       [0, 0, 0]

 -- Function: ident (N)
     Returns an N by N identity matrix.

 -- Function: innerproduct (X, Y)
 -- Function: inprod (X, Y)
     Returns the inner product (also called the scalar product or dot product) of X and Y, which are lists of equal length, or both 1-column or 1-row matrices of equal length. The return value is conjugate (x) . y, where . is the noncommutative multiplication operator.

     load ("eigen") loads this function.

     inprod is a synonym for innerproduct.

 -- Function: invert (M)
     Returns the inverse of the matrix M. The inverse is computed by the adjoint method.

     This allows a user to compute the inverse of a matrix with bfloat entries or polynomials with floating pt. coefficients without converting to cre-form.

     Cofactors are computed by the determinant function, so if ratmx is false the inverse is computed without changing the representation of the elements.

     The current implementation is inefficient for matrices of high order.

     When detout is true, the determinant is factored out of the inverse.

     The elements of the inverse are not automatically expanded. If M has polynomial elements, better appearing output can be generated by expand (invert (m)), detout. If it is desirable to then divide through by the determinant this can be accomplished by xthru (%) or alternatively from scratch by

          expand (adjoint (m)) / expand (determinant (m))
          invert (m) := adjoint (m) / determinant (m)

     See ^^ (noncommutative exponent) for another method of inverting a matrix.

 -- Option variable: lmxchar
     Default value: [

     lmxchar is the character displayed as the left delimiter of a matrix. See also rmxchar.

     Example:

          (%i1) lmxchar: "|"$
          (%i2) matrix ([a, b, c], [d, e, f], [g, h, i]);
                                     | a  b  c ]
                                     |         ]
          (%o2)                      | d  e  f ]
                                     |         ]
                                     | g  h  i ]

 -- Function: matrix (ROW_1, ..., ROW_N)
     Returns a rectangular matrix which has the rows ROW_1, ..., ROW_N. Each row is a list of expressions. All rows must be the same length.

     The operations + (addition), - (subtraction), * (multiplication), and / (division), are carried out element by element when the operands are two matrices, a scalar and a matrix, or a matrix and a scalar. The operation ^ (exponentiation, equivalently **) is carried out element by element if the operands are a scalar and a matrix or a matrix and a scalar, but not if the operands are two matrices. All operations are normally carried out in full, including . (noncommutative multiplication).

     Matrix multiplication is represented by the noncommutative multiplication operator .. The corresponding noncommutative exponentiation operator is ^^. For a matrix A, A.A = A^^2 and A^^-1 is the inverse of A, if it exists.

     There are switches for controlling simplification of expressions involving dot and matrix-list operations. These are doallmxops, domxexpt domxmxops, doscmxops, and doscmxplus.

     There are additional options which are related to matrices. These are: lmxchar, rmxchar, ratmx, listarith, detout, scalarmatrix, and sparse.

     There are a number of functions which take matrices as arguments or yield matrices as return values. See eigenvalues, eigenvectors, determinant, charpoly, genmatrix, addcol, addrow, copymatrix, transpose, echelon, and rank.

     Examples:

       * Construction of matrices from lists.

          (%i1) x: matrix ([17, 3], [-8, 11]);
                                     [ 17   3  ]
          (%o1)                      [         ]
                                     [ - 8  11 ]
          (%i2) y: matrix ([%pi, %e], [a, b]);
                                     [ %pi  %e ]
          (%o2)                      [         ]
                                     [  a   b  ]

       * Addition, element by element.

          (%i3) x + y;
                                [ %pi + 17  %e + 3 ]
          (%o3)                 [                  ]
                                [  a - 8    b + 11 ]

       * Subtraction, element by element.

          (%i4) x - y;
                                [ 17 - %pi  3 - %e ]
          (%o4)                 [                  ]
                                [ - a - 8   11 - b ]

       * Multiplication, element by element.

          (%i5) x * y;
                                  [ 17 %pi  3 %e ]
          (%o5)                   [              ]
                                  [ - 8 a   11 b ]

       * Division, element by element.

          (%i6) x / y;
                                  [ 17       - 1 ]
                                  [ ---  3 %e    ]
                                  [ %pi          ]
          (%o6)                   [              ]
                                  [   8    11    ]
                                  [ - -    --    ]
                                  [   a    b     ]

       * Matrix to a scalar exponent, element by element.

          (%i7) x ^ 3;
                                   [ 4913    27  ]
          (%o7)                    [             ]
                                   [ - 512  1331 ]

       * Scalar base to a matrix exponent, element by element.

          (%i8) exp(y);
                                   [   %pi    %e ]
                                   [ %e     %e   ]
          (%o8)                    [             ]
                                   [    a     b  ]
                                   [  %e    %e   ]

       * Matrix base to a matrix exponent. This is not carried out element by element.

          (%i9) x ^ y;
                                          [ %pi  %e ]
                                          [         ]
                                          [  a   b  ]
                               [ 17   3  ]
          (%o9)                [         ]
                               [ - 8  11 ]

       * Noncommutative matrix multiplication.

          (%i10) x . y;
                            [ 3 a + 17 %pi  3 b + 17 %e ]
          (%o10)            [                           ]
                            [ 11 a - 8 %pi  11 b - 8 %e ]
          (%i11) y . x;
                          [ 17 %pi - 8 %e  3 %pi + 11 %e ]
          (%o11)          [                              ]
                          [  17 a - 8 b     11 b + 3 a   ]

       * Noncommutative matrix exponentiation. A scalar base B to a matrix power M is carried out element by element and so b^^m is the same as b^m.

          (%i12) x ^^ 3;
                                  [  3833   1719 ]
          (%o12)                  [              ]
                                  [ - 4584  395  ]
          (%i13) %e ^^ y;
                                   [   %pi    %e ]
                                   [ %e     %e   ]
          (%o13)                   [             ]
                                   [    a     b  ]
                                   [  %e    %e   ]

       * A matrix raised to a -1 exponent with noncommutative exponentiation is the matrix inverse, if it exists.

          (%i14) x ^^ -1;
                                   [ 11      3  ]
                                   [ ---  - --- ]
                                   [ 211    211 ]
          (%o14)                   [            ]
                                   [  8    17   ]
                                   [ ---   ---  ]
                                   [ 211   211  ]
          (%i15) x . (x ^^ -1);
                                      [ 1  0 ]
          (%o15)                      [      ]
                                      [ 0  1 ]

 -- Function: matrixmap (F, M)
     Returns a matrix with element i,j equal to F(M[i,j]).

     See also map, fullmap, fullmapl, and apply.

 -- Function: matrixp (EXPR)
     Returns true if EXPR is a matrix, otherwise false.

 -- Option variable: matrix_element_add
     Default value: +

     matrix_element_add is the operation invoked in place of addition in a matrix multiplication. matrix_element_add can be assigned any n-ary operator (that is, a function which handles any number of arguments). The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

     See also matrix_element_mult and matrix_element_transpose.

     Example:

          (%i1) matrix_element_add: "*"$
          (%i2) matrix_element_mult: "^"$
          (%i3) aa: matrix ([a, b, c], [d, e, f]);
                                     [ a  b  c ]
          (%o3)                      [         ]
                                     [ d  e  f ]
          (%i4) bb: matrix ([u, v, w], [x, y, z]);
                                     [ u  v  w ]
          (%o4)                      [         ]
                                     [ x  y  z ]
          (%i5) aa . transpose (bb);
                               [  u  v  w   x  y  z ]
                               [ a  b  c   a  b  c  ]
          (%o5)                [                    ]
                               [  u  v  w   x  y  z ]
                               [ d  e  f   d  e  f  ]

 -- Option variable: matrix_element_mult
     Default value: *

     matrix_element_mult is the operation invoked in place of multiplication in a matrix multiplication. matrix_element_mult can be assigned any binary operator. The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

     The dot operator . is a useful choice in some contexts.

     See also matrix_element_add and matrix_element_transpose.

     Example:

          (%i1) matrix_element_add: lambda ([[x]], sqrt (apply ("+", x)))$
          (%i2) matrix_element_mult: lambda ([x, y], (x - y)^2)$
          (%i3) [a, b, c] . [x, y, z];
                                    2          2          2
          (%o3)         sqrt((c - z)  + (b - y)  + (a - x) )
          (%i4) aa: matrix ([a, b, c], [d, e, f]);
                                     [ a  b  c ]
          (%o4)                      [         ]
                                     [ d  e  f ]
          (%i5) bb: matrix ([u, v, w], [x, y, z]);
                                     [ u  v  w ]
          (%o5)                      [         ]
                                     [ x  y  z ]
          (%i6) aa . transpose (bb);
                         [             2          2          2  ]
                         [ sqrt((c - w)  + (b - v)  + (a - u) ) ]
          (%o6)  Col 1 = [                                      ]
                         [             2          2          2  ]
                         [ sqrt((f - w)  + (e - v)  + (d - u) ) ]

                                   [             2          2          2  ]
                                   [ sqrt((c - z)  + (b - y)  + (a - x) ) ]
                           Col 2 = [                                      ]
                                   [             2          2          2  ]
                                   [ sqrt((f - z)  + (e - y)  + (d - x) ) ]

 -- Option variable: matrix_element_transpose
     Default value: false

     matrix_element_transpose is the operation applied to each element of a matrix when it is transposed. matrix_element_mult can be assigned any unary operator. The assigned value may be the name of an operator enclosed in quote marks, the name of a function, or a lambda expression.

     When matrix_element_transpose equals transpose, the transpose function is applied to every element. When matrix_element_transpose equals nonscalars, the transpose function is applied to every nonscalar element. If some element is an atom, the nonscalars option applies transpose only if the atom is declared nonscalar, while the transpose option always applies transpose.

     The default value, false, means no operation is applied.

     See also matrix_element_add and matrix_element_mult.

     Examples:

          (%i1) declare (a, nonscalar)$
          (%i2) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o2)                   [              ]
                                  [      b       ]
          (%i3) matrix_element_transpose: nonscalars$
          (%i4) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o4)                   [              ]
                                  [      b       ]
          (%i5) matrix_element_transpose: transpose$
          (%i6) transpose ([a, b]);
                                  [ transpose(a) ]
          (%o6)                   [              ]
                                  [ transpose(b) ]
          (%i7) matrix_element_transpose: lambda ([x], realpart(x)
                - %i*imagpart(x))$
          (%i8) m: matrix ([1 + 5*%i, 3 - 2*%i], [7*%i, 11]);
                               [ 5 %i + 1  3 - 2 %i ]
          (%o8)                [                    ]
                               [   7 %i       11    ]
          (%i9) transpose (m);
                                [ 1 - 5 %i  - 7 %i ]
          (%o9)                 [                  ]
                                [ 2 %i + 3    11   ]

 -- Function: mattrace (M)
     Returns the trace (that is, the sum of the elements on the main diagonal) of the square matrix M.

     mattrace is called by ncharpoly, an alternative to Maxima's charpoly.

     load ("nchrpl") loads this function.

 -- Function: minor (M, I, J)
     Returns the I, J minor of the matrix M.  That is, M with row I and column J removed.

 -- Function: ncexpt (A, B)
     If a non-commutative exponential expression is too wide to be displayed as A^^B it appears as ncexpt (A,B).

     ncexpt is not the name of a function or operator; the name only appears in output, and is not recognized in input.

 -- Function: ncharpoly (M, X)
     Returns the characteristic polynomial of the matrix M with respect to X.  This is an alternative to Maxima's charpoly.

     ncharpoly works by computing traces of powers of the given matrix, which are known to be equal to sums of powers of the roots of the characteristic polynomial.  From these quantities the symmetric functions of the roots can be calculated, which are nothing more than the coefficients of the characteristic polynomial.  charpoly works by forming the determinant of X * ident [n] - a.  Thus ncharpoly wins, for example, in the case of large dense matrices filled with integers, since it avoids polynomial arithmetic altogether.

     load ("nchrpl") loads this file.

 -- Function: newdet (M, N)
     Computes the determinant of the matrix or array M by the Johnson-Gentleman tree minor algorithm. The argument N is the order; it is optional if M is a matrix.

 -- Declaration: nonscalar
     Makes atoms behave as does a list or matrix with respect to the dot operator.

 -- Function: nonscalarp (EXPR)
     Returns true if EXPR is a non-scalar, i.e., it contains atoms declared as non-scalars, lists, or matrices.

 -- Function: permanent (M, N)
     Computes the permanent of the matrix M.  A permanent is like a determinant but with no sign changes.

 -- Function: rank (M)
     Computes the rank of the matrix M.  That is, the order of the largest non-singular subdeterminant of M.

     RANK may return the wrong answer if it cannot determine that a matrix element that is equivalent to zero is indeed so.

 -- Option variable: ratmx
     Default value: false

     When ratmx is false, determinant and matrix addition, subtraction, and multiplication are performed in the representation of the matrix elements and cause the result of matrix inversion to be left in general representation.

     When ratmx is true, the 4 operations mentioned above are performed in CRE form and the result of matrix inverse is in CRE form.  Note that this may cause the elements to be expanded (depending on the setting of ratfac) which might not always be desired.

 -- Function: row (M, I)
     Returns the I'th row of the matrix M. The return value is a matrix.

 -- Option variable: scalarmatrixp
     Default value: true

     When scalarmatrixp is true, then whenever a 1 x 1 matrix is produced as a result of computing the dot product of matrices it is simplified to a scalar, namely the sole element of the matrix.

     When scalarmatrixp is all, then all 1 x 1 matrices are simplified to scalars.

     When scalarmatrixp is false, 1 x 1 matrices are not simplified to scalars.

 -- Function: scalefactors (COORDINATETRANSFORM)
     Here coordinatetransform evaluates to the form [[expression1, expression2, ...], indeterminate1, indeterminat2, ...], where indeterminate1, indeterminate2, etc. are the curvilinear coordinate variables and where a set of rectangular Cartesian components is given in terms of the curvilinear coordinates by [expression1, expression2, ...]. coordinates is set to the vector [indeterminate1, indeterminate2,...], and dimension is set to the length of this vector.  SF[1], SF[2], ..., SF[DIMENSION] are set to the coordinate scale factors, and sfprod is set to the product of these scale factors.  Initially, coordinates is [X, Y, Z], dimension is 3, and SF[1]=SF[2]=SF[3]=SFPROD=1, corresponding to 3-dimensional rectangular Cartesian coordinates. To expand an expression into physical components in the current coordinate system, there is a function with usage of the form

 -- Function: setelmx (X, I, J, M)
     Assigns X to the (I, J)'th element of the matrix M, and returns the altered matrix.

     M [I, J]: X has the same effect, but returns X instead of M.

 -- Function: similaritytransform (M)
 -- Function: simtran (M)
     similaritytransform computes a similarity transform of the matrix M. It returns a list which is the output of the uniteigenvectors command.  In addition if the flag nondiagonalizable is false two global matrices leftmatrix and rightmatrix are computed. These matrices have the property that leftmatrix . M . rightmatrix is a diagonal matrix with the eigenvalues of M on the diagonal.  If nondiagonalizable is true the left and right matrices are not computed.

     If the flag hermitianmatrix is true then leftmatrix is the complex conjugate of the transpose of rightmatrix.  Otherwise leftmatrix is the inverse of rightmatrix.

     rightmatrix is the matrix the columns of which are the unit eigenvectors of M.  The other flags (see eigenvalues and eigenvectors) have the same effects since similaritytransform calls the other functions in the package in order to be able to form rightmatrix.

     load ("eigen") loads this function.

     simtran is a synonym for similaritytransform.

 -- Option variable: sparse
     Default value: false

     When sparse is true, and if ratmx is true, then determinant will use special routines for computing sparse determinants.

 -- Function: submatrix (I_1, ..., I_M, M, J_1, ..., J_N)
 -- Function: submatrix (I_1, ..., I_M, M)
 -- Function: submatrix (M, J_1, ..., J_N)
     Returns a new matrix composed of the matrix M with rows I_1, ..., I_M deleted, and columns J_1, ..., J_N deleted.

 -- Function: transpose (M)
     Returns the transpose of M.

     If M is a matrix, the return value is another matrix N such that N[i,j] = M[j,i].

     If M is a list, the return value is a matrix N of length (m) rows and 1 column, such that N[i,1] = M[i].

     Otherwise M is a symbol, and the return value is a noun expression 'transpose (M).

 -- Function: triangularize (M)
     Returns the upper triangular form of the matrix M, as produced by Gaussian elimination. The return value is the same as echelon, except that the leading nonzero coefficient in each row is not normalized to 1.

     lu_factor and cholesky are other functions which yield triangularized matrices.

          (%i1) M: matrix ([3, 7, aa, bb], [-1, 8, 5, 2], [9, 2, 11, 4]);
                                 [  3   7  aa  bb ]
                                 [                ]
          (%o1)                  [ - 1  8  5   2  ]
                                 [                ]
                                 [  9   2  11  4  ]
          (%i2) triangularize (M);
                       [ - 1   8         5            2      ]
                       [                                     ]
          (%o2)        [  0   - 74     - 56         - 22     ]
                       [                                     ]
                       [  0    0    626 - 74 aa  238 - 74 bb ]

 -- Function: uniteigenvectors (M)
 -- Function: ueivects (M)
     Computes unit eigenvectors of the matrix M. The return value is a list of lists, the first sublist of which is the output of the eigenvalues command, and the other sublists of which are the unit eigenvectors of the matrix corresponding to those eigenvalues respectively.

     The flags mentioned in the description of the eigenvectors command have the same effects in this one as well.

     When knowneigvects is true, the eigen package assumes that the eigenvectors of the matrix are known to the user and are stored under the global name listeigvects.  listeigvects should be set to a list similar to the output of the eigenvectors command.

     If knowneigvects is set to true and the list of eigenvectors is given the setting of the flag nondiagonalizable may not be correct.  If that is the case please set it to the correct value.  The author assumes that the user knows what he is doing and will not try to diagonalize a matrix the eigenvectors of which do not span the vector space of the appropriate dimension.

     load ("eigen") loads this function.

     ueivects is a synonym for uniteigenvectors.

 -- Function: unitvector (X)
 -- Function: uvect (X)
     Returns X/norm(X); this is a unit vector in the same direction as X.

     load ("eigen") loads this function.

     uvect is a synonym for unitvector.

 -- Function: vectorsimp (EXPR)
     Applies simplifications and expansions according to the following global flags:

     expandall, expanddot, expanddotplus, expandcross, expandcrossplus, expandcrosscross, expandgrad, expandgradplus, expandgradprod, expanddiv, expanddivplus, expanddivprod, expandcurl, expandcurlplus, expandcurlcurl, expandlaplacian, expandlaplacianplus, and expandlaplacianprod.

     All these flags have default value false. The plus suffix refers to employing additivity or distributivity.  The prod suffix refers to the expansion for an operand that is any kind of product.

     expandcrosscross
          Simplifies p ~ (q ~ r) to (p . r)*q - (p . q)*r.

     expandcurlcurl
          Simplifies curl curl p to grad div p + div grad p.

     expandlaplaciantodivgrad
          Simplifies laplacian p to div grad p.

     expandcross
          Enables expandcrossplus and expandcrosscross.

     expandplus
          Enables expanddotplus, expandcrossplus, expandgradplus, expanddivplus, expandcurlplus, and expandlaplacianplus.

     expandprod
          Enables expandgradprod, expanddivprod, and expandlaplacianprod.

     These flags have all been declared evflag.

 -- Option variable: vect_cross
     Default value: false

     When vect_cross is true, it allows DIFF(X~Y,T) to work where ~ is defined in SHARE;VECT (where VECT_CROSS is set to true, anyway.)

 -- Function: zeromatrix (M, N)
     Returns an M by N matrix, all elements of which are zero.

 -- Special symbol: [
 -- Special symbol: ]
     [ and ] mark the beginning and end, respectively, of a list.

     [ and ] also enclose the subscripts of a list, array, hash array, or array function.

     Examples:

          (%i1) x: [a, b, c];
          (%o1)                       [a, b, c]
          (%i2) x[3];
          (%o2)                           c
          (%i3) array (y, fixnum, 3);
          (%o3)                           y
          (%i4) y[2]: %pi;
          (%o4)                          %pi
          (%i5) y[2];
          (%o5)                          %pi
          (%i6) z['foo]: 'bar;
          (%o6)                          bar
          (%i7) z['foo];
          (%o7)                          bar
          (%i8) g[k] := 1/(k^2+1);
                                            1
          (%o8)                     g  := ------
                                     k     2
                                          k  + 1
          (%i9) g[10];
                                          1
          (%o9)                          ---
                                         101