Tpetra_Details_shortSort.hpp

Go to the documentation of this file.

/*
// @HEADER
// ***********************************************************************
//
//          Tpetra: Templated Linear Algebra Services Package
//                 Copyright (2008) Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
// @HEADER
*/

#ifndef TPETRA_DETAILS_SHORTSORT_HPP
#define TPETRA_DETAILS_SHORTSORT_HPP


#include "TpetraCore_config.h"
#include "Kokkos_Macros.hpp"
#include <type_traits>

namespace Tpetra {
namespace Details {

// Make sure that the macro defined below wasn't defined somewhere else.
#ifdef TPETRA_DETAILS_SWAP_KEYSANDVALUES
#  error "The TPETRA_DETAILS_SWAP_KEYSANDVALUES macro is already defined."
#endif // TPETRA_DETAILS_SWAP_KEYSANDVALUES

#define TPETRA_DETAILS_SWAP_KEYSANDVALUES( i, j )  \
  if (keys[i] > keys[j]) { \
    const KeyType tmpKey (keys[i]); \
    keys[i] = keys[j]; \
    keys[j] = tmpKey; \
    const ValueType tmpVal (values[i]); \
    values[i] = values[j]; \
    values[j] = tmpVal; \
  }

// Make sure that the macro defined below wasn't defined somewhere else.
#ifdef TPETRA_DETAILS_SWAP_KEYS
#  error "The TPETRA_DETAILS_SWAP_KEYS macro is already defined."
#endif // TPETRA_DETAILS_SWAP_KEYSANDVALUES

#define TPETRA_DETAILS_SWAP_KEYS( i, j )  \
  if (keys[i] > keys[j]) { \
    const KeyType tmpKey (keys[i]); \
    keys[i] = keys[j]; \
    keys[j] = tmpKey; \
  }

template<class KeyType, class ValueType>
KOKKOS_FUNCTION void
shortSortKeysAndValues_2 (KeyType keys[2], ValueType values[2])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  For 2 entries, the sorting network is
  // nearly trivial.
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 1);
}

template<class KeyType>
KOKKOS_FUNCTION void
shortSortKeys_2 (KeyType keys[2])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  For 2 entries, the sorting network is
  // nearly trivial.
  TPETRA_DETAILS_SWAP_KEYS(0, 1);
}

template<class KeyType, class ValueType>
KOKKOS_FUNCTION void
shortSortKeysAndValues_3 (KeyType keys[3], ValueType values[3])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 2);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 2);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 1);
}

template<class KeyType>
KOKKOS_FUNCTION void
shortSortKeys_3 (KeyType keys[3])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYS(1, 2);
  TPETRA_DETAILS_SWAP_KEYS(0, 2);
  TPETRA_DETAILS_SWAP_KEYS(0, 1);
}

template<class KeyType, class ValueType>
KOKKOS_FUNCTION void
shortSortKeysAndValues_4 (KeyType keys[4], ValueType values[4])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 1);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(2, 3);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 2);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 3);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 2);
}

template<class KeyType>
KOKKOS_FUNCTION void
shortSortKeys_4 (KeyType keys[4])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYS(0, 1);
  TPETRA_DETAILS_SWAP_KEYS(2, 3);
  TPETRA_DETAILS_SWAP_KEYS(0, 2);
  TPETRA_DETAILS_SWAP_KEYS(1, 3);
  TPETRA_DETAILS_SWAP_KEYS(1, 2);
}

template<class KeyType, class ValueType>
KOKKOS_FUNCTION void
shortSortKeysAndValues_8 (KeyType keys[8], ValueType values[8])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 1);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(2, 3);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 2);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 3);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 2);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(4, 5);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(6, 7);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(4, 6);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(5, 7);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(5, 6);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(0, 4);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 5);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(1, 4);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(2, 6);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(3, 7);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(3, 6);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(2, 4);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(3, 5);
  TPETRA_DETAILS_SWAP_KEYSANDVALUES(3, 4);
}

template<class KeyType>
KOKKOS_FUNCTION void
shortSortKeys_8 (KeyType keys[8])
{
  // Since this function takes a constant number of entries, I use a
  // sorting network here.  To make the network, I used the generator
  // at
  //
  // http://pages.ripco.net/~jgamble/nw.html
  //
  // "Best" algorithm in this case is the Bose-Nelson Algorithm.
  TPETRA_DETAILS_SWAP_KEYS(0, 1);
  TPETRA_DETAILS_SWAP_KEYS(2, 3);
  TPETRA_DETAILS_SWAP_KEYS(0, 2);
  TPETRA_DETAILS_SWAP_KEYS(1, 3);
  TPETRA_DETAILS_SWAP_KEYS(1, 2);
  TPETRA_DETAILS_SWAP_KEYS(4, 5);
  TPETRA_DETAILS_SWAP_KEYS(6, 7);
  TPETRA_DETAILS_SWAP_KEYS(4, 6);
  TPETRA_DETAILS_SWAP_KEYS(5, 7);
  TPETRA_DETAILS_SWAP_KEYS(5, 6);
  TPETRA_DETAILS_SWAP_KEYS(0, 4);
  TPETRA_DETAILS_SWAP_KEYS(1, 5);
  TPETRA_DETAILS_SWAP_KEYS(1, 4);
  TPETRA_DETAILS_SWAP_KEYS(2, 6);
  TPETRA_DETAILS_SWAP_KEYS(3, 7);
  TPETRA_DETAILS_SWAP_KEYS(3, 6);
  TPETRA_DETAILS_SWAP_KEYS(2, 4);
  TPETRA_DETAILS_SWAP_KEYS(3, 5);
  TPETRA_DETAILS_SWAP_KEYS(3, 4);
}

