text:
/////////////////////////////////////////////////////////////////////
// = NMatrix
//
// A linear algebra library for scientific computation in Ruby.
// NMatrix is part of SciRuby.
//
// NMatrix was originally inspired by and derived from NArray, by
// Masahiro Tanaka: http://narray.rubyforge.org
//
// == Copyright Information
//
// SciRuby is Copyright (c) 2010 - 2014, Ruby Science Foundation
// NMatrix is Copyright (c) 2012 - 2014, John Woods and the Ruby Science Foundation
//
// Please see LICENSE.txt for additional copyright notices.
//
// == Contributing
//
// By contributing source code to SciRuby, you agree to be bound by
// our Contributor Agreement:
//
// * https://github.com/SciRuby/sciruby/wiki/Contributor-Agreement
//
// == getrs.h
//
// getrs function in native C++.
//

/*
 *             Automatically Tuned Linear Algebra Software v3.8.4
 *                    (C) Copyright 1999 R. Clint Whaley
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *   1. Redistributions of source code must retain the above copyright
 *      notice, this list of conditions and the following disclaimer.
 *   2. Redistributions in binary form must reproduce the above copyright
 *      notice, this list of conditions, and the following disclaimer in the
 *      documentation and/or other materials provided with the distribution.
 *   3. The name of the ATLAS group or the names of its contributers may
 *      not be used to endorse or promote products derived from this
 *      software without specific written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE ATLAS GROUP OR ITS CONTRIBUTORS
 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef GETRS_H
#define GETRS_H

extern "C" {
#if defined HAVE_CBLAS_H
  #include <cblas.h>
#elif defined HAVE_ATLAS_CBLAS_H
  #include <atlas/cblas.h>
#endif
}

namespace nm { namespace math {

/*
 * Solves a system of linear equations A*X = B with a general NxN matrix A using the LU factorization computed by GETRF.
 *
 * From ATLAS 3.8.0.
 */
template <typename DType>
int getrs(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE Trans, const int N, const int NRHS, const DType* A,
           const int lda, const int* ipiv, DType* B, const int ldb)
{
  // enum CBLAS_DIAG Lunit, Uunit; // These aren't used. Not sure why they're declared in ATLAS' src.

  if (!N || !NRHS) return 0;

  const DType ONE = 1;

  if (Order == CblasColMajor) {
    if (Trans == CblasNoTrans) {
      nm::math::laswp<DType>(NRHS, B, ldb, 0, N, ipiv, 1);
      nm::math::trsm<DType>(Order, CblasLeft, CblasLower, CblasNoTrans, CblasUnit, N, NRHS, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
    } else {
      nm::math::trsm<DType>(Order, CblasLeft, CblasUpper, Trans, CblasNonUnit, N, NRHS, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasLeft, CblasLower, Trans, CblasUnit, N, NRHS, ONE, A, lda, B, ldb);
      nm::math::laswp<DType>(NRHS, B, ldb, 0, N, ipiv, -1);
    }
  } else {
    if (Trans == CblasNoTrans) {
      nm::math::trsm<DType>(Order, CblasRight, CblasLower, CblasTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasRight, CblasUpper, CblasTrans, CblasUnit, NRHS, N, ONE, A, lda, B, ldb);
      nm::math::laswp<DType>(NRHS, B, ldb, 0, N, ipiv, -1);
    } else {
      nm::math::laswp<DType>(NRHS, B, ldb, 0, N, ipiv, 1);
      nm::math::trsm<DType>(Order, CblasRight, CblasUpper, CblasNoTrans, CblasUnit, NRHS, N, ONE, A, lda, B, ldb);
      nm::math::trsm<DType>(Order, CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, NRHS, N, ONE, A, lda, B, ldb);
    }
  }
  return 0;
}

/*
* Function signature conversion for calling LAPACK's getrs functions as directly as possible.
*
* For documentation: http://www.netlib.org/lapack/double/dgetrs.f
*
* This function should normally go in math.cpp, but we need it to be available to nmatrix.cpp.
*/
template <typename DType>
inline int clapack_getrs(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE trans, const int n, const int nrhs,
                         const void* a, const int lda, const int* ipiv, void* b, const int ldb) {
  return getrs<DType>(order, trans, n, nrhs, reinterpret_cast<const DType*>(a), lda, ipiv, reinterpret_cast<DType*>(b), ldb);
}

} } // end nm::math

#endif // GETRS_H

max_stars_repo_name: blackwinter-attic/nmatrix

max_issues_repo_name: blackwinter-attic/nmatrix

max_forks_repo_name: blackwinter-attic/nmatrix


text:
#ifndef FILE_PROCESSINGUNIT
#define FILE_PROCESSINGUNIT

#include "../devicedata/memorymanager.h"
#include "timer.h"

#ifndef CPU_ONLY
    #include "gpu_handle.h"

    #ifdef MAGMA
        #include "GpuMagma_queue.h"
        #include "magma_v2.h"
    #endif
#else
    #pragma message("Cpu only mode is activated. Gpu suppport is disabled. If you wish to activate Gpu support, remove compiler flag CPU_ONLY")
#endif

#ifdef __INTEL_CLANG_COMPILER
  #pragma message("Intel Math Kernel Library is used for acceleration of linear algebra operations.")
  #include <mkl.h>
#else
  #pragma message("CBLAS/LAPACK/LAPACKE are used for acceleration of linear algebra operations.")
  #include <cblas.h>
  #include <lapacke.h>
#endif

#include <map>
#include <type_traits>
#include <iostream>
#include <memory>
#include <string>

template<class floating>
class ProcessingUnitDevice;

template<class floating>
using ProcessingUnit = std::shared_ptr<ProcessingUnitDevice<floating>>;

/**
 * Class OperationType is used for choosing the matrix operation
 * type in several operations of interface ProcessingUnitDevice.
 */
enum class OperationType {
    Identical, //!< taking the matrix without any changes
    Transposed, //!< taking the transposition of the matrix
    Hermitian //!< taking the complex conjugate of the matrix
};

/**
 * Class ProcessingUnitDevice is an abstract interface for easy calls to
 * BLAS/LAPACK/LAPACKE-like function calls on a certain device,
 * e.g. for CUDA/cuSOLVER on a Gpu, or for MAGMA on a Gpu.
 *
 * The floating point type is passed as a template parameter.
 *
 * @tparam floating floating point type
 */
template<class floating>
class ProcessingUnitDevice {
public:
    /**
     * Constructor.
     */
    ProcessingUnitDevice();

    /**
     * Destructor.
     */
    virtual ~ProcessingUnitDevice() = default;

    /**
     * Creates a new instance of a Timer stopwatch on the current device.
     * @return new Timer instance
     */
    virtual Timer create_timer() const = 0;

    /**
     * Returns the MemoryManager associated with the current device.
     * @return MemoryManager
     */
    virtual MemoryManager get_memory_manager() const = 0;

    /**
     * Returns the subclass's name in a human-readable format as a string
     * @return subclass's human-readable name
     */
    virtual std::string display() const = 0;

    /**
     * Function ixamax is an abstraction of calls isamax/idamax/etc.
     * It returns the index of the first element having maximum absolute value.
     *
     * See also
     * https://www.netlib.org/lapack/explore-html/d0/d73/group__aux__blas_ga285793254ff0adaf58c605682efb880c.html
     *
     * @param[in] n number of elements of input vector
     * @param[in] x a C-style array containing (1 + (@param n-1) * abs(@param incx) elements of type @tparam floating
     * @param[in] incx storage spacing between elements of @param x
     * @return index of first element having maximum absolute value
     */
    virtual int ixamax(const int n, const floating * const x, const int incx) const = 0;

    /**
     * Function xaxpy is an abstraction of calls saxpy/daxpy/etc.
     * It performs the calculation y := alpha * x + y
     *
     * See also
     * https://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga8f99d6a644d3396aa32db472e0cfc91c.html
     *
     * @param[in] n number of elements in input vectors
     * @param[in] alpha scalar alpha
     * @param[in] x a C-style array containing (1 + (@param n-1) * abs(@param incx) elements of type @tparam floating
     * @param[in] incx storage spacing between elements of @param x
     * @param[in,out] y a C-style array containing (1 + (@param n-1) * abs(@param incy) elements of type @tparam floating
     * @param[in] incy storage spacing between elements of @param y
     */
    virtual void xaxpy(const int n, const floating alpha, const floating * const x, const int incx, floating * const y, const int incy) const = 0;

    /**
     * Function xcopy is an abstraction of calls scopy/dcopy/etc.
     * It performs the copy operation dest := source
     *
     * See also
     * https://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga21cdaae1732dea58194c279fca30126d.html
     * @param[in] n number of elements in input vectors
     * @param[in] source a C-style array containing (1 + (@param n-1) * abs(@param inc_source) elements of type @tparam floating
     * @param[in] inc_source storage spacing between elements of @param source
     * @param[out] dest a C-style array containing (1 + (@param n-1) * abs(@param inc_dest) elements of type @tparam floating
     * @param[in] inc_dest storage spacing between elements of @param dest
     */
    virtual void xcopy(const int n, const floating * const source, const int inc_source, floating * const dest, const int inc_dest) const = 0;

    /**
     * Function xdot is an abstraction of calls sdot/ddot/etc.
     * It performs the dot product <x, y>
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga75066c4825cb6ff1c8ec4403ef8c843a.html
     *
     * @param[in] n number of elements in input vectors
     * @param[in] x a C-style array containing (1 + (@param n-1) * abs(@param incx) elements of type @tparam floating
     * @param[in] incx spacing between elements of @param x
     * @param[in,out] y a C-style array containing (1 + (@param n-1) * abs(@param incy) elements of type @tparam floating
     * @param[in] incy spacing between elements of @param y
     * @return the value of the dot product <x, y>
     */
    virtual floating xdot(const int n, const floating * const x, const int incx, const floating * const y, const int incy) const = 0;

