Skip to content

Instantly share code, notes, and snippets.

@TruncatedDinoSour

TruncatedDinoSour/matrix.c

Last active August 6, 2024 17:05
Show Gist options
  • Save TruncatedDinoSour/5b2e939740a97e6942208a8b6a4dcc46 to your computer and use it in GitHub Desktop.
Save TruncatedDinoSour/5b2e939740a97e6942208a8b6a4dcc46 to your computer and use it in GitHub Desktop.
Download ZIP
Matrix arithmatic in C
Raw
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');
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment