text:
/*
 * This file is part of the ProVANT simulator project.
 * Licensed under the terms of the MIT open source license. More details at
 * https://github.com/Guiraffo/ProVANT-Simulator/blob/master/LICENSE.md
 */
/**
 * @file sdf_parser.h
 * @brief This file contains the declaration of the SDFParser class.
 *
 * @author Júnio Eduardo de Morais Aquino
 */

#ifndef PROVANT_SDF_READER_H
#define PROVANT_SDF_READER_H

#include <sdf/sdf.hh>

#include <gsl/gsl>

#include "sdf_status.h"

/**
 * @brief The SDFParser class is a helper class to read properties from a SDF
 * Element object.
 *
 * This class is specially usefull during the loading of Gazebo plugins where
 * is necessary to read parameters from the plugin SDF description.
 *
 * A SDF (Simulation Description Format) Element, is a specialized XML element
 * that may contain attributes and child elements, as any other XML element
 * can.
 *
 * Gazebo provides a SDF Reader, but many times it has a confusing syntax and
 * very lax documentation, which is very hard to understand and discourages the
 * usage of the SDF class to read its own elements.
 *
 * As there were many instances in the ProVANT Simulator source code that
 * users and developers were using the now depreacted XMLRead class to read the
 * values of SDF elements, this class was designed as a substitution for the
 * XMLRead class, that wraps the methods provided by gazebo into well documented
 * and tested methods.
 *
 * So far, this class provides methods to check the existence, and returns the
 * string values of attributes and elements. In the future, it is very likely
 * that this class will also contain methods to read the values of elements in
 * other types, checking for adequate conversion and validity.
 *
 * @todo It turns out that is really hard to create SDF elements in memory
 * in order to do the testing of this class, so by now THIS CLASS IS UNTESTED!
 * It is imperative that we develop tests for this class. The tests are not
 * hard to create, the problem is that it creating SDF descriptions in memory
 * throws a series of undocumented errors and leads the test executable to a
 * segmentation fault, rendering the test procedure unfeasible.
 */
class SDFParser
{
public:
  /**
   * @brief Construct a new SDFParser object and initializes the SDF pointer.
   *
   * @param sdf SDF element to read data from.
   */
  SDFParser(sdf::ElementPtr sdf);
  /**
   * @brief Destroy the SDFParser object.
   * Releases the shared pointer to the SDF element.
   */
  virtual ~SDFParser();

  /**
   * @brief Get the SDF ptr from which this class reads data from.
   * @return sdf::ElementPtr
   */
  virtual sdf::ElementPtr GetElementPtr() const;

  /**
   * @brief Checks if the SDF contain an attribute with the specified name.
   *
   * Attributes are data specific to a element, that are usually strings.
   * For example, the following XML element:
   * @code{.xml}
   * <test type="example">
   *  <other_data/>
   * </test>
   * @endcode
   *
   * Contains an attribute named type and with value example.
   *
   * Usage example:
   * Suppose a plugin SDF description as follows:
   * @code{.xml}
   * <plugin name="example" filename="provant_simulator_step_plugin">
   *  <plugin_data/>
   * </plugin>
   * @endcode
   *
   * Assuming this plugin was read, the following code has res=true, because
   * the plugin contains an attribute "name".
   * @code{.cpp}
   * SDFParser reader(_sdf);
   * bool res = reader.HasAttribute("name");
   * @endcode
   *
   * And the following code results in res=false, because the plugin does not
   * contain an attribute "type".
   * @code{.cpp}
   * SDFParser reader(_sdf);
   * bool res = reader.HasAttribute("type");
   * @endcode
   *
   * @param name
   * @return true If the SDF element contains an attribute with the specified
   * name.
   * @return false Otherwise.
   */
  virtual bool HasAttribute(const std::string& name) const;
  /**
   * @brief Checks if the SDF has an element with the specified name.
   *
   * Usage example:
   * Suppose a plugin SDF description as follows:
   * @code{.xml}
   * <plugin name="example" filename="provant_simulator_step_plugin">
   *  <plugin_data/>
   * </plugin>
   * @endcode
   *
   * Assuming this plugin was read, the following code has res=true, because
   * the plugin contains an element "plugin_data".
   * @code{.cpp}
   * SDFParser reader(_sdf);
   * bool res = reader.HasElement("plugin_data");
   * @endcode
   *
   * And the following code results in res=false, because the plugin does not
   * contain an attribute "type".
   * @code{.cpp}
   * SDFParser reader(_sdf);
   * bool res = reader.HasElement("type");
   * @endcode
   *
   * @param name Name of the element to look for.
   * @return true If the SDF contains an element with the specified name.
   * @return false Otherwise.
   */
  virtual bool HasElement(const std::string& name) const;

  /**
   * @brief Get an attribute value.
   *
   * This checks if an attribute with the specified value exists, and if it
   * exists, copies its value to the value parameter.
   *
   * If the attribute does not exists, or any other occurs, an SDFStatus with an
   * error and a message describing the problem is returned, and the value is
   * left emtpy.
   *
   * @param name Name of the attribute to return.
   * @param value Buffer string that will contain the value of the attribute with the specified name.
   * @return SDFStatus that informs the status of the operation. If the
   * attribute is found and returned, a OkStatus object is returned, in the case
   * the attribute does not exist an AttributeNotFoundError object is returned,
   * containing an error status and a message informing the problem.
   */
  virtual SDFStatus GetAttributeValue(const std::string& name, gsl::not_null<std::string*> value) const noexcept;
  /**
   * @brief Get an attribute value.
   *
   * This method checks if an attribute with the specified name exists, and if
   * it exists, returns its value.
   *
   * @param name Name of the attribute to locate.
   * @return std::string Value of the attribute.
   * @throw AttributeNotFoundError If an attribute with the specifed name does
   * not exist.
   */
  virtual std::string GetAttributeValue(const std::string& name) const;
  /**
   * @brief Get the text value of a SDF element.
   *
   * Check if an element with the specified name exists, and if it exists,
   * copies its value to the value parameter.
   *
   * If an element with the specified name does not exist, or any other
   * error occurs, the value string is left empty, and a SDFStatus error is
   * returned with a message identifying the problem.
   *
   * @param name Name of the element to read the value.
   * @param value Buffer string that will contain the value of the element
   * @return SDFStatus that informs the status of the operation. If at least one
   * element with the specified name exists, a OkStatus object is returned.
   * In case an error occurs, an object identifying the problem is returned
   * with an errors status and a message describing the problem.
   */
  virtual SDFStatus GetElementText(const std::string& name, gsl::not_null<std::string*> value) const noexcept;
  /**
   * @brief Get the text value of a SDF element.
   *
   * Check if an element with the specified name exists, and return is value.
   *
   * @sa GetElementText(const std::string& name, std::string* value)
   *
   * @param name Name of the element to read the value.
   * @return std::string Content of the element.
   * @throw ElementNotFoundError if the element cannot be found as a child of
   * the SDF element this class reads from.
   */
  virtual std::string GetElementText(const std::string& name) const;

  /**
   * @brief Get the value of a boolean SDF element.
   *
   * The valid options are:
   * For true elements: true in any case combination, or 1.
   * For false elements: false in any case combinatation, or 0.
   *
   * @sa GetElementBool(const std::string&)
   *
   * @param name Name of the element to read.
   * @param value Content of the element.
   * @return SDFStatus Status of the parse operation. If the element was found
   * and had a valid value, an OKStatus is returned, otherwise an object
   * containing and error status and error message is returned.
   */
  virtual SDFStatus GetElementBool(const std::string& name, gsl::not_null<bool*> value) const noexcept;
  /**
   * @brief Get the value of a boolean SDF element.
   *
   * The valid options are:
   * For true elements: true in any case combination, or 1.
   * For false elements: false in any case combinatation, or 0.
   *
   * @sa GetElementBool(const std::string&, gsl::not_null<int*>)
   *
   * @param name Name of the element to read.
   * @return Value of the element.
   * @throw SDFStatus exception in case the element is not found or has
   * invalid content, ie. content that cannot be parsed as a boolean.
   */
  virtual bool GetElementBool(const std::string& name) const;

  /**
   * @brief Get the value of a int SDF element.
   *
   * This method can read integer values in the binary (0b), octal (0), decimal
   * (no prefix), or hexadecimal (0x) formats, identified by their prefixes.
   * It can also parse signed values with either a - or + signal.
   *
   * @sa GetElementInt(const std::string&)
   *
   * @param name Name of the element to read.
   * @param value Value of the element.
   * @return SDFStatus Status of the parse operation. If the element was found
   * and had a valid value, an OKStatus is returned, otherwise an object
   * containing and error status and error message is returned.
   */
  virtual SDFStatus GetElementInt(const std::string& name, gsl::not_null<int*> value) const noexcept;
  /**
   * @brief Get the value of a int SDF element.
   *
   * This method can read integer values in the binary (0b), octal (0), decimal
   * (no prefix), or hexadecimal (0x) bases, identified by their prefixes.
   * It can also parse signed values with either a - or + signal.
   *
   * @sa GetElementInt(const std::string&, gsl::not_null<int*>)
   *
   * @param name Name of the element to read.
   * @return int Value of the element.
   * @throw SDFStatus exception in case the element is not found or has
   * invalid content, ie. content that cannot be parsed as an integer such as
   * string, or if the value of the element is beyond the current platform limits
   * for an integer value.
   */
  virtual int GetElementInt(const std::string& name) const;

  /**
   * @brief Get the value of an unsigned int SDF element.
   *
   * This method can read unsigned integer values in the binary (0b), octal(0),
   * decimal (no prefix) or hexadecimal (0x) bases.
   * In case this methods find a signed integer, it returns an error.
   *
   * @sa GetElementUnsignedInt(const std::string&)
   *
   * @param name Name of the element to read.
   * @param value Value of the element.
   * @return SDFStatus Status of the parse operation. If the element was found
   * and had a valid value, an OKStatus is returned, otherwise an object
   * containing and error status and error message is returned.
   */
  virtual SDFStatus GetElementUnsignedInt(const std::string& name, gsl::not_null<unsigned int*> value) const noexcept;
  /**
   * @brief Get the value of an unsigned int SDF element.
   *
   * This method can read unsigned integer values in the binary (0b), octal(0),
   * decimal (no prefix) or hexadecimal (0x) bases.
   * In case this methods find a signed integer, it throws an exception.
   *
   * @sa GetElementUnsignedInt(const std::string&, gsl::not_null<unsigned int*>)
   *
   * @param name Name of the element to read.
   * @return unsigned int Value of the element.
   * @throws SDFStatus exception in case the element is not found or has
   * invalid content, ie. content that cannot be parsed as an integer such as
   * string, if the value is signed, or if the value is beyond the limits of an
   * unsigned integer for the current platform.
   */
  virtual unsigned int GetElementUnsignedInt(const std::string& name) const;

  /**
   * @brief Get the value of a single precision (float) SDF Element.
   *
   * This method can read real values as single precision floating point values
   * (float). It is also possible to read floats specified in the hexadecimal
   * format with the 0x prefix.
   *
   * @sa GetElementFloat(const std::string&)
   *
   * @param name Name of the element to read.
   * @param value Value of the element.
   * @return SDFStatus Status of the parse operation. If the element was found
   * and had a valid value, an OKStatus is returned, otherwise an object
   * containing and error status and error message is returned.
   */
  virtual SDFStatus GetElementFloat(const std::string& name, gsl::not_null<float*> value) const noexcept;
  /**
   * @brief Get the value of a single precision (float) SDF Element.
   *
   * This method can read real values as single precision floating point values
   * (float). It is also possible to read floats specified in the hexadecimal
   * format with the 0x prefix.
   *
   * @sa GetElementFloat(const std::string&, gsl::not_null<float*>)
   *
   * @param name Name of the element to read.
   * @return float Value of the element.
   * @throws SDFStatus exception in case the element is not found or has
   * invalid content, ie. content that cannot be parsed as a float or is beyond
   * the limits of a float for the current platform.
   */
  virtual float GetElementFloat(const std::string& name) const;

  /**
   * @brief Get the value of a double precision floating point (double) SDF Element.
   *
   * This method can read real values as double precision floating point values
   * (double). It is also possible to read doubles specified in the hexadecimal
   * format with the 0x prefix.
   *
   * @sa GetElementDouble(const std::string&)
   *
   * @param name Name of the element to read.
   * @param value Value of the element.
   * @return SDFStatus Status of the parse operation. If the element was found
   * and had a valid value, an OKStatus is returned, otherwise an object
   * containing and error status and error message is returned.
   */
  virtual SDFStatus GetElementDouble(const std::string& name, gsl::not_null<double*> value) const noexcept;
  /**
   * @brief Get the value of a double precision floating point (double) SDF Element.
   *
   * This method can read real values as double precision floating point values
   * (double). It is also possible to read doubles specified in the hexadecimal
   * format with the 0x prefix.
   *
   * @sa GetElementDouble(const std::string&, gsl::not_null<double*>)
   *
   * @param name Name of the element to read.
   * @return double Value of the element.
   * @throws SDFStatus exception in case the element is not found or has
   * invalid content, ie. content that cannot be parsed as a double or is beyond
   * the limits of a double for the current platform.
   */
  virtual double GetElementDouble(const std::string& name) const;

private:
  //! Stores the pointer of the SDF from which this class reads the values of attributes and elements.
  sdf::ElementPtr _sdf;
};

