finish upgrading test driver to test all array classes with

both their own and noncontiguous/passed-in storage arrays


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@385 f88db872-0e4f-0410-b76b-b9085cfa78c5

1 files changed, 308 insertions, 222 deletions

diff --git a/src/jtutil/test_array.cc b/src/jtutil/test_array.cc
index c34c0cc..e9fd235 100644
--- a/src/jtutil/test_array.cc
+++ b/src/jtutil/test_array.cc
@@ -19,11 +19,10 @@ using jtutil::array4d;
 //
 // prototypes for non-class functions defined in this file
 //
-void test_array1d_double(bool print_flag);
-void test_array2d_float(bool print_flag);
-void test_array2d_float(bool print_flag);
-void test_array3d_double(bool print_flag);
-void test_array4d_double(bool print_flag);
+double test_array1d_double(bool fancy_storage, bool print_flag);
+double test_array2d_float(bool fancy_storage, bool print_flag);
+double test_array3d_double(bool fancy_storage, bool print_flag);
+double test_array4d_double(bool fancy_storage, bool print_flag);
 
 //******************************************************************************
 
@@ -36,236 +35,256 @@ void test_array4d_double(bool print_flag);
 int main(int argc, const char* const argv[])
 {
 const bool print_flag = (argc > 1);
-test_array1d_double(print_flag);
-test_array2d_float(print_flag);
-test_array3d_double(print_flag);
-test_array4d_double(print_flag);
+
+double sum_RMS_diff = 0.0;
+
+sum_RMS_diff += test_array1d_double(false, print_flag);
+sum_RMS_diff += test_array1d_double(true,  print_flag);
+
+sum_RMS_diff += test_array2d_float(false, print_flag);
+sum_RMS_diff += test_array2d_float(true,  print_flag);
+
+sum_RMS_diff += test_array3d_double(false, print_flag);
+sum_RMS_diff += test_array3d_double(true,  print_flag);
+
+sum_RMS_diff += test_array4d_double(false, print_flag);
+sum_RMS_diff += test_array4d_double(true,  print_flag);
+
+printf("\n");
+printf("grand total RMS_diff = %g\n", sum_RMS_diff);
+if (sum_RMS_diff == 0)
+   then {
+	printf("==> all ok!\n");
+	return 0;
+	}
+   else {
+	printf("==> ***** test(s) failed *****\n");
+	return 1;
+	}
 }
 
 //******************************************************************************
 
 //
-// This function tests the  array1d<double>  class twice, once with
-// its own storage array, and once with a noncontiguous storage array
-// which array1d<double> doesn't itself own.
+// This function tests the  array1d<double>  class.  It returns the
+// "RMS difference" error measure (should be 0.0 if all is working ok).
 //
