diff --git a/clients/calibrator.c b/clients/calibrator.c
index be65d9bdc..b98a676b3 100644
--- a/clients/calibrator.c
+++ b/clients/calibrator.c
@@ -117,11 +117,11 @@ finish_calibration (struct calibrator *calibrator)
 	 */
 	memset(&m, 0, sizeof(m));
 	for (i = 0; i < (int)ARRAY_LENGTH(test_ratios); i++) {
-		m.d[i] = calibrator->tests[i].clicked_x;
-		m.d[i + 4] = calibrator->tests[i].clicked_y;
-		m.d[i + 8] = 1;
+		m.M.col[0].el[i] = calibrator->tests[i].clicked_x;
+		m.M.col[1].el[i] = calibrator->tests[i].clicked_y;
+		m.M.col[2].el[i] = 1;
 	}
-	m.d[15] = 1;
+	m.M.col[3].el[3] = 1;
 
 	weston_matrix_invert(&inverse, &m);
 
@@ -129,8 +129,8 @@ finish_calibration (struct calibrator *calibrator)
 	memset(&y_calib, 0, sizeof(y_calib));
 
 	for (i = 0; i < (int)ARRAY_LENGTH(test_ratios); i++) {
-		x_calib.f[i] = calibrator->tests[i].drawn_x;
-		y_calib.f[i] = calibrator->tests[i].drawn_y;
+		x_calib.v.el[i] = calibrator->tests[i].drawn_x;
+		y_calib.v.el[i] = calibrator->tests[i].drawn_y;
 	}
 
 	/* Multiples into the vector */
@@ -138,8 +138,8 @@ finish_calibration (struct calibrator *calibrator)
 	weston_matrix_transform(&inverse, &y_calib);
 
 	printf ("Calibration values: %f %f %f %f %f %f\n",
-		x_calib.f[0], x_calib.f[1], x_calib.f[2],
-		y_calib.f[0], y_calib.f[1], y_calib.f[2]);
+		x_calib.v.el[0], x_calib.v.el[1], x_calib.v.el[2],
+		y_calib.v.el[0], y_calib.v.el[1], y_calib.v.el[2]);
 
 	exit(0);
 }
diff --git a/clients/simple-dmabuf-egl.c b/clients/simple-dmabuf-egl.c
index 9989497e7..92c63ac87 100644
--- a/clients/simple-dmabuf-egl.c
+++ b/clients/simple-dmabuf-egl.c
@@ -815,7 +815,7 @@ render(struct window *window, struct buffer *buffer)
 
 	glUniform1f(window->gl.offset_uniform, offset);
 	glUniformMatrix4fv(window->gl.reflection_uniform, 1, GL_FALSE,
-			   (GLfloat *) reflection.d);
+			   (GLfloat *) reflection.M.colmaj);
 
 	glClearColor(0.0,0.0, 0.0, 1.0);
 	glClear(GL_COLOR_BUFFER_BIT);
@@ -864,7 +864,7 @@ render_mandelbrot(struct window *window, struct buffer *buffer)
 
 	glUniform1f(window->gl.offset_uniform, offset);
 	glUniformMatrix4fv(window->gl.reflection_uniform, 1, GL_FALSE,
-			   (GLfloat *) reflection.d);
+			   (GLfloat *) reflection.M.colmaj);
 
 	glClearColor(0.6, 0.6, 0.6, 1.0);
 	glClear(GL_COLOR_BUFFER_BIT);
diff --git a/clients/simple-egl.c b/clients/simple-egl.c
index fadeac34c..8308fe6a5 100644
--- a/clients/simple-egl.c
+++ b/clients/simple-egl.c
@@ -940,10 +940,10 @@ redraw(struct window *window)
 		angle = ((time - window->initial_frame_time) / speed_div)
 			% 360 * M_PI / 180.0;
 	}
-	rotation.d[0] =   cos(angle);
-	rotation.d[2] =   sin(angle);
-	rotation.d[8] =  -sin(angle);
-	rotation.d[10] =  cos(angle);
+	rotation.M.col[0].el[0] =  cos(angle);
+	rotation.M.col[0].el[2] =  sin(angle);
+	rotation.M.col[2].el[0] = -sin(angle);
+	rotation.M.col[2].el[2] =  cos(angle);
 
 	switch (window->buffer_transform) {
 	case WL_OUTPUT_TRANSFORM_FLIPPED:
@@ -982,7 +982,7 @@ redraw(struct window *window)
 	glViewport(0, 0, window->buffer_size.width, window->buffer_size.height);
 
 	glUniformMatrix4fv(window->gl.rotation_uniform, 1, GL_FALSE,
-			   (GLfloat *) rotation.d);
+			   (GLfloat *) rotation.M.colmaj);
 
 	if (window->opaque || window->fullscreen)
 		glClearColor(0.0, 0.0, 0.0, 1);
diff --git a/clients/touch-calibrator.c b/clients/touch-calibrator.c
index 49dd9206b..374ae492e 100644
--- a/clients/touch-calibrator.c
+++ b/clients/touch-calibrator.c
@@ -330,9 +330,9 @@ compute_calibration(struct calibrator *cal, float *result)
 	 */
 	weston_matrix_init(&m);
 	for (i = 0; i < 3; i++) {
-		m.d[i + 0] = cal->samples[i].touched.x;
-		m.d[i + 4] = cal->samples[i].touched.y;
-		m.d[i + 8] = 1.0f;
+		m.M.col[0].el[i] = cal->samples[i].touched.x;
+		m.M.col[1].el[i] = cal->samples[i].touched.y;
+		m.M.col[2].el[i] = 1.0f;
 	}
 	m.type = WESTON_MATRIX_TRANSFORM_OTHER;
 
@@ -342,20 +342,20 @@ compute_calibration(struct calibrator *cal, float *result)
 	}
 
 	for (i = 0; i < 3; i++) {
-		x_calib.f[i] = cal->samples[i].drawn_cal.x;
-		y_calib.f[i] = cal->samples[i].drawn_cal.y;
+		x_calib.v.el[i] = cal->samples[i].drawn_cal.x;
+		y_calib.v.el[i] = cal->samples[i].drawn_cal.y;
 	}
