text:
/*
  Copyright [2017-2021] [IBM Corporation]
  Licensed under the Apache License, Version 2.0 (the "License");
  you may not use this file except in compliance with the License.
  You may obtain a copy of the License at
  http://www.apache.org/licenses/LICENSE-2.0
  Unless required by applicable law or agreed to in writing, software
  distributed under the License is distributed on an "AS IS" BASIS,
  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  See the License for the specific language governing permissions and
  limitations under the License.
*/

#ifndef _API_KVSTORE_H_
#define _API_KVSTORE_H_

#include <api/components.h>
#include <api/registrar_memory_direct.h>
#include <common/byte.h>
#include <common/byte_span.h>
#include <common/errors.h> /* ERROR_BASE */
#include <common/string_view.h>
#include <common/time.h>
#include <gsl/span>
#include <sys/uio.h> /* iovec */

#include <cinttypes> /* PRIx64 */
#include <cstdlib>
#include <functional>
#include <map>
#include <list>
#include <mutex>
#include <vector>

namespace nupm
{
  struct region_descriptor;
}

/* print format for the pool type */
#define PRIxIKVSTORE_POOL_T PRIx64

namespace component
{

#define DECLARE_OPAQUE_TYPE(NAME) \
  struct Opaque_##NAME {          \
    virtual ~Opaque_##NAME() {}   \
  }

/**
 * Key-value interface for pluggable backend (e.g. mapstore, hstore, hstore-cc)
 */
class KVStore : public Registrar_memory_direct {

protected:
  ~KVStore() {}
 private:
  DECLARE_OPAQUE_TYPE(lock_handle);
  template <typename E, typename ... Args>
    static E error_value(
      E e, Args && ...
    ) { return e; }

 public:
  DECLARE_OPAQUE_TYPE(memory_region); /* Buffer_manager::buffer_t need this */
  DECLARE_OPAQUE_TYPE(key);
  DECLARE_OPAQUE_TYPE(pool_iterator);

 public:
  using pool_t          = uint64_t;
  // using memory_handle_t = Opaque_memory_region*;
  using key_t           = Opaque_key*;
  using pool_lock_t     = Opaque_lock_handle*;
  using pool_iterator_t = Opaque_pool_iterator*;
  using byte            = common::byte;

 template <typename C>
    using basic_string_view = common::basic_string_view<C>;
  using string_view = common::string_view;
  using string_view_byte = basic_string_view<byte>;
  using string_view_key = string_view_byte;
  using string_view_value = string_view_byte;

  static constexpr memory_handle_t HANDLE_NONE = nullptr; /* old name */
  static constexpr memory_handle_t MEMORY_HANDLE_NONE = nullptr; /* better name */
  static constexpr key_t           KEY_NONE    = nullptr;

  struct Addr {
    explicit Addr(addr_t addr_) : addr(addr_) {}
    Addr() = delete;
    addr_t addr;
  };

  enum {
    THREAD_MODEL_UNSAFE,
    THREAD_MODEL_SINGLE_PER_POOL,
    THREAD_MODEL_RWLOCK_PER_POOL,
    THREAD_MODEL_MULTI_PER_POOL,
  };

  using flags_t = std::uint32_t ;
  static constexpr flags_t FLAGS_NONE          = 0x0;
  static constexpr flags_t FLAGS_READ_ONLY     = 0x1; /* lock read-only */
  static constexpr flags_t FLAGS_SET_SIZE      = 0x2;
  static constexpr flags_t FLAGS_CREATE_ONLY   = 0x4;  /* only succeed if no existing k-v pair exist */
  static constexpr flags_t FLAGS_DONT_STOMP    = 0x8;  /* do not overwrite existing k-v pair */
  static constexpr flags_t FLAGS_NO_RESIZE     = 0x10; /* if size < existing size, do not resize */
  static constexpr flags_t FLAGS_MAX_VALUE     = 0x10;

  using unlock_flags_t = std::uint32_t;
  static constexpr unlock_flags_t UNLOCK_FLAGS_NONE = 0x0;
  static constexpr unlock_flags_t UNLOCK_FLAGS_FLUSH = 0x1; /* indicates for PM backends to flush */

  static constexpr pool_t POOL_ERROR = 0;

  enum class Capability {
    POOL_DELETE_CHECK, /*< checks if pool is open before allowing delete */
    RWLOCK_PER_POOL,   /*< pools are locked with RW-lock */
    POOL_THREAD_SAFE,  /*< pools can be shared across multiple client threads */
    WRITE_TIMESTAMPS,  /*< support for write timestamping */
  };

  enum class Op_type {
    WRITE, /* copy bytes into memory region */
    ZERO,  /* zero the memory region */
    INCREMENT_UINT64,
    CAS_UINT64,
  };

  enum Attribute : std::uint32_t {
    VALUE_LEN                = 1, /* length of a value associated with key */
    COUNT                    = 2, /* number of objects */
    CRC32                    = 3, /* get CRC32 of a value */
    AUTO_HASHTABLE_EXPANSION = 4, /* set to true if the hash table should expand */
    PERCENT_USED             = 5, /* get percent used pool capacity at current size */
    WRITE_EPOCH_TIME         = 6, /* epoch time at which the key-value pair was last
                                     written or locked with STORE_LOCK_WRITE */
    MEMORY_TYPE              = 7, /* type of memory */
    MEMORY_SIZE              = 8, /* size of pool or store in bytes */
    NUMA_MASK                = 9, /* mask of first 64 numa nodes eligible for mapstore allocation */
  };

  enum {
    MEMORY_TYPE_DRAM        = 0x1,
    MEMORY_TYPE_PMEM_DEVDAX = 0x2,
    MEMORY_TYPE_PMEM_FSDAX  = 0x3,
    MEMORY_TYPE_UNKNOWN     = 0xFF,
  };

  enum lock_type_t {
    STORE_LOCK_NONE  = 0,
    STORE_LOCK_READ  = 1,
    STORE_LOCK_WRITE = 2,
  };

  enum {
    /* see common/errors.h */
    S_MORE           = 2,
    E_KEY_EXISTS     = E_ERROR_BASE - 1,
    E_KEY_NOT_FOUND  = E_ERROR_BASE - 2,
    E_POOL_NOT_FOUND = E_ERROR_BASE - 3,
    E_BAD_ALIGNMENT  = E_ERROR_BASE - 4,
    E_TOO_LARGE      = E_ERROR_BASE - 5, /* -55 */
    E_ALREADY_OPEN   = E_ERROR_BASE - 6,
  };

