text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>

double f (double px, void * par) {

            double ep=-.9, ed=.0, tsp=1.63, tpd=1.13, tpp=.2;

            struct my_f_params {  double x; double y; double z; double v;};
            struct my_f_params * params  = (struct my_f_params *)par;

             double T2 = (params->x);
             double es = (params->y);
             double eF = (params->z);
             double py = (params->v);

                    double xi,argument;
                    double sx=2.*sin(px/2.),sy=2.*sin(py/2.);

                    double data[] = {   ed ,  0.   ,  tpd*sx   ,  -tpd*sy,
                                    0., es    ,  tsp*sx   ,   tsp*sy,
                                tpd*sx, tsp*sx,  ep       ,  -tpp*sx*sy,
                               -tpd*sy, tsp*sy, -tpp*sx*sy,   ep };

               gsl_matrix_view m = gsl_matrix_view_array (data, 4, 4);
               gsl_vector *eval = gsl_vector_alloc (4);
               gsl_matrix *evec = gsl_matrix_alloc (4, 4);
               gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (4);
               gsl_eigen_symmv (&m.matrix, eval, evec, w);
               gsl_eigen_symmv_free (w);
               gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);

                     double e = gsl_vector_get (eval, 2);
                     gsl_vector_view evec_i   = gsl_matrix_column (evec, 2);
                     xi = gsl_vector_get(&evec_i.vector,0)*gsl_vector_get(&evec_i.vector,1);

                gsl_vector_free (eval);
                gsl_matrix_free (evec);

                double eeF=e-eF;
                if ( eeF == 0 )
                    argument=1./T2;
                else
                    argument=(tanh(eeF/T2))/eeF;
                    argument*=xi*xi;

	return argument;
}

double g(double py, void * p)
{
                struct my_f_params { double a; double b; double c; };
                struct my_f_params * params  = (struct my_f_params *)p;

                double T2 = (params->a);
                double es = (params->b);
                double eF = (params->c);

                //printf("%g %g %g\n",a,b,c);
                struct your_f_params { double x; double y; double z; double v; };

                struct your_f_params parametri = { T2, es , eF , py };

                double result, error;
                gsl_integration_workspace * w = gsl_integration_workspace_alloc (5000000);
                gsl_function F;

                F.function = &f;
                F.params = &parametri;

                gsl_integration_qags (&F, 0, 2*M_PI, 0, 1e-7, 5000000, w, &result, &error);
                gsl_integration_workspace_free (w);

return result;
}

double dvoenintegral(double T2, double es, double eF)
{
                double dvoenintegral, error;
                struct my_f_params { double a; double b; double c; };

                gsl_integration_workspace * v = gsl_integration_workspace_alloc (5000000);

                gsl_function F;
                struct my_f_params params = {T2, es , eF};

                F.function = &g;
                F.params = &params;

                gsl_integration_qags (&F, 0, 2*M_PI, 0, 1e-7, 5000000, v, &dvoenintegral, &error);
                gsl_integration_workspace_free (v);

return (dvoenintegral)/(2*M_PI*M_PI);
}

double optim(double L, void * pp)
{

                struct my_f_params { double a; double b; double c; };
                struct my_f_params * params  = (struct my_f_params *)pp;

                double Jsd = (params->a);
                double es = (params->b);
                double eF = (params->c);

        double TT=2.*exp(-L);
        double tmp=dvoenintegral(TT,es,eF);

return Jsd*tmp-1;
}

double nameri_nula( double Jsd, double es, double eF)
{
       int status;
       double Kelvin=1/11604.;
       struct my_f_params { double a; double b; double c; };
       struct my_f_params params = { Jsd, es , eF};
       int iter = 0, max_iter = 100;
       const gsl_root_fsolver_type *T;
       gsl_root_fsolver *s;
       double r = 0;
       double x_lo = 1E-6, x_hi = 7;
       gsl_function F;

       F.function = &optim;
       F.params = &params;

       T = gsl_root_fsolver_brent;
       s = gsl_root_fsolver_alloc (T);
       gsl_root_fsolver_set (s, &F, x_lo, x_hi);

       //printf ("using %s method\n", gsl_root_fsolver_name (s));

       //printf ("%5s [%9s, %9s] %9s %10s\n", "iter", "lower", "upper", "root", "err");

       do
         {
           iter++;
           status = gsl_root_fsolver_iterate (s);
           r = gsl_root_fsolver_root (s);
           x_lo = gsl_root_fsolver_x_lower (s);
           x_hi = gsl_root_fsolver_x_upper (s);
           status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001);

           //if (status == GSL_SUCCESS)
           //printf ("Converged:\n");

           //printf ("%5d [%.7f, %.7f] %.7f %+.7f\n", iter, x_lo, x_hi, r, x_hi - x_lo);
         }
       while (status == GSL_CONTINUE && iter < max_iter);

	//printf("T_c=%g\n",exp(-r)/Kelvin);

       gsl_root_fsolver_free (s);

       return (2*exp(-r)/Kelvin);
}

double granichni_tochki(double a, double b,double es)
{

    double sx=2.*sin(a/2.);
    double sy=2.*sin(b/2.);

    double ep=-.9, ed=.0, tsp=1.63, tpd=1.13, tpp=.2;

    double data[] = {   ed ,  0.   ,  tpd*sx   ,  -tpd*sy,
                            0., es    ,  tsp*sx   ,   tsp*sy,
                        tpd*sx, tsp*sx,  ep       ,  -tpp*sx*sy,
                       -tpd*sy, tsp*sy, -tpp*sx*sy,   ep };

       gsl_matrix_view m = gsl_matrix_view_array (data, 4, 4);
       gsl_vector *eval = gsl_vector_alloc (4);
       gsl_matrix *evec = gsl_matrix_alloc (4, 4);
       gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (4);
       gsl_eigen_symmv (&m.matrix, eval, evec, w);
       gsl_eigen_symmv_free (w);
       gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_VAL_ASC);

         gsl_vector_free (eval);
         gsl_matrix_free (evec);

         return ( gsl_vector_get (eval, 2) );
}

double pd(double e, double es)
{
    double ep=-.9, ed=.0, tsp=1.63, tpd=1.13, tpp=.2;
    double A,B,C,xd,eep,eed,ees,ppd;

eep=e-ep;
eed=e-ed;
ees=e-es;

A=16.*(4.*tpd*tpd*tsp*tsp + 2.*tsp*tsp*tpp*eed-2.*tpd*tpd*tpp*ees-tpp*tpp*eed*ees);
B=-4.*eep*(tsp*tsp*eed+tpd*tpd*ees);
C=ees*eep*eep*eed;
xd=(-B+sqrt(B*B-A*C))/A;
//printf("xd=%g\n",xd);
ppd=2*asin(sqrt(xd));
//printf("xd=%g ppd=%g e=%g es=%g\n",xd,ppd,e,es);
return ppd;
}

double fermi(double px, void * p)
{
                struct my_f_params { double a; double b; };
                struct my_f_params * params  = (struct my_f_params *)p;

                double eF = (params->a);
                double es = (params->b);

    double A,B,C,eep,eed,ees,e,F;
    double ep=-.9, ed=.0, tsp=1.63, tpd=1.13, tpp=.2;

    e=eF;
eep=e-ep;
eed=e-ed;
ees=e-es;

A=16.*(4.*tpd*tpd*tsp*tsp + 2.*tsp*tsp*tpp*eed-2.*tpd*tpd*tpp*ees-tpp*tpp*eed*ees);
B=-4.*eep*(tsp*tsp*eed+tpd*tpd*ees);
C=ees*eep*eep*eed;

double x=sin(px/2.)*sin(px/2.);
F=(px-2*asin(sqrt(- (C+B*x)/(A*x+B)  )));

return F;
}

double ffactor( double eF , double es )
{
                double res, error;
                struct my_f_params { double a; double b; };

                gsl_integration_workspace * v = gsl_integration_workspace_alloc (1000);

                gsl_function F;
                struct my_f_params params = {eF , es};

                F.function = &fermi;
                F.params = &params;

                gsl_integration_qag (&F, pd(eF,es) , M_PI , 0, 1e-7, 1000, 6 , v, &res, &error);
                gsl_integration_workspace_free (v);

return  (8*res)/(4*M_PI*M_PI) ;
}

double ffactornula( double eF , void * p )
{
                double res, error;
                double es = *(double *) p;

                struct my_f_params { double a; double b; };

                gsl_integration_workspace * v = gsl_integration_workspace_alloc (100000);

                gsl_function F;
                struct my_f_params params = {eF , es};

                F.function = &fermi;
                F.params = &params;

                gsl_integration_qags (&F, pd(eF,es) , M_PI , 0, 1e-7, 100000,  v, &res, &error);
                gsl_integration_workspace_free (v);

return  (8*res)/(4*M_PI*M_PI) - 0.66;
}

double ferminula( double evhs, double etop, double es)
{
    int status;

       int iter = 0, max_iter = 100;
       const gsl_root_fsolver_type *T;
       gsl_root_fsolver *s;
       double r = 0;
       double x_lo = evhs+1E-7, x_hi = etop-1E-7;
       gsl_function F;

       F.function = &ffactornula;
       F.params = &es;

       T = gsl_root_fsolver_brent;
       s = gsl_root_fsolver_alloc (T);
       gsl_root_fsolver_set (s, &F, x_lo, x_hi);

       //printf ("using %s method\n", gsl_root_fsolver_name (s));

       //printf ("%5s [%9s, %9s] %9s %10s\n", "iter", "lower", "upper", "root", "err");

       do
         {
           iter++;
           status = gsl_root_fsolver_iterate (s);
           r = gsl_root_fsolver_root (s);
           x_lo = gsl_root_fsolver_x_lower (s);
           x_hi = gsl_root_fsolver_x_upper (s);
           status = gsl_root_test_interval (x_lo, x_hi, 0, 0.001);

           //if (status == GSL_SUCCESS)
             //printf ("Converged:\n");

           //printf ("%5d [%.7f, %.7f] %.7f %+.7f\n", iter, x_lo, x_hi, r, x_hi - x_lo);
         }
       while (status == GSL_CONTINUE && iter < max_iter);

       gsl_root_fsolver_free (s);

       return r;
}

int main(void)
{
    double Kelvin=1/11604.;
    double es=6.0, eF=1.301531002,tsp=1.63,ep=-.9,tpd=1.13;
    double Jsd,Tc,evhs,etop,s,ss;
    //double tmp;
/*
    Jsd=1/dvoenintegral(90*Kelvin,es,eF);
    printf ("Jsd= %g\n",Jsd);
    Tc=nameri_nula(Jsd,es,eF);
    printf ("Tc= %g\n",Tc);
    evhs=granichni_tochki(0,M_PI,es);
    printf("evhs=%g\n",evhs);
    etop=granichni_tochki(M_PI,M_PI,es);
    printf("etop=%g\n",etop);
    ff=ffactor(evhs, es);
    printf("filling factor(evhs)=%g\n",ff);
    printf("%g\n",ferminula(evhs,etop,es));
*/
    Jsd=1/dvoenintegral(90*Kelvin,es,eF);

    for(es=5;es<8.5;es+=0.1)
    {
    evhs=granichni_tochki(0,M_PI,es);
    etop=granichni_tochki(M_PI,M_PI,es);
    eF=ferminula(evhs,etop,es);
    Tc=nameri_nula(Jsd,es,eF);
    //s=((es-eF)*(eF-ep))/(4*tsp*tsp);
    ss=((es-eF)*(eF-ep))/(4*tsp*tpd);
    //printf("Tc=%g es=%g eF=%g\n",Tc,es,eF);
    printf("%g %g\n",ss,log(Tc));
    //printf("%g %g\n",1/(2*(1+s)),Tc );
    }

    return 0;
}