#endif  // PROVANT_SDF_READER_H


text:
/*
    -- MAGMA (version 2.5.4) --
       Univ. of Tennessee, Knoxville
       Univ. of California, Berkeley
       Univ. of Colorado, Denver
       @date October 2020
*/

#ifndef MAGMA_TYPES_H
#define MAGMA_TYPES_H

#include <stdint.h>
#include <assert.h>

// for backwards compatability
#ifdef HAVE_clAmdBlas
#define HAVE_clBLAS
#endif

// each implementation of MAGMA defines HAVE_* appropriately.
#if ! defined(HAVE_CUBLAS) && ! defined(HAVE_clBLAS) && ! defined(HAVE_MIC)
#define HAVE_CUBLAS
#endif

// =============================================================================
// C99 standard defines __func__. Some older compilers use __FUNCTION__.
// Note __func__ in C99 is not a macro, so ifndef __func__ doesn't work.
#if __STDC_VERSION__ < 199901L
  #ifndef __func__
    #if __GNUC__ >= 2 || _MSC_VER >= 1300
      #define __func__ __FUNCTION__
    #else
      #define __func__ "<unknown>"
    #endif
  #endif
#endif

// =============================================================================
// To use int64_t, link with mkl_intel_ilp64 or similar (instead of mkl_intel_lp64).
// Similar to magma_int_t we declare magma_index_t used for row/column indices in sparse
#if defined(MAGMA_ILP64) || defined(MKL_ILP64)
typedef long long int magma_int_t;  // MKL uses long long int, not int64_t
#else
typedef int magma_int_t;
#endif

typedef int magma_index_t;
typedef unsigned int magma_uindex_t;

// Define new type that the precision generator will not change (matches PLASMA)
typedef double real_Double_t;

// =============================================================================
// define types specific to implementation (CUDA, OpenCL, MIC)
// define macros to deal with complex numbers
#if defined(HAVE_CUBLAS)
    // include cublas_v2.h, unless cublas.h has already been included, e.g., via magma.h
    #ifndef CUBLAS_H_
    #include <cuda.h>    // for CUDA_VERSION
    #include <cublas_v2.h>
    #endif

    #include <cusparse_v2.h>

    #ifdef __cplusplus
    extern "C" {
    #endif

    // opaque queue structure
    struct magma_queue;
    typedef struct magma_queue* magma_queue_t;
    typedef cudaEvent_t    magma_event_t;
    typedef magma_int_t    magma_device_t;

    // Half precision in CUDA
    #if defined(__cplusplus) && CUDA_VERSION >= 7500
    #include <cuda_fp16.h>
    typedef __half           magmaHalf;
    #else
    // use short for cuda older than 7.5
    // corresponding routines would not work anyway since there is no half precision
    typedef short            magmaHalf;
    #endif    // CUDA_VERSION >= 7500

    typedef cuDoubleComplex magmaDoubleComplex;
    typedef cuFloatComplex  magmaFloatComplex;

    cudaStream_t     magma_queue_get_cuda_stream    ( magma_queue_t queue );
    cublasHandle_t   magma_queue_get_cublas_handle  ( magma_queue_t queue );
    cusparseHandle_t magma_queue_get_cusparse_handle( magma_queue_t queue );

    /// @addtogroup magma_complex
    /// @{

    #define MAGMA_Z_MAKE(r,i)     make_cuDoubleComplex(r, i)    ///< @return complex number r + i*sqrt(-1).
    #define MAGMA_Z_REAL(a)       (a).x                         ///< @return real component of a.
    #define MAGMA_Z_IMAG(a)       (a).y                         ///< @return imaginary component of a.
    #define MAGMA_Z_ADD(a, b)     cuCadd(a, b)                  ///< @return (a + b).
    #define MAGMA_Z_SUB(a, b)     cuCsub(a, b)                  ///< @return (a - b).
    #define MAGMA_Z_MUL(a, b)     cuCmul(a, b)                  ///< @return (a * b).
    #define MAGMA_Z_DIV(a, b)     cuCdiv(a, b)                  ///< @return (a / b).
    #define MAGMA_Z_ABS(a)        cuCabs(a)                     ///< @return absolute value, |a| = sqrt( real(a)^2 + imag(a)^2 ).
    #define MAGMA_Z_ABS1(a)       (fabs((a).x) + fabs((a).y))   ///< @return 1-norm absolute value, | real(a) | + | imag(a) |.
    #define MAGMA_Z_CONJ(a)       cuConj(a)                     ///< @return conjugate of a.

    #define MAGMA_C_MAKE(r,i)     make_cuFloatComplex(r, i)
    #define MAGMA_C_REAL(a)       (a).x
    #define MAGMA_C_IMAG(a)       (a).y
    #define MAGMA_C_ADD(a, b)     cuCaddf(a, b)
    #define MAGMA_C_SUB(a, b)     cuCsubf(a, b)
    #define MAGMA_C_MUL(a, b)     cuCmulf(a, b)
    #define MAGMA_C_DIV(a, b)     cuCdivf(a, b)
    #define MAGMA_C_ABS(a)        cuCabsf(a)
    #define MAGMA_C_ABS1(a)       (fabsf((a).x) + fabsf((a).y))
    #define MAGMA_C_CONJ(a)       cuConjf(a)

    /// @}
    // end group magma_complex

    #ifdef __cplusplus
    }
    #endif
#elif defined(HAVE_clBLAS)
    #include <clBLAS.h>

    #ifdef __cplusplus
    extern "C" {
    #endif

    typedef cl_command_queue  magma_queue_t;
    typedef cl_event          magma_event_t;
    typedef cl_device_id      magma_device_t;

    typedef short         magmaHalf;    // placeholder until FP16 is supported
    typedef DoubleComplex magmaDoubleComplex;
    typedef FloatComplex  magmaFloatComplex;

    cl_command_queue magma_queue_get_cl_queue( magma_queue_t queue );

    #define MAGMA_Z_MAKE(r,i)     doubleComplex(r,i)
    #define MAGMA_Z_REAL(a)       (a).s[0]
    #define MAGMA_Z_IMAG(a)       (a).s[1]
    #define MAGMA_Z_ADD(a, b)     MAGMA_Z_MAKE((a).s[0] + (b).s[0], (a).s[1] + (b).s[1])
    #define MAGMA_Z_SUB(a, b)     MAGMA_Z_MAKE((a).s[0] - (b).s[0], (a).s[1] - (b).s[1])
    #define MAGMA_Z_MUL(a, b)     ((a) * (b))
    #define MAGMA_Z_DIV(a, b)     ((a) / (b))
    #define MAGMA_Z_ABS(a)        magma_cabs(a)
    #define MAGMA_Z_ABS1(a)       (fabs((a).s[0]) + fabs((a).s[1]))
    #define MAGMA_Z_CONJ(a)       MAGMA_Z_MAKE((a).s[0], -(a).s[1])

    #define MAGMA_C_MAKE(r,i)     floatComplex(r,i)
    #define MAGMA_C_REAL(a)       (a).s[0]
    #define MAGMA_C_IMAG(a)       (a).s[1]
    #define MAGMA_C_ADD(a, b)     MAGMA_C_MAKE((a).s[0] + (b).s[0], (a).s[1] + (b).s[1])
    #define MAGMA_C_SUB(a, b)     MAGMA_C_MAKE((a).s[0] - (b).s[0], (a).s[1] - (b).s[1])
    #define MAGMA_C_MUL(a, b)     ((a) * (b))
    #define MAGMA_C_DIV(a, b)     ((a) / (b))
    #define MAGMA_C_ABS(a)        magma_cabsf(a)
    #define MAGMA_C_ABS1(a)       (fabsf((a).s[0]) + fabsf((a).s[1]))
    #define MAGMA_C_CONJ(a)       MAGMA_C_MAKE((a).s[0], -(a).s[1])

    #ifdef __cplusplus
    }
    #endif
#elif defined(HAVE_MIC)
    #include <complex>

    #ifdef __cplusplus
    extern "C" {
    #endif

    typedef int   magma_queue_t;
    typedef int   magma_event_t;
    typedef int   magma_device_t;

    typedef short                 magmaHalf;    // placeholder until FP16 is supported
    typedef std::complex<float>   magmaFloatComplex;
    typedef std::complex<double>  magmaDoubleComplex;

    #define MAGMA_Z_MAKE(r, i)    std::complex<double>(r,i)
    #define MAGMA_Z_REAL(x)       (x).real()
    #define MAGMA_Z_IMAG(x)       (x).imag()
    #define MAGMA_Z_ADD(a, b)     ((a)+(b))
    #define MAGMA_Z_SUB(a, b)     ((a)-(b))
    #define MAGMA_Z_MUL(a, b)     ((a)*(b))
    #define MAGMA_Z_DIV(a, b)     ((a)/(b))
    #define MAGMA_Z_ABS(a)        abs(a)
    #define MAGMA_Z_ABS1(a)       (fabs((a).real()) + fabs((a).imag()))
    #define MAGMA_Z_CONJ(a)       conj(a)

    #define MAGMA_C_MAKE(r, i)    std::complex<float> (r,i)
    #define MAGMA_C_REAL(x)       (x).real()
    #define MAGMA_C_IMAG(x)       (x).imag()
    #define MAGMA_C_ADD(a, b)     ((a)+(b))
    #define MAGMA_C_SUB(a, b)     ((a)-(b))
    #define MAGMA_C_MUL(a, b)     ((a)*(b))
    #define MAGMA_C_DIV(a, b)     ((a)/(b))
    #define MAGMA_C_ABS(a)        abs(a)
    #define MAGMA_C_ABS1(a)       (fabs((a).real()) + fabs((a).imag()))
    #define MAGMA_C_CONJ(a)       conj(a)

    #ifdef __cplusplus
    }
    #endif
#else
    #error "One of HAVE_CUBLAS, HAVE_clBLAS, or HAVE_MIC must be defined. For example, add -DHAVE_CUBLAS to CFLAGS, or #define HAVE_CUBLAS before #include <magma.h>. In MAGMA, this happens in Makefile."
#endif