-void test_array1d_double(bool print_flag)
+double test_array1d_double(bool fancy_storage, bool print_flag)
 {
 const int min_i = 11;		const int max_i = 13;
 const double delta = 1.0/7.0;
 const double change_value = 144.025036146;
 
-	for (int test = 0 ; test < 2 ; ++test)
+double *array;
+array1d<double> *Aptr;
+int stride;
+
+// create the array
+printf("\n");
+printf("testing array1d<double>");
+if (fancy_storage)
+   then {
+	printf(" with noncontiguous non-owned storage array...\n");
+	const int N_i = jtutil::how_many_in_range(min_i, max_i);
+	stride = 4;
+	const int N = N_i*stride;
+	printf(
+   "allocating %d * stride=%d = %d underlying storage array...\n",
+	       N_i, stride, N);
+	array = new double[N];
+	printf("constructing...\n");
+	Aptr = new array1d<double>(min_i, max_i,
+				   array,
+				   stride);
+	}
+   else {
+	printf(" with its own storage array...\n");
+	Aptr = new array1d<double>(min_i, max_i);
+	stride = 1;
+	}
+array1d<double>& A = *Aptr;
+
+assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
+assert(A.N_i() == jtutil::how_many_in_range(min_i, max_i));
+
+printf("assigning values...\n");
+#define VALUE_1D(ix)	\
+		double(1000.0 + ix + delta)
+	for (int i = min_i ; i <= max_i ; ++i)
 	{
-	double *array;
-	array1d<double> *Aptr;
-	int stride;
-
-	// create the array
-	printf("\ntesting array1d<double>");
-	if (test == 0)
-	   then {
-		printf(" with its own storage array...\n");
-		Aptr = new array1d<double>(min_i,max_i);
-		stride = 1;
-		}
-	   else {
-		printf(" with noncontiguous non-owned storage array...\n");
-		const int N_i = jtutil::how_many_in_range(min_i, max_i);
-		stride = 4;
-		const int N = N_i*stride;
-		printf(
-		   "allocating N_i=%d * stride=%d = %d underlying storage array...\n",
-		       N_i, stride, N);
-		array = new double[N];
-		printf("constructing...\n");
-		Aptr = new array1d<double>(min_i,max_i,
-					   array,
-					   stride);
-		}
-	array1d<double>& A = *Aptr;
+	if (print_flag)
+	   then printf("   assigning %d\n", i);
+	A(i) = VALUE_1D(i);
+	}
 
-	assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
-	assert(A.N_i() == jtutil::how_many_in_range(min_i,max_i));
+printf("checking N-D vs 1-D addresses...\n");
+const double* array_ptr = A.get_array();
+int sub = 0;
+	for (int i = min_i ; i <= max_i ; ++i)
+	{
+	if (print_flag)
+	   then printf("   checking %d\n", i);
+	assert( &A(i) == &array_ptr[sub] );
+	sub += stride;
+	}
+assert( A.N_array() == sub );
 
-	printf("assigning values...\n");
-	#define VALUE_1D(ix)	\
-			double(1000.0 + ix + delta)
-		for (int i = min_i ; i <= max_i ; ++i)
-		{
-		if (print_flag)
-		   then printf("   assigning %d\n", i);
-		A(i) = VALUE_1D(i);
-		}
+printf("checking what changes when we assign an array element...\n");
+double sumsq_of_diff = 0.0;
+	for (int i = min_i ; i <= max_i ; ++i)
+	{
+	if (print_flag)
+	   then printf("   changing %d\n", i);
+	double save_a = A(i);
+	A(i) = change_value;
 
-	printf("checking N-D vs 1-D addresses...\n");
-	int sub = 0;
-	const double* array_ptr = A.get_array();
-		for (int i = min_i ; i <= max_i ; ++i)
+		for (int ii = min_i ; ii <= max_i ; ++ii)
 		{
 		if (print_flag)
-		   then printf("   checking %d\n", i);
-		assert( &A(i) == &array_ptr[sub] );
-		sub += stride;
-		}
-	assert( A.N_array() == sub );
+		   then printf("      checking %d\n", ii);
 
-	printf("checking what changes when we assign an array element...\n");
-	double sumsq_of_diff = 0.0;
-		for (int i = min_i ; i <= max_i ; ++i)
-		{
-		if (print_flag)
-		   then printf("   changing %d\n", i);
-		double save_a = A(i);
-		A(i) = change_value;
-
-			for (int ii = min_i ; ii <= max_i ; ++ii)
-			{
-			if (print_flag)
-			   then printf("      checking %d\n", ii);
-
-			double should_be = (ii == i)
-					   ? change_value : VALUE_1D(ii);
-			assert( fuzzy<double>::EQ(A(ii), should_be) );
-
-			double diff = A(ii) - should_be;
-			sumsq_of_diff += diff*diff;
-			}
-		A(i) = save_a;
+		double should_be = (ii == i)
+				   ? change_value : VALUE_1D(ii);
+		assert( fuzzy<double>::EQ(A(ii), should_be) );
+
+		double diff = A(ii) - should_be;
+		sumsq_of_diff += diff*diff;
 		}
-	double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
-	printf("==> everything looks ok (RMS_diff=%.3g)\n", double(RMS_diff));
+	A(i) = save_a;
+	}
+double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
+printf("==> everything looks ok (RMS_diff=%.3g)\n", double(RMS_diff));
 
-	delete Aptr;
+delete Aptr;
 
-	if (test == 1)
-	   then {
-		printf("deleting underlying storage array...\n");
-		delete[] array;
-		}
+if (fancy_storage)
+   then {
+	printf("deleting underlying storage array...\n");
+	delete[] array;
 	}
+
+return RMS_diff;
 }
 
 //******************************************************************************
 
 //
