MatrixOps 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.
M[i,j] read the (i,j)-entry of matrix M
SetEntry(M,i,j, val) set the (i,j)-entry of matrix M
GetRow(M,i) return the i-th row of M as a vector<RingElem>
GetCol(M,j) return the j-th col of M as a vector<RingElem>
GetRows(M) return the rows of M as a vector<vector<RingElem>>
GetCols(M) return the cols of M as a vector<vector<RingElem>>
FlattenByRows(M) return entries of M in a vector<RingElem> in order 1st row, 2nd row, etc
FlattenByCols(M) return entries of M in a vector<RingElem> in order 1st col, 2nd col, etc
Note that GetRow, GetCol, GetRows, GetCols, FlattenByRows and FlattenByCols make copies of the matrix entries.
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 FrobeniusNormSq gives the
square of the Frobenius norm (so that the value surely lies in the
same ring).
FrobeniusNormSq(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)
IsZeroDet(M) -- equivalent to IsZero(det(M)) (but may be faster)
HadamardBoundSq(M) -- computes row and column bounds in a struct (fields myRowBoundSq and myColBoundSq)
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".
rref(M) -- compute a reduced row echelon form of M (orig. matrix is unchanged); matrix must be over a field
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)
LinKerZZ(M) -- solve for x the linear system M*x = 0; returns a DenseMatrix whose columns are a ZZ-base for integer points in 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.
det2x2(M) -- only for 2x2 matrices
det3x3(M) -- only for 3x3 matrices
det4x4(M) -- only for 4x4 matrices
det5x5(M) -- only for 5x5 matrices
DetByGauss(M) -- matrix must be over an integral domain
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