max_stars_repo_name: zlatan/Superconductivity

max_issues_repo_name: zlatan/Superconductivity

max_forks_repo_name: zlatan/Superconductivity


text:
#ifndef UTILS_H
#define UTILS_H

#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <dirent.h>
#include <memory.h>
#include <time.h>

#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sf_psi.h>
#include <gsl/gsl_sf_gamma.h>

#include <vector>
using namespace std;

double log_sum(double log_a, double log_b);
double log_normalize(double * array, int nlen);
double log_normalize(vector<double> & vec, int nlen);
double log_subtract(double log_a, double log_b);
double log_factorial(int n, double a);
double similarity(const int* v1, const int* v2, int n);

bool   file_exists(const char * filename);
bool   dir_exists(const char * directory);

template <typename T> void free_vec_ptr(vector <T* > & v)
{
    int size = v.size();
    T* p = NULL;
    for (int i = 0; i < size; i ++)
    {
        p = v[i];
        delete [] p;
    }
    v.clear();
}

// find the max and argmax in an array
template <typename T> T max(const T * x, int n, int* argmax)
{
    *argmax = 0;
    T max_val = x[0];
    for (int i = 1; i < n; i++)
        if (x[i] > max_val)
        {
            max_val = x[i];
            *argmax = i;
        }
    return max_val;
}

// find the max and argmax in an vector
template <typename T> T max_vec(vector<T> & v, int n, int* argmax)
{
    *argmax = 0;
    T max_val = v[0];
    for (int i = 1; i < n; i++)
        if (v[i] > max_val)
        {
            max_val = v[i];
            *argmax = i;
        }
    return max_val;
}

// set a value to the entire array
template <typename T> void set_array(T * array, int size, T value)
{
    for (int i = 0; i < size; i++) array[i] = value;
}

// swap two elements in an array
template < typename T > void swap_array_element(T * array, int i, int j)
{
    if (i == j) return;
    T a = array[i];
    array[i] = array[j];
    array[j] = a;
}

// set a value to a entrie vector
template <typename T> void set_vector(vector<T> & v, T value)
{
    int size = v.size();
    for (int i = 0; i < size; i++) v[i] = value;
}

// swap two elements in vector
template < typename T > void swap_vec_element(vector<T> & v, int i, int j)
{
    if (i == j) return; // no need to swap
    T a = v[i];
    v[i] = v[j];
    v[j] = a;
}

/// gsl_wrappers
//double lgamma(double x); // linux has this
unsigned int rmultinomial(const double* p, int n, double tot_p=-1);
double rgamma(double a, double b);
double rbeta(double a, double b);
unsigned int rbernoulli(double p);
double runiform();
void rshuffle (void* base, size_t n, size_t size);
unsigned long int runiform_int(unsigned long int n);

#endif

ext: h

max_stars_repo_name: marziehf/Topic-Sensitive-Embeddings

max_issues_repo_name: osmanbaskaya/mapping-impact

max_forks_repo_name: osmanbaskaya/mapping-impact


text:
#pragma once

#include <gsl/gsl_assert>
#include <src/util/Logging.h>

#ifndef NDEBUG
    //#define RNR_EXPECTS(x) Expects(x)
    //#define RNR_ENSURES(x) Ensures(x)

/**
 * specifies a precondition that is expected to be true
 * terminates the program if the condition is not met and the program is compiled in debug mode
 */
#define RNR_EXPECTS(x) \
  do { \
    if (!(x)) { \
        LOG(ERROR) << "assertion failed: " << #x; \
        std::terminate(); \
    } \
  } while (0)

/**
* specifies a postcondition that is expected to be made true by the code block above
* terminates the program if the condition is not met and the program is compiled in debug mode
*/
#define RNR_ENSURES(x) \
  do { \
    if (!(x)) { \
        LOG(ERROR) << "assertion failed: " << #x; \
        std::terminate(); \
    } \
  } while (0)
#else
    #define RNR_EXPECTS(x)
    #define RNR_ENSURES(x)
#endif

max_stars_repo_name: DavHau/Riner

max_issues_repo_name: DavHau/Riner

max_forks_repo_name: DavHau/Riner


text:
// @file utils.h
//
//  \date Created on: Sep 23, 2017
//  \author Gopalakrishna Hegde
//
//   Description:
//
//
//

#ifndef INC_UTILS_H_
#define INC_UTILS_H_
#include "common_types.h"
#include <cblas.h>

void PrintMat(char *name, const float *ptr, int H, int W, CBLAS_LAYOUT layout);
void RandInitF32(float *p_data, int N);
void SeqInitF32(float *p_data, int N);
void PrintTensor(const float *data, TensorDim dim);
int TensorSize(TensorDim dim);
bool TensorCompare(const float *t1, const float *t2, TensorDim dim);
#endif  // INC_UTILS_H_

max_stars_repo_name: gplhegde/convolution-flavors

max_issues_repo_name: chayitw/convolution-flavors

max_forks_repo_name: chayitw/convolution-flavors


text:
#ifndef __FFT_DECON_OPERATOR_H__
#define __FFT_DECON_OPERATOR_H__
#include <string>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include "mspass/seismic/TimeWindow.h"
#include "mspass/utility/Metadata.h"
#include "mspass/seismic/CoreTimeSeries.h"
#include "mspass/algorithms/deconvolution/ComplexArray.h"
namespace mspass::algorithms::deconvolution{
/*! \brief Object to hold components needed in all fft based decon algorithms.

The fft based algorithms implemented here us the GNU Scientific Library
prime factorization fft algorithm.  Those methods require initialization
given length of the fft to load and store the factorization data.  This
object holds these for all such methods and recomputes them only when
needed for efficiency.  */
class FFTDeconOperator
{
public:
    FFTDeconOperator();
    FFTDeconOperator(const mspass::utility::Metadata& md);
    FFTDeconOperator(const FFTDeconOperator& parent);
    ~FFTDeconOperator();
    FFTDeconOperator& operator=(const FFTDeconOperator& parent);
    void changeparameter(const mspass::utility::Metadata& md);
    void change_size(const int nfft_new);
    void change_shift(const int shift) {
        sample_shift=shift;
    };
    int get_size(){return nfft;};
    int get_shift(){return sample_shift;};
    int operator_size() {
        return static_cast<int>(nfft);
    };
    int operator_shift() {
        return sample_shift;
    };
    double df(const double dt){
        double period;
        period=static_cast<double>(nfft)*dt;
        return 1.0/period;
    };
    /*! \brief Return inverse wavelet for Fourier methods.

    This is a helper to be used by the inverse_wavelet method for all
    Fourier based deconvolution methods.   It avoids repetitious code that
    would be required otherwise.  inverse_wavelet methods are only
    wrappers for this generic method.  See documentation for inverse_wavelet
    for description of tshift and t0parent.  */
    mspass::seismic::CoreTimeSeries FourierInverse(const ComplexArray& winv, const ComplexArray& sw,
   	const double dt, const double t0parent);

protected:
    int nfft;
    int sample_shift;
    gsl_fft_complex_wavetable *wavetable;
    gsl_fft_complex_workspace *workspace;
    ComplexArray winv;
};

/* This helper is best referenced here */

/*! \brief Circular buffer procedure.

 In the Fourier world circular vectors are an important thing
 to deal with because the fft is intimately connected with circlular
 things.   This routine can be used, for example, to time shift the
 time domain version of a signal after it was processed with an fft.

 \param d - is the input vector to be shifted
 \param i0 - is the wrap point.   On exit sample i0 of d will be sample 0.
   (Warning - this is C convention sample number)
*/

std::vector<double>  circular_shift(const std::vector<double>& d,const int i0);
/*! Derive fft length from a time window.

All deconvlution methods using an fft need to define nfft based on the
length of the working time series.   This procedure returns the size from
an input window and sample interval. */
int ComputeFFTLength(const mspass::seismic::TimeWindow w, const double dt);
/*! Derive fft length using parameters in a metadata object.

This procedure is basically a higher level version of the function of
the same name with a time window and sample interval argument.
This procedure extracts these using three parameter keys to extract
the real numbers form md:  deconvolution_data_window_start,
decon_window_end, and target dt. */
int ComputeFFTLength(const mspass::utility::Metadata& md);
/*! Returns next power of 2 larger than n.
 *
 * Some FFT implementations require the size of the input data vector be a power
 * of 2.  This routine can be used to define a buffer size satisfying that constraint.*/
extern "C" {
    unsigned int nextPowerOf2(unsigned int n);
}
}
#endif

max_stars_repo_name: seisman/mspass

max_issues_repo_name: seisman/mspass

max_forks_repo_name: seisman/mspass


text:
/**
 * @file coroutine/net.h
 * @brief Async I/O operation support with system socket functions
 * @author github.com/luncliff (luncliff@gmail.com)
 * @copyright CC BY 4.0
 */
#pragma once
#ifndef COROUTINE_NET_IO_H
#define COROUTINE_NET_IO_H
#include <gsl/gsl>

#include <coroutine/return.h>

/**
 * @defgroup Network
 * Helper types to apply `co_await` for socket operations + Name resolution utilities
 */

#if __has_include(<WinSock2.h>) // use winsock
#include <WS2tcpip.h>
#include <WinSock2.h>
#include <ws2def.h>

/** @brief indicates current import: using Windows Socket API? */
/** @brief indicates current import: using <netinet/*.h>? */
static constexpr bool is_winsock = true;
static constexpr bool is_netinet = false;

using io_control_block = OVERLAPPED;

#elif __has_include(<netinet/in.h>) // use netinet
#include <fcntl.h>
#include <netdb.h>
#include <netinet/in.h>
#include <netinet/tcp.h>
#include <sys/socket.h>
#include <unistd.h>

static constexpr bool is_winsock = false;
static constexpr bool is_netinet = true;

/**
 * @brief Follow the definition of Windows `OVERLAPPED`
 * @see https://docs.microsoft.com/en-us/windows/win32/sync/synchronization-and-overlapped-input-and-output
 * @see https://docs.microsoft.com/en-us/windows/win32/api/minwinbase/ns-minwinbase-overlapped
 * @ingroup Network
 */
struct io_control_block {
    uint64_t internal;      // uint32_t errc, int32_t flag
    uint64_t internal_high; // int64_t len, socklen_t addrlen
    union {
        struct {
            int32_t offset;
            int32_t offset_high;
        };
        void* ptr; // sockaddr* addr;
    };
    int64_t handle; // int64_t sd;
};

#endif // winsock || netinet

namespace coro {

/**
 * @brief This is simply a view to storage. Be aware that it doesn't have ownership
 * @ingroup Network
 */
using io_buffer_t = gsl::span<std::byte>;
static_assert(sizeof(io_buffer_t) <= sizeof(void*) * 2);

/**
 * @brief A struct to describe "1 I/O request" to system API.
 * @ingroup Network
 * When I/O request is submitted, an I/O task becomes 1 coroutine handle
 */
class io_work_t : public io_control_block {
  public:
    coroutine_handle<void> task{};
    io_buffer_t buffer{};

  protected:
    /**
     * @see await_ready
     * @return true  The given socket can be use for non-blocking operations
     * @return false For Windows, the return is always `false`
     */
    bool ready() const noexcept;

