Large scale scientific datasets are generally modeled as k-dimensional arrays, since this model is amenable to the form of analyses and visualization of the scientific phenomenon often investigated. In recent years, organizations have adopted the use of on-line analytical processing (OLAP), methods and statistical analyses to make strategic business decisions using enterprise data that are modeled as multi-dimensional arrays as well. In both of these domains, the datasets have the propensity to gradually grow, reaching orders of terabytes. However, the storage schemes used for these arrays correspond to those where the array elements are allocated in a sequence of consecutive locations according to an ordering of array mapping functions that map k-dimensional indices one-to-one onto the linear locations. Such schemes limit the degree of extendibility of the array to one dimension only. We present a method of allocating storage for the elements of a dense multidimensional extendible array such that the bounds on the indices of the respective dimensions can be arbitrarily extended without reorganizing previously allocated elements. We give a mapping function F ∗(), and its inverse F ∗−1(), for computing the linear address of an array element given its k-dimensional index. The technique adopts the mapping function, for realizing an extendible array with arbitrary extendibility in main memory, to implement such array files. We show how the extendible array file implementation gives an efficient storage scheme for both scientific and OLAP multi-dimensional datasets that are allowed to incrementally grow without incurring the prohibitive costs of reorganizations. Datasets used in large scale scientific applications, are generally modeled as multidimensional arrays and matrices. Matrices are 2-dimensional rectangular array of elements (or entries) laid out in rows and columns. Although the logical view of a rectangular array of elements need not be the same as the actual physical storage, the elements of the arrays are mapped into storage locations according to some linearization function of the indices. An array file is simply a file of the elements of an array in which the k-dimensional indices of the elements map into linear consecutive record locations in the file. The mapping of k-dimensional indices to linear addresses may either be by a computed access function, an indexing mechanism, or a mixture of both. Operations on these array files involve accessing and manipulating sub-arrays (or array chunks) of one or more of these array files. Such models for manipulating datasets are typical in large scale computing performed by various Scientific and Engineering Simulations, Climate Modeling, High Energy and Nuclear Physics, Astrophysics, Computational Fluid Dynamics, Scientific Visualization, etc.
On-line analytical processing (OLAP) and decision support systems depend highly on efficiently computing statistical summary information from multi-dimensional databases. A view of the dataset as a multi-dimensional array model forms the basis for efficient derivation of summary information. The literature distinguishes OLAP on relational model of data called ROLAP and that on multidimensional model termed MOLAP.
Consider a dataset that maintains monitored values of
temperature at defined locations given by the latitude and
longitude at various time-steps. Figure 1 illustrates a
simple multi-dimensional view of a sample data as a
3dimensional array A[
Each cell stores the temperature value (not shown) as the measure. Suppose we wish to maintain information on sales of some enterprise that has stores at different locations. The locations are defined by their spatial coordinates of latitude and longitude often abbreviated as Lat/Long. The sales information recorded for say three consecutive months can be represented by Figure 1 which gives the same view of the data as before but with values for sales as the measures. A multi-dimensional array correctly models the dataset in both domains of scientific applications and MOLAP. Further common characteristics associated with multi-dimensional model of data under these domains are that: 1. The data incrementally grow over time by appending new data elements. These may reach orders of terabytes, as in data warehousing. 2. The datasets are mainly read-only. However, they may be subject to expansions in the bounds of the dimensions, e.g., as new monitoring locations are enabled or new stores are opened. 3. The number of dimensions (i.e., the rank) of the array may be extended. 4. The array representation can be either dense or sparse.
In the illustration of Figure 1, the array may grow by appending data for new time-steps. When new locations are added, the bounds of the Lat/Long dimensions may be required to be extended. The mapping function depicted in the above figure (represented by the label inside each cell) corresponds to that of conventional array mapping that allows extendibility in one dimension only; namely the dimension in the mapping scheme that is least varying. We desire an array addressing scheme that allows extensions in both the bound and rank of the array file and also efficiently manages sparse arrays. In this paper we address primarily the storage scheme for managing dense extendible array efficiently. We discuss how the technique is tailored for handling sparse arrays as well without detailed experimentation due to space limitation.
Over the years special file formats have been extensively
studied and developed for storing multi-dimensional array
files that support sub-array accesses in high performance
computing. These include NetCDF [
Let A[N0][N1] . . . [Nk−1], denote a k-dimensional array where N j, 0 ≤ j < k − 1, is the bound on the indices of dimension j. An element, denoted by Ahi0, i1 . . . , ik−1i, is referenced by a k-dimensional index hi0, i1 . . . , ik−1i. Let L = {ℓ0, ℓ1 . . . , ℓM−1}, be a sequence of consecutive storage locations where M = ∏kj−=10 N j. An allocation (or mapping) function F (), maps the k-dimensional indices one-to-one, onto the sequence of consecutive indices {0, 1, . . . , M − 1}, i.e., F : N0 × N1 × · · · × Nk−1 → {0, 1, . . . , M − 1}. Given a location index j, the inverse mapping function F −1( j), computes the k-dimensional index hi0, i1 . . . , ik−1i that corresponds to j. The significance of inverse mapping functions have not been much appreciated in the past but has important use in computing the addresses of the neighbors of an element given its linear storage address and also for managing a general k-dimensional sparse array as opposed to sparse matrices which is 2-dimensional.
Modern programming languages provide native support for multidimensional arrays. The mapping function F () is normally defined so that elements are allocated either in row-major or column major order. We refer to arrays whose elements are allocated in this manner as conventional arrays. Conventional arrays limit their growth to only one dimension. We use the term row-major order in a general sense, beyond row-and-column matrices, to mean an order in which the leftmost index of a k-dimensional index is the least varying. This is also sometimes referred to as the lexicographic order. A column-major order refers to an ordering in which the rightmost index varies the slowest.
In conventional k-dimensional arrays, efficient computation of the mapping function is carried out with the aid of a vector that holds k multiplicative coefficients of the respective indices. Consequently, an N-element array in kdimensions, can be maintained in O(N) storage locations and the computation of the linear address of an element, given its k-dimensional coordinate index, is done in time O(k) since this is done by k − 1 multiplications and k − 1 additions. We achieve similar efficiency in the management of extendible arrays by maintaining vectors of multiplying coefficients each time the array expands. The computation of the linear address, corresponding to a k-dimensional index, uses these stored vectors of coefficients. The vectors of multiplicative coefficients maintain records of the history of the expansions and are organized into Axial-Vectors. There is one Axial-Vector for each dimension and holds entries that we refer to as expansion records. The fields of an expansion record are described later.
Consider a k-dimensional array of N elements and after some arbitrary extensions dimension j maintains E j records of the history of expansion for dimension j. Our approaches computes the linear address from the k-dimensional index in time O(k + ∑kj−=01 E j) using O(k ∑kj−=10 E j) additional space.
The technique being presented in this paper can be
adopted by compiler developers for generating mapping
functions for dense multidimensional extendible arrays in
memory as well. The Standard Template Library (STL) in
C++ also provides support for resizable vectors. However,
the general technique for allowing array extensions without
the extensive cost of reallocating storage is what we desire.
The approach we propose may be used in implementing
special libraries such as the Global Array (GA) library [
To handle sparse arrays, we adopt the technique of
array chunking [
The main contributions in this paper are: • A definition of a mapping function for extendible arrays that is applicable for both dense in-core arrays and out-of-core arrays. We call the approach the The Axial-Vector method. The general principle is not new. The idea was first proposed with the use of an auxiliary array. However, the storage overhead using an auxiliary array could be prohibitive. The new Axial-Vector method obviates the storage overhead and also gives a new method for computing the linear addresses. • When extremely large array files are generated and stored, any dimension can still be expanded by appending new array elements without the need to reorganize the already allocated storage which could be of the order of terabytes. • The mapping function proposed in this paper incurs very little storage overhead and can be used as a replacement for the access function for array files such as NetCDF and the data chunks in the HDF5 file format. • We discuss an extension of the technique to handle sparse extendible arrays for both in-core and out-ofcore arrays.
The organization of this paper is as follows. In the next section we present some details of the definition of the mapping function for dense multidimensional extendible arrays in main memory since the same mapping function is adopted for accessing elements of extendible array files. In section 3 we describe how an array file is implemented. We describe how sparse multidimensional array files are managed in section 4. We give some experimental comparison of the array mapping functions for extendible and conventional arrays in section 5. We conclude in section 6 and give some directions for our future work. 2
2.1
First we explore some details of how conventional arrays are mapped onto consecutive storage locations in memory. Consider a k-dimensional array A[N0][N1] . . . [Nk−1], where N j, 0 ≤ j < k − 1, denote the bounds of the index ranges of the respective dimensions. Suppose the elements of this array are allocated in the linear consecutive storage locations L = {ℓ0, ℓ1 . . . , ℓM−1}, in row-major order according to the mapping function F (). In the rest of this paper, we will always assume row-major ordering and for simplicity we will assume an array element occupies a unit of storage. A unit of storage could be one word of 4 bytes, a double word of 8 bytes, etc. An element Ahi0, i1, . . . ik−1i is assigned to location ℓq, where q is computed by the mapping function defined as q = F (hi0, i1, . . . ik−1i) = i0 ∗ C0 + i1 ∗ C1 + · · · + ik−1 ∗ Ck−1 where Cj = k−1 ∏ Nr, 0 ≤ j ≤ k − 1. r= j+1 (1) and Ah0, 0, . . . 0i assigned to ℓ0.
In most programming languages, since the bounds of the arrays are known at compilation time, the coefficients C0,C1, . . . ,Ck−1 are computed and stored during code generation. Consequently, given any k-dimensional index, the computation of the corresponding linear address using Equation 1, takes time O(k).
Suppose we know the linear address q of an array element, the k-dimensional index hi0, i1, . . . ik−1i corresponding to q can be computed by repeated modulus arithmetic with the coefficients Ck−2,Ck−3, . . . ,C1 in turn, i.e., F −1(q) → hi0, i1, . . . ik−1i. When allocating elements of a dense multidimensional array in a file, the same mapping function is used where the linear address q gives the displacement relative to the location of the first element in the file, in units of the size of the array elements. The limitation imposed by F is that the array can only be extended along dimension 0 since the evaluation of the function does not involve the bound of N0.
We can still circumvent the constraint imposed by F
by shuffling the order of the indices whenever the array is
extended along any dimension. The idea of extending the
index range of a dimension is simply to adjoin a block (or a
hyperslab) of array elements whose sizes on all dimensions
remain the same except for the dimension being extended.
The question of organizing an extendible array in
memory, such that the linear addresses can be computed in the
manner similar to those of a static array described above,
is a long standing one [
Consider Figure 2a that shows a map of the storage allocation of 3-dimensional extendible array. We denote this in general as A[N∗0][N∗1][N∗2], where N∗j represents the fact that the bound has the propensity to grow. In this paper we address only the problem of allowing extendibility in the array bounds but not its rank. The labels shown in the array cells represent the linear addresses of the respective elements, as a displacement from the location of the first element.
Suppose initially the array is allocated as A[
The labels shown in the array cells represent the linear
addresses of the respective elements. For example, in
Figure 2a the element Ah2, 1, 0i is assigned to location 7 and
element Ah3, 1, 2i is assigned to location 34. The array was
subsequently extended along the longitude axis by one
index, then along the latitude axis by 2 indices and then along
the time axis by one time-step. A hyperslab of array
elements can be perceived as an array chunk (a term used
in [
Consider now that we have a k-dimensional extendible array A[N∗0][N∗1] . . . [Nk∗−1], for which dimension l is extended by λl, so that the index range increases from Nl∗ to Nl∗ + λl. The strategy is to allocate a hyperslab of array elements such that addresses within the hyperslab are computed as displacements from the location of the first element of the hyperslab. Let the first element of a hyperslab of dimension l be denoted by Ah0, 0, . . . , Nl∗, . . . , 0i. Address calculation is computed in row-major order as before, except that now dimension l is the least varying dimension in the allocation scheme but all other dimensions retain their relative order. Denote the location of Ah0, 0, . . . , Nl∗, . . . , 0i as ℓMl∗ where Ml∗ = ∏rk=−01(Nr∗). Then the desired mapping function F ∗() that computes the address q∗ of a new element Ahi0, i1, . . . ik−1i during the allocation is given by: k−1 q∗ = F ∗(hi0, i1, . . . ik−1i) = Ml∗ + (il − Nl∗)Cl∗ + ∑ i jC∗j j=0 j6=l k−1 where Cl∗ = ∏ N∗
j j=0 j6=l and C∗j = k−1 ∏ Nr∗ r= j+1 r6=l (2)
We need to retain for dimension l the values of Ml∗ - the location of the first element of the hyperslab, Nl∗ - the fisrt index of the adjoined hyperslab, and Cr∗, 0 ≤ r < k - the multiplicative coefficients, in some data structure so that these can be easily retrieved for computating an element’s address within the adjoined hyperslab. The axial-vectors denoted by Γ j[E j], 0 ≤ j < k, and shown in Figure 2b, are used to retain the required information. E j is the number of stored records for axial-vector Γ j. Note that the number of elements in each axial-vector is always less than or equal to the number of indices of the corresponding dimension. It is exactly the number of uninterrupted expansions. In the example of Figure 2b, E 0 = 2, E 1 = 2, and E 2 = 3.
The information of each expansion record of a dimension is a record comprised of four fields. For dimension l, the ith entry denoted by Γl hii consists of Γl hii.Nl∗; Γlhii.Ml∗; Γlhii.C[k] - the stored multiplying coefficients for computing the displacement values within the hyperslab; and Γlhii.Sil - the memory address where the hyperslab is stored. Note however that for computing record addresses of array files, this last field is not required, since new records are always allocated by appending to the existing array file. In main memory an extendible array may be formed as a collection of disjoint hyperslabs since a block of memory acquired for each new hyperslab may not necessarily be contiguous to a previously allocated one. Contiguiety in memory allocation is only guranteed for unintrrupted expansion of a dimension, i.e., when the same dimension is repeatedly expanded.
Given a k-dimensional index hi0, i1, . . . , ik−1i, the main idea in correctly computing the linear address is in determining which of the records Γ0hz0i, Γ1hz1i . . . Γk−1hzk−1i, has the first maximum starting address of its hyperslab. The index z j is given by a modified binary search algorithm that always gives the highest index of the axial-vector where the expansion record has a maximum starting address of its hyperslab less than or equal to i j.
For example, suppose we desire the linear address of the element Ah4, 2, 2i, we first note that z0 = 1, z1 = 0, and Conventional from which we deduce that Ml∗ = 48, l = 0, and Nl∗ = N∗0 = 4. The computation F ∗(h4, 2, 2i) = 48 + 12 × (4 − 4) + 3 × 2 + 1 × 2 = 48 + 0 + 6 + 2 = 56. The value 56 is the linear address relative to the starting address of 0. The main characteristics of the extendible array realization is summarized in the following theorem Theorem 2.1. Suppose that in a k-dimensional extendible array, dimension j undergoes E j, uninterrupted expansions. Then if E = ∑kj=−01 E j, the complexity of computing the function F ∗() for an extendible array using axialvectors is O(k + k log(k + E ) − k logk), using O(kE ) worst case space.
Proof. The worst case sizes of the axial-vectors occur if each dimension has the same number of uninterrupted expansions; i.e., E j = E /k. The evaluation of F ∗() involves k log E j followed by k multiplications k additions and 1 subtraction, giving a total of O(k + k(log(1 + E /k))) = O(k + k log(k + E ) − k logk).
The additional space requirement for the k axial-vectors is O((k + 3) ∑kj=−01 E j) = O(kE ).
Given a linear address q∗, it is easy to find, an entry in the axial vectors that has the maximum starting address less than or equal to q∗ using a modified binary search algorithm as in the address computation. The repeated modulus arithmetic as described in section 2.1 is then used to extract the k-dimensional indices. We state without a formal proof the following.
Theorem 2.2. Given a linear address q of an element of an extendible array realized with the aid of k axial-vectors, the the k-dimensional index of an element is computable by the function F ∗−1 in time O(k + log E ).
The main results of the extendible array realization are summarized in Table 1.
An array file is simply a persistent storage of the corresponding main memory resident array in a file, augmented with some meta-data information, either as a header in the same file or in a separated file. We will consider array files as formed in pairs: the primary file Fp and the meta-data file Fm. For extremely large files, the content in memory at any time is a subset, or a subarray of the entire disk resident array file. Scientific applications that process these arrays consider array chunks as the unit of data access from secondary storage. Most high performance computations are executed as a parallel program either on a cluster of workstations or on massively parallel machines. The model of the data is a large global array from which subarrays are accessed into individual workstations. The subarrays of individual nodes together constitute tiles of the global array. Actual physical access of array elements is carried out in units of array chunks for both dense and sparse arrays. We discuss this in some detail in the next section.
For the subsequent discussions, we will ignore most of the details of the organization of the array files and also the details of the structure of the array elements. For simplicity, we consider the array files as being composed of fixed size elements, where each element is composed of a fixed number of attributes. Each attribute value has a corresponding ordinal number that serves as the index value. Mapping of attribute values to ordinal numbers is easily done for array files. The primary file Fp, contains elements of the multidimensional array that continuously grows by appending new data elements whenever a dimension is extended. For example, a dimension corresponding to time may be extended by the addition of data from new time-steps. A meta-data file Fm stores the records that correspond to the axial vectors. The contents of the meta-data file Fm are used to reconstruct the memory resident axial vectors. Each expansion of a dimension results in an update of an axial vector and consequently the meta-data file as well.
Suppose an application has already constructed the memory resident axial vectors from the meta-data file, then the linear address of an element (i.e., a record), of the array file given its k-dimensional coordinates, is computed using the mapping function F ∗(). Essentially, F ∗() serves as a hash function for the elements of the array file. Conversely, if the linear address q∗, of an element is given and one desires the neighbor element that lies some units of coordinate distances along some specified dimensions from the current, such an element can be easily retrieved with the aid of the inverse function F ∗−1(). First we compute the k-dimensional coordinate values from F ∗−1(q), adjust the coordinate values along the specified dimensions and compute the address of the desired element by F ∗(). The relevant algorithms for these operations are easily derived from the definitions of the functions presented in the preceding sections. One of the main features of our scheme is that when extremely large array files are generated, each dimensions can still be expanded by appending new array elements without the need to reorganize the allocated storage of the files.
Two popular data organization schemes for large scale
scientific datasets are NetCDF [
Besides the characteristics that multi-dimensional
databases incrementally grow over time, they also
have the unkind property of being sparse. Techniques
for managing large scale storage of multi-dimensional
data have either addressed the sparsity problem of the
array model [
Figure 3a shows a 3-dimensional array partitioned into chunks using index-intervals of 3. Addressing elements of the array is computed in two levels. The first level address computation gives the chunk address of the element. The second level address is computed as the displacement of the array element within the array chunk. The extendible array addressing method maps the k-dimensional index of each chunk into a Table Map. The table map of the chunks contains pointers to the physical addresses of the data blocks that hold chunks. An array chunk forms the unit of data transfer between secondary storage and main memory. A table map pointer is set to null if the chunk holds no array elements. As in extendible hashing, the use of table map to address chunks guarantees at most 2 disk accesses to locate a chunk. Array chunks that have less than some defined number of array elements can be compressed further.
The table map is one of the simplest techniques for
handling sparse extendible array but suffers from the
problem of a potential exponential growth for high
dimensional datasets. Instead of a table map for maintaining
array chunks, we use a dynamically constructed PATRICIA
trie [
The use of a table map for locating the physical locations of array chunks relaxes the need for directly addressing array elements with an extendible array mapping function. One can further relax this constraint by doing away with the extendible array mapping function entirely. Rather, a method for constructing a unique address Ihi0,...ik0 i of an array chunk from the k-dimensional index hi0, . . . ik0 i is all that is required. One such method is given by concatenating the binary representation of the coordinate indices of an array chunk. The unique address generated is then used in an index scheme such as a B+-tree, to locate the physical chunk where an array element is stored. This approach is actually implemented in the HDF5 storage format. There are two problems with this approach; 1. Either the k-dimensional index or the generated identifier for the chunk must be stored with the chunk. For the latter case, an inverse function for computing the k-dimensional chunk index from the chunk identifier is needed but is less space consuming. 2. It does not handle both memory resident array and disk resident arrays uniformly.
In general the table map can be replaced with any indexing scheme that maintains O(N) chunk address for exactly N non-empty chunks. The use of a PATRICIA trie guarantees that. Using extendible array mapping function for computing the linear address of a chunk has the advantage of: 1. giving us a uniform manner of managing extendible arrays resident both in-core and out-of-core. 2. allowing the replacement of the global k-dimensional address of an array element by one which only defines its linear location within an array chunk and yet enables us to compute the global k-dimensional index from the linear address. 4.2
Besides the creation of the multi-dimensional array files,
and operations for reading, writing (i.e., accessing ) and
appending new element, the application domain dictates
the type of operations the array files are subjected to.
While multi-dimensional OLAP applications see the
efficient computation of the Cube operator [
The mapping function provides direct addressing for array chunks and subsequently to the array elements. A naive approach to performing sub-array extraction operation would be to iterate over the k-dimensional coordinates of the elements to be selected and retrieve each array elements independently. Suppose the cardinality of the response set of the first selection class is ℜ. The naive approach performs ℜ independent disk accesses. But one can do better than this worst case number of disk accesses. An efficient method for processing any of the above queries is to compute the chunk identifies of the element; and for each chunk retrieved, extract all elements the chunk that satisfies the request. Due to space limitation, detailed discussions on the structures, algorithms and experiments on extendible sparse multidimensional arrays is left out in this paper. 5
The theoretical analysis of the mapping function for extendible arrays indicates that it is nearly of the same order of computational complexity as that of conventional arrays. The main difference being the additonal time requied by the mapping function for an extendible array to perform binary searches in the axial-vectors. We experimentally tested this by computing the average access times of both the conventional array and extendible array for an array size of approximately 108 elements of double data types. We varied the rank of the array from 2 to 8 while keeping the size of the array about the same. We plotted the average time over 10000 random element access for static arrays. These experiments were run on a 1.5GHz AMD Athlon processor running Centos-4 Linux with 2GByte memory. Figure 4a show the graphs of the access times averaged over 10000 element access.
The graphs indicate that for static arrays, the cost of computing the extendible array access function varies more significantly with the rank of the array than for the conventional array access function. Even though the complexity of computing the access functions are both O(k), the extendible array access function shows strong dependence on the array’s rank k due to the fact that k binary searches are done in the axial-vectors.
We also considered the impact on the average access cost when the arrays undergo interleaved expansions. The experiment considered arrays of ranks 2 and 3 where the initial array size grew from about 10000 elements to about 106 elements. For each sequence of about 10000 accesses the array is allowed to undergo up to 16 expansions. A dimension selected at random, is extended by a random integer amount of between 1 and 10. Each time the conventional array is extended, the storage allocation is reorganized. The average cost of accessing array elements with interleaved expansions is shown in Figure 4b. Under this model, we find that the average cost of accessing array elements for extendible arrays is significantly less than for conventional arrays that incur the additional cost of reorganization.
Similar experiments were conducted, for accessing elements of array files instead of memory resident arrays. Figure 5a compares the average time to access elements for 2, 3 and 4 dimensional static files. There is very little variation in the times of the conventional and extendible array functions. However, when these times are computed with interleaved expansions, the extendible array clearly outperforms the conventional array methods by several orders of magnitude. Figure 5b shows the graphs for 2, 3 and 4dimensional extendible array files only. The extra time and storage required to reorganize the conventional array files with interleaved expansions become prohibitive as the array becomes large. One can infer from these results that, in handling large scale dense multi-dimensional dataset that incrementally grow by appending elements, the extendible array mapping function should be the choice for addressing storage.
We should mention that a competitive model for comparisons of multi-dimensional array file implementations should be with the HDF5 data schemes. Our implementation currently does not include array caching of data pages which the HDF5 implementation uses. Work is still ongoing to add main memory buffer pools for data caching at which time a fair comparison would be made. 6
We have shown how a k-dimensional extendible array file can be used to implement multi-dimensional databases. The technique applies to extendible arrays in-core just as much as for out-of-core extendible arrays that can be either dense or sparse. The method relies on a mapping function that uses information retained in axial-vectors to compute the linear storage addresses from the k-dimensional indices. Given the characteristics of multi-dimensional databases that they incrementally grow into terabytes of data, developing a mapping function that does not require reorganization of the array file as the file grows is a desirable future.
The method proposed is highly appropriate for most scientific datasets where the model of the data is perceived typical as large global array files. The mapping function developed can be used to enhance current implementations of array files such as NetCDF, HDF5 and Global Arrays. Work is still on-going to incorporate our proposed solution to multidimensional array libraries and to extend the technique for multi-dimensional datasets whose dimensions or ranks are allowed to expand. We are also conducting comparative studies on the different techniques for managing sparse multi-dimensional extendible arrays.
This work is supported by the Director, Office of Laboratory Policy and Infrastructure Management of the U. S. Department of Energy under Contract No. DE-AC0376SF00098. This research used resources of the National Energy Research Scientific Computing (NERSC), which is supported by the Office of Science of the U.S. Department of Energy. 1 3 2 6 24 27 30 33 25 28 31 34 26 29 32 35 a 3-D extendible array (a) A 3-dimensional array partitioned into chunks 2-D Array File 3-D Array File 4-D Array File 500000 1e+06 (b) Element access costs with interleaved extensions (b) Access cost from a file with interleaved extensions type double