matlib

A header-only C++20 matrix library with both host (CPU) and device (CUDA GPU) implementations, designed for neural network applications.

Features

• Dual Backend: Seamless host (CPU) and device (CUDA GPU) support
• CRTP Architecture: Efficient Matrix and MatrixView classes sharing a common base
• Column-Major Storage: Compatible with MATLAB/FORTRAN conventions
• Strided Views: Zero-copy slicing, transposing, reshaping operations
• Rvalue Optimization: Overloads for & (lvalue) and && (rvalue) enable efficient chaining without unnecessary copies
• CUDA Integration: cuBLAS-accelerated matrix multiplication, custom CUDA kernels
• RAII GPU Memory: GPUBuffer/GPUSpan hide CUDA allocation/deallocation details
• Neural Network Utilities: Xavier/He/LeCun initialization, activation functors with gradients, im2col/col2im for convolution

Quick Start

Build

cmake -B build -DUSE_CUDA=ON  # or OFF for CPU-only
cmake --build build

Basic Usage

#include "matlib.hh"

using namespace matlib;

// Create matrices
host::Matrix<float> A({128, 64});   // 128×64 matrix
host::Matrix<float> B({64, 32});    // 64×32 matrix

// Initialize with random values
A.SetRandXavier();  // Xavier initialization
B.SetRandHe();      // He initialization

// Matrix operations
host::Matrix<float> C = A * B;      // Matrix multiplication
host::Matrix<float> D = A + B;      // Element-wise addition (broadcasting)
host::Matrix<float> E = A.transpose();  // Transpose (zero-copy view)

// Slicing and reshaping
auto slice = A(0, 10);              // First 10 rows (view)
auto flat = A.flatten(2);           // Flatten first 2 dimensions

// Reductions
float sum = A.sum();
float avg = A.avg();
host::Matrix<float> row_sum = A.sum(1);  // Sum across rows

// Transfer to GPU (when USE_CUDA=ON)
device::Matrix<float> A_gpu = A.toDevice();
device::Matrix<float> B_gpu = B.toDevice();
device::Matrix<float> C_gpu = A_gpu * B_gpu;  // cuBLAS matmul
host::Matrix<float> C_result = C_gpu.toHost();

Architecture

Core Classes

Class	Description
host::Matrix<T>	Owning CPU matrix stored in std::vector<T>
host::MatrixView<T>	Non-owning CPU view with custom strides
device::Matrix<T>	Owning GPU matrix with CUDA memory
device::MatrixView<T>	Non-owning GPU view

All classes inherit from MatrixBase<Derived, T> using CRTP, sharing operations like:

• Element-wise arithmetic (+, -, *, /)
• Scalar operations (+=, *=, etc.)
• Transformations (map, abs, exp, log, sqrt)
• Reductions (sum, avg, max, min)
• Shape manipulation (reshape, flatten, squeeze, unsqueeze, transpose)
• Indexing and slicing

Memory Layout

Uses column-major order (first dimension is contiguous):

• Shape {M, N} means M rows, N columns
• Element (i, j) is at offset i + j*M
• Compatible with cuBLAS and traditional linear algebra

Views and Strides

Views enable zero-copy operations:

Matrix<float> A({4, 3});
auto T = A.transpose();      // View with swapped strides
auto S = A(0, 2);             // Slice first 2 columns
auto R = A.reshape({6, 2});   // Reshape (contiguous only)

Views track dim(), stride(), and is_contiguous(). Non-contiguous views use custom iterators; operations like copy() or to_contiguous() compact data when needed.

Rvalue Optimization

Many operations have separate overloads for lvalue and rvalue references:

// Lvalue: returns a View (zero-copy)
auto T = A.transpose();           // MatrixView<T>

// Rvalue: modifies and returns the moved object
auto T = std::move(A).transpose(); // Matrix<T> (modified in-place)

This enables efficient chaining without intermediate allocations:

// Single allocation, no copies
auto result = std::move(A)
    .reshape({128, 64})
    .transpose()
    .unsqueeze(0, {1});

Operations with rvalue optimization: transpose(), reshape(), flatten(), unsqueeze(), map(), scalarOperator().

In-place Operators with Auto Broadcasting

In-place operators (+=, -=, hadamard_()) automatically broadcast and cycle the smaller operand:

Matrix<float> A({100, 50});  // 5000 elements
Matrix<float> B({50});       // 100 elements (row vector)

A += B;  // B cycles 50 times to cover all 5000 elements

Important: This behavior is powerful but can be surprising:

• If A.size() > B.size(): B repeats until A is fully covered
• If A.size() < B.size(): A's iterator restarts mid-operation