  public:
    /**
     * @brief Multiple retrieving won't be a matter
     * @return uint32_t error code from the system
     */
    uint32_t error() const noexcept;
};
static_assert(sizeof(io_work_t) <= 56);

/**
 * @brief Awaitable type to perform `sendto` I/O request
 * @see sendto
 * @see WSASendTo
 * @ingroup Network
 */
class io_send_to final : public io_work_t {
  private:
    /**
     * @brief makes an I/O request with given context(`coroutine_handle<void>`)
     * @throw std::system_error
     */
    void suspend(coroutine_handle<void> t) noexcept(false);
    /**
     * @brief Fetch I/O result/error
     * @return int64_t return of `sendto`
     *
     * This function must be used through `co_await`.
     * Multiple invoke of this will lead to malfunction.
     */
    int64_t resume() noexcept;

  public:
    bool await_ready() const noexcept {
        return this->ready();
    }
    /**
     * @throw std::system_error
     */
    void await_suspend(coroutine_handle<void> t) noexcept(false) {
        return this->suspend(t);
    }
    int64_t await_resume() noexcept {
        return this->resume();
    }
};
static_assert(sizeof(io_send_to) == sizeof(io_work_t));

/**
 * @brief Awaitable type to perform `recvfrom` I/O request
 * @see recvfrom
 * @see WSARecvFrom
 * @ingroup Network
 */
class io_recv_from final : public io_work_t {
  private:
    /**
     * @brief makes an I/O request with given context(`coroutine_handle<void>`)
     * @throw std::system_error
     */
    void suspend(coroutine_handle<void> t) noexcept(false);
    /**
     * @brief Fetch I/O result/error
     * @return int64_t return of `recvfrom`
     *
     * This function must be used through `co_await`.
     * Multiple invoke of this will lead to malfunction.
     */
    int64_t resume() noexcept;

  public:
    bool await_ready() const noexcept {
        return this->ready();
    }
    /**
     * @throw std::system_error
     */
    void await_suspend(coroutine_handle<void> t) noexcept(false) {
        return this->suspend(t);
    }
    int64_t await_resume() noexcept {
        return this->resume();
    }
};
static_assert(sizeof(io_recv_from) == sizeof(io_work_t));

/**
 * @brief Awaitable type to perform `send` I/O request
 * @see send
 * @see WSASend
 * @ingroup Network
 */
class io_send final : public io_work_t {
  private:
    /**
     * @brief makes an I/O request with given context(`coroutine_handle<void>`)
     * @throw std::system_error
     */
    void suspend(coroutine_handle<void> t) noexcept(false);
    /**
     * @brief Fetch I/O result/error
     * @return int64_t return of `send`
     *
     * This function must be used through `co_await`.
     * Multiple invoke of this will lead to malfunction.
     */
    int64_t resume() noexcept;

  public:
    bool await_ready() const noexcept {
        return this->ready();
    }
    void await_suspend(coroutine_handle<void> t) noexcept(false) {
        return this->suspend(t);
    }
    int64_t await_resume() noexcept {
        return this->resume();
    }
};
static_assert(sizeof(io_send) == sizeof(io_work_t));

/**
 * @brief Awaitable type to perform `recv` I/O request
 * @see recv
 * @see WSARecv
 * @ingroup Network
 */
class io_recv final : public io_work_t {
  private:
    /**
     * @brief makes an I/O request with given context(`coroutine_handle<void>`)
     * @throw std::system_error
     */
    void suspend(coroutine_handle<void> t) noexcept(false);

    /**
     * @brief Fetch I/O result/error
     * @return int64_t return of `recv`
     *
     * This function must be used through `co_await`.
     * Multiple invoke of this will lead to malfunction.
     */
    int64_t resume() noexcept;

  public:
    bool await_ready() const noexcept {
        return this->ready();
    }
    void await_suspend(coroutine_handle<void> t) noexcept(false) {
        return this->suspend(t);
    }
    int64_t await_resume() noexcept {
        return this->resume();
    }
};
static_assert(sizeof(io_recv) == sizeof(io_work_t));

/**
 * @brief Constructs `io_send_to` awaitable with the given parameters
 * @param sd
 * @param remote
 * @param buf
 * @param work
 * @return io_send_to&
 *
 * @ingroup Network
 */
auto send_to(uint64_t sd, const sockaddr_in& remote, io_buffer_t buf,
             io_work_t& work) noexcept(false) -> io_send_to&;

/**
 * @brief Constructs `io_send_to` awaitable with the given parameters
 * @param sd
 * @param remote
 * @param buf
 * @param work
 * @return io_send_to&
 *
 * @ingroup Network
 */
auto send_to(uint64_t sd, const sockaddr_in6& remote, io_buffer_t buf,
             io_work_t& work) noexcept(false) -> io_send_to&;

/**
 * @brief Constructs `io_recv_from` awaitable with the given parameters
 * @param sd
 * @param remote
 * @param buf
 * @param work
 * @return io_recv_from&
 *
 * @ingroup Network
 */
auto recv_from(uint64_t sd, sockaddr_in& remote, io_buffer_t buf,
               io_work_t& work) noexcept(false) -> io_recv_from&;

/**
 * @brief Constructs `io_recv_from` awaitable with the given parameters
 * @param sd
 * @param remote
 * @param buf
 * @param work
 * @return io_recv_from&
 *
 * @ingroup Network
 */
auto recv_from(uint64_t sd, sockaddr_in6& remote, io_buffer_t buf,
               io_work_t& work) noexcept(false) -> io_recv_from&;

/**
 * @brief Constructs `io_send` awaitable with the given parameters
 * @param sd
 * @param buf
 * @param flag
 * @param work
 * @return io_send&
 *
 * @ingroup Network
 */
auto send_stream(uint64_t sd, io_buffer_t buf, uint32_t flag,
                 io_work_t& work) noexcept(false) -> io_send&;

/**
 * @brief Constructs `io_recv` awaitable with the given parameters
 * @param sd
 * @param buf
 * @param flag
 * @param work
 * @return io_recv&
 *
 * @ingroup Network
 */
auto recv_stream(uint64_t sd, io_buffer_t buf, uint32_t flag,
                 io_work_t& work) noexcept(false) -> io_recv&;

/**
 * @brief Poll internal I/O works and invoke user callback
 * @param nano timeout in nanoseconds
 * @throw std::system_error
 *
 * @ingroup Network
 */
void poll_net_tasks(uint64_t nano) noexcept(false);

/**
 * @brief Thin wrapper of `getaddrinfo` for IPv4
 *
 * @param hint
 * @param host
 * @param serv
 * @param output
 * @return uint32_t Error code from the `getaddrinfo` that can be the argument of `gai_strerror`
 * @see getaddrinfo
 * @see gai_strerror
 *
 * @ingroup Network
 */
uint32_t get_address(const addrinfo& hint, //
                     gsl::czstring host, gsl::czstring serv,
                     gsl::span<sockaddr_in> output) noexcept;

/**
 * @brief Thin wrapper of `getaddrinfo` for IPv6
 *
 * @param hint
 * @param host
 * @param serv
 * @param output
 * @return uint32_t Error code from the `getaddrinfo` that can be the argument of `gai_strerror`
 * @see getaddrinfo
 * @see gai_strerror
 *
 * @ingroup Network
 */
uint32_t get_address(const addrinfo& hint, //
                     gsl::czstring host, gsl::czstring serv,
                     gsl::span<sockaddr_in6> output) noexcept;

/**
 * @brief Thin wrapper of `getnameinfo`
 *
 * @param addr
 * @param name
 * @param serv can be `nullptr`
 * @param flags
 * @return uint32_t EAI_AGAIN ...
 * @see getnameinfo
 *
 * @ingroup Network
 */
uint32_t get_name(const sockaddr_in& addr, //
                  gsl::basic_zstring<char, NI_MAXHOST> name, gsl::basic_zstring<char, NI_MAXSERV> serv,
                  int32_t flags = NI_NUMERICHOST | NI_NUMERICSERV) noexcept;

/**
 * @brief Thin wrapper of `getnameinfo`
 * @param addr
 * @param name
 * @param serv can be `nullptr`
 * @param flags
 * @return uint32_t EAI_AGAIN ...
 * @see getnameinfo
 *
 * @ingroup Network
 */
uint32_t get_name(const sockaddr_in6& addr, //
                  gsl::basic_zstring<char, NI_MAXHOST> name, gsl::basic_zstring<char, NI_MAXSERV> serv,
                  int32_t flags = NI_NUMERICHOST | NI_NUMERICSERV) noexcept;

} // namespace coro

#endif // COROUTINE_NET_IO_H

max_stars_repo_name: dmitrykobets-msft/coroutine

max_issues_repo_name: dmitrykobets-msft/coroutine

max_forks_repo_name: dmitrykobets-msft/coroutine


text:
/*! @copyright (c) 2017 King Abdullah University of Science and
 *                      Technology (KAUST). All rights reserved.
 *
 * STARS-H is a software package, provided by King Abdullah
 *             University of Science and Technology (KAUST)
 *
 * @file testing/rndtiled.c
 * @version 1.3.0
 * @author Aleksandr Mikhalev
 * @date 2017-11-07
 * */

#ifdef MKL
    #include <mkl.h>
#else
    #include <cblas.h>
    #include <lapacke.h>
#endif
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <starsh.h>
#include <starsh-randtlr.h>

int main(int argc, char **argv)
{
    if(argc != 6)
    {
        printf("%d arguments provided, but 5 are needed\n", argc-1);
        printf("randtlr N NB decay maxrank tol\n");
        return 1;
    }
    int N = atoi(argv[1]);
    int block_size = atoi(argv[2]);
    double decay = atof(argv[3]);
    int maxrank = atoi(argv[4]);
    double tol = atof(argv[5]);
    int onfly = 0;
    char symm = 'N', dtype = 'd';
    STARSH_int shape[2] = {N, N};
    int info;
    srand(0);
    // Init STARS-H
    info = starsh_init();
    if(info != 0)
        return info;
    // Generate problem by random matrices
    STARSH_randtlr *data;
    STARSH_kernel *kernel;
    info = starsh_application((void **)&data, &kernel, N, dtype,
            STARSH_RANDTLR, STARSH_RANDTLR_KERNEL1,
            STARSH_RANDTLR_NB, block_size,
            STARSH_RANDTLR_DECAY, decay, STARSH_RANDTLR_DIAG, (double)N, 0);
    if(info != 0)
    {
        printf("Problem was NOT generated (wrong parameters)\n");
        return info;
    }
    // Init problem with given data and kernel and print short info
    STARSH_problem *P;
    info = starsh_problem_new(&P, 2, shape, symm, dtype, data, data, kernel,
            "Randomly generated tiled blr-matrix");
    if(info != 0)
        return info;
    starsh_problem_info(P);
    // Init plain clusterization and print info
    STARSH_cluster *C;
    info = starsh_cluster_new_plain(&C, data, N, block_size);
    if(info != 0)
        return info;
    starsh_cluster_info(C);
    // Init tlr division into admissible blocks and print short info
    STARSH_blrf *F;
    info = starsh_blrf_new_tlr(&F, P, symm, C, C);
    if(info != 0)
        return info;
    starsh_blrf_info(F);
    // Approximate each admissible block
    STARSH_blrm *M;
    double time1 = omp_get_wtime();
    info = starsh_blrm_approximate(&M, F, maxrank, tol, onfly);
    if(info != 0)
        return info;
    time1 = omp_get_wtime()-time1;
    // Print info about updated format and approximation
    starsh_blrf_info(F);
    starsh_blrm_info(M);
    printf("TIME TO APPROXIMATE: %e secs\n", time1);
    // Measure approximation error
    time1 = omp_get_wtime();
    double rel_err = starsh_blrm__dfe_omp(M);
    time1 = omp_get_wtime()-time1;
    printf("TIME TO MEASURE ERROR: %e secs\nRELATIVE ERROR: %e\n",
            time1, rel_err);
    if(rel_err/tol > 10.)
    {
        printf("Resulting relative error is too big\n");
        return 1;
    }
    // Free block low-rank matrix
    starsh_blrm_free(M);
    // Free block low-rank format
    starsh_blrf_free(F);
    // Free clusterization info
    starsh_cluster_free(C);
    // Free STARS_Problem instance
    starsh_problem_free(P);
    // Free randomly generated data
    starsh_randtlr_free(data);
    return 0;
}

