
template<class T>
inline void blockMaxval(T* res, T* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    if(*src > *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*src > *res) *res = *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMaxval(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

template<class T>
inline void blockMaxval(T* res,
                        T* src, int* srcStrs, int* msk, int* mskStrs,
                        int* exts, int rank) {
  if(rank == 0) {
    if(*msk && *src > *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *srcStrs ;
    const int step2 = *mskStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk && *src > *res) *res = *src ;
      src += step1 ;
      msk += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = srcStrs [rank] ;
    const int step2 = mskStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMaxval(res, src, srcStrs, msk, mskStrs, exts, rank) ;
      src += step1 ;
      msk += step2;
    }
  }
}


template<class T>
inline void blockMinval(T* res, T* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    if(*src < *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*src < *res) *res = *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMinval(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

template<class T>
inline void blockMinval(T* res,
                        T* src, int* srcStrs, int* msk, int* mskStrs,
                        int* exts, int rank) {
  if(rank == 0) {
    if(*msk && *src < *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *srcStrs ;
    const int step2 = *mskStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk && *src < *res) *res = *src ;
      src += step1 ;
      msk += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = srcStrs [rank] ;
    const int step2 = mskStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMinval(res, src, srcStrs, msk, mskStrs, exts, rank) ;
      src += step1 ;
      msk += step2;
    }
  }
}


template<class T>
inline void blockSum(T* res, T* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    *res += *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res += *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockSum(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

template<class T>
inline void blockSum(T* res,
                     T* src, int* srcStrs, int* msk, int* mskStrs,
                     int* exts, int rank) {
  if(rank == 0) {
    if(*msk) *res += *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *srcStrs ;
    const int step2 = *mskStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk) *res += *src ;
      src += step1 ;
      msk += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = srcStrs [rank] ;
    const int step2 = mskStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockSum(res, src, srcStrs, msk, mskStrs, exts, rank) ;
      src += step1 ;
      msk += step2;
    }
  }
}

template<class T>
inline void blockProduct(T* res, T* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    *res *= *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res *= *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockProduct(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

template<class T>
inline void blockProduct(T* res,
                         T* src, int* srcStrs, int* msk, int* mskStrs,
                         int* exts, int rank) {
  if(rank == 0) {
    if(*msk) *res *= *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *srcStrs ;
    const int step2 = *mskStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk) *res *= *src ;
      src += step1 ;
      msk += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = srcStrs [rank] ;
    const int step2 = mskStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockProduct(res, src, srcStrs, msk, mskStrs, exts, rank) ;
      src += step1 ;
      msk += step2;
    }
  }
}

inline void blockAnd(int* res, int* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    *res += *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res = *res && *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockAnd(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

inline void blockOr(int* res, int* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    *res += *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res = *res || *src ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockOr(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}

inline void blockCount(int* res, int* src, int* srcStrs, int* exts, int rank) {
  if(rank == 0) {
    *res += *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*src) (*res)++ ;
      src += step ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockCount(res, src, srcStrs, exts, rank) ;
      src += step ;
    }
  }
}


template<class S, class T, class U>
inline void blockDotProduct(S* res,
                            T* src1, int* src1Strs, U* src2, int* src2Strs,
                            int* exts, int rank) {
  if(rank == 0) {
    *res += *src1 * *src2 ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *src1Strs ;
    const int step2 = *src2Strs ;
    for(int i = 0 ; i < count ; i++) {
      *res += *src1 * *src2 ;
      src1 += step1 ;
      src2 += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = src1Strs [rank] ;
    const int step2 = src2Strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockDotProduct(res, src1, src1Strs, src2, src2Strs, exts, rank) ;
      src1 += step1 ;
      src2 += step2;
    }
  }
}

inline void blockBoolDotProduct(int* res,
                                int* src1, int* src1Strs,
                                int* src2, int* src2Strs,
                                int* exts, int rank) {
  if(rank == 0) {
    *res = *res || (*src1 && *src2) ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step1 = *src1Strs ;
    const int step2 = *src2Strs ;
    for(int i = 0 ; i < count ; i++) {
      *res = *res || (*src1 && *src2) ;
      src1 += step1 ;
      src2 += step2;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step1 = src1Strs [rank] ;
    const int step2 = src2Strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockBoolDotProduct(res, src1, src1Strs, src2, src2Strs, exts, rank) ;
      src1 += step1 ;
      src2 += step2;
    }
  }
}


template<class T>
inline void blockCopy(T* dst, int* dstStrs, T* src, int* srcStrs,
                      int* exts, int rank) {
  if(rank == 0)
    *dst = *src ;
  else if(rank == 1) {
    const int count   = *exts ;
    const int dstStep = *dstStrs ;
    const int srcStep = *srcStrs ;
    for(int i = 0 ; i < count ; i++) {
      *dst = *src ;
      dst += dstStep ;
      src += srcStep ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int dstStep = dstStrs [rank] ;
    const int srcStep = srcStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockCopy(dst, dstStrs, src, srcStrs, exts, rank) ;
      dst += dstStep ;
      src += srcStep ;
    }
  }
}

template <class T>
inline T* blockGather(T* dat, int* strs, int* exts, int rank, T* buffer) {
  if(rank == 0) {
    *buffer = *dat ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *buffer = *dat ;
      dat    += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockGather(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

template <class T>
inline T* blockScatter(T* dat, int* strs, int* exts, int rank, T* buffer) {
  if(rank == 0) {
    *dat = *buffer ;
    buffer++ ;
  }
  if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *dat   = *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatter(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

template<class T>
inline void blockVecSum(T* res, int* resStrs, T* src, int* srcStrs,
		        int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecSum(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      *res += *src ;
      src  += srcStep ;
    }
  }
}

template<class T>
inline void blockVecSum(T* res, int* resStrs,
                        T* src, int* srcStrs, int* msk, int* mskStrs,
		        int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecSum(res, resStrs, src, srcStrs, msk, mskStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*msk) *res += *src ;
      src  += srcStep ;
      msk  += mskStep ;
    }
  }
}

template <class T>
inline T* blockScatterSum(T* dat, int* strs, int* exts, int rank, T* buffer) {
  if(rank == 0) {
    *dat  += *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *dat  += *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterSum(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}


template<class W>
inline void blockZero(W* res, int* resStrs, int* exts, int rank) {
  if(rank == 0)
    *res = 0 ;
  else if(rank == 1) {
    const int count   = *exts ;
    const int resStep = *resStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res = 0 ;
      res += resStep ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int resStep = resStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockZero(res, resStrs, exts, rank) ;
      res += resStep ;
    }
  }
}

template<class T>
inline void blockVecProduct(T* res, int* resStrs, T* src, int* srcStrs,
		            int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecProduct(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      *res *= *src ;
      src  += srcStep ;
    }
  }
}

template<class T>
inline void blockVecProduct(T* res, int* resStrs,
                            T* src, int* srcStrs, int* msk, int* mskStrs,
		            int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecProduct(res, resStrs, src, srcStrs, msk, mskStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*msk) *res *= *src ;
      src  += srcStep ;
      msk  += mskStep ;
    }
  }
}


template <class T>
inline T* blockScatterProduct(T* dat, int* strs, int* exts, int rank,
                              T* buffer) {
  if(rank == 0) {
    *dat  *= *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *dat  *= *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterProduct(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

template<class T>
inline void blockOne(T* res, int* resStrs, int* exts, int rank) {
  if(rank == 0)
    *res = 1 ;
  else if(rank == 1) {
    const int count   = *exts ;
    const int resStep = *resStrs ;
    for(int i = 0 ; i < count ; i++) {
      *res = 1 ;
      res += resStep ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int resStep = resStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockOne(res, resStrs, exts, rank) ;
      res += resStep ;
    }
  }
}


template<class T>
inline void blockVecMaxval(T* res, int* resStrs, T* src, int* srcStrs,
		           int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMaxval(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*src > *res) *res = *src ;
      src  += srcStep ;
    }
  }
}

template<class T>
inline void blockVecMaxval(T* res, int* resStrs,
                           T* src, int* srcStrs, int* msk, int* mskStrs,
		           int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMaxval(res, resStrs, src, srcStrs, msk, mskStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*msk && *src > *res) *res = *src ;
      src  += srcStep ;
      msk  += mskStep ;
    }
  }
}

template <class T>
inline T* blockScatterMaxval(T* dat, int* strs, int* exts, int rank, T* buffer) {
  if(rank == 0) {
    if(*buffer > *dat) *dat = *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      if(*buffer > *dat) *dat = *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterMaxval(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}


template<class T>
inline void blockVecMinval(T* res, int* resStrs, T* src, int* srcStrs,
		           int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMinval(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*src < *res) *res = *src ;
      src  += srcStep ;
    }
  }
}

template<class T>
inline void blockVecMinval(T* res, int* resStrs,
                           T* src, int* srcStrs, int* msk, int* mskStrs,
		           int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMinval(res, resStrs, src, srcStrs, msk, mskStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*msk && *src < *res) *res = *src ;
      src  += srcStep ;
      msk  += mskStep ;
    }
  }
}