#ifdef __cplusplus
extern "C" {
#endif

#define MAGMA_Z_EQUAL(a,b)        (MAGMA_Z_REAL(a)==MAGMA_Z_REAL(b) && MAGMA_Z_IMAG(a)==MAGMA_Z_IMAG(b))
#define MAGMA_Z_NEGATE(a)         MAGMA_Z_MAKE( -MAGMA_Z_REAL(a), -MAGMA_Z_IMAG(a))

#define MAGMA_C_EQUAL(a,b)        (MAGMA_C_REAL(a)==MAGMA_C_REAL(b) && MAGMA_C_IMAG(a)==MAGMA_C_IMAG(b))
#define MAGMA_C_NEGATE(a)         MAGMA_C_MAKE( -MAGMA_C_REAL(a), -MAGMA_C_IMAG(a))

#define MAGMA_D_MAKE(r,i)         (r)
#define MAGMA_D_REAL(x)           (x)
#define MAGMA_D_IMAG(x)           (0.0)
#define MAGMA_D_ADD(a, b)         ((a) + (b))
#define MAGMA_D_SUB(a, b)         ((a) - (b))
#define MAGMA_D_MUL(a, b)         ((a) * (b))
#define MAGMA_D_DIV(a, b)         ((a) / (b))
#define MAGMA_D_ABS(a)            ((a)>0 ? (a) : -(a))
#define MAGMA_D_ABS1(a)           ((a)>0 ? (a) : -(a))
#define MAGMA_D_CONJ(a)           (a)
#define MAGMA_D_EQUAL(a,b)        ((a) == (b))
#define MAGMA_D_NEGATE(a)         (-a)

#define MAGMA_S_MAKE(r,i)         (r)
#define MAGMA_S_REAL(x)           (x)
#define MAGMA_S_IMAG(x)           (0.0)
#define MAGMA_S_ADD(a, b)         ((a) + (b))
#define MAGMA_S_SUB(a, b)         ((a) - (b))
#define MAGMA_S_MUL(a, b)         ((a) * (b))
#define MAGMA_S_DIV(a, b)         ((a) / (b))
#define MAGMA_S_ABS(a)            ((a)>0 ? (a) : -(a))
#define MAGMA_S_ABS1(a)           ((a)>0 ? (a) : -(a))
#define MAGMA_S_CONJ(a)           (a)
#define MAGMA_S_EQUAL(a,b)        ((a) == (b))
#define MAGMA_S_NEGATE(a)         (-a)

#define MAGMA_Z_ZERO              MAGMA_Z_MAKE( 0.0, 0.0)
#define MAGMA_Z_ONE               MAGMA_Z_MAKE( 1.0, 0.0)
#define MAGMA_Z_HALF              MAGMA_Z_MAKE( 0.5, 0.0)
#define MAGMA_Z_NEG_ONE           MAGMA_Z_MAKE(-1.0, 0.0)
#define MAGMA_Z_NEG_HALF          MAGMA_Z_MAKE(-0.5, 0.0)

#define MAGMA_C_ZERO              MAGMA_C_MAKE( 0.0, 0.0)
#define MAGMA_C_ONE               MAGMA_C_MAKE( 1.0, 0.0)
#define MAGMA_C_HALF              MAGMA_C_MAKE( 0.5, 0.0)
#define MAGMA_C_NEG_ONE           MAGMA_C_MAKE(-1.0, 0.0)
#define MAGMA_C_NEG_HALF          MAGMA_C_MAKE(-0.5, 0.0)

#define MAGMA_D_ZERO              ( 0.0)
#define MAGMA_D_ONE               ( 1.0)
#define MAGMA_D_HALF              ( 0.5)
#define MAGMA_D_NEG_ONE           (-1.0)
#define MAGMA_D_NEG_HALF          (-0.5)

#define MAGMA_S_ZERO              ( 0.0)
#define MAGMA_S_ONE               ( 1.0)
#define MAGMA_S_HALF              ( 0.5)
#define MAGMA_S_NEG_ONE           (-1.0)
#define MAGMA_S_NEG_HALF          (-0.5)

#ifndef CBLAS_SADDR
#define CBLAS_SADDR(a)  &(a)
#endif

// for MAGMA_[CZ]_ABS
double magma_cabs ( magmaDoubleComplex x );
float  magma_cabsf( magmaFloatComplex  x );

#if defined(HAVE_clBLAS)
    // OpenCL uses opaque memory references on GPU
    typedef cl_mem magma_ptr;
    typedef cl_mem magmaInt_ptr;
    typedef cl_mem magmaIndex_ptr;
    typedef cl_mem magmaFloat_ptr;
    typedef cl_mem magmaDouble_ptr;
    typedef cl_mem magmaFloatComplex_ptr;
    typedef cl_mem magmaDoubleComplex_ptr;

    typedef cl_mem magma_const_ptr;
    typedef cl_mem magmaInt_const_ptr;
    typedef cl_mem magmaIndex_const_ptr;
    typedef cl_mem magmaFloat_const_ptr;
    typedef cl_mem magmaDouble_const_ptr;
    typedef cl_mem magmaFloatComplex_const_ptr;
    typedef cl_mem magmaDoubleComplex_const_ptr;
#else
    // MIC and CUDA use regular pointers on GPU
    typedef void               *magma_ptr;
    typedef magma_int_t        *magmaInt_ptr;
    typedef magma_index_t      *magmaIndex_ptr;
    typedef magma_uindex_t     *magmaUIndex_ptr;
    typedef float              *magmaFloat_ptr;
    typedef double             *magmaDouble_ptr;
    typedef magmaFloatComplex  *magmaFloatComplex_ptr;
    typedef magmaDoubleComplex *magmaDoubleComplex_ptr;
    typedef magmaHalf          *magmaHalf_ptr;

    typedef void               const *magma_const_ptr;
    typedef magma_int_t        const *magmaInt_const_ptr;
    typedef magma_index_t      const *magmaIndex_const_ptr;
    typedef magma_uindex_t     const *magmaUIndex_const_ptr;
    typedef float              const *magmaFloat_const_ptr;
    typedef double             const *magmaDouble_const_ptr;
    typedef magmaFloatComplex  const *magmaFloatComplex_const_ptr;
    typedef magmaDoubleComplex const *magmaDoubleComplex_const_ptr;
    typedef magmaHalf          const *magmaHalf_const_ptr;
#endif

// =============================================================================
// MAGMA constants

// -----------------------------------------------------------------------------
#define MAGMA_VERSION_MAJOR 2
#define MAGMA_VERSION_MINOR 5
#define MAGMA_VERSION_MICRO 4

// stage is "svn", "beta#", "rc#" (release candidate), or blank ("") for final release
#define MAGMA_VERSION_STAGE ""

#define MagmaMaxGPUs 8
#define MagmaMaxAccelerators 8
#define MagmaMaxSubs 16

// trsv template parameter
#define MagmaBigTileSize 1000000

// -----------------------------------------------------------------------------
// Return codes
// LAPACK argument errors are < 0 but > MAGMA_ERR.
// MAGMA errors are < MAGMA_ERR.
/// @addtogroup magma_error_codes
/// @{

#define MAGMA_SUCCESS               0       ///< operation was successful
#define MAGMA_ERR                  -100     ///< unspecified error
#define MAGMA_ERR_NOT_INITIALIZED  -101     ///< magma_init() was not called
#define MAGMA_ERR_REINITIALIZED    -102     // unused
#define MAGMA_ERR_NOT_SUPPORTED    -103     ///< not supported on this GPU
#define MAGMA_ERR_ILLEGAL_VALUE    -104     // unused
#define MAGMA_ERR_NOT_FOUND        -105     ///< file not found
#define MAGMA_ERR_ALLOCATION       -106     // unused
#define MAGMA_ERR_INTERNAL_LIMIT   -107     // unused
#define MAGMA_ERR_UNALLOCATED      -108     // unused
#define MAGMA_ERR_FILESYSTEM       -109     // unused
#define MAGMA_ERR_UNEXPECTED       -110     // unused
#define MAGMA_ERR_SEQUENCE_FLUSHED -111     // unused
#define MAGMA_ERR_HOST_ALLOC       -112     ///< could not malloc CPU host memory
#define MAGMA_ERR_DEVICE_ALLOC     -113     ///< could not malloc GPU device memory
#define MAGMA_ERR_CUDASTREAM       -114     // unused
#define MAGMA_ERR_INVALID_PTR      -115     ///< can't free invalid pointer
#define MAGMA_ERR_UNKNOWN          -116     ///< unspecified error
#define MAGMA_ERR_NOT_IMPLEMENTED  -117     ///< not implemented yet
#define MAGMA_ERR_NAN              -118     ///< NaN (not-a-number) detected

// some MAGMA-sparse errors
#define MAGMA_SLOW_CONVERGENCE     -201
#define MAGMA_DIVERGENCE           -202
#define MAGMA_NONSPD               -203
#define MAGMA_ERR_BADPRECOND       -204
#define MAGMA_NOTCONVERGED         -205

// When adding error codes, please add to interface_cuda/error.cpp

// map cusparse errors to magma errors
#define MAGMA_ERR_CUSPARSE                            -3000
#define MAGMA_ERR_CUSPARSE_NOT_INITIALIZED            -3001
#define MAGMA_ERR_CUSPARSE_ALLOC_FAILED               -3002
#define MAGMA_ERR_CUSPARSE_INVALID_VALUE              -3003
#define MAGMA_ERR_CUSPARSE_ARCH_MISMATCH              -3004
#define MAGMA_ERR_CUSPARSE_MAPPING_ERROR              -3005
#define MAGMA_ERR_CUSPARSE_EXECUTION_FAILED           -3006
#define MAGMA_ERR_CUSPARSE_INTERNAL_ERROR             -3007
#define MAGMA_ERR_CUSPARSE_MATRIX_TYPE_NOT_SUPPORTED  -3008
#define MAGMA_ERR_CUSPARSE_ZERO_PIVOT                 -3009

/// @}
// end group magma_error_codes

// -----------------------------------------------------------------------------
// parameter constants
// numbering is consistent with CBLAS and PLASMA; see plasma/include/plasma.h
// also with lapack_cwrapper/include/lapack_enum.h
// see http://www.netlib.org/lapack/lapwrapc/
typedef enum {
    MagmaFalse         = 0,
    MagmaTrue          = 1
} magma_bool_t;

typedef enum {
    MagmaRowMajor      = 101,
    MagmaColMajor      = 102
} magma_order_t;

// Magma_ConjTrans is an alias for those rare occasions (zlarfb, zun*, zher*k)
// where we want Magma_ConjTrans to convert to MagmaTrans in precision generation.
typedef enum {
    MagmaNoTrans       = 111,
    MagmaTrans         = 112,
    MagmaConjTrans     = 113,
    Magma_ConjTrans    = MagmaConjTrans
} magma_trans_t;

typedef enum {
    MagmaUpper         = 121,
    MagmaLower         = 122,
    MagmaFull          = 123,  /* lascl, laset */
    MagmaHessenberg    = 124   /* lascl */
} magma_uplo_t;

typedef magma_uplo_t magma_type_t;  /* lascl */

typedef enum {
    MagmaNonUnit       = 131,
    MagmaUnit          = 132
} magma_diag_t;

typedef enum {
    MagmaLeft          = 141,
    MagmaRight         = 142,
    MagmaBothSides     = 143   /* trevc */
} magma_side_t;

typedef enum {
    MagmaOneNorm       = 171,  /* lange, lanhe */
    MagmaRealOneNorm   = 172,
    MagmaTwoNorm       = 173,
    MagmaFrobeniusNorm = 174,
    MagmaInfNorm       = 175,
    MagmaRealInfNorm   = 176,
    MagmaMaxNorm       = 177,
    MagmaRealMaxNorm   = 178
} magma_norm_t;

typedef enum {
    MagmaDistUniform   = 201,  /* latms */
    MagmaDistSymmetric = 202,
    MagmaDistNormal    = 203
} magma_dist_t;

typedef enum {
    MagmaHermGeev      = 241,  /* latms */
    MagmaHermPoev      = 242,
    MagmaNonsymPosv    = 243,
    MagmaSymPosv       = 244
} magma_sym_t;

typedef enum {
    MagmaNoPacking     = 291,  /* latms */
    MagmaPackSubdiag   = 292,
    MagmaPackSupdiag   = 293,
    MagmaPackColumn    = 294,
    MagmaPackRow       = 295,
    MagmaPackLowerBand = 296,
    MagmaPackUpeprBand = 297,
    MagmaPackAll       = 298
} magma_pack_t;

typedef enum {
    MagmaNoVec         = 301,  /* geev, syev, gesvd */
    MagmaVec           = 302,  /* geev, syev */
    MagmaIVec          = 303,  /* stedc */
    MagmaAllVec        = 304,  /* gesvd, trevc */
    MagmaSomeVec       = 305,  /* gesvd, trevc */
    MagmaOverwriteVec  = 306,  /* gesvd */
    MagmaBacktransVec  = 307   /* trevc */
} magma_vec_t;

typedef enum {
    MagmaRangeAll      = 311,  /* syevx, etc. */
    MagmaRangeV        = 312,
    MagmaRangeI        = 313
} magma_range_t;

typedef enum {
    MagmaQ             = 322,  /* unmbr, ungbr */
    MagmaP             = 323
} magma_vect_t;

typedef enum {
    MagmaForward       = 391,  /* larfb */
    MagmaBackward      = 392
} magma_direct_t;

typedef enum {
    MagmaColumnwise    = 401,  /* larfb */
    MagmaRowwise       = 402
} magma_storev_t;

typedef enum {
    MagmaHybrid        = 701,
    MagmaNative        = 702
} magma_mode_t;
// -----------------------------------------------------------------------------
// sparse
typedef enum {
    Magma_CSR          = 611,
    Magma_ELLPACKT     = 612,
    Magma_ELL          = 613,
    Magma_DENSE        = 614,
    Magma_BCSR         = 615,
    Magma_CSC          = 616,
    Magma_HYB          = 617,
    Magma_COO          = 618,
    Magma_ELLRT        = 619,
    Magma_SPMVFUNCTION = 620,
    Magma_SELLP        = 621,
    Magma_ELLD         = 622,
    Magma_CSRLIST      = 623,
    Magma_CSRD         = 624,
    Magma_CSRL         = 627,
    Magma_CSRU         = 628,
    Magma_CSRCOO       = 629,
    Magma_CUCSR        = 630,
    Magma_COOLIST      = 631,
    Magma_CSR5         = 632
} magma_storage_t;

typedef enum {
    Magma_CG           = 431,
    Magma_CGMERGE      = 432,
    Magma_GMRES        = 433,
    Magma_BICGSTAB     = 434,
  Magma_BICGSTABMERGE  = 435,
  Magma_BICGSTABMERGE2 = 436,
    Magma_JACOBI       = 437,
    Magma_GS           = 438,
    Magma_ITERREF      = 439,
    Magma_BCSRLU       = 440,
    Magma_PCG          = 441,
    Magma_PGMRES       = 442,
    Magma_PBICGSTAB    = 443,
    Magma_PASTIX       = 444,
    Magma_ILU          = 445,
    Magma_ICC          = 446,
    Magma_PARILU       = 447,
    Magma_PARIC        = 448,
    Magma_BAITER       = 449,
    Magma_LOBPCG       = 450,
    Magma_NONE         = 451,
    Magma_FUNCTION     = 452,
    Magma_IDR          = 453,
    Magma_PIDR         = 454,
    Magma_CGS          = 455,
    Magma_PCGS         = 456,
    Magma_CGSMERGE     = 457,
    Magma_PCGSMERGE    = 458,
    Magma_TFQMR        = 459,
    Magma_PTFQMR       = 460,
    Magma_TFQMRMERGE   = 461,
    Magma_PTFQMRMERGE  = 462,
    Magma_QMR          = 463,
    Magma_PQMR         = 464,
    Magma_QMRMERGE     = 465,
    Magma_PQMRMERGE    = 466,
    Magma_BOMBARD      = 490,
    Magma_BOMBARDMERGE = 491,
    Magma_PCGMERGE     = 492,
    Magma_BAITERO      = 493,
    Magma_IDRMERGE     = 494,
  Magma_PBICGSTABMERGE = 495,
    Magma_PARICT       = 496,
    Magma_CUSTOMIC     = 497,
    Magma_CUSTOMILU    = 498,
    Magma_PIDRMERGE    = 499,
    Magma_BICG         = 500,
    Magma_BICGMERGE    = 501,
    Magma_PBICG        = 502,
    Magma_PBICGMERGE   = 503,
    Magma_LSQR         = 504,
    Magma_PARILUT      = 505,
    Magma_ISAI         = 506,
    Magma_CUSOLVE      = 507,
    Magma_VBJACOBI     = 508,
    Magma_PARDISO      = 509,
    Magma_SYNCFREESOLVE= 510,
    Magma_ILUT         = 511
} magma_solver_type;

typedef enum {
    Magma_CGSO         = 561,
    Magma_FUSED_CGSO   = 562,
    Magma_MGSO         = 563
} magma_ortho_t;

typedef enum {
    Magma_CPU          = 571,
    Magma_DEV          = 572
} magma_location_t;

typedef enum {
    Magma_GENERAL      = 581,
    Magma_SYMMETRIC    = 582
} magma_symmetry_t;

typedef enum {
    Magma_ORDERED      = 591,
    Magma_DIAGFIRST    = 592,
    Magma_UNITY        = 593,
    Magma_VALUE        = 594
} magma_diagorder_t;

typedef enum {
    Magma_DCOMPLEX     = 501,
    Magma_FCOMPLEX     = 502,
    Magma_DOUBLE       = 503,
    Magma_FLOAT        = 504
} magma_precision;

typedef enum {
    Magma_NOSCALE      = 511,
    Magma_UNITROW      = 512,
    Magma_UNITDIAG     = 513,
    Magma_UNITCOL      = 514,
    Magma_UNITROWCOL   = 515, // to be deprecated
    Magma_UNITDIAGCOL  = 516, // to be deprecated
} magma_scale_t;

typedef enum {
    Magma_SOLVE        = 801,
    Magma_SETUPSOLVE   = 802,
    Magma_APPLYSOLVE   = 803,
    Magma_DESTROYSOLVE = 804,
    Magma_INFOSOLVE    = 805,
    Magma_GENERATEPREC = 806,
    Magma_PRECONDLEFT  = 807,
    Magma_PRECONDRIGHT = 808,
    Magma_TRANSPOSE    = 809,
    Magma_SPMV         = 810
} magma_operation_t;

typedef enum {
    Magma_PREC_SS           = 900,
    Magma_PREC_SST          = 901,
    Magma_PREC_HS           = 902,
    Magma_PREC_HST          = 903,
    Magma_PREC_SH           = 904,
    Magma_PREC_SHT          = 905,

    Magma_PREC_XHS_H        = 910,
    Magma_PREC_XHS_HTC      = 911,
    Magma_PREC_XHS_161616   = 912,
    Magma_PREC_XHS_161616TC = 913,
    Magma_PREC_XHS_161632TC = 914,
    Magma_PREC_XSH_S        = 915,
    Magma_PREC_XSH_STC      = 916,
    Magma_PREC_XSH_163232TC = 917,
    Magma_PREC_XSH_323232TC = 918,

    Magma_REFINE_IRSTRS   = 920,
    Magma_REFINE_IRDTRS   = 921,
    Magma_REFINE_IRGMSTRS = 922,
    Magma_REFINE_IRGMDTRS = 923,
    Magma_REFINE_GMSTRS   = 924,
    Magma_REFINE_GMDTRS   = 925,
    Magma_REFINE_GMGMSTRS = 926,
    Magma_REFINE_GMGMDTRS = 927,

    Magma_PREC_HD         = 930,
} magma_refinement_t;

typedef enum {
    Magma_MP_BASE_SS              = 950,
    Magma_MP_BASE_DD              = 951,
    Magma_MP_BASE_XHS             = 952,
    Magma_MP_BASE_XSH             = 953,
    Magma_MP_BASE_XHD             = 954,
    Magma_MP_BASE_XDH             = 955,

    Magma_MP_ENABLE_DFLT_MATH     = 960,
    Magma_MP_ENABLE_TC_MATH       = 961,
    Magma_MP_SGEMM                = 962,
    Magma_MP_HGEMM                = 963,
    Magma_MP_GEMEX_I32_O32_C32    = 964,
    Magma_MP_GEMEX_I16_O32_C32    = 965,
    Magma_MP_GEMEX_I16_O16_C32    = 966,
    Magma_MP_GEMEX_I16_O16_C16    = 967,

    Magma_MP_TC_SGEMM             = 968,
    Magma_MP_TC_HGEMM             = 969,
    Magma_MP_TC_GEMEX_I32_O32_C32 = 970,
    Magma_MP_TC_GEMEX_I16_O32_C32 = 971,
    Magma_MP_TC_GEMEX_I16_O16_C32 = 972,
    Magma_MP_TC_GEMEX_I16_O16_C16 = 973,

} magma_mp_type_t;

// When adding constants, remember to do these steps as appropriate:
// 1)  add magma_xxxx_const()  converter below and in control/constants.cpp
// 2a) add to magma2lapack_constants[] in control/constants.cpp
// 2b) update min & max here, which are used to check bounds for magma2lapack_constants[]
// 2c) add lapack_xxxx_const() converter below and in control/constants.cpp
#define Magma2lapack_Min  MagmaFalse     // 0
#define Magma2lapack_Max  MagmaRowwise   // 402

// -----------------------------------------------------------------------------
// string constants for calling Fortran BLAS and LAPACK
// todo: use translators instead? lapack_const_str( MagmaUpper )
#define MagmaRowMajorStr      "Row"
#define MagmaColMajorStr      "Col"

#define MagmaNoTransStr       "NoTrans"
#define MagmaTransStr         "Trans"
#define MagmaConjTransStr     "ConjTrans"
#define Magma_ConjTransStr    "ConjTrans"

#define MagmaUpperStr         "Upper"
#define MagmaLowerStr         "Lower"
#define MagmaFullStr          "Full"

#define MagmaNonUnitStr       "NonUnit"
#define MagmaUnitStr          "Unit"

#define MagmaLeftStr          "Left"
#define MagmaRightStr         "Right"
#define MagmaBothSidesStr     "Both"

#define MagmaOneNormStr       "1"
#define MagmaTwoNormStr       "2"
#define MagmaFrobeniusNormStr "Fro"
#define MagmaInfNormStr       "Inf"
#define MagmaMaxNormStr       "Max"

#define MagmaForwardStr       "Forward"
#define MagmaBackwardStr      "Backward"

#define MagmaColumnwiseStr    "Columnwise"
#define MagmaRowwiseStr       "Rowwise"

#define MagmaNoVecStr         "NoVec"
#define MagmaVecStr           "Vec"
#define MagmaIVecStr          "IVec"
#define MagmaAllVecStr        "All"
#define MagmaSomeVecStr       "Some"
#define MagmaOverwriteVecStr  "Overwrite"

// -----------------------------------------------------------------------------
// Convert LAPACK character constants to MAGMA constants.
// This is a one-to-many mapping, requiring multiple translators
// (e.g., "N" can be NoTrans or NonUnit or NoVec).
magma_bool_t   magma_bool_const  ( char lapack_char );
magma_order_t  magma_order_const ( char lapack_char );
magma_trans_t  magma_trans_const ( char lapack_char );
magma_uplo_t   magma_uplo_const  ( char lapack_char );
magma_diag_t   magma_diag_const  ( char lapack_char );
magma_side_t   magma_side_const  ( char lapack_char );
magma_norm_t   magma_norm_const  ( char lapack_char );
magma_dist_t   magma_dist_const  ( char lapack_char );
magma_sym_t    magma_sym_const   ( char lapack_char );
magma_pack_t   magma_pack_const  ( char lapack_char );
magma_vec_t    magma_vec_const   ( char lapack_char );
magma_range_t  magma_range_const ( char lapack_char );
magma_vect_t   magma_vect_const  ( char lapack_char );
magma_direct_t magma_direct_const( char lapack_char );
magma_storev_t magma_storev_const( char lapack_char );

// -----------------------------------------------------------------------------
// Convert MAGMA constants to LAPACK(E) constants.
// The generic lapack_const_str works for all cases, but the specific routines
// (e.g., lapack_trans_const) do better error checking.

// magma  defines lapack_const_str, which returns char* to call lapack (Fortran interface).
// plasma defines lapack_const, which is roughly the same as MAGMA's lapacke_const
// (returns a char instead of char*) to call lapacke (C interface).

const char* lapack_const_str   ( int            magma_const );
const char* lapack_bool_const  ( magma_bool_t   magma_const );
const char* lapack_order_const ( magma_order_t  magma_const );
const char* lapack_trans_const ( magma_trans_t  magma_const );
const char* lapack_uplo_const  ( magma_uplo_t   magma_const );
const char* lapack_diag_const  ( magma_diag_t   magma_const );
const char* lapack_side_const  ( magma_side_t   magma_const );
const char* lapack_norm_const  ( magma_norm_t   magma_const );
const char* lapack_dist_const  ( magma_dist_t   magma_const );
const char* lapack_sym_const   ( magma_sym_t    magma_const );
const char* lapack_pack_const  ( magma_pack_t   magma_const );
const char* lapack_vec_const   ( magma_vec_t    magma_const );
const char* lapack_range_const ( magma_range_t  magma_const );
const char* lapack_vect_const  ( magma_vect_t   magma_const );
const char* lapack_direct_const( magma_direct_t magma_const );
const char* lapack_storev_const( magma_storev_t magma_const );

static inline char lapacke_const       ( int magma_const            ) { return *lapack_const_str   ( magma_const ); }
static inline char lapacke_bool_const  ( magma_bool_t   magma_const ) { return *lapack_bool_const  ( magma_const ); }
static inline char lapacke_order_const ( magma_order_t  magma_const ) { return *lapack_order_const ( magma_const ); }
static inline char lapacke_trans_const ( magma_trans_t  magma_const ) { return *lapack_trans_const ( magma_const ); }
static inline char lapacke_uplo_const  ( magma_uplo_t   magma_const ) { return *lapack_uplo_const  ( magma_const ); }
static inline char lapacke_diag_const  ( magma_diag_t   magma_const ) { return *lapack_diag_const  ( magma_const ); }
static inline char lapacke_side_const  ( magma_side_t   magma_const ) { return *lapack_side_const  ( magma_const ); }
static inline char lapacke_norm_const  ( magma_norm_t   magma_const ) { return *lapack_norm_const  ( magma_const ); }
static inline char lapacke_dist_const  ( magma_dist_t   magma_const ) { return *lapack_dist_const  ( magma_const ); }
static inline char lapacke_sym_const   ( magma_sym_t    magma_const ) { return *lapack_sym_const   ( magma_const ); }
static inline char lapacke_pack_const  ( magma_pack_t   magma_const ) { return *lapack_pack_const  ( magma_const ); }
static inline char lapacke_vec_const   ( magma_vec_t    magma_const ) { return *lapack_vec_const   ( magma_const ); }
static inline char lapacke_range_const ( magma_range_t  magma_const ) { return *lapack_range_const ( magma_const ); }
static inline char lapacke_vect_const  ( magma_vect_t   magma_const ) { return *lapack_vect_const  ( magma_const ); }
static inline char lapacke_direct_const( magma_direct_t magma_const ) { return *lapack_direct_const( magma_const ); }
static inline char lapacke_storev_const( magma_storev_t magma_const ) { return *lapack_storev_const( magma_const ); }

// -----------------------------------------------------------------------------
// Convert MAGMA constants to clBLAS constants.
#if defined(HAVE_clBLAS)
clblasOrder          clblas_order_const( magma_order_t order );
clblasTranspose      clblas_trans_const( magma_trans_t trans );
clblasUplo           clblas_uplo_const ( magma_uplo_t  uplo  );
clblasDiag           clblas_diag_const ( magma_diag_t  diag  );
clblasSide           clblas_side_const ( magma_side_t  side  );
#endif

// -----------------------------------------------------------------------------
// Convert MAGMA constants to CUBLAS constants.
#if defined(CUBLAS_V2_H_)
cublasOperation_t    cublas_trans_const ( magma_trans_t trans );
cublasFillMode_t     cublas_uplo_const  ( magma_uplo_t  uplo  );
cublasDiagType_t     cublas_diag_const  ( magma_diag_t  diag  );
cublasSideMode_t     cublas_side_const  ( magma_side_t  side  );
#endif

// -----------------------------------------------------------------------------
// Convert MAGMA constants to CBLAS constants.
#if defined(HAVE_CBLAS)
#include <cblas.h>
enum CBLAS_ORDER     cblas_order_const  ( magma_order_t order );
enum CBLAS_TRANSPOSE cblas_trans_const  ( magma_trans_t trans );
enum CBLAS_UPLO      cblas_uplo_const   ( magma_uplo_t  uplo  );
enum CBLAS_DIAG      cblas_diag_const   ( magma_diag_t  diag  );
enum CBLAS_SIDE      cblas_side_const   ( magma_side_t  side  );
#endif

#ifdef __cplusplus
}
#endif

