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

/***************************************************************************
 * blitz/array/methods.cc   General array class methods.
 *
 * $Id$
 *
 * Copyright (C) 1997-2011 Todd Veldhuizen <tveldhui@acm.org>
 *
 * This file is a part of Blitz.
 *
 * Blitz is free software: you can redistribute it and/or modify 
 * it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation, either version 3
 * of the License, or (at your option) any later version.
 *
 * Blitz 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public 
 * License along with Blitz.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Suggestions:          blitz-devel@lists.sourceforge.net
 * Bugs:                 blitz-support@lists.sourceforge.net    
 *
 * For more information, please see the Blitz++ Home Page:
 *    https://sourceforge.net/projects/blitz/
 *
 ****************************************************************************/
#ifndef BZ_ARRAYMETHODS_CC
#define BZ_ARRAYMETHODS_CC

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

BZ_NAMESPACE(blitz)

template<typename P_numtype, int N_rank> template<typename T_expr>
Array<P_numtype,N_rank>::Array(_bz_ArrayExpr<T_expr> expr)
{
    // Determine extent of the array expression

    TinyVector<int,N_rank> lbound, extent, ordering;
    TinyVector<bool,N_rank> ascendingFlag;
    TinyVector<bool,N_rank> in_ordering;
    in_ordering = false;

    int j = 0;
    for (int i=0; i < N_rank; ++i)
    {
        lbound(i) = expr.lbound(i);
        int ubound = expr.ubound(i);
        extent(i) = ubound - lbound(i) + 1;
        int orderingj = expr.ordering(i);
        if (orderingj != INT_MIN && orderingj < N_rank &&
            !in_ordering( orderingj )) { // unique value in ordering array
            in_ordering( orderingj ) = true;
            ordering(j++) = orderingj;
        }
        int ascending = expr.ascending(i);
        ascendingFlag(i) = (ascending == 1);

#ifdef BZ_DEBUG
        if ((lbound(i) == INT_MIN) || (ubound == INT_MAX) 
          || (ordering(i) == INT_MIN) || (ascending == INT_MIN))
        {
          BZPRECHECK(0,
           "Attempted to construct an array from an expression " << endl
           << "which does not have a shape.  To use this constructor, "
           << endl 
           << "the expression must contain at least one array operand.");
          return;
        }
#endif
    }

    // It is possible that ordering is not a permutation of 0,...,N_rank-1.
    // In that case j will be less than N_rank. We fill in ordering with the
    // usused values in decreasing order.
    for (int i = N_rank-1; j < N_rank; ++j) {
        while (in_ordering(i))
          --i;
        ordering(j) = i--;
    }

    Array<T_numtype,N_rank> A(lbound,extent,
        GeneralArrayStorage<N_rank>(ordering,ascendingFlag));
    A = expr;
    reference(A);
}

template<typename P_numtype, int N_rank>
Array<P_numtype,N_rank>::Array(const TinyVector<int, N_rank>& lbounds,
    const TinyVector<int, N_rank>& extent,
    const GeneralArrayStorage<N_rank>& storage)
    : storage_(storage)
{
    length_ = extent;
    storage_.setBase(lbounds);
    setupStorage(N_rank - 1);
}


/** This routine takes the storage information for the array
    (ascendingFlag_[], base_[], and ordering_[]) and the size of the
    array (length_[]) and computes the stride vector (stride_[]) and
    the zero offset (see explanation in array.h).
 */
template<typename P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::computeStrides()
{
    if (N_rank > 1)
    {
      diffType stride = 1;

      // This flag simplifies the code in the loop, encouraging
      // compile-time computation of strides through constant folding.
      bool allAscending = storage_.allRanksStoredAscending();

      // BZ_OLD_FOR_SCOPING
      int n;
      for (n=0; n < N_rank; ++n)
      {
          int strideSign = +1;

          // If this rank is stored in descending order, then the stride
          // will be negative.
          if (!allAscending)
          {
            if (!isRankStoredAscending(ordering(n)))
                strideSign = -1;
          }

          // The stride for this rank is the product of the lengths of
          // the ranks minor to it.
          stride_[ordering(n)] = stride * strideSign;

          if((storage_.padding()==paddedData)&&(n==0)) {
            // The lowest rank dimension is padded to vecWidth, so this
            // needs to be accounted for in the stride
            stride *= simdTypes<T_numtype>::paddedLength(length_[ordering(0)]);
          }
          else
            stride *= length_[ordering(n)];
      }
    }
    else {
        // Specialization for N_rank == 1
        // This simpler calculation makes it easier for the compiler
        // to propagate stride values.

        if (isRankStoredAscending(0))
            stride_[0] = 1;
        else
            stride_[0] = -1;
    }

    calculateZeroOffset();
}

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::calculateZeroOffset()
{
    // Calculate the offset of (0,0,...,0)
    zeroOffset_ = 0;

    // zeroOffset_ = - sum(where(ascendingFlag_, stride_ * base_,
    //     (length_ - 1 + base_) * stride_))
    for (int n=0; n < N_rank; ++n)
    {
        if (!isRankStoredAscending(n))
            zeroOffset_ -= (length_[n] - 1 + base(n)) * stride_[n];
        else
            zeroOffset_ -= stride_[n] * base(n);
    }
}

