interlace.cc

数值计算工具库,C语言编写的,可以直接调用.

/***************************************************************************
 * blitz/array/interlace.cc
 *
 * $Id$
 *
 * Copyright (C) 1997,1998 Todd Veldhuizen <[e-mail]>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * Suggestions:          [e-mail]
 * Bugs:                 [e-mail]
 *
 * For more information, please see the Blitz++ Home Page:
 *    http://seurat.uwaterloo.ca/blitz/
 *
 ***************************************************************************
 * $Log$
 *
 */

#ifndef BZ_ARRAYINTERLACE_CC
#define BZ_ARRAYINTERLACE_CC

#ifndef BZ_ARRAY_H
 #error <blitz/array/interlace.cc> must be included via <blitz/array.h>
#endif

#ifndef BZ_ARRAYSHAPE_H
 #include <blitz/array/shape.h>
#endif

BZ_NAMESPACE(blitz)

/*
 * This header provides two collections of routines:
 *
 * interlaceArrays(shape, A1, A2, ...);
 * allocateArrays(shape, A1, A2, ...);
 *
 * interlaceArrays allocates a set of arrays so that their data is
 * interlaced.  For example,
 *
 * Array<int,2> A, B;
 * interlaceArrays(shape(10,10), A, B);
 *
 * sets up the array storage so that A(0,0) is followed by B(0,0) in
 * memory; then A(0,1) and B(0,1), and so on.
 *
 * The allocateArrays() routines may or may not interlace the data,
 * depending on whether interlacing is advantageous for the architecture.
 * This is controlled by the setting of BZ_INTERLACE_ARRAYS in
 * <blitz/tuning.h>.
 */

// Warning: Can't instantiate TinyVector<Range,N> because causes
// conflict between TinyVector<T,N>::operator=(T) and
// TinyVector<T,N>::operator=(Range)

// NEEDS_WORK -- also shape for up to N=11
// NEEDS_WORK -- shape(Range r1, Range r2, ...) (return TinyVector<Range,n>)
//               maybe use Domain objects
// NEEDS_WORK -- doesn't make a lot of sense for user to provide a
//               GeneralArrayStorage<N_rank+1>

template<class T_numtype>
void makeInterlacedArray(Array<T_numtype,2>& mainArray,
    Array<T_numtype,1>& subarray, int slice)
{
    Array<T_numtype,1> tmp = mainArray(Range::all(), slice);
    subarray.reference(tmp);
}

template<class T_numtype>
void makeInterlacedArray(Array<T_numtype,3>& mainArray,
    Array<T_numtype,2>& subarray, int slice)
{
    Array<T_numtype,2> tmp = mainArray(Range::all(), Range::all(), 
        slice);
    subarray.reference(tmp);
}

template<class T_numtype>
void makeInterlacedArray(Array<T_numtype,4>& mainArray,
    Array<T_numtype,3>& subarray, int slice)
{
    Array<T_numtype,3> tmp = mainArray(Range::all(), Range::all(), 
        Range::all(), slice);
    subarray.reference(tmp);
}

// These routines always allocate interlaced arrays
template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 2, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 3, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 4, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 5, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 6, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
    Array<T_numtype,N_rank>& a7)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 7, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
    makeInterlacedArray(array, a7, 6);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 8, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
    makeInterlacedArray(array, a7, 6);
    makeInterlacedArray(array, a8, 7);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
    Array<T_numtype,N_rank>& a9)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 9, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
    makeInterlacedArray(array, a7, 6);
    makeInterlacedArray(array, a8, 7);
    makeInterlacedArray(array, a9, 8);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 10, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
    makeInterlacedArray(array, a7, 6);
    makeInterlacedArray(array, a8, 7);
    makeInterlacedArray(array, a9, 8);
    makeInterlacedArray(array, a10, 9);
}

template<class T_numtype, int N_rank>
void interlaceArrays(const TinyVector<int,N_rank>& shape,
    Array<T_numtype,N_rank>& a1, Array<T_numtype,N_rank>& a2,
    Array<T_numtype,N_rank>& a3, Array<T_numtype,N_rank>& a4,
    Array<T_numtype,N_rank>& a5, Array<T_numtype,N_rank>& a6,
    Array<T_numtype,N_rank>& a7, Array<T_numtype,N_rank>& a8,
    Array<T_numtype,N_rank>& a9, Array<T_numtype,N_rank>& a10,
    Array<T_numtype,N_rank>& a11)
{
    GeneralArrayStorage<N_rank+1> storage;
    Array<T_numtype, N_rank+1> array(shape, 11, storage);
    makeInterlacedArray(array, a1, 0);
    makeInterlacedArray(array, a2, 1);
    makeInterlacedArray(array, a3, 2);
    makeInterlacedArray(array, a4, 3);
    makeInterlacedArray(array, a5, 4);
    makeInterlacedArray(array, a6, 5);
    makeInterlacedArray(array, a7, 6);
    makeInterlacedArray(array, a8, 7);
    makeInterlacedArray(array, a9, 8);
    makeInterlacedArray(array, a10, 9);
    makeInterlacedArray(array, a11, 10);
}

// NEEDS_WORK -- make `storage' a parameter in these routines
//  Will be tricky: have to convert GeneralArrayStorage<N_rank> to
//  GeneralArrayStorage<N_rank+1>

// These routines may or may not interlace arrays, depending on
// whether it is advantageous for this platform.