    /**
     * Function xgemm is an abstraction of calls sgemm/dgemm/etc.
     * It performs the matrix product operation C := alpha * trans_A(A) * trans_B(B) + beta * C
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/d1/d54/group__double__blas__level3_gaeda3cbd99c8fb834a60a6412878226e1.html
     *
     * @param[in] trans_A defines the transformation applied to matrix A (see enum class OperationType)
     * @param[in] trans_B defines the transformation applied to matrix B (see enum class OperationType)
     * @param[in] m the number of rows of matrix trans_A(A) and C
     * @param[in] n the number of columns of matrix trans_B(B) and C
     * @param[in] k the number of columns of matrix C
     * @param[in] alpha scalar applied to A
     * @param[in] A matrix stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] lda the first dimension of A
     * @param[in] B matrix stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] ldb the first dimension of B
     * @param[in] beta scalar applied to C
     * @param[in,out] C matrix stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] ldc the first dimension of C
     */
    virtual void xgemm(const OperationType trans_A, const OperationType trans_B,
                       const int m, const int n, const int k, const floating alpha, const floating * const A,
                       const int lda, const floating * const B, const int ldb, const floating beta, floating * const C, const int ldc) const = 0;

    /**
     * Function xgemv is an abstraction of calls sgemv/dgemv/etc.
     * It performs the matrix times vector operation y := alpha * trans(A) * x + beta * y
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/d7/d15/group__double__blas__level2_gadd421a107a488d524859b4a64c1901a9.html
     *
     * @param[in] trans defines the transformation applied to matrix A (see enum class OperationType)
     * @param[in] m the number of rows of matrix A
     * @param[in] n the number of columns of matrix A
     * @param[in] alpha scalar applied to A
     * @param[in] A matrix stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] lda the first dimension of matrix A
     * @param[in] x A C-style array containing (1 + (@param DIM-1) * abs(@param incx) elements of type @tparam floating. DIM is @param n if trans is OperationType::Identical, and @param m else.
     * @param[in] incx storage spacing between elements of @param x
     * @param[in] beta scalar applied to @param y
     * @param[in,out] y A C-style array containing (1 + (@param DIM-1) * abs(@param incy) elements of type @tparam floating. DIM is @param n if trans is OperationType::Identical, and @param m else.
     * @param[in] incy storage spacing between elements of @param y
     */
    virtual void xgemv(const OperationType trans, const int m, const int n,
                       const floating alpha, floating const * const A, const int lda, floating const * const x, const int incx,
                       const floating beta, floating * const y, const int incy) const = 0;

    /**
     * Function xgetrf is an abstraction of calls sgetrf/dgetrf/etc.
     * It computes the PLU factorization of A general M-by-N matrix A. A is overwritten by L and U,
     * and L has unit diagonal elements.
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html
     *
     * @param[in] m the number of rows of matrix A
     * @param[in] n the number of columns of matrix A
     * @param[in,out] A matrix stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] lda the first dimension of matrix A
     * @param[out] ipiv permutation array, row i of matrix A is interchanged with row ipiv[i]
     * @param[out] info If == 0, successful exit. If < 0, the i-th argument had illegal value. If > 0, U(i, i) is exactly 0 and cannot be used to solve A system of equations.
     */
    virtual void xgetrf(int * const m, int * const n, floating * const A, int * const lda,
                        int * const ipiv, int * const info) const = 0;

    /**
     * Function xgetri is an abstraction of calls sgetri/dgetri/etc.
     * It computes the inverse of A square matrix using its PLU factorization as calculated by xgetrf.
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga56d9c860ce4ce42ded7f914fdb0683ff.html
     *
     * @param[in] n order of square matrix A
     * @param[in,out] A LU-factorization matrix (as calculated by xgetrf) stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] lda the first dimension of matrix A
     * @param[in] ipiv permutation array from xgetrf
     * @param[out] work is temporary working space in form of A C-style array containing elements of type floating. On exit, if info == 0, work[1] returns optimal @param lwork
     * @param[in] lwork dimension of array @param work. If == -1, the routine calculates the optimal size of array @param work.
     * @param[out] info If == 0, successful exit. If < 0, the i-th argument had illegal value. If > 0, U(i, i) is exactly 0 and cannot be used to solve A system of equations.
     */
    virtual void xgetri(const int * const n, floating * const A, const int * const lda, const int * const ipiv,
                        floating * const work, const int * const lwork, int * const info) const = 0;

    /**
     * Function xgetrs is an abstraction of calls sgetrs/dgetrs/etc.
     * It solves A system of linear equations trans(A) * X = B using its PLU factorization as calculated by xgetrf.
     * Matrix A has to be square, and the result is stored in B.
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga58e332cb1b8ab770270843221a48296d.html
     *
     * @param[in] trans defines the transformation applied to matrix A (see enum class OperationType)
     * @param[in] n order of square matrix A
     * @param[in] nrhs number of right hand sides, i.e. the number of columns of B
     * @param[in] A LU-factorization matrix (as calculated by xgetrf) stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] lda the first dimension of matrix A
     * @param[in] ipiv permutation array from xgetrf
     * @param[in,out] B matrix of right hand sides stored as A column-wise C-style array containing elements of type @tparam floating
     * @param[in] ldb the first dimension of matrix B
     * @param[out] info If == 0, successful exit. If < 0, the i-th argument had illegal value.
     */
    virtual void xgetrs(const OperationType trans, const int * const n, const int * const nrhs,
                        const floating * const A, const int * const lda, const int * const ipiv,
                        floating * const B, const int * const ldb, int * const info) const = 0;

    /**
     * Function xscal is an abstraction of calls sscal/dscal/etc.
     * It performs the scaling operation x := alpha * x.
     *
     * See also
     * http://www.netlib.org/lapack/explore-html/de/da4/group__double__blas__level1_ga793bdd0739bbd0e0ec8655a0df08981a.html
     *
     * @param n number of elements in input vectors
     * @param alpha scalar alpha
     * @param x a C-style array containing (1 + (@param n-1) * abs(@param incx) elements of type @tparam floating
     * @param incx storage spacing between elements of @param x
     */
    virtual void xscal(const int n, const floating alpha, floating * const x, const int incx) const = 0;

    /**
     * Returns true if @tparam floating is float.
     * @return true if @tparam floating is float
     */
    constexpr static bool isFloat();

    /**
     * Returns true if @tparam floating is double.
     * @return true if @tparam floating is double
     */
    constexpr static bool isDouble();
};

/**
 * Class Cpu provides linear algebra operations on the Cpu using
 * BLAS/LAPACK/LAPACKE. All memory allocations also take place in
 * ordinary RAM.
 *
 * @tparam floating floating point type, either float or double
 */
template<class floating>
class Cpu : public ProcessingUnitDevice<floating> {
public:
    using ProcessingUnitDevice<floating>::isFloat;
    using ProcessingUnitDevice<floating>::isDouble;

    /**
     * Constructor
     */
    Cpu() = default;

    /**
     * Destructor
     */
    ~Cpu() override = default;

    Timer create_timer() const override;
    MemoryManager get_memory_manager() const override;
    std::string display() const override;

    int ixamax(const int n, const floating * const x, const int incx) const override;
    void xaxpy(const int n, const floating alpha, const floating * const x, const int incx, floating * const y, const int incy) const override;
    void xcopy(const int n, const floating * const source, const int inc_source, floating * const dest, const int inc_dest) const override;
    floating xdot(const int n, const floating * const x, const int incx, const floating * const y, const int incy) const override;
    void xgemm(const OperationType trans_A, const OperationType trans_B,
               const int M, const int N, const int K, const floating alpha, const floating * const A,
               const int lda, const floating * const B, const int ldb, const floating beta, floating * const C, const int ldc) const override;
    void xgemv(const OperationType trans, const int m, const int n,
               const floating alpha, floating const * const A, const int lda, floating const * const x, const int incx,
               const floating beta, floating * const y, const int incy) const override;
    void xgetrf(int * const m, int * const n, floating * const A, int * const lda,
                int * const ipiv, int * const info) const override;
    void xgetri(const int * const n, floating * const A, const int * const lda, const int * const ipiv,
                floating * const work, const int * const lwork, int * const info) const override;
    void xgetrs(const OperationType trans, const int * const n, const int * const nrhs,
                const floating * const A, const int * const lda, const int * const ipiv,
                floating * const B, const int * const ldb, int * const info) const override;
    void xscal(const int N, const floating alpha, floating * const x, const int incx) const override;

private:
    static std::map<OperationType, CBLAS_TRANSPOSE> to_internal_operation_blas; //!< Converts the OperationType to the corresponding CBLAS_TRANSPOSE type
    static std::map<OperationType, const char * const> to_internal_operation_lapack; //!< Converts the OperationType to the corresponding LAPACK transposition type
    static MemoryManager _device_manager; //!< Memory manager associated with Cpu memory
};

#ifndef CPU_ONLY

/**
 * Class Gpu provides linear algebra operations on the Gpu using
 * cuBLAS/cuSOLVER. All memory allocations also take place in
 * Gpu RAM.
 *
 * @tparam floating floating point type, either float or double
 */
template<class floating>
class Gpu : public ProcessingUnitDevice<floating> {
public:
    using ProcessingUnitDevice<floating>::isFloat;
    using ProcessingUnitDevice<floating>::isDouble;

    /**
     * Constructor
     */
    Gpu() = default;

    /**
     * Destructor
     */
    ~Gpu() override = default;

    Timer create_timer() const override;
    MemoryManager get_memory_manager() const override;
    std::string display() const override;

    /**
     * Returns the instance's cuBLAS handle.
     * @return cuBLAS handle
     */
    auto get_cublas_handle() const;

    /**
     * Returns the instance's cuSOLVER handle.
     * @return cuSOLVER handle
     */
    auto get_cusolver_handle() const;

