diff --git a/src/spasm.h b/src/spasm.h
index 17ce048..4bfff96 100644
--- a/src/spasm.h
+++ b/src/spasm.h
@@ -161,6 +161,8 @@ u32 spasm_prng_u32(spasm_prng_ctx *ctx);
 spasm_ZZp spasm_prng_ZZp(spasm_prng_ctx *ctx);
 
 /* spasm_util.c */
+extern int (*logcallback)(char *);
+int logprintf(char *format, ...);
 double spasm_wtime();
 i64 spasm_nnz(const struct spasm_csr * A);
 void *spasm_malloc(i64 size);
diff --git a/src/spasm_echelonize.c b/src/spasm_echelonize.c
index 8c0778c..8e42886 100644
--- a/src/spasm_echelonize.c
+++ b/src/spasm_echelonize.c
@@ -40,7 +40,7 @@ bool spasm_echelonize_test_completion(const struct spasm_csr *A, const int *p, i
 	void *S = spasm_malloc(Sn * Sm * spasm_datatype_size(datatype));
 	int *q = spasm_malloc(Sm * sizeof(*q));
 	size_t *Sp = spasm_malloc(Sm * sizeof(*Sp));       /* for FFPACK */
-	fprintf(stderr, "[echelonize/completion] Testing completion with %" PRId64" random linear combinations (rank %d)\n", Sn, U->n);
+	logprintf("[echelonize/completion] Testing completion with %" PRId64" random linear combinations (rank %d)\n", Sn, U->n);
 	fflush(stderr);
 	spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, datatype, q, Sn, 0);
 	int rr = spasm_ffpack_rref(prime, Sn, Sm, S, Sm, datatype, Sp);
@@ -58,7 +58,7 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 	int m = A->m;
 	int r = spasm_min(A->n, m);  /* upper-bound on rank */
 	int verbose_step = spasm_max(1, n / 1000);
-	fprintf(stderr, "[echelonize/GPLU] processing matrix of dimension %d x %d\n", n, m);
+	logprintf("[echelonize/GPLU] processing matrix of dimension %d x %d\n", n, m);
 	
 	struct spasm_csr *U = fact->U;
 	struct spasm_triplet *L = fact->Ltmp;