#endif // MAGMA_TYPES_H

ext: h

max_stars_repo_name: klast/magma

max_issues_repo_name: klast/magma

max_forks_repo_name: klast/magma


text:
#pragma once

#include <gsl/span>

namespace miner {

    template<class T, std::ptrdiff_t Extent = gsl::dynamic_extent>
    using span = gsl::span<T, Extent>;

    template<std::ptrdiff_t Extent = gsl::dynamic_extent>
    using ByteSpan = span<uint8_t, Extent>;

    template<std::ptrdiff_t Extent = gsl::dynamic_extent>
    using cByteSpan = span<const uint8_t, Extent>;
}


text:
/*
   Copyright [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 __MM_PLUGIN_ITF_H__
#define __MM_PLUGIN_ITF_H__

#pragma GCC diagnostic push
#pragma GCC diagnostic ignored "-Wunused-variable"
#pragma GCC diagnostic ignored "-Wunused-parameter"
#pragma GCC diagnostic ignored "-Wunused-function"

#define PUBLIC __attribute__((__visibility__("default")))

#if defined(__cplusplus)
#pragma GCC diagnostic ignored "-Weffc++"
#include <common/byte_span.h>
#include <gsl/span>
#include <functional>
#endif

#include <stdlib.h>
#include <assert.h>

#if defined(__cplusplus)
extern "C"
{
#endif

  typedef int status_t;
  typedef void * mm_plugin_heap_t;
  typedef void (*request_memory_callback_t)(void * param, size_t alignment, size_t size_hint, void * addr_hint);
  typedef void (*release_memory_callback_t)(void * param, void * addr, size_t size);

  /**
   * Initialize mm library
   *
   * @return S_OK, E_FAIL
   */
  status_t mm_plugin_init();

  /**
   * Create a heap instance
   *
   * @param params Constructor parameters (e.g., JSON)
   * @param root_ptr Root point for persistent heaps
   * @param out_heap Heap context
   *
   * @return S_OK, E_FAIL
   */
  status_t mm_plugin_create(const char * params,
                            void * root_ptr,
                            mm_plugin_heap_t * out_heap);

  /**
   * Delete a heap instance
   *
   * @param heap Heap context to destroy
   *
   * @return S_OK, E_INVAL, E_FAIL
   */
  status_t mm_plugin_destroy(mm_plugin_heap_t heap);

  /**
   * Add (slab) region of memory to fuel the allocator.  The passthru
   * allocator does not need this as it takes directly from the OS.
   *
   * @param heap Heap context
   * @param region_base Pointer to beginning of region
   * @param region_size Region length in bytes
   *
   * @return E_NOT_IMPL, S_OK, E_FAIL, E_INVAL
   */
  status_t mm_plugin_add_managed_region(mm_plugin_heap_t heap,
                                        void * region_base,
                                        size_t region_size);

  /**
   * Query memory regions being used by the allocator
   *
   * @param heap Heap context
   * @param region_id Index counting from 0
   * @param [out] Base address of region
   * @param [out] Size of region in bytes
   *
   * @return S_MORE (continue to next region id), S_OK (done), E_INVAL
   */
  status_t mm_plugin_query_managed_region(mm_plugin_heap_t heap,
                                          unsigned region_id,
                                          void** out_region_base,
                                          size_t* out_region_size);

  /**
   * Register callback the allocator can use to request more memory
   *
   * @param heap Heap context
   * @param callback Call back function pointer
   * @param param Optional parameter which will be pass to callback function
   *
   * @return E_NOT_IMPL, S_OK
   */
  status_t mm_plugin_register_callback_request_memory(mm_plugin_heap_t heap,
                                                      request_memory_callback_t callback,
                                                      void * param);

  /**
   * Allocate a region of memory without alignment or hint
   *
   * @param heap Heap context
   * @param Length in bytes
   * @param [out] Pointer to allocated region
   *
   * @return S_OK or E_FAIL
   */
  status_t mm_plugin_allocate(mm_plugin_heap_t heap, size_t n, void ** out_ptr);

  /**
   * Allocation region of memory that is aligned
   *
   * @param heap Heap context
   * @param n Size of region in bytes
   * @param alignment Alignment in bytes
   * @param out_ptr Address if [in] nullptr, [out] pointer to allocated region
   *
   * @return S_OK, E_FAIL, E_INVAL (depending on implementation)
   */
  status_t mm_plugin_aligned_allocate(mm_plugin_heap_t heap, size_t n, size_t alignment, void ** out_ptr);

  /**
   * Special case for EASTL
   *
   * @param heap Heap context
   * @param n Size of region in bytes
   * @param alignment Alignment in bytes
   * @param offset Offset from start of region to location within region which satisfies alignment
   * @param out_ptr Address if [in] nullptr, [out] pointer to allocated region
   *
   * @return
   */
  status_t mm_plugin_aligned_allocate_offset(mm_plugin_heap_t heap, size_t n, size_t alignment, size_t offset, void ** out_ptr);

  /**
   * Free a previously allocated region of memory with length known
   *
   * @param heap Heap context
   * @param ptr Address of [in] pointer to previously allocated region, [out] nullptr
   * @param size Length of region in bytes
   *
   * @return S_OK or E_INVAL;
   */
  status_t mm_plugin_deallocate(mm_plugin_heap_t heap, void ** ptr, size_t size);

  /**
   * Free previously allocated region without known length
   *
   * @param heap Heap context
   * @param ptr Address of [in] pointer to previously allocated region, [out] nullptr
   *
   * @return S_OK
   */
  status_t mm_plugin_deallocate_without_size(mm_plugin_heap_t heap, void ** ptr);

  /**
   * Allocate region and zero memory
   *
   * @param heap Heap context
   * @param size Size of region in bytes
   * @param ptr [out] Pointer to allocated region
   *
   * @return S_OK
   */
  status_t mm_plugin_callocate(mm_plugin_heap_t heap, size_t n, void ** out_ptr);

  /*
    The POSIX realloc() function changes the size of the memory block pointed
    to by ptr to size bytes.  The contents will be unchanged in the
    range from the start of the region up to the minimum of the old and
    new sizes.  If the new size is larger than the old size, the added
    memory will not be initialized.  If ptr is NULL, then the call is
    equivalent to malloc(size), for all values of size; if size is equal
    to zero, and ptr is not NULL, then the call is equivalent to
    free(ptr).  Unless ptr is NULL, it must have been returned by an
    earlier call to malloc(), calloc(), or realloc().  If the area
    pointed to was moved, a free(ptr) is done.
  */

  /**
   * Resize an existing allocation
   *
   * @param heap Heap context
   * @param in_out_ptr Address of pointer to [in] existing allocated region, [out] new reallocated region or null if unable to reallocate
   * @param size New size in bytes
   *
   *
   * @return S_OK
   */
  status_t mm_plugin_reallocate(mm_plugin_heap_t heap, void ** in_out_ptr, size_t size);

  /**
   * Get the number of usable bytes in block pointed to by ptr.  The
   * allocator *may* not have this information and should then return
   * E_NOT_IMPL. Returned size may be larger than requested allocation
   *
   * @param heap Heap context
   * @param ptr Pointer to block base
   * @param out_ptr [out] Number of bytes in allocated block
   *
   * @return S_OK, E_NOT_IMPL
   */
  status_t mm_plugin_usable_size(mm_plugin_heap_t heap, void * ptr, size_t * out_size);

  /**
   * Inject an allocation back into the allocator (reconstituing)
   *
   * @param heap Heap context
   * @param ptr Pointer to region of memory to mark allocated
   * @param size Size of region in bytes
   *
   * @return S_OK, E_NOT_IMPL
   */
  status_t mm_plugin_inject_allocation(mm_plugin_heap_t heap, void * ptr, size_t size);

  /**
   * Report bytes remaining
   *
   * @param bytes_remaining Bytes remaining for allocation
   *
   * @return S_OK, E_NOT_IMPL
   */
  status_t mm_plugin_bytes_remaining(mm_plugin_heap_t heap, size_t *bytes_remaining);

  /**
   * Get debugging information
   *
   * @param heap Heap context
   */
  void mm_plugin_debug(mm_plugin_heap_t heap);

  /**
   * Check for crash consistency libccpm behavior)
   *
   * @return non-zero iff the root_ptr parameter of mm_plugin_create is interpreted as a struct ccpm_params *.
   */
  int mm_plugin_is_crash_consistent(mm_plugin_heap_t heap);

  /**
   * Check for inject_allocation capability
   *
   * @return non-zero iff mm_plugin_inject_allocation is implemented
   */
  int mm_plugin_can_inject_allocation(mm_plugin_heap_t heap);

  /**
   * Function pointer table for all methods
   *
   */
  typedef struct tag_mm_plugin_function_table_t
  {
    status_t (*mm_plugin_init)();
    status_t (*mm_plugin_create)(const char * params, void * root_ptr
			, mm_plugin_heap_t * out_heap);
    status_t (*mm_plugin_destroy)(mm_plugin_heap_t heap);
    status_t (*mm_plugin_add_managed_region)(mm_plugin_heap_t heap,
                                             void * region_base,
                                             size_t region_size);
    status_t (*mm_plugin_query_managed_region)(mm_plugin_heap_t heap,
                                               unsigned region_id,
                                               void** out_region_base,
                                               size_t* out_region_size);
    status_t (*mm_plugin_register_callback_request_memory)(mm_plugin_heap_t heap,
                                                           request_memory_callback_t callback,
                                                           void * param);
    status_t (*mm_plugin_allocate)(mm_plugin_heap_t heap, size_t n, void ** out_ptr);
    status_t (*mm_plugin_aligned_allocate)(mm_plugin_heap_t heap, size_t n, size_t alignment, void ** out_ptr);
    status_t (*mm_plugin_aligned_allocate_offset)(mm_plugin_heap_t heap, size_t n, size_t alignment, size_t offset, void ** out_ptr);
    status_t (*mm_plugin_deallocate)(mm_plugin_heap_t heap, void ** ptr, size_t size);
    status_t (*mm_plugin_deallocate_without_size)(mm_plugin_heap_t heap, void ** ptr);
    status_t (*mm_plugin_callocate)(mm_plugin_heap_t heap, size_t n, void ** out_ptr);
    status_t (*mm_plugin_reallocate)(mm_plugin_heap_t heap, void ** in_out_ptr, size_t size);
    status_t (*mm_plugin_usable_size)(mm_plugin_heap_t heap, void * ptr, size_t * out_size);
    status_t (*mm_plugin_bytes_remaining)(mm_plugin_heap_t heap, size_t * bytes_remaining);
    void     (*mm_plugin_debug)(mm_plugin_heap_t heap);
    status_t (*mm_plugin_inject_allocation)(mm_plugin_heap_t heap, void * ptr, size_t size);
    int (*mm_plugin_is_crash_consistent)(mm_plugin_heap_t heap);
    int (*mm_plugin_can_inject_allocation)(mm_plugin_heap_t heap);
  } mm_plugin_function_table_t;