-// This function tests the  array2d<float>  class twice, once with
-// its own storage array, and once with a noncontiguous storage array
-// which array2d<float> doesn't itself own.
+// This function tests the  array2d<float>  class.  It returns the
+// "RMS difference" error measure (should be 0.0 if all is working ok).
 //
-void test_array2d_float(bool print_flag)
+double test_array2d_float(bool fancy_storage, bool print_flag)
 {
 const int min_i = 11;		const int max_i = 13;
 const int min_j = 21;		const int max_j = 24;
 const float delta = 1.0/7.0;
 const float change_value = 144.025036146;
 
-	for (int test = 0 ; test < 2 ; ++test)
+float *array;
+array2d<float> *Aptr;
+int stride;
+
+// create the array
+printf("\n");
+printf("testing array2d<float>");
+if (fancy_storage)
+   then {
+	printf(" with noncontiguous non-owned storage array...\n");
+	const int N_i = jtutil::how_many_in_range(min_i, max_i);
+	const int N_j = jtutil::how_many_in_range(min_j, max_j);
+	stride = 3;
+	const int N = N_i*N_j*stride;
+	printf(
+   "allocating %d*%d * stride=%d = %d underlying storage array...\n",
+	       N_i,N_j, stride, N);
+	array = new float[N];
+	printf("constructing...\n");
+	Aptr = new array2d<float>(min_i, max_i,
+				  min_j, max_j,
+				  array,
+				  N_j*stride, stride);
+	}
+   else {
+	printf(" with its own storage array...\n");
+	Aptr = new array2d<float>(min_i, max_i,
+				  min_j, max_j);
+	stride = 1;
+	}
+array2d<float>& A = *Aptr;
+
+assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
+assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
+assert(A.N_i() == jtutil::how_many_in_range(min_i, max_i));
+assert(A.N_j() == jtutil::how_many_in_range(min_j, max_j));
+
+printf("assigning values...\n");
+#define VALUE_2D(ix,jx)	\
+		float(1000.0 + 10.0*jx + ix + delta)
+	for (int i = min_i ; i <= max_i ; ++i)
 	{
-	float *array;
-	array2d<float> *Aptr;
-	int stride;
-
-	// create the array
-	printf("\ntesting array2d<float>");
-	if (test == 0)
-	   then {
-		printf(" with its own storage array...\n");
-		Aptr = new array2d<float>(min_i,max_i,
-					  min_j,max_j);
-		stride = 1;
-		}
-	   else {
-		printf(" with noncontiguous non-owned storage array...\n");
-		const int N_i = jtutil::how_many_in_range(min_i, max_i);
-		const int N_j = jtutil::how_many_in_range(min_j, max_j);
-		stride = 3;
-		const int N = N_i*N_j*stride;
-		printf(
-		   "allocating N_i=%d * N_j=%d * stride=%d = %d underlying storage array...\n",
-		       N_i, N_j, stride, N);
-		array = new float[N];
-		printf("constructing...\n");
-		Aptr = new array2d<float>(min_i,max_i,
-					  min_j,max_j,
-					  array,
-					  N_j*stride, stride);
-		}
-	array2d<float>& A = *Aptr;
+	for (int j = min_j ; j <= max_j ; ++j)
+	{
+	if (print_flag)
+	   then printf("   assigning %d %d\n", i, j);
+	A(i,j) = VALUE_2D(i,j);
+	}
+	}
 
-	assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
-	assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
-	assert(A.N_i() == jtutil::how_many_in_range(min_i,max_i));
-	assert(A.N_j() == jtutil::how_many_in_range(min_j,max_j));
+printf("checking N-D vs 1-D addresses...\n");
+const float* array_ptr = A.get_array();
+int sub = 0;
+	for (int i = min_i ; i <= max_i ; ++i)
+	{
+	for (int j = min_j ; j <= max_j ; ++j)
+	{
+	if (print_flag)
+	   then printf("   checking %d %d\n", i, j);
+	assert( &A(i,j) == &array_ptr[sub] );
+	sub += stride;
+	}
+	}
+assert( A.N_array() == sub );
 
-	printf("assigning values...\n");
-	#define VALUE_2D(ix,jx)	\
-			float(1000.0 + 10.0*jx + ix + delta)
-		for (int i = min_i ; i <= max_i ; ++i)
-		{
-		for (int j = min_j ; j <= max_j ; ++j)
-		{
-		if (print_flag)
-		   then printf("   assigning %d %d\n", i, j);
-		A(i,j) = VALUE_2D(i,j);
-		}
-		}
+printf("checking what changes when we assign an array element...\n");
+double sumsq_of_diff = 0.0;
+	for (int i = min_i ; i <= max_i ; ++i)
+	{
+	for (int j = min_j ; j <= max_j ; ++j)
+	{
+	if (print_flag)
+	   then printf("   changing %d %d\n", i, j);
+	float save_a = A(i,j);
+	A(i,j) = change_value;
 
-	printf("checking N-D vs 1-D addresses...\n");
-	int sub = 0;
-	const float* array_ptr = A.get_array();
-		for (int i = min_i ; i <= max_i ; ++i)
+		for (int ii = min_i ; ii <= max_i ; ++ii)
 		{
-		for (int j = min_j ; j <= max_j ; ++j)
+		for (int jj = min_j ; jj <= max_j ; ++jj)
 		{
 		if (print_flag)
-		   then printf("   checking %d %d\n", i, j);
-		assert( &A(i,j) == &array_ptr[sub] );
-		sub += stride;
-		}
-		}
-	assert( A.N_array() == sub );
+		   then printf("      checking %d %d\n", ii, jj);
 
-	printf("checking what changes when we assign an array element...\n");
-	double sumsq_of_diff = 0.0;
-		for (int i = min_i ; i <= max_i ; ++i)
-		{
-		for (int j = min_j ; j <= max_j ; ++j)
-		{
-		if (print_flag)
-		   then printf("   changing %d %d\n", i, j);
-		float save_a = A(i,j);
-		A(i,j) = change_value;
-
-			for (int ii = min_i ; ii <= max_i ; ++ii)
-			{
-			for (int jj = min_j ; jj <= max_j ; ++jj)
-			{
-			if (print_flag)
-			   then printf("      checking %d %d\n", ii, jj);
-
-			float should_be = ((ii == i) && (jj == j))
-					  ? change_value : VALUE_2D(ii,jj);
-			assert( fuzzy<float>::EQ(A(ii,jj), should_be) );
-
-			double diff = A(ii,jj) - should_be;
-			sumsq_of_diff += diff*diff;
-			}
-			}
-		A(i,j) = save_a;
+		float should_be = ((ii == i) && (jj == j))
+				  ? change_value : VALUE_2D(ii,jj);
+		assert( fuzzy<float>::EQ(A(ii,jj), should_be) );
+
+		double diff = A(ii,jj) - should_be;
+		sumsq_of_diff += diff*diff;
 		}
 		}
-	float RMS_diff = sqrt(sumsq_of_diff / A.N_array());
-	printf("==> everything looks ok (RMS_diff=%.3g)\n", double(RMS_diff));
+	A(i,j) = save_a;
+	}
+	}
+double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
+printf("==> everything looks ok (RMS_diff=%.3g)\n", double(RMS_diff));
 
-	delete Aptr;
+delete Aptr;
 
-	if (test == 1)
-	   then {
-		printf("deleting underlying storage array...\n");
-		delete[] array;
-		}
+if (fancy_storage)
+   then {
+	printf("deleting underlying storage array...\n");
+	delete[] array;
 	}
+
+return RMS_diff;
 }
 
 //******************************************************************************
