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

#include <config.h>
#include <time.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_errno.h>

void benchmark (const gsl_rng_type * T);

#define N  1000000
int isum;
double dsum;

int
main (void)
{
  benchmark(gsl_rng_ranlux);
  benchmark(gsl_rng_ranlux389);
  benchmark(gsl_rng_ranlxs0);
  benchmark(gsl_rng_ranlxs1);
  benchmark(gsl_rng_ranlxs2);
  benchmark(gsl_rng_ranlxd1);
  benchmark(gsl_rng_ranlxd2);

  benchmark(gsl_rng_slatec);
  benchmark(gsl_rng_gfsr4);
  benchmark(gsl_rng_cmrg);
  benchmark(gsl_rng_minstd);
  benchmark(gsl_rng_mrg);
  benchmark(gsl_rng_mt19937);
  benchmark(gsl_rng_r250);
  benchmark(gsl_rng_ran0);
  benchmark(gsl_rng_ran1);
  benchmark(gsl_rng_ran2);
  benchmark(gsl_rng_ran3);
  benchmark(gsl_rng_rand48);
  benchmark(gsl_rng_rand);
  benchmark(gsl_rng_random8_bsd);
  benchmark(gsl_rng_random8_glibc2);
  benchmark(gsl_rng_random8_libc5);
  benchmark(gsl_rng_random128_bsd);
  benchmark(gsl_rng_random128_glibc2);
  benchmark(gsl_rng_random128_libc5);
  benchmark(gsl_rng_random256_bsd);
  benchmark(gsl_rng_random256_glibc2);
  benchmark(gsl_rng_random256_libc5);
  benchmark(gsl_rng_random32_bsd);
  benchmark(gsl_rng_random32_glibc2);
  benchmark(gsl_rng_random32_libc5);
  benchmark(gsl_rng_random64_bsd);
  benchmark(gsl_rng_random64_glibc2);
  benchmark(gsl_rng_random64_libc5);
  benchmark(gsl_rng_random_bsd);
  benchmark(gsl_rng_random_glibc2);
  benchmark(gsl_rng_random_libc5);
  benchmark(gsl_rng_randu);
  benchmark(gsl_rng_ranf);

  benchmark(gsl_rng_ranmar);
  benchmark(gsl_rng_taus);
  benchmark(gsl_rng_transputer);
  benchmark(gsl_rng_tt800);
  benchmark(gsl_rng_uni32);
  benchmark(gsl_rng_uni);
  benchmark(gsl_rng_vax);
  benchmark(gsl_rng_zuf);
  return 0;
}

void
benchmark (const gsl_rng_type * T)
{
  int start, end;
  int i = 0, d = 0 ;
  double t1, t2;

  gsl_rng *r = gsl_rng_alloc (T);

  start = clock ();
  do
    {
      int j;
      for (j = 0; j < N; j++)
	isum += gsl_rng_get (r);

      i += N;
      end = clock ();
    }
  while (end < start + CLOCKS_PER_SEC/10);

  t1 = (end - start) / (double) CLOCKS_PER_SEC;

  start = clock ();
  do
    {
      int j;
      for (j = 0; j < N; j++)
	dsum += gsl_rng_uniform (r);

      d += N;
      end = clock ();
    }
  while (end < start + CLOCKS_PER_SEC/10);

  t2 = (end - start) / (double) CLOCKS_PER_SEC;

  printf ("%6.0f k ints/sec, %6.0f k doubles/sec, %s\n",
	  i / t1 / 1000.0, d / t2 / 1000.0, gsl_rng_name (r));

  gsl_rng_free (r);
}


text:
/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */

// $Revision: 1.2 $
// $Date: 2007-02-14 20:12:32 $
// $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/petsc/ShadowPetscSOE.h,v $

// File: ~/system_of_eqn/linearSOE/petsc/ShadowPetscSOE.h
//
// Written: fmk & om
// Created: 7/98
// Revision: A
//
// Description: This file contains the class definition for ShadowPetscSOE
// ShadowPetscSOE is a subclass of LinearSOE.

// What: "@(#) ShadowPetscSOE.h, revA"

#ifndef ShadowPetscSOE_h
#define ShadowPetscSOE_h

#include <LinearSOE.h>
#include <Vector.h>
#include <PetscSOE.h>

extern "C" {
#include <petsc.h>
}

class PetscSolver;

class ShadowPetscSOE : public LinearSOE
{
  public:
    ShadowPetscSOE(PetscSolver &theSolver, int blockSize);

    ~ShadowPetscSOE();

    int solve(void);

    int getNumEqn(void) const;
    int setSize(Graph &theGraph);

    int addA(const Matrix &, const ID &, double fact = 1.0);
    int addB(const Vector &, const ID &, double fact = 1.0);
    int setB(const Vector &, double fact = 1.0);

    void zeroA(void);
    void zeroB(void);

    const Vector &getX(void);
    const Vector &getB(void);
    double normRHS(void);

    void setX(int loc, double value);

    int setSolver(PetscSolver &newSolver);

    int sendSelf(int commitTag, Channel &theChannel);
    int recvSelf(int commitTag, Channel &theChannel,
		 FEM_ObjectBroker &theBroker);

  protected:

  private:
    MPI_Comm theComm; // a comm for communicating to the ActorPetscSOE's
                      // without using PETSC_COMM_WORLD
    PetscSOE theSOE;  // the local portion of the SOE
    PetscSolver *theSolver; // created by the user
    int myRank;
    int numProcessors;
    int sendData[3];
    void *sendBuffer;
    int blockSize;
};

#endif

ext: h

max_stars_repo_name: kuanshi/ductile-fracture


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

#include "outcome/outcome.hpp"

#include <gsl/span>

