Switch to using dsyrk instead of dsyr for speed.

Also added a workspace (GenWork) structure for to hold working matrices
for the gensvm_get_update() and gensvm_get_loss() functions

6 files changed, 282 insertions, 65 deletions

diff --git a/include/gensvm_base.h b/include/gensvm_base.h
index efeaa4d..03b7ffa 100644
--- a/include/gensvm_base.h
+++ b/include/gensvm_base.h
@@ -99,12 +99,43 @@ struct GenModel {
 	///< array of kernel parameters, size depends on kernel type
 };
 
+/**
+ * @brief A structure to hold the GenSVM workspace
+ *
+ */
+struct GenWork {
+	long n;
+	///< number of instances for the workspace
+	long m;
+	///< number of features for the workspace
+	long K;
+	///< number of classes for the workspace
+
+	double *LZ;
+	///< n x (m+1) working matrix for the Z'*A*Z calculation
+	double *ZB;
+	///< (m+1) x (K-1) working matrix for the Z'*B calculation
+	double *ZBc;
+	///< (K-1) x (m+1) working matrix for the Z'*B calculation
+	double *ZAZ;
+	///< (m+1) x (m+1) working matrix for the Z'*A*Z calculation
+	double *ZV;
+	///< n x (K-1) working matrix for the Z * V calculation
+	double *beta;
+	///< K-1 working vector for a row of the B matrix
+};
+
 // function declarations
 struct GenModel *gensvm_init_model();
 void gensvm_allocate_model(struct GenModel *model);
 void gensvm_reallocate_model(struct GenModel *model, long n, long m);
 void gensvm_free_model(struct GenModel *model);
+
 struct GenData *gensvm_init_data();
 void gensvm_free_data(struct GenData *data);
 
+struct GenWork *gensvm_init_work(struct GenModel *model);
+void gensvm_free_work(struct GenWork *work);
+void gensvm_reset_work(struct GenWork *work);
+
 #endif
diff --git a/include/gensvm_optimize.h b/include/gensvm_optimize.h
index bbdf4c8..dec8914 100644
--- a/include/gensvm_optimize.h
+++ b/include/gensvm_optimize.h
@@ -19,8 +19,8 @@
 
 // function declarations
 void gensvm_optimize(struct GenModel *model, struct GenData *data);
-double gensvm_get_loss(struct GenModel *model, struct GenData *data,
-		double *ZV);
+double gensvm_get_loss(struct GenModel *model, struct GenData *data, 
+		struct GenWork *work);
 
 double gensvm_calculate_omega(struct GenModel *model, struct GenData *data,
 		long i);
@@ -33,7 +33,8 @@ void gensvm_calculate_ab_simple(struct GenModel *model, long i, long j,
 double gensvm_get_alpha_beta(struct GenModel *model, struct GenData *data,
 		long i, double *beta);
 
-void gensvm_get_update(struct GenModel *model, struct GenData *data);
+void gensvm_get_update(struct GenModel *model, struct GenData *data, 
+		struct GenWork *work);
 void gensvm_calculate_errors(struct GenModel *model, struct GenData *data,
 		double *ZV);
 void gensvm_calculate_huber(struct GenModel *model);
diff --git a/src/gensvm_base.c b/src/gensvm_base.c
index c77e704..0f3ac5d 100644
--- a/src/gensvm_base.c
+++ b/src/gensvm_base.c
@@ -195,3 +195,86 @@ void gensvm_free_model(struct GenModel *model)
 	free(model);
 }
 