#if defined(__cplusplus)
}
#endif

#if defined(__cplusplus)

#include <dlfcn.h>
#include <string>
#include <stdio.h>
#include <stdexcept>

#define LOAD_SYMBOL(X) _ft.X = reinterpret_cast<decltype(_ft.X)>(dlsym(_module, # X)); assert(_ft.X)

/**
 * C++ wrapper on C-based plugin API
 *
 */
class MM_plugin_wrapper
{
public:

  MM_plugin_wrapper(const std::string& plugin_path,
                    const std::string& config = "",
                    void * root_ptr = nullptr)
  {
    assert(plugin_path.empty() == false);
    _module = dlopen(plugin_path.c_str(), RTLD_NOW | RTLD_NODELETE); // RTLD_DEEPBIND |
    if(_module == nullptr) {
      char err[1024];
      sprintf(err, "%s\n", dlerror());
      throw std::invalid_argument(err);
    }

    LOAD_SYMBOL(mm_plugin_init);
    LOAD_SYMBOL(mm_plugin_create);
    LOAD_SYMBOL(mm_plugin_add_managed_region);
    LOAD_SYMBOL(mm_plugin_query_managed_region);
    LOAD_SYMBOL(mm_plugin_register_callback_request_memory);
    LOAD_SYMBOL(mm_plugin_allocate);
    LOAD_SYMBOL(mm_plugin_aligned_allocate);
    LOAD_SYMBOL(mm_plugin_aligned_allocate_offset);
    LOAD_SYMBOL(mm_plugin_deallocate);
    LOAD_SYMBOL(mm_plugin_deallocate_without_size);
    LOAD_SYMBOL(mm_plugin_callocate);
    LOAD_SYMBOL(mm_plugin_reallocate);
    LOAD_SYMBOL(mm_plugin_usable_size);
    LOAD_SYMBOL(mm_plugin_debug);
    LOAD_SYMBOL(mm_plugin_destroy);
    LOAD_SYMBOL(mm_plugin_inject_allocation);
    LOAD_SYMBOL(mm_plugin_bytes_remaining);
    LOAD_SYMBOL(mm_plugin_is_crash_consistent);
    LOAD_SYMBOL(mm_plugin_can_inject_allocation);

    //      dlclose(_module);

    /* create heap instance */
    _ft.mm_plugin_create(config.c_str(), root_ptr, &_heap);
  }

