diff --git a/demos/gtk-demo/four_point_transform.c b/demos/gtk-demo/four_point_transform.c
index 9df10ef163..beffc174ae 100644
--- a/demos/gtk-demo/four_point_transform.c
+++ b/demos/gtk-demo/four_point_transform.c
@@ -63,11 +63,15 @@ unit_to (graphene_point3d_t *p1,
   graphene_matrix_multiply (&s, &u, m);
 }
 
-/* Make a 4x4 matrix that maps
+/* Compute a 4x4 matrix m that maps
  * p1 -> q1
  * p2 -> q2
  * p3 -> q3
  * p4 -> q4
+ *
+ * This is not in general possible, because projective
+ * transforms preserve coplanarity. But in the cases we
+ * care about here, both sets of points are always coplanar.
  */
 void
 perspective_3d (graphene_point3d_t *p1,
diff --git a/demos/gtk-demo/singular_value_decomposition.c b/demos/gtk-demo/singular_value_decomposition.c
index 05d04d0667..e6d063ab4b 100644
--- a/demos/gtk-demo/singular_value_decomposition.c
+++ b/demos/gtk-demo/singular_value_decomposition.c
@@ -11,6 +11,20 @@
 
 #define MAX_ITERATION_COUNT 30
 
+/* Perform Householder reduction to bidiagonal form
+ *
+ * Input: Matrix A of size nrows x ncols
+ *
+ * Output: Matrices and vectors such that
+ * A = U*Bidiag(diagonal, superdiagonal)*Vt
+ *
+ * All matrices are allocated by the caller
+ *
+ * Sizes:
+ *  A, U: nrows x ncols
+ *  diagonal, superdiagonal: ncols
+ *  V: ncols x ncols
+ */
 static void
 householder_reduction (double *A,
                        int     nrows,
@@ -160,6 +174,20 @@ householder_reduction (double *A,
     }
 }
 
+/* Perform Givens reduction
+ *
+ * Input: Matrices such that
+ * A = U*Bidiag(diagonal,superdiagonal)*Vt
+ *
+ * Output: The same, with superdiagonal = 0
+ *
+ * All matrices are allocated by the caller
+ *
+ * Sizes:
+ *  U: nrows x ncols
+ *  diagonal, superdiagonal: ncols
+ *  V: ncols x ncols
+ */
 static int
 givens_reduction (int nrows,
                   int ncols,
@@ -298,6 +326,11 @@ givens_reduction (int nrows,
   return 0;
 }
 
+/* Given a singular value decomposition
+ * of an nrows x ncols matrix A = U*Diag(S)*Vt,
+ * sort the values of S by decreasing value,
+ * permuting V to match.
+ */
 static void
 sort_singular_values (int     nrows,
                       int     ncols,
@@ -339,6 +372,16 @@ sort_singular_values (int     nrows,
     }
 }
 
+/* Compute a singular value decomposition of A,
+ * A = U*Diag(S)*Vt
+ *
+ * All matrices are allocated by the caller
+ *
+ * Sizes:
+ *  A, U: nrows x ncols
+ *  S: ncols
+ *  V: ncols x ncols
+ */
 int
 singular_value_decomposition (double *A,
                               int     nrows,
@@ -364,6 +407,18 @@ singular_value_decomposition (double *A,
   return 0;
 }
 
+/*
+ * Given a singular value decomposition of A = U*Diag(S)*Vt,
+ * compute the best approximation x to A*x = B.
+ *
+ * All matrices are allocated by the caller
+ *
+ * Sizes:
+ *  U: nrows x ncols
+ *  S: ncols
+ *  V: ncols x ncols
+ *  B, x: ncols
+ */
 void
 singular_value_decomposition_solve (double *U,
                                     double *S,