+/**
+ * @brief Initialize the workspace structure
+ *
+ * @details
+ * During the computations in gensvm_get_update(), a number of matrices need 
+ * to be allocated which are used to store intermediate results. To avoid 
+ * reallocating these matrices at every call to gensvm_get_update(), we create 
+ * here a workspace with these matrices. This workspace is created once by 
+ * gensvm_optimize(), and is passed to gensvm_get_update() and 
+ * gensvm_get_loss(). See the documentation of the GenWork structure for 
+ * information on each allocated field.
+ *
+ * @param[in] 	model 	a GenModel with the dimensionality of the problem
+ * @returns 		an allocated GenWork instance
+ *
+ */
+struct GenWork *gensvm_init_work(struct GenModel *model)
+{
+	long n = model->n;
+	long m = model->m;
+	long K = model->K;
+
+	struct GenWork *work = Malloc(struct GenWork, 1);
+	work->n = n;
+	work->m = m;
+	work->K = K;
+
+	work->LZ = Calloc(double, n*(m+1));
+	work->ZB = Calloc(double, (m+1)*(K-1)),
+	work->ZBc = Calloc(double, (m+1)*(K-1)),
+	work->ZAZ = Calloc(double, (m+1)*(m+1)),
+	work->ZV = Calloc(double, n*(K-1));
+	work->beta = Calloc(double, K-1);
+
+	return work;
+}
+
+/**
+ * @brief Free an allocated GenWork instance
+ *
+ * @details
+ * This function simply frees every matrix allocated for in the GenWork 
+ * workspace.
+ *
+ * @param[in] 	work 	a pointer to an allocated GenWork instance
+ *
+ */
+void gensvm_free_work(struct GenWork *work)
+{
+	free(work->LZ);
+	free(work->ZB);
+	free(work->ZBc);
+	free(work->ZAZ);
+	free(work->ZV);
+	free(work->beta);
+	free(work);
+}
+
+/**
+ * @brief Reset all matrices of a GenWork instance
+ *
+ * @details
+ * In the gensvm_get_update() function, it is expected for some matrices that 
+ * all their entries are initialized to 0. This function sets the memory for 
+ * each of the matrices in the GenWork workspace to 0 to facilitate this.
+ *
+ * @param[in,out] 	work 	an initialized GenWork instance. On exit, the
+ * 				matrices of the GenWork instance are all 
+ * 				initialized to 0.
+ */
+void gensvm_reset_work(struct GenWork *work)
+{
+	long n = work->n;
+	long m = work->m;
+	long K = work->K;
+
+	Memset(work->LZ, double, n*(m+1));
+	Memset(work->ZB, double, (m+1)*(K-1)),
+	Memset(work->ZBc, double, (m+1)*(K-1)),
+	Memset(work->ZAZ, double, (m+1)*(m+1)),
+	Memset(work->ZV, double, n*(K-1));
+	Memset(work->beta, double, K-1);
+}
diff --git a/src/gensvm_optimize.c b/src/gensvm_optimize.c
index 1022ed4..184c45f 100644
--- a/src/gensvm_optimize.c
+++ b/src/gensvm_optimize.c
@@ -47,8 +47,10 @@ void gensvm_optimize(struct GenModel *model, struct GenData *data)
 	long m = model->m;
 	long K = model->K;
 
-	double *ZV = Calloc(double, n*(K-1));
+	// initialize the workspace
+	struct GenWork *work = gensvm_init_work(model);
 
+	// print some info on the dataset and model configuration
 	note("Starting main loop.\n");
 	note("Dataset:\n");
 	note("\tn = %i\n", n);
@@ -61,28 +63,33 @@ void gensvm_optimize(struct GenModel *model, struct GenData *data)
 	note("\tepsilon = %g\n", model->epsilon);
 	note("\n");
 
+	// compute necessary simplex vectors
 	gensvm_simplex(model);
 	gensvm_simplex_diff(model);
 
-	L = gensvm_get_loss(model, data, ZV);
+	// get initial loss
+	L = gensvm_get_loss(model, data, work);
 	Lbar = L + 2.0*model->epsilon*L;
 
+	// run main loop
 	while ((it < MAX_ITER) && (Lbar - L)/L > model->epsilon)
 	{
 		// ensures V contains newest V and Vbar contains V from
 		// previous
-		gensvm_get_update(model, data);
+		gensvm_get_update(model, data, work);
 		if (it > 50)
 			gensvm_step_doubling(model);
 
 		Lbar = L;
-		L = gensvm_get_loss(model, data, ZV);
+		L = gensvm_get_loss(model, data, work);
 
 		if (it%PRINT_ITER == 0)
 			note("iter = %li, L = %15.16f, Lbar = %15.16f, "
 			     "reldiff = %15.16f\n", it, L, Lbar, (Lbar - L)/L);
 		it++;
 	}
+
+	// print warnings if necessary
 	if (L > Lbar)
 		err("[GenSVM Warning]: Negative step occurred in "
 				"majorization.\n");
@@ -90,13 +97,19 @@ void gensvm_optimize(struct GenModel *model, struct GenData *data)
 		err("[GenSVM Warning]: maximum number of iterations "
 				"reached.\n");
 
