+@node Matrices, Indexed objects, Relations, Basic Concepts
+@c node-name, next, previous, up
+@section Matrices
+@cindex @code{matrix} (class)
+
+A @dfn{matrix} is a two-dimensional array of expressions. The elements of a
+matrix with @math{m} rows and @math{n} columns are accessed with two
+@code{unsigned} indices, the first one in the range 0@dots{}@math{m-1}, the
+second one in the range 0@dots{}@math{n-1}.
+
+There are a couple of ways to construct matrices, with or without preset
+elements:
+
+@example
+matrix::matrix(unsigned r, unsigned c);
+matrix::matrix(unsigned r, unsigned c, const lst & l);
+ex lst_to_matrix(const lst & l);
+ex diag_matrix(const lst & l);
+@end example
+
+The first two functions are @code{matrix} constructors which create a matrix
+with @samp{r} rows and @samp{c} columns. The matrix elements can be
+initialized from a (flat) list of expressions @samp{l}. Otherwise they are
+all set to zero. The @code{lst_to_matrix()} function constructs a matrix
+from a list of lists, each list representing a matrix row. Finally,
+@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
+elements. Note that the last two functions return expressions, not matrix
+objects.
+
+Matrix elements can be accessed and set using the parenthesis (function call)
+operator:
+
+@example
+const ex & matrix::operator()(unsigned r, unsigned c) const;
+ex & matrix::operator()(unsigned r, unsigned c);
+@end example
+
+It is also possible to access the matrix elements in a linear fashion with
+the @code{op()} method. But C++-style subscripting with square brackets
+@samp{[]} is not available.
+
+Here are a couple of examples that all construct the same 2x2 diagonal
+matrix:
+
+@example
+@{
+ symbol a("a"), b("b");
+ ex e;
+
+ matrix M(2, 2);
+ M(0, 0) = a;
+ M(1, 1) = b;
+ e = M;
+
+ e = matrix(2, 2, lst(a, 0, 0, b));
+
+ e = lst_to_matrix(lst(lst(a, 0), lst(0, b)));
+
+ e = diag_matrix(lst(a, b));
+
+ cout << e << endl;
+ // -> [[a,0],[0,b]]
+@}
+@end example
+
+@cindex @code{transpose()}
+@cindex @code{inverse()}
+There are three ways to do arithmetic with matrices. The first (and most
+efficient one) is to use the methods provided by the @code{matrix} class:
+
+@example
+matrix matrix::add(const matrix & other) const;
+matrix matrix::sub(const matrix & other) const;
+matrix matrix::mul(const matrix & other) const;
+matrix matrix::mul_scalar(const ex & other) const;
+matrix matrix::pow(const ex & expn) const;
+matrix matrix::transpose(void) const;
+matrix matrix::inverse(void) const;
+@end example
+
+All of these methods return the result as a new matrix object. Here is an
+example that calculates @math{A*B-2*C} for three matrices @math{A}, @math{B}
+and @math{C}:
+
+@example
+@{
+ matrix A(2, 2, lst(1, 2, 3, 4));
+ matrix B(2, 2, lst(-1, 0, 2, 1));
+ matrix C(2, 2, lst(8, 4, 2, 1));
+
+ matrix result = A.mul(B).sub(C.mul_scalar(2));
+ cout << result << endl;
+ // -> [[-13,-6],[1,2]]
+ ...
+@}
+@end example
+
+@cindex @code{evalm()}
+The second (and probably the most natural) way is to construct an expression
+containing matrices with the usual arithmetic operators and @code{pow()}.
+For efficiency reasons, expressions with sums, products and powers of
+matrices are not automatically evaluated in GiNaC. You have to call the
+method
+
+@example
+ex ex::evalm() const;
+@end example
+
+to obtain the result:
+
+@example
+@{
+ ...
+ ex e = A*B - 2*C;
+ cout << e << endl;
+ // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
+ cout << e.evalm() << endl;
+ // -> [[-13,-6],[1,2]]
+ ...
+@}
+@end example
+
+The non-commutativity of the product @code{A*B} in this example is
+automatically recognized by GiNaC. There is no need to use a special
+operator here. @xref{Non-commutative objects}, for more information about
+dealing with non-commutative expressions.
+
+Finally, you can work with indexed matrices and call @code{simplify_indexed()}
+to perform the arithmetic:
+
+@example
+@{
+ ...
+ idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
+ e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
+ cout << e << endl;
+ // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
+ cout << e.simplify_indexed() << endl;
+ // -> [[-13,-6],[1,2]].i.j
+@}
+@end example
+
+Using indices is most useful when working with rectangular matrices and
+one-dimensional vectors because you don't have to worry about having to
+transpose matrices before multiplying them. @xref{Indexed objects}, for
+more information about using matrices with indices, and about indices in
+general.
+
+The @code{matrix} class provides a couple of additional methods for
+computing determinants, traces, and characteristic polynomials:
+
+@example
+ex matrix::determinant(unsigned algo = determinant_algo::automatic) const;
+ex matrix::trace(void) const;
+ex matrix::charpoly(const symbol & lambda) const;
+@end example
+
+The @samp{algo} argument of @code{determinant()} allows to select between
+different algorithms for calculating the determinant. The possible values
+are defined in the @file{flags.h} header file. By default, GiNaC uses a
+heuristic to automatically select an algorithm that is likely to give the
+result most quickly.
+
+
+@node Indexed objects, Non-commutative objects, Matrices, Basic Concepts