Matrix arithmetic in C

matrix.c

#include <time.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

#include "include/mem.h"
#include "include/matrix.h"

#ifndef INFINITY
#define INFINITY HUGE_VAL
#endif /* INFINITY */

Bool matinit(Matrix *mat,
             const uint64_t rows,
             const uint64_t cols,
             const uint64_t stride) {
    if (!mat)
        return False;

    if (cols == 0 || rows == 0)
        return False;

    mat->data = Calloc(rows * cols, sizeof(*mat->data));

    if (!mat->data)
        return False;

    mat->rows   = rows;
    mat->cols   = cols;
    mat->stride = stride;

    return True;
}

Bool matinitrand(Matrix *mat,
                 const uint64_t rows,
                 const uint64_t cols,
                 const uint64_t stride,
                 const double init_scale) {
    if (!matinit(mat, rows, cols, stride))
        return False;

    if (!matrand(mat, init_scale)) {
        matdestroy(mat);
        return False;
    }

    return True;
}

Bool matresize(Matrix *mat,
               const uint64_t new_rows,
               const uint64_t new_cols,
               const uint64_t new_stride) {
    if (!mat)
        return False;

    if (new_cols == 0 || new_rows == 0)
        return False;

    double *new_data =
        Realloc(mat->data, new_rows * new_stride * sizeof(*new_data));

    if (!new_data)
        return False;

    mat->data   = new_data;
    mat->rows   = new_rows;
    mat->cols   = new_cols;
    mat->stride = new_stride;

    return True;
}

Bool matadd(const Matrix *A, const Matrix *B, Matrix *C) {
    uint64_t idx, jdx;

    if (!A || !B || !C)
        return False;

    if (A->rows != B->rows || A->cols != B->cols)
        return False;

    if (A->rows != C->rows || A->cols != C->cols)
        return False;

    for (idx = 0; idx < A->rows; ++idx)
        for (jdx = 0; jdx < A->cols; ++jdx)
            matset(C, idx, jdx, matget(A, idx, jdx) + matget(B, idx, jdx));

    return True;
}

Bool matsub(const Matrix *A, const Matrix *B, Matrix *C) {
    uint64_t idx, jdx;

    if (!A || !B || !C)
        return False;

    if (A->rows != B->rows || A->cols != B->cols)
        return False;

    if (A->rows != C->rows || A->cols != C->cols)
        return False;

    for (idx = 0; idx < A->rows; ++idx)
        for (jdx = 0; jdx < A->cols; ++jdx)
            matset(C, idx, jdx, matget(A, idx, jdx) - matget(B, idx, jdx));

    return True;
}

Bool matmul(const Matrix *A, const Matrix *B, Matrix *C) {
    double sum;
    uint64_t idx, jdx, kdx;

    if (!A || !B || !C)
        return False;

    if (A->cols != B->rows)
        return False;

    if (A->rows != C->rows || B->cols != C->cols)
        return False;

    for (idx = 0; idx < A->rows; ++idx) {
        for (jdx = 0; jdx < B->cols; ++jdx) {
            sum = 0;

            for (kdx = 0; kdx < A->cols; ++kdx)
                sum += matget(A, idx, kdx) * matget(B, kdx, jdx);

            matset(C, idx, jdx, sum);
        }
    }

    return True;
}

Bool matinv(const Matrix *A, Matrix *A_inv) {
    uint64_t idx, jdx, kdx;

    if (!A || !A_inv)
        return False;

    /* I'm so sorry for this code. Forgive me for my sins. */

    if (A->rows != A->cols)
        return False; /* Must be square */

    if (A->rows != A_inv->rows)
        return False; /* We can just check if one dimention matches */

    const uint64_t n = A->rows;

    Matrix aug;

    if (!matinit(&aug, n, 2 * n, 2 * n))
        return False;

    /* Form the augmented matrix [A | I] */

    for (idx = 0; idx < n; ++idx) {
        for (jdx = 0; jdx < n; ++jdx) {
            matset(&aug, idx, jdx, matget(A, idx, jdx));
            matset(&aug, idx, jdx + n, (idx == jdx) ? 1.0 : 0.0);
        }
    }

    /* Perform Gaussian elimination */

    for (idx = 0; idx < n; ++idx) {
        /* Partial pivoting */

        uint64_t max_row = idx;

        for (kdx = idx + 1; kdx < n; ++kdx) {
            if (fabs(matget(&aug, kdx, idx)) >
                fabs(matget(&aug, max_row, idx))) {
                max_row = kdx;
            }
        }

        if (max_row != idx)
            matswprow(&aug, idx, max_row);

        /* Scale the pivot row */

        const double pivot = matget(&aug, idx, idx);

        /* Likely singular */

        if (fabs(pivot) < 1e-14) {
            matdestroy(&aug);
            return False;
        }

        for (jdx = 0; jdx < 2 * n; ++jdx)
            matset(&aug, idx, jdx, matget(&aug, idx, jdx) / pivot);

        /* Eliminate below */

        for (kdx = 0; kdx < n; ++kdx) {
            if (kdx != idx) {
                const double factor = matget(&aug, kdx, idx);

                for (jdx = 0; jdx < 2 * n; ++jdx) {
                    matset(&aug, kdx, jdx,
                           matget(&aug, kdx, jdx) -
                               factor * matget(&aug, idx, jdx));
                }
            }
        }
    }

    /* Extract the inverse from the augmented matrix */

    for (idx = 0; idx < n; ++idx)
        for (jdx = 0; jdx < n; ++jdx)
            matset(A_inv, idx, jdx, matget(&aug, idx, jdx + n));

    matdestroy(&aug);

    return True;
}