template <class T>
inline T* blockScatterMinval(T* dat, int* strs, int* exts, int rank, T* buffer) {
  if(rank == 0) {
    if(*buffer > *dat) *dat = *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      if(*buffer < *dat) *dat = *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterMinval(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

inline void blockVecAnd(int* res, int* resStrs, int* src, int* srcStrs,
		        int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecAnd(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      *res = *res && *src ;
      src  += srcStep ;
    }
  }
}

inline int* blockScatterAnd(int* dat, int* strs, int* exts, int rank,
                            int* buffer) {
  if(rank == 0) {
    *dat  = *dat && *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *dat  = *dat && *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterAnd(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

inline void blockVecOr(int* res, int* resStrs, int* src, int* srcStrs,
		        int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecOr(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      *res = *res || *src ;
      src  += srcStep ;
    }
  }
}

inline int* blockScatterOr(int* dat, int* strs, int* exts, int rank,
                            int* buffer) {
  if(rank == 0) {
    *dat  = *dat || *buffer ;
    buffer++ ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *strs ;
    for(int i = 0 ; i < count ; i++) {
      *dat  = *dat || *buffer ;
      dat   += step ;
      buffer++ ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = strs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      buffer = blockScatterOr(dat, strs, exts, rank, buffer) ;
      dat   += step ;
    }
  }
  return buffer ;
}

inline void blockVecCount(int* res, int* resStrs, int* src, int* srcStrs,
		          int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int resStep = *resStrs ;
    rank-- ;
    exts++ ;
    resStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecCount(res, resStrs, src, srcStrs, exts, rank) ;
      res += resStep ;
      src += srcStep ;
    }
  }
  else {
    for(int i = 0 ; i < count ; i++) {
      if(*src) (*res)++ ;
      src  += srcStep ;
    }
  }
}


template <class T>
inline void blockScatterMaxtag(T* val, int* valStrs, int* tag, int* tagStrs,
			       int* exts, int rank,
			       T* &bufVal, int* &bufTag) {
  if(rank == 0) {
    if(*bufVal > *val || (*bufVal == *val && *bufTag < *tag)) {
      *val = *bufVal ;
      *tag = *bufTag ;
    }
    bufVal++ ;
    bufTag++ ;
  }
  else if(rank == 1) {
    const int count   = *exts ;
    const int valStep = *valStrs ;
    const int tagStep = *tagStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*bufVal > *val || (*bufVal == *val && *bufTag < *tag)) {
	*val = *bufVal ;
	*tag = *bufTag ;
      }
      val += valStep ;
      tag += tagStep ;
      bufVal++ ;
      bufTag++ ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int valStep = valStrs [rank] ;
    const int tagStep = tagStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockScatterMaxtag(val, valStrs, tag, tagStrs,
                         exts, rank, bufVal, bufTag) ;
      val += valStep ;
      tag += tagStep ;
    }
  }
}