@@ -83,12 +83,12 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 	for (i = 0; i < n; i++) {
 		/* test for early abort (if L not needed) */
 		if (L == NULL && U->n == r) {
-			fprintf(stderr, "\n[echelonize/GPLU] full rank reached\n");
+			logprintf("\n[echelonize/GPLU] full rank reached\n");
 			break;
 		}
 		/* TODO: make these hard-coded values options */
 		if (L == NULL && !early_abort_done && rows_since_last_pivot > 10 && (rows_since_last_pivot > (n / 100))) {
-			fprintf(stderr, "\n[echelonize/GPLU] testing for early abort...\n");
+			logprintf("\n[echelonize/GPLU] testing for early abort...\n");
 			if (spasm_echelonize_test_completion(A, p, n, U, Uqinv))
 				break;
 			early_abort_done = 1;
@@ -145,13 +145,13 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 
 		/* store new pivotal row into U */
 		Uqinv[jpiv] = U->n;
-		// fprintf(stderr, "setting Uqinv[%d] <--- %d\n", jpiv, U->n);
+		// logprintf("setting Uqinv[%d] <--- %d\n", jpiv, U->n);
 
 		i64 old_unz = unz;
 		Uj[unz] = jpiv;
 		Ux[unz] = 1;
 		unz += 1;
-		// fprintf(stderr, "setting U[%d, %d] <--- 1\n", U->n, jpiv);
+		// logprintf("setting U[%d, %d] <--- 1\n", U->n, jpiv);
 		assert(x[jpiv] != 0);
 		spasm_ZZp beta = spasm_ZZp_inverse(A->field, x[jpiv]);
 		for (int px = top; px < m; px++) {
@@ -159,7 +159,7 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 			if (x[j] != 0 && Uqinv[j] < 0) {
 				Uj[unz] = j;
 				Ux[unz] = spasm_ZZp_mul(A->field, beta, x[j]);
-				// fprintf(stderr, "setting U[%d, %d] <--- %d\n", U->n, j, Ux[unz]);
+				// logprintf("setting U[%d, %d] <--- %d\n", U->n, j, Ux[unz]);
 				unz += 1;
 			}
 		}
@@ -171,7 +171,7 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 		early_abort_done = 0;
 
 		if ((i % verbose_step) == 0) {
-			fprintf(stderr, "\r[echelonize/GPLU] %d / %d [|U| = %" PRId64 " / |L| = %" PRId64"] -- current density= (%.3f vs %.3f) --- rank >= %d", 
+			logprintf("\r[echelonize/GPLU] %d / %d [|U| = %" PRId64 " / |L| = %" PRId64"] -- current density= (%.3f vs %.3f) --- rank >= %d", 
 				i, n, unz, lnz, 1.0 * (m - top) / (m), 1.0 * (unz - old_unz) / m, U->n);
 			fflush(stderr);
 		}
@@ -181,7 +181,7 @@ static void echelonize_GPLU(const struct spasm_csr *A, const int *p, int n, cons
 		L->nz = lnz;
 		L->m = U->n;
 	}
-	fprintf(stderr, "\n");
+	logprintf("\n");
 	free(x);
 	free(xj);
 }
@@ -196,7 +196,7 @@ static void update_U_after_rref(int rr, int Sm, const void *S, spasm_datatype da
 	int *Uqinv = fact->qinv;
 	i64 extra_nnz = ((i64) (1 + Sm - rr)) * rr;     /* maximum size increase */
 	i64 unz = spasm_nnz(U);
-	fprintf(stderr, "[dense update] enlarging U from %" PRId64 " to %" PRId64 " entries\n", unz, unz + extra_nnz);
+	logprintf("[dense update] enlarging U from %" PRId64 " to %" PRId64 " entries\n", unz, unz + extra_nnz);
 	spasm_csr_realloc(U, unz + extra_nnz);
 	i64 *Up = U->p;
 	int *Uj = U->j;
@@ -267,7 +267,7 @@ static void update_fact_after_LU(int n, int Sm, int r, const void *S, spasm_data
 			Lx[lnz] = x;
 			lnz += 1;
 		}
-		fprintf(stderr, "L : %" PRId64 " --> %" PRId64 " --> %" PRId64 " ---> ", lnz_before, L->nz, lnz);
+		logprintf("L : %" PRId64 " --> %" PRId64 " --> %" PRId64 " ---> ", lnz_before, L->nz, lnz);
 	}
 
 	/* add new entries from S */
@@ -287,7 +287,7 @@ static void update_fact_after_LU(int n, int Sm, int r, const void *S, spasm_data
 			Lp[U->n + i] = iorig;
 	}
 	L->nz = lnz;
-	fprintf(stderr, "%" PRId64 "\n", lnz);
+	logprintf("%" PRId64 "\n", lnz);
 
 	/* fill U */
 	for (i64 i = 0; i < r; i++) {
@@ -329,7 +329,7 @@ static void echelonize_dense_lowrank(const struct spasm_csr *A, const int *p, in
 	double start = spasm_wtime();
 	int old_un = U->n;
 	int round = 0;
-	fprintf(stderr, "[echelonize/dense/low-rank] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n", 
+	logprintf("[echelonize/dense/low-rank] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n", 
 		n, Sm, opts->dense_block_size, spasm_datatype_name(datatype));
 	
 	/* 
@@ -349,7 +349,7 @@ static void echelonize_dense_lowrank(const struct spasm_csr *A, const int *p, in
 		int Sn = spasm_min(rank_ub, opts->dense_block_size);
 		if (Sn <= 0)
 			break;		
-		fprintf(stderr, "[echelonize/dense/low-rank] Round %d. Weight %d. Processing chunk (%d x %d), |U| = %"PRId64"\n", 
+		logprintf("[echelonize/dense/low-rank] Round %d. Weight %d. Processing chunk (%d x %d), |U| = %"PRId64"\n", 
 			round, w, Sn, Sm, spasm_nnz(U));
 		spasm_schur_dense_randomized(A, p, n, U, Uqinv, S, datatype, q, Sn, w);
 		int rr = spasm_ffpack_rref(prime, Sn, Sm, S, Sm, datatype, Sp);
@@ -357,22 +357,22 @@ static void echelonize_dense_lowrank(const struct spasm_csr *A, const int *p, in
 		if (rr == 0) {
 			if (spasm_echelonize_test_completion(A, p, n, U, Uqinv))
 				break;
-			fprintf(stderr, "[echelonize/dense/low-rank] Failed termination test; switching to full linear combinations\n");
+			logprintf("[echelonize/dense/low-rank] Failed termination test; switching to full linear combinations\n");
 			w = 0;
 			Sn = omp_get_max_threads();
 		}
 		if (rr < 0.9 * Sn) {
 			w *= 2;
-			fprintf(stderr, "[echelonize/dense/low-rank] Not enough pivots, increasing weight to %d\n", w);
+			logprintf("[echelonize/dense/low-rank] Not enough pivots, increasing weight to %d\n", w);
 		}
 		update_U_after_rref(rr, Sm, S, datatype, Sp, q, fact);
 		n -= rr;
 		Sm -= rr;
 		rank_ub -= rr;
 		round += 1;
-		fprintf(stderr, "[echelonize/dense/low-rank] found %d new pivots (%d new since beginning)\n", rr,  U->n - old_un);
+		logprintf("[echelonize/dense/low-rank] found %d new pivots (%d new since beginning)\n", rr,  U->n - old_un);
 	}
-	fprintf(stderr, "[echelonize/dense/low-rank] completed in %.1fs. %d new pivots found\n", spasm_wtime() - start, U->n - old_un);
+	logprintf("[echelonize/dense/low-rank] completed in %.1fs. %d new pivots found\n", spasm_wtime() - start, U->n - old_un);
 	free(S);
 	free(q);
 	free(Sp);
@@ -404,7 +404,7 @@ static void echelonize_dense(const struct spasm_csr *A, const int *p, int n, con
 	double start = spasm_wtime();
 	int old_un = U->n;
 	int round = 0;
-	fprintf(stderr, "[echelonize/dense] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n", 
+	logprintf("[echelonize/dense] processing dense schur complement of dimension %d x %d; block size=%d, type %s\n", 
 		n, Sm, opts->dense_block_size, spasm_datatype_name(datatype));
 	bool lowrank_mode = 0;
 	int rank_ub = spasm_min(A->n - U->n, A->m - U->n);
@@ -415,7 +415,7 @@ static void echelonize_dense(const struct spasm_csr *A, const int *p, int n, con
 		if (Sn <= 0)
 			break;
 		
-		fprintf(stderr, "[echelonize/dense] Round %d. processing S[%d:%d] (%d x %d)\n", round, processed, processed + Sn, Sn, Sm);	
+		logprintf("[echelonize/dense] Round %d. processing S[%d:%d] (%d x %d)\n", round, processed, processed + Sn, Sn, Sm);	
 
 		i64 lnz_before = (opts->L) ? fact->Ltmp->nz : -1;
 		spasm_schur_dense(A, p, Sn, p_in, fact, S, datatype, q, p_out);
@@ -437,7 +437,7 @@ static void echelonize_dense(const struct spasm_csr *A, const int *p, int n, con
 		p += Sn;
 		Sm = m - U->n;
 		rank_ub = spasm_min(A->n - U->n, A->m - U->n);
-		fprintf(stderr, "[echelonize/dense] found %d new pivots\n", rr);
+		logprintf("[echelonize/dense] found %d new pivots\n", rr);
 
 		/* 
 		 * switch to low-rank mode if yield drops too much 
@@ -455,10 +455,10 @@ static void echelonize_dense(const struct spasm_csr *A, const int *p, int n, con
 	free(p_out);
 	free(pivotal);
 	if (rank_ub > 0 && n - processed > 0 && lowrank_mode) {
-		fprintf(stderr, "[echelonize/dense] Too few pivots; switching to low-rank mode\n");
+		logprintf("[echelonize/dense] Too few pivots; switching to low-rank mode\n");
 		echelonize_dense_lowrank(A, p, n - processed, fact, opts);
 	} else {
-		fprintf(stderr, "[echelonize/dense] completed in %.1fs. %d new pivots found\n", spasm_wtime() - start, U->n - old_un);
+		logprintf("[echelonize/dense] completed in %.1fs. %d new pivots found\n", spasm_wtime() - start, U->n - old_un);
 	}
 }
 
@@ -474,14 +474,14 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 {
 	struct echelonize_opts default_opts;
 	if (opts == NULL) {
-		fprintf(stderr, "[echelonize] using default settings\n");
+		logprintf("[echelonize] using default settings\n");
 		opts = &default_opts;
 		spasm_echelonize_init_opts(opts);
 	}
 	int n = A->n;
 	int m = A->m;
 	i64 prime = spasm_get_prime(A);
-	fprintf(stderr, "[echelonize] Start on %d x %d matrix with %" PRId64 " nnz\n", n, m, spasm_nnz(A));
+	logprintf("[echelonize] Start on %d x %d matrix with %" PRId64 " nnz\n", n, m, spasm_nnz(A));
 	
 	/* options sanity check */
 	if (opts->complete)
@@ -525,26 +525,26 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 	for (round = 0; round < opts->max_round; round++) {
 		/* decide whether to move on to the next iteration */
 		if (spasm_nnz(A) == 0) {
-			fprintf(stderr, "[echelonize] empty matrix\n");
+			logprintf("[echelonize] empty matrix\n");
 			status = 1;
 			break;
 		}
 		
-		fprintf(stderr, "[echelonize] round %d\n", round);
+		logprintf("[echelonize] round %d\n", round);
 		npiv = spasm_pivots_extract_structural(A, p_in, fact, p, opts);
 
 		if (npiv < opts->min_pivot_proportion * spasm_min(n, m - U->n)) {
-			fprintf(stderr, "[echelonize] not enough pivots found; stopping\n");
+			logprintf("[echelonize] not enough pivots found; stopping\n");
 			status = 2;
 			break;     /* not enough pivots found */
 		}		
 		// if (density > opts->sparsity_threshold && aspect_ratio > opts->tall_and_skinny_ratio) {
-		// 	fprintf(stderr, "Schur complement is dense, tall and skinny (#rows / #cols = %.1f)\n", aspect_ratio);
+		// 	logprintf("Schur complement is dense, tall and skinny (#rows / #cols = %.1f)\n", aspect_ratio);
 		// 	break;
 		// }
 		density = spasm_schur_estimate_density(A, p + npiv, n - npiv, U, Uqinv, 100);
 		if (density > opts->sparsity_threshold) {
-			fprintf(stderr, "[echelonize] Schur complement is dense (estimated %.2f%%)\n", 100 * density);
+			logprintf("[echelonize] Schur complement is dense (estimated %.2f%%)\n", 100 * density);
 			status = 2;
 			break;
 		}
@@ -553,7 +553,7 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 		i64 nnz = (density * (n - npiv)) * (m - U->n);
 		char tmp[8];
 		spasm_human_format(sizeof(int) * (n - npiv + nnz) + sizeof(spasm_ZZp) * nnz, tmp);
-		fprintf(stderr, "Schur complement is %d x %d, estimated density : %.2f (%s byte)\n", n - npiv, m - U->n, density, tmp);
+		logprintf("Schur complement is %d x %d, estimated density : %.2f (%s byte)\n", n - npiv, m - U->n, density, tmp);
 		int *p_out = spasm_malloc((n - npiv) * sizeof(*p_out));
 		struct spasm_csr *S = spasm_schur(A, p + npiv, n - npiv, fact, density, L, p_in, p_out);
 		if (round > 0)
@@ -579,14 +579,14 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 
 	/* finish */
 	if (!opts->enable_tall_and_skinny)
-		fprintf(stderr, "[echelonize] dense low-rank mode disabled\n");
+		logprintf("[echelonize] dense low-rank mode disabled\n");
 	if (!opts->enable_dense)
-		fprintf(stderr, "[echelonize] regular dense mode disabled\n");
+		logprintf("[echelonize] regular dense mode disabled\n");
 	if (!opts->enable_GPLU)
-		fprintf(stderr, "[echelonize] GPLU mode disabled\n");
+		logprintf("[echelonize] GPLU mode disabled\n");
 	
 	double aspect_ratio = (double) (n - npiv) / (m - U->n);
-	fprintf(stderr, "[echelonize] finishing; density = %.3f; aspect ratio = %.1f\n", density, aspect_ratio);
+	logprintf("[echelonize] finishing; density = %.3f; aspect ratio = %.1f\n", density, aspect_ratio);
 	if (opts->enable_tall_and_skinny && aspect_ratio > opts->tall_and_skinny_ratio)
 		echelonize_dense_lowrank(A, p + npiv, n - npiv, fact, opts);
 	else if (opts->enable_dense && density > opts->sparsity_threshold)
@@ -594,12 +594,12 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 	else if (opts->enable_GPLU)
 		echelonize_GPLU(A, p + npiv, n - npiv, p_in, fact, opts);
 	else
-		fprintf(stderr, "[echelonize] Cannot finish (no valid method enabled). Incomplete echelonization returned\n");
+		logprintf("[echelonize] Cannot finish (no valid method enabled). Incomplete echelonization returned\n");
 
 cleanup:
 	free(p);
 	free(p_in);
-	fprintf(stderr, "[echelonize] Done in %.1fs. Rank %d, %" PRId64 " nz in basis\n", spasm_wtime() - start, U->n, spasm_nnz(U));
+	logprintf("[echelonize] Done in %.1fs. Rank %d, %" PRId64 " nz in basis\n", spasm_wtime() - start, U->n, spasm_nnz(U));
 	spasm_csr_resize(U, U->n, m);
 	spasm_csr_realloc(U, -1);
 	if (round > 0)
@@ -614,4 +614,4 @@ struct spasm_lu * spasm_echelonize(const struct spasm_csr *A, struct echelonize_
 	}
 	fact->r = U->n;
 	return fact;
-}
\ No newline at end of file
+}
diff --git a/src/spasm_io.c b/src/spasm_io.c
index 162470d..96cf862 100644
--- a/src/spasm_io.c
+++ b/src/spasm_io.c
@@ -93,7 +93,7 @@ struct spasm_triplet *spasm_triplet_load(FILE * f, i64 prime, u8 *hash)
 		if (sscanf(buffer, "%d %d %" SCNd64 "\n", &i, &j, &nnz) != 3)
 			errx(1, "[spasm_triplet_load] bad MatrixMarking dimensions (line %" PRId64 ")\n", line);
 		spasm_human_format(nnz, hnnz);
-		fprintf(stderr, "[IO] loading %d x %d MatrixMarket matrix modulo %" PRId64 " with %s non-zero... ", 
+		logprintf("[IO] loading %d x %d MatrixMarket matrix modulo %" PRId64 " with %s non-zero... ", 
 			i, j, prime, hnnz);
 		fflush(stderr);
 	} else {
@@ -103,7 +103,7 @@ struct spasm_triplet *spasm_triplet_load(FILE * f, i64 prime, u8 *hash)
 			errx(1, "[spasm_triplet_load] bad SMS file (header)\n");
 		if (prime != -1 && type != 'M')
 			errx(1, "[spasm_triplet_load] only ``Modular'' type supported\n");
-		fprintf(stderr, "[IO] loading %d x %d SMS matrix modulo %" PRId64 "... ", i, j, prime);
+		logprintf("[IO] loading %d x %d SMS matrix modulo %" PRId64 "... ", i, j, prime);
 		fflush(stderr);
 	}
 
@@ -143,17 +143,17 @@ struct spasm_triplet *spasm_triplet_load(FILE * f, i64 prime, u8 *hash)
 	if (!mm) {
 		spasm_triplet_realloc(T, -1);
 		spasm_human_format(T->nz, hnnz);
-		fprintf(stderr, "%s non-zero [%.1fs]\n", hnnz, spasm_wtime() - start);
+		logprintf("%s non-zero [%.1fs]\n", hnnz, spasm_wtime() - start);
 	} else {
-		fprintf(stderr, "[%.1fs]\n", spasm_wtime() - start);
+		logprintf("[%.1fs]\n", spasm_wtime() - start);
 	}
 
 	if (ctx != NULL) {
 		spasm_SHA256_final(hash, ctx);
-		fprintf(stderr, "[spasm_triplet_load] sha256(matrix) = ");
+		logprintf("[spasm_triplet_load] sha256(matrix) = ");
 		for (int i = 0; i < 32; i++)
-			fprintf(stderr, "%02x", hash[i]);
-		fprintf(stderr, " / size = %" PRId64" bytes\n", (((i64) ctx->Nh) << 29) + ctx->Nl / 8);
+			logprintf("%02x", hash[i]);
+		logprintf(" / size = %" PRId64" bytes\n", (((i64) ctx->Nh) << 29) + ctx->Nl / 8);
 	}
 	return T;
 }
@@ -250,7 +250,7 @@ void spasm_save_pnm(const struct spasm_csr *A, FILE *f, int x, int y, int mode,
 		limits_h[1] = ((long long int) cc[3]) * ((long long int) x) / ((long long int) m);
 		limits_v[0] = ((long long int) rr[1]) * ((long long int) y) / ((long long int) n);
 		limits_v[1] = ((long long int) rr[2]) * ((long long int) y) / ((long long int) n);
-		fprintf(stderr, "limits_v = %d, %d\n", limits_v[0], limits_v[1]);
+		logprintf("limits_v = %d, %d\n", limits_v[0], limits_v[1]);
 		r = DM->r;
 		c = DM->c;
 	}
diff --git a/src/spasm_kernel.c b/src/spasm_kernel.c
index 7bc4a06..580d829 100644
--- a/src/spasm_kernel.c
+++ b/src/spasm_kernel.c
@@ -16,7 +16,7 @@ struct spasm_csr * spasm_kernel(const struct spasm_lu *fact)
 	assert(n <= m);
 	char hnnz[8];
 	spasm_human_format(spasm_nnz(U), hnnz);
-	fprintf(stderr, "[kernel] start. U is %d x %d (%s nnz). Transposing U\n", n, m, hnnz);
+	logprintf("[kernel] start. U is %d x %d (%s nnz). Transposing U\n", n, m, hnnz);
 	double start_time = spasm_wtime();
 
 	struct spasm_csr *Ut = spasm_transpose(U, true);
@@ -73,7 +73,7 @@ struct spasm_csr * spasm_kernel(const struct spasm_lu *fact)
 					while (writing > 0) {
 						#pragma omp flush(writing)
 					}
-					fprintf(stderr, "increasing K; %" PRId64 " ---> %" PRId64 "\n", K->nzmax, 2 * K->nzmax + row_nz);
+					logprintf("increasing K; %" PRId64 " ---> %" PRId64 "\n", K->nzmax, 2 * K->nzmax + row_nz);
 					spasm_csr_realloc(K, 2 * K->nzmax + row_nz);
 					Kj = K->j;
 					Kx = K->x;
@@ -110,7 +110,7 @@ struct spasm_csr * spasm_kernel(const struct spasm_lu *fact)
 			if (tid == 0) {
 				char hnnz[8];
 				spasm_human_format(nnz, hnnz);
-	  			fprintf(stderr, "\rkernel: %d/%d, |K| = %s    ", Kn, m-n, hnnz);
+	  			logprintf("\rkernel: %d/%d, |K| = %s    ", Kn, m-n, hnnz);
 	  			fflush(stderr);
 	  		}
 		}
@@ -118,11 +118,11 @@ struct spasm_csr * spasm_kernel(const struct spasm_lu *fact)
 		free(xj);
 	}
 
-	fprintf(stderr, "\n");
+	logprintf("\n");
 	free(Utqinv);
 	spasm_csr_free(Ut);
 	spasm_human_format(spasm_nnz(K), hnnz);
-	fprintf(stderr, "[kernel] done in %.1fs. NNZ(K) = %s\n", spasm_wtime() - start_time, hnnz);
+	logprintf("[kernel] done in %.1fs. NNZ(K) = %s\n", spasm_wtime() - start_time, hnnz);
 	return K;
 }
 
