extend array3d<> and array4d<> to allow noncontiguous passed-in storage arrays

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

2 files changed, 110 insertions, 52 deletions

diff --git a/src/jtutil/array.cc b/src/jtutil/array.cc
index 05a22e2..1648d6f 100644
--- a/src/jtutil/array.cc
+++ b/src/jtutil/array.cc
@@ -46,9 +46,9 @@ if (stride_i_ == 0)
    then stride_i_ = 1;
 
 // must use unchecked subscripting here since setup isn't done yet
-offset_ = - subscript_unchecked(min_i_);
-assert( subscript_unchecked(min_i()) == 0 );
-max_subscript_ = subscript_unchecked(max_i());
+offset_ = - subscript_unchecked(min_i_);		// RHS uses offset_ = 0
+assert( subscript_unchecked(min_i_) == 0 );
+max_subscript_ = subscript_unchecked(max_i_);
 
 if (we_own_array_)
    then {
@@ -107,9 +107,9 @@ if (stride_i_ == 0)
    then stride_i_ = N_j();
 
 // must use unchecked subscripting here since setup isn't done yet
-offset_ = - subscript_unchecked(min_i_, min_j_);
-assert( subscript_unchecked(min_i(), min_j()) == 0 );
-max_subscript_ = subscript_unchecked(max_i(), max_j());
+offset_ = - subscript_unchecked(min_i_,min_j_);		// RHS uses offset_ = 0
+assert( subscript_unchecked(min_i_,min_j_) == 0 );
+max_subscript_ = subscript_unchecked(max_i_,max_j_);
 
 if (we_own_array_)
    then {
@@ -157,26 +157,48 @@ template <typename fp>
 array3d<fp>::array3d(int min_i_in, int max_i_in,
 		     int min_j_in, int max_j_in,
 		     int min_k_in, int max_k_in,
-		     fp *array_in = NULL)
-	: min_i_(min_i_in), max_i_(max_i_in),
+		     fp *array_in = NULL,
+		     int stride_i_in = 0, int stride_j_in = 0,
+		     int stride_k_in = 0)
+	: array_(array_in),
+	  offset_(0),		// temp value, changed below
+	  stride_i_(stride_i_in), stride_j_(stride_j_in),
+	  stride_k_(stride_k_in),
+	  min_i_(min_i_in), max_i_(max_i_in),
 	  min_j_(min_j_in), max_j_(max_j_in),
 	  min_k_(min_k_in), max_k_(max_k_in),
 	  we_own_array_(array_in == NULL)
 {
-stride_j_ = N_k();
-stride_i_ = N_j() * stride_j_;
-N_array_ = N_i() * stride_i_;
-array_ = we_own_array_ ? new fp[N_array_] : array_in;
+if (stride_k_ == 0)
+   then stride_k_ = 1;
+if (stride_j_ == 0)
+   then stride_j_ = N_k();
+if (stride_i_ == 0)
+   then stride_i_ = N_j()*N_k();
 
-// need to explicitly initialize -- new[] *doesn't* do this automagically
-	for (int i = 0 ; i < N_array() ; ++i)
-	{
-	array_[i] = fp(0);
+// must use unchecked subscripting here since setup isn't done yet
+offset_ = - subscript_unchecked(min_i_,min_j_,min_k_);	// RHS uses offset_ = 0
+assert( subscript_unchecked(min_i_,min_j_,min_k_) == 0 );
+max_subscript_ = subscript_unchecked(max_i_,max_j_,max_k_);
+
+if (we_own_array_)
+   then {
+	// allocate it
+	const int N_allocate = N_i() * N_j() * N_k();
+	array_ = new fp[N_allocate];
 	}
 
-offset_ = 0;	// temp value, needed for following  subscript()  call
-offset_ = - subscript_unchecked(min_i_, min_j_, min_k_);
-		// must use unchecked form here since setup isn't done yet
+// explicitly initialize array (new[] *doesn't* do this automagically)
+	for (int i = min_i() ; i <= max_i() ; ++i)
+	{
+	for (int j = min_j() ; j <= max_j() ; ++j)
+	{
+	for (int k = min_k() ; k <= max_k() ; ++k)
+	{
+	operator()(i,j,k) = fp(0);
+	}
+	}
+	}
 }
 	  }	// namespace jtutil::
 
@@ -209,28 +231,55 @@ array4d<fp>::array4d(int min_i_in, int max_i_in,
 		     int min_j_in, int max_j_in,
 		     int min_k_in, int max_k_in,
 		     int min_l_in, int max_l_in,
-		     fp *array_in = NULL)
-	: min_i_(min_i_in), max_i_(max_i_in),
+		     fp *array_in = NULL,
+		     int stride_i_in = 0, int stride_j_in = 0,
+		     int stride_k_in = 0, int stride_l_in = 0)
+	: array_(array_in),
+	  offset_(0),		// temp value, changed below
+	  stride_i_(stride_i_in), stride_j_(stride_j_in),
+	  stride_k_(stride_k_in), stride_l_(stride_l_in),
+	  min_i_(min_i_in), max_i_(max_i_in),
 	  min_j_(min_j_in), max_j_(max_j_in),
 	  min_k_(min_k_in), max_k_(max_k_in),
 	  min_l_(min_l_in), max_l_(max_l_in),
 	  we_own_array_(array_in == NULL)
 {
-stride_k_ = N_l();
-stride_j_ = N_k() * stride_k_;
-stride_i_ = N_j() * stride_j_;
-N_array_ = N_i() * stride_i_;
-array_ = we_own_array_ ? new fp[N_array_] : array_in;
-
-// need to explicitly initialize -- new[] *doesn't* do this automagically
-	for (int i = 0 ; i < N_array() ; ++i)
-	{
-	array_[i] = fp(0);
+if (stride_l_ == 0)
+   then stride_l_ = 1;
+if (stride_k_ == 0)
+   then stride_k_ = N_l();
+if (stride_j_ == 0)
+   then stride_j_ = N_k()*N_l();
+if (stride_i_ == 0)
+   then stride_i_ = N_j()*N_k()*N_l();
+
+// must use unchecked subscripting here since setup isn't done yet
+offset_ = - subscript_unchecked(min_i_,min_j_,
+				min_k_,min_l_);		// RHS uses offset_ = 0
+assert( subscript_unchecked(min_i_,min_j_,min_k_,min_l_) == 0 );
+max_subscript_ = subscript_unchecked(max_i_,max_j_,max_k_,max_l_);
+
+if (we_own_array_)
+   then {
+	// allocate it
+	const int N_allocate = N_i() * N_j() * N_k() * N_l();
+	array_ = new fp[N_allocate];
 	}
 
-offset_ = 0;	// temp value, needed for following  subscript()  call
-offset_ = - subscript_unchecked(min_i_, min_j_, min_k_, min_l_);
-		// must use unchecked form here since setup isn't done yet
+// explicitly initialize array (new[] *doesn't* do this automagically)
+	for (int i = min_i() ; i <= max_i() ; ++i)
+	{
+	for (int j = min_j() ; j <= max_j() ; ++j)
+	{
+	for (int k = min_k() ; k <= max_k() ; ++k)
+	{
+	for (int l = min_l() ; l <= max_l() ; ++l)
+	{
+	operator()(i,j,k,l) = fp(0);
+	}
+	}
+	}
+	}
 }
 	  }	// namespace jtutil::
 
