Linear Algebra Functions

Overview

The DataFrame library provides a comprehensive set of linear algebra types and functions for working with vectors and matrices. These types are fixed-size, stack-allocated, and designed for performance in scientific computing. All algebra functions operate on Serie containers and have bind_* versions for pipeline usage.

Types

Vectors

Vector Types


// Generic fixed-size vector
template <typename T, size_t N>
using Vector = std::array<T, N>;

// Common aliases
using Vector2D = Vector<double, 2>;
using Vector3D = Vector<double, 3>;
using Vector4D = Vector<double, 4>;

Full Matrices

Full Matrix Types


// Generic NxN full matrix (stored as flat array of N*N elements)
template <typename T, size_t N>
using FullMatrix = std::array<T, N * N>;

// Common aliases
using Matrix2D = FullMatrix<double, 2>;   // 4 components
using Matrix3D = FullMatrix<double, 3>;   // 9 components
using Matrix4D = FullMatrix<double, 4>;   // 16 components

Symmetric Matrices

Symmetric Matrix Types


// Symmetric NxN matrix (stores only upper triangle)
// 2x2: 3 components, 3x3: 6 components, 4x4: 10 components
template <typename T, size_t N>
using SymmetricMatrix = std::array<T, N * (N + 1) / 2>;

// Common aliases
using SMatrix2D = SymmetricMatrix<double, 2>;  // 3 components
using SMatrix3D = SymmetricMatrix<double, 3>;  // 6 components (Stress3D, Strain3D)
using SMatrix4D = SymmetricMatrix<double, 4>;  // 10 components

// Stress and strain tensors
using Stress2D = SymmetricMatrix<double, 2>;
using Strain2D = SymmetricMatrix<double, 2>;
using Stress3D = SymmetricMatrix<double, 3>;
using Strain3D = SymmetricMatrix<double, 3>;

Vector Operations

Cross Product


// Cross product of two 3D vector Series
Serie<Vector3D> cross(const Serie<Vector3D> &a, const Serie<Vector3D> &b);

// Pipeline version
auto bind_cross(const Serie<Vector3D> &b);

Dot Product


// Dot product of two vector Series
template <typename T, size_t N>
Serie<T> dot(const Serie<Vector<T,N>> &a, const Serie<Vector<T,N>> &b);

// Pipeline version
template <typename T, size_t N>
auto bind_dot(const Serie<Vector<T,N>> &b);

Norm

Vector Norm (Magnitude)


// Compute the magnitude (Euclidean norm) of each vector
template <typename T, size_t N>
Serie<T> norm(const Serie<Vector<T,N>> &serie);

// Pipeline version
template <typename T, size_t N>
auto bind_norm();

Vector Operations Example


df::Serie<Vector3D> a{{1, 0, 0}, {0, 1, 0}, {1, 1, 0}};
df::Serie<Vector3D> b{{0, 1, 0}, {0, 0, 1}, {0, 1, 0}};

// Cross product
auto c = df::algebra::cross(a, b);
// c[0] = {0, 0, 1}, c[1] = {1, 0, 0}, c[2] = {0, 0, 1}

// Dot product
auto d = df::algebra::dot(a, b);
// d = {0, 0, 1}

// Norm
auto magnitudes = df::algebra::norm(a);
// magnitudes = {1.0, 1.0, 1.414...}

// Pipeline usage
auto result = a | df::algebra::bind_cross(b);
auto norms  = a | df::algebra::bind_norm<double, 3>();

Matrix Operations

Determinant


// Determinant of full matrices
template <typename T, size_t N>
Serie<T> det(const Serie<FullMatrix<T,N>> &matrices);

// Determinant of symmetric matrices
template <typename T, size_t N>
Serie<T> det(const Serie<SymmetricMatrix<T,N>> &matrices);

// Pipeline versions
auto bind_det();

Inverse

Matrix Inverse


// Inverse for 1x1, 2x2, and 3x3 full matrices
template <typename T, size_t N>
Serie<FullMatrix<T,N>> inv(const Serie<FullMatrix<T,N>> &matrices);

// Inverse for symmetric matrices
template <typename T, size_t N>
Serie<SymmetricMatrix<T,N>> inv(const Serie<SymmetricMatrix<T,N>> &matrices);

// Pipeline version
auto bind_inv();

Transpose

Matrix Transpose


// Transpose for 2x2, 3x3, and 4x4 full matrices
template <typename T, size_t N>
Serie<FullMatrix<T,N>> transpose(const Serie<FullMatrix<T,N>> &matrices);