ext: c

max_stars_repo_name: wawando/stars-h

max_issues_repo_name: wawando/stars-h

max_forks_repo_name: wawando/stars-h


text:
/**
 *
 * @file qwrapper_dlantr.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 Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:57 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

void
CORE_dlantr_quark(Quark *quark);
void
CORE_dlantr_f1_quark(Quark *quark);

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlantr(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                       const double *A, int LDA, int szeA,
                       int szeW, double *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;
    QUARK_Insert_Task(quark, CORE_dlantr_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(PLASMA_enum),                &diag,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(double)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(double)*szeW,                 NULL,          SCRATCH,
        sizeof(double),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlantr_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, PLASMA_enum diag, int M, int N,
                          const double *A, int LDA, int szeA,
                          int szeW, double *result,
                          double *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANGE;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_dlantr_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(double)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_dlantr_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(PLASMA_enum),                &diag,  VALUE,
            sizeof(int),                        &M,     VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(double)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT,
            sizeof(double)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlantr_quark = PCORE_dlantr_quark
#define CORE_dlantr_quark PCORE_dlantr_quark
#endif
void CORE_dlantr_quark(Quark *quark)
{
    double *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    double *A;
    int LDA;
    double *work;

    quark_unpack_args_9(quark, norm, uplo, diag, M, N, A, LDA, work, normA);
    *normA = LAPACKE_dlantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlantr_f1_quark = PCORE_dlantr_f1_quark
#define CORE_dlantr_f1_quark PCORE_dlantr_f1_quark
#endif
void CORE_dlantr_f1_quark(Quark *quark)
{
    double *normA;
    PLASMA_enum norm, uplo, diag;
    int M;
    int N;
    double *A;
    int LDA;
    double *work;
    double *fake;

    quark_unpack_args_10(quark, norm, uplo, diag, M, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_dlantr_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        lapack_const(diag),
        M, N, A, LDA, work);
}


text:
#ifndef S3D_DISPARITY_VIEWER_DEPTH_CONVERTER_H
#define S3D_DISPARITY_VIEWER_DEPTH_CONVERTER_H

#include <s3d/utilities/eigen.h>

#include <vector>

#include <gsl/gsl>

namespace s3d {

struct ViewerContext;

class ViewerDepthConverter {
 public:
  using Point = Eigen::Vector2d;
  using Pointf = Eigen::Vector2f;

  explicit ViewerDepthConverter(gsl::not_null<ViewerContext*> context);

  std::vector<float> computePerceivedDepth(const std::vector<float>& disparitiesPercent);

  // horizontal position and depth in meters
  std::vector<Pointf> computeDepthPositions(const std::vector<Point>& imagePoints,
                                            const std::vector<float>& disparities);

  float computePerceivedDepth(float disparityPercent);

  // 0 in the middle of the image
  float computeHorizontalPosition(float imageX);

  std::vector<float> computeHorizontalPositions(const std::vector<Point>& imagePoints);

  void setViewerContext(gsl::not_null<ViewerContext*> viewerContext);

 private:
  ViewerContext* viewerContext_;
};
}  // namespace s3d

#endif  // S3D_DISPARITY_VIEWER_DEPTH_CONVERTER_H


text:
/* rstat/rquantile.c
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * 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 <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_rstat.h>

/*
 * Running quantile calculation based on the paper
 *
 * [1] R. Jain and I. Chlamtac, "The P^2 algorithm for dynamic
 *     calculation of quantiles and histograms without storing
 *     observations", Communications of the ACM, October 1985
 */

static double calc_psq(const double qp1, const double q, const double qm1,
                       const double d, const double np1, const double n, const double nm1);

gsl_rstat_quantile_workspace *
gsl_rstat_quantile_alloc(const double p)
{
  gsl_rstat_quantile_workspace *w;

  w = calloc(1, sizeof(gsl_rstat_quantile_workspace));
  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->p = p;

  gsl_rstat_quantile_reset(w);

  return w;
} /* gsl_rstat_quantile_alloc() */

void
gsl_rstat_quantile_free(gsl_rstat_quantile_workspace *w)
{
  free(w);
} /* gsl_rstat_quantile_free() */

int
gsl_rstat_quantile_reset(gsl_rstat_quantile_workspace *w)
{
  const double p = w->p;
  size_t i;

  /* initialize positions n */
  for (i = 0; i < 5; ++i)
    w->npos[i] = i + 1;

  /* initialize n' */
  w->np[0] = 1.0;
  w->np[1] = 1.0 + 2.0 * p;
  w->np[2] = 1.0 + 4.0 * p;
  w->np[3] = 3.0 + 2.0 * p;
  w->np[4] = 5.0;

  /* initialize dn' */
  w->dnp[0] = 0.0;
  w->dnp[1] = 0.5 * p;
  w->dnp[2] = p;
  w->dnp[3] = 0.5 * (1.0 + p);
  w->dnp[4] = 1.0;

  w->n = 0;

  return GSL_SUCCESS;
}

int
gsl_rstat_quantile_add(const double x, gsl_rstat_quantile_workspace *w)
{
  if (w->n < 5)
    {
      w->q[w->n] = x;
    }
  else
    {
      int i;
      int k = -1;

      if (w->n == 5)
        {
          /* initialization: sort the first five heights */
          gsl_sort(w->q, 1, w->n);
        }

      /* step B1: find k such that q_k <= x < q_{k+1} */
      if (x < w->q[0])
        {
          w->q[0] = x;
          k = 0;
        }
      else if (x >= w->q[4])
        {
          w->q[4] = x;
          k = 3;
        }
      else
        {
          for (i = 0; i <= 3; ++i)
            {
              if (w->q[i] <= x && x < w->q[i + 1])
                {
                  k = i;
                  break;
                }
            }
        }

      if (k < 0)
        {
          /* we could get here if x is nan */
          GSL_ERROR ("invalid input argument x", GSL_EINVAL);
        }

      /* step B2(a): update n_i */
      for (i = k + 1; i <= 4; ++i)
        ++(w->npos[i]);

      /* step B2(b): update n_i' */
      for (i = 0; i < 5; ++i)
        w->np[i] += w->dnp[i];

      /* step B3: update heights */
      for (i = 1; i <= 3; ++i)
        {
          double ni = (double) w->npos[i];
          double d = w->np[i] - ni;

          if ((d >= 1.0 && (w->npos[i + 1] - w->npos[i] > 1)) ||
              (d <= -1.0 && (w->npos[i - 1] - w->npos[i] < -1)))
            {
              int dsign = (d > 0.0) ? 1 : -1;
              double qp1 = w->q[i + 1];
              double qi = w->q[i];
              double qm1 = w->q[i - 1];
              double np1 = (double) w->npos[i + 1];
              double nm1 = (double) w->npos[i - 1];
              double qp = calc_psq(qp1, qi, qm1, (double) dsign,
                                   np1, ni, nm1);

              if (qm1 < qp && qp < qp1)
                w->q[i] = qp;
              else
                {
                  /* use linear formula */
                  w->q[i] += dsign * (w->q[i + dsign] - qi) / ((double) w->npos[i + dsign] - ni);
                }

              w->npos[i] += dsign;
            }
        }
    }

  ++(w->n);

  return GSL_SUCCESS;
} /* gsl_rstat_quantile_add() */

double
gsl_rstat_quantile_get(gsl_rstat_quantile_workspace *w)
{
  if (w->n >= 5)
    {
      return w->q[2];
    }
  else
    {
      /* not yet initialized */
      gsl_sort(w->q, 1, w->n);
      return gsl_stats_quantile_from_sorted_data(w->q, 1, w->n, w->p);
    }
} /* gsl_rstat_quantile_get() */

static double
calc_psq(const double qp1, const double q, const double qm1,
         const double d, const double np1, const double n, const double nm1)
{
  double outer = d / (np1 - nm1);
  double inner_left = (n - nm1 + d) * (qp1 - q) / (np1 - n);
  double inner_right = (np1 - n - d) * (q - qm1) / (n - nm1);

  return q + outer * (inner_left + inner_right);
} /* calc_psq() */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
Matrix multiplication.

Based on the amazing:
https://github.com/HandsOnOpenCL/Exercises-Solutions/tree/a908ac3f0fadede29f2735eb1264b0db7f4311a0/Solutions/Exercise08

The most basic / useful application where OpenCL might be faster than CPU.

TODO: make a SERIOUS matrix implementation. Also compare with existing SERIOUS CPU and GPU implementations:

- http://stackoverflow.com/questions/1907557/optimized-matrix-multiplication-in-c
- http://stackoverflow.com/questions/12289235/simple-and-fast-matrix-vector-multiplication-in-c-c
- https://www.quora.com/What-is-the-best-way-to-multiply-two-matrices-in-C++
- http://stackoverflow.com/questions/25900312/optimizing-batched-matrix-multiplication-opencl-code

Serious CPU implementation means it considers:

- caching
- SIMD

Articles:

- http://www.netlib.org/utk/papers/autoblock/node2.html
- http://codereview.stackexchange.com/questions/101144/simd-matrix-multiplication
*/

#include "common.h"

#include <cblas.h>
#include <clBLAS.h>

#define NELEMS(a) sizeof(a)/sizeof(a[0])
#define UNUSED(x) (void)(x)

/* TODO how to auto tune this?
 * GCC 6 was not smart enough to use widest type automatically:
 * anything larger than the widest type just dropped perf drastically. */
#define VECTOR_NELEMS 4
#define VECTOR_SIZEOF (VECTOR_NELEMS * sizeof(F))

typedef cl_float F;
typedef F Vec __attribute__ ((vector_size(VECTOR_SIZEOF)));
/* Stuff reused across OpenCL calls.
 * It would not be fair to consider their alocation time for each call,
 * as the functions will be called millions of times.
 *
 * TODO: also factor out the kernel compile time.
 * Currently recompiled every time.
 * Implementations already cache compiled kernels, but we would need
 * to call each function multiple times.
 * */
typedef struct {
    Common common;
    /* We only factor out the allocation of those buffers.
     * Data transfer time with clEnqueueWriteBuffer must still be considered. */
    cl_mem buf_a, buf_b, buf_c;
} Cache;
typedef void (*MatMul)(const F *A, const F *B, F *C, size_t n, Cache *cache);

void cl_buf_init(Cache *cache, size_t mat_sizeof) {
    cache->buf_a = clCreateBuffer(cache->common.context, CL_MEM_READ_ONLY, mat_sizeof, NULL, NULL);
    cache->buf_b = clCreateBuffer(cache->common.context, CL_MEM_READ_ONLY, mat_sizeof, NULL, NULL);
    cache->buf_c = clCreateBuffer(cache->common.context, CL_MEM_WRITE_ONLY, mat_sizeof, NULL, NULL);
}