  std::string strerro(int e)
  {
    static std::map<int, std::string> errs {
      { S_MORE, "MORE" }
      , { E_KEY_EXISTS, "E_KEY_EXISTS" }
      , { E_KEY_NOT_FOUND, "E_KEY_NOT_FOUND" }
      , { E_POOL_NOT_FOUND, "E_POOL_NOT_FOUND" }
      , { E_BAD_ALIGNMENT, "E_BAD_ALIGNMENT" }
      , { E_TOO_LARGE, "E_TOO_LARGE" }
      , { E_ALREADY_OPEN, "E_ALREADY_OPEN" }
    };
    auto it = errs.find(e);
    return it == errs.end() ? ( "non-IKVStore error " + std::to_string(e) ) : it->second;
  }

  class Operation {
    Op_type _type;
    size_t  _offset;

   protected:
    Operation(Op_type type, size_t offset) : _type(type), _offset(offset) {}

   public:
    Op_type type() const noexcept { return _type; }
    size_t  offset() const noexcept { return _offset; }
  };

  class Operation_sized : public Operation {
    size_t _len;

   protected:
    Operation_sized(Op_type type, size_t offset_, size_t len) : Operation(type, offset_), _len(len) {}

   public:
    size_t size() const noexcept { return _len; }
  };

  class Operation_write : public Operation_sized {
    const void* _data;

   public:
    Operation_write(size_t offset, size_t len, const void* data)
        : Operation_sized(Op_type::WRITE, offset, len), _data(data)
    {
    }
    const void* data() const noexcept { return _data; }
  };

  class Operation_zero : public Operation_sized {

   public:
    Operation_zero(size_t offset, size_t len)
        : Operation_sized(Op_type::ZERO, offset, len)
    {
    }
  };

  /**
   * Determine thread safety of the component
   * Check capability of component
   *
   * @param cap Capability type
   *
   * @return THREAD_MODEL_XXXX
   */
  virtual int thread_safety() const = 0;

  /**
   * Check capability of component
   *
   * @param cap Capability type
   *
   * @return THREAD_MODEL_XXXX
   */
  virtual int get_capability(Capability cap) const { return error_value(-1, cap); }

  /**
   * Create an object pool. If the pool exists and the FLAGS_CREATE_ONLY
   * is not provided, then the existing pool will be opened.  If
   * FLAGS_CREATE_ONLY is specified and the pool exists, POOL ERROR will be
   * returned.
   *
   * @param pool_name Unique pool name
   * @param size Size of pool in bytes (for keys,values and metadata)
   * @param flags Creation flags
   * @param expected_obj_count Expected maximum object count (optimization)
   * @param base Optional base address
   *
   * @return Pool handle or POOL_ERROR
   */
  virtual pool_t create_pool(const std::string& pool_name,
                             const size_t       size,
                             flags_t            flags              = 0,
                             uint64_t           expected_obj_count = 0,
                             const Addr         base_addr = Addr{0})
  {
    PERR("create_pool not implemented");
    return error_value(POOL_ERROR, pool_name, size, flags, expected_obj_count, base_addr);
  }

  virtual pool_t create_pool(const std::string& path,
                             const std::string& name,
                             const size_t       size,
                             flags_t            flags              = 0,
                             uint64_t           expected_obj_count = 0,
                             const Addr         base_addr_unused = Addr{0}) __attribute__((deprecated))
  {
    return create_pool(path + name, size, flags, expected_obj_count, base_addr_unused);
  }

  /**
   * Open an existing pool.
   *
   * @param name Name of object pool
   * @param flags Optional flags e.g., FLAGS_READ_ONLY
   * @param base Optional base address
   *
   * @return Pool handle or POOL_ERROR if pool cannot be opened, or flags
   * unsupported
   */
  virtual pool_t open_pool(const std::string& pool_name,
                           flags_t flags = 0,
                           const Addr base_addr_unused = Addr{0})
  {
    return error_value(POOL_ERROR, pool_name, flags, base_addr_unused);
  }

  virtual pool_t open_pool(const std::string& path,
                           const std::string& name,
                           flags_t flags = 0,
                           const Addr base_addr_unused = Addr{0}) __attribute__((deprecated))
  {
    return open_pool(path + name, flags, base_addr_unused);
  }

  /**
   * Close pool handle
   *
   * @param pool Pool handle
   *
   * @return S_OK on success, E_POOL_NOT_FOUND, E_ALREADY_OPEN if pool cannot be
   * closed due to open session.
   */
  virtual status_t close_pool(pool_t pool) = 0;

  /**
   * Delete an existing pool
   *
   * @param name Name of object pool
   *
   * @return S_OK on success, E_POOL_NOT_FOUND, E_ALREADY_OPEN if pool cannot be
   * deleted
   */
  virtual status_t delete_pool(const std::string& name) = 0;

  /**
   * Get a list of pool names
   *
   * @param inout_pool_names [inout] List of pool names
   *
   * @return S_OK or E_NOT_IMPL;
   */
  virtual status_t get_pool_names(std::list<std::string>& inout_pool_names) = 0;

  /**
   * Get mapped memory regions for pool.  This is used for pre-registration with
   * DMA engines.
   *
   * @param pool Pool handle
   * @param out_regions Backing file name (if any), Mapped memory regions
   *
   * @return S_OK on success or E_POOL_NOT_FOUND.  Components that do not
   * support this return E_NOT_SUPPORTED.
   */
  virtual status_t get_pool_regions(const pool_t pool, nupm::region_descriptor & out_regions)
  {
    return error_value(E_NOT_SUPPORTED, pool, out_regions); /* not supported in FileStore */
  }

  /**
   * Dynamically expand a pool.  Typically, this will add to the regions
   * belonging to a pool.
   *
   * @param pool Pool handle
   * @param increment_size Size in bytes to expand by
   * @param reconfigured_size [out] new size of pool
   *
   * @return S_OK on success or E_POOL_NOT_FOUND. Components that do not support
   * this return E_NOT_SUPPORTED (e.g. MCAS client)
   */
  virtual status_t grow_pool(const pool_t pool,
                             const size_t increment_size,
                             size_t& reconfigured_size)
  {
    PERR("grow_pool: not supported");
    return error_value(E_NOT_SUPPORTED, pool, increment_size, reconfigured_size);
  }