// Pipeline version
auto bind_transpose();

Solve Linear System

Solve Ax = b


// Solve the linear system Ax = b for each element
template <typename T, size_t N>
Serie<Vector<T,N>> solve(
    const Serie<FullMatrix<T,N>> &A,
    const Serie<Vector<T,N>> &b
);

Matrix Operations Example


// 2x2 matrix operations
df::Serie<Matrix2D> mats{
    {1, 2, 3, 4},    // [[1, 2], [3, 4]]
    {5, 6, 7, 8}     // [[5, 6], [7, 8]]
};

auto determinants = df::algebra::det(mats);
// determinants = {-2, -2}

auto inverses = df::algebra::inv(mats);
auto transposed = df::algebra::transpose(mats);

// Symmetric matrix (stress tensor) eigen decomposition
df::Serie<SMatrix3D> stresses{
    {1, 0, 0, 1, 0, 1},   // identity-like
    {2, 1, 0, 3, 0, 1}
};
auto dets = df::algebra::det(stresses);

// Solve linear system: Ax = b
df::Serie<Matrix3D> A{
    {2, 1, 0, 1, 3, 1, 0, 1, 2}
};
df::Serie<Vector3D> b{{1, 2, 3}};
auto x = df::algebra::solve(A, b);

Eigen Decomposition

Eigen Functions (Symmetric Matrices Only)


// Eigenvalues of symmetric matrices (returned sorted)
template <typename T, size_t N>
Serie<Vector<T,N>> eigenValues(const Serie<SymmetricMatrix<T,N>> &matrices);

// Eigenvectors of symmetric matrices
template <typename T, size_t N>
Serie<FullMatrix<T,N>> eigenVectors(const Serie<SymmetricMatrix<T,N>> &matrices);

// Both eigenvalues and eigenvectors at once
template <typename T, size_t N>
std::pair<Serie<Vector<T,N>>, Serie<FullMatrix<T,N>>>
eigenSystem(const Serie<SymmetricMatrix<T,N>> &matrices);

// Pipeline versions
auto bind_eigenValues();
auto bind_eigenVectors();
auto bind_eigenSystem();

Eigen Decomposition Example


// Stress tensor eigen decomposition
df::Serie<SMatrix3D> stresses{
    {10.0, 0.0, 0.0, 5.0, 0.0, 3.0},
    {20.0, 1.0, 0.5, 15.0, 0.0, 8.0}
};

// Get eigenvalues (principal stresses)
auto principals = df::algebra::eigenValues(stresses);
// principals[i] = {sigma1, sigma2, sigma3} sorted

// Get eigenvectors (principal directions)
auto directions = df::algebra::eigenVectors(stresses);

// Get both at once (more efficient)
auto [values, vectors] = df::algebra::eigenSystem(stresses);

// Pipeline usage
auto evals = stresses | df::algebra::bind_eigenValues<double, 3>();

Complete Example

Stress Analysis Example


#include <dataframe/Serie.h>
#include <dataframe/algebra.h>
#include <iostream>

int main() {
    // Define stress tensors at mesh nodes (symmetric 3x3)
    // Components: {xx, xy, xz, yy, yz, zz}
    df::Serie<Stress3D> stresses{
        {100.0, 10.0, 0.0, 50.0, 5.0, 30.0},
        {200.0, 20.0, 5.0, 150.0, 10.0, 80.0},
        {80.0,  5.0,  0.0, 40.0, 2.0, 25.0}
    };

    // Compute principal stresses (eigenvalues)
    auto [principals, directions] = df::algebra::eigenSystem(stresses);

    for (size_t i = 0; i < stresses.size(); ++i) {
        std::cout << "Node " << i << ": "
                  << "sigma1=" << principals[i][0]
                  << " sigma2=" << principals[i][1]
                  << " sigma3=" << principals[i][2] << "\n";
    }

    // Compute determinant of each stress tensor
    auto dets = df::algebra::det(stresses);

    // Compute inverse
    auto inv_stresses = df::algebra::inv(stresses);

    // Vector operations: displacement field
    df::Serie<Vector3D> displacements{
        {0.1, 0.2, 0.0},
        {0.3, 0.1, 0.05},
        {0.05, 0.15, 0.0}
    };

    auto magnitudes = df::algebra::norm(displacements);

    std::cout << "Displacement magnitudes: ";
    for (size_t i = 0; i < magnitudes.size(); ++i) {
        std::cout << magnitudes[i] << " ";
    }
    std::cout << std::endl;

    return 0;
}