    int ixamax(const int n, const floating * const x, const int incx) const override;
    void xaxpy(const int n, const floating alpha, const floating * const x, const int incx, floating * const y, const int incy) const override;
    void xcopy(const int n, const floating * const source, const int inc_source, floating * const dest, const int inc_dest) const override;
    floating xdot(const int n, const floating * const x, const int incx, const floating * const y, const int incy) const override;
    void xgemm(const OperationType trans_A, const OperationType trans_B,
               const int M, const int N, const int K, const floating alpha, const floating * const A,
               const int lda, const floating * const B, const int ldb, const floating beta, floating * const C, const int ldc) const override;
    void xgemv(const OperationType trans, const int m, const int n,
               const floating alpha, floating const * const A, const int lda, floating const * const x, const int incx,
               const floating beta, floating * const y, const int incy) const override;
    void xgetrf(int * const m, int * const n, floating * const A, int * const lda,
                int * const ipiv, int * const info) const override;
    void xgetri(const int * const n, floating * const A, const int * const lda, const int * const ipiv,
                floating * const work, const int * const lwork, int * const info) const override;
    void xgetrs(const OperationType trans, const int * const n, const int * const nrhs,
                const floating * const A, const int * const lda, const int * const ipiv,
                floating * const B, const int * const ldb, int * const info) const override;
    void xscal(const int N, const floating alpha, floating * const x, const int incx) const override;

private:
    static std::map<OperationType, cublasOperation_t> to_internal_operation; //!< Converts the OperationType to the corresponding cuBLAS operation type
    static GpuHandle _handle; //!< Gpu handle for cuBLAS/cuSOLVER operations
    static MemoryManager _deviceManager; //!< Memory manager associated with CUDA memory

};

#ifdef MAGMA
/**
 * Class GpuMagma provides linear algebra operations on the Gpu using
 * MAGMA. All memory allocations also take place in Gpu RAM.
 *
 * @tparam floating floating point type, either float or double
 */
template<class floating>
class GpuMagma : public ProcessingUnitDevice<floating> {
public:
    using ProcessingUnitDevice<floating>::isFloat;
    using ProcessingUnitDevice<floating>::isDouble;

    /**
     * Constructor
     */
    GpuMagma() = default;

    /**
     * Destructor
     */
    ~GpuMagma() override = default;

    Timer create_timer() const override;
    MemoryManager get_memory_manager() const override;
    std::string display() const override;

    /**
     * Returns the instance's magma queue.
     */
    auto get_magma_queue() const;

    int ixamax(const int n, const floating * const x, const int incx) const override;
    void xaxpy(const int n, const floating alpha, const floating * const x, const int incx, floating * const y, const int incy) const override;
    void xcopy(const int n, const floating * const source, const int inc_source, floating * const dest, const int inc_dest) const override;
    floating xdot(const int n, const floating * const x, const int incx, const floating * const y, const int incy) const override;
    void xgemm(const OperationType trans_A, const OperationType trans_B,
               const int M, const int N, const int K, const floating alpha, const floating * const A,
               const int lda, const floating * const B, const int ldb, const floating beta, floating * const C, const int ldc) const override;
    void xgemv(const OperationType trans, const int m, const int n,
               const floating alpha, floating const * const A, const int lda, floating const * const x, const int incx,
               const floating beta, floating * const y, const int incy) const override;
    void xgetrf(int * const m, int * const n, floating * const A, int * const lda,
                int * const ipiv, int * const info) const override;
    void xgetri(const int * const n, floating * const A, const int * const lda, const int * const ipiv,
                floating * const work, const int * const lwork, int * const info) const override;
    void xgetrs(const OperationType trans, const int * const n, const int * const nrhs,
                const floating * const A, const int * const lda, const int * const ipiv,
                floating * const B, const int * const ldb, int * const info) const override;
    void xscal(const int N, const floating alpha, floating * const x, const int incx) const override;

private:
    static std::map<OperationType, magma_trans_t> to_internal_operation; //!< Converts the OperationType to the corresponding MAGMA transposition type
    static GpuMagmaQueue _queue; //!< MAGMA queue for MAGMA operations
    static MemoryManager _device_manager; //!< Memory manager associated with CUDA memory
};

/**
 * Class ReplacementNumber is used for passing the operation to GpuMixed, which is then calculated via
 * GpuMagma instead of ordinary Gpu processing unit.
 */
enum class ReplacementNumber {
    IXAMAX,
    XAXPY,
    XCOPY,
    XDOT,
    XGEMM,
    XGEMV,
    XGETRF,
    XGETRS,
    XSCAL
};

/**
 * Class GpuMixed provides linear algebra operations on the Gpu using
 * cuBLAS/cuSOLVER, except for the one operation specified by replacementNumber,
 * which is calculated using MAGMA.
 * All memory allocations also take place in Gpu RAM.
 *
 * @tparam floating floating point type, either float or double
 */
template<class floating>
class GpuMixed : public ProcessingUnitDevice<floating> {
public:
    /**
     * Constructor
     */
    GpuMixed(const ReplacementNumber replacementNumber = ReplacementNumber::XGETRF);

    /**
     * Destructor
     */
    ~GpuMixed() override = default;

    Timer create_timer() const override;

    MemoryManager get_memory_manager() const override;

    std::string display() const override;

    int ixamax(const int n, const floating * const x, const int incx) const override;
    void xaxpy(const int n, const floating alpha, const floating * const x, const int incx, floating * const y, const int incy) const override;
    void xcopy(const int n, const floating * const source, const int inc_source, floating * const dest, const int inc_dest) const override;
    floating xdot(const int n, const floating * const x, const int incx, const floating * const y, const int incy) const override;
    void xgemm(const OperationType trans_A, const OperationType trans_B,
               const int M, const int N, const int K, const floating alpha, const floating * const A,
               const int lda, const floating * const B, const int ldb, const floating beta, floating * const C, const int ldc) const override;
    void xgemv(const OperationType trans, const int m, const int n,
               const floating alpha, floating const * const A, const int lda, floating const * const x, const int incx,
               const floating beta, floating * const y, const int incy) const override;
    void xgetrf(int * const m, int * const n, floating * const A, int * const lda,
                int * const ipiv, int * const info) const override;
    void xgetri(const int * const n, floating * const A, const int * const lda, const int * const ipiv,
                floating * const work, const int * const lwork, int * const info) const override;
    void xgetrs(const OperationType trans, const int * const n, const int * const nrhs,
                const floating * const A, const int * const lda, const int * const ipiv,
                floating * const B, const int * const ldb, int * const info) const override;
    void xscal(const int N, const floating alpha, floating * const x, const int incx) const override;

    ReplacementNumber get_replacement_number() const;

private:
    static Gpu<floating> _gpu; //!< Gpu processing unit
    static GpuMagma<floating> _gpu_magma; //!< MAGMA processing unit
    ReplacementNumber _replacementNumber; //!< current instance's operation which should be exectued on MAGMA instead of cuBLAS/cuSOLVER
};
#endif
#endif

#include "processingunit.hpp"
#include "processingunit_cpu.hpp"

#ifndef CPU_ONLY
    #include "processingunit_gpu.hpp"
    #ifdef MAGMA
        #include "processingunit_gpu_magma.hpp"
        #include "processingunit_gpu_mixed.hpp"
    #endif
#endif

#endif

ext: h

max_stars_repo_name: meteoritenhagel/fractional-pde

max_issues_repo_name: meteoritenhagel/fractional-pde

max_forks_repo_name: meteoritenhagel/fractional-pde


text:
#pragma once

#include <list>
#include <vector>
#include <optional>
#include <memory>
#include <functional>

#include <gsl/gsl-lite.hpp>

#include "UtilMacros.Namespace.h"
#include "util.std_aliases.h"
#include "zstring_view.h"

NAMESPACE_BEGIN( cpp )

using vlr::string;
using vlr::wstring;
using vlr::tstring;
using vlr::zstring_view;
using vlr::wzstring_view;
using vlr::tzstring_view;

using std::shared_ptr;
using std::weak_ptr;
using std::make_shared;

using std::function;

using gsl::span;

NAMESPACE_END //( cpp )


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  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.
 *
 * @brief interface of Table
 * @file Table.h
 * @author: xingqiangbai
 * @date: 2021-04-07
 */
#pragma once

#include "../../interfaces/protocol/ProtocolTypeDef.h"
#include "../../libutilities/Error.h"
#include "boost/algorithm/string.hpp"
#include "tbb/spin_mutex.h"
#include "tbb/spin_rw_mutex.h"
#include <boost/throw_exception.hpp>
#include <algorithm>
#include <any>
#include <cstdlib>
#include <gsl/span>
#include <map>
#include <memory>
#include <string>
#include <utility>
#include <vector>

#define STORAGE_LOG(LEVEL) BCOS_LOG(LEVEL) << "[STORAGE]"

namespace bcos
{
namespace storage
{
enum StorageError
{
    UnknownError = -60000,
    TableNotExists,
    SystemTableNotExists,
    TableExists,
    UnknownEntryType,
    ReadError,
    WriteError,
    EmptyStorage,
    ReadOnly,
};

struct Condition
{
    Condition() = default;
    ~Condition() = default;
    void NE(const std::string& value) { m_conditions.emplace_back(Comparator::NE, value); }
    // string compare, "2" > "12"
    void GT(const std::string& value) { m_conditions.emplace_back(Comparator::GT, value); }
    void GE(const std::string& value) { m_conditions.emplace_back(Comparator::GE, value); }
    // string compare, "12" < "2"
    void LT(const std::string& value) { m_conditions.emplace_back(Comparator::LT, value); }
    void LE(const std::string& value) { m_conditions.emplace_back(Comparator::LE, value); }
    void limit(size_t start, size_t end) { m_limit = std::pair<size_t, size_t>(start, end); }

    std::pair<size_t, size_t> getLimit() const { return m_limit; }

    bool isValid(const std::string_view& key) const
    {  // all conditions must be satisfied
        for (auto& cond : m_conditions)
        {  // conditions should few, so not parallel check for now
            switch (cond.cmp)
            {
            case Comparator::NE:
                if (key == cond.value)
                {
                    return false;
                }
                break;
            case Comparator::GT:
                if (key <= cond.value)
                {
                    return false;
                }
                break;
            case Comparator::GE:
                if (key < cond.value)
                {
                    return false;
                }
                break;
            case Comparator::LT:
                if (key >= cond.value)
                {
                    return false;
                }
                break;
            case Comparator::LE:
                if (key > cond.value)
                {
                    return false;
                }
                break;
            default:
                // undefined Comparator
                break;
            }
        }
        return true;
    }