template<class KeyType, class ValueType, class IndexType>
KOKKOS_FUNCTION void
shellSortKeysAndValues (KeyType keys[],
                        ValueType values[],
                        const IndexType n)
{
  static_assert (std::is_integral<IndexType>::value,
                 "IndexType must be a signed integer type.");
  static_assert (std::is_signed<IndexType>::value,
                 "IndexType must be a signed integer type.  "
                 "This implementation does a count-down loop, "
                 "and may thus loop forever "
                 "if one attempts to use it with unsigned IndexType.");
  constexpr IndexType ZERO = 0;
  IndexType midpoint = n / static_cast<IndexType> (2);

  while (midpoint > ZERO) {
    // Avoid names like "max" in case they collide with macros.
    const IndexType theMax = n - midpoint;
    for (IndexType j = 0; j < theMax; ++j) {
      // IndexType is signed, so it's legit to compare >= 0.
      for (IndexType k = j; k >= 0; k -= midpoint) {
        if (keys[k + midpoint] >= keys[k]) {
          break;
        }
        const KeyType tmpKey = keys[k + midpoint];
        keys[k + midpoint] = keys[k];
        keys[k] = tmpKey;
        const ValueType tmpVal = values[k + midpoint];
        values[k + midpoint] = values[k];
        values[k] = tmpVal;
      }
    }
    midpoint = midpoint / 2;
  }
}

template<class KeyType, class IndexType>
KOKKOS_FUNCTION void
shellSortKeys (KeyType keys[], const IndexType n)
{
  static_assert (std::is_integral<IndexType>::value,
                 "IndexType must be a signed integer type.");
  static_assert (std::is_signed<IndexType>::value,
                 "IndexType must be a signed integer type.  "
                 "This implementation does a count-down loop, "
                 "and may thus loop forever "
                 "if one attempts to use it with unsigned IndexType.");
  constexpr IndexType ZERO = 0;
  IndexType midpoint = n / static_cast<IndexType> (2);

  while (midpoint > ZERO) {
    // Avoid names like "max" in case they collide with macros.
    const IndexType theMax = n - midpoint;
    for (int j = 0; j < theMax; ++j) {
      // IndexType must be signed, so it's legit to compare >= 0.
      for (int k = j; k >= 0; k -= midpoint) {
        if (keys[k + midpoint] >= keys[k]) {
          break;
        }
        const KeyType tmpKey = keys[k + midpoint];
        keys[k + midpoint] = keys[k];
        keys[k] = tmpKey;
      }
    }
    midpoint = midpoint / 2;
  }
}

template<typename KeyType, typename ValueType, typename IndexType>
KOKKOS_FUNCTION void
shellSortKeysAndValues2(KeyType* keys, ValueType* values, IndexType n)
{
  IndexType ngaps = 10;
  //Use Ciura's gap choices
  IndexType gaps[] = {3548, 1577, 701, 301, 132, 57, 23, 10, 4, 1};
  for(IndexType gapIndex = 0; gapIndex < ngaps; gapIndex++)
  {
    auto gap = gaps[gapIndex];
    if(n < gap)
      continue;
    //insertion sort the array {keys[0*gap], keys[1*gap], keys[2*gap], ...}
    for(IndexType gapOffset = 0; gapOffset < gap; gapOffset++)
    {
      for(IndexType i = gap + gapOffset; i < n; i += gap)
      {
        //avoid extra swaps: scan for the final position of keys[i]
        if(keys[i - gap] > keys[i])
        {
          IndexType finalPos = i - gap;
          while(finalPos - gap >= 0 && keys[finalPos - gap] > keys[i])
          {
            finalPos -= gap;
          }
          //save keys/values [i], then shift up all keys/values between finalPos and i-gap (inclusive)
          auto tempKey = keys[i];
          auto tempVal = values[i];
          for(IndexType j = i - gap; j >= finalPos; j -= gap)
          {
            keys[j + gap] = keys[j];
            values[j + gap] = values[j];
          }
          keys[finalPos] = tempKey;
          values[finalPos] = tempVal;
        }
      }
    }
  }
#undef SHELL_SWAP
}

} // namespace Details
} // namespace Tpetra

#endif // TPETRA_DETAILS_SHORTSORT_HPP