@@ -273,7 +292,7 @@ const float change_value = 144.025036146;
 //
 // This function tests the  array3d<double>  class.
 //
-void test_array3d_double(bool print_flag)
+double test_array3d_double(bool fancy_storage, bool print_flag)
 {
 const int min_i = 11;		const int max_i = 13;
 const int min_j = 21;		const int max_j = 24;
@@ -281,26 +300,47 @@ const int min_k = 31;		const int max_k = 35;
 const double delta = 1.0/7.0;
 const double change_value = 144.025036146;
 
-printf("\ntesting array3d<double> with passed-in array...\n");
-
-const int N_i = jtutil::how_many_in_range(min_i,max_i);
-const int N_j = jtutil::how_many_in_range(min_j,max_j);
-const int N_k = jtutil::how_many_in_range(min_k,max_k);
-const int N = N_i*N_j*N_k;
-printf("allocating %d*%d*%d = %d underlying storage array...\n",
-       N_i,N_j,N_k, N);
-double *array = new double[N];
+double *array;
+array3d<double> *Aptr;
+int stride;
+
+// create the array
+printf("\n");
+printf("testing array3d<double>");
+if (fancy_storage)
+   then {
+	printf(" with noncontiguous non-owned storage array...\n");
+	const int N_i = jtutil::how_many_in_range(min_i, max_i);
+	const int N_j = jtutil::how_many_in_range(min_j, max_j);
+	const int N_k = jtutil::how_many_in_range(min_k, max_k);
+	stride = 4;
+	const int N = N_i*N_j*N_k*stride;
+	printf(
+   "allocating %d*%d*%d * stride=%d = %d   underlying storage array...\n",
+	       N_i,N_j,N_k, stride, N);
+	array = new double[N];
+	printf("constructing...\n");
+	Aptr = new array3d<double>(min_i, max_i,
+				   min_j, max_j,
+				   min_k, max_k,
+				   array,
+				   N_j*N_k*stride,N_k*stride, stride);
+	}
+   else {
+	printf(" with its own storage array...\n");
+	Aptr = new array3d<double>(min_i, max_i,
+				   min_j, max_j,
+				   min_k, max_k);
+	stride = 1;
+	}
+array3d<double>& A = *Aptr;
 
-  // local scope for array
-  {
-printf("constructing array3d<double> object...\n");
-array3d<double> A(min_i,max_i, min_j,max_j, min_k,max_k, array);
 assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
 assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
 assert(A.min_k() == min_k);	assert(A.max_k() == max_k);
-assert(A.N_i() == jtutil::how_many_in_range(min_i,max_i));
-assert(A.N_j() == jtutil::how_many_in_range(min_j,max_j));
-assert(A.N_k() == jtutil::how_many_in_range(min_k,max_k));
+assert(A.N_i() == jtutil::how_many_in_range(min_i, max_i));
+assert(A.N_j() == jtutil::how_many_in_range(min_j, max_j));
+assert(A.N_k() == jtutil::how_many_in_range(min_k, max_k));
 
 printf("assigning values...\n");
 #define VALUE_3D(ix,jx,kx)	\
@@ -319,8 +359,8 @@ printf("assigning values...\n");
 	}
 
 printf("checking N-D vs 1-D addresses...\n");
-int sub = 0;
 const double* array_ptr = A.get_array();
+int sub = 0;
 	for (int i = min_i ; i <= max_i ; ++i)
 	{
 	for (int j = min_j ; j <= max_j ; ++j)
@@ -329,7 +369,8 @@ const double* array_ptr = A.get_array();
 	{
 	if (print_flag)
 	   then printf("   checking %d %d %d\n", i, j, k);
-	assert( &A(i,j,k) == &array_ptr[sub++] );
+	assert( &A(i,j,k) == &array_ptr[sub] );
+	sub += stride;
 	}
 	}
 	}