-	x_calib.f[3] = 0.0f;
-	y_calib.f[3] = 0.0f;
+	x_calib.v.el[3] = 0.0f;
+	y_calib.v.el[3] = 0.0f;
 
 	/* Multiples into the vector */
 	weston_matrix_transform(&inverse, &x_calib);
 	weston_matrix_transform(&inverse, &y_calib);
 
 	for (i = 0; i < 3; i++)
-		result[i] = x_calib.f[i];
+		result[i] = x_calib.v.el[i];
 	for (i = 0; i < 3; i++)
-		result[i + 3] = y_calib.f[i];
+		result[i + 3] = y_calib.v.el[i];
 
 	return 0;
 }
diff --git a/clients/window.c b/clients/window.c
index 9d0aaf964..104bb3991 100644
--- a/clients/window.c
+++ b/clients/window.c
@@ -1842,8 +1842,9 @@ widget_cairo_update_transform(struct widget *widget, cairo_t *cr)
 				     surface->allocation.width,
 				     surface->allocation.height,
 				     surface->buffer_scale);
-	cairo_matrix_init(&m, matrix.d[0], matrix.d[1], matrix.d[4],
-			  matrix.d[5], matrix.d[12], matrix.d[13]);
+	cairo_matrix_init(&m, matrix.M.col[0].x, matrix.M.col[0].y,
+			      matrix.M.col[1].x, matrix.M.col[1].y,
+			      matrix.M.col[3].x, matrix.M.col[3].y);
 	cairo_transform(cr, &m);
 }
 
diff --git a/include/libweston/matrix.h b/include/libweston/matrix.h
index c309db48d..ce4af304b 100644
--- a/include/libweston/matrix.h
+++ b/include/libweston/matrix.h
@@ -31,6 +31,7 @@
 #include <stdbool.h>
 
 #include <wayland-server-protocol.h>
+#include <libweston/linalg-types.h>
 
 #ifdef  __cplusplus
 extern "C" {
@@ -44,12 +45,12 @@ enum weston_matrix_transform_type {
 };
 
 struct weston_matrix {
-	float d[16];
+	struct weston_mat4f M;
 	unsigned int type;
 };
 
 struct weston_vector {
-	float f[4];
+	struct weston_vec4f v;
 };
 
 /** Arbitrary coordinates in any space */
diff --git a/ivi-shell/ivi-layout.c b/ivi-shell/ivi-layout.c
index a50ed5a9e..20388ea55 100644
--- a/ivi-shell/ivi-layout.c
+++ b/ivi-shell/ivi-layout.c
@@ -463,33 +463,33 @@ calc_inverse_matrix_transform(const struct weston_matrix *matrix,
 	}
 
 	/* The vectors and matrices involved will always produce f[3] == 1.0. */
-	top_left.f[0] = rect_input->x;
-	top_left.f[1] = rect_input->y;
-	top_left.f[2] = 0.0f;
-	top_left.f[3] = 1.0f;
+	top_left.v.el[0] = rect_input->x;
+	top_left.v.el[1] = rect_input->y;
+	top_left.v.el[2] = 0.0f;
+	top_left.v.el[3] = 1.0f;
 
-	bottom_right.f[0] = rect_input->x + rect_input->width;
-	bottom_right.f[1] = rect_input->y + rect_input->height;
-	bottom_right.f[2] = 0.0f;
-	bottom_right.f[3] = 1.0f;
+	bottom_right.v.el[0] = rect_input->x + rect_input->width;
+	bottom_right.v.el[1] = rect_input->y + rect_input->height;
+	bottom_right.v.el[2] = 0.0f;
+	bottom_right.v.el[3] = 1.0f;
 
 	weston_matrix_transform(&m, &top_left);
 	weston_matrix_transform(&m, &bottom_right);
 
-	if (top_left.f[0] < bottom_right.f[0]) {
-		rect_output->x = floorf(top_left.f[0]);
-		rect_output->width = ceilf(bottom_right.f[0] - rect_output->x);
+	if (top_left.v.el[0] < bottom_right.v.el[0]) {
+		rect_output->x = floor(top_left.v.el[0]);
+		rect_output->width = ceil(bottom_right.v.el[0] - rect_output->x);
 	} else {
-		rect_output->x = floorf(bottom_right.f[0]);
-		rect_output->width = ceilf(top_left.f[0] - rect_output->x);
+		rect_output->x = floor(bottom_right.v.el[0]);
+		rect_output->width = ceil(top_left.v.el[0] - rect_output->x);
 	}
 
-	if (top_left.f[1] < bottom_right.f[1]) {
-		rect_output->y = floorf(top_left.f[1]);
-		rect_output->height = ceilf(bottom_right.f[1] - rect_output->y);
+	if (top_left.v.el[1] < bottom_right.v.el[1]) {
+		rect_output->y = floor(top_left.v.el[1]);
+		rect_output->height = ceil(bottom_right.v.el[1] - rect_output->y);
 	} else {
-		rect_output->y = floorf(bottom_right.f[1]);
-		rect_output->height = ceilf(top_left.f[1] - rect_output->y);
+		rect_output->y = floor(bottom_right.v.el[1]);
+		rect_output->height = ceil(top_left.v.el[1] - rect_output->y);
 	}
 
 	ivi_rectangle_intersect(rect_output, boundingbox, rect_output);
diff --git a/libweston/compositor.c b/libweston/compositor.c
index 819fb9530..2fcbf6bf0 100644
--- a/libweston/compositor.c
+++ b/libweston/compositor.c
@@ -1798,14 +1798,14 @@ weston_view_update_transform_disable(struct weston_view *view)
 
 	/* Otherwise identity matrix, but with x and y translation. */
 	view->transform.position.matrix.type = WESTON_MATRIX_TRANSFORM_TRANSLATE;
-	view->transform.position.matrix.d[12] = view->geometry.pos_offset.x;
-	view->transform.position.matrix.d[13] = view->geometry.pos_offset.y;
+	view->transform.position.matrix.M.col[3].x = view->geometry.pos_offset.x;
+	view->transform.position.matrix.M.col[3].y = view->geometry.pos_offset.y;
 
 	view->transform.matrix = view->transform.position.matrix;
 
 	view->transform.inverse = view->transform.position.matrix;
-	view->transform.inverse.d[12] = -view->geometry.pos_offset.x;
-	view->transform.inverse.d[13] = -view->geometry.pos_offset.y;
+	view->transform.inverse.M.col[3].x = -view->geometry.pos_offset.x;
+	view->transform.inverse.M.col[3].y = -view->geometry.pos_offset.y;
 
 	pixman_region32_init_rect(&view->transform.boundingbox,
 				  0, 0,
@@ -1850,8 +1850,8 @@ weston_view_update_transform_enable(struct weston_view *view)
 
 	/* Otherwise identity matrix, but with x and y translation. */
 	view->transform.position.matrix.type = WESTON_MATRIX_TRANSFORM_TRANSLATE;
-	view->transform.position.matrix.d[12] = view->geometry.pos_offset.x;
-	view->transform.position.matrix.d[13] = view->geometry.pos_offset.y;
+	view->transform.position.matrix.M.col[3].x = view->geometry.pos_offset.x;
+	view->transform.position.matrix.M.col[3].y = view->geometry.pos_offset.y;
 
 	weston_matrix_init(matrix);
 	wl_list_for_each(tform, &view->geometry.transformation_list, link)
@@ -1889,8 +1889,8 @@ weston_view_update_transform_enable(struct weston_view *view)
 							  &view->transform.opaque,
 							  &view->geometry.scissor);
 			pixman_region32_translate(&view->transform.opaque,
-						  matrix->d[12],
-						  matrix->d[13]);
+						  matrix->M.col[3].x,
+						  matrix->M.col[3].y);
 		}
 	} else if (view->alpha == 1.0 &&
 		 matrix->type < WESTON_MATRIX_TRANSFORM_ROTATE &&
@@ -10751,7 +10751,7 @@ weston_output_finish_frame_from_timer(struct weston_output *output)
 }
 
 /** Retrieve the backend type of as described in enum
- * weston_compositor_backend. 
+ * weston_compositor_backend.
  *
  * Note that the backend must be loaded, with weston_compositor_load_backend
  *
diff --git a/libweston/pixman-renderer.c b/libweston/pixman-renderer.c
index ec7ddf5e6..84031b703 100644
--- a/libweston/pixman-renderer.c
+++ b/libweston/pixman-renderer.c
@@ -152,15 +152,15 @@ weston_matrix_to_pixman_transform(pixman_transform_t *pt,
 {
 	/* Pixman supports only 2D transform matrix, but Weston uses 3D, *
 	 * so we're omitting Z coordinate here. */