    enum class Comparator
    {
        EQ,
        NE,
        GT,
        GE,
        LT,
        LE,
    };
    struct cond
    {
        cond(Comparator _cmp, const std::string& _value) : cmp(_cmp), value(_value) {}
        Comparator cmp;
        std::string value;
    };

    std::vector<cond> m_conditions;
    std::pair<size_t, size_t> m_limit;
};

class TableInfo
{
public:
    using Ptr = std::shared_ptr<TableInfo>;
    using ConstPtr = std::shared_ptr<const TableInfo>;

    TableInfo(std::string name, std::vector<std::string> fields)
      : m_name(std::move(name)), m_fields(std::move(fields))
    {
        m_order.reserve(m_fields.size());
        for (size_t i = 0; i < m_fields.size(); ++i)
        {
            m_order.push_back({m_fields[i], i});
        }
        std::sort(m_order.begin(), m_order.end(),
            [](auto&& lhs, auto&& rhs) { return std::get<0>(lhs) < std::get<0>(rhs); });
    }

    std::string_view name() const { return m_name; }

    const std::vector<std::string>& fields() const { return m_fields; }

    size_t fieldIndex(const std::string_view& field) const
    {
        auto it = std::lower_bound(m_order.begin(), m_order.end(), field,
            [](auto&& lhs, auto&& rhs) { return std::get<0>(lhs) < rhs; });
        if (it != m_order.end() && std::get<0>(*it) == field)
        {
            return std::get<1>(*it);
        }
        else
        {
            BOOST_THROW_EXCEPTION(
                BCOS_ERROR(-1, std::string("Can't find field: ") + std::string(field)));
        }
    }

private:
    std::string m_name;
    std::vector<std::string> m_fields;
    std::vector<std::tuple<std::string_view, size_t>> m_order;

private:
    void* operator new(size_t s) { return malloc(s); };
    void operator delete(void* p) { free(p); };
};

}  // namespace storage
}  // namespace bcos


text:
#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_ntuple.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

struct data
{
  int num;
  double x;
  double y;
  double z;
};
int sel_func (void *ntuple_data, void * params);
double val_func (void *ntuple_data, void * params);

int
main (void)
{
  struct data ntuple_row;
  int i;

  double x[1000], y[1000], z[1000], f[100];

  gsl_ntuple_select_fn S;
  gsl_ntuple_value_fn V;

  double scale = 1.5;

  gsl_ieee_env_setup ();

  /* zero struct including padding bytes to avoid valgrind errors */
  memset(&ntuple_row, 0, sizeof(struct data));

  S.function = &sel_func;
  S.params = &scale;

  V.function = &val_func;
  V.params = &scale;

  {
    gsl_ntuple *ntuple = gsl_ntuple_create ("test.dat", &ntuple_row,
                                            sizeof (ntuple_row));

    int status = 0;

    for (i = 0; i < 100; i++) f[i] = 0;

    for (i = 0; i < 1000; i++)
      {
        double xi = 1.0 / (i + 1.5);
        double yi = xi * xi ;
        double zi = xi * xi * xi;

        ntuple_row.x = xi;
        ntuple_row.y = yi;
        ntuple_row.z = zi;
        ntuple_row.num = i;

        x[i] = xi; y[i] = yi; z[i] = zi;

        if (xi * scale < 0.1)
          {
            double v = xi + yi + zi;
            int k = (int)(100.0*v*scale);
            f[k]++;
          }

        /* printf ("x,y,z = %f,%f,%f; n=%x \n", ntuple_row.x,
           ntuple_row.y, ntuple_row.z, ntuple_row.num); */

        {
          int s = gsl_ntuple_bookdata (ntuple);

          if (s != GSL_SUCCESS)
            {
              status = 1;
            }
        }
      }

    gsl_ntuple_close (ntuple);

    gsl_test (status, "writing ntuples");
  }

  {
    gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row,
                                          sizeof (ntuple_row));
    int status = 0;

    for (i = 0; i < 1000; i++)
      {
        gsl_ntuple_read (ntuple);

        status = (ntuple_row.num != i);
        status |= (ntuple_row.x != x[i]);
        status |= (ntuple_row.y != y[i]);
        status |= (ntuple_row.z != z[i]);

        /* printf ("x,y,z = %f,%f,%f; n=%d\n", ntuple_row.x,
                ntuple_row.y, ntuple_row.z, ntuple_row.num); */
      }
    gsl_ntuple_close (ntuple);

    gsl_test (status, "reading ntuples");
  }

  {
    int status = 0;

    gsl_ntuple *ntuple = gsl_ntuple_open ("test.dat", &ntuple_row,
                                          sizeof (ntuple_row));

    gsl_histogram *h = gsl_histogram_calloc_uniform (100, 0., 1.);

    gsl_ntuple_project (h, ntuple, &V, &S);

    gsl_ntuple_close (ntuple);

    /* gsl_histogram_fprintf (stdout, h, "%f", "%f"); */

    for (i = 0; i < 100; i++)
      {
        /* printf ("h  %g f  %g\n", h->bin[i], f[i]); */

        if (h->bin[i] != f[i])
          {
            status = 1;
          }
      }

    gsl_test (status, "histogramming ntuples");

    gsl_histogram_free (h);
  }

  exit (gsl_test_summary());
}

int
sel_func (void *ntuple_data, void * params)
{
  double x, y, z, scale;
  scale = *(double *)params;

  x = ((struct data *) ntuple_data)->x;
  y = ((struct data *) ntuple_data)->y;
  z = ((struct data *) ntuple_data)->z;

  return (x*scale < 0.1);
}

double
val_func (void *ntuple_data, void * params)
{
  double x, y, z, scale;
  scale = *(double *)params;

  x = ((struct data *) ntuple_data)->x;
  y = ((struct data *) ntuple_data)->y;
  z = ((struct data *) ntuple_data)->z;

  return (x + y + z) * scale;
}

ext: c

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#include <math.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>

double t_test(const double *data1, int n1, const double *data2, int n2)
{
    double mean1, mean2, var1, var2, num, sp, denom, t, dof, p;
    mean1 = gsl_stats_mean(data1, 1, n1);
    mean2 = gsl_stats_mean(data2, 1, n2);
    var1 = gsl_stats_variance(data1, 1, n1);
    var2 = gsl_stats_variance(data2, 1, n2);
    num = mean1 - mean2;
    sp = sqrt((((n1 - 1) * var1) + ((n2 - 1) * var2)) / (n1 + n2 - 2));
    denom = sp * sqrt((1 / (double) n1) + (1 / (double) n2));
    t = fabs(num / denom);
    dof = (double) n1 + n2 - 2;
    p = (1 - gsl_cdf_tdist_P(t, dof)) * 2;
    return p;
}

max_stars_repo_name: kgoettler/python-c-ext

max_issues_repo_name: kgoettler/python-c-ext

max_forks_repo_name: kgoettler/python-c-ext


text:
/* siman/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_siman.h>
#include <gsl/gsl_ieee_utils.h>
#include <stdio.h>

/* set up parameters for this simulated annealing run */
#define N_TRIES 200             /* how many points do we try before stepping */
#define ITERS_FIXED_T 1000      /* how many iterations for each T? */
#define STEP_SIZE 1.0           /* max step size in random walk */
#define K 1.0                   /* Boltzmann constant */
#define T_INITIAL 0.008         /* initial temperature */
#define MU_T 1.003              /* damping factor for temperature */
#define T_MIN 2.0e-6

gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
                             K, T_INITIAL, MU_T, T_MIN};

double square (double x) ;
double square (double x) { return x * x ; }

double E1(void *xp);
double M1(void *xp, void *yp);
void S1(const gsl_rng * r, void *xp, double step_size);
void P1(void *xp);

/* now some functions to test in one dimension */
double E1(void *xp)
{
  double x = * ((double *) xp);

  return exp(-square(x-1))*sin(8*x) - exp(-square(x-1000))*0.89;
}

double M1(void *xp, void *yp)
{
  double x = *((double *) xp);
  double y = *((double *) yp);

  return fabs(x - y);
}

void S1(const gsl_rng * r, void *xp, double step_size)
{
  double old_x = *((double *) xp);
  double new_x;

  new_x = gsl_rng_uniform(r)*2*step_size - step_size + old_x;

  memcpy(xp, &new_x, sizeof(new_x));
}

void P1(void *xp)
{
  printf(" %12g ", *((double *) xp));
}

int main(void)
{
  double x_min = 1.36312999455315182 ;
  double x ;

  gsl_rng * r = gsl_rng_alloc (gsl_rng_env_setup()) ;

  gsl_ieee_env_setup ();

  /* The function tested here has multiple mimima.
     The global minimum is at    x = 1.36312999, (f = -0.87287)
     There is a local minimum at x = 0.60146196, (f = -0.84893) */

  x = -10.0 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
                  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=-10") ;

  x = +10.0 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
                  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=10") ;

  /* Start at the false minimum */

  x = +0.6 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
                  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.6") ;

  x = +0.5 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
                  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.5") ;

  x = +0.4 ;
  gsl_siman_solve(r, &x, E1, S1, M1, NULL, NULL, NULL, NULL,
                  sizeof(double), params);
  gsl_test_rel(x, x_min, 1e-3, "f(x)= exp(-(x-1)^2) sin(8x), x0=0.4") ;

  gsl_rng_free(r);
  exit (gsl_test_summary ());

#ifdef JUNK
  x0.D1 = 12.0;
  printf("#one dimensional problem, x0 = %f\n", x0.D1);
  gsl_siman_Usolve(r, &x0, test_E_1D, test_step_1D, distance_1D,
                   print_pos_1D, params);

  x0.D2[0] = 12.0;
  x0.D2[1] = 5.5;
  printf("#two dimensional problem, (x0,y0) = (%f,%f)\n",
         x0.D2[0], x0.D2[1]);
  gsl_siman_Usolve(r, &x0, test_E_2D, test_step_2D, distance_2D,
                   print_pos_2D, params);

  x0.D3[0] = 12.2;
  x0.D3[1] = 5.5;
  x0.D3[2] = -15.5;
  printf("#three dimensional problem, (x0,y0,z0) = (%f,%f,%f)\n",
         x0.D3[0], x0.D3[1], x0.D3[2]);
  gsl_siman_Usolve(r, &x0, test_E_3D, test_step_3D, distance_3D,
                   print_pos_3D, params);

  x0.D2[0] = 12.2;
  x0.D2[1] = 5.5;

  gsl_siman_solve(r, &x0, test_E_2D, test_step_2D, distance_2D, print_pos_2D, params);

  x0.D3[0] = 12.2;
  x0.D3[1] = 5.5;
  x0.D3[2] = -15.5;

  gsl_siman_solve(r, &x0, test_E_3D, test_step_3D, distance_3D, print_pos_3D, params);

  return 0;
#endif
}


