source: XMLIO_V2/external/include/blitz/array/interlace.cc @ 80

Last change on this file since 80 was 80, checked in by ymipsl, 14 years ago

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