source: Sophya/trunk/Poubelle/DPC:FitsIOServer/Blitz/blitz/array/methods.cc@ 2878

/*
 * $Id: methods.cc,v 1.1.1.1 1999-11-26 16:37:07 ansari Exp $
 *
 * $Log: not supported by cvs2svn $
// Revision 1.1.1.1  1999/04/09  17:59:03  ansari
// Creation module DPC/Blitz (blitz 0.4) Reza 09/04/99
//
 * Revision 1.4  1998/03/14 00:04:47  tveldhui
 * 0.2-alpha-05
 *
 * Revision 1.3  1997/08/18 19:13:08  tveldhui
 * Just prior to implementing fastRead() optimization for array
 * expression evaluation.
 *
 * Revision 1.2  1997/08/15 21:14:10  tveldhui
 * Just prior to loop-collapse change
 *
 */

#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

#include <blitz/minmax.h>  // Needed for resizeAndPreserve()

BZ_NAMESPACE(blitz)

template<class P_numtype, int N_rank> template<class T_expr>
Array<P_numtype,N_rank>::Array(_bz_ArrayExpr<T_expr> expr)
{
    BZ_NOT_IMPLEMENTED();

    // Obtain storage order from an operand in the expression
    // (if possible).  Probably best to assume C-style storage,
    // then pass the storage object to the expression for possible
    // modification.

    // Obtain ubounds/lbounds from array operands.  Precondition
    // failure if any bounds missing.

    // Size array.

    // assignment of expression.
}

template<class T_numtype, int N_rank>
Array<T_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<class P_numtype, int N_rank>
_bz_inline2 void Array<P_numtype, N_rank>::computeStrides()
{
    if (N_rank > 1)
    {
      int stride = 1;

      // This flag simplifies the code in the loop, encouraging
      // compile-time computation of strides through constant folding.
      _bz_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;

          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<class T_numtype, int N_rank>
void Array<T_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<class 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
       << "storageContiguous = " << storageContiguous_ << endl;
}

/*
 * Make this array a view of another array's data.
 */
template<class P_numtype, int N_rank>
void Array<P_numtype, N_rank>::reference(Array<P_numtype, N_rank>& array)
{
    storage_ = array.storage_;
    length_ = array.length_;
    stride_ = array.stride_;
    zeroOffset_ = array.zeroOffset_;
    storageContiguous_ = array.storageContiguous_;

    MemoryBlockReference<P_numtype>::changeBlock(array, array.zeroOffset_);

    data_ = array.data_;
}

/*
 * This method is called to allocate memory for a new array.  
 */
template<class 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
    MemoryBlockReference<P_numtype>::newBlock(numElements());

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

    // A new array will always have contiguous storage
    storageContiguous_ = _bz_true;
}

template<class T_numtype, int N_rank>
Array<T_numtype, N_rank> Array<T_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;
    }
}

template<class T_numtype, int N_rank>
void Array<T_numtype, N_rank>::makeUnique()
{
    if (numReferences() > 1)
    {
        T_array tmp = copy();
        reference(tmp);
    }
}

template<class T_numtype, int N_rank>
Array<T_numtype, N_rank> Array<T_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)
{
    T_array B(*this);
    B.transposeSelf(r0,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10);
    return B;
}

template<class T_numtype, int N_rank>
void Array<T_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<class T_numtype, int N_rank>
void Array<T_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<class T_numtype, int N_rank>
void Array<T_numtype, N_rank>::reverseSelf(int rank)
{
    BZPRECONDITION(rank < N_rank);

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

    int adjustment = stride_[rank] * (length_[rank] - 1);
    zeroOffset_ += adjustment;
    data_ += adjustment;
    stride_[rank] *= -1;
}

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

template<class T_numtype, int N_rank> template<class T_numtype2>
Array<T_numtype2,N_rank> Array<T_numtype,N_rank>::extractComponent(T_numtype2, 
    int componentNumber, int numComponents)
{
    BZPRECONDITION((componentNumber >= 0) && (componentNumber < numComponents));

    TinyVector<int,N_rank> stride2;
    stride2 = stride_ * numComponents;
    T_numtype2* dataFirst2 = ((T_numtype2*)dataFirst()) + componentNumber;
    return Array<T_numtype2,N_rank>(dataFirst2, length_, stride2, storage_);
}

BZ_NAMESPACE_END

#endif // BZ_ARRAY_CC