-	pt->matrix[0][0] = pixman_double_to_fixed(wm->d[0]);
-	pt->matrix[0][1] = pixman_double_to_fixed(wm->d[4]);
-	pt->matrix[0][2] = pixman_double_to_fixed(wm->d[12]);
-	pt->matrix[1][0] = pixman_double_to_fixed(wm->d[1]);
-	pt->matrix[1][1] = pixman_double_to_fixed(wm->d[5]);
-	pt->matrix[1][2] = pixman_double_to_fixed(wm->d[13]);
-	pt->matrix[2][0] = pixman_double_to_fixed(wm->d[3]);
-	pt->matrix[2][1] = pixman_double_to_fixed(wm->d[7]);
-	pt->matrix[2][2] = pixman_double_to_fixed(wm->d[15]);
+	pt->matrix[0][0] = pixman_double_to_fixed(wm->M.col[0].x);
+	pt->matrix[0][1] = pixman_double_to_fixed(wm->M.col[1].x);
+	pt->matrix[0][2] = pixman_double_to_fixed(wm->M.col[3].x);
+	pt->matrix[1][0] = pixman_double_to_fixed(wm->M.col[0].y);
+	pt->matrix[1][1] = pixman_double_to_fixed(wm->M.col[1].y);
+	pt->matrix[1][2] = pixman_double_to_fixed(wm->M.col[3].y);
+	pt->matrix[2][0] = pixman_double_to_fixed(wm->M.col[0].w);
+	pt->matrix[2][1] = pixman_double_to_fixed(wm->M.col[1].w);
+	pt->matrix[2][2] = pixman_double_to_fixed(wm->M.col[3].w);
 }
 
 static bool
diff --git a/libweston/renderer-gl/gl-renderer.c b/libweston/renderer-gl/gl-renderer.c
index 96155b98f..d19ca5d04 100644
--- a/libweston/renderer-gl/gl-renderer.c
+++ b/libweston/renderer-gl/gl-renderer.c
@@ -43,6 +43,8 @@
 #include <gbm.h>
 #endif
 
+#include <libweston/linalg-4.h>
+
 #include "linux-sync-file.h"
 #include "timeline.h"
 
@@ -2343,12 +2345,12 @@ blit_shadow_to_output(struct weston_output *output,
 			.input_is_premult = true,
 		},
 		.projection = {
-			.d = { /* transpose */
-				 2.0f,	0.0f,              0.0f, 0.0f,
-				 0.0f,  go->y_flip * 2.0f, 0.0f, 0.0f,
-				 0.0f,  0.0f,              1.0f, 0.0f,
-				-1.0f, -go->y_flip,        0.0f, 1.0f
-			},
+			.M = WESTON_MAT4F(
+				2.0,              0.0, 0.0,        -1.0,
+				0.0, go->y_flip * 2.0, 0.0, -go->y_flip,
+				0.0,              0.0, 1.0,         0.0,
+				0.0,              0.0, 0.0,         1.0
+			),
 			.type = WESTON_MATRIX_TRANSFORM_SCALE |
 				WESTON_MATRIX_TRANSFORM_TRANSLATE,
 		},