namespace katla {

enum class MqttQos {AtMostOnce = 0, AtLeastOnce = 1, ExactlyOnce = 2};

struct MqttMessage {
    int messageId {};
    std::string topic;
    gsl::span<std::byte> payload;
    MqttQos qos;
    bool retain;
};

} // namespace katla

#endif

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla


text:
#pragma once
#include <gsl/multi_span>
#include <droidCrypto/Defines.h>
#include <droidCrypto/MatrixView.h>
#include <cstring>

namespace droidCrypto
{
    enum class AllocType
    {
        Uninitialized,
        Zeroed
    };

    template<typename T>
    class Matrix : public MatrixView<T>
    {
        uint64_t mCapacity = 0;
    public:
        Matrix()
        {}

        Matrix(uint64_t rows, uint64_t columns, AllocType t = AllocType::Zeroed)
        {
            resize(rows, columns, t);
        }

        Matrix(const MatrixView<T>& copy)
            : MatrixView<T>(new T[copy.size()], copy.bounds()[0], copy.stride())
            , mCapacity(copy.size())
        {
            memcpy(MatrixView<T>::mView.data(), copy.data(), copy.mView.size_bytes());
        }

        Matrix(Matrix<T>&& copy)
            : MatrixView<T>(copy.data(), copy.bounds()[0], copy.stride())
            , mCapacity(copy.mCapacity)
        {
            copy.mView = span<T>();
            copy.mStride = 0;
            copy.mCapacity = 0;
        }

        ~Matrix()
        {
            delete[] MatrixView<T>::mView.data();
        }

        const Matrix<T>& operator=(const Matrix<T>& copy)
        {
            resize(copy.rows(), copy.stride());
            memcpy(MatrixView<T>::mView.data(), copy.data(), copy.mView.size_bytes());
            return copy;
        }

        void resize(uint64_t rows, uint64_t columns, AllocType type = AllocType::Zeroed)
        {
            if (rows * columns > mCapacity)
            {
                mCapacity = rows * columns;
                auto old = MatrixView<T>::mView;

                if (type == AllocType::Zeroed)
                    MatrixView<T>::mView = span<T>(new T[mCapacity](), mCapacity);
                else
                    MatrixView<T>::mView = span<T>(new T[mCapacity], mCapacity);

                auto min = std::min<uint64_t>(old.size(), mCapacity) * sizeof(T);
                memcpy(MatrixView<T>::mView.data(), old.data(), min);

                delete[] old.data();

            }
            else
            {
                auto newSize = rows * columns;
                if (newSize > MatrixView<T>::size() && type == AllocType::Zeroed)
                {
                    memset(MatrixView<T>::data() + MatrixView<T>::size(), 0, newSize - MatrixView<T>::size());
                }

                MatrixView<T>::mView = span<T>(MatrixView<T>::data(), newSize);
            }

            MatrixView<T>::mStride = columns;
        }

        // return the internal memory, stop managing its lifetime, and set the current container to null.
        T* release()
        {
            auto ret = MatrixView<T>::mView.data();
            MatrixView<T>::mView = {};
            mCapacity = 0;
            return ret;
        }
    };

}


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <time.h>

/*
  gsl_matrix_printf prints a matrix as a column vector.  This function
  prints a matrix in block form.
*/
void pretty_print(const gsl_matrix * M)
{
  // Get the dimension of the matrix.
  int rows = M->size1;
  int cols = M->size2;
  // Now print out the data in a square format.
  for(int i = 0; i < rows; i++){
    for(int j = 0; j < cols; j++){
      printf("%f ", gsl_matrix_get(M, i, j));
    }
    printf("\n");
  }
  printf("\n");
}

int run_svd(const gsl_matrix * a) {
  // Need to transpose the input
  gsl_matrix *a_transpose = gsl_matrix_alloc(a->size2, a->size1);
  gsl_matrix_transpose_memcpy(a_transpose, a);
  printf("a_matrix'\n");
  pretty_print(a_transpose);

  int m = a->size1;
  gsl_matrix * V = gsl_matrix_alloc(m, m);
  gsl_vector * S = gsl_vector_alloc(m);
  gsl_vector * work = gsl_vector_alloc(m);

  gsl_linalg_SV_decomp(a_transpose, V, S, work);
  printf("U\n");
  pretty_print(V); printf("\n");
  printf("V\n");
  pretty_print(a_transpose); printf("\n");
  printf("S\n");
  gsl_vector_fprintf(stdout, S, "%g");

  gsl_matrix_free(V);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return 0;
}

ext: c

max_stars_repo_name: paul-reiners/getting-smaller

max_issues_repo_name: paul-reiners/getting-smaller

max_forks_repo_name: paul-reiners/getting-smaller


text:
#ifndef DOG_H
#define DOG_H

/* from local library */
#include "colorhistogram.h"

/* from OpenCV library */
#include <opencv2/core.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/imgproc.hpp>

/* from gsl library */
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* from gflags library */
#include <gflags/gflags.h>

// use console flag
DECLARE_double(std);

class Dog {
public:
    Dog();

    Dog(const int fw, const int fh, const cv::Rect &rect, cv::Mat *hist);

    void run(cv::Mat &frame, gsl_rng *rng);
    //bool operator<(const Dog &p) { return weight < p.weight; }

    /* coefficient for normalizing */
    static const int LAMBDA;

    /* autoregressive dynamics parameters for transition model */
    static const float A1;
    static const float A2;
    static const float B0;

    /* standard deviations for gaussian sampling in transition model */
    static const double TRANS_X_STD;
    static const double TRANS_Y_STD;
    static const double TRANS_S_STD;

    int fw;
    int fh;
    float x0;
    float y0;
    float xp;
    float yp;
    float x;
    float y;
    float sp;
    float s;
    float width;
    float height;
    cv::Mat *hist;
    float weight;

private:
    static ColorHistogram *pHist;