template <class T>
inline void blockScatterMintag(T* val, int* valStrs, int* tag, int* tagStrs,
			       int* exts, int rank,
			       T* &bufVal, int* &bufTag) {
  if(rank == 0) {
    if(*bufVal < *val || (*bufVal == *val && *bufTag < *tag)) {
      *val = *bufVal ;
      *tag = *bufTag ;
    }
    bufVal++ ;
    bufTag++ ;
  }
  else if(rank == 1) {
    const int count   = *exts ;
    const int valStep = *valStrs ;
    const int tagStep = *tagStrs ;
    for(int i = 0 ; i < count ; i++) {
      if(*bufVal < *val || (*bufVal == *val && *bufTag < *tag)) {
	*val = *bufVal ;
	*tag = *bufTag ;
      }
      val += valStep ;
      tag += tagStep ;
      bufVal++ ;
      bufTag++ ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int valStep = valStrs [rank] ;
    const int tagStep = tagStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockScatterMintag(val, valStrs, tag, tagStrs,
                         exts, rank, bufVal, bufTag) ;
      val += valStep ;
      tag += tagStep ;
    }
  }
}


template<class T>
inline void blockMostneg(T* res, int* resStrs, int* exts, int rank) {
  if(rank == 0)
    mostneg(res) ;
  else if(rank == 1) {
    const int count   = *exts ;
    const int resStep = *resStrs ;
    for(int i = 0 ; i < count ; i++) {
      mostneg(res) ;
      res += resStep ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int resStep = resStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMostneg(res, resStrs, exts, rank) ;
      res += resStep ;
    }
  }
}