text:
#ifndef SYNAPSE_H_
#define SYNAPSE_H_

#include "spike_train.h"
#include <gsl/gsl_rng.h>

namespace neurophys {

class DynamicSynapse {
public:
	virtual SpikeTrain filter(const SpikeTrain& in) = 0;
};

class FDSynapse: public DynamicSynapse {
public:
	static const int STATIC = 0;
	static const int DEPRESSING = 1;
	static const int FACILITATING = 2;
	static const int GENERAL = DEPRESSING | FACILITATING;
	static const int KEEP_STATE_VARIABLES = 4; // no reset between trials

	FDSynapse(double F0, double tau_f, double tau_d, double delta, int flags):
		F0_(F0), tau_f_(tau_f), tau_d_(tau_d), delta_(delta), flags_(flags),
        FC0_(0), t_(0) {};

protected:
	const double F0_;
	const double tau_f_;
	const double tau_d_;
	const double delta_;
	const int flags_;
    double FC0_;
    double t_;
};

class DeterministicFDSynapse: public FDSynapse {
public:
	DeterministicFDSynapse(double F0, double tau_f, double tau_d, double delta,
           int flags):
		FDSynapse(F0, tau_f, tau_d, delta, flags), D0_(0) {};

	virtual SpikeTrain filter(const SpikeTrain& in);
private:
    double D0_;
};

// rng in constructor is not so nice -- make explicit as param of filter and
// screw poplymorphism?
class StochasticFDSynapse: public FDSynapse {
public:
	StochasticFDSynapse(double F0, double tau_f, double tau_d, double delta,
            int flags, gsl_rng* rng, unsigned int n_release_sites):
		FDSynapse(F0, tau_f, tau_d, delta, flags), rng_(rng), n_release_sites_(n_release_sites), n_ves_present_(n_release_sites) {};

	virtual SpikeTrain filter(const SpikeTrain& in);
private:
	gsl_rng* rng_;
	const unsigned int n_release_sites_;
	unsigned int n_ves_present_;
};

}

#endif // SYNAPSE_H_


text:
static char help[] = "Broker for DYNHELICS\n";

#include <stdio.h>
#include <unistd.h>
#include <helics.h>
#include <petsc.h>

int main(int argc,char **argv)
{
  helics_broker broker;
  const char*   helicsversion;
  int            isconnected;
  char           initstring[PETSC_MAX_PATH_LEN];
  PetscInt       nfederates=0;
  PetscErrorCode ierr;

  PetscInitialize(&argc,&argv,"petscopt",help);

  helicsversion = helicsGetVersion();

  printf("BROKER: Helics version = %s\n",helicsversion);
  printf("%s",help);

  ierr = PetscOptionsGetInt(NULL,NULL,"-nfeds",&nfederates,NULL);CHKERRQ(ierr);
  if(!nfederates) {
    SETERRQ(PETSC_COMM_SELF,0,"Number of federates need to be given with option -nfeds");
  }
  ierr = PetscSNPrintf(initstring,PETSC_MAX_PATH_LEN-1,"%d",nfederates);
  ierr = PetscStrcat(initstring," --name=mainbroker");

  /* Create broker */
  broker = helicsCreateBroker("zmq","",initstring);

  isconnected = helicsBrokerIsConnected(broker);

  if(isconnected) {
    printf("BROKER: Created and connected\n");
  }

  while(helicsBrokerIsConnected(broker)) {
    usleep(1000); /* Sleep for 1 millisecond */
  }
  printf("BROKER: disconnected\n");
  helicsBrokerFree(broker);
  helicsCloseLibrary();

  printf("Helics library closed\n");
  PetscFinalize();
  return 0;
}

max_stars_repo_name: GMLC-TDC/Use-Cases

max_issues_repo_name: GMLC-TDC/Use-Cases

max_forks_repo_name: GMLC-TDC/Use-Cases


text:
/**
 *
 * @file qwrapper_cgemm.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:56 2014
 *
 **/
#include <cblas.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum transA, int transB,
                      int m, int n, int k, int nb,
                      PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                                                const PLASMA_Complex32_t *B, int ldb,
                      PLASMA_Complex32_t beta, PLASMA_Complex32_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    C,                 INOUT,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm2( Quark *quark, Quark_Task_Flags *task_flags,
                        PLASMA_enum transA, int transB,
                        int m, int n, int k, int nb,
                        PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                        const PLASMA_Complex32_t *B, int ldb,
                        PLASMA_Complex32_t beta, PLASMA_Complex32_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    C,                 INOUT | LOCALITY | GATHERV,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemm_quark = PCORE_cgemm_quark
#define CORE_cgemm_quark PCORE_cgemm_quark
#endif
void CORE_cgemm_quark(Quark *quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int m;
    int n;
    int k;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int lda;
    PLASMA_Complex32_t *B;
    int ldb;
    PLASMA_Complex32_t beta;
    PLASMA_Complex32_t *C;
    int ldc;

    quark_unpack_args_13(quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
    cblas_cgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        CBLAS_SADDR(alpha), A, lda,
        B, ldb,
        CBLAS_SADDR(beta), C, ldc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm_f2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                                                   const PLASMA_Complex32_t *B, int ldb,
                         PLASMA_Complex32_t beta, PLASMA_Complex32_t *C, int ldc,
                         PLASMA_Complex32_t *fake1, int szefake1, int flag1,
                         PLASMA_Complex32_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_f2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(PLASMA_Complex32_t)*szefake1, fake1,             flag1,
        sizeof(PLASMA_Complex32_t)*szefake2, fake2,             flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemm_f2_quark = PCORE_cgemm_f2_quark
#define CORE_cgemm_f2_quark PCORE_cgemm_f2_quark
#endif
void CORE_cgemm_f2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int LDA;
    PLASMA_Complex32_t *B;
    int LDB;
    PLASMA_Complex32_t beta;
    PLASMA_Complex32_t *C;
    int LDC;
    void *fake1, *fake2;

    quark_unpack_args_15(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1, fake2);
    cblas_cgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm_p2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                         const PLASMA_Complex32_t **B, int ldb,
                         PLASMA_Complex32_t beta, PLASMA_Complex32_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_p2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemm_p2_quark = PCORE_cgemm_p2_quark
#define CORE_cgemm_p2_quark PCORE_cgemm_p2_quark
#endif
void CORE_cgemm_p2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int LDA;
    PLASMA_Complex32_t **B;
    int LDB;
    PLASMA_Complex32_t beta;
    PLASMA_Complex32_t *C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_cgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        *B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm_p3(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                           const PLASMA_Complex32_t *B, int ldb,
                           PLASMA_Complex32_t beta, PLASMA_Complex32_t **C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_p3_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t)*ldb*nb,   B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t*),         C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemm_p3_quark = PCORE_cgemm_p3_quark
#define CORE_cgemm_p3_quark PCORE_cgemm_p3_quark
#endif
void CORE_cgemm_p3_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int LDA;
    PLASMA_Complex32_t *B;
    int LDB;
    PLASMA_Complex32_t beta;
    PLASMA_Complex32_t **C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_cgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        B, LDB,
        CBLAS_SADDR(beta), *C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cgemm_p2f1(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           PLASMA_Complex32_t alpha, const PLASMA_Complex32_t *A, int lda,
                           const PLASMA_Complex32_t **B, int ldb,
                           PLASMA_Complex32_t beta, PLASMA_Complex32_t *C, int ldc,
                           PLASMA_Complex32_t *fake1, int szefake1, int flag1)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_cgemm_p2f1_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex32_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex32_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex32_t*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex32_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex32_t)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(PLASMA_Complex32_t)*szefake1, fake1,             flag1,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemm_p2f1_quark = PCORE_cgemm_p2f1_quark
#define CORE_cgemm_p2f1_quark PCORE_cgemm_p2f1_quark
#endif
void CORE_cgemm_p2f1_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex32_t alpha;
    PLASMA_Complex32_t *A;
    int LDA;
    PLASMA_Complex32_t **B;
    int LDB;
    PLASMA_Complex32_t beta;
    PLASMA_Complex32_t *C;
    int LDC;
    void *fake1;

    quark_unpack_args_14(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1);
    cblas_cgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        *B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}


text:
#ifndef TETRA_DOS_ELHAMILTONIAN_H
#define TETRA_DOS_ELHAMILTONIAN_H

#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "environment.h"

void ElHamiltonian_Recip(Environment *env, double k[3], gsl_matrix_complex *H);

void ElHamiltonian(Environment *env, double k[3], gsl_matrix_complex *H);

gsl_complex EpsilonAE(Environment *env, double k[3]);

gsl_complex EpsilonBE(Environment *env, double k[3]);

gsl_complex EpsilonAO(Environment *env, double k[3]);

gsl_complex EpsilonBO(Environment *env, double k[3]);

#endif // TETRA_DOS_ELHAMILTONIAN_H