    float _likelihood(const cv::Mat &frame);

    float _histo_dist_sq(cv::Mat *h1, cv::Mat *h2);
};

#endif

max_stars_repo_name: Redoblue/find-the-bombs

max_issues_repo_name: Redoblue/find-the-bombs

max_forks_repo_name: Redoblue/find-the-bombs


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

#ifndef __GSL_FFT_REAL_FLOAT_H__
#define __GSL_FFT_REAL_FLOAT_H__

#include <stddef.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_fft.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

int gsl_fft_real_float_radix2_transform (float data[], const size_t stride, const size_t n) ;

typedef struct
  {
    size_t n;
    size_t nf;
    size_t factor[64];
    gsl_complex_float *twiddle[64];
    gsl_complex_float *trig;
  }
gsl_fft_real_wavetable_float;

typedef struct
  {
    size_t n;
    float *scratch;
  }
gsl_fft_real_workspace_float;

gsl_fft_real_wavetable_float * gsl_fft_real_wavetable_float_alloc (size_t n);

void  gsl_fft_real_wavetable_float_free (gsl_fft_real_wavetable_float * wavetable);

gsl_fft_real_workspace_float * gsl_fft_real_workspace_float_alloc (size_t n);

void  gsl_fft_real_workspace_float_free (gsl_fft_real_workspace_float * workspace);

int gsl_fft_real_float_transform (float data[], const size_t stride, const size_t n,
                                  const gsl_fft_real_wavetable_float * wavetable,
                                  gsl_fft_real_workspace_float * work);

int gsl_fft_real_float_unpack (const float real_float_coefficient[],
			       float complex_coefficient[],
			       const size_t stride, const size_t n);

__END_DECLS

#endif /* __GSL_FFT_REAL_FLOAT_H__ */


text:
#include <stdio.h>
#include <gsl/gsl_rng.h>