@@ -176,4 +176,4 @@ struct spasm_csr * spasm_kernel_from_rref(const struct spasm_csr *R, const int *
 	free(p);
 	spasm_csr_free(Rt);
 	return K;
-}
\ No newline at end of file
+}
diff --git a/src/spasm_matching.c b/src/spasm_matching.c
index 1163286..652883b 100644
--- a/src/spasm_matching.c
+++ b/src/spasm_matching.c
@@ -124,10 +124,10 @@ int spasm_maximum_matching(const struct spasm_csr *A, int *p, int *qinv)
 	for (int i = 0; (i < n) && (k < r); i++) {
 		if (p[i] < 0)
 			k += spasm_augmenting_path(A, i, istack, jstack, pstack, marks, plookahead, p, qinv);
-		fprintf(stderr, "\r[matching] %d / %d, size %d", i, n, k);
+		logprintf("\r[matching] %d / %d, size %d", i, n, k);
 		fflush(stderr);
 	}
-	fprintf(stderr, " [%.1f s]\n", spasm_wtime() - start);
+	logprintf(" [%.1f s]\n", spasm_wtime() - start);
 	free(istack);
 	free(jstack);
 	free(pstack);
@@ -190,4 +190,4 @@ int *spasm_submatching(const int *match, int a, int b, int c, int d) {
 		}
 	}
 	return pmatch;