  /**
   * Write or overwrite an object value. If there already exists an
   * object with matching key, then it should be replaced
   * (i.e. reallocated) or overwritten.
   *
   * @param pool Pool handle
   * @param key Object key
   * @param value Value data
   * @param value_len Size of value in bytes
   *
   * @return S_OK or E_POOL_NOT_FOUND, E_KEY_EXISTS
   */
  virtual status_t put(const pool_t       pool,
                       const std::string& key,
                       const void*        value,
                       const size_t       value_len,
                       flags_t            flags = FLAGS_NONE)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, value, value_len, flags);
  }

  /**
   * Zero-copy put operation.  If there does not exist an object
   * with matching key, then an error E_KEY_EXISTS should be returned.
   *
   * @param pool Pool handle
   * @param key Object key
   * @param value Value
   * @param value_len Value length in bytes
   * @param handle Memory registration handle
   * @param flags Optional flags
   *
   * @return S_OK or E_POOL_NOT_FOUND, E_KEY_EXISTS
   */
  virtual status_t put_direct(const pool_t   pool,
                              const std::string&    key,
                              gsl::span<const common::const_byte_span> values,
                              gsl::span<const memory_handle_t> handles = gsl::span<const memory_handle_t>(),
                              flags_t        flags  = FLAGS_NONE)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, values, handles, flags);
  }

 /**
   * Zero-copy put operation (see exceptions above).
   *
   * @param pool Pool handle
   * @param key Object key
   * @param value Value
   * @param value_len Value length in bytes
   * @param handle Memory registration handle
   * @param flags Optional flags
   *
   * @return S_OK or error code
   */
  virtual status_t put_direct(const pool_t   pool,
                              const std::string&    key,
                              const void*           value,
                              const size_t          value_len,
                              const memory_handle_t handle = MEMORY_HANDLE_NONE,
                              flags_t        flags  = FLAGS_NONE)
  {
    return put_direct(
      pool
      , key
      , std::array<const common::const_byte_span,1>{common::make_const_byte_span(value, value_len)}
      , std::array<memory_handle_t,1>{handle}
      , flags
    );
  }

  /**
   * Resize memory for a value
   *
   * @param pool Pool handle
   * @param key Object key (should be unlocked)
   * @param new_size New size of value in bytes (can be more or less)
   *
   * @return S_OK on success, E_BAD_ALIGNMENT, E_POOL_NOT_FOUND,
   * E_KEY_NOT_FOUND, E_TOO_LARGE, E_ALREADY(?)
   */
  virtual status_t resize_value(const pool_t       pool,
                                const std::string& key,
                                const size_t       new_size,
                                const size_t       alignment)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, new_size, alignment);
  }

  /**
   * Read an object value
   *
   * @param pool Pool handle
   * @param key Object key
   * @param out_value Value data (if null, component will allocate memory)
   * @param out_value_len Size of value in bytes
   *
   * @return S_OK or E_POOL_NOT_FOUND, E_KEY_NOT_FOUND if key not found
   */
  virtual status_t get(const pool_t       pool,
                       const std::string& key,
                       void*&             out_value, /* release with free_memory() API */
                       size_t&            out_value_len) = 0;

  /**
   * Read an object value directly into client-provided memory.
   *
   * @param pool Pool handle
   * @param key Object key
   * @param out_value Client provided buffer for value
   * @param out_value_len [in] size of value memory in bytes [out] size of value
   * @param handle Memory registration handle
   *
   * @return S_OK, S_MORE if only a portion of value is read,
   * E_BAD_ALIGNMENT on invalid alignment, E_POOL_NOT_FOUND, E_KEY_NOT_FOUND, or other
   * error code
   *
   * Note: S_MORE is reduncant, it could have been inferred from S_OK and
   * out_value_len [in] < out_value_len [out].
   */
  virtual status_t get_direct(pool_t             pool,
                              const std::string& key,
                              void*              out_value,
                              size_t&            out_value_len,
                              memory_handle_t    handle = HANDLE_NONE)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, out_value, out_value_len, handle);
  }

  /**
   * Get attribute for key or pool (see enum Attribute)
   *
   * @param pool Pool handle
   * @param attr Attribute to retrieve
   * @param out_value Vector of attribute values
   * @param key [optional] Key
   *
   * @return S_OK on success, E_POOL_NOT_FOUND, E_INVALID_ARG, E_KEY_NOT_FOUND
   */
  virtual status_t get_attribute(pool_t                 pool,
                                 Attribute              attr,
                                 std::vector<uint64_t>& out_value,
                                 const std::string*     key = nullptr) = 0;

  /**
   * Atomically (crash-consistent for pmem) swap keys (K,V)(K',V') -->
   * (K,V')(K',V).  Before calling this API, both KV-pairs must be
   * unlocked.
   *
   * @param pool Pool handle
   * @param key0 First key
   * @param key1 Second key
   *
   * @return S_OK on success, E_POOL_NOT_FOUND, E_KEY_NOT_FOUND, E_LOCKED
   * (unable to take both locks)
   */
  virtual status_t swap_keys(const pool_t pool,
                             const std::string key0,
                             const std::string key1)
  {
    return error_value(E_NOT_SUPPORTED, pool, key0, key1);
  }

  /**
   * Set attribute on a pool.
   *
   * @param pool Pool handle
   * @param attr Attribute to set
   * @param value Vector of values to set (for boolean 0=false, 1=true)
   * @param key [optional] key
   *
   * @return S_OK, E_INVALID_ARG (e.g. key==nullptr), E_POOL_NOT_FOUND
   */
  virtual status_t set_attribute(const pool_t                 pool,
                                 const Attribute              attr,
                                 const std::vector<uint64_t>& value,
                                 const std::string*           key = nullptr)
  {
    return error_value(E_NOT_SUPPORTED, pool, attr, value, key);
  }

  /**
   * Allocate memory for zero copy DMA
   *
   * @param vaddr [out] allocated memory buffer
   * @param len [in] length of memory buffer in bytes
   * @param handle memory handle
   *
   */
  virtual status_t allocate_direct_memory(void*& vaddr,
                                          size_t len,
                                          memory_handle_t& handle)
  {
    return error_value(E_NOT_SUPPORTED, vaddr, len, handle);
  }

  /**
   * Free memory for zero copy DMA
   *
   * @param handle handle to memory region to free
   *
   * @return S_OK on success
   */
  virtual status_t free_direct_memory(memory_handle_t handle)
  {
    return error_value(E_NOT_SUPPORTED, handle);
  }

  /**
   * Register memory for zero copy DMA
   *
   * @param vaddr Appropriately aligned memory buffer
   *
   * @return Memory handle or NULL on not supported.
   */

  memory_handle_t register_direct_memory(common::const_byte_span bytes) override
  {
    return error_value(nullptr, bytes);
  }

  using Registrar_memory_direct::register_direct_memory;

  /**
   * Direct memory regions should be unregistered before the memory is released
   * on the client side.
   *
   * @param handle (as returned by register_direct_memory) of region to deregister.
   *
   * @return S_OK on success
   */
  status_t unregister_direct_memory(memory_handle_t handle) override
  {
    return error_value(E_NOT_SUPPORTED, handle);
  }

  /**
   * Take a lock on an object. If the object does not exist and inout_value_len
   * is non-zero, create it with value space according to out_value_len (this
   * is very important for mcas context). If the object does not exist and
   * inout_value_len is zero, return E_KEY_NOT_FOUND.
   *
   * @param pool Pool handle
   * @param key Key
   * @param type STORE_LOCK_READ | STORE_LOCK_WRITE
   * @param out_value [out] Pointer to data
   * @param inout_value_len [in-out] Size of data in bytes
   * @param alignment [in] Alignment of new value space in bytes.
   * @param out_key [out]  Handle to key for unlock
   * @param out_key_ptr [out]  Optional request for key-string pointer (set to
   * nullptr if not required)
   *
   * @return S_OK, S_CREATED_OK (if created on demand), E_KEY_NOT_FOUND,
   * E_LOCKED (already locked), E_INVAL (e.g., no key & no length),
   * E_TOO_LARGE (cannot allocate space for lock), E_NOT_SUPPORTED
   * if unable to take lock or other error
   */
  virtual status_t lock(const pool_t       pool,
                        const std::string& key,
                        const lock_type_t  type,
                        void*&             out_value,
                        size_t&            inout_value_len,
                        size_t             alignment,
                        key_t&             out_key_handle,
                        const char**       out_key_ptr = nullptr)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, type, out_value, inout_value_len, alignment,
                       out_key_handle, out_key_ptr);
  }

  /**
   * Unlock a key-value pair
   *
   * @param pool Pool handle
   * @param key_handle Handle (opaque) for key used to unlock
   * @param flags Optional unlock flags, UNLOCK_FLAGS_FLUSH
   *
   * @return S_OK, S_MORE (for async), E_INVAL or other error
   */
  virtual status_t unlock(const pool_t pool,
                          const key_t key_handle,
                          const unlock_flags_t flags = UNLOCK_FLAGS_NONE)
  {
    return error_value(E_NOT_SUPPORTED, pool, key_handle, flags);
  }

  /**
   * Update an existing value by applying a series of operations.
   * Together the set of operations make up an atomic transaction.
   * If the operation requires a result the operation type may provide
   * a method to accept the result. No operation currently requires
   * a result, but compare and swap probably would.
   *
   * @param pool Pool handle
   * @param key Object key
   * @param op_vector Operation vector
   * @param take_lock Set to true for automatic locking of object
   *
   * @return S_OK or error code
   */
  virtual status_t atomic_update(const pool_t                   pool,
                                 const std::string&             key,
                                 const std::vector<Operation*>& op_vector,
                                 bool                           take_lock = true)
  {
    return error_value(E_NOT_SUPPORTED, pool, key, op_vector, take_lock);
  }

  /**
   * Erase an object
   *
   * @param pool Pool handle
   * @param key Object key
   *
   * @return S_OK or error code (e.g. E_LOCKED)
   */
  virtual status_t erase(pool_t pool, const std::string& key) = 0;

  /**
   * Return number of objects in the pool
   *
   * @param pool Pool handle
   *
   * @return Number of objects
   */
  virtual size_t count(pool_t pool) = 0;

  /**
   * Apply functor to all objects in the pool
   *
   * @param pool Pool handle
   * @param function Functor to apply
   *
   * @return S_OK, E_POOL_NOT_FOUND
   */
  virtual status_t map(const pool_t pool,
                       std::function<int(const void* key,
                                         const size_t key_len,
                                         const void* value,
                                         const size_t value_len)> function)
  {
    return error_value(E_NOT_SUPPORTED, pool, function);
  }

  /**
   * Apply functor to all objects in the pool according
   * to given time constraints
   *
   * @param pool Pool handle
   * @param function Functor to apply (not in time order). If
   *                 functor returns < 0, then map aborts
   * @param t_begin Time must be after or equal. If set to zero, no constraint.
   * @param t_end Time must be before or equal. If set to zero, no constraint.
   *
   * @return S_OK, E_POOL_NOT_FOUND
   */
  virtual status_t map(const pool_t pool,
                       std::function<int(const void*              key,
                                         const size_t             key_len,
                                         const void*              value,
                                         const size_t             value_len,
                                         const common::tsc_time_t timestamp)> function,
                       const common::epoch_time_t t_begin,
                       const common::epoch_time_t t_end)
  {
    return error_value(E_NOT_SUPPORTED, pool, function, t_begin, t_end);
  }

  /**
   * Apply functor to all keys only. Useful for file_store (now deprecated)
   *
   * @param pool Pool handle
   * @param function Functor
   *
   * @return S_OK, E_POOL_NOT_FOUND
   */
  virtual status_t map_keys(const pool_t pool, std::function<int(const std::string& key)> function)
  {
    return error_value(E_NOT_SUPPORTED, pool, function);
  }

  /*
     auto iter = open_pool_iterator(pool);

     while(deref_pool_iterator(iter, ref, true) == S_OK)
       process_record(ref);

     close_pool_iterator(iter);
  */

  struct pool_reference_t {
  public:
    pool_reference_t()
      : key(nullptr), key_len(0), value(nullptr), value_len(0), timestamp() {}

    const void*          key;
    size_t               key_len;
    const void*          value;
    size_t               value_len;
    common::epoch_time_t timestamp; /* zero if not supported */

    inline std::string get_key() const {
      std::string k(static_cast<const char*>(key), key_len);
      return k;
    }

  };

  /**
   * Open pool iterator to iterate over objects in pool.
   *
   * @param pool Pool handle
   *
   * @return Pool iterator or nullptr
   */
  virtual pool_iterator_t open_pool_iterator(const pool_t pool)
  {
    return error_value(nullptr, pool);
  }

  /**
   * Deference pool iterator position and optionally increment
   *
   * @param pool Pool handle
   * @param iter Pool iterator
   * @param t_begin Time must be after or equal. If set to zero, no constraint.
   * @param t_end Time must be before or equal. If set to zero, no constraint.
   * @param ref [out] Output reference record
   * @param ref [out] Set to true if within time bounds
   * @param increment Move iterator forward one position
   *
   * @return S_OK on success and valid reference, E_INVAL (bad iterator),
   *   E_OUT_OF_BOUNDS (when attempting to dereference out of bounds)
   *   E_ITERATOR_DISTURBED (when writes have been made since last iteration)
   */
  virtual status_t deref_pool_iterator(const pool_t       pool,
                                       pool_iterator_t    iter,
                                       const common::epoch_time_t t_begin,
                                       const common::epoch_time_t t_end,
                                       pool_reference_t&  ref,
                                       bool&              time_match,
                                       bool               increment = true)
  {
    /* Not sure of the difference between "not supported" and "not implemented",
     * but they are separated codes.
     */
    return error_value(E_NOT_IMPL, pool, iter, t_begin, t_end, ref, time_match, increment);
  }

  /**
   * Unlock pool, release iterator and associated resources
   *
   * @param pool Pool handle
   * @param iter Pool iterator
   *
   * @return S_OK on success, E_INVAL (bad iterator)
   */
  virtual status_t close_pool_iterator(const pool_t pool,
                                       pool_iterator_t iter)
  {
    return error_value(E_NOT_IMPL, pool, iter);
  }

  /**
   * Free server-side allocated memory
   *
   * @param p Pointer to memory allocated through a get call
   *
   * @return S_OK on success
   */
  virtual status_t free_memory(void* p)
  {
    ::free(p);
    return S_OK;
  }

  /**
   * Allocate memory from pool
   *
   * @param pool Pool handle
   * @param size Size in bytes
   * @param alignment Alignment hint in bytes, 0 if no alignment is needed
   * @param out_addr Pointer to allocated region
   *
   * @return S_OK on success, E_BAD_ALIGNMENT, E_POOL_NOT_FOUND, E_NOT_SUPPORTED
   */
  virtual status_t allocate_pool_memory(const pool_t pool,
                                        const size_t size,
                                        const size_t alignment_hint,
                                        void*&       out_addr)
  {
    return error_value(E_NOT_SUPPORTED, pool, size, alignment_hint, out_addr);
  }

  /**
   * Free memory from pool
   *
   * @param pool Pool handle
   * @param addr Address of memory to free
   * @param size Size in bytes of allocation; if provided this accelerates
   * release
   *
   * @return S_OK on success, E_INVAL, E_POOL_NOT_FOUND, E_NOT_SUPPORTED
   */
  virtual status_t free_pool_memory(pool_t pool, const void* addr, size_t size = 0)
  {
    return error_value(E_NOT_SUPPORTED, pool, addr, size);
  }

  /**
   * Flush memory from pool
   *
   * @param pool Pool handle
   * @param addr Address of memory to flush
   * @param size Size in bytes to flush
   *
   * @return S_OK on success, E_INVAL, E_POOL_NOT_FOUND, E_NOT_SUPPORTED
   */
  virtual status_t flush_pool_memory(pool_t pool, const void* addr, size_t size)
  {
    return error_value(E_NOT_SUPPORTED, pool, addr, size);
  }

  /**
   * Perform control invocation on component
   *
   * @param command String representation of command (component-interpreted)
   *
   * @return S_OK on success or error otherwise
   */
  virtual status_t ioctl(const std::string& command) { return error_value(E_NOT_SUPPORTED, command); }

  /**
   * Debug routine
   *
   * @param pool Pool handle
   * @param cmd Debug command
   * @param arg Parameter for debug operation
   */
  virtual void debug(pool_t pool, unsigned cmd, uint64_t arg) = 0;
};

}  // namespace component
#endif