void cl_buf_deinit(Cache *cache) {
    clReleaseMemObject(cache->buf_a);
    clReleaseMemObject(cache->buf_b);
    clReleaseMemObject(cache->buf_c);
}

/* No, this was not created for debugging, my code is flawless from the first try. */
void mat_print(const F *A, size_t n) {
    size_t i, j;
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            printf("%f ", A[i*n+j]);
        }
        puts("");
    }
}

/* Zero a matrix. */
void mat_zero(F *A, size_t n) {
    size_t i, n2;
    n2 = n*n;
    for (i = 0; i < n2; ++i) {
        A[i] = 0.0;
    }
}

/* Initialize a random matrix. */
void mat_rand(F *A, size_t n) {
    size_t i, n2;
    n2 = n*n;
    for (i = 0; i < n2; ++i) {
        A[i] = ((float)rand()) / ((float)RAND_MAX);
    }
}

/* Initialize a random matrix. */
void mat_trans(F *A, size_t n) {
    size_t i, j, i1, i2;
    F tmp;
    for (i = 0; i < n; ++i) {
        for (j = 0; j < i; ++j) {
            i1 = i*n+j;
            i2 = j*n+i;
            tmp = A[i1];
            A[i1] = A[i2];
            A[i2] = tmp;
        }
    }
}

/* Check if two matrices are equal with given mean squared err_maxor. */
int mat_eq(const F *A, const F *B, size_t n) {
    const F err_max = 10e-3;
    F err, diff, a, b;
    size_t i, i_max;

    err = 0.0;
    i_max = n*n;
    for (i = 0; i < i_max; ++i) {
        a = A[i];
        b = B[i];
        diff = a - b;
        err += diff * diff;
    }
    return (sqrt(err) / i_max) < err_max;
}

void mat_assert_eq(const F *A, const F *B, size_t n) {
    if (!mat_eq(A, B, n)) {
        puts("");
        mat_print(A, n);
        puts("");
        mat_print(B, n);
        assert(0);
    }
}

size_t zmin(size_t x, size_t y) {
    return (x < y) ? x : y;
}

/* C = A*B, width n, naive. */
void mat_mul_cpu(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    F tmp;
    size_t i, j, k;
    UNUSED(cache);

    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            tmp = 0.0;
            for (k = 0; k < n; ++k) {
                tmp += A[i*n+k] * B[k*n+j];
            }
            C[i*n+j] = tmp;
        }
    }
}

/* Transpose matrix B to increase cache hits. */
void mat_mul_cpu_trans(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    F tmp;
    size_t i, j, k;
    UNUSED(cache);

    mat_trans((F*)B, n);
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            tmp = 0.0;
            for (k = 0; k < n; ++k) {
                tmp += A[i*n+k] * B[j*n+k];
            }
            C[i*n+j] = tmp;
        }
    }
    mat_trans((F*)B, n);
}

/* Zero vector. */
void vec_zero(Vec *vec, size_t vec_n) {
    size_t i;
    for (i = 0; i < vec_n; ++i) {
        (*vec)[i] = 0.0;
    }
}

/* Load vector from array. */
void vec_load(Vec *vec, size_t vec_n, const F* A, size_t A_base) {
    size_t i;
    for (i = 0; i < vec_n; ++i) {
        (*vec)[i] = A[A_base + i];
    }
}

/* Sum elements of the vector. */
F vec_sum(Vec vec, size_t vec_n) {
    size_t i;
    F sum;
    sum = 0.0;
    for (i = 0; i < vec_n; ++i) {
        sum += vec[i];
    }
    return sum;
}

/* Transpose matrix B to both:
 *
 * - increase cache hits
 * - simd GCC vector extensions which is made possible.
 *   by the transposition, to increase likelyhood of SIMDs.
 *
 * Note that GCC 6 O=3 is smart enough to use SIMD
 * even for the naive CPU method. However this was still way faster.
 * */
void mat_mul_cpu_trans_vec(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    F tmpf;
    size_t i, j, k, k_max, ai, bi;
    Vec tmp, a, b;
    UNUSED(cache);

    mat_trans((F*)B, n);
    k_max = (n / VECTOR_NELEMS) * VECTOR_NELEMS;
    for (i = 0; i < n; ++i) {
        for (j = 0; j < n; ++j) {
            vec_zero(&tmp, VECTOR_NELEMS);
            for (k = 0; k < k_max; k += VECTOR_NELEMS) {
                ai = i * n + k;
                bi = j * n + k;
                vec_load(&a, VECTOR_NELEMS, A, ai);
                vec_load(&b, VECTOR_NELEMS, B, bi);
                tmp += a * b;
            }
            tmpf = 0.0;
            for (; k < n; ++k) {
                tmpf += A[i*n+k] * B[j*n+k];
            }
            C[i*n+j] = vec_sum(tmp, VECTOR_NELEMS) + tmpf;
        }
    }
    mat_trans((F*)B, n);
}

/* Blocked matrix multiplication.
 * TODO slower than transpose, sometimes similar timing to naive.
 * Why do they say that this is better?
 * */
void mat_mul_cpu_block(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    size_t ib, jb, kb, i, j, k, i_max, j_max, k_max, nb ;
    F tmp;
    UNUSED(cache);

    nb = lround(sqrt(n));
    for (ib = 0; ib < n; ib += nb) {
        i_max = zmin(ib + nb, n);
        for (jb = 0; jb < n; jb += nb) {
            k_max = zmin(jb + nb, n);
            for (kb = 0; kb < n; kb += nb) {
                j_max = zmin(kb + nb, n);

                /* TODO would be cool to recurse here, but would require more offset parameters,
                 * likely similar to tose of CBLAS. */
                for (i = ib; i < i_max; ++i) {
                    for (j = kb; j < j_max; ++j) {
                        tmp = 0.0;
                        for (k = jb; k < k_max; ++k) {
                            tmp += A[i*n+k] * B[k*n+j];
                        }
                        C[i*n+j] += tmp;
                    }
                }
            }
        }
    }
}

/* The golden single thread CPU standard. */
void mat_mul_cpu_cblas(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    UNUSED(cache);
    cblas_sgemm(
        CblasRowMajor,
        CblasNoTrans,
        CblasNoTrans,
        n,
        n,
        n,
        1.0,
        A,
        n,
        B,
        n,
        0.0,
        C,
        n
    );
}