  MM_plugin_wrapper(const MM_plugin_wrapper &) = delete;

  MM_plugin_wrapper(MM_plugin_wrapper && other) noexcept
    : _module(std::move(other._module))
    , _ft(std::move(other._ft))
    , _heap(std::move(other._heap))
  {
    other._heap = nullptr;
  }

  virtual ~MM_plugin_wrapper() noexcept {
    if ( _heap )
    {
      _ft.mm_plugin_destroy(_heap);
    }
  }

  /* forwarding in liners */
  inline status_t init() noexcept {
    return _ft.mm_plugin_init();
  }

  inline status_t add_managed_region(void * region_base, size_t region_size) noexcept {
    return _ft.mm_plugin_add_managed_region(_heap, region_base, region_size);
  }

  inline status_t query_managed_region(unsigned region_id, void** out_region_base, size_t* out_region_size) noexcept {
    return _ft.mm_plugin_query_managed_region(_heap, region_id, out_region_base, out_region_size);
  }

  inline status_t register_callback_request_memory(request_memory_callback_t callback, void * param) noexcept {
    return _ft.mm_plugin_register_callback_request_memory(_heap, callback, param);
  }

  inline status_t allocate(size_t n, void ** out_ptr) noexcept {
    return _ft.mm_plugin_allocate(_heap, n, out_ptr);
  }

  inline status_t aligned_allocate(size_t n, size_t alignment, void ** out_ptr) noexcept {
    return _ft.mm_plugin_aligned_allocate(_heap, n, alignment, out_ptr);
  }

  inline status_t aligned_allocate_offset(size_t n, size_t alignment, size_t offset, void ** out_ptr) noexcept {
    return _ft.mm_plugin_aligned_allocate_offset(_heap, n, alignment, offset, out_ptr);
  }

  inline status_t deallocate(void ** ptr, size_t size) noexcept {
    return _ft.mm_plugin_deallocate(_heap, ptr, size);
  }

  inline status_t deallocate_without_size(void ** ptr) noexcept {
    return _ft.mm_plugin_deallocate_without_size(_heap, ptr);
  }

  inline status_t callocate(size_t n, void ** out_ptr) noexcept {
    return _ft.mm_plugin_callocate(_heap, n, out_ptr);
  }

  inline status_t reallocate(void ** in_out_ptr, size_t size) noexcept {
    return _ft.mm_plugin_reallocate(_heap, in_out_ptr, size);
  }

  inline status_t usable_size(void * ptr, size_t * out_size) noexcept {
    return _ft.mm_plugin_usable_size(_heap, ptr, out_size);
  }

  inline status_t bytes_remaining(size_t *bytes_remaining) const noexcept {
    return _ft.mm_plugin_bytes_remaining(_heap, bytes_remaining);
  }

  inline void debug(mm_plugin_heap_t heap) noexcept {
    return _ft.mm_plugin_debug(_heap);
  }

  inline status_t inject_allocation(void * ptr, size_t size) noexcept {
    return _ft.mm_plugin_inject_allocation(_heap, ptr, size);
  }

  inline int is_crash_consistent() noexcept {
    return _ft.mm_plugin_is_crash_consistent(_heap);
  }

  inline int can_inject_allocation() noexcept {
    return _ft.mm_plugin_can_inject_allocation(_heap);
  }

private:
  void *                     _module;
  mm_plugin_function_table_t _ft;
#if 0
  /* intention, but we avoid appealing to moveable_ptr */
  common::moveable_ptr<void> _heap;
#else
  void *                     _heap;
#endif
};

#include <limits>
#include <new>

/**
 * Standard C++ allocator wrapper
 *
 */
template <class T>
class MM_plugin_cxx_allocator
{
public:
  using value_type    = T;
  using pointer       = value_type*;
  using const_pointer = typename std::pointer_traits<pointer>::template
    rebind<value_type const>;
  using void_pointer       = typename std::pointer_traits<pointer>::template
    rebind<void>;
  using const_void_pointer = typename std::pointer_traits<pointer>::template
    rebind<const void>;

  using difference_type = typename std::pointer_traits<pointer>::difference_type;
  using size_type       = std::make_unsigned_t<difference_type>;

  template <class U> struct rebind {typedef MM_plugin_cxx_allocator<U> other;};

  MM_plugin_cxx_allocator(MM_plugin_wrapper& wrapper) noexcept : _wrapper(wrapper)
  {
  }

  template <class U>
  MM_plugin_cxx_allocator(MM_plugin_cxx_allocator<U> const& a_) noexcept : _wrapper(a_._wrapper)  {}

  pointer allocate(std::size_t n)
  {
    pointer p = nullptr;
    auto status = _wrapper.allocate(n*sizeof(value_type), reinterpret_cast<void**>(&p));
    if(status != 0) throw std::bad_alloc();
    return p;
  }

  void deallocate(pointer p, std::size_t n) noexcept
  {
    _wrapper.deallocate(reinterpret_cast<void**>(&p), n);
  }

  pointer allocate(std::size_t n, const_void_pointer)
  {
    return allocate(n);
  }

  template <class U, class ...Args>
  void construct(U* p, Args&& ...args)
  {
    ::new(p) U(std::forward<Args>(args)...);
  }

  template <class U>
  void destroy(U* p) noexcept
  {
    p->~U();
  }

  std::size_t max_size() const noexcept
  {
    return std::numeric_limits<size_type>::max();
  }

  MM_plugin_cxx_allocator select_on_container_copy_construction() const
  {
    return *this;
  }

  using propagate_on_container_copy_assignment = std::false_type;
  using propagate_on_container_move_assignment = std::false_type;
  using propagate_on_container_swap            = std::false_type;
  using is_always_equal                        = std::is_empty<MM_plugin_cxx_allocator>;

  MM_plugin_wrapper& _wrapper;
};

template <class T, class U>
bool
operator==(MM_plugin_cxx_allocator<T> const&, MM_plugin_cxx_allocator<U> const&) noexcept
{
  return true;
}

template <class T, class U>
bool
operator!=(MM_plugin_cxx_allocator<T> const& x, MM_plugin_cxx_allocator<U> const& y) noexcept
{
  return !(x == y);
}

#undef LOAD_SYMBOL
#endif

#pragma GCC diagnostic pop

#endif // __MM_PLUGIN_ITF_H__

max_stars_repo_name: omriarad/mcas

max_issues_repo_name: omriarad/mcas

max_forks_repo_name: omriarad/mcas


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"

void
cblas_ccopy (const int N, const void *X, const int incX, void *Y,
             const int incY)
{
#define BASE float
#include "source_copy_c.h"
#undef BASE
}

ext: c

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


text:
/**
 * @file iftCommon.h
 * @brief Common definitions and functions to remaining modules.
 */

#ifndef IFT_COMMON_H
#define IFT_COMMON_H