-}
\ No newline at end of file
+}
diff --git a/src/spasm_pivots.c b/src/spasm_pivots.c
index 9068718..a8eaa24 100644
--- a/src/spasm_pivots.c
+++ b/src/spasm_pivots.c
@@ -61,7 +61,7 @@ static int spasm_find_FL_pivots(const struct spasm_csr *A, int *p, int *qinv)
 		if (qinv[j] == -1 || spasm_row_weight(A, i) < spasm_row_weight(A, qinv[j]))
 			npiv += register_pivot(i, j, p, qinv);
 	}
-	fprintf(stderr, "[pivots] Faugère-Lachartre: %d pivots found [%.1fs]\n", npiv, spasm_wtime() - start);
+	logprintf("[pivots] Faugère-Lachartre: %d pivots found [%.1fs]\n", npiv, spasm_wtime() - start);
 	return npiv;
 }
 
@@ -116,7 +116,7 @@ static int spasm_find_FL_column_pivots(const struct spasm_csr *A, int *pinv, int
 		}
 	}
 	free(w);
-	fprintf(stderr, "[pivots] ``Faugère-Lachartre on columns'': %d pivots found [%.1fs]\n", 
+	logprintf("[pivots] ``Faugère-Lachartre on columns'': %d pivots found [%.1fs]\n", 
 		npiv, spasm_wtime() - start);
 	return npiv;
 }