template<typename P_numtype, int N_rank>
bool Array<P_numtype, N_rank>::isStorageContiguous() const
{
    // The storage is contiguous if for the set
    // { | stride[i] * extent[i] | }, i = 0..N_rank-1,
    // there is only one value which is not in the set
    // of strides; and if there is one stride which is 1.

    // This algorithm is quadratic in the rank.  It is hard
    // to imagine this being a serious problem.

    int numStridesMissing = 0;
    bool haveUnitStride = false;

    for (int i=0; i < N_rank; ++i)
    {
      diffType stride = BZ_MATHFN_SCOPE(abs)(stride_[i]);
        if (stride == 1)
            haveUnitStride = true;

        diffType vi = stride * length_[i];

        int j = 0;
        for (j=0; j < N_rank; ++j)
            if (BZ_MATHFN_SCOPE(abs)(stride_[j]) == vi)
                break;

        if (j == N_rank)
        {
            ++numStridesMissing;
            if (numStridesMissing == 2)
                return false;
        }
    }

    return haveUnitStride;
}

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::dumpStructureInformation(ostream& os) const
{
    os << "Dump of Array<" << BZ_DEBUG_TEMPLATE_AS_STRING_LITERAL(P_numtype) 
       << ", " << N_rank << ">:" << endl
       << "ordering_      = " << storage_.ordering() << endl
       << "ascendingFlag_ = " << storage_.ascendingFlag() << endl
       << "base_          = " << storage_.base() << endl
       << "length_        = " << length_ << endl
       << "stride_        = " << stride_ << endl
       << "zeroOffset_    = " << zeroOffset_ << endl
       << "numElements()  = " << numElements() << endl
       << "isStorageContiguous() = " << isStorageContiguous() << endl;
}

/**
  Make this array a view of another array's data. This overrides the
  current storage of the array.
 */
template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::reference(const Array<P_numtype, N_rank>& array)
{
    storage_ = array.storage_;
    length_ = array.length_;
    stride_ = array.stride_;
    zeroOffset_ = array.zeroOffset_;

    T_base::changeBlock(array.noConst());
}

/** This method makes the array reference another, but it does it as a
    "weak" reference that is not counted. If you can guarantee that
    the array memory block containing the data is persistent, this
    will allow reference counting to be bypassed for this array, which
    if mutex-locking is involved is a significant overhead. */
template<typename P_numtype, int N_rank>
void 
Array<P_numtype, N_rank>::weakReference(const Array<P_numtype, N_rank>& array)
{
    storage_ = array.storage_;
    length_ = array.length_;
    stride_ = array.stride_;
    zeroOffset_ = array.zeroOffset_;

    T_base::changeToNullBlock();
    data_ = array.data_;
}


/**
   Modify the Array storage.  Array must be unallocated.
 */
template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::setStorage(GeneralArrayStorage<N_rank> x)
{
#ifdef BZ_DEBUG
    if (size() != 0) {
        BZPRECHECK(0,"Cannot modify storage format of an Array that has already been allocated!" << endl);
        return;
    }
#endif
    storage_ = x;
    return;
}

/**
   This method is called to allocate memory for a new array. It
   assumes the storage_ and length_ members are already initialized.
 */
template<typename P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::setupStorage(int lastRankInitialized)
{
    TAU_TYPE_STRING(p1, "Array<T,N>::setupStorage() [T="
        + CT(P_numtype) + ",N=" + CT(N_rank) + "]");
    TAU_PROFILE(" ", p1, TAU_BLITZ);

    /*
     * If the length of some of the ranks was unspecified, fill these
     * in using the last specified value.
     *
     * e.g. Array<int,3> A(40) results in a 40x40x40 array.
     */
    for (int i=lastRankInitialized + 1; i < N_rank; ++i)
    {
        storage_.setBase(i, storage_.base(lastRankInitialized));
        length_[i] = length_[lastRankInitialized];
    }

    // Compute strides
    computeStrides();

    // Allocate a block of memory.
    TinyVector<int, N_rank> alloc_length = length();
    if(storage_.padding()==paddedData) {
      // The size of the block is NOT equal to numelements, because the
      // lowest rank dimension is padded to vecWidth
      alloc_length[ordering(0)] = 
        simdTypes<T_numtype>::paddedLength(alloc_length[ordering(0)]);
    }
    sizeType numElem = _bz_returntype<sizeType>::product(alloc_length);
    if (numElem==0)
        T_base::changeToNullBlock();
    else
        T_base::newBlock(numElem);

    // Adjust the base of the array to account for non-zero base
    // indices and reversals
    data_ += zeroOffset_;
}

/** Return a deep copy of an array (as opposed to the reference copy
    done by the copy constructor. */
template<typename P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype, N_rank>::copy() const
{
    if (numElements())
    {
        Array<T_numtype, N_rank> z(length_, storage_);
        z = *this;
        return z;
    }
    else {
        // Null array-- don't bother allocating an empty block.
        return *this;
    }
}