template<class T>
inline void blockMostpos(T* res, int* resStrs, int* exts, int rank) {
  if(rank == 0)
    mostpos(res) ;
  else if(rank == 1) {
    const int count   = *exts ;
    const int resStep = *resStrs ;
    for(int i = 0 ; i < count ; i++) {
      mostpos(res) ;
      res += resStep ;
    }
  }
  else {
    rank-- ;
    const int count   = exts [rank] ;
    const int resStep = resStrs [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMostpos(res, resStrs, exts, rank) ;
      res += resStep ;
    }
  }
}

inline int le(int* a, int* b, const int n) {

  // Vectorwise a <= b, reverse-lexicographic ordering.

  for(int i = n - 1 ; i >= 0 ; i--)
    if(a [i] > b [i]) return 0 ;
  return 1 ;
}

template<class T>
inline void blockMaxloc(T* res, int* tags, const int tagsRank,
                        T* src, int* srcStrs,
                        int* blockBases, int* blockSteps,
			int* exts, int rank) {
  if(rank == 0) {
    if(*src > *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    const int idxStep = *blockSteps ;
    int idx = *blockBases ;

    for(int i = 0 ; i < count ; i++) {
      if(*src > *res ||
            (*src == *res &&
             le(blockBases + 1, tags + 1, tagsRank - 1) && idx < *tags)) {
        *res  = *src ;
        *tags = idx  ;
        for(int r = 1 ; r < tagsRank ; r++)
          tags [r] = blockBases [r] ;
              // Higher elements of `blockBases' hold higher `idxs' values.
      }
      src += step ;
      idx += idxStep ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    const int idxStep = blockSteps [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMaxloc(res, tags, tagsRank, src, srcStrs, blockBases, blockSteps,
		  exts, rank) ;
      src += step ;
      blockBases [rank] += idxStep ;        // Overwrite with current `idx'
    }
    blockBases [rank] -= count * idxStep ;  // Restore to original value.
  }
}

template<class T>
inline void blockMaxloc(T* res, int* tags, const int tagsRank,
                        T* src, int* srcStrs,
                        int* msk, int* mskStrs,
                        int* blockBases, int* blockSteps,
			int* exts, int rank) {
  if(rank == 0) {
    if(*msk && *src > *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int srcStep = *srcStrs ;
    const int mskStep = *srcStrs ;
    const int idxStep = *blockSteps ;
    int idx = *blockBases ;

    for(int i = 0 ; i < count ; i++) {
      if(*msk) {
	if(*src > *res ||
	      (*src == *res &&
	       le(blockBases + 1, tags + 1, tagsRank - 1) && idx < *tags)) {
	  *res  = *src ;
	  *tags = idx  ;
	  for(int r = 1 ; r < tagsRank ; r++)
	    tags [r] = blockBases [r] ;
		// Higher elements of `blockBases' hold higher `idxs' values.
	}
      }
      src += srcStep ;
      msk += mskStep ;
      idx += idxStep ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int srcStep = srcStrs [rank] ;
    const int mskStep = mskStrs [rank] ;
    const int idxStep = blockSteps [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMaxloc(res, tags, tagsRank, src, srcStrs, msk, mskStrs,
                  blockBases, blockSteps, exts, rank) ;
      src += srcStep ;
      msk += mskStep ;
      blockBases [rank] += idxStep ;        // Overwrite with current `idx'
    }
    blockBases [rank] -= count * idxStep ;  // Restore to original value.
  }
}

template<class T>
inline void blockMinloc(T* res, int* tags, const int tagsRank,
                        T* src, int* srcStrs,
                        int* blockBases, int* blockSteps,
			int* exts, int rank) {
  if(rank == 0) {
    if(*src < *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int step  = *srcStrs ;
    const int idxStep = *blockSteps ;
    int idx = *blockBases ;

    for(int i = 0 ; i < count ; i++) {
      if(*src < *res ||
            (*src == *res &&
             le(blockBases + 1, tags + 1, tagsRank - 1) && idx < *tags)) {
        *res  = *src ;
        *tags = idx  ;
        for(int r = 1 ; r < tagsRank ; r++)
          tags [r] = blockBases [r] ;
              // Higher elements of `blockBases' hold higher `idxs' values.
      }
      src += step ;
      idx += idxStep ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int step  = srcStrs [rank] ;
    const int idxStep = blockSteps [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMinloc(res, tags, tagsRank, src, srcStrs, blockBases, blockSteps,
		  exts, rank) ;
      src += step ;
      blockBases [rank] += idxStep ;        // Overwrite with current `idx'
    }
    blockBases [rank] -= count * idxStep ;  // Restore to original value.
  }
}

template<class T>
inline void blockMinloc(T* res, int* tags, const int tagsRank,
                        T* src, int* srcStrs,
                        int* msk, int* mskStrs,
                        int* blockBases, int* blockSteps,
			int* exts, int rank) {
  if(rank == 0) {
    if(*msk && *src < *res) *res = *src ;
  }
  else if(rank == 1) {
    const int count = *exts ;
    const int srcStep = *srcStrs ;
    const int mskStep = *srcStrs ;
    const int idxStep = *blockSteps ;
    int idx = *blockBases ;

    for(int i = 0 ; i < count ; i++) {
      if(*msk) {
	if(*src < *res ||
	      (*src == *res &&
	       le(blockBases + 1, tags + 1, tagsRank - 1) && idx < *tags)) {
	  *res  = *src ;
	  *tags = idx  ;
	  for(int r = 1 ; r < tagsRank ; r++)
	    tags [r] = blockBases [r] ;
		// Higher elements of `blockBases' hold higher `idxs' values.
	}
      }
      src += srcStep ;
      msk += mskStep ;
      idx += idxStep ;
    }
  }
  else {
    rank-- ;
    const int count = exts [rank] ;
    const int srcStep = srcStrs [rank] ;
    const int mskStep = mskStrs [rank] ;
    const int idxStep = blockSteps [rank] ;
    for(int i = 0 ; i < count ; i++) {
      blockMinloc(res, tags, tagsRank, src, srcStrs, msk, mskStrs,
                  blockBases, blockSteps, exts, rank) ;
      src += srcStep ;
      msk += mskStep ;
      blockBases [rank] += idxStep ;        // Overwrite with current `idx'
    }
    blockBases [rank] -= count * idxStep ;  // Restore to original value.
  }
}

template <class T>
inline void blockVecMaxloc(T* val, int* valStrs, int* tag, int* tagStrs,
                           T* src, int* srcStrs,
                           const int blockBase, const int blockStep,
			   int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int valStep = *valStrs ;
    int tagStep = *tagStrs ;
    rank-- ;
    exts++ ;
    valStrs++ ;
    tagStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMaxloc(val, valStrs, tag, tagStrs,
                     src, srcStrs, blockBase, blockStep,
                     exts, rank) ;
      val += valStep ;
      tag += tagStep ;
      src += srcStep ;
    }
  }
  else {
    int loc = blockBase ;
    for(int i = 0 ; i < count ; i++) {
      if(*src > *val || (*src == *val && loc < *tag)) {
	*val = *src ;
	*tag = loc ;
      }
      src += srcStep ;
      loc += blockStep ;
    }
  }
}

template <class T>
inline void blockVecMaxloc(T* val, int* valStrs, int* tag, int* tagStrs,
                           T* src, int* srcStrs, int* msk, int* mskStrs,
                           const int blockBase, const int blockStep,
			   int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int valStep = *valStrs ;
    int tagStep = *tagStrs ;
    rank-- ;
    exts++ ;
    valStrs++ ;
    tagStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMaxloc(val, valStrs, tag, tagStrs,
                     src, srcStrs, msk, mskStrs, blockBase, blockStep,
                     exts, rank) ;
      val += valStep ;
      tag += tagStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    int loc = blockBase ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk && (*src > *val || (*src == *val && loc < *tag))) {
	*val = *src ;
	*tag = loc ;
      }
      src += srcStep ;
      msk += mskStep ;
      loc += blockStep ;
    }
  }
}

template <class T>
inline void blockVecMinloc(T* val, int* valStrs, int* tag, int* tagStrs,
                           T* src, int* srcStrs,
                           const int blockBase, const int blockStep,
			   int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs ;
  if(rank) {
    int valStep = *valStrs ;
    int tagStep = *tagStrs ;
    rank-- ;
    exts++ ;
    valStrs++ ;
    tagStrs++ ;
    srcStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMinloc(val, valStrs, tag, tagStrs,
                     src, srcStrs, blockBase, blockStep,
                     exts, rank) ;
      val += valStep ;
      tag += tagStep ;
      src += srcStep ;
    }
  }
  else {
    int loc = blockBase ;
    for(int i = 0 ; i < count ; i++) {
      if(*src < *val || (*src == *val && loc < *tag)) {
	*val = *src ;
	*tag = loc ;
      }
      src += srcStep ;
      loc += blockStep ;
    }
  }
}