diff --git a/src/jtutil/array.hh b/src/jtutil/array.hh
index 492c1d9..0185b8d 100644
--- a/src/jtutil/array.hh
+++ b/src/jtutil/array.hh
@@ -261,7 +261,7 @@ public:
 	// FIXME: should we also provide the reverse mapping, i.e.
 	//        subscript --> (i,j,k) ?
 	int subscript_unchecked(int i, int j, int k) const
-		{ return offset_ + stride_i_*i + stride_j_*j + k; }
+		{ return offset_ + stride_i_*i + stride_j_*j + stride_k_*k; }
 	int subscript(int i, int j, int k) const
 		{
 		assert( is_valid_subscript(i,j,k) );
@@ -270,7 +270,7 @@ public:
 		//	source line, so an assert() failure message can
 		//	pinpoint *which* index is bad
 		assert(posn >= 0);
-		assert(posn <= N_array_);
+		assert(posn <= max_subscript_);
 		return posn;
 		}
 
@@ -284,16 +284,20 @@ public:
 
 	// get access to internal 0-origin 1D storage array
 	// (low-level, dangerous, use with caution!)
-	int N_array() const { return N_array_; }
+	// ... semantics of N_array() may not be what you want
+	//     if strides specify noncontiguous storage
+	int N_array() const { return max_subscript_+stride_k_; }
 	fp* get_array() const { return array_; }
 
 	// constructor, destructor
-	// constructor initializes all array elements to fp(0.0)
+	// ... constructor initializes all array elements to fp(0.0)
+	// ... omitted strides default to C storage order
 	array3d(int min_i_in, int max_i_in,
 		int min_j_in, int max_j_in,
 		int min_k_in, int max_k_in,
-		fp *array_in = NULL);	// caller-provided storage array
+		fp *array_in = NULL,	// caller-provided storage array
 					// if non-NULL