text:
/* rng/ranlux.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is a lagged fibonacci generator with skipping developed by Luescher.
   The sequence is a series of 24-bit integers, x_n,

   x_n = d_n + b_n

   where d_n = x_{n-10} - x_{n-24} - c_{n-1}, b_n = 0 if d_n >= 0 and
   b_n = 2^24 if d_n < 0, c_n = 0 if d_n >= 0 and c_n = 1 if d_n < 0,
   where after 24 samples a group of p integers are "skipped", to
   reduce correlations. By default p = 199, but can be increased to
   365.

   The period of the generator is around 10^171.

   From: M. Luescher, "A portable high-quality random number generator
   for lattice field theory calculations", Computer Physics
   Communications, 79 (1994) 100-110.

   Available on the net as hep-lat/9309020 at http://xxx.lanl.gov/

   See also,

   F. James, "RANLUX: A Fortran implementation of the high-quality
   pseudo-random number generator of Luscher", Computer Physics
   Communications, 79 (1994) 111-114

   Kenneth G. Hamilton, F. James, "Acceleration of RANLUX", Computer
   Physics Communications, 101 (1997) 241-248

   Kenneth G. Hamilton, "Assembler RANLUX for PCs", Computer Physics
   Communications, 101 (1997) 249-253  */

static inline unsigned long int ranlux_get (void *vstate);
static double ranlux_get_double (void *vstate);
static void ranlux_set_lux (void *state, unsigned long int s, unsigned int luxury);
static void ranlux_set (void *state, unsigned long int s);
static void ranlux389_set (void *state, unsigned long int s);

static const unsigned long int mask_lo = 0x00ffffffUL;	/* 2^24 - 1 */
static const unsigned long int mask_hi = ~0x00ffffffUL;
static const unsigned long int two24 = 16777216;	/* 2^24 */

typedef struct
  {
    unsigned int i;
    unsigned int j;
    unsigned int n;
    unsigned int skip;
    unsigned int carry;
    unsigned long int u[24];
  }
ranlux_state_t;

static inline unsigned long int increment_state (ranlux_state_t * state);

static inline unsigned long int
increment_state (ranlux_state_t * state)
{
  unsigned int i = state->i;
  unsigned int j = state->j;
  long int delta = state->u[j] - state->u[i] - state->carry;

  if (delta & mask_hi)
    {
      state->carry = 1;
      delta &= mask_lo;
    }
  else
    {
      state->carry = 0;
    }

  state->u[i] = delta;

  if (i == 0)
    {
      i = 23;
    }
  else
    {
      i--;
    }

  state->i = i;

  if (j == 0)
    {
      j = 23;
    }
  else
    {
      j--;
    }

  state->j = j;

  return delta;
}

static inline unsigned long int
ranlux_get (void *vstate)
{
  ranlux_state_t *state = (ranlux_state_t *) vstate;
  const unsigned int skip = state->skip;
  unsigned long int r = increment_state (state);

  state->n++;

  if (state->n == 24)
    {
      unsigned int i;
      state->n = 0;
      for (i = 0; i < skip; i++)
	increment_state (state);
    }

  return r;
}

static double
ranlux_get_double (void *vstate)
{
  return ranlux_get (vstate) / 16777216.0;
}

static void
ranlux_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
{
  ranlux_state_t *state = (ranlux_state_t *) vstate;
  int i;

  long int seed;

  if (s == 0)
    s = 314159265;	/* default seed is 314159265 */

  seed = s;

  /* This is the initialization algorithm of F. James, widely in use
     for RANLUX. */

  for (i = 0; i < 24; i++)
    {
      unsigned long int k = seed / 53668;
      seed = 40014 * (seed - k * 53668) - k * 12211;
      if (seed < 0)
	{
	  seed += 2147483563;
	}
      state->u[i] = seed % two24;
    }

  state->i = 23;
  state->j = 9;
  state->n = 0;
  state->skip = luxury - 24;

  if (state->u[23] & mask_hi)
    {
      state->carry = 1;
    }
  else
    {
      state->carry = 0;
    }
}

static void
ranlux_set (void *vstate, unsigned long int s)
{
  ranlux_set_lux (vstate, s, 223);
}

static void
ranlux389_set (void *vstate, unsigned long int s)
{
  ranlux_set_lux (vstate, s, 389);
}

static const gsl_rng_type ranlux_type =
{"ranlux",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlux_state_t),
 &ranlux_set,
 &ranlux_get,
 &ranlux_get_double};

static const gsl_rng_type ranlux389_type =
{"ranlux389",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlux_state_t),
 &ranlux389_set,
 &ranlux_get,
 &ranlux_get_double};

const gsl_rng_type *gsl_rng_ranlux = &ranlux_type;
const gsl_rng_type *gsl_rng_ranlux389 = &ranlux389_type;


text:
#ifndef HYPER_CONST
#define HYPER_CONST

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

class Hyper_Const
{
public :
    static const int dim0 = 120;           // max_clause
    static const int dim1 = 20;            // max_var
    static const int dim2 = 2;             // nc
    static const int nact = 40;            // nact
    static const float c_act;     	       // c_act is a hyperparameter for MCTS (decide the level of exploration)

    static const gsl_rng *r;              // random generater
    static const double alpha[nact];      // alpha parameter
    static void generate_dirichlet(double*);     // function used to generate dirichlet noise
    static const int MCTS_size_lim; // the size of MCT we want to achieve.
};

#endif

max_stars_repo_name: dmeoli/MCTSminisat

max_issues_repo_name: dmeoli/MCTSminisat

max_forks_repo_name: dmeoli/MCTSminisat


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
/*#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit_nlin.h>*/
#include "ral_nlls.h"

/* define the usertype */
struct usertype {
  double y_data[67];
  double x_data[67];
};

void generate_data_example ( double *x_data, double *y_data, const int m); // prototype
int eval_F  ( const int n, const int m,  void *params,
      const double *X, double *f);
int eval_J  ( const int n, const int m,  void *params,
      const double *X, double *f);
int eval_HF ( const int n, const int m,  void *params,
	       const double *X, const double *f, double *hf);

/* A c driver for the ral_nlls program */
int main(void) {

  /* Problem data */
  const int n = 2;
  const int m = 67;

  /* Derived types */
  struct ral_nlls_options options;
  struct ral_nlls_inform status;
  struct usertype params;

  printf("===============\n");
  printf("RAL NLLS driver\n");
  printf("~  C version  ~\n");
  printf("===============\n");

  /* Generate the data... */
  generate_data_example(params.x_data,params.y_data,m);

  double x[n];
  x[0] = 1.0;
  x[1] = 2.0;

  ral_nlls_default_options(&options);

  options.print_level = 3;

  nlls_solve(n, m, x,
	   eval_F, eval_J, eval_HF, &params,
	     &options, &status, NULL, NULL, NULL, NULL );

  int i;
  printf("\nX = \n");
  for(i=0; i < n; i++) {
    printf("  %5.4f \n",x[i]);
  }

  return 0; /* success */
}

/* Do a function evaluation */
int eval_F(int n, int m, void *params,
	      const double *X, double *f){

  struct usertype *myparams = (struct usertype *) params;

  int i;

  for(i=0; i<m; i++) {
    f[i] = myparams->y_data[i] - exp( X[0] * myparams->x_data[i] + X[1] );
  }

  return 0;
}

/* Evaluate the Jacobian */
int eval_J( const int n, const int m, void *params,
	     const double *X, double *J){

  struct usertype *myparams = (struct usertype *) params;

  int i;

  for(i=0; i<m; i++) {
    J[i] = -myparams->x_data[i] * exp( X[0] * myparams->x_data[i] + X[1] );
    J[m + i] = - exp( X[0] * myparams->x_data[i] + X[1] );
  }

  return 0;
}

/* Evaluate the Hessian */
int eval_HF( const int n, const int m, void *params,
	     const double *X, const double *f,
	     double *hf){

  struct usertype *myparams = (struct usertype *) params;

  int i;

  for(i=0; i<n*n; i++) {
    hf[i] = 0.0;
  }

  return 0;
}

/* Generate some example data... */
void generate_data_example( double *x_data, double *y_data, const int m ) {

  int i;
  /* Note the 67 here needs to be hard-coded, and you can't initialize
     an array with a variable length in C */
  double tempx[67] =  { 0.0,
			0.075000000000000,
			0.150000000000000,
			0.225000000000000,
			0.300000000000000,
			0.375000000000000,
			0.450000000000000,
			0.525000000000000,
			0.600000000000000,
			0.675000000000000,
			0.750000000000000,
			0.825000000000000,
			0.900000000000000,
			0.975000000000000,
			1.050000000000000,
			1.125000000000000,
			1.200000000000000,
			1.275000000000000,
			1.350000000000000,
			1.425000000000000,
			1.500000000000000,
			1.575000000000000,
			1.650000000000000,
			1.725000000000000,
			1.800000000000000,
			1.875000000000000,
			1.950000000000000,
			2.025000000000000,
			2.100000000000000,
			2.175000000000000,
			2.250000000000000,
			2.325000000000000,
			2.400000000000000,
			2.475000000000000,
			2.550000000000000,
			2.625000000000000,
			2.700000000000000,
			2.775000000000000,
			2.850000000000000,
			2.925000000000000,
			3.000000000000000,
			3.075000000000000,
			3.150000000000000,
			3.225000000000001,
			3.300000000000000,
			3.375000000000000,
			3.450000000000000,
			3.525000000000000,
			3.600000000000001,
			3.675000000000000,
			3.750000000000000,
			3.825000000000000,
			3.900000000000000,
			3.975000000000000,
			4.050000000000001,
			4.125000000000000,
			4.200000000000000,
			4.275000000000000,
			4.350000000000001,
			4.425000000000000,
			4.500000000000000,
			4.575000000000000,
			4.650000000000000,
			4.725000000000001,
			4.800000000000000,
			4.875000000000000,
			4.950000000000000};

  double tempy[67] =  {0.907946872110432,
		       1.199579396036134,
		       1.060092431384317,
		       1.298370500472354,
		       0.952768858414788,
		       1.209665290655204,
		       1.256912538155493,
		       1.163922146095987,
		       1.004877938808100,
		       1.205944250961060,
		       0.952693297695969,
		       1.449662692280761,
		       1.402015259144406,
		       1.378094012325746,
		       1.560882147577552,
		       1.437185539058121,
		       1.559853079888265,
		       1.877814947316832,
		       1.818781749024682,
		       1.375546045112591,
		       1.233967904388409,
		       1.887793124397751,
		       1.610237096463521,
		       1.787032484792262,
		       1.850015127982676,
		       2.120553361509177,
		       1.942913663511919,
		       2.106517132599766,
		       2.271787117356578,
		       1.727554346001754,
		       2.002909500898113,
		       1.975837413903495,
		       2.337446525801909,
		       1.960190841677278,
		       2.447097025572309,
		       2.161663720225506,
		       2.748798529374621,
		       2.507814238594416,
		       2.423769408403069,
		       2.578119353028746,
		       2.460310096221557,
		       2.638362783992324,
		       2.765540456237868,
		       2.837165966564409,
		       3.179711963042789,
		       3.245315453091675,
		       3.289631922410174,
		       3.360995198615834,
		       3.470489725998371,
		       3.169513520153466,
		       3.363740517933189,
		       3.665288099084969,
		       3.620334359722351,
		       4.018911445550667,
		       3.512715166706162,
		       3.874661411575566,
		       4.197746303653517,
		       3.703511523106007,
		       4.076351488309604,
		       4.056340365649961,
		       4.297751562451419,
		       4.373076571153739,
		       4.577093065941748,
		       4.856619059058190,
		       4.927350280596274,
		       4.703122139742729,
		       4.870205182453842};

  for(i=0;i<m;i++){
    x_data[i] = tempx[i];
    y_data[i] = tempy[i];
  }

}