template <class T>
inline void blockVecMinloc(T* val, int* valStrs, int* tag, int* tagStrs,
                           T* src, int* srcStrs, int* msk, int* mskStrs,
                           const int blockBase, const int blockStep,
			   int* exts, int rank) {

  int count = *exts, srcStep = *srcStrs, mskStep = *mskStrs ;
  if(rank) {
    int valStep = *valStrs ;
    int tagStep = *tagStrs ;
    rank-- ;
    exts++ ;
    valStrs++ ;
    tagStrs++ ;
    srcStrs++ ;
    mskStrs++ ;
    for(int j = 0 ; j < count ; j++) {
      blockVecMinloc(val, valStrs, tag, tagStrs,
                     src, srcStrs, msk, mskStrs, blockBase, blockStep,
                     exts, rank) ;
      val += valStep ;
      tag += tagStep ;
      src += srcStep ;
      msk += mskStep ;
    }
  }
  else {
    int loc = blockBase ;
    for(int i = 0 ; i < count ; i++) {
      if(*msk && (*src < *val || (*src == *val && loc < *tag))) {
	*val = *src ;
	*tag = loc ;
      }
      src += srcStep ;
      msk += mskStep ;
      loc += blockStep ;
    }
  }
}


template<class S, class T, class U>
inline void blockMatmul(S* res,  int* resStrs,
                        T* src1, int* src1Strs, U* src2, int* src2Strs,
                        const int m, const int n, const int l) {

  // res : m by n, src1 : m by l, src2 : l by n.

  // Should really use an optimized sequential matrix multiplier here.

  const int resStep0  = resStrs  [0] ;
  const int resStep1  = resStrs  [1] ;
  const int src1Step0 = src1Strs [0] ;
  const int src1Step1 = src1Strs [1] ;
  const int src2Step0 = src2Strs [0] ;
  const int src2Step1 = src2Strs [1] ;

  for(int i = 0 ; i < m ; i++) {
    S* resPtr   = res  ;
    U* src2Ptr1 = src2 ;
    for(int j = 0 ; j < n ; j++) {
      T* src1Ptr  = src1 ;
      U* src2Ptr2 = src2Ptr1 ;
      for(int k = 0 ; k < l ; k++) {
        *resPtr  += *src1Ptr * *src2Ptr2 ;
        src1Ptr  += src1Step1 ;
        src2Ptr2 += src2Step0 ;
      }
      resPtr   += resStep1  ;
      src2Ptr1 += src2Step1 ;
    }
    res  += resStep0  ;
    src1 += src1Step0 ;
  }
}

// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //
//                                                                     //
// `blockCopy', `blockGather', `blockScatter': Untyped versions.       //
//                                                                     //
// - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //


void blockCopy(char* dst, int* dstStrs, char* src, int* srcStrs,
               int* exts, int rank, int elementLen) ;

char* blockGather(char* dat, int* strs, int* exts, int rank,
                  char* buffer, int elementLen) ;

char* blockScatter(char* dat, int* strs, int* exts, int rank,
                   char* buffer, int elementLen) ;



// -------------------------------------------------------------
// Things to do
//
//   Eliminate return value of `blockGather'...  Make `buffer' argument
//   a reference?
//
//   Find out why order of loops in, eg, `blockSum' is different to
//   order in, eg, `blockVecSum'.  Does it matter?  Which is best?
//
//   In general, these functions need organizing/rationalizing/optimizing.
//   Probably this file should be generated by a macro preprocessor.