#ifdef __cplusplus
extern "C" {
#endif

#include <stdarg.h>
#include <assert.h>
#include <cblas.h>
#include <ctype.h>
#include <dirent.h>
#include <float.h>
#include <libgen.h>
#include <limits.h>
#include <math.h>
#include <mm_malloc.h>
#include <regex.h>
#include <stdarg.h>
#include <stdbool.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>

#ifdef __linux
#include <dirent.h>
#include <omp.h>
#endif

#if !defined(__APPLE__)
    #include <malloc.h>
#endif
#if defined(__linux)
#endif

#include "iftBasicDataTypes.h"
#include "iftMemory.h"
#include "iftDialog.h"
#include "iftString.h"

/**
 * Common definitions
 */

#ifndef PI
#define PI IFT_PI
#endif

typedef enum ift_axis_order {
  XYZ, XZY, YXZ, YZX, ZXY, ZYX
} iftAxisOrder;

#define IFT_RANDOM_SEED (unsigned int) 213344
#define IFT_MAXWEIGHT     4095.0
#define IFT_AXIS_X 0
#define IFT_AXIS_Y 1
#define IFT_AXIS_Z 2
#define IFT_PI      3.14159265358979323846
#define IFT_INTERIOR    0
#define IFT_EXTERIOR    1
#define IFT_BOTH        2
#define IFT_WHITE       0
#define IFT_GRAY        1
#define IFT_BLACK       2
#define IFT_NIL        -1
#define IFT_INCREASING  1
#define IFT_DECREASING  0
#define IFT_EPSILON     1E-07

#ifndef FALSE
#define FALSE 0
#endif
#ifndef TRUE
#define TRUE 1
#endif

/**
 * Common operations
 */

#ifndef iftMax
#define iftMax(x,y) (((x) > (y))?(x):(y))
#endif

#ifndef iftMin
#define iftMin(x,y) (((x) < (y))?(x):(y))
#endif

#define iftRound(x) ((x < 0)?(int)(x-0.5):(int)(x+0.5))
#define iftSign(x) ((x >= 0)?1:-1)

/** @brief Euclidean Distance between voxels or points */
#define iftPointDistance(u,v) (sqrtf((u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y)+(u.z-v.z)*(u.z-v.z)))

/** @brief Euclidean Distance between voxels or porints */
#define iftVoxelDistance(u,v) iftPointDistance(u,v)

/** @brief Computes the Squared Distance between two voxels or points */
#define iftSquaredVoxelDistance(u,v) ((u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y)+(u.z-v.z)*(u.z-v.z))

/* @brief Computes the Smooth Euclidean Distances from a pre-computed Squared Euclidean distance among 26-neighbors
 * @param squared_dist The pre-computed Squared Euclidean Distance among 26-neighbors. The details are in N. Kiryati and G. Sze kely, "Estimating shortest paths and minimal distances on digitized three-dimensional surfaces," Pattern Recognition, vol. 26, pp. 1623-1637, 1993.
 * @return The resulting Smooth Euclidean Distance
 */
#define iftSmoothEuclideanDistance(squared_dist) ((fabs(squared_dist-1)<1.0e-6)? 0.9016 : (fabs(squared_dist-2)<1.0e-6)? 1.289 : 1.615 )

/**
 * @brief Converts a Point to a Voxel
 * @author Samuel
 * @date May 24, 2016
 */
#define iftPointToVoxel(point) ((iftVoxel) {(int) (point).x, (int) (point).y, (int) (point).z})

/**
 * @brief Converts a Voxel to a Point
 * @author Samuel
 * @date May 24, 2016
 */
#define iftVoxelToPoint(v) ((iftPoint) {(v).x, (v).y, (v).z})

#define iftImageCenter(img) ((iftVoxel){((img)->xsize - 1)/2, ((img)->ysize - 1)/2, ((img)->zsize - 1)/2})

/** @brief Constant used to align memory */
#define IFT_MEMORY_ALIGNMENT 16

/**
 * @brief Init the GPU devices.
 * @warning If not called, the GPU devices will be started in the first call for a GPU function.
 * @author Peixinho
 */
void iftStartGPU();

/**
 * @brief Stop the GPU devices.
 * @author Peixinho
 */
void iftStopGPU();

/**
 * @brief Checks if the long positive number is prime or not.
 * @author Samuel Martins
 * @date Jan 17th, 2016
 * @note Based on https://en.wikipedia.org/wiki/Primality_test
 */
bool iftIsPrime(long n);

/**
 * @brief Computes a fast power for small integer exponents.
 * @param base
 * @param b
 * @return Returns base^b
 * @author Peixinho
 * @date May, 2016
 * @warning The funcion is faster upto b=10, after that the regular function pow is called.
 */
double iftFastPow(double base, unsigned b);

/**
 * @brief Prints a Float array on stdout of lentgh <n>.
 */
void      iftPrintFloatArray(float* v, int n);

/**
 * @brief Prints a Double array on stdout of lentgh <n>.
 */
void      iftPrintDoubleArray(double* v, int n);

/**
 * @brief Prints an Integer array on stdout of lentgh <n>.
 */
void      iftPrintIntArray(int* v, int n);

/**
 * @brief Returns a random integer number between low and high.
 */
int iftRandomInteger (int low, int high);

/**
 * @brief Randomly selects nelems of the set [low, high]
 */
int *iftRandomIntegers (int low, int high, int nelems);

/**
 * @brief Randomly selects a normal distributed (N(0,1)) float number
 */
float iftRandomNormalFloat(void);

/**
 * @brief Returns initial time
 */
timer *iftTic(void);

/**
 * @brief Returns final time
 */
timer *iftToc(void);

/**
 * @brief Prints the computational time from tstart to tend, as obtained by iftTic() and iftToc() functions.
 */
void iftPrintCompTime(timer *tstart, timer *tend, const char *msg, ...);

/**
 * @brief Computes the difference in ms from the initial time to the final time.
 *
 * @param tic The initial time.
 * @param toc The final time.
 *
 * @return The computed time difference in milliseconds.
 *
 * @warning The memory allocated for <tic> and <toc> is freed inside this function.
 */
float iftCompTime(timer *tic, timer *toc);

/**
 * @brief Returns a string with the formatted time: days hours mins secs ms.
 *
 * @author Samuel
 *
 * Given a time in ms, this function returns a formatted time: days hours mins secs ms.
 * The runtime in MILISECONDS can be obtained with the function iftCompTime(timer *tic, timer *toc).
 *
 * @param runtime Time in ms.
 * @return The formatted time.
 */
char *iftFormattedTime(float runtime);

/**
 * @brief Prints the runtime in the following format: days hours mins secs ms.
 *
 * @author Samuel Martins
 * @date October 23, 2015
 * @sa iftFormattedTime()
 *
 * Given a time in ms, this function prints the runtime following the format: days hours mins secs ms.
 * The runtime in MILISECONDS can be obtained with the function iftCompTime(timer *tic, timer *toc).
 *
 * @param runtime Time in ms.
 */
void iftPrintFormattedTime(float runtime);

/**
 * @brief Writes a timer to a given file, using the following format %s: %f, where %s is the
 * given information corresponding to the time and %f is the current time in milliseconds.
 *
 * @param tic is freed by the function.
 */
void iftWriteTimerToFile(const char *filename, const char *information, timer *tic);

/**
 * @brief Generates seed for rand(), used in iftRandomInteger.
 */
void iftRandomSeed(unsigned int);

/**
 * @brief Returns the factorial of a number or NIL in case of overflow
 */
long double iftFactorial(int n);

/**
 * @brief Returns the limit to avoid overflow in factorial computation
 */
int iftFactorialLimit(void);

void iftUnitNorm(float *feats, int nelems);
void iftNormalizeFeatures(float *feats, int nelems);

float iftSquaredFeatDistance(float *A, float *B, int n);
float iftFeatDistance(float *A, float *B, int n);

/**
 * @brief Computes the Manhattan Distance between two arrays of float considering the sinal, ie.,
 * wihout the module between the feature difference: dist = (b[0]-a[0]) + (b[1]-a[1]) + ...
 */
float iftSignedManhattanDistance(float *A, float *B, int n);

/**
 * @brief Computes the log(val) in the specified base.
 * @author Peixinho
 */
double iftLog(double val, double base);

/**
 * @brief Returns the Central Voxel from the Bounding Box
 * @author Samuel Martins
 * @date Mar 15, 2016
 * @ingroup Geometry
 */
iftVoxel iftCenterFromBoundingBox(iftBoundingBox bb);

/**
 * @brief Vector operations defined as macros so that iftVoxels, iftVectors, and iftPoints may be used interchangeably.
 */
/**
 * @brief Subtracts point vec1 from vec2 (it works for iftVectors/iftVoxels/iftPoints).
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorSub(vec1, vec2) {(vec1).x - (vec2).x, (vec1).y - (vec2).y, (vec1).z - (vec2).z}
/**
 * @brief Adds points vec1 and vec2 (it works for iftVectors/iftVoxels/iftPoints).
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorSum(vec1, vec2) {(vec1).x + (vec2).x, (vec1).y + (vec2).y, (vec1).z + (vec2).z}
/**
 * @brief Computes the cross product between two voxels, points, or vectors.
 *
 * @author Thiago Vallin Spina (changed the macro function's name)
 * @date Mar 04, 2016
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorCrossProd(a, b) {(a).y*(b).z - (a).z*(b).y, (a).z*(b).x - (a).x*(b).z, (a).x*(b).y - (a).y*(b).x}
/**
 * @brief Computes the inner product between two voxels, points, or vectors.
 *
 * @author Thiago Vallin Spina (changed the macro function's name)
 * @date Mar 04, 2016
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorInnerProduct(a, b) (((a).x*(b).x + (a).y*(b).y + (a).z*(b).z))

/**
 * @brief Computes the Vector Magnitude from a voxel or point
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 */
#define iftVectorMagnitude(v) (sqrtf(iftVectorInnerProduct((v),(v))))

/**
 * @brief Returns the XY angle of a 2D vector
 * @author Peixinho
 * @date April, 2016
 * @ingroup Geometry
 */
#define iftVectorXYAngle(v) ( (atan2(v.y, v.x)+IFT_PI) * (180.0/IFT_PI) )

/**
 * @brief Rounds a vector to integer coordinates.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorRound(vec1) {iftRound((vec1).x), iftRound((vec1).y), iftRound((vec1).z)}
/**
 * @brief Multiplies a vector by an scalar.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorScalarProd(vec, s) {(vec).x*(s), (vec).y*(s), (vec).z*(s)}
/**
 * @brief Adds an scalar to a vector.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
#define iftVectorScalarSum(vec, s) {(vec).x+(s), (vec).y+(s), (vec).z+(s)}
/**
 * @brief Divides a vector by an scalar.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @return The vector divided by the scalar, or the vector itself if the scalar is close to 0.0
 *
 * @warning The result must be cast as an iftVoxel, iftVector, or iftPoint.
 */
iftVector iftVectorScalarDiv(iftVector vec, double s);
/**
 * @brief Verifies if a vector has coordinates (0, 0, 0).
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
#define iftIsNullVector(vec) (iftAlmostZero((vec).x) && iftAlmostZero((vec).y) && iftAlmostZero((vec).z))
/**
 * @brief Verifies if three points, voxels, or vectors are collinear.
 *
 * This is achieved by testing if the cross product between the vectors P1->P2 and P1->P3 is Null (0,0,0).
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
#define iftCollinearPoints(P1, P2, P3) (iftIsNullVector((iftVector)iftVectorCrossProd((iftVector)iftVectorSub((P2), (P1)), (iftVector)iftVectorSub((P3), (P1)))))

/**
 * @brief Tests if two iftVoxels are equal
 *
 * @author Thiago Vallin Spina
 * @date Apr 19, 2016
 *
 * @param u1 The first voxel
 * @param u2 The second voxel
 * @return True if their coordinates are the same
 */
bool iftVoxelsAreEqual(iftVoxel u1, iftVoxel u2);

/**
 * @brief Computes the mean voxel
 *
 * @author Thiago Vallin Spina
 * @date May 1, 2016
 *
 * @param u1 The first voxel
 * @param u2 The second voxel
 * @return The mean voxel (u1 + u2) / 2
 */
iftVoxel iftMeanVoxel(iftVoxel u1, iftVoxel u2);

/**
 * @brief Tests if two iftPoints/iftVectors are equal
 *
 * @author Thiago Vallin Spina
 * @date Apr 19, 2016
 *
 * @param u1 The first point/vector
 * @param u2 The second point/vector
 * @return True if their coordinates are (almost) the same
 */
#define iftVectorsAreEqual(u1, u2) (iftAlmostZero((u1).x-(u2).x) && iftAlmostZero((u1).y-(u2).y) && \
                                      iftAlmostZero((u1).z-(u2).z))

/**
 * @brief Normalizes a vector if its norm is greater than 0.0, otherwise returns the vector itself.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
iftVector iftNormalizeVector(iftVector v);

/**
 * @brief Projects vector U onto vector V.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 *
 * @warning If vector V has norm close to 0.0 the function issues and iftError.
 */
iftVector iftProjectVector(iftVector U, iftVector V);

/**
 * @brief Returns the distance between P0 and the line from P1 to P2, whose size is P1P2
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
double iftVoxelLineDist2D(iftVoxel P0, iftVoxel P1, iftVoxel P2, double P1P2);

/**
 * @brief Given the vector from P1 to P0, returns the voxel projected on the line between P1 and P2.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
iftVoxel iftClosestVoxelOnLine2D(iftVoxel P0, iftVoxel P1, iftVoxel P2);

/**
 * @brief Returns the position of P0 with respect to the line from P1 to P2.
 *
 * Negative values indicate left side, 0 indicates on the line,
 * and positive values indicate right side.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 * @ingroup Geometry
 */
int iftVoxelLinePosition2D(iftVoxel P0, iftVoxel P1, iftVoxel P2);

/**
 * @brief Computes the circle's area given a radius.
 *
 * @author Thiago Vallin Spina
 * @date Mar 04, 2016
 *
 * @ingroup Geometry
 */
static inline double iftCircleArea(double r) {
    return (r*r)*IFT_PI;
}

/* @brief Gets rid of the carriage return and the line feed characteres introduced by DOS systems when reading
 * strings from ASCII files.
 *
 * @author Alexandre Xavier Falcao
 *
 * @param line The string in which all carriage return characters are removed.
 */
void iftRemoveCarriageReturn(char *line);

/**
 * @brief Evaluates the sigmoid function, with x = value.
 *
 * @param alfa Controls the decay of the function.
 * @author Renzo Phellan
 */
float iftSigmoidalValue(float value, float alfa);

void iftVerifyToken(FILE *fp, char *token, char *function);
void iftReadIntValue(FILE *fp, int *value, char *token, char *function);
void iftReadIntValues(FILE *fp, int **value, int nvalues, char *token, char *function);
void iftWriteIntValue(FILE *fp, int value, char *token);
void iftWriteIntValues(FILE *fp, int *value, int nvalues, char *token);
void iftReadFloatValue(FILE *fp, float *value, char *token, char *function);
void iftReadFloatValues(FILE *fp, float **value, int nvalues, char *token, char *function);
void iftWriteFloatValue(FILE *fp, float value, char *token);
void iftWriteFloatValues(FILE *fp, float *value, int nvalues, char *token);
void iftReadDoubleValue(FILE *fp, double *value, char *token, char *function);
void iftReadDoubleValues(FILE *fp, double **value, int nvalues, char *token, char *function);
void iftWriteDoubleValue(FILE *fp, double value, char *token);
void iftWriteDoubleValues(FILE *fp, double *value, int nvalues, char *token);
void iftSkipComments(FILE *fp);

/* These functions are currently used to communicate with numpy */
void iftWriteRawIntArray(char *filename, int *array, int n);
int* iftReadRawIntArray(char *filename, int n);

/**
 * @brief Computes the mean value of a float array.
 */
float iftMean(float *x, int n);
/**
 * @brief Computes the variance of a float array.
 */
float iftVar(float *x, int n);
/**
 * @brief Computes the covariance of float arrays.
 */
float iftCov(float *x, float *y, int n);

/**
 * @brief Tests if a real (float/double) number is zero. It should be used for comparisons if two floats are equal.
 *
 * @author Thiago Vallin Spina
 */
int iftAlmostZero(double x);
/**
 * @brief Computes the mod operation of <a> by <n>. If the mod value is negative, it is placed within [0,n-1].
 */
int iftSafeMod(int a, int n);

/**
 * @brief Gets an Integer Array with the Unique Elements from an Input Array.
 *
 * Ex: Given the array = [1, 2, 3, 4, 4, 3, 2, 1], with size n = 8, the resulting array of unique elements will be:
 * [1, 2, 3, 4]
 *
 * @param  array Input Int Array to be scanned.
 * @param  n     Number of Elements from the Input Int Array.
 * @return       An Integer Array with the Unique Elements from <b>array</b>.
 */
iftIntArray *iftIntArrayUnique(const int *array, int n);

/**
 * @brief Counts the Number of Unique Elements from an array.
 *
 * @param array Integer Array whose Number of Unique Elements will be counted.
 * @param n     Array Size.
 * @return      The number of Unique Elements.
 */
int iftCountUniqueIntElems(const int *array, int n);

/**
 * @brief This function finds the most suitable normalization value to be used during image color/feature conversion/normalization.
 *
 * It is afunction of the most common number of bits in imaging sensors. It assumes that these values are in {1, 8, 10,
 * 12, 16, 24, 32} bits per voxel.
 *
 * @author Alexandre Xavier Falcao.
 *
 * @param maxval The maximum value of an image, for instance.
 */
int iftNormalizationValue(int maxval);

/**
 * @brief Converts a number in radians to degrees.
 *
 * @author Thiago Vallin Spina
 */
static inline double iftDegrees(double rad) {
  return rad * 180.0 / IFT_PI;
}

/**
 * @brief Converts a number in degrees to radians.
 *
 * @author Thiago Vallin Spina
 */
static inline double iftRadians(double deg) {
  return deg * IFT_PI / 180.0;
}

/**
 * @brief Converts a negative radian number to positive.
 *
 * @author Thiago Vallin Spina
 */
static inline double iftNegativeToPositiveRadians(double rad) {
  return (rad < 0) ? 2 * IFT_PI + rad : rad;
}