max_stars_repo_name: ralna/RALFit

max_issues_repo_name: ralna/RALFit

max_forks_repo_name: ralna/RALFit


text:
/* chii - CHIsq test integers belonging to a given contiguous range. C
 * not R because a) a second opinion is good and b) I'm too lazy to
 * figure out how to make R honor sparse input, where a particular
 * bucket might have no values. otherwise, use the "chisq.test" in R or
 * the "equichisq" arlet of my r-fu script
 *
 * the following should fail, as there are six buckets (implicit 0 for
 * minimum through 5, inclusive) and the Perl only produces numbers for
 * five of those buckets (0..4)
 *
 *   perl -E 'for (1..1000) { say int rand 5}' | chii -M 5
 *
 * (it could also fail for reasons unrelated to the missing bucket, or
 * it could pass if there are not enough trials to reveal the lack)
 *
 * if you are actually testing the fitness of non-cryptographic hash
 * functions, consider instead the "smhasher" project */

#include <ctype.h>
#include <err.h>
#include <errno.h>
#include <getopt.h>
#include <locale.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sysexits.h>
#include <unistd.h>

// https://github.com/thrig/goptfoo
#include <goptfoo.h>

// GNU Science Library
#include <gsl/gsl_cdf.h>

double Flag_Fitness = 0.05;     // -F @ 95% confidence by default
long long Flag_Min;             // -m
long long Flag_Max;             // -M

double chisq(long long valuec, long long *values, long long minv,
             long long maxv, long long sum);
void emit_help(void);

int main(int argc, char *argv[])
{
    int ch;

    char *ep;
    long long number_of_nums = 0;
    long long range, *numstore, value;

    char *line = NULL;
    char *lp;
    FILE *fh;
    size_t linesize = 0;
    ssize_t linelen;
    ssize_t lineno = 1;

    double pvalue;

#ifdef __OpenBSD__
    if (pledge("stdio", NULL) == -1)
        err(1, "pledge failed");
#endif

    setlocale(LC_ALL, "C");

    while ((ch = getopt(argc, argv, "F:h?m:M:")) != -1) {
        switch (ch) {
        case 'F':
            Flag_Fitness = flagtod(ch, optarg, 0.0, 1.0);
            break;
        case 'm':
            Flag_Min = flagtoll(ch, optarg, LLONG_MIN, LLONG_MAX);
            break;
        case 'M':
            Flag_Max = flagtoll(ch, optarg, LLONG_MIN, LLONG_MAX);
            break;
        case 'h':
        case '?':
        default:
            emit_help();
            /* NOTREACHED */
        }
    }
    argc -= optind;
    argv += optind;

    if ((Flag_Min == 0 && Flag_Max == 0) || Flag_Min == Flag_Max
        || Flag_Min > Flag_Max)
        emit_help();

    if (argc == 0 || strncmp(*argv, "-", (size_t) 2) == 0) {
        fh = stdin;
    } else {
        if (!(fh = fopen(*argv, "r")))
            err(EX_IOERR, "could not open '%s'", *argv);
    }

    // inclusive, need buckets for all possible numbers
    range = Flag_Max - Flag_Min + 1;

    if (!(numstore = calloc((size_t) range, sizeof(long long))))
        err(EX_OSERR, "could not calloc() for %lld long longs", range);

    while ((linelen = getline(&line, &linesize, fh)) != -1) {
        // so can parse multiple numbers per line, assuming they are
        // isspace(3) delimited (downside: complicates the parsing)
        lp = line;
        while (*lp != '\0') {
            errno = 0;
            value = strtoll(lp, &ep, 10);
            if (*lp == '\0')
                errx(EX_DATAERR, "could not parse long long at line %ld",
                     lineno);
            if (errno == ERANGE && (value == LLONG_MIN || value == LLONG_MAX))
                errx(EX_DATAERR, "value out of range at line %ld", lineno);
            if (value < Flag_Min || value > Flag_Max)
                errx(EX_DATAERR, "value beyond bounds at line %ld", lineno);

            numstore[value - Flag_Min]++;
            number_of_nums++;

            // nom up not-strtoll() material so parser does not get wedged
            // on trailing material such as `echo -n 0 1 2 3 4" "`
            while (*ep != '+' && *ep != '-' && !isdigit(*ep)) {
                if (*ep == '\0' || *ep == '\n') {
                    *ep = '\0';
                    break;
                }
                ep++;
            }
            lp = ep;
        }
        lineno++;
    }
    if (number_of_nums == 0)
        errx(EX_DATAERR, "could not parse any numbers");
    if (ferror(fh))
        err(EX_IOERR, "ferror() reading input");

    pvalue = chisq(range, numstore, Flag_Min, Flag_Max, number_of_nums);
    printf("p-value %.6g%s\n", pvalue, pvalue <= Flag_Fitness ? " FAIL" : "");

    exit(EXIT_SUCCESS);
}

double chisq(long long valuec, long long *values, long long minv,
             long long maxv, long long sum)
{
    unsigned long delta, v;
    double expected, sum_of_squares;

    delta = maxv - minv + 1;
    expected = sum / (double) delta;

    sum_of_squares = 0.0;
    for (v = 0; v < valuec; v++) {
        sum_of_squares += pow(values[v] - expected, 2.0) / expected;
    }
    return gsl_cdf_chisq_Q(sum_of_squares, (double) delta - 1.0);
}

void emit_help(void)
{
    fputs("Usage: chii [-F fit] -m min -M max [file|-]\n", stderr);
    exit(EX_USAGE);
}

max_stars_repo_name: thrig/scripts

max_issues_repo_name: thrig/scripts

max_forks_repo_name: thrig/scripts


text:
// 2019-05-16 by Liya Ding

#include <string>
#include <math.h>
#include <gsl/gsl_poly.h>
#include <algorithm>
#include <stdlib.h>
#include <iostream>
//#include "convolver.h"

static double approxZero(double n);
static double opposite(double theta);
static int getTemplateN(int M);
double* getWeights(int M, double sigma);
void computeBaseTemplates(double* input, int nx, int ny, int M, int borderCondition, double sigma, double** templates);
double pointRespM1(int i, double angle, double* alpha, double** templates);
double pointRespM2(int i, double angle, double* alpha, double** templates);
double pointRespM3(int i, double angle, double* alpha, double** templates);
double pointRespM4(int i, double angle, double* alpha, double** templates);
double pointRespM5(int i, double angle, double* alpha, double** templates);

int getRealRoots(double* z, int nz, double* roots);
void filterM1(double** templates, int nx, int ny, double* alpha, double* response, double* orientation);
void filterM2(double** templates, int nx, int ny, double* alpha, double* response, double* orientation);
void filterM3(double** templates, int nx, int ny, double* alpha, double* response, double* orientation);
void filterM4(double** templates, int nx, int ny, double* alpha, double* response, double* orientation);
void filterM5(double** templates, int nx, int ny, double* alpha, double* response, double* orientation);
int mirror(int x, int nx);
double interp(double* image, int nx, int ny, double x, double y);
void computeNMS(double* response, double* orientation, double* nms, int nx, int ny);
void steerablefilter2Dcore(double * input, long* in_sz, int M, double sigma,double* &response, double* &orientation, double* &nms);


text:
#ifndef __STATS_H__
#define __STATS_H__

#include <stdint.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <inttypes.h>

#include <gsl/gsl_rstat.h>

typedef uint64_t word;

#define TIMER_STACK 8
enum Timer{MutTimer, Gen0Timer, Gen1Timer, MiscTimer, MaxTimer};

typedef struct {
  uint64_t allocated;
  uint64_t nursery_n_collections;
  uint64_t nursery_time_max;
  uint64_t nursery_copied;

  uint64_t gen1_collections;
  uint64_t gen1_copied;

  uint64_t max_heap;
  uint64_t max_residency;

  uint64_t start_time;
  enum Timer active_timer[TIMER_STACK];
  uint64_t n_timers;
  uint64_t timers[MaxTimer];

  gsl_rstat_quantile_workspace *latency_50;
  gsl_rstat_quantile_workspace *latency_90;
  gsl_rstat_quantile_workspace *latency_99;
} Stats;

void stats_init(Stats *s);
void stats_pprint(Stats *s);

static uint64_t clock_gettime_nsec(clockid_t clk_id) {
  struct timespec t;
  clock_gettime(clk_id,&t);
  return (uint64_t)t.tv_sec*1e9 + t.tv_nsec;
}