max_stars_repo_name: IBM/artemis

max_issues_repo_name: IBM/artemis

max_forks_repo_name: IBM/artemis


text:
#ifndef CHEVALIER_H
#define CHEVALIER_H

#include <math.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_multiroots.h>

#define parsec_in_cm 3.085678e+18 //parsec in cm
#define parsec_in_km 3.085678e+13 //parsec in km
#define msun_in_g    1.988920e+33 //msun in g
#define year_in_sec  3.155760e+07 //year in s

typedef struct Chevalier
{
    double M_dot;               //mass input rate
    double E_dot;               //energy input rate
    double R;                   //limit of mass and energy input
    double gamma;               //adiabatic index
    double V;                   //Volume of input area
    double q;                   //M_dot / V
    double Q;                   //E_dot / V
    gsl_root_fsolver *grsolve;  //GSL root solver
} Chevalier;

// Function declarations for Chevalier struct
void SetChevalier(Chevalier *c, double M_dot, double E_dot, double gamma, double R);
double MachNumber(Chevalier *c, double r);      //mach number vs. radius
double MomentumDensity(Chevalier *c, double r); //momentum density vs. radius in Msun/pc^3 * km/s
double Pressure(Chevalier *c, double r);        //Pressure in dyn cm^-2
double Density(Chevalier *c, double r);         //density in g cm^-3
double Density_Msunpc3(Chevalier *c, double r); //density in Msun/pc3
double WindVelocity(Chevalier *c, double r);    //wind velocity in km/s
double EnergyIntegral(Chevalier *c, double r);  //msun/yr (km/s)^2 / pc^3
double u_star(Chevalier *c, double r);          //normalized wind velocity
double rho_star(Chevalier *c, double r);        //normalized density
double P_star(Chevalier *c, double r);          //normalized pressure

