source: XIOS/dev/dev_olga/src/extern/blitz/include/blitz/array/interlace.cc @ 1022

Last change on this file since 1022 was 1022, checked in by mhnguyen, 7 years ago
File size: 17.9 KB
Line 
1/***************************************************************************
2 * blitz/array/interlace.cc   Interlaced array storage.
3 *
4 * $Id$
5 *
6 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
7 *
8 * This file is a part of Blitz.
9 *
10 * Blitz is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation, either version 3
13 * of the License, or (at your option) any later version.
14 *
15 * Blitz is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 * GNU Lesser General Public License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with Blitz.  If not, see <http://www.gnu.org/licenses/>.
22 *
23 * Suggestions:          blitz-devel@lists.sourceforge.net
24 * Bugs:                 blitz-support@lists.sourceforge.net   
25 *
26 * For more information, please see the Blitz++ Home Page:
27 *    https://sourceforge.net/projects/blitz/
28 *
29 ****************************************************************************/
30#ifndef BZ_ARRAYINTERLACE_CC
31#define BZ_ARRAYINTERLACE_CC
32
33#ifndef BZ_ARRAY_H
34 #error <blitz/array/interlace.cc> must be included via <blitz/array.h>
35#endif
36
37#ifndef BZ_ARRAYSHAPE_H
38 #include <blitz/array/shape.h>
39#endif
40
41BZ_NAMESPACE(blitz)
42
43/*
44 * This header provides two collections of routines:
45 *
46 * interlaceArrays(shape, A1, A2, ...);
47 * allocateArrays(shape, A1, A2, ...);
48 *
49 * interlaceArrays allocates a set of arrays so that their data is
50 * interlaced.  For example,
51 *
52 * Array<int,2> A, B;
53 * interlaceArrays(shape(10,10), A, B);
54 *
55 * sets up the array storage so that A(0,0) is followed by B(0,0) in
56 * memory; then A(0,1) and B(0,1), and so on.
57 *
58 * The allocateArrays() routines may or may not interlace the data,
59 * depending on whether interlacing is advantageous for the architecture.
60 * This is controlled by the setting of BZ_INTERLACE_ARRAYS in
61 * <blitz/tuning.h>.
62 */
63
64// Warning: Can't instantiate TinyVector<Range,N> because causes
65// conflict between TinyVector<T,N>::operator=(T) and
66// TinyVector<T,N>::operator=(Range)
67
68// NEEDS_WORK -- also shape for up to N=11
69// NEEDS_WORK -- shape(Range r1, Range r2, ...) (return TinyVector<Range,n>)
70//               maybe use Domain objects
71// NEEDS_WORK -- doesn't make a lot of sense for user to provide a
72//               GeneralArrayStorage<N_rank+1>
73
74template<typename T_numtype>
75void makeInterlacedArray(Array<T_numtype,2>& mainArray,
76    Array<T_numtype,1>& subarray, int slice)
77{
78    Array<T_numtype,1> tmp = mainArray(Range::all(), slice);
79    subarray.reference(tmp);
80}
81
82template<typename T_numtype>
83void makeInterlacedArray(Array<T_numtype,3>& mainArray,
84    Array<T_numtype,2>& subarray, int slice)
85{
86    Array<T_numtype,2> tmp = mainArray(Range::all(), Range::all(), 
87        slice);
88    subarray.reference(tmp);
89}
90
91template<typename T_numtype>
92void makeInterlacedArray(Array<T_numtype,4>& mainArray,
93    Array<T_numtype,3>& subarray, int slice)
94{
95    Array<T_numtype,3> tmp = mainArray(Range::all(), Range::all(), 
96        Range::all(), slice);
97    subarray.reference(tmp);
98}
99
100// These routines always allocate interlaced arrays
101template<typename T_numtype, int N_rank>
102void interlaceArrays(const TinyVector<int,N_rank>& shape,
103    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2)
104{
105    GeneralArrayStorage<N_rank+1> storage(contiguousData);
106    Array<T_numtype, N_rank+1> array(shape, 2, storage);
107    makeInterlacedArray(array, a1, 0);
108    makeInterlacedArray(array, a2, 1);
109}
110
111template<typename T_numtype, int N_rank>
112void interlaceArrays(const TinyVector<int,N_rank>& shape,
113    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
114    Array<T_numtype,N_rank>& a3)
115{
116    GeneralArrayStorage<N_rank+1> storage(contiguousData);
117    Array<T_numtype, N_rank+1> array(shape, 3, storage);
118    makeInterlacedArray(array, a1, 0);
119    makeInterlacedArray(array, a2, 1);
120    makeInterlacedArray(array, a3, 2);
121}
122
123template<typename T_numtype, int N_rank>
124void interlaceArrays(const TinyVector<int,N_rank>& shape,
125    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
126    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4)
127{
128    GeneralArrayStorage<N_rank+1> storage(contiguousData);
129    Array<T_numtype, N_rank+1> array(shape, 4, storage);
130    makeInterlacedArray(array, a1, 0);
131    makeInterlacedArray(array, a2, 1);
132    makeInterlacedArray(array, a3, 2);
133    makeInterlacedArray(array, a4, 3);
134}
135
136template<typename T_numtype, int N_rank>
137void interlaceArrays(const TinyVector<int,N_rank>& shape,
138    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
139    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
140    Array<T_numtype,N_rank>& a5)
141{
142    GeneralArrayStorage<N_rank+1> storage(contiguousData);
143    Array<T_numtype, N_rank+1> array(shape, 5, storage);
144    makeInterlacedArray(array, a1, 0);
145    makeInterlacedArray(array, a2, 1);
146    makeInterlacedArray(array, a3, 2);
147    makeInterlacedArray(array, a4, 3);
148    makeInterlacedArray(array, a5, 4);
149}
150
151template<typename T_numtype, int N_rank>
152void interlaceArrays(const TinyVector<int,N_rank>& shape,
153    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
154    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
155    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6)
156{
157    GeneralArrayStorage<N_rank+1> storage(contiguousData);
158    Array<T_numtype, N_rank+1> array(shape, 6, storage);
159    makeInterlacedArray(array, a1, 0);
160    makeInterlacedArray(array, a2, 1);
161    makeInterlacedArray(array, a3, 2);
162    makeInterlacedArray(array, a4, 3);
163    makeInterlacedArray(array, a5, 4);
164    makeInterlacedArray(array, a6, 5);
165}
166
167template<typename T_numtype, int N_rank>
168void interlaceArrays(const TinyVector<int,N_rank>& shape,
169    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
170    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
171    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
172    Array<T_numtype,N_rank>& a7)
173{
174    GeneralArrayStorage<N_rank+1> storage(contiguousData);
175    Array<T_numtype, N_rank+1> array(shape, 7, storage);
176    makeInterlacedArray(array, a1, 0);
177    makeInterlacedArray(array, a2, 1);
178    makeInterlacedArray(array, a3, 2);
179    makeInterlacedArray(array, a4, 3);
180    makeInterlacedArray(array, a5, 4);
181    makeInterlacedArray(array, a6, 5);
182    makeInterlacedArray(array, a7, 6);
183}
184
185template<typename T_numtype, int N_rank>
186void interlaceArrays(const TinyVector<int,N_rank>& shape,
187    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
188    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
189    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
190    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8)
191{
192    GeneralArrayStorage<N_rank+1> storage(contiguousData);
193    Array<T_numtype, N_rank+1> array(shape, 8, storage);
194    makeInterlacedArray(array, a1, 0);
195    makeInterlacedArray(array, a2, 1);
196    makeInterlacedArray(array, a3, 2);
197    makeInterlacedArray(array, a4, 3);
198    makeInterlacedArray(array, a5, 4);
199    makeInterlacedArray(array, a6, 5);
200    makeInterlacedArray(array, a7, 6);
201    makeInterlacedArray(array, a8, 7);
202}
203
204template<typename T_numtype, int N_rank>
205void interlaceArrays(const TinyVector<int,N_rank>& shape,
206    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
207    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
208    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
209    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
210    Array<T_numtype,N_rank>& a9)
211{
212    GeneralArrayStorage<N_rank+1> storage(contiguousData);
213    Array<T_numtype, N_rank+1> array(shape, 9, storage);
214    makeInterlacedArray(array, a1, 0);
215    makeInterlacedArray(array, a2, 1);
216    makeInterlacedArray(array, a3, 2);
217    makeInterlacedArray(array, a4, 3);
218    makeInterlacedArray(array, a5, 4);
219    makeInterlacedArray(array, a6, 5);
220    makeInterlacedArray(array, a7, 6);
221    makeInterlacedArray(array, a8, 7);
222    makeInterlacedArray(array, a9, 8);
223}
224
225template<typename T_numtype, int N_rank>
226void interlaceArrays(const TinyVector<int,N_rank>& shape,
227    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
228    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
229    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
230    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
231    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10)
232{
233    GeneralArrayStorage<N_rank+1> storage(contiguousData);
234    Array<T_numtype, N_rank+1> array(shape, 10, storage);
235    makeInterlacedArray(array, a1, 0);
236    makeInterlacedArray(array, a2, 1);
237    makeInterlacedArray(array, a3, 2);
238    makeInterlacedArray(array, a4, 3);
239    makeInterlacedArray(array, a5, 4);
240    makeInterlacedArray(array, a6, 5);
241    makeInterlacedArray(array, a7, 6);
242    makeInterlacedArray(array, a8, 7);
243    makeInterlacedArray(array, a9, 8);
244    makeInterlacedArray(array, a10, 9);
245}
246
247template<typename T_numtype, int N_rank>
248void interlaceArrays(const TinyVector<int,N_rank>& shape,
249    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
250    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
251    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
252    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
253    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10,
254    Array<T_numtype,N_rank>& a11)
255{
256    GeneralArrayStorage<N_rank+1> storage(contiguousData);
257    Array<T_numtype, N_rank+1> array(shape, 11, storage);
258    makeInterlacedArray(array, a1, 0);
259    makeInterlacedArray(array, a2, 1);
260    makeInterlacedArray(array, a3, 2);
261    makeInterlacedArray(array, a4, 3);
262    makeInterlacedArray(array, a5, 4);
263    makeInterlacedArray(array, a6, 5);
264    makeInterlacedArray(array, a7, 6);
265    makeInterlacedArray(array, a8, 7);
266    makeInterlacedArray(array, a9, 8);
267    makeInterlacedArray(array, a10, 9);
268    makeInterlacedArray(array, a11, 10);
269}
270
271// NEEDS_WORK -- make `storage' a parameter in these routines
272//  Will be tricky: have to convert GeneralArrayStorage<N_rank> to
273//  GeneralArrayStorage<N_rank+1>
274
275// These routines may or may not interlace arrays, depending on
276// whether it is advantageous for this platform.
277
278template<typename T_numtype, int N_rank>
279void allocateArrays(const TinyVector<int,N_rank>& shape,
280    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2)
281{
282#ifdef BZ_INTERLACE_ARRAYS
283    interlaceArrays(shape, a1, a2);
284#else
285    a1.resize(shape);
286    a2.resize(shape);
287#endif
288}
289
290template<typename T_numtype, int N_rank>
291void allocateArrays(const TinyVector<int,N_rank>& shape,
292    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
293    Array<T_numtype,N_rank>& a3)
294{
295#ifdef BZ_INTERLACE_ARRAYS
296    interlaceArrays(shape, a1, a2, a3);
297#else
298    a1.resize(shape);
299    a2.resize(shape);
300    a3.resize(shape);
301#endif
302}
303
304template<typename T_numtype, int N_rank>
305void allocateArrays(const TinyVector<int,N_rank>& shape,
306    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
307    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4)
308{
309#ifdef BZ_INTERLACE_ARRAYS
310    interlaceArrays(shape, a1, a2, a3, a4);
311#else
312    a1.resize(shape);
313    a2.resize(shape);
314    a3.resize(shape);
315    a4.resize(shape);
316#endif
317}
318
319template<typename T_numtype, int N_rank>
320void allocateArrays(const TinyVector<int,N_rank>& shape,
321    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
322    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
323    Array<T_numtype,N_rank>& a5)
324{
325#ifdef BZ_INTERLACE_ARRAYS
326    interlaceArrays(shape, a1, a2, a3, a4, a5);
327#else
328    a1.resize(shape);
329    a2.resize(shape);
330    a3.resize(shape);
331    a4.resize(shape);
332    a5.resize(shape);
333#endif
334}
335
336template<typename T_numtype, int N_rank>
337void allocateArrays(const TinyVector<int,N_rank>& shape,
338    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
339    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
340    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6)
341{
342#ifdef BZ_INTERLACE_ARRAYS
343    interlaceArrays(shape, a1, a2, a3, a4, a5, a6);
344#else
345    a1.resize(shape);
346    a2.resize(shape);
347    a3.resize(shape);
348    a4.resize(shape);
349    a5.resize(shape);
350    a6.resize(shape);
351#endif
352}
353
354template<typename T_numtype, int N_rank>
355void allocateArrays(const TinyVector<int,N_rank>& shape,
356    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
357    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
358    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
359    Array<T_numtype,N_rank>& a7)
360{
361#ifdef BZ_INTERLACE_ARRAYS
362    interlaceArrays(shape, a1, a2, a3, a4, a5, a6, a7);
363#else
364    a1.resize(shape);
365    a2.resize(shape);
366    a3.resize(shape);
367    a4.resize(shape);
368    a5.resize(shape);
369    a6.resize(shape);
370    a7.resize(shape);
371#endif
372}
373
374template<typename T_numtype, int N_rank>
375void allocateArrays(const TinyVector<int,N_rank>& shape,
376    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
377    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
378    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
379    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8)
380{
381#ifdef BZ_INTERLACE_ARRAYS
382    interlaceArrays(shape, a1, a2, a3, a4, a5, a6, a7, a8);
383#else
384    a1.resize(shape);
385    a2.resize(shape);
386    a3.resize(shape);
387    a4.resize(shape);
388    a5.resize(shape);
389    a6.resize(shape);
390    a7.resize(shape);
391    a8.resize(shape);
392#endif
393}
394
395template<typename T_numtype, int N_rank>
396void allocateArrays(const TinyVector<int,N_rank>& shape,
397    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
398    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
399    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
400    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
401    Array<T_numtype,N_rank>& a9)
402{
403#ifdef BZ_INTERLACE_ARRAYS
404    interlaceArrays(shape, a1, a2, a3, a4, a5, a6, a7, a8, a9);
405#else
406    a1.resize(shape);
407    a2.resize(shape);
408    a3.resize(shape);
409    a4.resize(shape);
410    a5.resize(shape);
411    a6.resize(shape);
412    a7.resize(shape);
413    a8.resize(shape);
414    a9.resize(shape);
415#endif
416}
417
418template<typename T_numtype, int N_rank>
419void allocateArrays(const TinyVector<int,N_rank>& shape,
420    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
421    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
422    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
423    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
424    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10)
425{
426#ifdef BZ_INTERLACE_ARRAYS
427    interlaceArrays(shape, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10);
428#else
429    a1.resize(shape);
430    a2.resize(shape);
431    a3.resize(shape);
432    a4.resize(shape);
433    a5.resize(shape);
434    a6.resize(shape);
435    a7.resize(shape);
436    a8.resize(shape);
437    a9.resize(shape);
438    a10.resize(shape);
439#endif
440}
441
442template<typename T_numtype, int N_rank>
443void allocateArrays(const TinyVector<int,N_rank>& shape,
444    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
445    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
446    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
447    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
448    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10,
449    Array<T_numtype,N_rank>& a11)
450{
451#ifdef BZ_INTERLACE_ARRAYS
452    interlaceArrays(shape, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11);
453#else
454    a1.resize(shape);
455    a2.resize(shape);
456    a3.resize(shape);
457    a4.resize(shape);
458    a5.resize(shape);
459    a6.resize(shape);
460    a7.resize(shape);
461    a8.resize(shape);
462    a9.resize(shape);
463    a10.resize(shape);
464    a11.resize(shape);
465#endif
466}
467
468// NEEDS_WORK -- allocateArrays for TinyVector<Range,N_rank>
469
470/** This constructor is used to create interlaced arrays. */
471template<typename T_numtype, int N_rank>
472Array<T_numtype,N_rank>::Array(const TinyVector<int,N_rank-1>& shape,
473    int lastExtent, const GeneralArrayStorage<N_rank>& storage)
474    : storage_(storage)
475{
476    // Create an array with the given shape, plus an extra dimension
477    // for the number of arrays being allocated.  This extra dimension
478    // must have minor storage order.
479
480    if (ordering(0) == 0)
481    {
482        // Column major storage order (or something like it)
483        length_[0] = lastExtent;
484        storage_.setBase(0,0);
485        for (int i=1; i < N_rank; ++i)
486            length_[i] = shape[i-1];
487    }
488    else if (ordering(0) == N_rank-1)
489    {
490        // Row major storage order (or something like it)
491        for (int i=0; i < N_rank-1; ++i)
492            length_[i] = shape[i];
493        length_[N_rank-1] = lastExtent;
494        storage_.setBase(N_rank-1, 0);
495    }
496    else {
497        BZPRECHECK(0, "Used allocateArrays() with a peculiar storage format");
498    }
499
500    setupStorage(N_rank-1);
501}
502
503// NEEDS_WORK -- see note about TinyVector<Range,N> in <blitz/arrayshape.h>
504#if 0
505template<typename T_numtype, int N_rank>
506Array<T_numtype,N_rank>::Array(const TinyVector<Range,N_rank-1>& shape,
507    int lastExtent, const GeneralArrayStorage<N_rank>& storage)
508    : storage_(storage)
509{
510#ifdef BZ_DEBUG
511    for (int i=0; i < N_rank; ++i)
512      BZPRECHECK(shape[i].isAscendingContiguous(),
513        "In call to allocateArrays(), a Range object is not ascending" << endl
514        << "contiguous: " << shape[i] << endl);
515#endif
516
517    if (ordering(0) == 0)
518    {
519        // Column major storage order (or something like it)
520        length_[0] = lastExtent;
521        storage_.setBase(0,0);
522        for (int i=1; i < N_rank; ++i)
523        {
524            length_[i] = shape[i-1].length();
525            storage_.setBase(i, shape[i-1].first());
526        }
527    }
528    else if (ordering(0) == N_rank-1)
529    {
530        // Row major storage order (or something like it)
531        for (int i=0; i < N_rank-1; ++i)
532        {
533            length_[i] = shape[i];
534            storage_.setBase(i, shape[i].first());
535        }
536        length_[N_rank-1] = lastExtent;
537        storage_.setBase(N_rank-1, 0);
538    }
539    else {
540        BZPRECHECK(0, "Used allocateArrays() with a peculiar storage format");
541    }
542
543    setupStorage(N_rank-1);
544}
545#endif
546
547BZ_NAMESPACE_END
548
549#endif // BZ_ARRAYINTER_CC
550
Note: See TracBrowser for help on using the repository browser.