static char *pp_time(uint64_t time) {
  static int i=0;
  static char buffer[10][128];
  i = (i+1) % 10;
  if(time < 1000000) {
    snprintf(buffer[i], 128, "%4.0f us", ((double)time)/1000);
  } else if(time < 10000000) {
    snprintf(buffer[i], 128, "%4.2f ms", ((double)time)/1000000);
  } else if(time < 100000000) {
    snprintf(buffer[i], 128, "%4.1f ms", ((double)time)/1000000);
  } else if(time < 1000000000) {
    snprintf(buffer[i], 128, "%4.0f ms", ((double)time)/1000000);
  } else {
    snprintf(buffer[i], 128, "%4.2f s", ((double)time)/1000000000);
  }
  return buffer[i];
}

static char *pp_speed(uint64_t words, uint64_t time) {
  static int i=0;
  static char buffer[10][128];
  uint64_t rate = words*8/1024 / (time/1e9);
  i = (i+1) % 10;
  if(rate < 1024) {
    snprintf(buffer[i], 128, "%4" PRIu64 " KB/s", rate);
  } if(rate < 1024*1024*10) { // Less than 10000 MB/s
    snprintf(buffer[i], 128, "%4" PRIu64 " MB/s", rate/1024);
  } else {
    snprintf(buffer[i], 128, "%4.1f GB/s", (double)rate/1024/1024);
  }
  return buffer[i];
}

static char *pp_bytes(uint64_t words) {
  static int i=0;
  static char buffer[10][128];
  uint64_t bytes = words * 8;
  i = (i+1) % 10;
  if(bytes < 1024) {
    snprintf(buffer[i], 128, "%4" PRIu64 " B", bytes);
  } else if(bytes < 1024*1024) {
    snprintf(buffer[i], 128, "%4" PRIu64 " KB", bytes/1024);
  } else if(bytes < (uint64_t)1024*1024*1024*10) {
    snprintf(buffer[i], 128, "%4" PRIu64 " MB", bytes/1024/1024);
  } else {
    snprintf(buffer[i], 128, "%4.1f GB", ((double)bytes)/1024/1024/1024);
  }
  return buffer[i];
}

#define CLOCK_ID CLOCK_PROCESS_CPUTIME_ID

static void stats_timer_begin(Stats *s, enum Timer t) {
  uint64_t now = clock_gettime_nsec(CLOCK_ID);
  assert(s->n_timers < TIMER_STACK);
  s->active_timer[s->n_timers] = t;
  if(s->n_timers > 0) {
    s->timers[s->active_timer[s->n_timers-1]] += now - s->start_time;
  }

  s->n_timers++;
  s->start_time = now;
}
static void stats_timer_end(Stats *s) {
  uint64_t now = clock_gettime_nsec(CLOCK_ID);
  uint64_t t = now - s->start_time;
  assert(s->n_timers > 0);
  if(s->active_timer[s->n_timers-1] == Gen0Timer) {
    gsl_rstat_quantile_add((double)t, s->latency_50);
    gsl_rstat_quantile_add((double)t, s->latency_90);
    gsl_rstat_quantile_add((double)t, s->latency_99);
    if(t > s->nursery_time_max) {
      s->nursery_time_max = t;
    }
  }
  s->timers[s->active_timer[s->n_timers-1]] += t;
  s->start_time = now;
  s->n_timers--;
}

#endif

max_stars_repo_name: Lemmih/lhc

max_issues_repo_name: lemmih/lhc

max_forks_repo_name: lemmih/lhc


text:
/*This program evaluates the differential equations for a photon's geodesic in Minkowski spacetime.
The equations are written in the form $\frac{d(x or p)^{\alpha}}{d\lambda}=f(x^{\alpha},p^{\alpha})$.
Where $p^{\alpha}={\dot{x}}^{\alpha}$*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>          //GSL error management module

#define C 299792.458     //Speed of light in velocity units
#define DLAMBDA 0.001   //Geodesics parameter step

/*Function for solving the differential equations of Minkowski geodesics.*/
void euler1(double *x0, double *x1, double *x2, double *x3, double *p0, double *p1, double *p2, double *p3, double *lambda)
{
  /*Increment in the variables of the differential equation we want to solve*/
  double dx0, dx1, dx2, dx3, dp0, dp1, dp2, dp3;

  /*Calculation of the increments*/
  dx0 = *p0*DLAMBDA; dx1 = *p1*DLAMBDA; dx2 = *p2*DLAMBDA; dx3 = *p3*DLAMBDA;
  dp0 = 0.0; dp1 = 0.0; dp2 = 0.0; dp3 = 0.0;

  /*New values of the variables of the differential equation. Since we are using pointers, when called the routine the value of variable change.*/
  *x0 = *x0 + dx0; *x1 = *x1 + dx1; *x2 = *x2 + dx2; *x3 = *x3 + dx3;
  *p0 = *p0 + dp0; *p1 = *p1 + dp1; *p2 = *p2 + dp2; *p3 = *p3 + dp3;

  /*Increment of parameter of geodesics*/
  *lambda = *lambda + DLAMBDA;
}

int main(void)
{
  int i;            //For array manipulation

  /*Initial conditions*/
  double x0 = 0.0, x1 = 0.0, x2 = 0.0, x3 = 0.0, p1 = 0.1*C, p2 = 0.9*C, p3 = 0.0, lambda = 0.0;
  double p0 = sqrt(pow(p1,2) + pow(p2,2) + pow(p3,2))/C;

  /*Pointer to file where solution of differential equation will be saved.*/
  FILE *geodesic;
  geodesic = fopen("geodesic_solution.dat","w");

  /*Write line of initial values in file*/
  fprintf(geodesic, "%.12lf %.12lf %.12lf %.12lf %.12lf %.12lf %.12lf %.12lf %.12lf\n", lambda, x0, x1, x2, x3, p0, p1, p2, p3);

  /*Solution of the differential equation*/
  for(i=0; i<1000; i++)
    {
      euler1(&x0, &x1, &x2, &x3, &p0, &p1, &p2, &p3, &lambda);
      fprintf(geodesic, "%12lf %12lf %12lf %12lf %12lf %12lf %12lf %12lf %12lf\n", lambda, x0, x1, x2, x3, p0, p1, p2, p3);
    }

  /** Releasing all used space in memory **/
  fclose(geodesic); //Close file storing the results

  /*GNUPLOT*/

  FILE *plot;
  plot = popen("gnuplot -persist","w");
  fprintf(plot, "set terminal x11 0\n");
  fprintf(plot, "set multiplot layout 1,3\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:2 not\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:3 not\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:4 not\n");
  fprintf(plot, "unset multiplot\n");
  fprintf(plot, "reset\n");
  fprintf(plot, "set terminal x11 1\n");
  fprintf(plot, "set multiplot layout 1,3\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:6 not\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:7 not\n");
  fprintf(plot, "plot 'geodesic_solution.dat' using 1:8 not\n");
  fprintf(plot, "unset multiplot\n");
  fprintf(plot, "set terminal x11 2\n");
  fprintf(plot, "splot 'geodesic_solution.dat' using 2:3:1\n");

  pclose(plot);
  //system("rm geodesic_solution.dat");
}


text:
static char help[] = "Solve a 2-variable polynomial optimization problem.\n"
"Implements an objective function and its gradient (the residual).  The \n"
"Jacobian (Hessian) is not implemented.  Usage is either objective-only\n"
"    ./cartoon -snes_[fd|mf] -snes_fd_function\n"
"or with gradient\n"
"    ./cartoon -snes_[fd|mf]\n\n";

#include <petsc.h>

PetscErrorCode FormObjective(SNES snes, Vec x, PetscReal *Phi, void *ctx) {
    PetscErrorCode  ierr;
    const PetscReal *ax;

    ierr = VecGetArrayRead(x,&ax); CHKERRQ(ierr);
    *Phi = 0.25 * (pow(ax[0],4.0) + pow(ax[1],4.0)) - 2.0 * ax[0] + 2.0 * ax[1];
    ierr = VecRestoreArrayRead(x,&ax); CHKERRQ(ierr);
    return 0;
}

// residual F = grad Phi
PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void *ctx) {
    PetscErrorCode  ierr;
    const PetscReal *ax;
    PetscReal       *aF;

    ierr = VecGetArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecGetArray(F,&aF);CHKERRQ(ierr);
    aF[0] = ax[0]*ax[0]*ax[0] - 2.0;
    aF[1] = ax[1]*ax[1]*ax[1] + 2.0;
    ierr = VecRestoreArrayRead(x,&ax);CHKERRQ(ierr);
    ierr = VecRestoreArray(F,&aF);CHKERRQ(ierr);
    return 0;
}

int main(int argc,char **argv) {
    PetscErrorCode ierr;
    SNES           snes;
    Vec            x, r, xexact;
    const PetscInt ix[2] = {0,1};
    PetscReal      iv[2], err;

    PetscInitialize(&argc,&argv,NULL,help);

    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,2); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);

    iv[0] = 0.5;   iv[1] = 0.5;    // initial iterate corresponds to nonsingular Hessian
    ierr = VecSetValues(x,2,ix,iv,INSERT_VALUES); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(x); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(x); CHKERRQ(ierr);

    ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
    ierr = SNESSetObjective(snes,FormObjective,NULL); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&r); CHKERRQ(ierr);
    ierr = SNESSetFunction(snes,r,FormFunction,NULL); CHKERRQ(ierr);
    ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
    ierr = SNESSolve(snes,NULL,x); CHKERRQ(ierr);

    iv[0] = pow(2.0,1.0/3.0);  iv[1] = - iv[0];  // exact soln
    ierr = VecDuplicate(x,&xexact); CHKERRQ(ierr);
    ierr = VecSetValues(xexact,2,ix,iv,ADD_VALUES); CHKERRQ(ierr); // note ADD_VALUES
    ierr = VecAssemblyBegin(xexact); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(xexact); CHKERRQ(ierr);
    ierr = VecAXPY(x,-1.0,xexact); CHKERRQ(ierr);    // x <-- x + (-1.0) xexact
    ierr = VecNorm(x,NORM_INFINITY,&err); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"numerical error |x-x_exact|_inf = %g\n",err); CHKERRQ(ierr);

    VecDestroy(&x);  VecDestroy(&r);  VecDestroy(&xexact);
    SNESDestroy(&snes);
    PetscFinalize();
    return 0;
}