double mach_crossing_A(double M, void *fp);     //functions for Mach number root finding
double mach_crossing_B(double M, void *fp);

// Function definitions
void SetChevalier(Chevalier *c, double M_dot, double E_dot, double gamma, double R) {
    // M_dot = Mass input in Msun/yr
    // E_dot = Energy input in erg / s
    // gamma = Adiabatic index
    // R     = Radius of the wind region in parsec

    //Set model parameters
    c->M_dot = M_dot;
    c->E_dot = E_dot;
    c->gamma = gamma;
    c->R     = R;
    c->V     = (4.*M_PI/3. * R*R*R); //volume in pc^3
    //compute input density rate (msun / yr / pc^3)
    c->q     = c->M_dot / c->V;
    //compute input energy density rate (msun/yr (km/s)^2 / pc^3)
    c->Q     = (c->E_dot/msun_in_g*year_in_sec / 1.0e10) / c->V;
    //define the GSL root solver to use
    const gsl_root_fsolver_type *grsolve_T = gsl_root_fsolver_brent;
    c->grsolve = gsl_root_fsolver_alloc(grsolve_T);
}

double u_star(Chevalier *c, double r) {
    // Returns dimensionless wind velocity
    double u = WindVelocity(c, r);
    double u_prime = sqrt(c->E_dot)/sqrt(c->M_dot*msun_in_g/year_in_sec)/1.0e5;
    return u/u_prime;
}