@@ -3926,9 +3928,9 @@ gl_renderer_surface_copy_content(struct weston_surface *surface,
 	glViewport(0, 0, cw, ch);
 	set_blend_state(gr, false);
 	if (buffer->buffer_origin == ORIGIN_TOP_LEFT)
-		ARRAY_COPY(sconf.projection.d, projmat_normal);
+		ARRAY_COPY(sconf.projection.M.colmaj, projmat_normal);
 	else
-		ARRAY_COPY(sconf.projection.d, projmat_yinvert);
+		ARRAY_COPY(sconf.projection.M.colmaj, projmat_yinvert);
 	sconf.projection.type = WESTON_MATRIX_TRANSFORM_SCALE |
 				WESTON_MATRIX_TRANSFORM_TRANSLATE;
 
diff --git a/libweston/renderer-gl/gl-shaders.c b/libweston/renderer-gl/gl-shaders.c
index eeaad86e2..5332d1fb3 100644
--- a/libweston/renderer-gl/gl-shaders.c
+++ b/libweston/renderer-gl/gl-shaders.c
@@ -663,11 +663,11 @@ gl_shader_load_config(struct gl_renderer *gr,
 	int i, j;
 
 	glUniformMatrix4fv(shader->proj_uniform,
-			   1, GL_FALSE, sconf->projection.d);
+			   1, GL_FALSE, sconf->projection.M.colmaj);
 
 	if (shader->surface_to_buffer_uniform != -1)
 		glUniformMatrix4fv(shader->surface_to_buffer_uniform,
-			           1, GL_FALSE, sconf->surface_to_buffer.d);
+			           1, GL_FALSE, sconf->surface_to_buffer.M.colmaj);
 
 	if (shader->color_uniform != -1)
 		glUniform4fv(shader->color_uniform, 1, sconf->unicolor);
diff --git a/shared/matrix.c b/shared/matrix.c
index b78ff38ae..9d1819cbb 100644
--- a/shared/matrix.c
+++ b/shared/matrix.c
@@ -48,66 +48,37 @@
 WL_EXPORT void
 weston_matrix_init(struct weston_matrix *matrix)
 {
-	static const struct weston_matrix identity = {
-		.d = { 1, 0, 0, 0,  0, 1, 0, 0,  0, 0, 1, 0,  0, 0, 0, 1 },
-		.type = 0,
-	};
-
-	memcpy(matrix, &identity, sizeof identity);
+	matrix->M = WESTON_MAT4F_IDENTITY;
+	matrix->type = 0;
 }
 
 /* m <- n * m, that is, m is multiplied on the LEFT. */
 WL_EXPORT void
 weston_matrix_multiply(struct weston_matrix *m, const struct weston_matrix *n)
 {
-	struct weston_matrix tmp;
-	const float *row, *column;
-	int i, j, k;
-
-	for (i = 0; i < 4; i++) {
-		row = m->d + i * 4;
-		for (j = 0; j < 4; j++) {
-			tmp.d[4 * i + j] = 0;
-			column = n->d + j;
-			for (k = 0; k < 4; k++)
-				tmp.d[4 * i + j] += row[k] * column[k * 4];
-		}
-	}
-	tmp.type = m->type | n->type;
-	memcpy(m, &tmp, sizeof tmp);
+	m->M = weston_m4f_mul_m4f(n->M, m->M);
+	m->type |= n->type;
 }
 
 WL_EXPORT void
 weston_matrix_translate(struct weston_matrix *matrix, float x, float y, float z)
 {
-	struct weston_matrix translate = {
-		.d = { 1, 0, 0, 0,  0, 1, 0, 0,  0, 0, 1, 0,  x, y, z, 1 },
-		.type = WESTON_MATRIX_TRANSFORM_TRANSLATE,
-	};
-
-	weston_matrix_multiply(matrix, &translate);
+	matrix->M = weston_m4f_mul_m4f(weston_m4f_translation(x, y, z), matrix->M);
+	matrix->type = WESTON_MATRIX_TRANSFORM_TRANSLATE;
 }
 
 WL_EXPORT void
 weston_matrix_scale(struct weston_matrix *matrix, float x, float y,float z)
 {
-	struct weston_matrix scale = {
-		.d = { x, 0, 0, 0,  0, y, 0, 0,  0, 0, z, 0,  0, 0, 0, 1 },
-		.type = WESTON_MATRIX_TRANSFORM_SCALE,
-	};
-
-	weston_matrix_multiply(matrix, &scale);
+	matrix->M = weston_m4f_mul_m4f(weston_m4f_scaling(x, y, z), matrix->M);
+	matrix->type = WESTON_MATRIX_TRANSFORM_SCALE;
 }
 
 WL_EXPORT void
 weston_matrix_rotate_xy(struct weston_matrix *matrix, float cos, float sin)
 {
-	struct weston_matrix translate = {
-		.d = { cos, sin, 0, 0,  -sin, cos, 0, 0,  0, 0, 1, 0,  0, 0, 0, 1 },
-		.type = WESTON_MATRIX_TRANSFORM_ROTATE,
-	};
-
-	weston_matrix_multiply(matrix, &translate);
+	matrix->M = weston_m4f_mul_m4f(weston_m4f_rotation_xy(cos, sin), matrix->M);
+	matrix->type = WESTON_MATRIX_TRANSFORM_ROTATE;
 }
 
 /* v <- m * v */
@@ -115,16 +86,7 @@ WL_EXPORT void
 weston_matrix_transform(const struct weston_matrix *matrix,
 			struct weston_vector *v)
 {
-	int i, j;
-	struct weston_vector t;
-
-	for (i = 0; i < 4; i++) {
-		t.f[i] = 0;
-		for (j = 0; j < 4; j++)
-			t.f[i] += v->f[j] * matrix->d[i + j * 4];
-	}
-
-	*v = t;
+	v->v = weston_m4f_mul_v4f(matrix->M, v->v);
 }
 
 WL_EXPORT struct weston_coord
@@ -132,17 +94,29 @@ weston_matrix_transform_coord(const struct weston_matrix *matrix,
 			      struct weston_coord c)
 {
 	struct weston_coord out;
-	struct weston_vector t = { { c.x, c.y, 0.0, 1.0 } };
+	struct weston_vector t = { .v.el = { c.x, c.y, 0.0, 1.0 } };
 
 	weston_matrix_transform(matrix, &t);
 
-	assert(fabsf(t.f[3]) > 1e-6);
+	assert(fabsf(t.v.el[3]) > 1e-6);
 
-	out.x = t.f[0] / t.f[3];
-	out.y = t.f[1] / t.f[3];
+	out.x = t.v.el[0] / t.v.el[3];
+	out.y = t.v.el[1] / t.v.el[3];
 	return out;
 }
 