Use carefully when dimensions don't match—ensure intentional broadcasting.

Error Handling Philosophy

This library is exception-free by design. All errors trigger std::terminate() immediately:

// error.hh
[[noreturn]] void error(const std::string& message) {
    std::cerr << "\033[31m[ERROR]\033[0m " << message << std::endl;
    std::terminate();
}

void check(bool assertion, const std::string& message) {
    if (!assertion) error(message);
}

Design rationale:

1. All errors are mathematical errors — dimension mismatches, invalid shapes, out-of-bounds access
2. Mathematical errors are fatal, not recoverable — they indicate bugs in algorithm design or implementation
3. No exception overhead — enables aggressive optimization, no stack unwinding, no try/catch cost

This philosophy assumes:

• Correct matrix dimensions are known at design time
• Invalid operations indicate programmer error, not runtime conditions
• Immediate termination is preferable to silent corruption or unexpected behavior

Typical error triggers:

• A * B with mismatched inner dimensions
• reshape() to incompatible size
• flatten() on non-contiguous view
• Slicing out of bounds
• dim and vector size mismatch in construction

If you need graceful error handling for dynamic dimensions, validate externally before calling library operations.

Device Matrix Contiguity Requirement

Warning: Unlike host matrices, many device (device::Matrix/device::MatrixView) operations require contiguity and will terminate if violated:

device::Matrix<float> A({64, 32});
auto T = A.transpose();  // Creates non-contiguous view

T += B;  // TERMINATES: "non contiguous += not implemented"
T *= 2.0f;  // TERMINATES: "non contiguous *= not implemented"
T = 0.0f;  // TERMINATES: "non contiguous =T not implemented"

Operations requiring contiguity on device:

Operation	Contiguity Required
A += B, A -= B	Both A and B
A.hadamard_(B)	A
A += c, A -= c, A *= c, A /= c	A
A = c (scalar assignment)	A
A = B (matrix assignment)	A (B can be non-contiguous)
A.flatten()	A

Workaround: Convert to contiguous before operating:

auto T = A.transpose();
device::Matrix<float> T_contiguous = T.copy();  // Or: Matrix<float>(T)
T_contiguous += B;  // OK

This limitation exists because GPU kernels for strided access are more complex and often less efficient. For performance-critical code, maintain contiguous layout whenever possible.

Matrix Multiplication

host::Matrix<float> C = A * B;           // Standard matmul
host::Matrix<float> C = MatMul_AtB(A, B); // Aᵀ × B
host::Matrix<float> C = MatMul_ABt(A, B); // A × Bᵀ

device::Matrix<float> C = A_gpu * B_gpu; // cuBLAS batched GEMM

Supports batched multiplication with broadcasting:

• (M×N, batch=A) * (N×P, batch=B) → (M×P, batch=max(A,B))

Element-wise Operations

auto C = A + B;      // Addition (broadcasting)
auto C = A - B;      // Subtraction
auto C = A.hadamard(B);  // Element-wise multiplication
auto C = A.abs();    // Absolute value
auto C = A.exp();    // Exponential
auto C = A.log();    // Natural log
auto C = A.sqrt();   // Square root
auto C = -A;         // Negation

Reductions

float s = A.sum();           // Sum all elements
float m = A.max();           // Maximum element
auto row_sum = A.sum(1);     // Sum along dimension 0 (rows)
auto batch_avg = A.avg(2);   // Average first 2 dimensions

Shape Operations

auto flat = A.flatten();          // 1D view
auto flat = A.flatten(2);         // Flatten first 2 dims
auto sq = A.squeeze();            // Remove size-1 dimensions
auto exp = A.unsqueeze(1, {3});   // Insert dimension
auto T = A.transpose(0, 1);       // Swap dimensions
auto R = A.reshape({M, N});       // New shape

Initialization

A.SetRand();           // Uniform [-1, 1]
A.SetRandXavier();     // Xavier/Glorot normal
A.SetRandHe();         // He/Kaiming normal
A.SetRandLecun();      // LeCun normal
A.SetRandBinary(0.5);  // Random ±scale

Convolution Helpers

// im2col: Extract patches for convolution
auto col = im2col(image, {kW, kH});  // [kW*kH*Cin, Wout*Hout, Batch]

// col2im: Accumulate patches back to image
auto img = col2im(col, {kW, kH}, {W, H, Cin});

CUDA Backend

When USE_CUDA=ON, the device namespace provides GPU-accelerated operations:

• cuBLAS GEMM: Batched matrix multiplication with strided batching
• Custom Kernels: Element-wise ops, reductions, broadcasting
• Device Management: Automatic cuBLAS handle management