double rho_star(Chevalier *c, double r) {
    // Returns dimensionless density
    double rho = Density(c,r); //g cm^-3
    double M_dot_prime = c->M_dot * msun_in_g/year_in_sec; // g/s
    double R_prime = c->R * parsec_in_cm; //cm
    double rho_prime = pow(M_dot_prime,1.5)/(sqrt(c->E_dot)*R_prime*R_prime);
    return rho/rho_prime;
}

double P_star(Chevalier *c, double r) {
    // Returns dimensionless pressure
    double P = Pressure(c,r); //dyn cm^-2
    double M_dot_prime = c->M_dot * msun_in_g/year_in_sec; // g/s
    double R_prime = c->R * parsec_in_cm; //cm
    double P_prime = sqrt(M_dot_prime)*sqrt(c->E_dot)/(R_prime*R_prime);
    return P/P_prime;
}

double Pressure(Chevalier *c, double r) {
    // Returns pressure in dyn cm^-2
    double Mach = MachNumber(c,r);          //mach number
    double u    = WindVelocity(c,r);        //km/s
    double cs   = u/Mach * 1.0e5;           //sound speed in cm/s
    u *= year_in_sec / parsec_in_km;        //to pc/yr
    double rho  = MomentumDensity(c,r)/u;   //Msun / pc^3
    rho *= msun_in_g;                       //g / pc^3
    rho /= pow(parsec_in_cm, 3);            //g / cm^3
    double P = cs*cs*rho/c->gamma;          //pressure in cgs
    return P;                               //in dyn cm^-2
}

double Density(Chevalier *c, double r) {
    // Returns density in g cm^-3
    double rho = Density_Msunpc3(c, r);     //Msun/pc^3
    rho *= msun_in_g;                       //g / pc^3
    rho /= pow(parsec_in_cm, 3);            //g / cm^3
    return rho;                             // in g/cm^3
}

double Density_Msunpc3(Chevalier *c, double r) {
    // Returns density in g cm^-3
    double u = WindVelocity(c, r);          //km/s
    u *= year_in_sec / parsec_in_km;        //to pc/yr
    return MomentumDensity(c,r)/u;          //Msun/pc^3;
}

double WindVelocity(Chevalier *c, double r) {
    // Returns wind velocity in km/s

    //First, find integrated energy input density
    double Qint = EnergyIntegral(c, r);     //msun/yr (km/s)^2 / pc^3
    //Second, find momentum density
    double rhou = MomentumDensity(c, r);     //msun/yr/pc^2
    //find sq of velocity (modulo gamma + Mach correction)
    double usq = Qint/rhou;                 //1/2 u^2 + (gamma/(gamma-1))*P/rho
                                            //1/2 u^2 + c^2 / (gamma-1)
                                            //1/2 u^2 + u^2 / (M^2 (gamma-1))
                                            //u^2 * ( 1/2 + 1/(M^2 (gamma-1)) )
    //get the mach number
    double Mach = MachNumber(c,r);
    //find the adjustment factor
    double fac = (0.5 + 1./(Mach*Mach*(c->gamma-1.)));
    return sqrt(usq/fac);	//km/s
}

double EnergyIntegral(Chevalier *c, double r) {
    double Qint = r<c->R ? 1./3.*c->Q*r : 1./3.*c->Q*c->R*c->R*c->R/(r*r);
    return Qint;    //msun/yr (km/s)^2 / pc^3
}

double MomentumDensity(Chevalier *c, double r) {
    double rhou = r<c->R ? 1./3.*c->q*r : 1./3.*c->q*c->R*c->R*c->R/(r*r) ;
    return rhou;    // in Msun/yr/pc^2
}

double MachNumber(Chevalier *c, double r) {

    int    status, iter=0, max_iter=100;
    double M_lo = 1.0e-5, M_hi = 5.;
    double Mx, answer = 0;
    double x = r/c->R;
    gsl_function func;
    double fp[2] = {c->gamma, x};

    //choose which solution to use
    if(x<=1.0) {
        M_lo = 1.0e-5;
        M_hi = 1.0;
        func.function = &mach_crossing_A;
    }
    else {
        M_lo = 1.0;
        M_hi = 10000.0;
        func.function = &mach_crossing_B;
    }
    func.params   = &fp[0];

    gsl_root_fsolver_set(c->grsolve, &func, M_lo, M_hi);

    while(++iter<max_iter){
        status = gsl_root_fsolver_iterate(c->grsolve);
        Mx = gsl_root_fsolver_root(c->grsolve);
        M_lo = gsl_root_fsolver_x_lower(c->grsolve);
        M_hi = gsl_root_fsolver_x_upper(c->grsolve);
        status = gsl_root_test_interval(M_lo,M_hi,0,1.0e-5);
        if(status==GSL_SUCCESS) {answer = Mx; break;}
    }

    return answer;
}