+	// print final iteration count and loss
 	note("Optimization finished, iter = %li, loss = %15.16f, "
 			"rel. diff. = %15.16f\n", it-1, L,
 			(Lbar - L)/L);
+
+	// compute and print the number of SVs in the model
 	note("Number of support vectors: %li\n", gensvm_num_sv(model));
 
+	// store the training error in the model
 	model->training_error = (Lbar - L)/L;
-	free(ZV);
+
+	// free the workspace
+	gensvm_free_work(work);
 }
 
 /**
@@ -110,12 +123,11 @@ void gensvm_optimize(struct GenModel *model, struct GenData *data)
  * @param[in] 		model 	GenModel structure which holds the current
  * 				estimate V
  * @param[in]  		data 	GenData structure
- * @param[in,out] 	ZV 	pre-allocated matrix ZV which is updated on
- * 				output
+ * @param[in] 		work 	allocated workspace with the ZV matrix to use
  * @returns 			the current value of the loss function
  */
-double gensvm_get_loss(struct GenModel *model, struct GenData *data,
-		double *ZV)
+double gensvm_get_loss(struct GenModel *model, struct GenData *data, 
+		struct GenWork *work)
 {
 	long i, j;
 	long n = model->n;
@@ -124,7 +136,7 @@ double gensvm_get_loss(struct GenModel *model, struct GenData *data,
 
 	double value, rowvalue, loss = 0.0;
 
-	gensvm_calculate_errors(model, data, ZV);
+	gensvm_calculate_errors(model, data, work->ZV);
 	gensvm_calculate_huber(model);
 
 	for (i=0; i<n; i++) {
@@ -434,10 +446,12 @@ double gensvm_get_alpha_beta(struct GenModel *model, struct GenData *data,
  * @todo
  * Consider using CblasColMajor everywhere
  *
- * @param [in,out] 	model 	model to be updated
- * @param [in] 		data 	data used in model
+ * @param[in,out] 	model 	model to be updated
+ * @param[in] 		data 	data used in model
+ * @param[in] 		work 	allocated workspace to use
  */
-void gensvm_get_update(struct GenModel *model, struct GenData *data)
+void gensvm_get_update(struct GenModel *model, struct GenData *data, 
+		struct GenWork *work)
 {
 	int status;
 	long i, j;
@@ -447,52 +461,53 @@ void gensvm_get_update(struct GenModel *model, struct GenData *data)
 	long m = model->m;
 	long K = model->K;
 
-	// allocate working matrices and vectors
-	double *ZB = Calloc(double, (m+1)*(K-1)),
-	       *ZAZ = Calloc(double, (m+1)*(m+1)),
-	       *ZBc = Calloc(double, (m+1)*(K-1)),
-	       *beta = Calloc(double, K-1);
+	gensvm_reset_work(work);
 
 	// generate Z'*A*Z and Z'*B by rank 1 operations
 	for (i=0; i<n; i++) {
-		alpha = gensvm_get_alpha_beta(model, data, i, beta);
+		alpha = gensvm_get_alpha_beta(model, data, i, work->beta);
 
-		// rank 1 update of symmetric matrix Z'*A*Z (A is diagonal)
-		// Note: only the upper triangular part is calculated
-		cblas_dsyr(CblasRowMajor, CblasUpper, m+1, alpha,
-				&data->Z[i*(m+1)], 1, ZAZ, m+1);
+		// calculate row of matrix LZ, which is a scalar 
+		// multiplication of sqrt(alpha_i) and row z_i' of Z
+		cblas_daxpy(m+1, sqrt(alpha), &data->Z[i*(m+1)], 1, 
+				&work->LZ[i*(m+1)], 1);
 
 		// rank 1 update of matrix Z'*B
 		// Note: LDA is the second dimension of ZB because of
 		// Row-Major order
 		cblas_dger(CblasRowMajor, m+1, K-1, 1, &data->Z[i*(m+1)], 1,
-				beta, 1, ZB, K-1);
+				work->beta, 1, work->ZB, K-1);
 	}
 
+	// calculate Z'*A*Z by symmetric multiplication of LZ with itself 
+	// (ZAZ = (LZ)' * (LZ)
+	cblas_dsyrk(CblasRowMajor, CblasUpper, CblasTrans, m+1, n, 1.0,
+			work->LZ, m+1, 0.0, work->ZAZ, m+1);
+
 	// Calculate right-hand side of system we want to solve
 	// dsymm performs ZB := 1.0 * (ZAZ) * Vbar + 1.0 * ZB
 	// the right-hand side is thus stored in ZB after this call
 	// Note: LDB and LDC are second dimensions of the matrices due to
 	// Row-Major order
-	cblas_dsymm(CblasRowMajor, CblasLeft, CblasUpper, m+1, K-1, 1, ZAZ,
-			m+1, model->V, K-1, 1.0, ZB, K-1);
+	cblas_dsymm(CblasRowMajor, CblasLeft, CblasUpper, m+1, K-1, 1, 
+			work->ZAZ, m+1, model->V, K-1, 1.0, work->ZB, K-1);
 
 	// Calculate left-hand side of system we want to solve
 	// Add lambda to all diagonal elements except the first one. Recall
 	// that ZAZ is of size m+1 and is symmetric.
 	for (i=m+2; i<=m*(m+2); i+=m+2)
-		ZAZ[i] += model->lambda;
+		work->ZAZ[i] += model->lambda;
 
 	// Lapack uses column-major order, so we transform the ZB matrix to
 	// correspond to this.
 	for (i=0; i<m+1; i++)
 		for (j=0; j<K-1; j++)
-			ZBc[j*(m+1)+i] = ZB[i*(K-1)+j];
+			work->ZBc[j*(m+1)+i] = work->ZB[i*(K-1)+j];
 
 	// Solve the system using dposv. Note that above the upper triangular
 	// part has always been used in row-major order for ZAZ. This
 	// corresponds to the lower triangular part in column-major order.
-	status = dposv('L', m+1, K-1, ZAZ, m+1, ZBc, m+1);
+	status = dposv('L', m+1, K-1, work->ZAZ, m+1, work->ZBc, m+1);
 
 	// Use dsysv as fallback, for when the ZAZ matrix is not positive
 	// semi-definite for some reason (perhaps due to rounding errors).
@@ -502,13 +517,13 @@ void gensvm_get_update(struct GenModel *model, struct GenData *data)
 				"dposv: %i\n", status);
 		int *IPIV = Malloc(int, m+1);
 		double *WORK = Malloc(double, 1);
-		status = dsysv('L', m+1, K-1, ZAZ, m+1, IPIV, ZBc, m+1, WORK,
-				-1);
+		status = dsysv('L', m+1, K-1, work->ZAZ, m+1, IPIV, work->ZBc, 
+				m+1, WORK, -1);
 
 		int LWORK = WORK[0];
 		WORK = Realloc(WORK, double, LWORK);
-		status = dsysv('L', m+1, K-1, ZAZ, m+1, IPIV, ZBc, m+1, WORK,
-				LWORK);
+		status = dsysv('L', m+1, K-1, work->ZAZ, m+1, IPIV, work->ZBc, 
+				m+1, WORK, LWORK);
 		if (status != 0)
 			err("[GenSVM Warning]: Received nonzero "
 					"status from dsysv: %i\n", status);
@@ -520,7 +535,7 @@ void gensvm_get_update(struct GenModel *model, struct GenData *data)
 	// convert this back to row-major order
 	for (i=0; i<m+1; i++)
 		for (j=0; j<K-1; j++)
-			ZB[i*(K-1)+j] = ZBc[j*(m+1)+i];
+			work->ZB[i*(K-1)+j] = work->ZBc[j*(m+1)+i];
 
 	// copy the old V to Vbar and the new solution to V
 	for (i=0; i<m+1; i++) {
@@ -528,14 +543,9 @@ void gensvm_get_update(struct GenModel *model, struct GenData *data)
 			matrix_set(model->Vbar, K-1, i, j,
 					matrix_get(model->V, K-1, i, j));
 			matrix_set(model->V, K-1, i, j,
-					matrix_get(ZB, K-1, i, j));
+					matrix_get(work->ZB, K-1, i, j));
 		}
 	}
-
-	free(ZB);
-	free(ZAZ);
-	free(ZBc);
-	free(beta);
 }
 
 /**
diff --git a/tests/src/test_gensvm_base.c b/tests/src/test_gensvm_base.c
index 69ba236..89319c5 100644
--- a/tests/src/test_gensvm_base.c
+++ b/tests/src/test_gensvm_base.c
@@ -73,6 +73,95 @@ char *test_init_free_data_3()
 	return NULL;
 }
 
+char *test_init_free_work()
+{
+	struct GenModel *model = gensvm_init_model();
+	model->n = 10;
+	model->m = 4;
+	model->K = 3;
+
+	struct GenWork *work = gensvm_init_work(model);
+
+	mu_assert(model->n == work->n, "n variable copied incorrectly");
+	mu_assert(model->m == work->m, "m variable copied incorrectly");
+	mu_assert(model->K == work->K, "K variable copied incorrectly");
+
+	mu_assert(work->LZ != NULL, "LZ variable is NULL");
+	mu_assert(work->ZB != NULL, "ZB variable is NULL");
+	mu_assert(work->ZBc != NULL, "ZBc variable is NULL");
+	mu_assert(work->ZAZ != NULL, "ZAZ variable is NULL");
+	mu_assert(work->ZV != NULL, "ZV variable is NULL");
+	mu_assert(work->beta != NULL, "beta variable is NULL");
+
+	gensvm_free_model(model);
+	gensvm_free_work(work);
+
+	return NULL;
+}
+
+void fill_with_noise(double *ptr, int elem)
+{
+	int i;
+	for (i=0; i<elem; i++) {
+		ptr[i] = ((double) rand())/((double) RAND_MAX);
+	}
+}
+
+bool all_elements_zero(double *ptr, int elem)
+{
+	int i;
+	for (i=0; i<elem; i++) {
+		if (ptr[i] != 0)
+			return false;
+	}
+	return true;
+}
+
+char *test_reset_work()
+{
+	struct GenModel *model = gensvm_init_model();
+	int n = 10,
+	    m = 4,
+	    K = 3;
+
+	model->n = n;
+	model->m = m;
+	model->K = K;
+
+	struct GenWork *work = gensvm_init_work(model);
+
+	// start test code //
+	fill_with_noise(work->LZ, n*(m+1));
+	fill_with_noise(work->ZB, (m+1)*(K-1));
+	fill_with_noise(work->ZBc, (m+1)*(K-1));
+	fill_with_noise(work->ZAZ, (m+1)*(m+1));
+	fill_with_noise(work->ZV, n*(K-1));
+	fill_with_noise(work->beta, K-1);
+
+	gensvm_reset_work(work);
+
+	mu_assert(all_elements_zero(work->LZ, n*(m+1)), 
+			"Not all elements of LZ are zero");
+	mu_assert(all_elements_zero(work->ZB, (m+1)*(K-1)),
+			"Not all elements of ZB are zero");
+	mu_assert(all_elements_zero(work->ZBc, (m+1)*(K-1)),
+				"Not all elements of ZBc are zero");
+	mu_assert(all_elements_zero(work->ZAZ, (m+1)*(m+1)),
+				"Not all elements of ZAZ are zero");
+	mu_assert(all_elements_zero(work->ZV, n*(K-1)),
+			"Not all elements of ZV are zero");
+	mu_assert(all_elements_zero(work->beta, K-1),
+			"Not all elements of beta are zero");
+
+	// end test code //
+
+	gensvm_free_model(model);
+	gensvm_free_work(work);
+
+	return NULL;
+}
+
+
 char *all_tests()
 {
 	mu_suite_start();
@@ -80,10 +169,14 @@ char *all_tests()
 	mu_run_test(test_init_free_model_2);
 	mu_run_test(test_allocate_free_model);
 	mu_run_test(test_reallocate_free_model);
+
 	mu_run_test(test_init_free_data_1);
 	mu_run_test(test_init_free_data_2);
 	mu_run_test(test_init_free_data_3);
 
+	mu_run_test(test_init_free_work);
+	mu_run_test(test_reset_work);
+
 	return NULL;
 }
 
diff --git a/tests/src/test_gensvm_optimize.c b/tests/src/test_gensvm_optimize.c
index de7e074..d7326f7 100644
--- a/tests/src/test_gensvm_optimize.c
+++ b/tests/src/test_gensvm_optimize.c
@@ -137,12 +137,17 @@ char *test_gensvm_get_loss_1()
 	int n = 8,
 	    m = 3,
 	    K = 3;
+	model->n = n;
+	model->m = m;
+	model->K = K;
+	struct GenWork *work = gensvm_init_work(model);
 
 	// initialize the data
 	data->n = n;
 	data->K = K;
 	data->m = m;
 
+
 	data->y = Calloc(long, data->n);
 	data->y[0] = 2;
 	data->y[1] = 1;
@@ -188,9 +193,6 @@ char *test_gensvm_get_loss_1()
 	matrix_set(data->Z, data->m+1, 7, 3, 0.0540516162042732);
 
 	// initialize the model
-	model->n = n;
-	model->m = m;
-	model->K = K;
 	model->weight_idx = 1;
 	model->kappa = 0.5;
 	model->p = 1.5;
@@ -211,17 +213,16 @@ char *test_gensvm_get_loss_1()
 	matrix_set(model->V, model->K-1, 3, 1, 0.5912336906588613);
 
 	// start test code //
-	double *ZV = Calloc(double, (data->n)*(data->K-1));
-	double loss = gensvm_get_loss(model, data, ZV);
+	double loss = gensvm_get_loss(model, data, work);
 
 	mu_assert(fabs(loss - 0.903071383013108) < 1e-14,
 			"Incorrect value of the loss");
 
-	free(ZV);
 	// end test code //
 
 	gensvm_free_model(model);
 	gensvm_free_data(data);
+	gensvm_free_work(work);
 
 
 	return NULL;
@@ -235,6 +236,10 @@ char *test_gensvm_get_loss_2()
 	int n = 8,
 	    m = 3,
 	    K = 3;
+	model->n = n;
+	model->m = m;
+	model->K = K;
+	struct GenWork *work = gensvm_init_work(model);
 
 	// initialize the data
 	data->n = n;
@@ -286,9 +291,6 @@ char *test_gensvm_get_loss_2()
 	matrix_set(data->Z, data->m+1, 7, 3, 0.0540516162042732);
 
 	// initialize the model
-	model->n = n;
-	model->m = m;
-	model->K = K;
 	model->weight_idx = 2;
 	model->kappa = 0.5;
 	model->p = 1.5;
@@ -309,17 +311,15 @@ char *test_gensvm_get_loss_2()
 	matrix_set(model->V, model->K-1, 3, 1, 0.5912336906588613);
 
 	// start test code //
-	double *ZV = Calloc(double, (data->n)*(data->K-1));
-	double loss = gensvm_get_loss(model, data, ZV);
+	double loss = gensvm_get_loss(model, data, work);
 
 	mu_assert(fabs(loss - 0.972847045993281) < 1e-14,
 			"Incorrect value of the loss");
-
-	free(ZV);
 	// end test code //
 
 	gensvm_free_model(model);
 	gensvm_free_data(data);
+	gensvm_free_work(work);
 
 	return NULL;
 }
@@ -652,6 +652,11 @@ char *test_gensvm_get_update()
 	    m = 3,
 	    K = 3;
 
+	model->n = n;
+	model->m = m;
+	model->K = K;
+	struct GenWork *work = gensvm_init_work(model);
+
 	// initialize data
 	data->n = n;
 	data->m = m;
@@ -702,9 +707,6 @@ char *test_gensvm_get_update()
 	matrix_set(data->Z, data->m+1, 7, 3, -0.7292593770500276);
 
 	// initialize model
-	model->n = n;
-	model->m = m;
-	model->K = K;
 	model->p = 1.1;
 	model->lambda = 0.123;
 	model->weight_idx = 1;
@@ -727,14 +729,13 @@ char *test_gensvm_get_update()
 	matrix_set(model->V, model->K-1, 3, 1, 0.7134997072555367);
 
 	// start test code //
-	double *ZV = Calloc(double, n*(K-1));
 
 	// these need to be prepared for the update call
-	gensvm_calculate_errors(model, data, ZV);
+	gensvm_calculate_errors(model, data, work->ZV);
 	gensvm_calculate_huber(model);
 
 	// run the actual update call
-	gensvm_get_update(model, data);
+	gensvm_get_update(model, data, work);
 
 	// test values
 	mu_assert(fabs(matrix_get(model->V, model->K-1, 0, 0) -
@@ -761,13 +762,11 @@ char *test_gensvm_get_update()
 	mu_assert(fabs(matrix_get(model->V, model->K-1, 3, 1) -
 				0.4390030236354089) < 1e-14,
 			"Incorrect value of model->V at 3, 1");
-
-	free(ZV);
-
 	// end test code //
 
 	gensvm_free_model(model);
 	gensvm_free_data(data);
+	gensvm_free_work(work);
 
 	return NULL;
 }