RAII GPU Memory

GPUBuffer<T> and GPUSpan<T> provide safe GPU memory handling:

// GPUBuffer: owning RAII wrapper
GPUBuffer<float> buf(1024);           // cudaMalloc on construction
buf.to_device(host_span);             // H2D copy
auto host_data = buf.to_host();       // D2H copy
// cudaFree called automatically on destruction

// GPUSpan: non-owning view (like std::span for GPU)
GPUSpan<float> view = buf.span();     // No allocation, just pointer+size
view.to_host();                       // Can still download data

Key features:

• Automatic cudaMalloc/cudaFree in constructors/destructors
• Move semantics for efficient transfers
• Implicit conversion from GPUBuffer to GPUSpan
• Handles both contiguous and strided copies

Neural Network Activation Functors

extra_functors.hh provides activation functions with gradient support for both host and device:

using namespace matlib::host::functors;

// Forward pass
auto activated = A.map(SigmoidOp());

// Backward pass (gradient w.r.t. output y)
auto grad_y = activated.map(DSigmoidYOp());

// Each activation has its gradient as a nested type
// SigmoidOp::GradY = DSigmoidYOp
// TanhOp::GradY = DTanhYOp
// ReLUOp::GradY = DReLUYOp

Available activations:

Forward	Gradient	Notes
SigmoidOp	DSigmoidYOp	σ(x) = 1/(1+e^(-x))
TanhOp	DTanhYOp	tanh(x)
ReLUOp	DReLUYOp	max(0, x)
LeakyReLUOp	DLeakyReLUYOp	x<0 → 0.01x
SoftplusOp	DSoftplusYOp	log(1+e^x)
SwishOp	DSwishXYOp	x·sigmoid(x), needs both x and y

Utility functors:

• ClipOp<Scale>: Clamp to [-Scale, Scale]
• InvSqrtStableOp<eps>: 1/sqrt(max(x, eps)) for safe normalization
• LogStableOp<eps>: log(max(x, eps)) for numerical stability

Device functors are identical but decorated with __device__ for CUDA kernels.

GPU Operations

device::Matrix<float> A_gpu = host_A.toDevice();
device::Matrix<float> B_gpu({64, 32});
B_gpu.SetRandHe();

// GPU operations
auto C_gpu = A_gpu * B_gpu;          // cuBLAS matmul
auto D_gpu = A_gpu + B_gpu;          // CUDA kernel (broadcasting)
auto E_gpu = A_gpu.sum(1);           // Group reduction
auto F_gpu = A_gpu.map(functors::ReLUOp());  // Element-wise activation

host::Matrix<float> result = C_gpu.toHost();

File Structure

matlib/
├── matlib.hh              # Main include
├── common.hh              # Shared utilities
├── error.hh               # Error handling
├── transfer.hh            # Host↔Device transfer
├── host/
│   ├── core.hh            # Host module include
│   ├── matrix.hh          # Matrix<T>
│   ├── matrix_base.hh     # CRTP base
│   ├── view.hh            # MatrixView<T>
│   ├── iterator.hh        # Strided iterator
│   ├── slice.hh           # Slicing operators
│   ├── binary_ops.hh      # Binary operations
│   ├── matmul_extension.hh # im2col, dot product
│   ├── op_functors.hh     # Operation functors
│   └── extra_functors.hh  # Neural network activations
└── device/
    ├── core.hh            # Device module include
    ├── matrix.hh          # device::Matrix<T>
    ├── matrix_base.hh     # CRTP base
    ├── view.hh            # device::MatrixView<T>
    ├── gpu_buffer.hh      # GPUBuffer<T>, GPUSpan<T> (RAII)
    ├── slice.hh           # Slicing operators
    ├── extra_functors.hh  # Neural network activations (__device__)
    ├── cuda_ops/          # CUDA implementations
    │   ├── gemm.hh        # cuBLAS GEMM wrapper
    │   ├── matmul.hh      # Batched matmul
    │   ├── reduce.hh      # Reduction kernels
    │   ├── scalar_op.hh   # Scalar operations
    │   └── init.hh        # Random initialization
    │   ├── group.hh       # Group reduce operations
    │   └── welford.hh     # Welford algorithm
    │   └── contiguity.hh  # Strided-to-contiguous copy
    └── kernels/           # CUDA kernel sources
        ├── elem_wise.cuh
        ├── group.cuh
        └── welford.cuh

Requirements

• C++20 compiler
• CUDA 12+ (optional, for GPU support)
• CMake 3.18+

License

MIT License - See LICENSE