+WL_EXPORT int
+weston_matrix_invert(struct weston_matrix *inverse,
+		     const struct weston_matrix *matrix)
+{
+	if (weston_m4f_invert(&inverse->M, matrix->M)) {
+		inverse->type = matrix->type;
+		return 0;
+	}
+
+	return -1;
+}
+
 static inline void
 swap_rows(double *a, double *b)
 {
@@ -177,118 +151,6 @@ find_pivot(double *column, unsigned k)
 	return p;
 }
 
-/*
- * reference: Gene H. Golub and Charles F. van Loan. Matrix computations.
- * 3rd ed. The Johns Hopkins University Press. 1996.
- * LU decomposition, forward and back substitution: Chapter 3.
- */
-
-static int
-matrix_invert(double *A, unsigned *p, const struct weston_matrix *matrix)
-{
-	unsigned i, j, k;
-	unsigned pivot;
-	double pv;
-
-	for (i = 0; i < 4; ++i)
-		p[i] = i;
-	for (i = 16; i--; )
-		A[i] = matrix->d[i];
-
-	/* LU decomposition with partial pivoting */
-	for (k = 0; k < 4; ++k) {
-		pivot = find_pivot(&A[k * 4], k);
-		if (pivot != k) {
-			swap_unsigned(&p[k], &p[pivot]);
-			swap_rows(&A[k], &A[pivot]);
-		}
-
-		pv = A[k * 4 + k];
-		if (fabs(pv) < 1e-9)
-			return -1; /* zero pivot, not invertible */
-
-		for (i = k + 1; i < 4; ++i) {
-			A[i + k * 4] /= pv;
-
-			for (j = k + 1; j < 4; ++j)
-				A[i + j * 4] -= A[i + k * 4] * A[k + j * 4];
-		}
-	}
-
-	return 0;
-}
-
-static void
-inverse_transform(const double *LU, const unsigned *p, float *v)
-{
-	/* Solve A * x = v, when we have P * A = L * U.
-	 * P * A * x = P * v  =>  L * U * x = P * v
-	 * Let U * x = b, then L * b = P * v.
-	 */
-	double b[4];
-	unsigned j;
-
-	/* Forward substitution, column version, solves L * b = P * v */
-	/* The diagonal of L is all ones, and not explicitly stored. */
-	b[0] = v[p[0]];
-	b[1] = (double)v[p[1]] - b[0] * LU[1 + 0 * 4];
-	b[2] = (double)v[p[2]] - b[0] * LU[2 + 0 * 4];
-	b[3] = (double)v[p[3]] - b[0] * LU[3 + 0 * 4];
-	b[2] -= b[1] * LU[2 + 1 * 4];
-	b[3] -= b[1] * LU[3 + 1 * 4];
-	b[3] -= b[2] * LU[3 + 2 * 4];
-
-	/* backward substitution, column version, solves U * y = b */
-#if 1
-	/* hand-unrolled, 25% faster for whole function */
-	b[3] /= LU[3 + 3 * 4];
-	b[0] -= b[3] * LU[0 + 3 * 4];
-	b[1] -= b[3] * LU[1 + 3 * 4];
-	b[2] -= b[3] * LU[2 + 3 * 4];
-
-	b[2] /= LU[2 + 2 * 4];
-	b[0] -= b[2] * LU[0 + 2 * 4];
-	b[1] -= b[2] * LU[1 + 2 * 4];
-
-	b[1] /= LU[1 + 1 * 4];
-	b[0] -= b[1] * LU[0 + 1 * 4];
-
-	b[0] /= LU[0 + 0 * 4];
-#else
-	for (j = 3; j > 0; --j) {
-		unsigned k;
-		b[j] /= LU[j + j * 4];
-		for (k = 0; k < j; ++k)
-			b[k] -= b[j] * LU[k + j * 4];
-	}
-
-	b[0] /= LU[0 + 0 * 4];
-#endif
-
-	/* the result */
-	for (j = 0; j < 4; ++j)
-		v[j] = b[j];
-}
-
-WL_EXPORT int
-weston_matrix_invert(struct weston_matrix *inverse,
-		     const struct weston_matrix *matrix)
-{
-	double LU[16];		/* column-major */
-	unsigned perm[4];	/* permutation */
-	unsigned c;
-
-	if (matrix_invert(LU, perm, matrix) < 0)
-		return -1;
-
-	weston_matrix_init(inverse);
-	for (c = 0; c < 4; ++c)
-		inverse_transform(LU, perm, &inverse->d[c * 4]);
-	inverse->type = matrix->type;
-
-	return 0;
-}
-
 static bool
 m4f_LU_decompose(double *restrict LU, unsigned *restrict p, struct weston_mat4f M)
 {
@@ -521,7 +383,7 @@ get_el(const struct weston_matrix *matrix, int row, int col)
 	assert(row >= 0 && row <= 3);
 	assert(col >= 0 && col <= 3);
 
-	return matrix->d[col * 4 + row];
+	return matrix->M.col[col].el[row];
 }
 
 static bool
diff --git a/tests/color_util.c b/tests/color_util.c
index d18800935..413da9860 100644
--- a/tests/color_util.c
+++ b/tests/color_util.c
@@ -366,7 +366,7 @@ weston_matrix_from_lcmsMAT3(struct weston_matrix *w, const struct lcmsMAT3 *m)
 
 	for (c = 0; c < 3; c++) {
 		for (r = 0; r < 3; r++)
-			w->d[c * 4 + r] = m->v[c].n[r];
+			w->M.col[c].el[r] = m->v[c].n[r];
 	}
 }
 
@@ -377,7 +377,7 @@ lcmsMAT3_from_weston_matrix(struct lcmsMAT3 *m, const struct weston_matrix *w)
 
 	for (c = 0; c < 3; c++) {
 		for (r = 0; r < 3; r++)
-			m->v[c].n[r] = w->d[c * 4 + r];
+			m->v[c].n[r] = w->M.col[c].el[r];
 	}
 }
 