Bool matdiv(const Matrix *A, const Matrix *B, Matrix *C) {
    Bool result;

    if (!A || !B || !C)
        return False;

    if (B->rows != B->cols || A->cols != B->rows)
        return False;

    Matrix B_inv;

    if (!matinit(&B_inv, B->rows, B->cols, B->stride))
        return False;

    if (!matinv(B, &B_inv)) {
        matdestroy(&B_inv);
        return False;
    }

    result = matmul(A, &B_inv, C);
    matdestroy(&B_inv);

    return result;
}

Bool matcpy(const Matrix *src, Matrix *dst) {
    uint64_t idx, jdx;

    if (!src || !dst)
        return False;

    if (src->rows != dst->rows || src->cols != dst->cols)
        return False;

    for (idx = 0; idx < src->rows; ++idx)
        for (jdx = 0; jdx < src->cols; ++jdx)
            matset(dst, idx, jdx, matget(src, idx, jdx));

    return True;
}

Bool matswprow(const Matrix *mat, const uint64_t r1, const uint64_t r2) {
    uint64_t col;

    if (!mat)
        return False;

    if (r1 == r2 || r1 > mat->rows || r2 > mat->rows)
        return False;

    for (col = 0; col < mat->cols; ++col) {
        const double temp = matget(mat, r1, col);
        matset(mat, r1, col, matget(mat, r2, col));
        matset(mat, r2, col, temp);
    }

    return True;
}

Bool mattanh(const Matrix *A, Matrix *H) {
    uint64_t idx, jdx;

    if (!A || !H)
        return False;

    if (A->rows != H->rows || A->cols != H->cols)
        return False;

    for (idx = 0; idx < A->rows; ++idx)
        for (jdx = 0; jdx < A->cols; ++jdx)
            matset(H, idx, jdx, tanh(matget(A, idx, jdx)));

    return True;
}

Bool matsoftmax(const Matrix *Z, Matrix *A) {
    uint64_t idx;
    double max_val = -INFINITY;

    if (!Z || !A)
        return False;

    if (Z->rows != A->rows || Z->cols != A->cols)
        return False;

    for (idx = 0; idx < Z->rows; ++idx) {
        const double val = matget(Z, idx, 0);

        if (val > max_val)
            max_val = val;
    }

    double sum_exp = 0.0;

    for (idx = 0; idx < Z->rows; ++idx) {
        const double exp_val = exp(matget(Z, idx, 0) - max_val);
        matset(A, idx, 0, exp_val);
        sum_exp += exp_val;
    }

    for (idx = 0; idx < A->rows; ++idx)
        matset(A, idx, 0, matget(A, idx, 0) / sum_exp);

    return True;
}

Bool matrelu(const Matrix *Z, Matrix *A) {
    uint64_t idx;

    if (!Z || !A)
        return False;

    if (Z->rows != A->rows || Z->cols != A->cols)
        return False;

    for (idx = 0; idx < Z->rows * Z->cols; ++idx)
        A->data[idx] = fmax(0.0, Z->data[idx]);

    return True;
}

Bool matrand(Matrix *mat, const double init_scale) {
    uint64_t idx, jdx;

    if (!mat)
        return False;

    for (idx = 0; idx < mat->rows; ++idx)
        for (jdx = 0; jdx < mat->cols; ++jdx)
            matset(mat, idx, jdx,
                   (double)rand() / (double)RAND_MAX * init_scale -
                       init_scale / (double)2.0f);

    return True;
}

Bool matclear(Matrix *mat) {
    if (!mat)
        return False;

    memset(mat->data, 0, mat->rows * mat->cols * sizeof(*mat->data));
    return True;
}

Bool matdestroy(Matrix *mat) {
    if (!mat)
        return False;

    Free(mat->data);
    return True;
}

void matprint(const Matrix *mat) {
    uint64_t idx, jdx;

    if (!mat)
        return;

    printf("Matrix [%zu x %zu]:\n", mat->rows, mat->cols);

    for (idx = 0; idx < mat->rows; ++idx) {
        for (jdx = 0; jdx < mat->cols; ++jdx)
            printf("%9g ", matget(mat, idx, jdx));

        putchar('\n');
    }
}