stencils.cc

著名的数学计算类库

template<typename T_stencil, typename T_numtype1, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8, typename T_array9, typename T_array10,
    class T_array11>
void applyStencil_imp(const T_stencil& stencil, Array<T_numtype1,2>& A,
    T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F, 
    T_array7& G, T_array8& H, T_array9& I, T_array10& J, T_array11& K)
{
    checkShapes(A,B,C,D,E,F,G,H,I,J,K);

    // Determine stencil extent
    TinyVector<int,2> minb, maxb;
    getStencilExtent(minb, maxb, stencil, A, B, C, D, E, F, G, H, I, J, K);

    // Now determine the subdomain over which the stencil
    // can be applied without worrying about overrunning the
    // boundaries of the array
    int stencil_lbound0 = minb(0);
    int stencil_lbound1 = minb(1);

    int stencil_ubound0 = maxb(0);
    int stencil_ubound1 = maxb(1);

    int lbound0 = minmax::max(A.lbound(0), A.lbound(0) - stencil_lbound0);
    int lbound1 = minmax::max(A.lbound(1), A.lbound(1) - stencil_lbound1);

    int ubound0 = minmax::min(A.ubound(0), A.ubound(0) - stencil_ubound0);
    int ubound1 = minmax::min(A.ubound(1), A.ubound(1) - stencil_ubound1);

#if 0
    cout << "Stencil bounds are:" << endl
     << lbound0 << '\t' << ubound0 << endl
     << lbound1 << '\t' << ubound1 << endl;
#endif 

    // Now do the actual loop
    FastArrayIterator<T_numtype1,2> Aiter(A);
    _bz_typename T_array2::T_iterator Biter(B);
    _bz_typename T_array3::T_iterator Citer(C);
    _bz_typename T_array4::T_iterator Diter(D);
    _bz_typename T_array5::T_iterator Eiter(E);
    _bz_typename T_array6::T_iterator Fiter(F);
    _bz_typename T_array7::T_iterator Giter(G);
    _bz_typename T_array8::T_iterator Hiter(H);
    _bz_typename T_array9::T_iterator Iiter(I);
    _bz_typename T_array10::T_iterator Jiter(J);
    _bz_typename T_array11::T_iterator Kiter(K);

    // Load the strides for the innermost loop
    Aiter.loadStride(1);
    Biter.loadStride(1);
    Citer.loadStride(1);
    Diter.loadStride(1);
    Eiter.loadStride(1);
    Fiter.loadStride(1);
    Giter.loadStride(1);
    Hiter.loadStride(1);
    Iiter.loadStride(1);
    Jiter.loadStride(1);
    Kiter.loadStride(1);

    for (int i=lbound0; i <= ubound0; ++i)
    {
        Aiter.moveTo(i,lbound1);
        Biter.moveTo(i,lbound1);
        Citer.moveTo(i,lbound1);
        Diter.moveTo(i,lbound1);
        Eiter.moveTo(i,lbound1);
        Fiter.moveTo(i,lbound1);
        Giter.moveTo(i,lbound1);
        Hiter.moveTo(i,lbound1);
        Iiter.moveTo(i,lbound1);
        Jiter.moveTo(i,lbound1);
        Kiter.moveTo(i,lbound1);

        for (int k=lbound1; k <= ubound1; ++k)
        {
            stencil.apply(Aiter, Biter, Citer, Diter, Eiter, Fiter, Giter,
                Hiter, Iiter, Jiter, Kiter);

            Aiter.advance();
            Biter.advance();
            Citer.advance();
            Diter.advance();
            Eiter.advance();
            Fiter.advance();
            Giter.advance();
            Hiter.advance();
            Iiter.advance();
            Jiter.advance();
            Kiter.advance();
        }
    }
}

/*
 * This version applies a stencil to a set of 1D arrays.  Up to 11 arrays
 * may be used.  Any unused arrays are turned into dummyArray objects.
 * Operations on dummyArray objects are translated into no-ops.
 */
template<typename T_stencil, typename T_numtype1, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8, typename T_array9, typename T_array10,
    class T_array11>
