20171009:Practical:C:Matrix
Sommaire
Introduction
This week we will work on matrix. The purpose of this session is to introduce good practices on matrix manipulation in C.\
Submission
The directory in your GIT repository will be 20171009_Matrix, content will be specified in each exercise.
- S3: As for other session, your code is due to next monday (2017-10-16) 2pm.
- S3#: It is due on Friday, March 2, 2018 at 2pm
Forewords: Matrix in C, two dimensions in one
There exists many ways to represent matrix in C, the language provides static (i.e. fixed size) multidimensional arrays or let you play with dynamic arrays.
Experience shows that, in most cases, the best to implement matrix is to use a single dimension array and use offset inside the array. There are many reasons, first the code will be much more simpler, second we keep the matrix in a contiguous memory area, offering better performances in lot of cases.
Let's explain how it works: if you consider a simple 3×3 matrix, like the following one:
| 1 | 2 | 3 |
| 4 | 5 | 6 |
| 7 | 8 | 9 |
We will store it in row-major mode (see [1]) i.e. line by line. This means that we have a simple array of the form:
| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
Accessing a cell can then be done by recomputing the offset of the line (line number times number of column) and add the column number.
The following function illustrates this technique:
/* matrix_get(m, lines, cols, i, j) * returns the content of the cell (i,j) of lines x cols matrix m */ int matrix_get(int m[], size_t lines, size_t cols, size_t i, size_t j) { assert(i < lines); assert(j < cols); return m[j + i * cols]; }
Note: we use assert to enforce bounded access to matrix cells.
The important part to remember is (j + i * cols) the classical formula to transform the 2D index (i,j) (where i is the line number) into a single dimension index.
Standard operations on matrices
Submission sub-directory: MatrixStdOperations
In this part we will implement classical matrices operations: transpose, addition and multiplication. In addition we will implement a display function.
The code will be organized in the following files:
- MatrixStdOperations/
- Makefile
- matrix.h: header
- matrix.c: code
- matrix_test.c: some tests
All operations (but printing) will produce their results in a matrix provided as a parameter (so you don't have to care about allocating memory.)
Here is the header matrix.h:
/* * matrix.h: standard matrix operations * * 2017 Marwan Burelle / EPITA */ #ifndef _MATRIXSTDOPERATIONS_MATRIX_H_ #define _MATRIXSTDOPERATIONS_MATRIX_H_ /* * transpose(mat, lines, cols, res) * transpose the lines X cols matrix mat into a cols X lines matrix * result will be stored in res */ void transpose(double mat[], size_t lines, size_t cols, double res[]); /* * add(mat1, mat2, lines, cols, res) * computes res = mat1 + mat2 * All 3 matrices have the form lines X cols */ void add(double mat1[], double mat2[], size_t lines, size_t cols, double res[]); /* * mul(m1, m2, n, m, p, res) * computes res = m1 * m2 * Where: * m1 is a n X m matrix * m2 is a m X p matrix * res is a n X p matrix */ void mul(double m1[], double m2[], size_t n, size_t m, size_t p, double res[]); /* * print_matrix(mat, lines, cols) * prints the lines X cols matrix mat * All coeficient will be printed using the format " %4g ", line by line */ void print_matrix(double mat[], size_t lines, size_t cols); #endif /* _MATRIXSTDOPERATIONS_MATRIX_H_ */
Implementations
All the implementation code will be provided in the file matrix.c. You can assume that all dimensions are correct and matrices are correctly allocated.
Transpose
This operation just exchanges dimensions, see [2].
Addition
This is the matrix addition, the two matrix must have the same dimensions, see [3].
Multiplication
This is the matrix multiplication for the general case, the matrix m1 and m2 must verify that the number of columns of m1 match the number of lines of m2, see [4].
You can implement the straightforward algorithm for multiplication: the cell (i,j) in the result is the sum of m1(i, k) × m2(k, j) for all between 0 and m.
There exists more advanced algorithms if you want to go further, see [5].
Printing
Printing will be done in the most straightforward way, line by line using printf(3) with the format "%4g " (4 characters wide generic float conversion.)
Here is an example of output for the following matrix:
| 1 | 2 | 3 | 4 |
| 5 | 6 | 7 | 8 |
| 9 | 10 | 11 | 12 |
1 2 3 4 5 6 7 8 9 10 11 12
Tests
You can now implement some tests, we provide a skeleton for this but no operations (but printing) are tested in this file, you need to implement your own tests.
/* * matrix_test.c: standard matrix operations, very quick tests * * 2017 Marwan Burelle / EPITA */ #include <stdio.h> #include <stdlib.h> #include "matrix.h" /* * some matrices and their dimensions */ #define lines1 3 #define cols1 4 double m1[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0 }; #define lines2 4 #define cols2 3 double m2[] = { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0 }; #define lines3 3 #define cols3 3 double m3[] = { 1, 2, 3, 4, 5, 6, 7, 8, 9 }; double id[] = { 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 }; /* * some empty matrices with correct dimensions for results */ /* * lines1 X cols1 or cols1 X lines1 matrix */ double res1[ lines1 * cols1 ]; /* * lines1 X cols2 matrix (for m1 * m2 results, or m3 * m3 or m3 * id) */ double res2[ lines1 * cols2 ]; int main() { printf("m1 =\n"); print_matrix(m1, lines1, cols1); printf("\nm2 =\n"); print_matrix(m2, lines2, cols2); printf("\nm3 =\n"); print_matrix(m3, lines3, cols3); printf("\nid =\n"); print_matrix(id, lines3, cols3); /* TEST YOUR OPERATIONS HERE */ return 0; }
The following Makefile can be useful in order to compile your code:
# Makefile CC = gcc -fsanitize=address CPPFLAGS = -MMD CFLAGS = -Wall -Wextra -std=c99 -O0 -g LDFLAGS = LDLIBS = OBJ = matrix.o matrix_test.o DEP = ${OBJ:.o=.d} all: matrix_test matrix_test: ${OBJ} .PHONY: clean clean: ${RM} ${OBJ} ${DEP} matrix_test -include ${DEP} # END
Some Algebra
Submissions sub-directory: SomeAlgebra
The purpose of this exercise is to experiment with some special cases of matrix and vector operations that are useful in linear algebra. For vectors, we will interpret, depending on the context, its dimension as a 1×n or n×1 without loss of correctness.
The content of the exercise will be organised as follow:
- SomeAlgebra/
- Makefile
- algebra.h
- algebra.c
- algebra_test.c
In a lot of applications of linear algebra weed need to compute operations of the form (with M a matrix, v and b vectors):
- M * v: application of a linear transformation to a vector (in 3D this how you compute geometric transformation like rotation)
- f(M * v + b): this computes the activated output of a dense layer in a neural network with f the activation function (like sigmoid), here each line of M represents the weights of each neuron in the current layer, v is the output of the previous layer and each cell in b is the bias for the corresponding neuron.
Theoretically, M * v is just a multiplication between a m×n matrix and a n×1 matrix (yielding a m×1 matrix), but since we know that we work with vectors, we can simplify the code (and remove some transpositions in more complex cases).
Here is the header you have to implement:
/* * algebra.h: some useful linear algebra operations * * 2017 Marwan Burelle / EPITA */ #ifndef _SOMEALGEBRA_ALGEBRA_H_ #define _SOMEALGEBRA_ALGEBRA_H_ /* * func_t: function pointer type */ typedef double (*func_t)(double); /* * vector_apply(f, vect, size, res) * applies a function to all element on a vector of size n * for i in 0..size, res[i] = f[vect[i]] */ void vector_apply(func_t f, double vect[], size_t size, double res[]); /* * product(M, v, m, n, res) * computes res = M * v * where M is a m X n matrix, v is of length n and res of length m */ void product(double M[], double v[], size_t m, size_t n, double res[]); /* * vector_sum(v1, v2, size, res) * computes res = v1 + v2 * with res, v1 and v2 three vector of lenght size */ void vector_sum(double v1[], double v2[], size_t size, double res[]); #endif /* _SOMEALGEBRA_ALGEBRA_H_ */
Warning: in order to simplify the usage of the functions, you must assume that the result parameter can also be used as one of the parameter. For example, we must be able to do vector_apply(f, r, size, r) with r used as result and as operand. In most cases this is not a problem, but take of the way you write product.
Applying a function to a vector
The purpose of the function vector_apply is to apply a function to all elements in a vector.
Note on function pointer: a function pointer permits a bit of genericity, you can write a generic application function and provide to it any function (with the correct type). The only things that you need to know is that the parameter f can be use as normal function (i.e. you can call it like f(x)) in your function and you can call your function with the name of any function of the correct type (see examples in the tests part).
Product and sums
This two operations can be implemented using the generic matrix code of the previous exercise, but what is interesting here is to take benefit of the form of the vector.
Multiplication becomes: r[i] is the sum for all j of M(i,j) * v[j]. Sum is even more simpler: just sums the cells of the two vectors.
Testing
As for previous exercise, we provides a skeleton for tests, mainly to demonstrate how to use these functions. The tests are incomplete and you must implement your own. You must also adapt the previous Makefile for this exercise.
/* * algebra_test.c: some useful linear algebra operations, tests * * 2017 Marwan Burelle / EPITA */ #include <math.h> #include <stdio.h> #include <stdlib.h> #include "algebra.h" /* * sigmoid(x) the classical sigmoid function */ double sigmoid(double x) { return 1 / (1 + exp(-x)); } /* * 2D rotation of v by theta centered on the origin */ void rotation(double theta, double v[], double r[]) { double M[] = { cos(theta), -sin(theta), sin(theta), cos(theta) }; product(M, v, 2, 2, r); } void layer(double W[], double b[], double a[], size_t m, size_t n, double r[]) { product(W, a, m, n, r); vector_sum(b, r, m, r); vector_apply(sigmoid, r, m, r); } double W[] = { 6.36634445, 6.36769295 }; double b[] = { -9.00631428 }; void test(void) { double a[] = {1,1}; double r[] = {0}; layer(W, b, a, 1, 2, r); printf("1, 1 -> %g\n", r[0]); a[0] = 0; a[1] = 0; layer(W, b, a, 1, 2, r); printf("0, 0 -> %g\n", r[0]); a[0] = 1; a[1] = 0; layer(W, b, a, 1, 2, r); printf("1, 0 -> %g\n", r[0]); a[0] = 0; a[1] = 1; layer(W, b, a, 1, 2, r); printf("0, 1 -> %g\n", r[0]); } int main() { test(); return 0; }