/**
 * @brief Converts a negative degree number to positive.
 *
 * @author Thiago Vallin Spina
 */
static inline double iftNegativeToPositiveDegrees(double deg) {
  return (deg < 0) ? 360.0+deg : deg;
}

/**
 * @brief Compute the index of the element in vector x with minimum value
 *
 * @param x Float vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftArgmin(const int *x, int n);

/**
 * @brief Compute the index of the element in vector x with maximum value
 *
 * @param x Float vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftArgmax(const int *x, int n);

/**
 * @brief Compute the index of the element in vector x with minimum value
 *
 * @param x Float vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftFArgmin(const float *x, int n);

/**
 * @brief Compute the index of the element in vector x with maximum value
 *
 * @param x Float vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftFArgmax(const float *x, int n);

/**
 * @brief Compute the index of the element in vector x with minimum value
 *
 * @param x Double vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftDArgmin(const double *x, int n);

/**
 * @brief Compute the index of the element in vector x with maximum value
 *
 * @param x Double vector.
 * @param n Size of the vectors.
 * @return index.
*/
int iftDArgmax(const double *x, int n);

/**
 * @brief Finds the index of <elem> in array <x> and returns it.
 *
 * @author Peixinho
 *
 * @param x Unsorted array of elements.
 * @param n The size of the array.
 * @param elem The elem to be found.
 *
 * @return The element's index in <x> or NIL if not found
 */
int iftFindIntArrayElementIndex(int *x, int n, int elem);

/**
 * @brief Compute the Gaussian Probability Density Function for a set of values.
 * @author Samuel Martins
 * @date May 12, 2016
 *
 * @param  vals   Array of Values to compute the Gaussian PDF.
 * @param  n_vals Number of values.
 * @param  mean   Mean used to compute the Gaussian PDF.
 * @param  stdev  Standard Deviation used to compute the Gaussian PDF.
 * @return        Array of probabilities, one for each input value.
 */
double *iftGaussPDF(double *vals, size_t n_vals, double mean, double stdev);

/**
 * @brief Compute the Cumulative Density of the (Gaussian) PDF (Normal Distribution) in a Random Variable with value x.
 * @author Samuel Martins
 * @date May 12, 2016
 *
 * It gives the Area under the Probability Density Function (Normal Distribution) from minus infinity to x. \n
 * Then, it also describes the probability that a random variable X takes on a value less than or
 * equal to a number x. \n
 *
 * @param  x      Value from the Random Variable to compute the CDF.
 * @param  mean   Mean used to compute the Gaussian PDF.
 * @param  stdev  Standard Deviation used to compute the Gaussian PDF.
 * @return        The area under the PDF from minus to x.
 */
double iftCDF(double x, double mean, double stdev);

#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: Fabio-Kubo/processamento-imagens

max_issues_repo_name: Fabio-Kubo/processamento-imagens

max_forks_repo_name: Fabio-Kubo/processamento-imagens


text:
#ifndef IFT_COMMON_H_
#define IFT_COMMON_H_

#include <omp.h>
#include <stdlib.h>
#include <stdio.h>
#if !defined(__APPLE__)
	#include <malloc.h>
#endif
#include <math.h>
#include <string.h>
#include <limits.h>
#include <float.h>
#include <sys/time.h>
#include <time.h>
#include <unistd.h>
#include <assert.h>
#include <cblas.h>
#include <dirent.h>
#include <regex.h>

/*
 * Common data types
 */

#define INFINITY_INT  INT_MAX
#define INFINITY_FLT  FLT_MAX
#define INFINITY_DBL  DBL_MAX
#define INFINITY_LDBL LDBL_MAX

typedef struct timeval timer;

typedef unsigned char  uchar;
typedef unsigned short ushort;
typedef unsigned int   uint;
typedef unsigned long long ullong;

typedef struct ift_band {
  float *val;
} iftBand;

typedef struct ift_vector {
  float x,y,z;
} iftVector, iftPoint;

typedef struct ift_voxel {
  int x,y,z;
} iftVoxel;

typedef struct ift_dcomplex
{
  double r;
  double i;
} iftComplex;

typedef struct file_list {
  char **filesRoutes;
  char **filesNames;
  int    n;
} fileList;

typedef struct ift_name_metric_pair {
  char *name;
  float metric;
} iftNameMetricPair;

/**
 * Common definitions
 */

#define IFT_RANDOM_SEED (unsigned int) 213344
#define MAXWEIGHT     4095.0
#define AXIS_X 0
#define AXIS_Y 1
#define AXIS_Z 2
#define PI          3.1415926536
#define INTERIOR    0
#define EXTERIOR    1
#define BOTH        2
#define WHITE       0
#define GRAY        1
#define BLACK       2
#define NIL        -1
#define INCREASING  1
#define DECREASING  0
#define Epsilon     1E-05

/**
 * Common operations
 */

#ifndef MAX
#define MAX(x,y) (((x) > (y))?(x):(y))
#endif

#ifndef MIN
#define MIN(x,y) (((x) < (y))?(x):(y))
#endif

#define ROUND(x) ((x < 0)?(int)(x-0.5):(int)(x+0.5))
#define SIGN(x) ((x >= 0)?1:-1)

/**
 * Common functions to allocate memory
 */

char     *iftAllocCharArray(int n);
uchar    *iftAllocUCharArray(int n);
short    *iftAllocShortArray(int n);
ushort   *iftAllocUShortArray(int n);
uint     *iftAllocUIntArray(int n);
ullong   *iftAllocULLongArray(int n);
int      *iftAllocIntArray(int n);
float    *iftAllocFloatArray(int n);
double   *iftAllocDoubleArray(int n);
iftComplex  *iftAllocComplexArray(int n);
long double *iftAllocLongDoubleArray(int n);

void      iftPrintFloatArray(float* v, int n);
/**
 * Error messages
 */

#define MSG1  "Cannot allocate memory space"
#define MSG2  "Cannot open file"

/**
 *  Error message msg is printed in function func and the program
 *  exits abnormally.
 */

void iftError(char *msg,char *func);

/**
 *  Warning message msg is printed in function func and the program
 *  continues.
 */

void iftWarning(char *msg,char *func);

/**
 *  The contents of a and b are interchanged.
 */

void iftSwitchValues(int *a, int *b);
void iftSSwitchValues(char *a, char *b, int size);
void iftFSwitchValues(float *a, float *b);
void iftDSwitchValues(double *a, double *b);
void iftSwitchVoxels(iftVoxel *u, iftVoxel *v);

/**
 * Returns a random integer number between low and high.
 */

int iftRandomInteger (int low, int high);

/*
 * Randomly selects nelems of the set [low, high]
 */
int *iftRandomIntegers (int low, int high, int nelems);

/**
 * Randomly selects a normal distributed (N(0,1)) float number
 */

float iftRandomNormalFloat(void);

/**
 * Returns the distance between P0 and the line from P1 to P2, whose
 * size is P1P2
 */

float iftVoxelLineDist2D(iftVoxel P0, iftVoxel P1, iftVoxel P2, float P1P2);

/**
 * Returns the position of P0 with respect to the line from P1 to
 * P2. Negative values indicate left side, 0 indicates on the line,
 * and positive values indicate right side.
 */

int iftVoxelLinePosition2D (iftVoxel P0, iftVoxel P1, iftVoxel P2);

/**
 * Returns initial time
 */

timer *iftTic(void);

/**
 * Returns final time
 */

timer *iftToc(void);

/**
 * Computes the difference in ms from the initial time to the final time
 */

float iftCompTime(timer *tic, timer *toc);

/**
 * Generates seed for rand(), used in iftRandomInteger.
 */

void iftRandomSeed(unsigned int);

/**
 * Returns the factorial of a number or NIL in case of overflow
 */

long double iftFactorial(int n);

/**
 * Returns the limit to avoid overflow in factorial computation
 */

int iftFactorialLimit(void);

float iftVoxelDistance(iftVoxel u, iftVoxel v);
int iftSquaredVoxelDistance(iftVoxel u, iftVoxel v);

float iftPointDistance(iftPoint u, iftPoint v);

int   iftVoxelSquareDistance(iftVoxel u, iftVoxel v);

float iftInnerProduct(iftVector a, iftVector b);

iftVector iftCrossProduct(iftVector a, iftVector b);

char iftCollinearPoints(iftPoint P1, iftPoint P2, iftPoint P3);

char iftCollinearVoxels(iftVoxel P1, iftVoxel P2, iftVoxel P3);

iftVector iftNormalizeVector(iftVector v);

void iftNormalizeFloatArray(float *array, int nelems);

float iftVectorMagnitude(iftVector v);

void iftRemoveCarriageReturn(char *token); /* useful to get rid of the
					      carriage return and the
					      line feed characteres
					      introduced by DOS
					      systems when reading
					      strings from ASCII
					      files */

void   iftWriteFloatArray(float *v, int size, char *filename);
float *iftReadFloatArray(char *filename, int *size);

/**
 * Evaluates the sigmoid function, with x = value.
 * Alfa controls the decay of the function.
 */

float iftSigmoidalValue(float value, float alfa);

void iftVerifyToken(FILE *fp, char *token, char *function);
void iftReadIntValue(FILE *fp, int *value, char *token, char *function);
void iftReadIntValues(FILE *fp, int **value, int nvalues, char *token, char *function);
void iftWriteIntValue(FILE *fp, int value, char *token);
void iftWriteIntValues(FILE *fp, int *value, int nvalues, char *token);
void iftReadFloatValue(FILE *fp, float *value, char *token, char *function);
void iftReadFloatValues(FILE *fp, float **value, int nvalues, char *token, char *function);
void iftWriteFloatValue(FILE *fp, float value, char *token);
void iftWriteFloatValues(FILE *fp, float *value, int nvalues, char *token);
void iftReadDoubleValue(FILE *fp, double *value, char *token, char *function);
void iftReadDoubleValues(FILE *fp, double **value, int nvalues, char *token, char *function);
void iftWriteDoubleValue(FILE *fp, double value, char *token);
void iftWriteDoubleValues(FILE *fp, double *value, int nvalues, char *token);
void iftSkipComments(FILE *fp);
char iftVoxelsAreEqual(iftVoxel u1, iftVoxel u2);
char iftPointsAreEqual(iftPoint u1, iftPoint u2);

/**
 * Common function to handle arrays
 */

void iftCopyFloatArray(float *array1, float *array2, int nelems);
void iftCopyIntArray(int *array1, int *array2, int nelems);
int *iftMergeIntArray(int *array1, int n1, int *array2, int n2, int *nelems);
int *iftIntArrayOfUniqueElemsTransform(int *array, int *n);
fileList             *iftCreateFileList(void);
void                 iftDestroyFileList(fileList **list);
fileList             *iftGetFiles(char *dirname, char *type);

/* These functions are currently used to communicate with numpy */
void iftWriteRawIntArray(char *filename, int *array, int n);
int* iftReadRawIntArray(char *filename, int n);

float iftMean(float *x, int n);
float iftVar(float *x, int n);
float iftCov(float *x, float *y, int n);

int iftAlmostZero(float x);

#endif


text:
#pragma once

#include <nextalign/nextalign.h>

#include <gsl/string_span>
#include <string>

#include "../nextalign_private.h"
#include "../utils/to_underlying.h"

using NucleotideSequenceSpan = SequenceSpan<Nucleotide>;

Nucleotide toNucleotide(char nuc);

char nucToChar(Nucleotide nuc);

inline std::ostream& operator<<(std::ostream& os, const Nucleotide& nucleotide) {
  os << std::string{to_underlying(nucleotide)};
  return os;
}

max_stars_repo_name: davidcroll/nextclade

max_issues_repo_name: davidcroll/nextclade

max_forks_repo_name: davidcroll/nextclade


text:
#ifndef QR_H
#define QR_H

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>

/**
 * Computes, stores, and maintains a QR decomposition of a gsl_matrix.
 *
 * A matrix-pointer is passed to the QR on construction. This pointer is kept through
 * the life-time of the QR object (but ownership is not assumed). If the contents
 * of the matrix is changed, the QR can be update to reflect the change using the
 * abstract UpdateFromMatrix function. This is not automatically done on construction.
 */
class QR{
 protected:
  const int m, n; ///< Dimensions of matrix

  /** Decomposes M into U*S*V^t.  */
  QR(gsl_matrix* M);
 public:

  /**
   * Constructs a QR object. If MKL is available it will be an MKLQR and if GSL
   * is available it will be GSLQR
   */
  static QR* createQR(gsl_matrix* M);

  virtual ~QR();

  /** Update Q and R to reflect the decomposition of the matrix. */
  virtual void updateFromMatrix() = 0;

  /** Print the QR to standard out */
  void print() const;

  gsl_matrix* getMatrix() const;
  gsl_matrix* getQ() const;
  gsl_matrix* getR() const;
 protected:
  gsl_matrix * const m_matrix;
  gsl_matrix * const m_Q;
  gsl_matrix * const m_R;

};

#endif