MatrixOperations gathers together a number of operations on matrices; in
most cases these operations are happy to accept a MatrixView
(see MatrixView) as argument.
When not specified, a matrix argument is of type ConstMatrixView.
There are two ways of multiplying two matrices together. The infix
operators return a DenseMatrix; the procedural version may be
slightly faster than the infix operator.
mul(matrix& lhs, M1, M2) -- a procedure equivalent to lhs = M1*M2;, note that lhs might be a SparseMatrix (not yet implemented)
operator*(M1, M2) -- the product M1*M2
operator+(M1, M2) -- the sum M1+M2
operator-(M1, M2) -- the difference M1-M2
power(M, n) compute n-th power of M; if n is negative then M must be invertible
operator*(n, M1) -- scalar multiple of M1 by n (integer or RingElem)
operator*(M1, n) -- scalar multiple of M1 by n (integer or RingElem)
operator/(M1, n) -- scalar multiple of M1 by 1/n (where n is integer or RingElem)
operator-(M1) -- scalar multiple of M1 by -1
Here are some matrix norms. The result is an element of the ring
containing the matrix elements. Note that FrobeniusNorm2 gives the
square of the Frobenius norm (so that the value surely lies in the
same ring).
FrobeniusNorm2(M) -- the square of the Frobenius norm
OperatorNormInfinity(M) -- the infinity norm, ring must be ordered
OperatorNorm1(M) -- the one norm, ring must be ordered
Here are some fairly standard functions on matrices.
det(M) -- determinant of M (M must be square)
rk(M) -- rank of M (the base ring must be an integral domain)
inverse(M) -- inverse of M as a DenseMatrix
adj(M) -- classical adjoint of M as a DenseMatrix; sometimes called "adjugate".
PseudoInverse(M) -- PseudoInverse of M as a DenseMatrix.
I suspect that it requires that the matrix be of full rank.
LinSolve(M,rhs) -- solve for x the linear system M*x = rhs; result is a DenseMatrix; if no soln exists, result is the 0-by-0 matrix
LinKer(M) -- solve for x the linear system M*x = 0; returns a DenseMatrix whose columns are a base for ker(M)
Here are some standard operations where the method used is specified explicitly. It would usually be better to use the generic operations above, as those should automatically select the most appropriate method for the given matrix.
void det2x2(RingElem& d, M) -- for 2x2 matrices
void det3x3(RingElem& d, M) -- for 3x3 matrices
void DetByGauss(RingElem& d, M)
RankByGauss(std::vector<long>& IndepRows, M)
InverseByGauss(M) -- some restrictions (needs gcd)
AdjointByDetOfMinors(M)
AdjointByInverse(M) -- base ring must be integral domain
LinSolveByGauss(M,rhs) -- solve a linear system using gaussian elimination
(base ring must be a field), result is a DenseMatrix
LinSolveByHNF(M,rhs) -- solve a linear system using Hermite NormalForm
(base ring must be a PID), result is a DenseMatrix
LinSolveByModuleRepr(M,rhs) -- solve a linear system using module element representation, result is a DenseMatrix
void GrammSchmidtRows(MatrixView& M) -- NYI
void GrammSchmidtRows(MatrixView& M, long row) -- NYI
Most impls are quite straightforward.
power is slightly clever with its iterative impl of binary powering.
LinSolveByGauss is a little complicated because it tries to handle all
cases (e.g. full rank or not, square or more rows than cols or more cols than rows)
Can we make a common "gaussian elimination" impl which is called by the various algorithms needing it, rather than having several separate implementations?
Is the procedure mul really any faster than the infix operator?
2012
2011