double mach_crossing_A(double M, void *fp) {
    double *g     = (double *) fp;
    double gamma  = g[0];
    double x      = g[1];
    double alpha  = -1.*(3.*gamma+1.)/(5.*gamma+1.);
    double beta   = (gamma+1.)/(2.*(5.*gamma+1.));
    double A      = pow( (3.*gamma + 1./(M*M))/(1.+3.*gamma), alpha);
    double B      = pow( (gamma-1.+2./(M*M))/(1.+gamma), beta);
    return A*B - x;
}

double mach_crossing_B(double M, void *fp) {
    double *g     = (double *) fp;
    double gamma = g[0];
    double x     = g[1];
    double alpha = 2./(gamma-1.);
    double beta  = (gamma+1.)/(2.*(gamma-1.));
    double A     = pow( M, alpha);
    double B     = pow( (gamma-1.+2./(M*M))/(1.+gamma), beta);
    return A*B - x*x;
}

#endif //CHEVALIER_H

ext: h

max_stars_repo_name: editeodoro/galpy

max_issues_repo_name: editeodoro/galpy

max_forks_repo_name: editeodoro/galpy


text:
/*
 * SPDX-License-Identifier:  BSD-3-Clause
 *
 * example_02-StructOfArrays-Naive-Omp.c :
 *      Example of SPH Density Calculation using a
 *      naive implementation of the main density loop,
 *      no neighbours earch, and Struct of Arrays (SoA)
 *      data layout and OpenMP parallelization.
 *
 * (C) Copyright 2021 José Hugo Elsas
 * Author: José Hugo Elsas <jhelsas@gmail.com>
 *
 * Command Line Options:
 *   -runs  <int>   : Set the number of repetitions (runs) for
 *                      calculating the density. The value of
 *                      the density is based on the last
 *                      iteration.
 *                    Default value: 1
 *   -run_seed <int>: Flag to set an alternative seed use for
 *                      for the PRNG. Instead of feeding seed
 *                      to the PRNG directly, it feeds
 *                      seed + iteration, as to generate different
 *                      configurations for each iteration.
 *                    Default value: 0 - (possible 0/1)
 *   -seed     <int>: Set the seed to use for the SPH particles
 *                      uniform position generation in the box
 *                    Default value: 123123123
 *
 *   -N        <int>: Set the number of SPH particles to be used
 *                    Default value: 1e5 = 100,000
 *   -h      <float>: Set the value of the smoothing kernel
 *                      parameter h, which corresponds to half
 *                      of the support of the kernel.
 *                    Default value: 0.05
 *
 *   -Nx       <int>: Set the number of Cells in the X direction
 *                    Default value: 10
 *   -Ny       <int>: Set the number of Cells in the Y direction
 *                    Default value: 10
 *   -Nz       <int>: Set the number of Cells in the Z direction
 *                    Default value: 10
 *
 *   -Xmin   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Y direction for
 *                    the Cell Linked List box
 *                      Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *
 *   -Xmax   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Ymax   <float>: Set the lower bound in the Y direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Zmax   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 */

#include <math.h>
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <stdbool.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_utils.h"

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#define COMPUTE_BLOCKS 1

int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr, linkedListBox *box, SPHparticle *lsph, double *times);

int compute_density_3d_naive_omp(int N,double h,
                                 double* restrict x, double* restrict y,
                                 double* restrict z, double* restrict nu,
                                 double* restrict rho);

double w_bspline_3d(double r,double h);

int main(int argc, char **argv){
  bool run_seed = false;       // By default the behavior is is to use the same seed
  int runs = 1,err;            // it only runs once
  long int seed = 123123123;   // The default seed is 123123123
  int64_t N = 100000;          // The default number of particles is N = 1e5 = 100,000
  double h=0.05;               // The default kernel smoothing length is h = 0.05
  linkedListBox *box;          // Uninitialized Box containing the cells for the cell linked list method
  SPHparticle *lsph;           // Uninitialized array of SPH particles

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox)); // Create a box

  // allow for command line customization of the run
  arg_parse(argc,argv,&N,&h,&seed,&runs,&run_seed,box);  // Parse the command
                                                         // line arguments and override default values

  err = SPHparticle_SoA_malloc(N,&lsph);                 // Create an array of N particles
  if(err)
    fprintf(stderr,"error in SPHparticle_SoA_malloc\n");

  void *swap_arr = malloc(N*sizeof(double));
  double times[runs*COMPUTE_BLOCKS];

  for(int run=0;run<runs;run+=1)
    main_loop(run,run_seed,N,h,seed,swap_arr,box,lsph,times);

  bool is_cll = false;
  const char *prefix = "ex02,naive,SoA,omp";
  print_time_stats(prefix,is_cll,N,h,seed,runs,lsph,box,times);
  print_sph_particles_density(prefix,is_cll,N,h,seed,runs,lsph,box);

  SPHparticleSOA_safe_free(N,&lsph);
  safe_free_box(box);
  free(swap_arr);

  return 0;
}

/*
 *  Function main_loop:
 *    Runs the main loop of the program, including the particle array generation,
 *    density calculation and the timings annotations.
 *
 *    Arguments:
 *       run <int>            : index (or value) or the present iteration
 *       run_seed <bool>      : boolean defining whether to use run index for seed or not
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       seed <long int>      : seed for GSL PRNG generator to generate particle positions
 *       box  <linkedListBox> : Box of linked list cells, encapsulating the 3d domain
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *       times <double>       : Array to store the computation timings to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 *       times <double>       : Times is updated by reference
 */
int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr, linkedListBox *box, SPHparticle *lsph, double *times)
{
  int err;

  // Initialize the particles' positions and other values
  if(run_seed)
    err = gen_unif_rdn_pos_box(N,seed+run,box,lsph);
  else
    err = gen_unif_rdn_pos_box(N,seed,box,lsph);

  if(err)
    fprintf(stderr,"error in gen_unif_rdn_pos\n");

  // ------------------------------------------------------ //

  double t0,t1;

  t0 = omp_get_wtime();

  compute_density_3d_naive_omp(N,h,lsph->x,lsph->y,lsph->z,
                               lsph->nu,lsph->rho);         // Compute the density for all particles

  t1 = omp_get_wtime();

  // ------------------------------------------------------ //

  times[COMPUTE_BLOCKS*run+0] = t1-t0;                      // Only one component to measure time

  return 0;
}

