Power $A^n$ (binary exponentiation). $n < 0$ computes $(A^{-1})^{|n|}$
Note: The ^ operator has low precedence in C++. Use parentheses in expressions: (A ^ 'T') * B
Matrix<double> A = {{1, 2}, {3, 4}};
Matrix<double> B = {{5, 6}, {7, 8}};
auto C = A * B; // {{19, 22}, {43, 50}}
Vector<double> v = {1, 2};
auto w = A * v; // {5, 11}
auto At = A ^ 'T'; // transpose
auto Ainv = A ^ (-1); // inverse
auto A3 = A ^ 3; // A*A*A
Matrix Operations
Matrix transpose() const;
T determinant() const;
[[nodiscard]] Matrix inverse() const;
T trace() const;
T norm() const;
void zero();
void fill(const T& value);
void identity();
static Matrix identity(size_type size);
Matrix<double> A = {{1, 2}, {3, 4}};
std::cout << A.determinant() << std::endl; // -2
auto Ainv = A.inverse();
auto I = A * Ainv; // approximately identity matrix
Low-Level Free Functions
Output-parameter form. Customizable via the ComputePolicy template parameter.
template<...> void add(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void subtract(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void scale(const MatrixA& a, const Scalar& alpha, MatrixB& b);
template<...> void multiply(const MatrixA& a, const MatrixB& b, MatrixC& c);
template<...> void multiply_vector(const MatrixA& a, const VectorX& x, VectorY& y);
template<...> void transpose(const MatrixA& a, MatrixB& b);
template<typename MatrixA> auto trace(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto frobenius_norm(const MatrixA& a) -> typename MatrixA::value_type;
template<typename MatrixA> auto determinant(const MatrixA& a) -> typename MatrixA::value_type;
template<typename T> Matrix<T> inverse(const Matrix<T>& a);
Multiplication between StaticMatrix instances is size-checked at compile time:
StaticMatrix<T, M, P> = StaticMatrix<T, M, N> * StaticMatrix<T, N, P>.
Examples
Basic Matrix Creation and Operations
#include <math/core/matrix.hpp>
#include <iostream>
using namespace calx;
int main() {
// initializer_list construction
Matrix<double> A = {{1, 2, 3},
{4, 5, 6},
{7, 8, 10}};
// transpose, trace, determinant
auto At = A.transpose();
std::cout << "trace(A) = " << A.trace() << std::endl; // 16
std::cout << "det(A) = " << A.determinant() << std::endl; // -3
// matrix multiplication
Matrix<double> B = {{1, 0}, {0, 1}, {1, 1}};
auto C = A * B; // 3x2 result
print_matrix(C);
// scalar operations
auto D = 2.0 * A - A; // == A
// factory functions
auto I = eye<double>(3);
auto R = rotationZ<double>(M_PI / 6); // Z-axis 30 deg rotation
}
Solving Linear Systems with solve()
#include <math/core/matrix.hpp>
#include <math/linalg/solve.hpp>
#include <iostream>
using namespace calx;
int main() {
// Solve Ax = b
Matrix<double> A = {{2, 1, -1},
{-3, -1, 2},
{-2, 1, 2}};
Vector<double> b = {8, -11, -3};
auto x = solve(A, b);
std::cout << "x = " << x << std::endl; // {2, 3, -1}
// verify: A * x == b
auto residual = A * x - b;
std::cout << "residual norm = " << norm(residual) << std::endl; // ~0
}
Exact Inverse of a Hilbert Matrix with Matrix<Rational>
#include <math/core/matrix.hpp>
#include <math/core/mp/Rational.hpp>
#include <iostream>
using namespace calx;
int main() {
const int n = 4;
Matrix<Rational> H(n, n);
for (int i = 0; i < n; ++i)
for (int j = 0; j < n; ++j)
H(i, j) = Rational(1, i + j + 1);
// exact inverse -- no floating-point error
auto Hinv = H.inverse();
auto I = H * Hinv;
// verify: every diagonal element is exactly 1
for (int i = 0; i < n; ++i)
std::cout << "I(" << i << "," << i << ") = " << I(i, i) << std::endl;
// Output: 1/1 for all
print_matrix(Hinv);
}
Related Mathematical Background
The following articles explain the mathematical concepts underlying matrix operations.
Matrix Norms — Definitions and properties of various norms