@@ -373,11 +414,13 @@ double sumsq_of_diff = 0.0;
 double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
 printf("==> everything looks ok (RMS_diff=%.3g)\n", RMS_diff);
 
-printf("exiting local scope ==> destroying array3d object...\n");
-  }
+if (fancy_storage)
+   then {
+	printf("deleting underlying storage array...\n");
+	delete[] array;
+	}
 
-printf("deleting underlying storage array...\n");
-delete[] array;
+return RMS_diff;
 }
 
 //******************************************************************************
@@ -385,7 +428,7 @@ delete[] array;
 //
 // This function tests the  array4d<double>  class.
 //
-void test_array4d_double(bool print_flag)
+double test_array4d_double(bool fancy_storage, bool print_flag)
 {
 const int min_i = 11;		const int max_i = 13;
 const int min_j = 21;		const int max_j = 24;
@@ -394,18 +437,53 @@ const int min_l = 41;		const int max_l = 46;
 const double delta = 1.0/7.0;
 const double change_value = 144.025036146;
 
-printf("\ntesting array4d<double>...\n");
+double *array;
+array4d<double> *Aptr;
+int stride;
+
+// create the array
+printf("\n");
+printf("testing array4d<double>");
+if (fancy_storage)
+   then {
+	printf(" with noncontiguous non-owned storage array...\n");
+	const int N_i = jtutil::how_many_in_range(min_i, max_i);
+	const int N_j = jtutil::how_many_in_range(min_j, max_j);
+	const int N_k = jtutil::how_many_in_range(min_k, max_k);
+	const int N_l = jtutil::how_many_in_range(min_l, max_l);
+	stride = 3;
+	const int N = N_i*N_j*N_k*N_l*stride;
+	printf(
+   "allocating %d*%d*%d*%d * stride=%d = %d   underlying storage array...\n",
+	       N_i,N_j,N_k,N_l, stride, N);
+	array = new double[N];
+	printf("constructing...\n");
+	Aptr = new array4d<double>(min_i, max_i,
+				   min_j, max_j,
+				   min_k, max_k,
+				   min_l, max_l,
+				   array,
+				   N_j*N_k*N_l*stride, N_k*N_l*stride,
+				   N_l*stride, stride);
+	}
+   else {
+	printf(" with its own storage array...\n");
+	Aptr = new array4d<double>(min_i, max_i,
+				   min_j, max_j,
+				   min_k, max_k,
+				   min_l, max_l);
+	stride = 1;
+	}
+array4d<double>& A = *Aptr;
 
-printf("constructing...\n");
-array4d<double> A(min_i,max_i, min_j,max_j, min_k,max_k, min_l,max_l);
 assert(A.min_i() == min_i);	assert(A.max_i() == max_i);
 assert(A.min_j() == min_j);	assert(A.max_j() == max_j);
 assert(A.min_k() == min_k);	assert(A.max_k() == max_k);
 assert(A.min_l() == min_l);	assert(A.max_l() == max_l);
-assert(A.N_i() == jtutil::how_many_in_range(min_i,max_i));
-assert(A.N_j() == jtutil::how_many_in_range(min_j,max_j));
-assert(A.N_k() == jtutil::how_many_in_range(min_k,max_k));
-assert(A.N_l() == jtutil::how_many_in_range(min_l,max_l));
+assert(A.N_i() == jtutil::how_many_in_range(min_i, max_i));
+assert(A.N_j() == jtutil::how_many_in_range(min_j, max_j));
+assert(A.N_k() == jtutil::how_many_in_range(min_k, max_k));
+assert(A.N_l() == jtutil::how_many_in_range(min_l, max_l));
 
 printf("assigning values...\n");
 #define VALUE_4D(ix,jx,kx,lx)	\
@@ -427,8 +505,8 @@ printf("assigning values...\n");
 	}
 
 printf("checking N-D vs 1-D addresses...\n");
-int sub = 0;
 const double* array_ptr = A.get_array();
+int sub = 0;
 	for (int i = min_i ; i <= max_i ; ++i)
 	{
 	for (int j = min_j ; j <= max_j ; ++j)
@@ -439,7 +517,8 @@ const double* array_ptr = A.get_array();
 	{
 	if (print_flag)
 	   then printf("   checking %d %d %d %d\n", i, j, k, l);
-	assert( &A(i,j,k,l) == &array_ptr[sub++] );
+	assert( &A(i,j,k,l) == &array_ptr[sub] );
+	sub += stride;
 	}
 	}
 	}
@@ -489,6 +568,13 @@ double sumsq_of_diff = 0.0;
 	}
 	}
 double RMS_diff = sqrt(sumsq_of_diff / A.N_array());
-
 printf("==> everything looks ok (RMS_diff=%.3g)\n", RMS_diff);
+
+if (fancy_storage)
+   then {
+	printf("deleting underlying storage array...\n");
+	delete[] array;
+	}
+
+return RMS_diff;
 }