/** Make the array have its own memory block by making a copy if the
    block has a reference count greater than one. */
template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::makeUnique()
{
    if (T_base::numReferences() > 1)
    {
        T_array tmp = copy();
        reference(tmp);
    }
}

template<typename P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype, N_rank>::transpose(int r0, int r1, 
    int r2, int r3, int r4, int r5, int r6, int r7, int r8, int r9, int r10) const
{
    T_array B(*this);
    B.transposeSelf(r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10);
    return B;
}

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::transposeSelf(int r0, int r1, int r2, int r3,
    int r4, int r5, int r6, int r7, int r8, int r9, int r10)
{
    BZPRECHECK(r0+r1+r2+r3+r4+r5+r6+r7+r8+r9+r10 == N_rank * (N_rank-1) / 2,
        "Invalid array transpose() arguments." << endl
        << "Arguments must be a permutation of the numerals (0,...,"
        << (N_rank - 1) << ")");

    // Create a temporary reference copy of this array
    Array<T_numtype, N_rank> x(*this);

    // Now reorder the dimensions using the supplied permutation
    doTranspose(0, r0, x);
    doTranspose(1, r1, x);
    doTranspose(2, r2, x);
    doTranspose(3, r3, x);
    doTranspose(4, r4, x);
    doTranspose(5, r5, x);
    doTranspose(6, r6, x);
    doTranspose(7, r7, x);
    doTranspose(8, r8, x);
    doTranspose(9, r9, x);
    doTranspose(10, r10, x);
}

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::doTranspose(int destRank, int sourceRank,
    Array<T_numtype, N_rank>& array)
{
    // BZ_NEEDS_WORK: precondition check

    if (destRank >= N_rank)
        return;

    length_[destRank] = array.length_[sourceRank];
    stride_[destRank] = array.stride_[sourceRank];
    storage_.setAscendingFlag(destRank, 
        array.isRankStoredAscending(sourceRank));
    storage_.setBase(destRank, array.base(sourceRank));

    // BZ_NEEDS_WORK: Handling the storage ordering is currently O(N^2)
    // but it can be done fairly easily in linear time by constructing
    // the appropriate permutation.

    // Find sourceRank in array.storage_.ordering_
    int i=0;
    for (; i < N_rank; ++i)
        if (array.storage_.ordering(i) == sourceRank)
            break;

    storage_.setOrdering(i, destRank);
}

template<typename P_numtype, int N_rank>
void Array<P_numtype, N_rank>::reverseSelf(int rank)
{
    BZPRECONDITION(rank < N_rank);

    storage_.setAscendingFlag(rank, !isRankStoredAscending(rank));

    diffType adjustment = static_cast<ptrdiff_t>(stride_[rank]) * (lbound(rank) + ubound(rank));
    zeroOffset_ += adjustment;
    data_ += adjustment;
    stride_[rank] *= -1;
}

template<typename P_numtype, int N_rank>
Array<P_numtype, N_rank> Array<P_numtype,N_rank>::reverse(int rank)
{
    T_array B(*this);
    B.reverseSelf(rank);
    return B;
}

template<typename P_numtype, int N_rank> template<typename P_numtype2>
Array<P_numtype2,N_rank> Array<P_numtype,N_rank>::extractComponent(P_numtype2, 
    int componentNumber, int numComponents) const
{
    BZPRECONDITION((componentNumber >= 0) 
        && (componentNumber < numComponents));

    // If P_numtype is a multicomponent type, it may have an alignment
    // setting. For this reason it is not correct to use
    // numComponents, we must use sizeof(P_numtype)/sizeof(P_numtype2)
    // instead.
    BZASSERT(sizeof(P_numtype)%sizeof(P_numtype2)==0);

    TinyVector<diffType, N_rank> stride2;
    for (int i=0; i < N_rank; ++i)
      stride2(i) = stride_(i) * sizeof(P_numtype)/sizeof(P_numtype2);
    const P_numtype2* dataFirst2 = 
        ((const P_numtype2*)dataFirst()) + componentNumber;
    return Array<P_numtype2,N_rank>(const_cast<P_numtype2*>(dataFirst2), 
        length_, stride2, storage_);
}

/* 
 * These routines reindex the current array to use a new base vector.
 * The first reindexes the array, the second just returns a reindex view
 * of the current array, leaving the current array unmodified.
 * (Contributed by Derrick Bass)
 */
template<typename P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::reindexSelf(const 
    TinyVector<int, N_rank>& newBase) 
{
  diffType delta = 0;
    for (int i=0; i < N_rank; ++i)
      delta += (base(i) - newBase(i)) * stride_(i);

    data_ += delta;

    // WAS: dot(base() - newBase, stride_);

    storage_.setBase(newBase);
    calculateZeroOffset();
}

template<typename P_numtype, int N_rank>
_bz_inline2 Array<P_numtype, N_rank> 
Array<P_numtype, N_rank>::reindex(const TinyVector<int, N_rank>& newBase) 
{
    T_array B(*this);
    B.reindexSelf(newBase);
    return B;
}

BZ_NAMESPACE_END

#endif // BZ_ARRAY_CC