@@ -186,7 +186,7 @@ static int spasm_find_cycle_free_pivots(const struct spasm_csr *A, int *pinv, in
 			 *   w[j] ==  0  column j was absent and is unreachable
 			 */
 			if ((tid == 0) && (i % v) == 0) {
-				fprintf(stderr, "\r[pivots] %d / %d --- found %d new", processed, n, npiv);
+				logprintf("\r[pivots] %d / %d --- found %d new", processed, n, npiv);
 				fflush(stderr);
 			}
 			if (pinv[i] >= 0)
@@ -288,7 +288,7 @@ static int spasm_find_cycle_free_pivots(const struct spasm_csr *A, int *pinv, in
 		free(queue);
 	}
 	free(journal);
-	fprintf(stderr, "\r[pivots] greedy alternating cycle-free search: %d pivots found [%.1fs]\n", 
+	logprintf("\r[pivots] greedy alternating cycle-free search: %d pivots found [%.1fs]\n", 
 		npiv, spasm_wtime() - start);
 	return npiv;
 }
@@ -314,7 +314,7 @@ static int spasm_pivots_find(const struct spasm_csr *A, int *pinv, int *qinv, st
 	npiv += spasm_find_FL_column_pivots(A, pinv, qinv);	
 	if (opts->enable_greedy_pivot_search)
 		npiv += spasm_find_cycle_free_pivots(A, pinv, qinv);
-	fprintf(stderr, "\r[pivots] %d pivots found\n", npiv);
+	logprintf("\r[pivots] %d pivots found\n", npiv);
 	return npiv;
 }
 
@@ -421,7 +421,7 @@ int spasm_pivots_extract_structural(const struct spasm_csr *A, const int *p_in,
 		if (L != NULL) {
 			int i_out = (p_in != NULL) ? p_in[i] : i;
 			spasm_add_entry(L, i_out, U->n, pivot);
-			// fprintf(stderr, "Adding L[%d, %d] = %d\n", i_out, U->n, pivot);
+			// logprintf("Adding L[%d, %d] = %d\n", i_out, U->n, pivot);
 			Lp[U->n] = i_out;
 		}
 
@@ -476,4 +476,4 @@ struct spasm_csr *spasm_permute_pivots(const struct spasm_csr *A, const int *p,
 		}
 	return spasm_permute(A, p, qinv, true);
 }
-#endif
\ No newline at end of file
+#endif
diff --git a/src/spasm_rref.c b/src/spasm_rref.c
index 21d7b9f..6c12a4a 100644
--- a/src/spasm_rref.c
+++ b/src/spasm_rref.c
@@ -29,7 +29,7 @@ struct spasm_csr * spasm_rref(const struct spasm_lu *fact, int *Rqinv)
 	i64 prime = spasm_get_prime(U);
 	char hnnz[8];
 	spasm_human_format(spasm_nnz(U), hnnz);
-	fprintf(stderr, "[rref] start. U is %d x %d (%s nnz)\n", n, m, hnnz);
+	logprintf("[rref] start. U is %d x %d (%s nnz)\n", n, m, hnnz);
 	double start_time = spasm_wtime();
 	struct spasm_csr *R = spasm_csr_alloc(n, m, spasm_nnz(U), prime, true);
 	i64 *Rp = R->p;
@@ -120,7 +120,7 @@ struct spasm_csr * spasm_rref(const struct spasm_lu *fact, int *Rqinv)
 			if (tid == 0) {
 				char hnnz[8];
 				spasm_human_format(nnz, hnnz);
-	  			fprintf(stderr, "\rRREF: %d/%d, |R| = %s    ", Rn, n, hnnz);
+	  			logprintf("\rRREF: %d/%d, |R| = %s    ", Rn, n, hnnz);
 	  			fflush(stderr);
 	  		}
 		}