+		int stride_i_in = 0, int stride_j_in = 0, int stride_k_in = 0);
 	~array3d();
 
 private:
@@ -304,7 +308,7 @@ private:
 	array3d<fp>& operator=(const array3d<fp>& rhs);
 
 private:
-	// note we declare the array pointer first in the class
+	// n.b. we declare the array pointer first in the class
 	// ==> it's probably at 0 offset
 	// ==> we may get slightly faster array access
 	fp* array_;		// --> new-allocated 1D storage array
@@ -312,13 +316,13 @@ private:
 	// subscripting info
 	// n.b. put this next in class so it should be in the same
 	//	cpu cache line as  array_  ==> faster array access
-	int offset_, stride_i_, stride_j_;
+	int offset_, stride_i_, stride_j_, stride_k_;
 
 	// min/max array bounds
-	int N_array_;
 	const int min_i_, max_i_;
 	const int min_j_, max_j_;
 	const int min_k_, max_k_;
+	int max_subscript_;
 
 	// n.b. put this at end of class since performance doesn't matter
 	bool we_own_array_;	// true ==> array_ --> new[] array which we own
@@ -351,8 +355,7 @@ public:
 	bool is_valid_j(int j) const { return (j >= min_j_) && (j <= max_j_); }
 	bool is_valid_k(int k) const { return (k >= min_k_) && (k <= max_k_); }
 	bool is_valid_l(int l) const { return (l >= min_l_) && (l <= max_l_); }
-	bool is_valid_subscript(int i, int j, int k, int l)
-		const
+	bool is_valid_subscript(int i, int j, int k, int l) const
 		{
 		return    is_valid_i(i) && is_valid_j(j)
 		       && is_valid_k(k) && is_valid_l(l);
@@ -364,7 +367,8 @@ public:
 	//        subscript --> (i,j,k,l) ?
 	int subscript_unchecked(int i, int j, int k, int l) const
 		{
-		return offset_ + stride_i_*i + stride_j_*j + stride_k_*k + l;
+		return offset_ + stride_i_*i + stride_j_*j
+			       + stride_k_*k + stride_l_*l;
 		}
 	int subscript(int i, int j, int k, int l) const
 		{
@@ -374,7 +378,7 @@ public:
 		//	source line, so an assert() failure message can
 		//	pinpoint *which* index is bad
 		assert(posn >= 0);
-		assert(posn <= N_array_);
+		assert(posn <= max_subscript_);
 		return posn;
 		}
 
@@ -388,17 +392,22 @@ public:
 
 	// get access to internal 0-origin 1D storage array
 	// (low-level, dangerous, use with caution!)
-	int N_array() const { return N_array_; }
+	// ... semantics of N_array() may not be what you want
+	//     if strides specify noncontiguous storage
+	int N_array() const { return max_subscript_+stride_l_; }
 	fp* get_array() const { return array_; }
 
 	// constructor, destructor
-	// constructor initializes all array elements to fp(0.0)
+	// ... constructor initializes all array elements to fp(0.0)
+	// ... omitted strides default to C storage order
 	array4d(int min_i_in, int max_i_in,
 		int min_j_in, int max_j_in,
 		int min_k_in, int max_k_in,
 		int min_l_in, int max_l_in,
-		fp *array_in = NULL);	// caller-provided storage array
+		fp *array_in = NULL,	// caller-provided storage array
 					// if non-NULL
+		int stride_i_in = 0, int stride_j_in = 0,
+		int stride_k_in = 0, int stride_l_in = 0);
 	~array4d();
 
 private:
@@ -409,7 +418,7 @@ private:
 	array4d<fp>& operator=(const array4d<fp>& rhs);
 
 private:
-	// note we declare the array pointer first in the class
+	// n.b. we declare the array pointer first in the class
 	// ==> it's probably at 0 offset
 	// ==> we may get slightly faster array access
 	fp* array_;		// --> new-allocated 1D storage array
@@ -417,14 +426,14 @@ private:
 	// subscripting info
 	// n.b. put this next in class so it should be in the same
 	//	cpu cache line as  array_  ==> faster array access
-	int offset_, stride_i_, stride_j_, stride_k_;
+	int offset_, stride_i_, stride_j_, stride_k_, stride_l_;
 
 	// min/max array bounds
-	int N_array_;
 	const int min_i_, max_i_;
 	const int min_j_, max_j_;
 	const int min_k_, max_k_;
 	const int min_l_, max_l_;
+	int max_subscript_;
 
 	// n.b. put this at end of class since performance doesn't matter
 	bool we_own_array_;	// true ==> array_ --> new[] array which we own