int
main (void)
{
  const gsl_rng_type * T;
  gsl_rng * r;

  int i, n = 10;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (i = 0; i < n; i++)
    {
      double u = gsl_rng_uniform (r);
      printf ("%.5f\n", u);
    }

  gsl_rng_free (r);

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <gsl/gsl_sf_bessel.h>

int
main (void)
{
  double x = 5.0;
  double expected = -0.17759677131433830434739701;

  double y = gsl_sf_bessel_J0 (x);

  printf ("J0(5.0) = %.18f\n", y);
  printf ("exact   = %.18f\n", expected);
  return 0;
}


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``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 Javier Burguete Tolosa 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.
*/

/**
 * \file rk_2_2.c
 * \brief Source file to optimize Runge-Kutta 3 steps 2nd order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_2_2.h"

#define DEBUG_RK_2_2 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 2 steps 2nd order Runge-Kutta
 * method.
 */
int
rk_tb_2_2 (Optimize * optimize) ///< Optimize struct.
{
  long double *tb, *r;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_tb_2_2: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t2 (tb) = 1.L;
  t1 (tb) = r[0];
  b21 (tb) = 0.5L / t1 (tb);
  rk_b_2 (tb);
#if DEBUG_RK_2_2
  rk_print_tb (optimize, "rk_tb_2_2", stderr);
  fprintf (stderr, "rk_tb_2_2: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 2 steps 2nd order, 3rd order in
 * equations depending only in time, Runge-Kutta method.
 */
int
rk_tb_2_2t (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_tb_2_2t: start\n");
#endif
  tb = optimize->coefficient;
  t2 (tb) = 1.L;
  t1 (tb) = 2.L / 3.L;
  b21 (tb) = 0.5L / t1 (tb);
  rk_b_2 (tb);
#if DEBUG_RK_2_2
  rk_print_tb (optimize, "rk_tb_2_2t", stderr);
  fprintf (stderr, "rk_tb_2_2t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 2 steps 1st-2nd order Runge-Kutta
 * pair.
 */
int
rk_tb_2_2p (Optimize * optimize)        ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_tb_2_2p: start\n");
#endif
  if (!rk_tb_2_2 (optimize))
    return 0;
  tb = optimize->coefficient;
  e20 (tb) = 1.L;
#if DEBUG_RK_2_2
  rk_print_tb (optimize, "rk_tb_2_2p", stderr);
  fprintf (stderr, "rk_tb_2_2p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 2 steps 1st-2nd order, 1st-3rd order
 * in equations depending only on time, Runge-Kutta pair.
 */
int
rk_tb_2_2tp (Optimize * optimize)       ///< Optimize struct.
{
  long double *tb;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_tb_2_2tp: start\n");
#endif
  if (!rk_tb_2_2t (optimize))
    return 0;
  tb = optimize->coefficient;
  e20 (tb) = 1.L;
#if DEBUG_RK_2_2
  rk_print_tb (optimize, "rk_tb_2_2tp", stderr);
  fprintf (stderr, "rk_tb_2_2tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 2 steps 2nd order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_2_2 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_objective_tb_2_2: start\n");
#endif
  tb = rk->tb->coefficient;
  if (b20 (tb) < 0.L)
    {
      o = 40.L - b20 (tb);
      goto end;
    }
  o = 30.L + fmaxl (1.L, t1 (tb));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_objective_tb_2_2: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_2_2: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 2 steps 2nd order, 3rd
 * order in equations depending only in time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_2_2t (RK * rk)  ///< RK struct.
{
  long double o;
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_objective_tb_2_2t: start\n");
#endif
  o = 31.L;
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
#if DEBUG_RK_2_2
  fprintf (stderr, "rk_objective_tb_2_2t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_2_2t: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
#ifndef CWANNIER_SPINORBIT_H
#define CWANNIER_SPINORBIT_H

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include "HTightBinding.h"

typedef struct {
    gsl_complex val;
    int row;
    int col;
} socElem;

HTightBinding* HamiltonianWithSOC(double soc_strength, double theta, double phi, HTightBinding *Hrs_up, HTightBinding *Hrs_dn);

void addSOCElems(gsl_matrix_complex *Hr_onsite, double soc_strength, double theta, double phi);

gsl_matrix_complex* onSiteSOC_SpinZ();

void setWithHermitianConjugate(gsl_matrix_complex *M, int row, int col, gsl_complex val);

gsl_matrix_complex* onSiteSOC(double theta, double phi);

gsl_matrix_complex *RotationMatrix(double theta, double phi);
#endif // CWANNIER_SPINORBIT_H

max_stars_repo_name: tflovorn/cwannier

max_issues_repo_name: tflovorn/cwannier

max_forks_repo_name: tflovorn/cwannier


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

#include <stdlib.h>
#include <fftw.h>

#include "mex.h"

/**************************************************************************/

/* MEX programs need to use special memory allocation routines,
   so we use the hooks provided by FFTW to ensure that MEX
   allocation is used: */

void *fftw_mex_malloc_hook(size_t n)
{
     void *buf;

     buf = mxMalloc(n);

     /* Call this routine so that we can retain allocations and
	data between calls to the FFTW MEX: */
     mexMakeMemoryPersistent(buf);

     return buf;
}

void fftw_mex_free_hook(void *buf)
{
     mxFree(buf);
}

void install_fftw_hooks(void)
{
     fftw_malloc_hook = fftw_mex_malloc_hook;
     fftw_free_hook = fftw_mex_free_hook;
}

/**************************************************************************/

/* We retain various information between calls to the FFTW MEX in
   order to maximize performance.  (Reusing plans, data, and
   allocated blocks where possible.)  This information is referenced
   by the following variables, which are initialized in the function
   initialize_fftw_mex_data. */

#define MAX_RANK 10

int first_call = 1; /* 1 if this is the first call to the FFTW MEX,
		       and nothing has been initialized yet.  0 otherwise. */

/* Keep track of the array dimensions that stored data (below) is for.
   When these dimensions changed, we have to recompute the plans,
   work arrays, etc. */
int  cur_rank = 0,         /* rank of the array */
     cur_dims[MAX_RANK],   /* dimensions */
     cur_N;                /* product of the dimensions */

/* Work arrays.  MATLAB stores complex numbers as separate real/imag.
   arrays, so we have to translate into our format before the FFT.
   In case allocation is slow, we retain these work arrays between
   calls so that they can be reused. */
fftw_complex *input_work = 0, *output_work = 0;

/* The number of floating point operations required for the FFT.
   (Starting with FFTW 1.3, an exact count is computed by the planner.)
   This is used to update MATLAB's flops count. */
int fftw_mex_flops = 0, ifftw_mex_flops = 0;

/* Plans.  These are computed once and then reused as long as the
   dimensions of the array don't changed.  At any point in time,
   at most two plans are cached: a forward and backwards plan,
   either for one- or multi-dimensional transforms. */
fftw_plan p = 0, ip = 0;
fftwnd_plan pnd = 0, ipnd = 0;

/**************************************************************************/

int compute_fftw_mex_flops(fftw_direction dir)
{
#ifdef FFTW_HAS_COUNT_PLAN_OPS /* this feature will be in FFTW 1.3 */
     fftw_op_count ops;

     if (dir == FFTW_FORWARD) {
	  if (cur_rank == 1)
	       fftw_count_plan_ops(p,&ops);
	  else
	       fftwnd_count_plan_ops(pnd,&ops);
     }
     else {
	  if (cur_rank == 1)
	       fftw_count_plan_ops(ip,&ops);
	  else
	       fftwnd_count_plan_ops(ipnd,&ops);
     }

     return (ops.fp_additions + ops.fp_multiplications);
#else
     return 0;
#endif
}

/**************************************************************************/

/* The following functions destroy and/or initialize the data that
   FFTW-MEX caches between calls. */

void destroy_fftw_mex_data(void) {
     if (output_work != input_work)
	  fftw_free(output_work);
     if (input_work)
	  fftw_free(input_work);
     if (p)
	  fftw_destroy_plan(p);
     if (pnd)
	  fftwnd_destroy_plan(pnd);
     if (ip)
	  fftw_destroy_plan(ip);
     if (ipnd)
	  fftwnd_destroy_plan(ipnd);

     cur_rank = 0;
     input_work = output_work = 0;
     ip = p = 0;
     ipnd = pnd = 0;
}

/* This function is called when MATLAB exits or the MEX file is
   cleared, in which case we want to dispose of all data and
   free any allocated blocks. */

void fftw_mex_exit_function(void)
{
     if (!first_call) {
	  destroy_fftw_mex_data();
	  fftw_forget_wisdom();
	  first_call = 1;
     }
}

#define MAGIC(x) #x
#define STRINGIZE(x) MAGIC(x)

/* Initialize the cached data each time the MEX file is called.  First,
   we check if we have previously computed plans and data for these
   array dimensions.  Only if the dimensions have changed since the
   last call must we recompute the plans, etc. */

void initialize_fftw_mex_data(int rank, const int *dims, fftw_direction dir)
{
     int new_plan = 0;

     if (first_call) {
	  /* The following things need only be done once: */
	  install_fftw_hooks();
	  mexAtExit(fftw_mex_exit_function);
	  first_call = 0;
     }

     if (rank == 1) {
	  if (cur_rank != 1 || cur_dims[0] != dims[0]) {
	       destroy_fftw_mex_data();

	       cur_rank = 1;
	       cur_dims[0] = cur_N = dims[0];

	       input_work = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * cur_N);
	       output_work = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * cur_N);

	       new_plan = 1;
	  }
	  else if (dir == FFTW_FORWARD && !p ||
		   dir == FFTW_BACKWARD && !ip)
	       new_plan = 1;

	  if (new_plan) {
	       if (dir == FFTW_FORWARD) {
		    p = fftw_create_plan(cur_N,dir,
					 FFTW_MEASURE | FFTW_USE_WISDOM);

		    fftw_mex_flops = compute_fftw_mex_flops(dir);
	       }
	       else {
		    ip = fftw_create_plan(cur_N,dir,
					  FFTW_MEASURE | FFTW_USE_WISDOM);

		    ifftw_mex_flops = compute_fftw_mex_flops(dir);
	       }
	  }
     }
     else {
	  int same_dims = 1, dim;

	  if (cur_rank == rank)
	       for (dim = 0; dim < rank && same_dims; ++dim)
		    same_dims = (cur_dims[dim] == dims[rank-1-dim]);
	  else
	       same_dims = 0;

	  if (!same_dims) {
	       if (rank > MAX_RANK)
		    mexErrMsgTxt("Sorry, dimensionality > " STRINGIZE(MAX_RANK)
				 " is not supported.");

	       destroy_fftw_mex_data();

	       cur_rank = rank;

	       cur_N = 1;
	       for (dim = 0; dim < rank; ++dim)
		    cur_N *= (cur_dims[dim] = dims[rank-1-dim]);

	       input_work = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * cur_N);
	       output_work = input_work;

	       new_plan = 1;
	  }
          else if (dir == FFTW_FORWARD && !pnd ||
                   dir == FFTW_BACKWARD && !ipnd)
               new_plan = 1;

	  if (new_plan) {
	       if (dir == FFTW_FORWARD) {
		    pnd = fftwnd_create_plan(rank,cur_dims,dir,
					     FFTW_IN_PLACE |
					     FFTW_MEASURE | FFTW_USE_WISDOM);

		    fftw_mex_flops = compute_fftw_mex_flops(dir);
	       }
	       else {
		    ipnd = fftwnd_create_plan(rank,cur_dims,dir,
					      FFTW_IN_PLACE |
					      FFTW_MEASURE | FFTW_USE_WISDOM);

		    ifftw_mex_flops = compute_fftw_mex_flops(dir);
	       }
	  }

     }
}

/**************************************************************************/

/* MATLAB stores complex numbers as separate arrays for real and
   imaginary parts.  The following functions take the data in
   this format and pack it into a fftw_complex work array, or
   unpack it, respectively. The globals input_work and output_work
   are used as the arrays to pack to/unpack from.*/

void pack_input_work(double *input_re, double *input_im)
{
     int i;

     if (input_im)
	  for (i = 0; i < cur_N; ++i) {
	       c_re(input_work[i]) = input_re[i];
	       c_im(input_work[i]) = input_im[i];
	  }
     else
	  for (i = 0; i < cur_N; ++i) {
	       c_re(input_work[i]) = input_re[i];
	       c_im(input_work[i]) = 0.0;
	  }
}

void unpack_output_work(double *output_re, double *output_im)
{
     int i;

     for (i = 0; i < cur_N; ++i) {
	  output_re[i] = c_re(output_work[i]);
	  output_im[i] = c_im(output_work[i]);
     }
}

/**************************************************************************/

/* The following function is called by MATLAB when the FFTW
   MEX is invoked from within the program.

   The rhs parameters are the list of arrays on the right-hand-side
   (rhs) of the MATLAB command--the arguments to FFTW.  The lhs
   parameters are the list of arrays on the left-hand-side (lhs) of
   the MATLAB command--these are what the output(s) of FFTW are
   assigned to.

   The syntax for the FFTW call in MATLAB is fftw(array,sign),
   as described in fftw.m  */

void mexFunction(int nlhs, mxArray *plhs[],
		 int nrhs, const mxArray *prhs[])
{
     int rank;
     const int *dims;
     int m, n;  /* Array is m x n, C-ordered */
     fftw_direction dir;

     if (nrhs != 2)
	  mexErrMsgTxt("Two input arguments are expected.");

     if (!mxIsDouble(prhs[0]))
	  mexErrMsgTxt("First input must be a double precision matrix.");
     if (mxIsSparse(prhs[0]))
	  mexErrMsgTxt("Sorry, sparse matrices are not currently supported.");

     if (mxGetM(prhs[1]) * mxGetN(prhs[1]) != 1)
	  mexErrMsgTxt("Second input must be a scalar (+/- 1).");

     if (mxGetScalar(prhs[1]) > 0.0)
	  dir = FFTW_BACKWARD;
     else
	  dir = FFTW_FORWARD;

     if ((rank = mxGetNumberOfDimensions(prhs[0])) == 2) {
	  int dims2[2];
	  m = mxGetM(prhs[0]);
	  n = mxGetN(prhs[0]);
	  if (m == 1 || n == 1) {
	       dims2[0] = m * n;
	       initialize_fftw_mex_data(1,dims2,dir);
	  }
	  else {
	       dims2[0] = m;
	       dims2[1] = n;
	       initialize_fftw_mex_data(2,dims2,dir);
	  }
     }
     else
	  initialize_fftw_mex_data(rank,dims = mxGetDimensions(prhs[0]),dir);

     pack_input_work(mxGetPr(prhs[0]),mxGetPi(prhs[0]));

     if (dir == FFTW_FORWARD) {
	  if (cur_rank == 1)
	       fftw(p,1, input_work,1,0, output_work,1,0);
	  else
	       fftwnd(pnd,1, input_work,1,0, 0,0,0);

	  mexAddFlops(fftw_mex_flops);
     }
     else {
	  if (cur_rank == 1)
	       fftw(ip,1, input_work,1,0, output_work,1,0);
	  else
	       fftwnd(ipnd,1, input_work,1,0, 0,0,0);

	  mexAddFlops(ifftw_mex_flops);
     }

     /* Create a matrix for the return argument. */
     if (cur_rank <= 2)
	  plhs[0] = mxCreateDoubleMatrix(m, n, mxCOMPLEX);
     else
	  plhs[0] = mxCreateNumericArray(rank,dims,
					 mxDOUBLE_CLASS,mxCOMPLEX);

     unpack_output_work(mxGetPr(plhs[0]),mxGetPi(plhs[0]));
}

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#ifndef __CAD_UTILC_HH__
#define __CAD_UTILC_HH__

/********************************* TRICK HEADER *******************************
PURPOSE:
      (CAD Utility functions C Version)
LIBRARY DEPENDENCY:
      ((../src/cad_utility_c.c)
        (../src/global_constants.c))
*******************************************************************************/

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <math.h>
#include <assert.h>
#include <stdlib.h>

#ifdef __cplusplus
extern "C" {
#endif
/** Global constants **/
extern const double __FLATTENING;
extern const double __SMAJOR_AXIS;
extern const double __SMALL;
extern const double __RAD;
extern const double __GM;

/** Functions **/
int cad_in_geo84(const double lon, const double lat
                    , const double alt, const gsl_matrix *TEI, gsl_vector *SBII);
int cad_tdi84(const double lon, const double lat, const double alt, const gsl_matrix *TEI, gsl_matrix *TDI);
int cad_geo84_in(const gsl_vector *SBII, const gsl_matrix *TEI, double *lon, double *lat, double *alt);
int cad_kepler(const gsl_vector *SBII, const gsl_vector *VBII, const double tgo
                , gsl_vector *SPII, gsl_vector *VPII);
#ifdef __cplusplus
}
#endif
#endif  //  __CAD_UTILC_HH__


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2, 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.

 $Id$
*/

#ifndef _LIB_OCT_H
#define _LIB_OCT_H

#include <gsl/gsl_complex.h>

int         parse_init   (char *file_out, int *mpiv_node);
int         parse_input  (char *file_in);
void        parse_end    (void);

int         parse_isdef  (char *name);
int         parse_int    (char *name, int def);
double      parse_double (char *name, double def);
gsl_complex parse_complex(char *name, gsl_complex def);
char       *parse_string (char *name, char *def);

/* Now comes stuff for the blocks */
typedef struct sym_block_line{
  int n;
  char **fields;
} sym_block_line;

typedef struct sym_block{
  int n;
  sym_block_line *lines;
} sym_block;

int parse_block        (char *name, sym_block **blk);
int parse_block_end    (sym_block **blk);
int parse_block_n      (sym_block *blk);
int parse_block_cols   (sym_block *blk, int l);
int parse_block_int    (sym_block *blk, int l, int col, int *r);
int parse_block_double (sym_block *blk, int l, int col, double *r);
int parse_block_complex(sym_block *blk, int l, int col, gsl_complex *r);
int parse_block_string (sym_block *blk, int l, int col, char **r);

/* from parse_exp.c */
enum pr_type {PR_NONE,PR_CMPLX, PR_STR};
typedef struct parse_result{
  union {
    gsl_complex c;
    char *s;
  } value;
  enum pr_type type;
} parse_result;

void parse_result_free(parse_result *t);

int parse_exp(char *exp, parse_result *t);

void parse_putsym_int(char *s, int i);
void parse_putsym_double(char *s, double d);
void parse_putsym_complex(char *s, gsl_complex c);

#endif

max_stars_repo_name: neelravi/octopus

max_issues_repo_name: neelravi/octopus

max_forks_repo_name: neelravi/octopus


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

#define N 8

int
main (void)
{
  int i;
  double xi, yi;
  double x[N] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
  double y[N] = { 0.0, 0.0, 2.0, 0.0,-1.0,-2.0,-1.0, 0.0 };

  printf ("#m=0,S=2\n");

  for (i = 0; i < N; i++)
    {
      printf ("%g %g\n", x[i], y[i]);
    }

  printf ("#m=1,S=0\n");

  {
    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    gsl_spline *spline = gsl_spline_alloc (gsl_interp_akima, N);
    gsl_spline_init (spline, x, y, N);

    for (xi = x[0]; xi < x[N-1]; xi += 0.01)
      {
        double yi = gsl_spline_eval (spline, xi, acc);
        printf ("%g %g\n", xi, yi);
      }
    gsl_spline_free (spline);
    gsl_interp_accel_free(acc);
  }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef CIMPLE_CIMPLE_AUXILIARY_FUNCTIONS_H
#define CIMPLE_CIMPLE_AUXILIARY_FUNCTIONS_H

/**
 * @brief Macro computes number of an array.
 *
 * Remark: It is defined as a Macro so that it can be used independently of the type of the array
 * e.g. int a[2]; =>N_ELEMS(a)=2
 *      polytope b[3]; =>N_ELEMS(a)=2 (this is an array of 3 polytopes that's why it is not "polytope*")
 */
#define N_ELEMS(x)  (sizeof(x) / sizeof((x)[0]))

#include <stdio.h>
#include <string.h>
#include <pthread.h>
#include <stddef.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/**
 * @brief Generates random numbers with normal distribution
 * @param mu
 * @param sigma
 * @return
 */
double randn (double mu,
              double sigma);

/**
 * @brief Timer simulating one time step in discrete control
 * @param arg
 * @return
 */
void * timer(void * arg);

/**
 * @brief Breaks infinte loop in a pthread
 * @param mtx
 * @return
 */
int needQuit(pthread_mutex_t *mtx);
#endif //CIMPLE_CIMPLE_AUXILIARY_FUNCTIONS_H

max_stars_repo_name: shaesaert/TuLiPXML

max_issues_repo_name: shaesaert/TuLiPXML

max_forks_repo_name: shaesaert/TuLiPXML


text:
/* linalg/test_tri.c
 *
 * Copyright (C) 2019 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 <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>

static int
test_symmtd_decomp_eps(const gsl_matrix * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = m->size1;
  size_t i, j;

  gsl_matrix * Q = gsl_matrix_alloc(N, N);
  gsl_matrix * T = gsl_matrix_calloc(N, N);
  gsl_matrix * A  = gsl_matrix_alloc(N, N);
  gsl_matrix * B  = gsl_matrix_alloc(N, N);
  gsl_vector * tau = gsl_vector_alloc(N - 1);
  gsl_vector_view diag = gsl_matrix_diagonal(T);
  gsl_vector_view subdiag = gsl_matrix_subdiagonal(T, 1);
  gsl_vector_view superdiag = gsl_matrix_superdiagonal(T, 1);

  gsl_matrix_memcpy(A, m);
  s += gsl_linalg_symmtd_decomp(A, tau);

  s += gsl_linalg_symmtd_unpack(A, tau, Q, &diag.vector, &subdiag.vector);
  gsl_vector_memcpy(&superdiag.vector, &subdiag.vector);

  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Q, T, 0.0, A); /* A := Q T */
  gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, A, Q, 0.0, B);   /* B := Q T Q^T */

  for (i = 0; i < N; i++)
    {
      for (j = 0; j <= i; j++)
        {
          double bij = gsl_matrix_get(B, i, j);
          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(bij, mij, eps, "%s (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i, j, bij, mij);
        }
    }

  gsl_matrix_free(T);
  gsl_matrix_free(A);
  gsl_matrix_free(Q);
  gsl_matrix_free(B);
  gsl_vector_free(tau);

  return s;
}

static int
test_symmtd_decomp(gsl_rng * r)
{
  int s = 0;
  size_t N;

  for (N = 2; N <= 50; ++N)
    {
      gsl_matrix * A = gsl_matrix_alloc(N, N);

      create_symm_matrix(A, r);
      s += test_symmtd_decomp_eps(A, 1.0e5 * N * GSL_DBL_EPSILON, "symmtd_decomp random");

      gsl_matrix_free(A);
    }

  return s;
}

static int
test_hermtd_decomp_eps(const gsl_matrix_complex * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = m->size1;
  size_t i, j;

  gsl_matrix_complex * Q = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * T = gsl_matrix_complex_calloc(N, N);
  gsl_matrix_complex * A  = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * B  = gsl_matrix_complex_alloc(N, N);
  gsl_vector_complex * tau = gsl_vector_complex_alloc(N - 1);
  gsl_vector_view diag, subdiag, superdiag;
  gsl_vector_complex_view v;

  v = gsl_matrix_complex_diagonal(T);
  diag = gsl_vector_complex_real(&v.vector);

  v = gsl_matrix_complex_subdiagonal(T, 1);
  subdiag = gsl_vector_complex_real(&v.vector);

  v = gsl_matrix_complex_superdiagonal(T, 1);
  superdiag = gsl_vector_complex_real(&v.vector);

  gsl_matrix_complex_memcpy(A, m);
  s += gsl_linalg_hermtd_decomp(A, tau);

  s += gsl_linalg_hermtd_unpack(A, tau, Q, &diag.vector, &subdiag.vector);
  gsl_vector_memcpy(&superdiag.vector, &subdiag.vector);

  gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, Q, T, GSL_COMPLEX_ZERO, A);    /* A := Q T */
  gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, GSL_COMPLEX_ONE, A, Q, GSL_COMPLEX_ZERO, B);  /* B := Q T Q^T */

  for (i = 0; i < N; i++)
    {
      for (j = 0; j <= i; j++)
        {
          gsl_complex bij = gsl_matrix_complex_get(B, i, j);
          gsl_complex mij = gsl_matrix_complex_get(m, i, j);

          gsl_test_rel(GSL_REAL(bij), GSL_REAL(mij), eps, "%s real (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i, j, GSL_REAL(bij), GSL_REAL(mij));
          gsl_test_rel(GSL_IMAG(bij), GSL_IMAG(mij), eps, "%s imag (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i, j, GSL_IMAG(bij), GSL_IMAG(mij));
        }
    }

  gsl_matrix_complex_free(T);
  gsl_matrix_complex_free(A);
  gsl_matrix_complex_free(Q);
  gsl_matrix_complex_free(B);
  gsl_vector_complex_free(tau);

  return s;
}

static int
test_hermtd_decomp(gsl_rng * r)
{
  int s = 0;
  size_t N;

  for (N = 2; N <= 50; ++N)
    {
      gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);

      create_herm_matrix(A, r);
      s += test_hermtd_decomp_eps(A, 1.0e5 * N * GSL_DBL_EPSILON, "hermtd_decomp random");

      gsl_matrix_complex_free(A);
    }

  return s;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdarg.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_mass_functions.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

// Compute cumulative comoving mass function
typedef struct mass_function_cumulative_params_struct mass_function_cumulative_params;
struct mass_function_cumulative_params_struct {
    cosmo_info **cosmo;
    double       z;
    int          mode;
    double *     P;
};
double mass_function_cumulative_integrand(double lM, void *params_in);
double mass_function_cumulative_integrand(double lM, void *params_in) {
    mass_function_cumulative_params *params;
    params       = (mass_function_cumulative_params *)params_in;
    double M     = take_alog10(lM);
    double r_val = mass_function(M, params->z, params->cosmo, params->mode, params->P);
    return (r_val);
}
double mass_function_cumulative(double M_interp, double z, cosmo_info **cosmo, int mode, ...) {
    va_list vargs;
    va_start(vargs, mode);

    // Get passed parameters (if mode specifies they are there).
    double *P = NULL;
    if(SID_CHECK_BITFIELD_SWITCH(mode, MF_PASS_PARAMS))
        P = (double *)va_arg(vargs, double *);

    // Initialize integral
    if(!ADaPS_exist((*cosmo), "lk_P"))
        init_power_spectrum_TF(cosmo);
    double *lk_P     = (double *)ADaPS_fetch((*cosmo), "lk_P");
    double  limit_lo = take_log10(M_interp);
    double  limit_hi = take_log10(M_of_k(take_alog10(lk_P[0]), (*cosmo)));

    // Perform integral
    int    n_int        = 5000;
    double rel_accuracy = 1e-4;
    double r_val        = 0.;
    if(limit_lo < limit_hi) {
        double                          abs_error;
        mass_function_cumulative_params params;
        gsl_function                    integrand;
        gsl_integration_workspace *     wspace;
        params.z           = z;
        params.cosmo       = cosmo;
        params.mode        = mode;
        params.P           = P;
        integrand.function = mass_function_cumulative_integrand;
        integrand.params   = (void *)(&params);
        wspace             = gsl_integration_workspace_alloc(n_int);

        // Integrate mass function
        gsl_integration_qag(&integrand, limit_lo, limit_hi, 0, rel_accuracy, n_int, GSL_INTEG_GAUSS61, wspace, &r_val, &abs_error);

        // Clean-up
        gsl_integration_workspace_free(wspace);
    }

    va_end(vargs);
    return (r_val);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#ifndef DETERMINISTIC_NETWORK_SIR_H
#define DETERMINISTIC_NETWORK_SIR_H

#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include <iostream>
#include <math.h>
#include "DiffEq_Sim.h"

using namespace std;

class Deterministic_Network_SIR_Sim : public DiffEq_Sim {

    private:
        const double r;
        const double mu;
        const vector<double> deg_dist;

    public:
        Deterministic_Network_SIR_Sim() : r(0.0), mu(0.0) { nbins=4;}
        Deterministic_Network_SIR_Sim(double r_param, double mu_param, vector<double> deg_dist_param):
            r(r_param),
            mu(mu_param),
            deg_dist(deg_dist_param) {
                nbins=4;
            }
        ~Deterministic_Network_SIR_Sim() {};

        void initialize( double theta, double pS, double pI, double I) {
            y = new double[nbins];
            y[0] = theta;
            y[1] = pS;
            y[2] = pI;
            y[3] = I;
        }

        double current_susceptible() { return g( y[0] ); }
        double current_infectious() { return y[3]; }
        double current_recovered() { return 1.0 - current_susceptible() - current_infectious(); }

        double g( double theta) {
            double val = 0;
            for (unsigned int i = 0; i<deg_dist.size(); i++) {
                val += deg_dist[i] * pow(theta,i);
            }
            return val;
        }

        double dg( double theta) {
            double val = 0;
            for (unsigned int i = 1; i<deg_dist.size(); i++) {
                val += (i) * deg_dist[i] * pow(theta,i-1);
            }
            return val;
        }

        double ddg( double theta) {
            double val = 0;
            for (unsigned int i = 2; i<deg_dist.size(); i++) {
                val += (i) * (i-1) *deg_dist[i] * pow(theta,i-2);
            }
            return val;
        }

        void derivative(double const y[], double dydt[]) {
            const double theta = y[0];
            const double pS = y[1];
            const double pI = y[2];
            const double I = y[3];

            dydt[0] = -r * pI * theta;                                                  // dtheta.dt
            dydt[1] = r *  pS * pI * (1 - theta * ddg(theta)/dg(theta));                // dpS.dt
            dydt[2] = r *  pI * pS * theta * ddg(theta)/dg(theta) - pI*(1-pI)*r- pI*mu; // dpI.dt

            dydt[3] = r * pI * theta * dg(theta) - mu*I;
        }

};

#endif


text:
/****************************************************************************
 *                                                                          *
 *  Author : lukasz.iwaszkiewicz@gmail.com                                  *
 *  ~~~~~~~~                                                                *
 *  License : see COPYING file for details.                                 *
 *  ~~~~~~~~~                                                               *
 ****************************************************************************/

#pragma once
#include "Condition.h"
#include <gsl/gsl>

/**
 * Sequence is like AndCondition, but inputs must fulfill conditions in order.
 */
template <typename EventT = LIB_STATE_MACHINE_DEFAULT_EVENT_TYPE> class SequenceCondition : public Condition<EventT> {
public:
        using EventType = EventT;

        SequenceCondition (Condition<EventT> *a, Condition<EventT> *b) : a (a), b (b) {}
        virtual ~SequenceCondition () = default;

        bool getResult () const override { return a->getResult () && b->getResult (); }
        void reset () override
        {
                a->reset ();
                b->reset ();
        }

        bool check (EventType  &event, EventType &retainedEvent) const override
        {
                if (!a->getResult ()) {
                        a->check (event, retainedEvent);
                }

                if (a->getResult () && !b->getResult ()) {
                        b->check (event, retainedEvent);
                }

                return getResult ();
        }

private:
        gsl::not_null <Condition<EventType>*> a;
        gsl::not_null <Condition<EventType>*> b;
};

template <typename EventT = LIB_STATE_MACHINE_DEFAULT_EVENT_TYPE> SequenceCondition<EventT> *seq (Condition<EventT> *a, Condition<EventT> *b)
{
        return new SequenceCondition<EventT> (a, b);
}

max_stars_repo_name: iwasz/libstatemachine

max_issues_repo_name: iwasz/libstatemachine

max_forks_repo_name: iwasz/libstatemachine