diff --git a/toys/mat.c b/toys/mat.c
new file mode 100644
index 0000000..bb729f5
--- /dev/null
+++ b/toys/mat.c
@@ -0,0 +1,124 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+int rowcmp(float *row0, float *row1, int n)
+{
+	for (int i = 0; i < n; ++i) {
+		if (row0[i])
+			return 0;
+		if (row1[i])
+			return 1;
+	}
+	return 0;
+}
+
+/* Upper-triangularise matrix mat of size mxn */
+void ut(float *mat, int m, int n)
+{
+	long i, j, k;
+	float *current, *other, *tmp = malloc(n * sizeof(float));
+
+	for (i = 0; i + 1 < m; ++i) {
+		current = mat + i * n;
+		for (j = i + 1; j < m; ++j) {
+			other = mat + j * n;
+			if (rowcmp(current, other, n)) {
+				memcpy(tmp, current, n);
+				memcpy(current, other, n);
+				memcpy(other, tmp, n);
+			}
+		}
+
+		for (j = i + 1; j < m; ++j) {
+			for (k = 0; k < n && !mat[i*n + k]; ++k);
+			if (k == n)
+				break;
+
+			float c = mat[j*n + k] / mat[i*n + k];
+			mat[j*n + k] = 0.0; /* floating point imprecision */
+			while (++k < n)
+				mat[j*n + k] -= c * mat[i*n + k];
+		}
+	}
+	free(tmp);
+}
+
+float *transpose(float *mat, int m, int n)
+{
+	float *result = malloc(m * n * sizeof(float));
+	for (long i = 0; i < m; ++i)
+		for (long j = 0; j < n; ++j)
+			result[j*m + i] = mat[i*n + j];
+	return result;
+}
+
+void printmatf(char *format, float *mat, int m, int n)
+{
+	for (long i = 0; i < m; ++i) {
+		for (long j = 0; j < n; ++j)
+			printf(format, mat[i*n + j]);
+		putchar(10);
+	}
+}
+
+int main()
+{
+	int n;
+	puts("Please input the size of the square matrix:");
+	fputs("n = ", stdout);
+	scanf("%d", &n);
+
+	long i, j, k, len = n * n * sizeof(float);
+	float sum = 0.0, *a = malloc(len);
+	puts("Now please input the matrix:");
+	for (i = 0; i < n * n; ++i) {
+		scanf("%f", a + i);
+		sum += a[i];
+	}
+	printf("The sum of all the matrix elements: %.2f\n", sum);
+
+	float *b = transpose(a, n, n);
+	puts("The transpose of the matrix:");
+	printmatf("%8.3f", b, n, n);
+	free(a);
+
+	b = realloc(b, len * 2);
+	for (i = 0; i < n; ++i)
+		for (j = 0; j < n; ++j)
+			b[n*n + i*n + j] = i == j;
+	a = transpose(b, n * 2, n);
+	ut(a, n, n * 2);
+	free(b);
+
+	float det = 1.0;
+	for (i = 0; i < n; ++i)
+		det *= a[i*n*2 + i];
+	printf("The determinant of the matrix: %.2f\n", det);
+
+	if (!det) {
+		puts("The given matrix is not invertible.");
+		free(a);
+		return 0;
+	}
+
+	for (i = n - 1; i >= 0; --i) {
+		float d = a[i*n*2 + i];
+		for (j = n; j < n * 2; ++j)
+			a[i*n*2 + j] /= d;
+		for (k = i - 1; k >= 0; --k) {
+			float c = a[k*n*2 + i];
+			for (j = n; j < n * 2; ++j)
+				a[k*n*2 + j] -= a[i*n*2 + j] * c;
+		}
+	}
+
+	b = transpose(a, n, n * 2);
+	a = transpose(b + n * n, n, n);
+	puts("The inverse of the matrix:");
+	printmatf("%8.3f", a, n, n);
+	free(a);
+	free(b);
+
+	return 0;
+}