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

2 files changed, 114 insertions, 22 deletions

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;
 }