/* Simplest possible CL implementation. No speedup. */
void mat_mul_cl(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    cl_uint ncl;
    size_t global_work_size[2], mat_sizeof;

    /* Setup variables. */
    global_work_size[0] = n;
    global_work_size[1] = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    common_create_kernel_file(&cache->common, "matmul.cl", NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, sizeof(ncl), &ncl);
    clEnqueueNDRangeKernel(cache->common.command_queue, cache->common.kernel, 2, NULL, global_work_size, NULL, 0, NULL, NULL);
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

/* Cache rows in private memory. Drastic speedups expected over naive CPU. */
void mat_mul_cl_row_priv(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    char options[256];
    cl_uint ncl;
    size_t global_work_size, mat_sizeof;

    /* Setup variables. */
    global_work_size = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    snprintf(options, sizeof(options), "-DPRIV_ROW_SIZE=%ju", n);
    common_create_kernel_file(&cache->common, "matmul_row_priv.cl", options);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, sizeof(ncl), &ncl);
    clEnqueueNDRangeKernel(cache->common.command_queue, cache->common.kernel, 1, NULL, &global_work_size, NULL, 0, NULL, NULL);
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

/* Let's see if this is any different from local memory.
 * Outcome: much slower than private memory, slower than naive method. */
void mat_mul_cl_row_local(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    cl_uint ncl;
    size_t global_work_size, local_work_size, mat_sizeof;

    /* Setup variables. */
    /* Cannot be larger than 1 on this example, otherwise memory conflicts
     * will happen between work items. */
    local_work_size = 1;
    global_work_size = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    common_create_kernel_file(&cache->common, "matmul_row_local.cl", NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, n * sizeof(F), NULL);
    clSetKernelArg(cache->common.kernel, 4, sizeof(ncl), &ncl);
    clEnqueueNDRangeKernel(cache->common.command_queue, cache->common.kernel, 1, NULL, &global_work_size, &local_work_size, 0, NULL, NULL);
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

/* Like row private, but to reduce global memory accesses,
 * we copy only once per work group to local memory.
 *
 * Each work group contains a few rows of A.
 *
 * We load the first column, multiply all rows by that column, synrhconize,
 * load another column, and so on.
 *
 * This leads to a thread blockage / memory access tradeoff.
 *
 * We make work groups as large as possible to reload memory less times. */
void mat_mul_cl_row_priv_col_local(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    char options[256];
    cl_uint ncl;
    size_t global_work_size, local_work_size, mat_sizeof;

    /* Setup variables. */
    global_work_size = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    snprintf(options, sizeof(options), "-DPRIV_ROW_SIZE=%ju", n);
    common_create_kernel_file(&cache->common, "matmul_row_priv_col_local.cl", options);
    local_work_size = 0;
    clGetDeviceInfo(cache->common.device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(local_work_size), &local_work_size, NULL);
    local_work_size = zmin(local_work_size, n);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, n * sizeof(F), NULL);
    clSetKernelArg(cache->common.kernel, 4, sizeof(ncl), &ncl);
    clEnqueueNDRangeKernel(cache->common.command_queue, cache->common.kernel, 1, NULL, &global_work_size, &local_work_size, 0, NULL, NULL);
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

/* Copy as many cols from B as possibl to the local memory, only then start multiplying.
 * This leads to less memory barrier hits.
 * How many rows we copy is limited by the local memory size, ideally the entire matrix will fit. */
void mat_mul_cl_row_priv_cols_local(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    char options[256];
    cl_uint ncl, n_local_cols;
    cl_ulong local_mem_size;
    size_t col_size, global_work_size, local_work_size, mat_sizeof;

    /* Setup variables. */
    col_size = n * sizeof(F);
    global_work_size = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    snprintf(options, sizeof(options), "-DPRIV_ROW_SIZE=%ju", n);
    common_create_kernel_file(&cache->common, "matmul_row_priv_cols_local.cl", options);
    local_work_size = 0;
    clGetDeviceInfo(cache->common.device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(local_work_size), &local_work_size, NULL);
    local_work_size = zmin(local_work_size, n);
    local_mem_size = 0;
    clGetDeviceInfo(cache->common.device, CL_DEVICE_LOCAL_MEM_SIZE, sizeof(local_mem_size), &local_mem_size, NULL);
    /* TODO: can blow up without that - 1. Why?
     * It only reaches the max without it, not crosses, right?
     * So bug in the kernel? */
    n_local_cols = zmin(local_mem_size / col_size, n) - 1;
    /*puts("");*/
    /*printf("max memory %llu\n", (unsigned long long)local_mem_size);*/
    /*printf("n_local_cols %llu\n", (unsigned long long)n_local_cols);*/
    /*printf("memory %llu\n", (unsigned long long)n_local_cols * n * sizeof(F));*/
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, n_local_cols * col_size, NULL);
    clSetKernelArg(cache->common.kernel, 4, sizeof(ncl), &ncl);
    clSetKernelArg(cache->common.kernel, 5, sizeof(n_local_cols), &n_local_cols);
    clEnqueueNDRangeKernel(cache->common.command_queue, cache->common.kernel, 1, NULL, &global_work_size, &local_work_size, 0, NULL, NULL);
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

void mat_mul_cl_block(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    cl_uint ncl, nblkcl;
    size_t global_work_size[2], local_work_size[2], mat_sizeof, nblk;

    /* Setup variables. */
    global_work_size[0] = n;
    global_work_size[1] = n;
    mat_sizeof = n * n * sizeof(F);
    ncl = n;

    /* Run kernel. */
    common_create_kernel_file(&cache->common, "matmul_block.cl", NULL);
    clGetDeviceInfo(cache->common.device, CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(nblk), &nblk, NULL);
    nblk = sqrt(zmin(nblk, n));
    nblk = zmin(nblk, 3);
    nblkcl = nblk;
    local_work_size[0] = nblk;
    local_work_size[1] = nblk;
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clSetKernelArg(cache->common.kernel, 0, sizeof(cache->buf_a), &cache->buf_a);
    clSetKernelArg(cache->common.kernel, 1, sizeof(cache->buf_b), &cache->buf_b);
    clSetKernelArg(cache->common.kernel, 2, sizeof(cache->buf_c), &cache->buf_c);
    clSetKernelArg(cache->common.kernel, 3, nblk * nblk * sizeof(F), NULL);
    printf("nblk = %llu\n", (unsigned long long)nblk);
    printf("local memory = %llu\n", (unsigned long long)2 * nblk * nblk * sizeof(F));
    clSetKernelArg(cache->common.kernel, 4, nblk * nblk * sizeof(F), NULL);
    clSetKernelArg(cache->common.kernel, 5, sizeof(ncl), &ncl);
    clSetKernelArg(cache->common.kernel, 6, sizeof(nblkcl), &nblkcl);
    clEnqueueNDRangeKernel(
        cache->common.command_queue, cache->common.kernel, 2, NULL,
        global_work_size, local_work_size, 0, NULL, NULL
    );
    clFlush(cache->common.command_queue);
    clFinish(cache->common.command_queue);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

void mat_mul_cl_clblas(const F *A, const F *B, F *C, size_t n, Cache *cache) {
    cl_event event;
    size_t mat_sizeof;

    mat_sizeof = n * n * sizeof(F);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_a, CL_TRUE, 0, mat_sizeof, (F*)A, 0, NULL, NULL);
    clEnqueueWriteBuffer(cache->common.command_queue, cache->buf_b, CL_TRUE, 0, mat_sizeof, (F*)B, 0, NULL, NULL);
    clblasSgemm(
        clblasRowMajor,
        clblasNoTrans,
        clblasNoTrans,

        n,
        n,
        n,
        1.0,

        cache->buf_a,
        0,
        n,

        cache->buf_b,
        0,
        n,

        0.0,
        cache->buf_c,
        0,
        n,

        1,
        &(cache->common.command_queue),
        0,
        NULL,
        &event
    );
    clWaitForEvents(1, &event);
    clEnqueueReadBuffer(cache->common.command_queue, cache->buf_c, CL_TRUE, 0, mat_sizeof, C, 0, NULL, NULL);
}

double bench(MatMul f, const F *A, const F *B, F *C, F *C_ref, size_t n, Cache *cache) {
    double dt, time;
    mat_zero(C, n);
    time = common_get_nanos();
    f(A, B, C, n, cache);
    dt = common_get_nanos() - time;
    if (C_ref != NULL)
        mat_assert_eq(C, C_ref, n);
    printf("%.3e ", dt);
    fflush(stdout);
    return dt;
}

int main(int argc, char **argv) {
    srand(time(NULL));
    Cache cache;
    double max_runtime;
    /* Overly slow ones commented out by default. */
    MatMul mat_mul_funcs[] = {
        /*mat_mul_cpu,*/
        mat_mul_cpu_trans,
        mat_mul_cpu_trans_vec,
        mat_mul_cpu_block,
        mat_mul_cpu_cblas,
        /*mat_mul_cl,*/
        mat_mul_cl_row_priv,
        mat_mul_cl_row_local,
        mat_mul_cl_row_priv_col_local,
        mat_mul_cl_row_priv_cols_local,
        /* TODO broken for larger matrics, some cells contain trash.
         * Likey some memory overflow problem. */
        /*mat_mul_cl_block,*/
        mat_mul_cl_clblas,
    };
    int first, func_done[NELEMS(mat_mul_funcs)] = {0};
    size_t f, i;
    size_t mat_sizeof;

    /* CLI args. */
    if (argc > 1) {
        max_runtime = strtod(argv[1], NULL);
    } else {
        max_runtime = 1.0;
    }

    common_init(&(cache.common), NULL);

    /* Unit test 2x2. */
    {
        const F A[] = {
            1.0, 2.0,
            3.0, 4.0
        };
        const F B[] = {
            5.0, 6.0,
            7.0, 8.0
        };
        enum N { n = 2 };
        F C[n*n];
        const F C_ref[] = {
            19.0, 22.0,
            43.0, 50.0
        };
        cl_buf_init(&cache, n * n * sizeof(F));
        for (f = 0; f < sizeof(mat_mul_funcs)/sizeof(mat_mul_funcs[0]); ++f) {
            mat_zero(C, n);
            mat_mul_funcs[f](A, B, C, n, &cache);
            mat_assert_eq(C, C_ref, n);
        }
        cl_buf_deinit(&cache);
    }

    /* Unit test 4x4. */
    {
        const F A[] = {
             1.0,  2.0,  3.0,  4.0,
             5.0,  6.0,  7.0,  8.0,
             9.0, 10.0, 11.0, 12.0,
            13.0, 14.0, 15.0, 16.0,
        };
        const F B[] = {
            17.0, 18.0, 19.0, 20.0,
            21.0, 22.0, 23.0, 24.0,
            25.0, 26.0, 27.0, 28.0,
            29.0, 30.0, 31.0, 32.0,
        };
        const F C_ref[] = {
             250.0,  260.0,  270.0,  280.0,
             618.0,  644.0,  670.0,  696.0,
             986.0, 1028.0, 1070.0, 1112.0,
            1354.0, 1412.0, 1470.0, 1528.0,
        };
        enum N { n = 4 };
        F C[n*n];
        cl_buf_init(&cache, n * n * sizeof(F));
        for (f = 0; f < NELEMS(mat_mul_funcs); ++f) {
            mat_zero(C, n);
            mat_mul_funcs[f](A, B, C, n, &cache);
            mat_assert_eq(C, C_ref, n);
        }
        cl_buf_deinit(&cache);
    }

    /* Benchmarks. */
    {
        double dt;
        F *A = NULL, *B = NULL, *C = NULL, *C_ref = NULL, *dst = NULL, *ref = NULL;
        int done;
        size_t n = 2;

        puts("#matmul");
        done = 0;
        while(1) {
            printf("%zu ", (size_t)log2(n));
            mat_sizeof = n * n * sizeof(F);

            /* CPU setup. */
            A = aligned_alloc(VECTOR_SIZEOF, mat_sizeof);
            B = aligned_alloc(VECTOR_SIZEOF, mat_sizeof);
            C = aligned_alloc(VECTOR_SIZEOF, mat_sizeof);
            C_ref = aligned_alloc(VECTOR_SIZEOF, mat_sizeof);
            if (NULL == A || NULL == B || NULL == C) {
                printf("error: could not allocate memory for n = %zu", n);
                break;
            }
            mat_rand(A, n);
            mat_rand(B, n);

            cl_buf_init(&cache, mat_sizeof);
            first = 1;
            for (f = 0; f < NELEMS(mat_mul_funcs); ++f) {
                if (func_done[f]) {
                    printf("%*s", 10, "");
                } else {
                    if (first) {
                        dst = C_ref;
                        ref = NULL;
                        first = 0;
                    } else {
                        dst = C;
                        ref = C_ref;
                    }
                    dt = bench(mat_mul_funcs[f], A, B, dst, ref, n, &cache);
                    if (dt > max_runtime)
                        func_done[f] = 1;
                }
            }
            puts("");
            done = 1;
            for (i = 0; i < NELEMS(mat_mul_funcs); ++i) {
                if (!func_done[i]) {
                    done = 0;
                    break;
                }
            }
            if (done)
                break;
            n *= 2;

            /* CPU deinit. */
            free(A);
            free(B);
            free(C);
            free(C_ref);

            cl_buf_deinit(&cache);
        }
        common_deinit(&cache.common);
    }

    return EXIT_SUCCESS;
}

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
#ifndef matops_h
#define matops_h

#include <ceed.h>
#include <petsc.h>
#include "../include/structs.h"

// This function uses libCEED to compute the local action of an operator
PetscErrorCode ApplyLocalCeedOp(Vec X, Vec Y, UserMult user);

// This function uses libCEED to compute the non-linear residual
PetscErrorCode FormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to apply the Jacobian for assembly via a SNES
PetscErrorCode ApplyJacobianCoarse_Ceed(SNES snes, Vec X, Vec Y, void *ctx);

// This function uses libCEED to compute the action of the Jacobian
PetscErrorCode ApplyJacobian_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the prolongation operator
PetscErrorCode Prolong_Ceed(Mat A, Vec X, Vec Y);

// This function uses libCEED to compute the action of the restriction operator
PetscErrorCode Restrict_Ceed(Mat A, Vec X, Vec Y);

// This function returns the computed diagonal of the operator
PetscErrorCode GetDiag_Ceed(Mat A, Vec D);

// This function calculates the strain energy in the final solution
PetscErrorCode ComputeStrainEnergy(DM dm_energy, UserMult user,
                                   CeedOperator op_energy, Vec X,
                                   PetscReal *energy);

// this function checks to see if the computed energy is close enough to reference file energy.
PetscErrorCode RegressionTests_solids(AppCtx app_ctx, PetscReal energy);

#endif // matopts_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
//
//  her2.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 12/23/11.
//  Copyright (c) 2011 University of Colorado Denver. All rights reserved.
//

#ifndef _her2_h
#define _her2_h

/// @file her2.h Performs complex outer product of two vectors.

#include <cctype>
#include "latl.h"

namespace LATL
{
   /// @brief Performs a vector outer product of two complex vectors.
   ///
   /// For a complex Hermitian matrix A, complex vectors x and y, and complex scalar alpha,
   ///
   ///        A := alpha * x * y' + conj(alpha) * y * x' + A
   ///
   /// is computed. The result is Hermitian and is stored as either upper or lower triangular in A.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param uplo Specifies whether A is stored as upper or lower triangular.
   ///
   ///        if uplo = 'U' or 'u' then A is upper triangular
   ///        if uplo = 'L' or 'l' then A is lower triangular
   ///
   /// @param n Specifies the order of the matrix A.  n>=0
   /// @param alpha Complex scalar.
   /// @param x Pointer to complex vector x.
   /// @param incx Increment of the vector x.  x!=0
   /// @param y Pointer to complex vector y.
   /// @param incy Increment of the vector y.  y!=0
   /// @param A Pointer to complex Hermitian n-by-n matrix A.If uplo = 'U' or 'u' then only the upper
   /// triangular part of A is referenced and the lower part is not referenced.  If uplo = 'L' or 'l' then
   /// only the lower triangular part of A is referenced and the upper part is not referenced.
   /// @param ldA Column length of matrix A.  ldA>=n.
   /// @ingroup BLAS

   template <typename real_t>
   int HER2(char uplo, int_t n, complex<real_t> alpha, complex<real_t> *x, int_t incx, complex<real_t> *y, int_t incy, complex<real_t> *A, int_t ldA)
   {
      using std::conj;
      using std::real;
      using std::toupper;
      const complex<real_t> zero(0.0,0.0);
      complex<real_t> tx,ty;
      int_t i,j,kx,ky,jx,ix,jy,iy;

      uplo=toupper(uplo);

      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(incy==0)
         return -7;
      else if(ldA<n)
         return -9;
      else if((n==0)||(alpha==zero))
         return 0;

      if((incx==1)&&(incy==1))
      {
         if(uplo=='U')
         {
            for(j=0;j<n;j++)
            {
               tx=conj(alpha*x[j]);
               ty=alpha*conj(y[j]);
               for(i=0;i<j;i++)
                  A[i]+=x[i]*ty+y[i]*tx;
               A[j]=real(A[j])+real(x[j]*ty+y[j]*tx);
               A+=ldA;
            }
         }
         else
         {
            for(j=0;j<n;j++)
            {
               tx=conj(alpha*x[j]);
               ty=alpha*conj(y[j]);
               A[j]=real(A[j])+real(x[j]*ty+y[j]*tx);
               for(i=j+1;i<n;i++)
                  A[i]+=x[i]*ty+y[i]*tx;
               A+=ldA;
            }
         }
      }
      else
      {
         kx=(incx>0)?0:(1-n)*incx;
         ky=(incy>0)?0:(1-n)*incy;
         if(uplo=='U')
         {
            jx=kx;
            jy=ky;
            for(j=0;j<n;j++)
            {
               tx=conj(alpha*x[jx]);
               ty=alpha*conj(y[jy]);
               ix=kx;
               iy=ky;
               for(i=0;i<j;i++)
               {
                  A[i]+=x[ix]*ty+y[iy]*tx;
                  ix+=incx;
                  iy+=incy;
               }
               A[j]=real(A[j])+real(x[jx]*ty+y[jy]*tx);
               A+=ldA;
               jx+=incx;
               jy+=incy;
            }
         }
         else
         {
            jx=kx;
            jy=ky;
            for(j=0;j<n;j++)
            {
               tx=conj(alpha*x[jx]);
               ty=alpha*conj(y[jy]);
               ix=jx;
               iy=jy;
               A[j]=real(A[j])+real(x[jx]*ty+y[jy]*tx);
               for(i=j+1;i<n;i++)
               {
                  ix+=incx;
                  iy+=incy;
                  A[i]+=x[ix]*ty+y[iy]*tx;
               }
               jx+=incx;
               jy+=incy;
               A+=ldA;
            }
         }
      }
      return 0;
   }
#ifdef __latl_cblas
#include <cblas.h>

   template <> int HER2<float>(char uplo, int_t n, complex<float> alpha, complex<float> *x, int_t incx, complex<float> *y, int_t incy, complex<float> *A, int_t ldA)
   {
      uplo=std::toupper(uplo);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(ldA<n)
         return -7;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_cher2(CblasColMajor,Uplo,n,&alpha,x,incx,y,incy,A,ldA);

      return 0;
   }

   template <> int HER2<double>(char uplo, int_t n, complex<double> alpha, complex<double> *x, int_t incx, complex<double> *y, int_t incy, complex<double> *A, int_t ldA)
   {
      uplo=std::toupper(uplo);
      if((uplo!='U')&&(uplo!='L'))
         return -1;
      else if(n<0)
         return -2;
      else if(incx==0)
         return -5;
      else if(ldA<n)
         return -7;

      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;

      cblas_zher2(CblasColMajor,Uplo,n,&alpha,x,incx,y,incy,A,ldA);

      return 0;
   }
#endif

}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl


text:
/* Copyright (c) 2011-2017, NVIDIA CORPORATION. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *  * Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *  * 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.
 *  * Neither the name of NVIDIA CORPORATION nor the names of its
 *    contributors may be used to endorse or promote products derived
 *    from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``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 COPYRIGHT OWNER OR
 * 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.
 */

#pragma once

namespace amgx
{
template <class Matrix>  typename Matrix::value_type estimate_largest_eigen_value(Matrix &A);
}

#include <norm.h>
#include <multiply.h>
#include <blas.h>

namespace amgx
{

template <class Matrix>
typename Matrix::value_type estimate_largest_eigen_value(Matrix &A)
{
    typedef typename Matrix::TConfig TConfig;
    typedef typename Matrix::value_type ValueTypeA;
    typedef typename TConfig::VecPrec ValueTypeB;
    typedef Vector<TConfig> VVector;
    VVector x(A.get_num_rows()), y(A.get_num_rows());
    fill(x, 1);

    for (int i = 0; i < 20; i++)
    {
        ValueTypeB Lmax = get_norm(A, x, LMAX);
        scal(x, ValueTypeB(1) / Lmax);
        multiply(A, x, y);
        x.swap(y);
    }

    ValueTypeB retval = get_norm(A, x, L2) / get_norm(A, y, L2);
    return retval;
}

} // namespace amgx

max_stars_repo_name: fizmat/AMGX

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
#include <stdbool.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include "numstates.h"

int main(int argc, char *argv[]) {
    int num_bands = 2;
    double t = 1.0;
    double E0 = 6.0;
    double deltaE = 14.0;
    void Efn(double k[3], gsl_vector *energies) {
        double cx = cos(2.0 * M_PI * k[0]);
        double cy = cos(2.0 * M_PI * k[1]);
        double cz = cos(2.0 * M_PI * k[2]);
        double tk = -2.0 * t * (cx + cy + cz);
        int i;
        double E0_band;
        for (i = 0; i < num_bands; i++) {
            E0_band = E0 + ((double)i) * deltaE;
            gsl_vector_set(energies, i, E0_band + tk);
        }
    }
    int na = 8;
    int nb = 8;
    int nc = 8;
    int G_order[3] = {0, 1, 2};
    int G_neg[3] = {1, 1, 1};
    bool use_cache = true;

    EnergyCache *Ecache = init_EnergyCache(na, nb, nc, num_bands, G_order, G_neg, Efn, use_cache);

    double E = 0.0;
    double count = NumStates(E, Ecache);
    double eps = 1e-9;
    double expected_E0 = 0.0;
    if (fabs(count - expected_E0) > eps) {
        printf("Incorrect occupation; got %f, expected %f\n", count, expected_E0);
        return 1;
    }

    E = 6.0;
    count = NumStates(E, Ecache);
    double expected_E6 = 0.5;
    if (fabs(count - expected_E6) > eps) {
        printf("Incorrect occupation; got %f, expected %f\n", count, expected_E6);
        return 1;
    }

    E = 12.0;
    count = NumStates(E, Ecache);
    double expected_E12 = 1.0;
    if (fabs(count - expected_E12) > eps) {
        printf("Incorrect occupation; got %f, expected %f\n", count, expected_E12);
        return 1;
    }

    E = 26.0;
    count = NumStates(E, Ecache);
    double expected_E24 = 2.0;
    if (fabs(count - expected_E24) > eps) {
        printf("Incorrect occupation; got %f, expected %f\n", count, expected_E24);
        return 1;
    }

    printf("Numstates test passed.\n");
    return 0;
}

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
/*
* -----------------------------------------------------------------
*  pmsr_lib.c
*  Pairwise Mixing Stirred Reactor Library
*  Version: 2.0
*  Last Update: Oct 1, 2019
*
*  This is a computational library with routines to implement
*  a Pairwise Mixing Stirred Reactor (PMSR) model.
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2019 by Americo Barbosa da Cunha Junior
* -----------------------------------------------------------------
*/

#include <math.h>
#include <time.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_errno.h>

#include "../include/pmsr_lib.h"

/*
* -----------------------------------------------------------------
*   pmsr_title
*
*   This function prints the program title on the screen.
*
*   last update: Oct 20, 2009
* -----------------------------------------------------------------
*/

void pmsr_title()
{
    time_t curtime = time (NULL);
    struct tm *loctime = localtime (&curtime);

    printf("\n");
    printf("=============================================");
    printf("\n Pairwise Mixing Stirred Reactor --- PMSR\n");
    printf("\n by");
    printf("\n Americo Barbosa da Cunha Junior");
    printf("\n americo.cunhajr@gmail.com\n");
    printf("\n Date:");
    printf("\n %s", asctime(loctime) );
    printf("=============================================");
    printf("\n");

    return;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   pmsr_input
*
*   This function receives PMSR parameters from user
*   and returns GSL_SUCCESS if there is no error.
*
*   Input:
*   seed         - RNG seed
*   Np           - # of particles
*   Ndt          - # of time steps
*   t0           - initial time
*   delta_t      - time step
*   tau_res      - residence time
*   tau_mix      - mixture time
*   tau_pair     - pairwise time
*
*   Output:
*   success or error
*
*   last update: Oct 1, 2019
* -----------------------------------------------------------------
*/

int pmsr_input(unsigned long int *seed,unsigned int *Np,unsigned int *Ndt,
	       double *t0,double *delta_t,double *tau_res,double *tau_mix,
	       double *tau_pair,double *Tmax,double *Tmin,
               double *atol,double *rtol)
{
    printf("\n Input PMSR parameters:\n");

    printf("\n RNG seed:");
    scanf("%lu", seed);
    printf("\n %lu\n", *seed);
    if( *seed <= 0 )
	*seed = time(NULL);

    printf("\n # of particles:");
    scanf("%d", Np);
    printf("\n %d\n", *Np);
    if( *Np <= 0 || *Np % 2 != 0 )
	GSL_ERROR(" Np must be a positive even integer",GSL_EINVAL);

    printf("\n # of time steps:");
    scanf("%d", Ndt);
    printf("\n %d\n", *Ndt);
    if( *Ndt <= 0 )
	GSL_ERROR(" Ndt must be a positive integer",GSL_EINVAL);

    printf("\n initial time:");
    scanf("%lf", t0);
    printf("\n %+.1e\n", *t0);
    if( *t0 < 0.0 )
	GSL_ERROR(" t0 must be a grather than zero",GSL_EINVAL);

    printf("\n time step:");
    scanf("%lf", delta_t);
    printf("\n %+.1e\n", *delta_t);
    if( *delta_t <= 0.0 )
	GSL_ERROR(" delta_t must be a grather than zero",GSL_EINVAL);

    printf("\n residence time:");
    scanf("%lf", tau_res);
    printf("\n %+.1e\n", *tau_res);
    if( *tau_res <= 0.0 )
	GSL_ERROR(" tau_res must be a grather than zero",GSL_EINVAL);

    printf("\n mixture time:");
    scanf("%lf", tau_mix);
    printf("\n %+.1e\n", *tau_mix);
    if( *tau_mix <= 0.0 )
	GSL_ERROR(" tau_mix must be a grather than zero",GSL_EINVAL);

    printf("\n pairwise time:");
    scanf("%lf", tau_pair);
    printf("\n %+.1e\n", *tau_pair);
    if( *tau_pair <= 0.0 )
	GSL_ERROR(" tau_pair must be a grather than zero",GSL_EINVAL);

    printf("\n temperature upper bound:");
    scanf("%lf", Tmax);
    printf("\n %+.1e\n", *Tmax);
    if( *Tmax <= 0.0 )
	GSL_ERROR(" Tmax must be a grather than zero",GSL_EINVAL);

    printf("\n temperature lower bound:");
    scanf("%lf", Tmin);
    printf("\n %+.1e\n", *Tmin);
    if( *Tmin < 0.0 || *Tmin > *Tmax )
	GSL_ERROR(" Tmin must be less than Tmax and a grather than zero",GSL_EINVAL);

    printf("\n absolute tolerance:");
    scanf("%lf", atol);
    printf("\n %+.1e\n", *atol);
    if( *atol <= 0.0 )
	GSL_ERROR(" atol must be a grather than zero",GSL_EINVAL);

    printf("\n relative tolerance:");
    scanf("%lf", rtol);
    printf("\n %+.1e\n", *rtol);
    if( *rtol <= 0.0 )
	GSL_ERROR(" rtol must be a grather than zero",GSL_EINVAL);

    return GSL_SUCCESS;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*  pmsr_alloc
*
*  This function allocates an internal memory block for
*  pmsr workspace struct. If successful, returns a pointer to
*  pmsr_wrk. If a startup error occurs returns NULL.
*
*   last update: Nov 5, 2009
* -----------------------------------------------------------------
*/

pmsr_wrk *pmsr_alloc()
{
    pmsr_wrk *pmsr = NULL;

    /* memory allocation for pmsr struct */
    pmsr = (pmsr_wrk *) malloc(sizeof(pmsr_wrk));
    if ( pmsr == NULL )
        return NULL;

    /* setting pmsr_wrk elements */
    pmsr->pmsr_idx = NULL;
    pmsr->pmsr_ps  = NULL;

    return pmsr;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*  pmsr_init
*
*  This function initiates the system of particles
*  with an initial composition phi.
*
*   Input:
*   Np   - # of particles
*   Neq  - # of equations
*   pmsr - pointer to a struct pmsr
*
*   Output:
*   success or error
*
*  last update: Sep 22, 2009
* -----------------------------------------------------------------
*/

int pmsr_init(int Np,int Neq,pmsr_wrk *pmsr)
{
    int i;

    /* checking if pmsr is NULL */
    if ( pmsr == NULL )
        return GSL_EINVAL;

    /* memory allocation for the system of particles */
    pmsr->pmsr_ps = (gsl_vector **) malloc(Np*sizeof(gsl_vector *));
    if ( pmsr->pmsr_ps == NULL )
	return GSL_ENOMEM;

    for ( i = 0; i < Np; i++ )
        pmsr->pmsr_ps[i] = gsl_vector_calloc(Neq);

    /* memory allocation for the vector of indices */
    pmsr->pmsr_idx = (int *) malloc(Np*sizeof(int));
    if ( pmsr->pmsr_idx == NULL )
    {
	for ( i = 0; i < Np; i++ )
	{
	    gsl_vector_free(pmsr->pmsr_ps[i]);
	    pmsr->pmsr_ps[i] = NULL;
	}

        free(pmsr->pmsr_ps);
	pmsr->pmsr_ps = NULL;

	return GSL_ENOMEM;
    }

    /* initial configuration of particles pairwise */
    for ( i = 0; i < Np; i++ )
        pmsr->pmsr_idx[i] = i;

    return GSL_SUCCESS;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   pmsr_free
*
*   This routine frees the memory allocated by pmsr_alloc.
*
*   Input:
*   pmsr - pointer to a struct pmsr
*
*   Output:
*   void
*
*   last update: May 10, 2009
* -----------------------------------------------------------------
*/

void pmsr_free(int Np,void **pmsr_bl)
{
    int i;
    pmsr_wrk *pmsr = NULL;

    /* checking if pmsr_bl is NULL */
    if ( *pmsr_bl == NULL )
        return;

    /* checking if Np is positive */
    if ( Np <= 0 )
        return;

    pmsr = (pmsr_wrk *) (*pmsr_bl);

    /* releasing pmsr elements */
    if ( pmsr->pmsr_ps != NULL )
    {
	for ( i = 0; i < Np; i++ )
	    if ( pmsr->pmsr_ps[i] != NULL )
	    {
	        gsl_vector_free(pmsr->pmsr_ps[i]);
	        pmsr->pmsr_ps[i] = NULL;
	    }
	/* releasing pmsr struct */
	free(pmsr->pmsr_ps);
	pmsr->pmsr_ps = NULL;
    }

    if ( pmsr->pmsr_idx != NULL )
    {
	free(pmsr->pmsr_idx);
	pmsr->pmsr_idx = NULL;
    }

    free(*pmsr_bl);
    *pmsr_bl = NULL;

    return;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*  pmsr_set_all
*
*  This function sets the system of particles with a given
*  composition phi.
*
*   Input:
*   pmsr - pointer to a struct pmsr
*   phi  - composition
*
*   Output:
*   success or error
*
*  last update: Mar 2, 2010
* -----------------------------------------------------------------
*/

void pmsr_set_all(int Np,gsl_vector *phi,pmsr_wrk *pmsr)
{
    int i;

    /* setting all particles with an initial phi */
    for ( i = 0; i < Np; i++ )
        gsl_vector_memcpy(pmsr->pmsr_ps[i],phi);

    return;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*  pmsr_Nexc
*
*  This function computes the number of particles to be exchanged,
*  Nexc, corresponding to outﬂow/inﬂow particles.
*
*  Nexc  = ceil( 1/2 Np delta_t/taus_res )
*
*   Input:
*   Np      - # of particles
*   delta_t - time step
*   tau_res - residence time
*
*   Output:
*   Nexc    - # of particles to be exchanged
*
*   last update: Feb 12, 2009
* -----------------------------------------------------------------
*/

int pmsr_Nexc(int Np,double delta_t,double tau_res)
{
    return ceil( (Np/2)*(delta_t/tau_res) );
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*  pmsr_Npair
*
*  This function computes the number of particles to be pairwised,
*  Npair, corresponding to pairing particles.
*
*  Npair = ceil( 1/2 Np delta_t/taus_pair )
*
*   Input:
*   Np       - # of particles
*   delta_t  - time step
*   tau_pair - pairwise time
*
*   Output:
*   Npair    - # of particles to be pairwised
*
*   last update: Feb 12, 2009
* -----------------------------------------------------------------
*/

int pmsr_Npair(int Np,double delta_t,double tau_pair)
{
    return ceil( (Np/2)*(delta_t/tau_pair) );
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   pmsr_meanvar
*
*   This function computes the mean and the variance
*   of a particular property in a system of particles
*   using the algorithm of Knuth/Welford.
*
*   Ref.:
*   D. E. Knuth (1998)
*   The Art of Computer Programming,
*   volume 2: Seminumerical Algorithms
*   3rd ed., Addison-Wesley.
*
*   B. P. Welford (1962).
*   "Note on a method for calculating corrected sums
*   of squares and products". Technometrics 4(3):419–420
*
*   Input:
*   Np - # of particles
*   k  - property index
*   ps - system of particle
*
*   Output:
*   mean - propertie mean value
*   var  - propertie variance
*
*   last update: Mar 4, 2010
* -----------------------------------------------------------------
*/

void pmsr_meanvar(int Np,int k,gsl_vector **ps,double *mean,double *var)
{
    int i;
    double M     = 0.0;
    double S     = 0.0;
    double delta = 0.0;

    for ( i = 0; i < Np; i++ )
    {
        delta = ps[i]->data[k] - M;
	M    += delta/(double) (i+1);
        S    += delta*(ps[i]->data[k] - M);
    }

    /* mean */
    *mean  = M;

    /* variance */
    *var   = S/(double) (i-1);

    return;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   pmsr_mixture
*
*   This function computes the mixture step in PMSR model.
*
*   Input:
*   Np      - # of particles
*   Neq     - # of equations
*   delta_t - time step
*   tau_mix - mixture time
*   idx     - particles index vector
*
*   Output:
*   ps      - system of particles
*
*   last update: Mar 9, 2009
* -----------------------------------------------------------------
*/

void pmsr_mixture(int Np,int Neq,double delta_t,double tau_mix,
					    int *idx,gsl_vector **ps)
{
    int k, i;
    double sum, dif;

    /* particle index */
    k = 0;

    /* mixture of the Np particles */
    while( k < Np/2 )
    {
        /* mixture of two particles */
        for ( i = 0; i < Neq; i++ )
        {
            /* sum = 1/2(phi_p0 + phi_q0) */
            sum = 0.5*( ps[idx[2*k]]->data[i] + ps[idx[2*k+1]]->data[i] );

            /* dif = 1/2(phi_p0 - phi_q0) */
            dif = 0.5*( ps[idx[2*k]]->data[i] - ps[idx[2*k+1]]->data[i] );

            /* phi_p =  dif.exp(-2/tau delta_t) + sum */
            ps[idx[2*k]]->data[i] =
                        dif*exp(-(2.0/tau_mix)*delta_t) + sum;

            /* phi_q = -dif.exp(-2/tau delta_t) + sum */
            ps[idx[2*k+1]]->data[i] =
                        -dif*exp(-(2.0/tau_mix)*delta_t) + sum;
        }
        k++;
    }

    return;
}
/*----------------------------------------------------------------*/

/*
* -----------------------------------------------------------------
*   pmsr_iops
*
*   This function executes the inflow/outflow,
*   pairwise and shuflle steps.
*
*   Input:
*   Np      - # of particles
*   Nexc    - # of particles to be exchanged
*   Npair   - # of particles to be pairwised
*   rand    - rng workspace
*   phi     - inflow composition
*
*   Output:
*   p       - system of particles
*   pairs   - pairs of particles index
*
*   last update: Mar 1, 2010
* -----------------------------------------------------------------
*/

void pmsr_iops(int Np,int Nexc,int Npair,gsl_rng *rand,
				    gsl_vector *phi,pmsr_wrk *pmsr)
{
    int i;
    int idx1[Nexc+Npair];
    int idx2[2*(Nexc+Npair)];
    int *aux = NULL;

    /* memory allocation */
    aux = (int *) malloc((Np/2)*sizeof(int));

    /* setting aux vector with integers from 0 to Np/2-1 */
    for ( i = 0; i < Np/2; i++ )
        aux[i] = i;

    /* selecting Nexc + Npair pairs at random */
    gsl_ran_choose(rand,idx1,Nexc+Npair,aux,Np/2,sizeof(int));
    for ( i = 0; i < Nexc + Npair; i++ )
    {
        idx2[2*i  ] = pmsr->pmsr_idx[2*idx1[i]  ];
        idx2[2*i+1] = pmsr->pmsr_idx[2*idx1[i]+1];
    }

    /* changing composition of Nexc pairs (inflow / outflow) */
    for ( i = 0; i < Nexc; i++ )
    {
        gsl_vector_memcpy(pmsr->pmsr_ps[idx2[2*i]  ],phi);
        gsl_vector_memcpy(pmsr->pmsr_ps[idx2[2*i+1]],phi);
    }

    /* shuffling the selected pairs */
    gsl_ran_shuffle(rand,idx2,2*(Nexc+Npair),sizeof(int));

    /* setting the new partners in the system */
    for ( i = 0; i < Nexc + Npair; i++ )
    {
        pmsr->pmsr_idx[2*idx1[i]  ] = idx2[2*i  ];
        pmsr->pmsr_idx[2*idx1[i]+1] = idx2[2*i+1];
    }

    /* releasing allocated memory */
    free(aux);
    aux = NULL;

    return;
}
/*----------------------------------------------------------------*/

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#ifndef _createK2_h
#define _createK2_h

#include <petsc.h>
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>

PetscErrorCode bsscr_buildK2(KSP ksp);
PetscErrorCode bsscr_DGMiGtD( Mat *_K2, Mat K, Mat G, Mat M);
PetscErrorCode bsscr_GMiGt( Mat *_K2, Mat K, Mat G, Mat M);
PetscErrorCode bsscr_GGt( Mat *_K2, Mat K, Mat G);
#endif