@@ -139,9 +139,9 @@ struct spasm_csr * spasm_rref(const struct spasm_lu *fact, int *Rqinv)
 			Rqinv[j] = i;
 		}
 	}
-	fprintf(stderr, "\n");
+	logprintf("\n");
 
 	spasm_human_format(spasm_nnz(R), hnnz);
-	fprintf(stderr, "[rref] done in %.1fs. NNZ(R) = %s\n", spasm_wtime() - start_time, hnnz);
+	logprintf("[rref] done in %.1fs. NNZ(R) = %s\n", spasm_wtime() - start_time, hnnz);
 	return R;
-}
\ No newline at end of file
+}
diff --git a/src/spasm_schur.c b/src/spasm_schur.c
index 80b3423..0baf82b 100644
--- a/src/spasm_schur.c
+++ b/src/spasm_schur.c
@@ -165,7 +165,7 @@ struct spasm_csr *spasm_schur(const struct spasm_csr *A, const int *p, int n, co
 					Li[local_lnz] = i_orig;
 					Lj[local_lnz] = qinv[j];
 					Lx[local_lnz] = x[j];
-					// fprintf(stderr, "Adding L[%d, %d] = %d\n", i_out, qinv[j], x[j]);
+					// logprintf("Adding L[%d, %d] = %d\n", i_out, qinv[j], x[j]);
 					local_lnz += 1;
 				}
 			}