/*
 *  Function compute_density_3d_naive_omp:
 *    Computes the SPH density from the particles naively (i.e. direct loop).
 *    The outer-most loop is parallelized using openMP.
 *
 *    Arguments:
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       x       <double*>    : Array of particles' X positions
 *       y       <double*>    : Array of particles' Y positions
 *       z       <double*>    : Array of particles' Z positions
 *       nu      <double*>    : Array of particles' density weights (i.e. masses)
 *    Returns:
 *       0                    : error code returned
 *       rho       <double*>  : Array of particles' densities
 */
int compute_density_3d_naive_omp(int N,double h,
                                 double* restrict x, double* restrict y,
                                 double* restrict z,double* restrict nu,
                                 double* restrict rho){
  memset(rho,(int)0,N*sizeof(double));        // Pre-initialize the density to zero

  #pragma omp parallel for                    // Iterate in parallel
  for(int64_t ii=0;ii<N;ii+=1){               // For every particle
    for(int64_t jj=0;jj<N;jj+=1){             // Run over every other particle
      double dist = 0.;                       // the contributions from each direction by

      dist += (x[ii]-x[jj])*(x[ii]-x[jj]);    // adding the contribution in the X direction
      dist += (y[ii]-y[jj])*(y[ii]-y[jj]);    // adding the contribution in the Y direction
      dist += (z[ii]-z[jj])*(z[ii]-z[jj]);    // adding the contribution in the Z direction

      dist = sqrt(dist);                      // take the sqrt to have the distance

      rho[ii] += nu[jj]*w_bspline_3d(dist,h); // and add the contribution to the density
    }
  }

  return 0;
}

/*
 *  Function w_bspline_3d:
 *    Returns the normalized value of the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       q <double>           : Distance between particles
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *    Returns:
 *       wq <double>          : Normalized value of the kernel
 */
double w_bspline_3d(double r,double h){
  const double A_d = 3./(2.*M_PI*h*h*h);      // The 3d normalization constant
  double q=0.;                                // normalized distance, initialized to zero

  if(r<0||h<=0.)                              // If either distance or smoothing length
    exit(10);                                 // are negative, declare an emergency

  q = r/h;                                    // Compute the normalized distance
  if(q<=1)                                    // If the distance is small
    return A_d*(2./3.-q*q + q*q*q/2.0);       // Compute this first polynomal
  else if((1.<=q)&&(q<2.))                    // If the distance is a bit larger
    return A_d*(1./6.)*(2.-q)*(2.-q)*(2.-q);  // Compute this other polynomial
  else                                        // Otherwise, if the distance is large
    return 0.;                                // The value of the kernel is 0
}

ext: c

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
/**
 *
 * @file core_dgeqp3_tntpiv.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:50 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************
 *
 * @ingroup CORE_double
 *
 *  CORE_dgeqp3_tntpiv computes a QR factorization with column pivoting of a
 *  matrix A:  A*P = Q*R  using Level 3 BLAS.
 *
 *  The matrix Q is represented as a product of elementary reflectors
 *
 *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
 *
 *  Each H(i) has the form
 *
 *     H(i) = I - tau * v * v**T
 *
 *  where tau is a complex scalar, and v is a real/complex vector
 *  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
 *  A(i+1:m,i), and tau in TAU(i).
 *
 *******************************************************************************
 *
 *  Arguments:
 *  ==========
 *
 * @param[in] m
 *          The number of rows of the matrix A. M >= 0.
 *
 * @param[in] n
 *          The number of columns of the matrix A.  N >= 0.
 *
 * @param[in,out] A
 *          A is COMPLEX*16 array, dimension (LDA,N)
 *          On entry, the M-by-N matrix A.
 *          On exit, the upper triangle of the array contains the
 *          min(M,N)-by-N upper trapezoidal matrix R; the elements below
 *          the diagonal, together with the array TAU, represent the
 *          unitary matrix Q as a product of min(M,N) elementary
 *          reflectors.
 *
 * @param[in] lda
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[out] IPIV
 *          IPIV is INTEGER array, dimension min(M,N)
 *          The pivot indices; for 1 <= j <= min(M,N), column j of the
 *          tile was interchanged with column IPIV(j).
 *
 * @param[out] TAU
 *          TAU is COMPLEX*16 array, dimension (min(M,N))
 *          The scalar factors of the elementary reflectors.
 *
 * @param[in,out] iwork
 *          iwork is INTEGER array, dimension (N)
 *          On entry, if iwork(J).ne.0, the J-th column of A is permuted
 *          to the front of A*P (a leading column); if iwork(J)=0,
 *          the J-th column of A is a free column.
 *          On exit, if iwork(J)=K, then the J-th column of A*P was the
 *          the K-th column of A.
 *
 * @param[out] INFO
 *          = 0: successful exit.
 *          < 0: if INFO = -i, the i-th argument had an illegal value.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dgeqp3_tntpiv = PCORE_dgeqp3_tntpiv
#define CORE_dgeqp3_tntpiv PCORE_dgeqp3_tntpiv
#endif
int CORE_dgeqp3_tntpiv(int m, int n,
                       double *A, int lda,
                       int *IPIV, double *tau,
                       int *iwork)
{
    int i, tmp, info;
    memset(iwork, 0, n*sizeof(int));
    info = LAPACKE_dgeqp3(LAPACK_COL_MAJOR, m, n, A, lda, iwork, tau );

    /* Convert IPIV from permutation array, to pivot array
     * WARNING: this is because this kernel is only used in
     * tournament pivoting with rank revealing QR */
    if (info == 0) {
        for(i=0; i<min(m,n); i++) {
            assert(iwork[i] != 0 );

            tmp = iwork[i]-1;
            while( tmp < i ) {
                tmp = IPIV[ tmp ] - 1;
            }
            IPIV[i] = tmp+1;
        }
    }
    return info;
}


text:
#include <cblas.h>
#include <stdio.h>

void main()
{
  int i=0;
  double A[6] = {1.0,2.0,1.0,-3.0,4.0,-1.0};
  double B[6] = {1.0,2.0,1.0,-3.0,4.0,-1.0};
  double C[9] = {.5,.5,.5,.5,.5,.5,.5,.5,.5};
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,3,3,2,1,A, 3, B, 3,2,C,3);

  for(i=0; i<9; i++)
    printf("%lf ", C[i]);
  printf("\n");
}

max_stars_repo_name: chelini/openBlas

max_issues_repo_name: chelini/openBlas

max_forks_repo_name: chelini/openBlas