diff --git a/tests/lcms_util.c b/tests/lcms_util.c
index 3fe606329..5f2060970 100644
--- a/tests/lcms_util.c
+++ b/tests/lcms_util.c
@@ -32,6 +32,7 @@
 #include <stdlib.h>
 
 #include <libweston/matrix.h>
+#include <libweston/linalg-4.h>
 #include "shared/helpers.h"
 #include "color_util.h"
 #include "lcms_util.h"
@@ -281,11 +282,9 @@ roundtrip_verification(cmsPipeline *DToB, cmsPipeline *BToD, float tolerance)
 	test_assert_f32_lt(stat.two_norm.max, tolerance);
 }
 
-static const struct weston_vector ZEROS = {
-	.f = { 0.0, 0.0, 0.0, 1.0 }
-};
+static const struct weston_vector ZEROS = { .v = WESTON_VEC4F_ZERO };
 static const struct weston_vector PCS_BLACK = {
-	.f = {
+	.v.el = {
 		cmsPERCEPTUAL_BLACK_X,
 		cmsPERCEPTUAL_BLACK_Y,
 		cmsPERCEPTUAL_BLACK_Z,
@@ -310,19 +309,19 @@ static cmsInt32Number
 transform_sampler(const float src[], float dst[], void *cargo)
 {
 	const struct transform_sampler_context *tsc = cargo;
-	struct weston_vector stmp = { .f = { src[0], src[1], src[2], 1.0 } };
-	struct weston_vector dtmp = { .f = { 0.0, 0.0, 0.0, 1.0 } };
+	struct weston_vector stmp = { .v.el = { src[0], src[1], src[2], 1.0 } };
+	struct weston_vector dtmp = { .v.el = { 0.0, 0.0, 0.0, 1.0 } };
 
 	if (tsc->dir == BPC_DIR_BTOD)
 		weston_matrix_transform(&tsc->bpc, &stmp);
 
-	cmsDoTransform(tsc->t, stmp.f, dtmp.f, 1);
+	cmsDoTransform(tsc->t, stmp.v.el, dtmp.v.el, 1);
 
 	if (tsc->dir == BPC_DIR_DTOB)
 		weston_matrix_transform(&tsc->bpc, &dtmp);
 
 	for (int i = 0; i < 3; i++)
-		dst[i] = dtmp.f[i];
+		dst[i] = dtmp.v.el[i];
 
 	return 1; /* Success. */
 }
@@ -346,17 +345,17 @@ ComputeBlackPointCompensation(struct weston_matrix *m,
 	// a = (bpout - D50) / (bpin - D50)
 	// b = - D50* (bpout - bpin) / (bpin - D50)
 
-	tx = src_bp->f[0] - cmsD50_XYZ()->X;
-	ty = src_bp->f[1] - cmsD50_XYZ()->Y;
-	tz = src_bp->f[2] - cmsD50_XYZ()->Z;
+	tx = src_bp->v.x - cmsD50_XYZ()->X;
+	ty = src_bp->v.y - cmsD50_XYZ()->Y;
+	tz = src_bp->v.z - cmsD50_XYZ()->Z;
 
-	ax = (dst_bp->f[0] - cmsD50_XYZ()->X) / tx;
-	ay = (dst_bp->f[1] - cmsD50_XYZ()->Y) / ty;
-	az = (dst_bp->f[2] - cmsD50_XYZ()->Z) / tz;
+	ax = (dst_bp->v.x - cmsD50_XYZ()->X) / tx;
+	ay = (dst_bp->v.y - cmsD50_XYZ()->Y) / ty;
+	az = (dst_bp->v.z - cmsD50_XYZ()->Z) / tz;
 
-	bx = - cmsD50_XYZ()-> X * (dst_bp->f[0] - src_bp->f[0]) / tx;
-	by = - cmsD50_XYZ()-> Y * (dst_bp->f[1] - src_bp->f[1]) / ty;
-	bz = - cmsD50_XYZ()-> Z * (dst_bp->f[2] - src_bp->f[2]) / tz;
+	bx = - cmsD50_XYZ()-> X * (dst_bp->v.x - src_bp->v.x) / tx;
+	by = - cmsD50_XYZ()-> Y * (dst_bp->v.y - src_bp->v.y) / ty;
+	bz = - cmsD50_XYZ()-> Z * (dst_bp->v.z - src_bp->v.z) / tz;
 
 	/*
 	 *     [ax,  0,  0, bx ]
diff --git a/tests/matrix-test.c b/tests/matrix-test.c
index 2b8bd30b2..1e93a8c51 100644
--- a/tests/matrix-test.c
+++ b/tests/matrix-test.c
@@ -27,51 +27,16 @@
 
 #include <math.h>
 #include <libweston/matrix.h>
+#include <libweston/linalg-4.h>
 #include "weston-test-client-helper.h"
 #include "weston-test-assert.h"
 
-/*
- * A helper to lay out a matrix in the natural writing order in code
- * instead of needing to transpose in your mind every time you read it.
- * The matrix is laid out as written:
- *     ⎡ a11 a12 a13 a14 ⎤
- *     ⎢ a21 a22 a23 a24 ⎥
- *     ⎢ a31 a32 a33 a34 ⎥
- *     ⎣ a41 a42 a43 a44 ⎦
- * where the first digit is row and the second digit is column.
- *
- * The type field is set to the most pessimistic case possible so that if
- * weston_matrix_invert() ever gets special-case code paths, we don't take
- * them.
- */
-#define MAT(a11, a12, a13, a14,						\
-	    a21, a22, a23, a24,						\
-	    a31, a32, a33, a34,						\
-	    a41, a42, a43, a44) ((struct weston_matrix)			\
-	{								\
-		.d[0] = a11, .d[4] = a12, .d[ 8] = a13, .d[12] = a14,	\
-		.d[1] = a21, .d[5] = a22, .d[ 9] = a23, .d[13] = a24,	\
-		.d[2] = a31, .d[6] = a32, .d[10] = a33, .d[14] = a34,	\
-		.d[3] = a41, .d[7] = a42, .d[11] = a43, .d[15] = a44,	\
-		.type = WESTON_MATRIX_TRANSFORM_TRANSLATE |		\
-			WESTON_MATRIX_TRANSFORM_SCALE |			\
-			WESTON_MATRIX_TRANSFORM_ROTATE |		\
-			WESTON_MATRIX_TRANSFORM_OTHER, \
-	})
-
-static const struct weston_matrix IDENTITY =
-	MAT(1, 0, 0, 0,
-	    0, 1, 0, 0,
-	    0, 0, 1, 0,
-	    0, 0, 0, 1);
+static const struct weston_matrix IDENTITY = { .M = WESTON_MAT4F_IDENTITY };
 
 static void
 subtract_matrix(struct weston_matrix *from, const struct weston_matrix *what)
 {
-	unsigned i;
-
-	for (i = 0; i < ARRAY_LENGTH(from->d); i++)
-		from->d[i] -= what->d[i];
+	from->M = weston_m4f_sub_m4f(from->M, what->M);
 }
 
 static void
@@ -81,7 +46,7 @@ print_matrix(const struct weston_matrix *m)
 
 	for (r = 0; r < 4; ++r) {
 		for (c = 0; c < 4; ++c)
-			testlog(" %14.6e", m->d[r + c * 4]);
+			testlog(" %14.6e", m->M.col[c].el[r]);
 		testlog("\n");
 	}
 }
@@ -93,21 +58,7 @@ print_matrix(const struct weston_matrix *m)
 static double
 matrix_inf_norm(const struct weston_matrix *mat)
 {
-	unsigned row;
-	double infnorm = -1.0;
-
-	for (row = 0; row < 4; row++) {
-		unsigned col;
-		double sum = 0.0;
-
-		for (col = 0; col < 4; col++)
-			sum += fabs(mat->d[col * 4 + row]);
-
-		if (infnorm < sum)
-			infnorm = sum;
-	}
-
-	return infnorm;
+	return weston_m4f_inf_norm(mat->M);
 }
 
 struct test_matrix {
@@ -128,7 +79,8 @@ struct test_matrix {
 static const struct test_matrix matrices[] = {
 	/* A very trivial case. */
 	{
-		.M = MAT(1, 0, 0, 0,
+		.M.M = WESTON_MAT4F(
+			 1, 0, 0, 0,
 		         0, 2, 0, 0,
 		         0, 0, 3, 0,
 		         0, 0, 0, 4),
@@ -179,7 +131,8 @@ static const struct test_matrix matrices[] = {
 	 * https://gitlab.freedesktop.org/pq/fourbyfour
 	 */
 	{
-		.M = MAT(1, 0, 0, 1980,
+		.M.M = WESTON_MAT4F(
+			 1, 0, 0, 1980,
 		         0, 1, 0, 1080,
 		         0, 0, 1, 0,
 		         0, 0, 0, 1),
@@ -199,7 +152,8 @@ static const struct test_matrix matrices[] = {
 
 	/* cond = 1e3 */
 	{
-		.M = MAT(-4.12798022231678357619e-02, -7.93301899046665176529e-02, 2.49367040174418935772e-01, -2.22400462135059429070e-01,
+		.M.M = WESTON_MAT4F(
+			 -4.12798022231678357619e-02, -7.93301899046665176529e-02, 2.49367040174418935772e-01, -2.22400462135059429070e-01,
 		          2.02416121867255743849e-01, -2.25754422240346010187e-02, -2.91283152417864787953e-01, 1.49354988316431153139e-01,
 		          6.18473094065821293874e-01, 5.81511312950217934548e-02, -1.18363610818063924590e+00, 8.00087538947595322547e-01,
 		          1.25723127083294305972e-01, 7.72723720984487272290e-02, -3.76023220287807879991e-01, 2.82473279931768073148e-01),
@@ -208,7 +162,8 @@ static const struct test_matrix matrices[] = {
 
 	/* cond = 1e3, abs det = 15 */
 	{
-		.M = MAT(6.84154939885726509630e+00, -6.87241565273813304060e+00, -2.56772939909334070308e+01, -2.52185055099662420730e+01,
+		.M.M = WESTON_MAT4F(
+			 6.84154939885726509630e+00, -6.87241565273813304060e+00, -2.56772939909334070308e+01, -2.52185055099662420730e+01,
 		         2.04511561406330022450e+00, -3.67551043874248994925e+00, -1.96421641406619129633e+00, -2.40644091603848320204e+00,
 		         5.83631095663641819016e+00, -9.31051765621826277197e+00, -1.80402129629135217215e+01, -1.78475057662460052654e+01,
 		        -9.88588496379959025262e+00, 1.49790516545410774540e+01, 2.64975800675967363418e+01, 2.65795891678410747261e+01),
@@ -217,7 +172,8 @@ static const struct test_matrix matrices[] = {
 
 	/* cond = 700, abs det = 1e-6, invertible regardless of det */
 	{
-		.M = MAT(1.32125189257677579449e-03, -1.67411409720826992453e-01, 1.07940907587735196449e-01, -1.22163309792902186057e-01,
+		.M.M = WESTON_MAT4F(
+			 1.32125189257677579449e-03, -1.67411409720826992453e-01, 1.07940907587735196449e-01, -1.22163309792902186057e-01,
 		        -5.42113793774764013422e-02, 5.30455105336593901733e-01, -2.59607412684229155175e-01, 4.36480803188117993940e-01,
 		         2.88175168292948129939e-03, -1.85262537685181277736e-01, 1.46265858042118279680e-01, -9.41398969709369287662e-02,
 		        -2.88900393087768159184e-03, 1.57987202530630227448e-01, -1.20781192010860280450e-01, 8.95194304475115387731e-02),
@@ -226,7 +182,8 @@ static const struct test_matrix matrices[] = {
 
 	/* cond = 1e6, this is a little more challenging */
 	{
-		.M = MAT(-4.41851445093878913983e-01, -5.16386185043831491548e-01, 2.86186055948129847160e-01, -5.79440137716940473211e-01,
+		.M.M = WESTON_MAT4F(
+			 -4.41851445093878913983e-01, -5.16386185043831491548e-01, 2.86186055948129847160e-01, -5.79440137716940473211e-01,
 		          2.49798696238173301154e-01, 2.84965614532234345901e-01, -1.65729639683955931595e-01, 3.12568045963485974248e-01,
 		          3.15253213984537428161e-01, 3.71270066781250074328e-01, -2.02675623845341434937e-01, 4.19969870491003371971e-01,
 		          5.60818677658178832424e-01, 6.45373659426444201692e-01, -3.68902466471524526082e-01, 7.13785795079988516498e-01),
@@ -235,7 +192,8 @@ static const struct test_matrix matrices[] = {
 
 	/* cond = 15, abs det = 1e-9, should be well invertible */
 	{
-		.M = MAT(-5.37536200142514660589e-05, 7.92552373388843642288e-03, -3.90554524958281433500e-03, 2.68892064500873568395e-03,
+		.M.M = WESTON_MAT4F(
+			 -5.37536200142514660589e-05, 7.92552373388843642288e-03, -3.90554524958281433500e-03, 2.68892064500873568395e-03,
 		         -9.72329428437283989350e-03, 8.32075145342783470404e-03, 6.52648485926096092596e-03, 1.06707947887298994737e-03,
 		          1.04453728969657322345e-02, -1.03627268579679666927e-02, -3.56835980207569763989e-03, -3.95935925157862422114e-03,
 		          5.37160838929722633805e-03, 6.13466744624343262009e-05, -1.23695935407398946090e-04, 8.21231194921675112380e-04),
diff --git a/tests/matrix-transform-test.c b/tests/matrix-transform-test.c
index 6c94e70d8..2a07ec6a8 100644
--- a/tests/matrix-transform-test.c
+++ b/tests/matrix-transform-test.c
@@ -34,6 +34,7 @@
 #include <wayland-client.h>
 #include "libweston-internal.h"
 #include "libweston/matrix.h"
+#include <libweston/linalg-4.h>
 
 #include "weston-test-client-helper.h"
 #include "weston-test-assert.h"
@@ -382,40 +383,40 @@ simple_transform_vector(struct weston_output *output, struct weston_vector in)
 
 	switch (output->transform) {
 	case WL_OUTPUT_TRANSFORM_NORMAL:
-		out.f[0] = (-output->pos.c.x + in.f[0]) * scale;
-		out.f[1] = (-output->pos.c.y + in.f[1]) * scale;
+		out.v.el[0] = (-output->pos.c.x + in.v.el[0]) * scale;
+		out.v.el[1] = (-output->pos.c.y + in.v.el[1]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_FLIPPED:
-		out.f[0] = (output->pos.c.x + output->width - in.f[0]) * scale;
-		out.f[1] = (-output->pos.c.y + in.f[1]) * scale;
+		out.v.el[0] = (output->pos.c.x + output->width - in.v.el[0]) * scale;
+		out.v.el[1] = (-output->pos.c.y + in.v.el[1]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_90:
-		out.f[0] = (-output->pos.c.y + in.f[1]) * scale;
-		out.f[1] = (output->pos.c.x + output->width - in.f[0]) * scale;
+		out.v.el[0] = (-output->pos.c.y + in.v.el[1]) * scale;
+		out.v.el[1] = (output->pos.c.x + output->width - in.v.el[0]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_FLIPPED_90:
-		out.f[0] = (-output->pos.c.y + in.f[1]) * scale;
-		out.f[1] = (-output->pos.c.x + in.f[0]) * scale;
+		out.v.el[0] = (-output->pos.c.y + in.v.el[1]) * scale;
+		out.v.el[1] = (-output->pos.c.x + in.v.el[0]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_180:
-		out.f[0] = (output->pos.c.x + output->width - in.f[0]) * scale;
-		out.f[1] = (output->pos.c.y + output->height - in.f[1]) * scale;
+		out.v.el[0] = (output->pos.c.x + output->width - in.v.el[0]) * scale;
+		out.v.el[1] = (output->pos.c.y + output->height - in.v.el[1]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_FLIPPED_180:
-		out.f[0] = (-output->pos.c.x + in.f[0]) * scale;
-		out.f[1] = (output->pos.c.y + output->height - in.f[1]) * scale;
+		out.v.el[0] = (-output->pos.c.x + in.v.el[0]) * scale;
+		out.v.el[1] = (output->pos.c.y + output->height - in.v.el[1]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_270:
-		out.f[0] = (output->pos.c.y + output->height - in.f[1]) * scale;
-		out.f[1] = (-output->pos.c.x + in.f[0]) * scale;
+		out.v.el[0] = (output->pos.c.y + output->height - in.v.el[1]) * scale;
+		out.v.el[1] = (-output->pos.c.x + in.v.el[0]) * scale;
 		break;
 	case WL_OUTPUT_TRANSFORM_FLIPPED_270:
-		out.f[0] = (output->pos.c.y + output->height - in.f[1]) * scale;
-		out.f[1] = (output->pos.c.x + output->width - in.f[0]) * scale;
+		out.v.el[0] = (output->pos.c.y + output->height - in.v.el[1]) * scale;
+		out.v.el[1] = (output->pos.c.x + output->width - in.v.el[0]) * scale;
 		break;
 	}
-	out.f[2] = 0;
-	out.f[3] = 1;
+	out.v.el[2] = 0;
+	out.v.el[3] = 1;
 
 	return out;
 }
@@ -426,7 +427,7 @@ output_test_all_transforms(struct weston_output *output,
 {
 	int i;
 	int transform;
-	struct weston_vector t = { { 7.0, 13.0, 0.0, 1.0 } };
+	struct weston_vector t = { .v = WESTON_VEC4F(7.0, 13.0, 0.0, 1.0) };
 	struct weston_vector v, sv;
 
 	for (transform = WL_OUTPUT_TRANSFORM_NORMAL;
@@ -443,7 +444,7 @@ output_test_all_transforms(struct weston_output *output,
 		weston_matrix_transform(&output->matrix, &v);
 		sv = simple_transform_vector(output, t);
 		for (i = 0; i < 4; i++)
-			test_assert_f32_eq (sv.f[i], v.f[i]);
+			test_assert_f32_eq(sv.v.el[i], v.v.el[i]);
 	}
 }