@@ -176,7 +176,7 @@ struct spasm_csr *spasm_schur(const struct spasm_csr *A, const int *p, int n, co
 
 			if (tid == 0 && (i % verbose_step) == 0) {
 				double density =  1.0 * snz / (1.0 * m * Sn);
-				fprintf(stderr, "\rSchur complement: %d/%d [%" PRId64 " nz / density= %.3f]", Sn, n, snz, density);
+				logprintf("\rSchur complement: %d/%d [%" PRId64 " nz / density= %.3f]", Sn, n, snz, density);
 				fflush(stderr);
 			}
 		}
@@ -188,7 +188,7 @@ struct spasm_csr *spasm_schur(const struct spasm_csr *A, const int *p, int n, co
 		L->nz = lnz;
 	spasm_csr_realloc(S, -1);
 	double density = 1.0 * snz / (1.0 * m * n);
-	fprintf(stderr, "\rSchur complement: %d * %d [%" PRId64 " nz / density= %.3f], %.1fs\n", n, m, snz, density, spasm_wtime() - start);
+	logprintf("\rSchur complement: %d * %d [%" PRId64 " nz / density= %.3f], %.1fs\n", n, m, snz, density, spasm_wtime() - start);
 	return S;
 }
 
@@ -263,7 +263,7 @@ void spasm_schur_dense(const struct spasm_csr *A, const int *p, int n, const int
 	int m = A->m;
 	int Sm = m - U->n;                                   /* #columns of S */
 	prepare_q(m, qinv, q);                               /* FIXME: useless if many invokations */
-	fprintf(stderr, "[schur/dense] dimension %d x %d...\n", n, Sm);
+	logprintf("[schur/dense] dimension %d x %d...\n", n, Sm);
 	double start = spasm_wtime();
 	int verbose_step = spasm_max(1, n / 1000);
 	int r = 0;
@@ -322,14 +322,14 @@ void spasm_schur_dense(const struct spasm_csr *A, const int *p, int n, const int
 			#pragma omp atomic update
 			r += 1;
 			if (tid == 0 && (r % verbose_step) == 0) {
-				fprintf(stderr, "\r[schur/dense] %d/%d", r, n);
+				logprintf("\r[schur/dense] %d/%d", r, n);
 				fflush(stderr);
 			}
 		}
 		free(x);
 		free(xj);
 	}