void applyStencil_imp(const T_stencil& stencil, Array<T_numtype1,1>& A,
    T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F, 
    T_array7& G, T_array8& H, T_array9& I, T_array10& J, T_array11& K)
{
    checkShapes(A,B,C,D,E,F,G,H,I,J,K);

    // Determine stencil extent
    TinyVector<int,1> minb, maxb;
    getStencilExtent(minb, maxb, stencil, A, B, C, D, E, F, G, H, I, J, K);

    // Now determine the subdomain over which the stencil
    // can be applied without worrying about overrunning the
    // boundaries of the array
    int stencil_lbound0 = minb(0);
    int stencil_ubound0 = maxb(0);

    int lbound0 = minmax::max(A.lbound(0), A.lbound(0) - stencil_lbound0);
    int ubound0 = minmax::min(A.ubound(0), A.ubound(0) - stencil_ubound0);

#if 0
    cout << "Stencil bounds are:" << endl
     << lbound0 << '\t' << ubound0 << endl;
#endif

    // Now do the actual loop
    FastArrayIterator<T_numtype1,1> Aiter(A);
    _bz_typename T_array2::T_iterator Biter(B);
    _bz_typename T_array3::T_iterator Citer(C);
    _bz_typename T_array4::T_iterator Diter(D);
    _bz_typename T_array5::T_iterator Eiter(E);
    _bz_typename T_array6::T_iterator Fiter(F);
    _bz_typename T_array7::T_iterator Giter(G);
    _bz_typename T_array8::T_iterator Hiter(H);
    _bz_typename T_array9::T_iterator Iiter(I);
    _bz_typename T_array10::T_iterator Jiter(J);
    _bz_typename T_array11::T_iterator Kiter(K);

    // Load the strides for the innermost loop
    Aiter.loadStride(0);
    Biter.loadStride(0);
    Citer.loadStride(0);
    Diter.loadStride(0);
    Eiter.loadStride(0);
    Fiter.loadStride(0);
    Giter.loadStride(0);
    Hiter.loadStride(0);
    Iiter.loadStride(0);
    Jiter.loadStride(0);
    Kiter.loadStride(0);

    // Patch from Derrick Bass
    Aiter.moveTo(lbound0);
    Biter.moveTo(lbound0);
    Citer.moveTo(lbound0);
    Diter.moveTo(lbound0);
    Eiter.moveTo(lbound0);
    Fiter.moveTo(lbound0);
    Giter.moveTo(lbound0);
    Hiter.moveTo(lbound0);
    Iiter.moveTo(lbound0);
    Jiter.moveTo(lbound0);
    Kiter.moveTo(lbound0);

    for (int i=lbound0; i <= ubound0; ++i)
    {
        stencil.apply(Aiter, Biter, Citer, Diter, Eiter, Fiter, Giter,
            Hiter, Iiter, Jiter, Kiter);

        Aiter.advance();
        Biter.advance();
        Citer.advance();
        Diter.advance();
        Eiter.advance();
        Fiter.advance();
        Giter.advance();
        Hiter.advance();
        Iiter.advance();
        Jiter.advance();
        Kiter.advance();
    }
}

/*
 * These 11 versions of applyStencil handle from 1 to 11 array parameters.
 * They pad their argument list with enough dummyArray objects to call
 * applyStencil_imp with 11 array parameters.
 */
template<typename T_stencil, typename T_numtype1, int N_rank>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A)
{
    applyStencil_imp(stencil, A, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
    T_array2& B)
{
    applyStencil_imp(stencil, A, B, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
    T_array2& B, T_array3& C)
{
    applyStencil_imp(stencil, A, B, C, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray, _dummyArray,
        _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
    T_array2& B, T_array3& C, T_array4& D)
{
    applyStencil_imp(stencil, A, B, C, D, _dummyArray, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E)
{
    applyStencil_imp(stencil, A, B, C, D, E, _dummyArray,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F)
{
    applyStencil_imp(stencil, A, B, C, D, E, F,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F,
   T_array7& G)
{
    applyStencil_imp(stencil, A, B, C, D, E, F, G,
        _dummyArray, _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F,
   T_array7& G, T_array8& H)
{
    applyStencil_imp(stencil, A, B, C, D, E, F, G, H,
        _dummyArray, _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8, typename T_array9>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F,
   T_array7& G, T_array8& H, T_array9& I)
{
    applyStencil_imp(stencil, A, B, C, D, E, F, G, H, I,
        _dummyArray, _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8, typename T_array9, typename T_array10>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F,
   T_array7& G, T_array8& H, T_array9& I, T_array10& J)
{
    applyStencil_imp(stencil, A, B, C, D, E, F, G, H, I, J,
        _dummyArray);
}

template<typename T_stencil, typename T_numtype1, int N_rank, typename T_array2,
    class T_array3, typename T_array4, typename T_array5, typename T_array6,
    class T_array7, typename T_array8, typename T_array9, typename T_array10,
    class T_array11>
inline void applyStencil(const T_stencil& stencil, Array<T_numtype1,N_rank>& A,
   T_array2& B, T_array3& C, T_array4& D, T_array5& E, T_array6& F,
   T_array7& G, T_array8& H, T_array9& I, T_array10& J, T_array11& K)
{
    applyStencil_imp(stencil, A, B, C, D, E, F, G, H, I, J, K);
}

BZ_NAMESPACE_END

#endif // BZ_ARRAYSTENCIL_CC