diff --git a/algebra/csc_tools/csc_math.h b/algebra/csc_tools/csc_math.h
index 20e10514c..3058aba55 100644
--- a/algebra/csc_tools/csc_math.h
+++ b/algebra/csc_tools/csc_math.h
@@ -36,16 +36,15 @@ void csc_update_values(csc           *M,
 /*****************************************************************************
 * CSC Algebraic Operations                                    ******************************************************************************/
 
-//DEBUG : ADD documentation
-
-/* matrix times scalar */
-
+// A = sc*A
 void csc_scale(csc* A, c_float sc);
 
+// A = diag(L)*A
 void csc_lmult_diag(csc* A, const c_float *L);
 
+// A = A*diag(R)
 void csc_rmult_diag(csc* A, const c_float *R);
-    
+
 //y = alpha*A*x + beta*y, where A is symmetric and only triu is stored
 void csc_Axpy_sym_triu(const csc   *A,
                        const c_float *x,
@@ -60,19 +59,23 @@ void csc_Axpy(const csc   *A,
                     c_float alpha,
                     c_float beta);
 
-
+//y = alpha*A^T*x + beta*y
 void csc_Atxpy(const csc *A,
                      const c_float *x,
                      c_float *y,
                      c_float alpha,
                      c_float beta);
 
+// returns 1/2 x'*P*x
 c_float csc_quad_form(const csc *P, const c_float *x);
 
+// E[i] = inf_norm(M(:,i))
 void csc_col_norm_inf(const csc *M, c_float *E);
 
+// E[i] = inf_norm(M(i,:))
 void csc_row_norm_inf(const csc *M, c_float *E);
 
+// E[i] = inf_norm(M(i,:)), where M stores triu part only
 void csc_row_norm_inf_sym_triu(const csc *M, c_float *E);
 
 
diff --git a/algebra/default/matrix.c b/algebra/default/matrix.c
index 7e0bc4cca..52ab167ed 100644
--- a/algebra/default/matrix.c
+++ b/algebra/default/matrix.c
@@ -4,7 +4,7 @@
 #include "csc_math.h"
 #include "csc_utils.h"
 
-/*  logical functions ------------------------------------------------------*/
+/*  logical test functions ----------------------------------------------------*/
 
 c_int OSQPMatrix_is_eq(OSQPMatrix *A, OSQPMatrix* B, c_float tol){
   return (A->symmetry == B->symmetry &&
@@ -53,7 +53,7 @@ c_int    OSQPMatrix_get_n(const OSQPMatrix *M){return M->csc->n;}
 c_float* OSQPMatrix_get_x(const OSQPMatrix *M){return M->csc->x;}
 c_int*   OSQPMatrix_get_i(const OSQPMatrix *M){return M->csc->i;}
 c_int*   OSQPMatrix_get_p(const OSQPMatrix *M){return M->csc->p;}
-c_int    OSQPMatrix_get_nnz(const OSQPMatrix *M){return M->csc->p[M->csc->n];}
+c_int    OSQPMatrix_get_nz(const OSQPMatrix *M){return M->csc->p[M->csc->n];}
 
 
 /* math functions ----------------------------------------------------------*/
diff --git a/include/algebra_matrix.h b/include/algebra_matrix.h
index 81ab63eee..ea54c7ad2 100644
--- a/include/algebra_matrix.h
+++ b/include/algebra_matrix.h
@@ -33,8 +33,8 @@ OSQPMatrix* OSQPMatrix_new_from_csc(const csc* A, c_int is_triu);
 
 /*  direct data access functions ---------------------------------------------*/
 
-/*  These functions allow getting/setting data
-*   in the OSQPMatrix type.   Data is passed in/out using bare
+/*  These functions allow getting data in csc format from the
+*   the OSQPMatrix type.   Data is passed in/out using bare
 *   pointers instead of OSQPVectors since these functions interface
 *   with user defined linear solvers and the user API
 */
@@ -44,12 +44,23 @@ void OSQPMatrix_update_values(OSQPMatrix    *M,
                             const c_int     *Mx_new_idx,
                             c_int           M_new_n);
 
+/* returns the row dimension */
 c_int    OSQPMatrix_get_m(const OSQPMatrix *M);
+
+/* returns the columns dimension */
 c_int    OSQPMatrix_get_n(const OSQPMatrix *M);
+
+/* returns a pointer to the array of data values */
 c_float* OSQPMatrix_get_x(const OSQPMatrix *M);
+
+/* returns a pointer to the array of row indices */
 c_int*   OSQPMatrix_get_i(const OSQPMatrix *M);
+
+/* returns a pointer to the array of col indices (csc format).  Should be n+1 long */
 c_int*   OSQPMatrix_get_p(const OSQPMatrix *M);
-c_int    OSQPMatrix_get_nnz(const OSQPMatrix *M);
+
+/* returns the number of nonzeros (length of x and i arrays) */
+c_int    OSQPMatrix_get_nz(const OSQPMatrix *M);
 
 
 /* math functions ----------------------------------------------------------*/
diff --git a/include/algebra_vector.h b/include/algebra_vector.h
index 05744086c..0a798c2c1 100644
--- a/include/algebra_vector.h
+++ b/include/algebra_vector.h
@@ -170,7 +170,7 @@ c_float OSQPVectorf_mean(const OSQPVectorf *a);
 c_float OSQPVectorf_dot_prod(const OSQPVectorf *a,
                             const OSQPVectorf *b);
 
-/* Inner product a'b, but using only the positive or Negative
+/* Inner product a'b, but using only the positive or negative
  * terms in b.   Use sign = 1 for positive terms, sign = -1 for
  * negative terms.   Setting any other value for sign will return
  * the normal dot product
@@ -184,7 +184,9 @@ void OSQPVectorf_ew_prod(OSQPVectorf       *c,
                          const OSQPVectorf *a,
                          const OSQPVectorf *b);
 
-/* check l <= u elementwise */
+/* check l <= u elementwise.  Returns 1 if inequality is true
+ * for every element pair in both vectors  
+ */
 c_int OSQPVectorf_all_leq(OSQPVectorf *l, OSQPVectorf* u);
 
 /* Elementwise bounding vectors x = min(max(z,l),u)
diff --git a/lin_sys/direct/pardiso/pardiso_interface.c b/lin_sys/direct/pardiso/pardiso_interface.c
index 7d4adb5e9..63f894e5a 100644
--- a/lin_sys/direct/pardiso/pardiso_interface.c
+++ b/lin_sys/direct/pardiso/pardiso_interface.c
@@ -137,8 +137,8 @@ c_int init_linsys_solver_pardiso(pardiso_solver ** sp, const OSQPMatrix * P, con
     else { // Called from ADMM algorithm
 
         // Allocate vectors of indices
-        s->PtoKKT = c_malloc(OSQPMatrix_get_nnz(P) * sizeof(c_int));
-        s->AtoKKT = c_malloc(OSQPMatrix_get_nnz(A) * sizeof(c_int));
+        s->PtoKKT = c_malloc(OSQPMatrix_get_nz(P) * sizeof(c_int));
+        s->AtoKKT = c_malloc(OSQPMatrix_get_nz(A) * sizeof(c_int));
         s->rhotoKKT = c_malloc(OSQPMatrix_get_m(A) * sizeof(c_int));
 
         // Use s->rho_inv_vec for storing param2 = rho_inv_vec
diff --git a/lin_sys/direct/qdldl/qdldl_interface.c b/lin_sys/direct/qdldl/qdldl_interface.c
index 27bfa9726..e369b9677 100644
--- a/lin_sys/direct/qdldl/qdldl_interface.c
+++ b/lin_sys/direct/qdldl/qdldl_interface.c
@@ -271,8 +271,8 @@ c_int init_linsys_solver_qdldl(qdldl_solver ** sp, const OSQPMatrix* P, const OS
     else { // Called from ADMM algorithm
 
         // Allocate vectors of indices
-        s->PtoKKT = c_malloc(OSQPMatrix_get_nnz(P) * sizeof(c_int));
-        s->AtoKKT = c_malloc(OSQPMatrix_get_nnz(A) * sizeof(c_int));
+        s->PtoKKT = c_malloc(OSQPMatrix_get_nz(P) * sizeof(c_int));
+        s->AtoKKT = c_malloc(OSQPMatrix_get_nz(A) * sizeof(c_int));
         s->rhotoKKT = c_malloc(OSQPMatrix_get_m(A) * sizeof(c_int));
 
         // Use p->rho_inv_vec for storing param2 = rho_inv_vec
@@ -295,7 +295,7 @@ c_int init_linsys_solver_qdldl(qdldl_solver ** sp, const OSQPMatrix* P, const OS
 
         // Permute matrix
         if (KKT_temp)
-            permute_KKT(&KKT_temp, s, OSQPMatrix_get_nnz(P), OSQPMatrix_get_nnz(A), s->m, s->PtoKKT, s->AtoKKT, s->rhotoKKT);
+            permute_KKT(&KKT_temp, s, OSQPMatrix_get_nz(P), OSQPMatrix_get_nz(A), s->m, s->PtoKKT, s->AtoKKT, s->rhotoKKT);
     }
 
     // Check if matrix has been created
diff --git a/src/osqp_api.c b/src/osqp_api.c
index 4128ac471..2fb2ba638 100644
--- a/src/osqp_api.c
+++ b/src/osqp_api.c
@@ -1038,7 +1038,7 @@ c_int osqp_update_P(OSQPSolver    *solver,
   osqp_tic(work->timer); // Start timer
 #endif /* ifdef PROFILING */
 
-  nnzP = OSQPMatrix_get_nnz(work->data->P);
+  nnzP = OSQPMatrix_get_nz(work->data->P);
 
   if (Px_new_idx) { // Passing the index of elements changed
     // Check if number of elements is less or equal than the total number of
@@ -1109,7 +1109,7 @@ c_int osqp_update_A(OSQPSolver *solver,
   osqp_tic(work->timer); // Start timer
 #endif /* ifdef PROFILING */
 
-  nnzA = OSQPMatrix_get_nnz(work->data->A);
+  nnzA = OSQPMatrix_get_nz(work->data->A);
 
   if (Ax_new_idx) { // Passing the index of elements changed
     // Check if number of elements is less or equal than the total number of
@@ -1184,8 +1184,8 @@ c_int osqp_update_P_A(OSQPSolver    *solver,
   osqp_tic(work->timer); // Start timer
 #endif /* ifdef PROFILING */
 
-  nnzP = OSQPMatrix_get_nnz(work->data->P);
-  nnzA = OSQPMatrix_get_nnz(work->data->A);
+  nnzP = OSQPMatrix_get_nz(work->data->P);
+  nnzA = OSQPMatrix_get_nz(work->data->A);
 
 
   if (Px_new_idx) { // Passing the index of elements changed
diff --git a/src/util.c b/src/util.c
index 25f68ce37..a68737467 100644
--- a/src/util.c
+++ b/src/util.c
@@ -70,7 +70,7 @@ void print_setup_header(const OSQPSolver *solver) {
   settings = solver->settings;
 
   // Number of nonzeros
-  nnz = OSQPMatrix_get_nnz(data->P) + OSQPMatrix_get_nnz(data->A);
+  nnz = OSQPMatrix_get_nz(data->P) + OSQPMatrix_get_nz(data->A);
 
   print_line();
   c_print("           OSQP v%s  -  Operator Splitting QP Solver\n"