-	fprintf(stderr, "\n[schur/dense] finished in %.1fs, rank <= %d\n", spasm_wtime() - start, r);
+	logprintf("\n[schur/dense] finished in %.1fs, rank <= %d\n", spasm_wtime() - start, r);
 }
 
 
@@ -354,7 +354,7 @@ void spasm_schur_dense_randomized(const struct spasm_csr *A, const int *p, int n
 	const i64 *Up = U->p;
 	const int *Uj = U->j;
 	prepare_q(m, qinv, q);
-	fprintf(stderr, "[schur/dense/random] dimension %d x %d, weight %d...\n", N, Sm, w);
+	logprintf("[schur/dense/random] dimension %d x %d, weight %d...\n", N, Sm, w);
 	double start = spasm_wtime();
 	int verbose_step = spasm_max(1, N / 1000);
 
@@ -403,11 +403,11 @@ void spasm_schur_dense_randomized(const struct spasm_csr *A, const int *p, int n
 
 			/* verbosity */
 			if ((k % verbose_step) == 0) {
-				fprintf(stderr, "\r[schur/dense/random] %" PRId64 "/%d", k, N);
+				logprintf("\r[schur/dense/random] %" PRId64 "/%d", k, N);
 				fflush(stderr);
 			}
 		}
 		free(y);
 	}
-	fprintf(stderr, "\n[schur/dense/random] finished in %.1fs\n", spasm_wtime() - start);
-}
\ No newline at end of file
+	logprintf("\n[schur/dense/random] finished in %.1fs\n", spasm_wtime() - start);
+}
diff --git a/src/spasm_triplet.c b/src/spasm_triplet.c
index 9c1cb58..14c568b 100644
--- a/src/spasm_triplet.c
+++ b/src/spasm_triplet.c
@@ -106,7 +106,7 @@ struct spasm_csr *spasm_compress(const struct spasm_triplet * T)
 	spasm_ZZp *Tx = T->x;
 	
 	double start = spasm_wtime();
-	fprintf(stderr, "[CSR] Compressing... ");
+	logprintf("[CSR] Compressing... ");
 	fflush(stderr);
 
 	/* allocate result */
@@ -153,6 +153,6 @@ struct spasm_csr *spasm_compress(const struct spasm_triplet * T)
 	char mem[16];
 	int size = sizeof(int) * (n + nz) + sizeof(spasm_ZZp) * ((Cx != NULL) ? nz : 0);
 	spasm_human_format(size, mem);
-	fprintf(stderr, "%" PRId64 " actual NZ, Mem usage = %sbyte [%.2fs]\n", spasm_nnz(C), mem, spasm_wtime() - start);
+	logprintf("%" PRId64 " actual NZ, Mem usage = %sbyte [%.2fs]\n", spasm_nnz(C), mem, spasm_wtime() - start);
 	return C;
-}
\ No newline at end of file
+}
diff --git a/src/spasm_util.c b/src/spasm_util.c
index b5710dc..ebb908e 100644
--- a/src/spasm_util.c
+++ b/src/spasm_util.c
@@ -4,9 +4,24 @@
 #include <sys/time.h>
 #include <err.h>
 #include <inttypes.h>
+#include <stdarg.h>
 
 #include "spasm.h"
 
+int (*logcallback)(char *) = NULL;
+
+int logprintf(char *format, ...) {
+  va_list args;
+  va_start(args, format);
+
+  if (logcallback) {
+    char s[1000];
+    vsnprintf(s, 1000, format, args);
+    return (*logcallback)(s);
+  } else
+    return vfprintf(stderr, format, args);
+}
+
 int spasm_get_num_threads() {
 #ifdef _OPENMP
 	return omp_get_num_threads();
@@ -140,12 +155,12 @@ void spasm_triplet_realloc(struct spasm_triplet *A, i64 nzmax)
 {
 	if (nzmax < 0)
 		nzmax = A->nz;
-	// fprintf(stderr, "[realloc] nzmax=%ld. before %px %px %px\n", nzmax, A->i, A->j, A->x);
+	// logprintf("[realloc] nzmax=%ld. before %px %px %px\n", nzmax, A->i, A->j, A->x);
 	A->i = spasm_realloc(A->i, nzmax * sizeof(int));
 	A->j = spasm_realloc(A->j, nzmax * sizeof(int));
 	if (A->x != NULL)
 		A->x = spasm_realloc(A->x, nzmax * sizeof(spasm_ZZp));
-	// fprintf(stderr, "[realloc] after %px %px %px\n", A->i, A->j, A->x);
+	// logprintf("[realloc] after %px %px %px\n", A->i, A->j, A->x);
 	A->nzmax = nzmax;
 }