text:
/* fft/gsl_fft_halfcomplex.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 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 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.
 */

#ifndef __GSL_FFT_HALFCOMPLEX_H__
#define __GSL_FFT_HALFCOMPLEX_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stddef.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_fft.h>
#include <gsl/gsl_fft_real.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

GSL_FUN int gsl_fft_halfcomplex_radix2_backward (double data[], const size_t stride, const size_t n);
GSL_FUN int gsl_fft_halfcomplex_radix2_inverse (double data[], const size_t stride, const size_t n);
GSL_FUN int gsl_fft_halfcomplex_radix2_transform (double data[], const size_t stride, const size_t n);

typedef struct
  {
    size_t n;
    size_t nf;
    size_t factor[64];
    gsl_complex *twiddle[64];
    gsl_complex *trig;
  }
gsl_fft_halfcomplex_wavetable;

GSL_FUN gsl_fft_halfcomplex_wavetable * gsl_fft_halfcomplex_wavetable_alloc (size_t n);

GSL_FUN void
gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * wavetable);

GSL_FUN int gsl_fft_halfcomplex_backward (double data[], const size_t stride, const size_t n,
                                  const gsl_fft_halfcomplex_wavetable * wavetable,
                                  gsl_fft_real_workspace * work);

GSL_FUN int gsl_fft_halfcomplex_inverse (double data[], const size_t stride, const size_t n,
                                 const gsl_fft_halfcomplex_wavetable * wavetable,
                                 gsl_fft_real_workspace * work);

GSL_FUN int gsl_fft_halfcomplex_transform (double data[], const size_t stride, const size_t n,
                                   const gsl_fft_halfcomplex_wavetable * wavetable,
                                   gsl_fft_real_workspace * work);

GSL_FUN int
gsl_fft_halfcomplex_unpack (const double halfcomplex_coefficient[],
                            double complex_coefficient[],
                            const size_t stride, const size_t n);

GSL_FUN int
gsl_fft_halfcomplex_radix2_unpack (const double halfcomplex_coefficient[],
                                   double complex_coefficient[],
                                   const size_t stride, const size_t n);

__END_DECLS

#endif /* __GSL_FFT_HALFCOMPLEX_H__ */

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:
#pragma once

#include <random>
#include <type_traits>

#include <gsl-lite/gsl-lite.hpp>

namespace thrustshift {

namespace random {

template <class Range, class Distribution, class Engine>
void generate(Range&& range, Distribution& dist, Engine& rng) {
	for (size_t i = 0; i < range.size(); ++i) {
		range[i] = dist(rng);
	}
}

template <class Range,
          typename T = typename std::remove_reference<Range>::type::value_type>
void generate_uniform_real(Range&& range,
                           T min = T{0},
                           T max = T{1},
                           unsigned long long seed = 0) {
	gsl_Expects(min <= max);
	std::default_random_engine rng(seed);
	std::uniform_real_distribution<T> dist(min, max);
	generate(std::forward<Range>(range), dist, rng);
}

template <class Range,
          typename T = typename std::remove_reference<
              Range>::type::value_type::value_type>
void generate_uniform_complex(Range&& range,
                              T min = T{0},
                              T max = T{1},
                              unsigned long long seed = 0) {
	gsl_Expects(min <= max);
	using Complex = typename std::remove_reference<Range>::type::value_type;
	std::default_random_engine rng(seed);
	std::uniform_real_distribution<T> dist(min, max);
	auto cdist = [&](auto& rng) -> Complex { return {dist(rng), dist(rng)}; };
	generate(std::forward<Range>(range), cdist, rng);
}

template <class Range,
          typename T = typename std::remove_reference<Range>::type::value_type>
void generate_normal_real(Range&& range,
                          T mean = T{0},
                          T stddev = T{1},
                          unsigned long long seed = 0) {
	gsl_Expects(stddev >= 0);
	std::default_random_engine rng(seed);
	std::normal_distribution<T> dist(mean, stddev);
	generate(std::forward<Range>(range), dist, rng);
}

template <class Range,
          typename T = typename std::remove_reference<
              Range>::type::value_type::value_type>
void generate_normal_complex(Range&& range,
                             T mean = T{0},
                             T stddev = T{1},
                             unsigned long long seed = 0) {
	gsl_Expects(stddev >= 0);
	using Complex = typename std::remove_reference<Range>::type::value_type;
	std::default_random_engine rng(seed);
	std::normal_distribution<T> dist(mean, stddev);
	auto cdist = [&](auto& rng) -> Complex { return {dist(rng), dist(rng)}; };
	generate(std::forward<Range>(range), cdist, rng);
}

} // namespace random

} // namespace thrustshift

ext: h

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/** @file kjg_fpca.h
 * @brief Runs fastPCA.
 * This module also has methods to multiply a genotype matrix against the GSL
 * matrices.
 */

#ifndef KJG_FPCA_H_
#define KJG_FPCA_H_

#include <gsl/gsl_matrix.h>

extern size_t KJG_FPCA_ROWS;	// number of rows to process at once

/** Performs a fast PCA
 * @param *eval eigenvalues
 * @param *evec eigenvectors
 * @param K number of eigenvalues/vectors
 * @param L width of projection matrix
 * @param I iterations to do exponentiation
 */

void kjg_fpca (size_t K, size_t L, size_t I, double *eval, double *evec);

/** Multiplies B=X*A1 and A2 = XT*B = XT*X*A1
 * @param *A1 some matrix
 * @param *B intermediate matrix
 * @param *A2 next matrix
 */

void kjg_fpca_XTXA (const gsl_matrix * A1, gsl_matrix * B, gsl_matrix * A2);

/** Multiplies B = X*A
 * @param *A some matrix
 * @param *B another matrix
 */

void kjg_fpca_XA (const gsl_matrix * A, gsl_matrix * B);

/** Multiplies A = XT*B
 * @param *B some matrix
 * @param *A another matrix
 */

void kjg_fpca_XTB (const gsl_matrix * B, gsl_matrix * A);

#endif /* KJG_FPCA_H_ */

max_stars_repo_name: nevrome/EIG

max_issues_repo_name: nevrome/EIG

max_forks_repo_name: nevrome/EIG


text:
/* randist/hyperg.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

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

/* The hypergeometric distribution has the form,

   prob(k) =  choose(n1,t) choose(n2, t-k) / choose(n1+n2,t)

   where choose(a,b) = a!/(b!(a-b)!)

   n1 + n2 is the total population (tagged plus untagged)
   n1 is the tagged population
   t is the number of samples taken (without replacement)
   k is the number of tagged samples found

*/

unsigned int
gsl_ran_hypergeometric (const gsl_rng * r, unsigned int n1, unsigned int n2,
                        unsigned int t)
{
  const unsigned int n = n1 + n2;

  unsigned int i = 0;
  unsigned int a = n1;
  unsigned int b = n1 + n2;
  unsigned int k = 0;

  if (t > n)
    {
      t = n ;
    }

  if (t < n / 2)
    {
      for (i = 0 ; i < t ; i++)
        {
          double u = gsl_rng_uniform(r) ;

          if (b * u < a)
            {
              k++ ;
              if (k == n1)
                return k ;
              a-- ;
            }
          b-- ;
        }
      return k;
    }
  else
    {
      for (i = 0 ; i < n - t ; i++)
        {
          double u = gsl_rng_uniform(r) ;

          if (b * u < a)
            {
              k++ ;
              if (k == n1)
                return n1 - k ;
              a-- ;
            }
          b-- ;
        }
      return n1 - k;
    }

}

double
gsl_ran_hypergeometric_pdf (const unsigned int k,
                            const unsigned int n1,
                            const unsigned int n2,
                            unsigned int t)
{
  if (t > n1 + n2)
    {
      t = n1 + n2 ;
    }

  if (k > n1 || k > t)
    {
      return 0 ;
    }
  else if (t > n2 && k + n2 < t )
    {
      return 0 ;
    }
  else
    {
      double p;

      double c1 = gsl_sf_lnchoose(n1,k);
      double c2 = gsl_sf_lnchoose(n2,t-k);
      double c3 = gsl_sf_lnchoose(n1+n2,t);

      p = exp(c1 + c2 - c3) ;

      return p;
    }
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
////////////////////////////////////////////////////////////
//
// Copyright (c) 2018 Jan Filipowicz, Filip Turobos
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
////////////////////////////////////////////////////////////

#ifndef GENETIC_ALGORITHM_LIBRARY_EVALUATED_SPECIMEN_H
#define GENETIC_ALGORITHM_LIBRARY_EVALUATED_SPECIMEN_H

#include <optional>
#include <utility>
#include <tuple>
#include <type_traits>
#include <gsl/gsl_assert>

template<class Specimen, class Rating>
class evaluated_specimen {
public:
	using value_type = Specimen;
	using rating_type = Rating;
	template<class... Args, class = std::enable_if_t<!std::is_same_v<std::tuple<std::decay_t<Args>...>, std::tuple<evaluated_specimen>>>>
	constexpr evaluated_specimen(Args&&... args) noexcept(sizeof...(Args) == 0);
	constexpr value_type& value() & noexcept;
	constexpr const value_type& value() const& noexcept;
	constexpr bool has_rating() const noexcept;
	constexpr rating_type rating() const;
	template<class Function>
	void evaluate(Function&& evaluator);
private:
	value_type specimen;
	std::optional<rating_type> grade;
};

template<class Specimen, class Rating>
template<class... Args, class>
constexpr evaluated_specimen<Specimen, Rating>::evaluated_specimen(Args&&... args) noexcept(sizeof...(Args) == 0)
	: specimen(std::forward<Args>(args)...) {}

template<class Specimen, class Rating>
constexpr auto evaluated_specimen<Specimen, Rating>::value() & noexcept -> value_type& {
	return specimen;
}

template<class Specimen, class Rating>
constexpr auto evaluated_specimen<Specimen, Rating>::value() const& noexcept -> const value_type& {
	return specimen;
}

template<class Specimen, class Rating>
constexpr bool evaluated_specimen<Specimen, Rating>::has_rating() const noexcept {
	return grade.has_value();
}

template<class Specimen, class Rating>
constexpr auto evaluated_specimen<Specimen, Rating>::rating() const -> rating_type {
	return grade.value();
}

template<class Specimen, class Rating>
template<class Function>
inline void evaluated_specimen<Specimen, Rating>::evaluate(Function&& evaluator) {
	grade.emplace(evaluator(std::as_const(specimen)));
	Ensures(has_rating());
}

#endif

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


text:
#pragma once

#include "BasicBox.h"
#include "BasicTypes.h"
#include "BoxTypes.h"
#include "General.h"
#include "NodeDim.h"
#include "NodeListDim.h"
#include <gsl.h>

namespace formulae{ namespace node {

/**
 *
 * @param width
 * @param node center node -> reference point, (must be Box0 or Rule !!!)
 * @param upList above center, enumerated from bottom to top
 * @param dnList below center, enumerated from top to bottom
 * @return
 */
template<typename VList1=vlist_t, typename VList2=vlist_t, typename Node=NodeRef>
box_t makeVBox(dist_t  width, Node && node,
               VList1 && upList, VList2 && dnList) {

  using general::revAppend;
  using general::transfer;

  auto h = vlistVsize(upList) + node.height();
  auto d = vlistVsize(dnList) + node.depth();
  std::remove_reference_t<VList1> emptyList;
  auto nodeList = revAppend(std::forward<VList1>(upList), emptyList);
  nodeList.push_back(std::forward<Node>(node));
  transfer(nodeList, std::forward<VList2>(dnList));

  return vbox(dim_t{width, h, d}, std::move(nodeList));
}
/**
 *
 * @param width
 * @param center center node -> reference point,
 * @param upList above center, enumerated from bottom to top
 * @return
 */
template <typename Box = box_t, typename VList = vlist_t>
box_t upVBox(dist_t  width, Box && center, VList&& upList) {
  std::remove_reference_t<VList> emptyList;
  return makeVBox(width, NodeRef{Box0(std::forward<Box>(center))},
                  std::forward<VList>(upList), emptyList);
}

/**
 *
 * @param width
 * @param center center node -> reference point
 * @param dnList below center, enumerated from top to bottom
 * @return
 */
template <typename Box = box_t, typename VList = vlist_t>
box_t dnVBox(dist_t  width, Box && center, VList && dnList) {
  std::remove_reference_t<VList> emptyList;
  return makeVBox(width, Box0(std::forward<Box>(center)), emptyList, std::forward<VList>(dnList));
}

/**
 *
 * @param n1
 * @param dist1 distance from bottom of n1 to baseline
 * @param dist  vertical distance between n1 and n2
 * @param n2
 * @return
 */
template <typename Node = NodeRef, typename DownNode=NodeRef>
NodeRef above(Node && n1, dist_t  dist1, dist_t  dist, DownNode && n2) {
  auto w = std::max(n1.vwidth(), n2.vwidth());
  auto h = n1.vsize() + dist1;
  auto d = n2.vsize() + dist - dist1;
  using E = std::remove_reference_t<Node>;
  list_t<E> nodeList;
  nodeList.push_back(std::forward<Node>(n1));
  nodeList.push_back(Kern_t::create(dist));
  nodeList.push_back(std::forward<DownNode>(n2));
  return Box0(vbox(dim_t{w, h, d}, std::move(nodeList)));
};

}
}


text:
#ifndef __COMPLEX_ARRAY_H__
#define __COMPLEX_ARRAY_H__

#include <complex>
#include <vector>
#include <iostream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>

#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
namespace mspass::algorithms::deconvolution{
typedef std::complex<double> Complex64;
typedef std::complex<float> Complex32;

typedef struct FortranComplex32 {
    float real;
    float imag;
} FortranComplex32;
typedef struct FortranComplex64 {
    double real;
    double imag;
} FortranComplex64;

/* \brief Interfacing object to ease conversion between FORTRAN and C++ complex.

   */
class ComplexArray
{
public:
    /*! Empty constructor. */
    ComplexArray();
    /*! Construct from stl vector container of complex. */
    ComplexArray(std::vector<Complex64> &d);
    /*! Similar for 32 bit version */
    ComplexArray(std::vector<Complex32> &d);
    /*! Construct from a FORTRAN data array.

      Fortran stores complex numbers in a mulitplexed array
      structure (real(1), imag(1), real(2), imag(2), etc.).
      The constructors below provide a mechanism for building
      this object from various permutations of this.
      \param nsamp is the number of elements in the C vector
      \param d is the pointer to the first compoment of the
         fortran vector.
      */
    ComplexArray(int nsamp, FortranComplex32 *d);
    ComplexArray(int nsamp, FortranComplex64 *d);
    ComplexArray(int nsamp, float *d);
    ComplexArray(int nsamp, double *d);
    ComplexArray(int nsamp);
    /*! Construct from different length of vector, adds zoeros to it
    	And construct a constant arrays
    	*/
    template<class T> ComplexArray(int nsamp, std::vector<T> d);
    template<class T> ComplexArray(int nsamp, T d);

    /*! Construct from magnitude and phase arrays.*/
    ComplexArray(std::vector<double> mag,std::vector<double> phase);

    /* These will need to be implemented.  Likely cannot
       depend on the compiler to generate them correctly */
    ComplexArray(const ComplexArray &parent);
    ComplexArray& operator=(const ComplexArray &parent);
    ~ComplexArray();

    /* These are kind of the inverse of the constructor.
    Independent of what the internal representation is they
    will return useful interface representations. */
    /*! Return a pointer to a fortran array containing
      the data vector.

      The array linked to the returned pointer should be
      created with the new operator and the caller should
      be sure to use delete [] to free this memory when
      finished. */
    template<class T> T *FortranData();
    /* This is same for what I think fortran calls
       double complex */
//        double *FortranData();
    /* C representation.  This can be templated easily.
    See below.  The syntax is weird and should probably
    be wrapped with a typedef */
    template<class T> std::vector<std::complex<T> > CPPData();

    /* Operators are the most important elements of this
       thing to make life easier. */
    /*! Index operator.
    	Cannot make it work by getting the address from reference.
    	Have to call the ptr() function to get the address.
      \param sample is the sample number to return.
      \return contents of vector at position sample.
      */
    Complex64 operator[](int sample);
    double *ptr();
    double *ptr(int sample);
    ComplexArray& operator +=(const ComplexArray& other) noexcept(false);
    ComplexArray& operator -=(const ComplexArray& other) noexcept(false);
    /* This actually is like .* in matlab - sample by sample multiply not
    a dot product */
    ComplexArray& operator *= (const ComplexArray& other)noexcept(false);
    /* This is like *= but complex divide element by element */
    ComplexArray& operator /= (const ComplexArray& other)noexcept(false);
    const ComplexArray operator +(const ComplexArray& other)const noexcept(false);
    //template<class T> ComplexArray operator +(const vector<T> &other);
    //template<class T> ComplexArray operator +(const T &other);
    const ComplexArray operator -(const ComplexArray& other)const noexcept(false);
    //template<class T> ComplexArray operator -(const vector<T> &other);
    //template<class T> ComplexArray operator -(const T &other);
    const ComplexArray operator *(const ComplexArray& other)const noexcept(false);
    const ComplexArray operator /(const ComplexArray& other)const noexcept(false);
    /*! product of complex and real vectors */
    //template<class T> ComplexArray operator *(const vector<T> &other);
    //template<class T> friend ComplexArray operator *(const vector<T> &lhs,const ComplexArray &rhs);
    /*! product of complex and a number */
    //template<class T> ComplexArray operator *(const T &other);
    //template<class T> friend ComplexArray operator *(const T &lhs,const ComplexArray &rhs);
    /*! Change vector to complex conjugates. */
    void conj();
    /* Return stl vector of amplitude spectrum.  */
    std::vector<double> abs() const;
    /* Return rms value.  */
    double rms() const;
    /* Return 2-norm value.  */
    double norm2() const;
    /* Return stl vector of phase */
    std::vector<double> phase() const;
    /* Return size of the array*/
    int size() const;
private:
    /* Here is an implementation detail.   There are three ways
       I can think to do this.  First, we could internally store
       data as fortran array of 32 bit floats.   That is probably
       the best because we can use BLAS routines (if you haven't
       heard of this - likely - I need to educate you.)  to do
       most of the numerics fast. Second, we could use stl
       vector container of std::complex.  The third is excessively
       messy but technically feasible - I would not recommend it.
       That is, one could store pointers to either representation
       and internally convert back and forth.  Ugly and dangerous
       I think.

       I suggest we store a FORTRAN 32 bit form since that is
       what standard numeric libraries (e.g. most fft routines)
       use.  */
    /*I decided to use 64 bit, since the GSL's fft routine is using that.*/
    FortranComplex64 *data;
    int nsamp;
};
/* This would normally be in the .h file and since I don't think
   you've used templates worth showing you how it would work. */
template <class T> std::vector<std::complex<T> > ComplexArray::CPPData()
{
    std::vector<std::complex<T> > result;
    result.reserve(nsamp);
    int i;
    for(i=0; i<nsamp; ++i)
    {
        std::complex<T> z(data[i].real, data[i].imag);
        result.push_back(z);
    }
    return result;
}

template<class T> T* ComplexArray::FortranData()
{
    T* result=new T[nsamp];
    for(int i=0; i<nsamp; i++)
        result[i]=data[i];
    return result;
}

template<class T> ComplexArray::ComplexArray(int n, std::vector<T> d)
{
    nsamp=n;
    if(nsamp>d.size())
    {
        data=new FortranComplex64[nsamp];
        for(int i=0; i<d.size(); i++)
        {
            data[i].real=d[i];
            data[i].imag=0.0;
        }
        for(int i=d.size(); i<nsamp; i++)
        {
            data[i].real=0.0;
            data[i].imag=0.0;
        }
    }
    else
    {
        data=new FortranComplex64[nsamp];
        for(int i=0; i<nsamp; i++)
        {
            data[i].real=d[i];
            data[i].imag=0.0;
        }
    }
}
template<class T> ComplexArray::ComplexArray(int n, T d)
{
    nsamp=n;
    data=new FortranComplex64[nsamp];
    for(int i=0; i<nsamp; i++)
    {
        data[i].real=d;
        data[i].imag=0.0;
    }
}
/*
template<class T> ComplexArray ComplexArray::operator +(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real+other[i];
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator +(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real+other;
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator -(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real-other[i];
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator -(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real-other;
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator *(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real*other[i];
        result.data[i].imag=data[i].imag*other[i];
    }
    return result;
}
template<class T> ComplexArray operator *(const vector<T>& lhs,const ComplexArray& rhs)
{
    return rhs*lhs;
}
template<class T> ComplexArray ComplexArray::operator *(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real*other;
        result.data[i].imag=data[i].imag*other;
    }
    return result;
}
template<class T> ComplexArray operator *(const T& lhs,const ComplexArray& rhs)
{
    return rhs*lhs;
}
*/
}
#endif

max_stars_repo_name: Yangzhengtang/mspass

max_issues_repo_name: Yangzhengtang/mspass

max_forks_repo_name: Yangzhengtang/mspass


text:
#pragma once
#include <SDL2/SDL.h>
#include <SDL2/SDL_audio.h>
#include <chrono>
#include <cstdint>
#include <gsl/span>
#include <queue>

class CSound
{
private:
	static constexpr std::int32_t Frequency{ 44100 };

	struct SBeep
	{
		double Frequency;
		std::uint32_t SamplesLeft;
	};

	SDL_AudioDeviceID mDeviceID;
	std::queue<SBeep> mBeeps;
	double mV;

public:
	CSound();
	~CSound();

	void Beep(double frequency, std::chrono::milliseconds duration);

private:
	void GenerateSamples(gsl::span<std::int16_t> stream);

	static void SdlCallback(CSound* self, std::uint8_t* stream, std::int32_t len);
};


text:
/* Author: G. Jungman
 */
#include <config.h>
#include <gsl/gsl_ieee_utils.h>

#include "gsl_qrng.h"

int test_sobol(void)
{
  int status = 0;
  double v[3];
  /* int i; */

  /* test in dimension 2 */
  gsl_qrng * g = gsl_qrng_alloc(gsl_qrng_sobol, 2);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.25 || v[1] != 0.75 );
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.375 || v[1] != 0.375 );
  gsl_qrng_free(g);

  /* test in dimension 3 */
  g = gsl_qrng_alloc(gsl_qrng_sobol, 3);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.25 );
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.375 || v[1] != 0.375 || v[2] != 0.625 );
  gsl_qrng_free(g);

  return status;
}

int test_nied2(void)
{
  int status = 0;
  double v[3];
  /* int i; */

  /* test in dimension 2 */
  gsl_qrng * g = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 2);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.75 || v[1] != 0.25 );
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.25 || v[1] != 0.75 );
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.625 || v[1] != 0.125 );
  gsl_qrng_free(g);

  /* test in dimension 3 */
  g = gsl_qrng_alloc(gsl_qrng_niederreiter_2, 3);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.75 || v[1] != 0.25 || v[2] != 0.3125 );
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.25 || v[1] != 0.75 || v[2] != 0.5625 );
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  gsl_qrng_get(g, v);
  status += ( v[0] != 0.625 || v[1] != 0.125 || v[2] != 0.6875 );
  gsl_qrng_free(g);

  return status;
}

int main()
{
  int status = 0;

  gsl_ieee_env_setup ();

  status += test_sobol();
  status += test_nied2();

  return status;
}


text:
/*
 * Copyright 2014 Marc Normandin
 *
 * 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 INC_PARTICLE_H
#define INC_PARTICLE_H

#include <iostream>
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
#include <cstdlib>
#include <gsl/gsl_rng.h>
#include <limits>
#include <fstream>
#include <cassert>

#include "rng.h"
#include "dim.h"

class Particle
{
public:
    typedef std::vector<dim_t> dvector;

    Particle();

    Particle(const std::vector<Dim>& dim);

    Particle(const dvector& pos, const prob_t fitness = -1.0 * std::numeric_limits<dim_t>::max());

    const dvector& getPosition() const;

    dvector getVelocity() const;

    double getFitness() const;

    dvector getPBestPosition() const;

    double getPBestFitness() const;

    void updateFitness(const prob_t newFitness);

    void updatePosition(const Particle& GBest, const std::vector<Dim>& dim, const double C1, const double C2,
                        const RandomNumberGenerator& rng, const double inertiaWeight);

    size_t size() const;

    friend bool operator<(const Particle& p1, const Particle& p2);

private:
    dvector     mPos;
    dvector     mVel;
    prob_t      mFitness;

    // The particles best position
    dvector		mPBestPos;
    prob_t 		mPBestFitness;
};

bool operator<(const Particle& p1, const Particle& p2);

#endif // #ifndef INC_PARTICLE_H

max_stars_repo_name: marcnormandin/ParticleSwarmOptimization

max_issues_repo_name: marcnormandin/ParticleSwarmOptimization

max_forks_repo_name: marcnormandin/ParticleSwarmOptimization


text:
#include <gsl/gsl_qrng.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char* argv[]) {
    if (argc < 3) {
        printf("Usage: ./sampling <dimensions> <samples>\n");
        return 0;
    }

    int dimensions = atoi(argv[1]);
    int samples = atoi(argv[2]);

    if (samples < 0) {
        printf("Invalid number of samples: %d\n", samples);
        return 0;
    }

    gsl_qrng * q;
    if (dimensions <= 0) {
        printf("Invalid number of dimensions\n");
        return 1;
    }
    else if (dimensions <= 12) q = gsl_qrng_alloc (gsl_qrng_niederreiter_2, dimensions);
    else if (dimensions <= 40) q = gsl_qrng_alloc (gsl_qrng_sobol, dimensions);
    else if (dimensions <= 1229) q = gsl_qrng_alloc (gsl_qrng_halton, dimensions);
    else {
        printf("Invalid number of dimensions (max 1229 supported): dimensions\n", dimensions);
        return 1;
    }

    for (int i = 0; i < samples; i++) {
        double v[dimensions];
        gsl_qrng_get(q, v);
        for (int d = 0; d < dimensions; d++) {
            printf("%f ", v[d]);
        }
        printf ("\n");
    }
    gsl_qrng_free(q);

    return 0;
}


text:
/*
FRACTAL - A program growing fractals to benchmark parallelization and drawing
libraries.

Copyright 2009-2021, 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 fractal.c
 * \brief Source file to define the fractal data and functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2009-2021, Javier Burguete Tolosa.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#if HAVE_SYSINFO
#include <sys/sysinfo.h>
#endif
#include <gsl/gsl_rng.h>
#include <libxml/parser.h>
#include <glib.h>
#if HAVE_GTOP
#include <glibtop.h>
#include <glibtop/close.h>
#endif
#ifdef G_OS_WIN32
#include <windows.h>
#endif
#include <libintl.h>
#include <GL/glew.h>
#include <ft2build.h>
#include FT_FREETYPE_H
#include <gtk/gtk.h>
#include "config.h"
#include "fractal.h"
#include "simulator.h"

extern void draw ();

unsigned int width = WIDTH;     ///< Medium width.
unsigned int height = HEIGHT;   ///< Medium height.
unsigned int length = LENGTH;   ///< Medium length.
unsigned int area;              ///< Medium area.
unsigned int medium_bytes;      ///< Number of bytes used by the medium.
unsigned int breaking = 0;      ///< 1 on breaking, 0 otherwise.
unsigned int simulating = 0;    ///< 1 on simulating, 0 otherwise.
unsigned int animating = 1;     ///< 1 on animating, 0 otherwise.

unsigned int fractal_type = FRACTAL_TYPE_TREE;  ///< Fractal type.
unsigned int fractal_3D = 0;    ///< 1 on 3D fractals, 0 on 2D fractals.
unsigned int fractal_diagonal = 0;
///< 1 on diagonal point movement, 0 otherwise.

unsigned long t0;               ///< Computational time.

unsigned int max_d = 0;         ///< Maximum fractal size.
unsigned int *medium = NULL;    ///< Array of fractal points.
Point3D *point = NULL;          ///< Array of 3D points.
unsigned int npoints = 0;       ///< Number of points.

unsigned int random_algorithm = 0;
///< Type of random numbers generator algorithm.
unsigned int random_seed_type = 1;      ///< Type of random seed.
unsigned long random_seed = SEED;       ///< Random seed.
static void *(*parallel_fractal) (gsl_rng * rng);
///< Pointer to the function to calculate the fractal.

static const float color3f[16][3] = {
  {0., 0., 0.},
  {1., 0., 0.},
  {0., 1., 0.},
  {0., 0., 1.},
  {0.5, 0.5, 0.},
  {0.5, 0., 0.5},
  {0., 0.5, 0.5},
  {0.5, 0.25, 0.25},
  {0.25, 0.5, 0.25},
  {0.25, 0.25, 0.5},
  {0.75, 0.25, 0.},
  {0.75, 0., 0.25},
  {0.25, 0.75, 0.},
  {0., 0.75, 0.25},
  {0.25, 0., 0.75},
  {0., 0.25, 0.75}
};                              ///< Array of colors.

DialogOptions dialog_options[1];
  ///< DialogOptions to set the fractal options.
DialogSimulator dialog_simulator[1];
  ///< DialogSimulator to show the main window.

// PARALLELIZING DATA

unsigned int nthreads;          ///< Threads number.
GMutex mutex[1];                ///< Mutex to lock memory saves.

// END

/**
 * Function to return the square of an unsigned int.

 * \return square.
 */
static inline unsigned int
sqr (int x)                     ///< unsigned int.
{
  return x * x;
}

/**
 * Function to get the number of cores of the CPU.
 *
 * \return CPU number of cores.
 */
int
threads_number ()
{
#if HAVE_GTOP
  int ncores;
  glibtop *top;
  top = glibtop_init ();
  ncores = top->ncpu + 1;
  glibtop_close ();
  return ncores;
#elif HAVE_GET_NPROCS
  return get_nprocs ();
#elif defined(G_OS_WIN32)
  SYSTEM_INFO sysinfo;
  GetSystemInfo (&sysinfo);
  return sysinfo.dwNumberOfProcessors;
#else
  return (int) sysconf (_SC_NPROCESSORS_ONLN);
#endif
}

/**
 * Function to add a point to the array.
 */
static inline void
points_add (int x,              ///< Point x-coordinate.
            int y,              ///< Point y-coordinate.
            int z,              ///< Point z-coordinate.
            unsigned int c)     ///< Point color.
{
  Point3D *p;
  ++npoints;
  point = (Point3D *) g_realloc (point, npoints * sizeof (Point3D));
  p = point + npoints - 1;
  p->r[0] = x;
  p->r[1] = y;
  p->r[2] = z;
  memcpy (p->c, color3f[c], 3 * sizeof (float));
}

/**
 * Function to make a random 2D movement on a point.
 */
static inline void
point_2D_move (int *x,          ///< Point x-coordinate.
               int *y,          ///< Point y-coordinate.
               gsl_rng * rng)   ///< Pseudo-random number generator.
{
  register unsigned int k;
  static const int mx[4] = { 0, 0, 1, -1 }, my[4] = { 1, -1, 0, 0 };
  k = gsl_rng_uniform_int (rng, 4);
  *x += mx[k];
  *y += my[k];
}

/**
 * Function to make a random 2D movement on a point enabling diagonals.
 */
static inline void
point_2D_move_diagonal (int *x, ///< Point x-coordinate.
                        int *y, ///< Point y-coordinate.
                        gsl_rng * rng)  ///< Pseudo-random number generator.
{
  register unsigned int k;
  static const int mx[8] = { 1, 1, 1, 0, -1, -1, -1, 0 },
    my[8] = { 1, 0, -1, -1, -1, 0, 1, 1 };
  k = gsl_rng_uniform_int (rng, 8);
  *x += mx[k];
  *y += my[k];
}

/**
 * Function to make a random 3D movement on a point.
 */
static inline void
point_3D_move (int *x,          ///< Point x-coordinate.
               int *y,          ///< Point y-coordinate.
               int *z,          ///< Point z-coordinate.
               gsl_rng * rng)   ///< Pseudo-random number generator.
{
  register unsigned int k;
  static const int mx[6] = { 0, 1, -1, 0, 0, 0 },
    my[6] = { 0, 0, 0, 1, -1, 0 }, mz[6] = { 1, 0, 0, 0, 0, -1 };
  k = gsl_rng_uniform_int (rng, 6);
  *x += mx[k];
  *y += my[k];
  *z += mz[k];
}

/**
 * Function to make a random 3D movement on a point enabling diagonals.
 */
static inline void
point_3D_move_diagonal (int *x, ///< Point x-coordinate.
                        int *y, ///< Point y-coordinate.
                        int *z, ///< Point z-coordinate.
                        gsl_rng * rng)  ///< Pseudo-random number generator.
{
  register int k;
  static const int mx[26] = {
    1, 1, 1, 0, 0, 0, -1, -1, -1,
    1, 1, 1, 0, 0, -1, -1, -1,
    1, 1, 1, 0, 0, 0, -1, -1, -1
  }, my[26] = {
    1, 0, -1, 1, 0, -1, 1, 0, -1,
    1, 0, -1, 1, -1, 1, 0, -1, 1, 0, -1, 1, 0, -1, 1, 0, -1
  }, mz[26] = {
    1, 1, 1, 1, 1, 1, 1, 1, 1,
    0, 0, 0, 0, 0, 0, 0, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1
  };
  k = gsl_rng_uniform_int (rng, 26);
  *x += mx[k];
  *y += my[k];
  *z += mz[k];
}

/**
 * Function to start a new 2D tree point.
 */
static inline void
tree_2D_point_new (int *x,      ///< Point x-coordinate.
                   int *y,      ///< Point y-coordinate.
                   gsl_rng * rng)       ///< Pseudo-random number generator.
{
  *x = gsl_rng_uniform_int (rng, width);
  *y = max_d;
#if DEBUG
  printf ("New point x %d y %d\n", *x, *y);
#endif
}

/**
 * Function to check the limits of a 2D tree point.
 */
static inline void
tree_2D_point_boundary (int *x, ///< Point x-coordinate.
                        int *y, ///< Point y-coordinate.
                        gsl_rng * rng)  ///< Pseudo-random number generator.
{
  if (*y < 0 || *y == (int) height)
    {
      tree_2D_point_new (x, y, rng);
      return;
    }
  if (*x < 0)
    *x = width - 1;
  else if (*x == (int) width)
    *x = 0;
#if DEBUG
  printf ("Boundary point x %d y %d\n", *x, *y);
#endif
}

/**
 * Function to fix a 2D tree point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
tree_2D_point_fix (int x,       ///< Point x-coordinate.
                   int y)       ///< Point y-coordinate.
{
  register unsigned int *point;
#if DEBUG
  printf ("x=%d y=%d max_d=%d width=%d height=%d\n", x, y, max_d, width,
          height);
#endif
  if (y > (int) max_d || x == 0 || y == 0 || x == (int) width - 1
      || y == (int) height - 1)
    return 0;
  point = medium + y * width + x;
  if (point[1] || point[-1] || point[width] || point[-(int) width])
    {
#if DEBUG
      printf ("fixing point\n");
#endif
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      point[0] = 2;
      points_add (x, y, 0, 2);
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to init a 2D tree.
 */
static inline void
tree_2D_init ()
{
  medium[width / 2] = 2;
}

/**
 * Function to check the end a 2D tree.
 *
 * \return 1 on ending, 0 on continuing.
 */
static inline unsigned int
tree_2D_end (int y)             ///< Point y-coordinate.
{
  if (y == (int) max_d)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      ++max_d;
      g_mutex_unlock (mutex);
// END
    }
  if (max_d >= height - 1)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      max_d = height - 1;
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to start a new 3D tree point.
 */
static inline void
tree_3D_point_new (int *x,      ///< Point x-coordinate.
                   int *y,      ///< Point y-coordinate.
                   int *z,      ///< Point z-coordinate.
                   gsl_rng * rng)       ///< Pseudo-random number generator.
{
  *x = gsl_rng_uniform_int (rng, length);
  *y = gsl_rng_uniform_int (rng, width);
  *z = max_d;
#if DEBUG
  printf ("New point x %d y %d z %d\n", *x, *y, *z);
#endif
}

/**
 * Function to check the limits of a 2D tree point.
 */
static inline void
tree_3D_point_boundary (int *x, ///< Point x-coordinate.
                        int *y, ///< Point y-coordinate.
                        int *z, ///< Point z-coordinate.
                        gsl_rng * rng)  ///< Pseudo-random number generator.
{
  if (*z < 0 || *z == (int) height)
    {
      tree_3D_point_new (x, y, z, rng);
      return;
    }
  if (*x < 0)
    *x = length - 1;
  else if (*x == (int) length)
    *x = 0;
  if (*y < 0)
    *y = width - 1;
  else if (*y == (int) width)
    *y = 0;
#if DEBUG
  printf ("New point x %d y %d z %d\n", *x, *y, *z);
#endif
}

/**
 * Function to fix a 2D tree point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
tree_3D_point_fix (int x,       ///< Point x-coordinate.
                   int y,       ///< Point y-coordinate.
                   int z)       ///< point z-coordinate.
{
  register unsigned int *point;
  if (z > (int) max_d || z == 0 || y == 0 || x == 0 || z == (int) height - 1
      || y == (int) width - 1 || x == (int) length - 1)
    return 0;
  point = medium + z * area + y * length + x;
  if (point[1] || point[-1] || point[length] || point[-(int) length]
      || point[area] || point[-(int) area])
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      point[0] = 2;
      points_add (x, y, z, 2);
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to init a 3D tree.
 */
static inline void
tree_3D_init ()
{
  medium[length * (width / 2) + length / 2] = 2;
}

/**
 * Function to check the end a 3D tree.
 *
 * \return 1 on ending, 0 on continuing.
 */
static inline unsigned int
tree_3D_end (int z)             ///< Point z-coordinate.
{
  if (z == (int) max_d)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      ++max_d;
      g_mutex_unlock (mutex);
// END
    }
  if (max_d >= height - 1)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      max_d = height - 1;
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to check the limits of a 2D forest point.
 */
static inline void
forest_2D_point_boundary (int *x,       ///< Point x-coordinate.
                          int *y,       ///< Point y-coordinate.
                          gsl_rng * rng)
///< Pseudo-random number generator.
{
  if (*y == (int) height || *y < 0)
    {
      tree_2D_point_new (x, y, rng);
      return;
    }
  if (*x < 0)
    *x = width - 1;
  else if (*x == (int) width)
    *x = 0;
#if DEBUG
  printf ("Boundary point x %d y %d\n", *x, *y);
#endif
}

/**
 * Function to fix a 2D forest point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
forest_2D_point_fix (int x,     ///< Point x-coordinate.
                     int y,     ///< Point y-coordinate.
                     gsl_rng * rng)     ///< Pseudo-random number generator.
{
  register unsigned int k, *point;
  if (y > (int) max_d || x == 0 || x == (int) width - 1
      || y == (int) height - 1)
    return 0;
  point = medium + y * width + x;
  if (y == 0)
    {
      k = 1 + gsl_rng_uniform_int (rng, 15);
      goto forest;
    }
  k = point[1];
  if (k)
    goto forest;
  k = point[-1];
  if (k)
    goto forest;
  k = point[width];
  if (k)
    goto forest;
  k = point[-(int) width];
  if (k)
    goto forest;
  return 0;

forest:
// PARALLELIZING MUTEX
  g_mutex_lock (mutex);
  point[0] = k;
  points_add (x, y, 0, k);
  g_mutex_unlock (mutex);
// END
  return k;
}

/**
 * Function to check the limits of a 3D forest point.
 */
static inline void
forest_3D_point_boundary (int *x,       ///< Point x-coordinate.
                          int *y,       ///< Point y-coordinate.
                          int *z,       ///< Point z-coordinate.
                          gsl_rng * rng)
///< Pseudo-random number generator.
{
  if (*z == (int) height || *z < 0)
    {
      tree_3D_point_new (x, y, z, rng);
      return;
    }
  if (*y < 0)
    *y = width - 1;
  else if (*y == (int) width)
    *y = 0;
  if (*x < 0)
    *x = length - 1;
  else if (*x == (int) length)
    *x = 0;
#if DEBUG
  printf ("Boundary point x %d y %d z %d\n", *x, *y, *z);
#endif
}

/**
 * Function to fix a 3D forest point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
forest_3D_point_fix (int x,     ///< Point x-coordinate.
                     int y,     ///< Point y-coordinate.
                     int z,     ///< Point z-coordinate.
                     gsl_rng * rng)     ///< Pseudo-random number generator.
{
  register unsigned int k, *point;
  if (z > (int) max_d || y == 0 || x == 0 || z == (int) height - 1
      || y == (int) width - 1 || x == (int) length - 1)
    return 0;
  point = medium + z * area + y * length + x;
  if (z == 0)
    {
      k = 1 + gsl_rng_uniform_int (rng, 15);
      goto forest;
    }
  k = point[1];
  if (k)
    goto forest;
  k = point[-1];
  if (k)
    goto forest;
  k = point[length];
  if (k)
    goto forest;
  k = point[-(int) length];
  if (k)
    goto forest;
  k = point[area];
  if (k)
    goto forest;
  k = point[-(int) area];
  if (k)
    goto forest;
  return 0;

forest:
// PARALLELIZING MUTEX
  g_mutex_lock (mutex);
  point[0] = k;
  points_add (x, y, z, k);
  g_mutex_unlock (mutex);
// END
  return k;
}

/**
 * Function to start a new 2D neuron point.
 */
static inline void
neuron_2D_point_new (int *x,    ///< Point x-coordinate.
                     int *y,    ///< Point y-coordinate.
                     gsl_rng * rng)     ///< Pseudo-random number generator.
{
  register double angle;
  angle = 2 * M_PI * gsl_rng_uniform (rng);
  *x = width / 2 + max_d * cos (angle);
  *y = height / 2 + max_d * sin (angle);
#if DEBUG
  printf ("New point x %d y %d\n", *x, *y);
#endif
}

/**
 * Function to check the limits of a 2D neuron point.
 */
static inline void
neuron_2D_point_boundary (int *x,       ///< Point x-coordinate.
                          int *y,       ///< Point y-coordinate.
                          gsl_rng * rng)
///< Pseudo-random number generator.
{
  if (*y < 0 || *y == (int) height || *x < 0 || *x == (int) width)
    {
      neuron_2D_point_new (x, y, rng);
#if DEBUG
      printf ("Boundary point x %d y %d\n", *x, *y);
#endif
    }
}

/**
 * Function to fix a 2D neuron point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
neuron_2D_point_fix (int x,     ///< Point x-coordinate.
                     int y)     ///< Point y-coordinate.
{
  register unsigned int *point;
  if (x == 0 || y == 0 || x == (int) width - 1 || y == (int) height - 1)
    return 0;
  point = medium + y * width + x;
  if (point[1] || point[-1] || point[width] || point[-(int) width])
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      point[0] = 2;
      points_add (x, y, 0, 2);
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to init a 2D neuron.
 */
static inline void
neuron_2D_init ()
{
  medium[(height / 2) * width + width / 2] = 2;
}

/**
 * Function to check the end a 2D neuron.
 *
 * \return 1 on ending, 0 on continuing.
 */
static inline unsigned int
neuron_2D_end (int x,           ///< Point x-coordinate.
               int y)           ///< Point y-coordinate.
{
  register int r, k;
  r = 1 + round (sqrt (sqr (x - width / 2) + sqr (y - height / 2)));
  if (r >= (int) max_d)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      ++max_d;
      g_mutex_unlock (mutex);
// END
    }
  if (height < width)
    k = height;
  else
    k = width;
  k = k / 2 - 1;
  if ((int) max_d >= k)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      max_d = k;
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to start a new 3D neuron point.
 */
static inline void
neuron_3D_point_new (int *x,    ///< Point x-coordinate.
                     int *y,    ///< Point y-coordinate.
                     int *z,    ///< Point z-coordinate.
                     gsl_rng * rng)     ///< Pseudo-random number generator.
{
  double c1, s1, c2, s2;
  sincos (2. * M_PI * gsl_rng_uniform (rng), &s1, &c1);
  sincos (asin (2. * gsl_rng_uniform (rng) - 1), &s2, &c2);
  *x = length / 2 + max_d * c1 * c2;
  *y = width / 2 + max_d * s1 * c2;
  *z = height / 2 + max_d * s2;
#if DEBUG
  printf ("New point x %d y %d z %d\n", *x, *y, *z);
#endif
}

/**
 * Function to check the limits of a 3D neuron point.
 */
static inline void
neuron_3D_point_boundary (int *x,       ///< Point x-coordinate.
                          int *y,       ///< Point y-coordinate.
                          int *z,       ///< Point z-coordinate.
                          gsl_rng * rng)
///< Pseudo-random number generator.
{
  if (*z < 0 || *y < 0 || *x < 0 || *z == (int) height || *y == (int) width
      || *x == (int) length)
    {
      neuron_3D_point_new (x, y, z, rng);
#if DEBUG
      printf ("Boundary point x %d y %d z %d\n", *x, *y, *z);
#endif
    }
}

/**
 * Function to fix a 3D neuron point.
 *
 * \return 1 on fixing point, 0 on otherwise.
 */
static inline unsigned int
neuron_3D_point_fix (int x,     ///< Point x-coordinate.
                     int y,     ///< Point y-coordinate.
                     int z)     ///< Point z-coordinate.
{
  register unsigned int *point;
  if (z == 0 || y == 0 || x == 0 || z == (int) height - 1
      || y == (int) width - 1 || x == (int) length - 1)
    return 0;
  point = medium + z * area + y * length + x;
  if (point[1] || point[-1] || point[length] || point[-(int) length] ||
      point[area] || point[-(int) area])
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      point[0] = 2;
      points_add (x, y, z, 2);
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to init a 3D neuron.
 */
static inline void
neuron_3D_init ()
{
  medium[area * (height / 2) + length * (width / 2) + length / 2] = 2;
}

/**
 * Function to check the end a 3D neuron.
 *
 * \return 1 on ending, 0 on continuing.
 */
static inline unsigned int
neuron_3D_end (int x,           ///< Point x-coordinate.
               int y,           ///< Point y-coordinate.
               int z)           ///< Point z-coordinate.
{
  register int r, k;
  r = 1 + sqrt (sqr (x - length / 2) + sqr (y - width / 2)
                + sqr (z - height / 2));
  if (r >= (int) max_d)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      ++max_d;
      g_mutex_unlock (mutex);
// END
    }
  k = length;
  if ((int) width < k)
    k = width;
  if ((int) height < k)
    k = height;
  k = k / 2 - 1;
  if ((int) max_d >= k)
    {
// PARALLELIZING MUTEX
      g_mutex_lock (mutex);
      max_d = k;
      g_mutex_unlock (mutex);
// END
      return 1;
    }
  return 0;
}

/**
 * Function to stop the fractal simulation.
 */
void
fractal_stop ()
{
// PARALLELIZING MUTEX
  g_mutex_lock (mutex);
  breaking = 1;
  g_mutex_unlock (mutex);
// END
}

// PARALLELIZED FUNCTIONS

/**
 * Function to create a 2D fractal tree.
 *
 * \return NULL.
 */
void *
parallel_fractal_tree_2D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
#if DEBUG
  printf ("parallel_fractal_tree_2D: start\n");
#endif
  t0 = time (NULL);
  do
    {
#if DEBUG
      printf ("creating point\n");
#endif
      tree_2D_point_new (&x, &y, rng);
#if DEBUG
      printf ("checking fix\n");
#endif
      while (!breaking && !tree_2D_point_fix (x, y))
        {
#if DEBUG
          printf ("moving point\n");
#endif
          point_2D_move (&x, &y, rng);
#if DEBUG
          printf ("checking boundary\n");
#endif
          tree_2D_point_boundary (&x, &y, rng);
        }
#if DEBUG
      printf ("checking end\n");
#endif
      if (animating && time (NULL) > t0)
        break;
      if (tree_2D_end (y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal tree.
 *
 * \return NULL.
 */
void *
parallel_fractal_tree_3D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !tree_3D_point_fix (x, y, z))
        {
          point_3D_move (&x, &y, &z, rng);
          tree_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_3D_end (z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 2D fractal forest.
 *
 * \return NULL.
 */
void *
parallel_fractal_forest_2D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_2D_point_new (&x, &y, rng);
      while (!breaking && !forest_2D_point_fix (x, y, rng))
        {
          point_2D_move (&x, &y, rng);
          forest_2D_point_boundary (&x, &y, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_2D_end (y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal forest.
 *
 * \return NULL.
 */
void *
parallel_fractal_forest_3D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !forest_3D_point_fix (x, y, z, rng))
        {
          point_3D_move (&x, &y, &z, rng);
          forest_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_3D_end (z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 2D fractal neuron.
 *
 * \return NULL.
 */
void *
parallel_fractal_neuron_2D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
  t0 = time (NULL);
  do
    {
      neuron_2D_point_new (&x, &y, rng);
      while (!breaking && !neuron_2D_point_fix (x, y))
        {
          point_2D_move (&x, &y, rng);
          neuron_2D_point_boundary (&x, &y, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (neuron_2D_end (x, y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal neuron.
 *
 * \return NULL.
 */
void *
parallel_fractal_neuron_3D (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      neuron_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !neuron_3D_point_fix (x, y, z))
        {
          point_3D_move (&x, &y, &z, rng);
          neuron_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (neuron_3D_end (x, y, z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 2D fractal tree with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_tree_2D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_2D_point_new (&x, &y, rng);
      while (!breaking && !tree_2D_point_fix (x, y))
        {
          point_2D_move (&x, &y, rng);
          tree_2D_point_boundary (&x, &y, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_2D_end (y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal tree with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_tree_3D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !tree_3D_point_fix (x, y, z))
        {
          point_3D_move (&x, &y, &z, rng);
          tree_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_3D_end (z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 2D fractal forest with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_forest_2D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_2D_point_new (&x, &y, rng);
      while (!breaking && !forest_2D_point_fix (x, y, rng))
        {
          point_2D_move_diagonal (&x, &y, rng);
          forest_2D_point_boundary (&x, &y, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_2D_end (y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal forest with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_forest_3D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      tree_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !forest_3D_point_fix (x, y, z, rng))
        {
          point_3D_move_diagonal (&x, &y, &z, rng);
          forest_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (tree_3D_end (z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 2D fractal neuron with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_neuron_2D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y;
  long t0;
  t0 = time (NULL);
  do
    {
      neuron_2D_point_new (&x, &y, rng);
      while (!breaking && !neuron_2D_point_fix (x, y))
        {
          point_2D_move_diagonal (&x, &y, rng);
          neuron_2D_point_boundary (&x, &y, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (neuron_2D_end (x, y))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

/**
 * Function to create a 3D fractal neuron with diagonal movements.
 *
 * \return NULL.
 */
void *
parallel_fractal_neuron_3D_diagonal (gsl_rng * rng)
///< Pseudo-random number generator.
{
  int x, y, z;
  long t0;
  t0 = time (NULL);
  do
    {
      neuron_3D_point_new (&x, &y, &z, rng);
      while (!breaking && !neuron_3D_point_fix (x, y, z))
        {
          point_3D_move_diagonal (&x, &y, &z, rng);
          neuron_3D_point_boundary (&x, &y, &z, rng);
        }
      if (animating && time (NULL) > t0)
        break;
      if (neuron_3D_end (x, y, z))
        fractal_stop ();
    }
  while (!breaking);
  g_thread_exit (NULL);
  return NULL;
}

//END OF PARALLELIZED FUNCTIONS

/**
 * Function to start the fractal functions and data.
 */
void
medium_start ()
{
  register int i, j;

#if DEBUG
  printf ("Deleting points\n");
#endif
  g_free (point);
  point = NULL;
  npoints = 0;

  j = area = width * length;
  if (fractal_3D)
    j *= length;
  medium_bytes = j * sizeof (unsigned int);
  medium = (unsigned int *) g_slice_alloc (medium_bytes);
  for (i = j; --i >= 0;)
    medium[i] = 0;
#if DEBUG
  printf ("Medium size=%d pointer=%ld\n", j, (size_t) medium);
#endif

#if DEBUG
  printf ("Setting functions\n");
#endif
  if (fractal_diagonal)
    {
      if (fractal_3D)
        {
          switch (fractal_type)
            {
            case FRACTAL_TYPE_TREE:
              tree_3D_init ();
              parallel_fractal = parallel_fractal_tree_3D_diagonal;
              max_d = 1;
              break;
            case FRACTAL_TYPE_FOREST:
              parallel_fractal = parallel_fractal_forest_3D_diagonal;
              max_d = 1;
              break;
            default:
              neuron_3D_init ();
              parallel_fractal = parallel_fractal_neuron_3D_diagonal;
              max_d = 2;
            }
        }
      else
        {
          switch (fractal_type)
            {
            case FRACTAL_TYPE_TREE:
              tree_2D_init ();
              parallel_fractal = parallel_fractal_tree_2D_diagonal;
              max_d = 1;
              break;
            case FRACTAL_TYPE_FOREST:
              parallel_fractal = parallel_fractal_forest_2D_diagonal;
              max_d = 1;
              break;
            default:
              neuron_2D_init ();
              parallel_fractal = parallel_fractal_neuron_2D_diagonal;
              max_d = 2;
            }
        }
    }
  else
    {
      if (fractal_3D)
        {
          switch (fractal_type)
            {
            case FRACTAL_TYPE_TREE:
              tree_3D_init ();
              parallel_fractal = parallel_fractal_tree_3D;
              max_d = 1;
              break;
            case FRACTAL_TYPE_FOREST:
              parallel_fractal = parallel_fractal_forest_3D;
              max_d = 1;
              break;
            default:
              neuron_3D_init ();
              parallel_fractal = parallel_fractal_neuron_3D;
              max_d = 2;
            }
        }
      else
        {
          switch (fractal_type)
            {
            case FRACTAL_TYPE_TREE:
              tree_2D_init ();
              parallel_fractal = parallel_fractal_tree_2D;
              max_d = 1;
              break;
            case FRACTAL_TYPE_FOREST:
              parallel_fractal = parallel_fractal_forest_2D;
              max_d = 1;
              break;
            default:
              neuron_2D_init ();
              parallel_fractal = parallel_fractal_neuron_2D;
              max_d = 2;
            }
        }
    }
}

/**
 * Function to get an unsigned integer number of a XML node property.
 *
 * \return Unsigned integer number value.
 */
static unsigned int
xml_node_get_uint_with_default (xmlNode * node, ///< XML node.
                                const xmlChar * prop,   ///< XML property.
                                unsigned int default_value,     ///< Default value.
                                int *error_code)        ///< Error code.
{
  unsigned int i = 0;
  xmlChar *buffer;
  buffer = xmlGetProp (node, prop);
  *error_code = 0;
  if (!buffer)
    return default_value;
  else
    {
      if (sscanf ((char *) buffer, "%u", &i) != 1)
        *error_code = 1;
      xmlFree (buffer);
    }
  return i;
}

/**
 * Function to get an unsigned long integer number of a XML node property.
 *
 * \return Unsigned long integer number value.
 */
static unsigned long
xml_node_get_ulong_with_default (xmlNode * node,        ///< XML node.
                                 const xmlChar * prop,  ///< XML property.
                                 unsigned long default_value,
                                 ///< Default value.
                                 int *error_code)       ///< Error code.
{
  unsigned long i = 0l;
  xmlChar *buffer;
  buffer = xmlGetProp (node, prop);
  *error_code = 0;
  if (!buffer)
    return default_value;
  else
    {
      if (sscanf ((char *) buffer, "%lu", &i) != 1)
        *error_code = 1;
      xmlFree (buffer);
    }
  return i;
}

/**
 * Function to open the fractal data on a file.
 *
 * \return 1 on success, 0 on error.
 */
int
fractal_input (char *filename)  ///< File name.
{
  xmlDoc *doc;
  xmlNode *node;
  xmlChar *buffer;
  const char *error_message;
  int error_code;
  buffer = NULL;
  xmlKeepBlanksDefault (0);
  doc = xmlParseFile ((const char *) filename);
  if (!doc)
    {
      error_message = _("Unable to parse the input file");
      goto exit_on_error;
    }
  node = xmlDocGetRootElement (doc);
  if (!node || xmlStrcmp (node->name, XML_FRACTAL))
    {
      error_message = _("Bad XML root node");
      goto exit_on_error;
    }
  width = xml_node_get_uint_with_default (node, XML_WIDTH, WIDTH, &error_code);
  if (error_code)
    {
      error_message = _("Bad width");
      goto exit_on_error;
    }
  height
    = xml_node_get_uint_with_default (node, XML_HEIGHT, HEIGHT, &error_code);
  if (error_code)
    {
      error_message = _("Bad height");
      goto exit_on_error;
    }
  length
    = xml_node_get_uint_with_default (node, XML_LENGTH, LENGTH, &error_code);
  if (error_code)
    {
      error_message = _("Bad length");
      goto exit_on_error;
    }
  random_seed
    = xml_node_get_ulong_with_default (node, XML_SEED, SEED, &error_code);
  if (error_code)
    {
      error_message = _("Bad random seed");
      goto exit_on_error;
    }
  nthreads
    = xml_node_get_uint_with_default (node, XML_THREADS, threads_number (),
                                      &error_code);
  if (!nthreads || error_code)
    {
      error_message = _("Bad threads number");
      goto exit_on_error;
    }
  buffer = xmlGetProp (node, XML_DIAGONAL);
  if (!buffer || !xmlStrcmp (buffer, XML_NO))
    fractal_diagonal = 0;
  else if (!xmlStrcmp (buffer, XML_YES))
    fractal_diagonal = 1;
  else
    {
      error_message = _("Bad diagonal movement");
      goto exit_on_error;
    }
  xmlFree (buffer);
  buffer = xmlGetProp (node, XML_3D);
  if (!buffer || !xmlStrcmp (buffer, XML_NO))
    fractal_3D = 0;
  else if (!xmlStrcmp (buffer, XML_YES))
    fractal_3D = 1;
  else
    {
      error_message = _("Bad 3D");
      goto exit_on_error;
    }
  xmlFree (buffer);
  buffer = xmlGetProp (node, XML_ANIMATE);
  if (!buffer || !xmlStrcmp (buffer, XML_YES))
    animating = 1;
  else if (!xmlStrcmp (buffer, XML_NO))
    animating = 0;
  else
    {
      error_message = _("Bad animation");
      goto exit_on_error;
    }
  xmlFree (buffer);
  buffer = xmlGetProp (node, XML_TYPE);
  if (!buffer || !xmlStrcmp (buffer, XML_TREE))
    fractal_type = FRACTAL_TYPE_TREE;
  else if (!xmlStrcmp (buffer, XML_FOREST))
    fractal_type = FRACTAL_TYPE_FOREST;
  else if (!xmlStrcmp (buffer, XML_NEURON))
    fractal_type = FRACTAL_TYPE_NEURON;
  else
    {
      error_message = _("Unknown fractal type");
      goto exit_on_error;
    }
  xmlFree (buffer);
  buffer = xmlGetProp (node, XML_RANDOM_SEED);
  if (!buffer || !xmlStrcmp (buffer, XML_CLOCK))
    random_seed_type = RANDOM_SEED_TYPE_CLOCK;
  else if (!xmlStrcmp (buffer, XML_DEFAULT))
    random_seed_type = RANDOM_SEED_TYPE_DEFAULT;
  else if (!xmlStrcmp (buffer, XML_FIXED))
    random_seed_type = RANDOM_SEED_TYPE_FIXED;
  else
    {
      error_message = _("Unknown random seed type");
      goto exit_on_error;
    }
  xmlFree (buffer);
  buffer = xmlGetProp (node, XML_RANDOM_TYPE);
  if (!buffer || !xmlStrcmp (buffer, XML_MT19937))
    random_algorithm = 0;
  else if (!xmlStrcmp (buffer, XML_RANLXS0))
    random_algorithm = 1;
  else if (!xmlStrcmp (buffer, XML_RANLXS1))
    random_algorithm = 2;
  else if (!xmlStrcmp (buffer, XML_RANLXS2))
    random_algorithm = 3;
  else if (!xmlStrcmp (buffer, XML_RANLXD1))
    random_algorithm = 4;
  else if (!xmlStrcmp (buffer, XML_RANLXD2))
    random_algorithm = 5;
  else if (!xmlStrcmp (buffer, XML_RANLUX))
    random_algorithm = 6;
  else if (!xmlStrcmp (buffer, XML_RANLUX389))
    random_algorithm = 7;
  else if (!xmlStrcmp (buffer, XML_CMRG))
    random_algorithm = 8;
  else if (!xmlStrcmp (buffer, XML_MRG))
    random_algorithm = 9;
  else if (!xmlStrcmp (buffer, XML_TAUS2))
    random_algorithm = 10;
  else if (!xmlStrcmp (buffer, XML_GFSR4))
    random_algorithm = 11;
  else
    {
      error_message = _("Unknown random algorithm");
      goto exit_on_error;
    }
  xmlFree (buffer);
  xmlFreeDoc (doc);
  return 1;

exit_on_error:
  show_error (error_message);
  xmlFree (buffer);
  xmlFreeDoc (doc);
  return 0;
}

/**
 * Function with the main bucle to draw the fractal.
 */
void
fractal ()
{
  const gsl_rng_type *random_type[N_RANDOM_TYPES] = {
    gsl_rng_mt19937,
    gsl_rng_ranlxs0,
    gsl_rng_ranlxs1,
    gsl_rng_ranlxs2,
    gsl_rng_ranlxd1,
    gsl_rng_ranlxd2,
    gsl_rng_ranlux,
    gsl_rng_ranlux389,
    gsl_rng_cmrg,
    gsl_rng_mrg,
    gsl_rng_taus2,
    gsl_rng_gfsr4
  };
  FILE *file;
  unsigned int i;

// PARALLELIZING DATA
  gsl_rng *rng[nthreads];
  GThread *thread[nthreads];

  t0 = time (NULL);
#if DEBUG
  printf ("t0=%lu\n", t0);
#endif

#if DEBUG
  printf ("Opening log file\n");
#endif
  file = fopen ("log", "w");

#if DEBUG
  printf ("Opening pseudo-random generators\n");
#endif
  for (i = 0; i < nthreads; ++i)
    {
      rng[i] = gsl_rng_alloc (random_type[random_algorithm]);
      switch (random_seed_type)
        {
        case RANDOM_SEED_TYPE_DEFAULT:
          break;
        case RANDOM_SEED_TYPE_CLOCK:
          gsl_rng_set (rng[i], (unsigned long) clock () + i);
          break;
        default:
          gsl_rng_set (rng[i], random_seed + i);
        }
    }

// END

  breaking = 0;
  simulating = 1;

#if DEBUG
  printf ("Updating simulator dialog\n");
#endif
  dialog_simulator_update ();

#if DEBUG
  printf ("Starting medium\n");
#endif
  medium_start ();

#if DEBUG
  printf ("Main bucle\n");
#endif
  do
    {
#if DEBUG
      printf ("Calculating fractal\n");
#endif
// PARALLELIZING CALLS
      for (i = 0; i < nthreads; ++i)
        thread[i] = g_thread_new (NULL, (void (*)) parallel_fractal, rng[i]);
      for (i = 0; i < nthreads; ++i)
        g_thread_join (thread[i]);
// END

#if DEBUG
      printf ("Updating simulator dialog\n");
#endif
      dialog_simulator_progress ();

// DISPLAYING DRAW
#if DEBUG
      printf ("Redisplaying draw\n");
#endif
      draw ();
// END

#if DEBUG
      printf ("Saving log data\n");
#endif
      fprintf (file, "%d %d\n", max_d, npoints);
    }
  while (!breaking);

#if DEBUG
  printf ("Closing log file\n");
#endif
  fclose (file);

#if DEBUG
  printf ("Updating simulator dialog\n");
#endif
  breaking = simulating = 0;
  dialog_simulator_update ();

#if DEBUG
  printf ("Freeing threads\n");
#endif
  for (i = 0; i < nthreads; ++i)
    gsl_rng_free (rng[i]);
  g_slice_free1 (medium_bytes, medium);
}

max_stars_repo_name: jburguete/fractal

max_issues_repo_name: jburguete/fractal

max_forks_repo_name: jburguete/fractal


text:
#ifndef bcfunctions_h
#define bcfunctions_h

#include <petsc.h>

PetscErrorCode BCsDiff(PetscInt dim, PetscReal time, const PetscReal x[],
                       PetscInt num_comp_u, PetscScalar *u, void *ctx);
PetscErrorCode BCsMass(PetscInt dim, PetscReal time, const PetscReal x[],
                       PetscInt num_comp_u, PetscScalar *u, void *ctx);

#endif // bcfunctions_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
//
//  trsm.h
//  Linear Algebra Template Library
//
//  Created by Rodney James on 1/5/12.
//  Copyright (c) 2012 University of Colorado Denver. All rights reserved.
//

#ifndef _trsm_h
#define _trsm_h

/// @file trsm.h Solves triangular matrix equations.

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

namespace LATL
{
   /// @brief Solves a real triangular matrix equation.
   ///
   /// For a real upper or lower triangular matrix A, real rectangular matrix B and real scalar alpha,
   ///
   ///        A*X=alpha*B  or  A'*X=alpha*B  or  X*A=alpha*B  or  X*A'=alpha*B
   /// is solved for the matrix X.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then A*X=alpha*B  or  A'*X=alpha*B is solved
   ///        if side = 'R' or 'r' then X*A=alpha*B  or  X*A'=alpha*B is solved
   /// @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 trans Specifies whether the transpose or conjugate transpose of A is to be used:
   ///
   ///        if trans = 'N' or 'n' then A*X=alpha*B or X*A=alpha*B is solved
   ///        if trans = 'T' or 't' then A'*X=alpha*B or X*A'=alpha*B is solved
   ///        if trans = 'C' or 'c' then A'*X=alpha*B or X*A'=alpha*B is solved
   /// @param diag specifies whether or not A is unit triangular as follows:
   ///
   ///        if diag = 'U' or 'u' then A is assumed to be unit triangular
   ///        if diag = 'N' or 'n' then A is not assumed to be unit triangular.
   /// If A is unit triangular, the diagonal elements are assumed to be unity and are not referenced.
   /// @param m Specifies the number of rows of the matrix B.  m>=0
   /// @param n Specifies the number of columns of the matrix B.  n>=0
   /// @param alpha Real scalar.
   /// @param A Pointer to real triangular matrix A.  The order of A is n if side = 'L' or 'l';
   /// the order of A is m is side = 'R' or 'r'.
   /// If uplo = 'U' or 'u, A is upper triangular and the lower triangular part is not referenced.
   /// If uplo = 'L' or 'l A is lower triangular and the upper triangular part is not referenced.
   /// @param ldA Column length of the matrix A.  If side='L' or 'l' then ldA>=n; if side='R' or 'r' then ldA>=m.
   /// @param B Pointer to real m-by-n matrix B.  On entry, B contains the right hand side matrix B.
   /// On exit, B is overwritten with the solution matrix X.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int TRSM(char side, char uplo, char trans, char diag, int_t m, int_t n, real_t alpha, real_t *A, int_t ldA, real_t *B, int_t ldB)
   {
      using std::toupper;

      const real_t zero(0.0);
      const real_t one(1.0);
      int_t i,j,k;
      real_t *a,*b,*bt;
      real_t t;
      bool nounit;

      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      nounit=(diag=='N')?1:0;

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      if(alpha==zero)
      {
         b=B;
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               b[i]=zero;
            b+=ldB;
         }
      }
      else
      {
         if(side=='L')
         {
            if(trans=='N')
            {
               if(uplo=='U')
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a=A+m*ldA;
                     for(k=m-1;k>=0;k--)
                     {
                        a-=ldA;
                        if(nounit)
                           b[k]=b[k]/a[k];
                        for(i=0;i<k;i++)
                           b[i]-=b[k]*a[i];
                     }
                     b+=ldB;
                  }
               }
               else
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a=A;
                     for(k=0;k<m;k++)
                     {
                        if(nounit)
                           b[k]=b[k]/a[k];
                        for(i=k+1;i<m;i++)
                           b[i]-=b[k]*a[i];
                        a+=ldA;
                     }
                     b+=ldB;
                  }
               }
            }
            else
            {
               if(uplo=='U')
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A;
                     for(i=0;i<m;i++)
                     {
                        t=alpha*b[i];
                        for(k=0;k<i;k++)
                           t-=a[k]*b[k];
                        if(nounit)
                           t=t/a[i];
                        b[i]=t;
                        a+=ldA;
                     }
                     b+=ldB;
                  }
               }
               else
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A+m*ldA;
                     for(i=m-1;i>=0;i--)
                     {
                        a-=ldA;
                        t=alpha*b[i];
                        for(k=i+1;k<m;k++)
                           t-=a[k]*b[k];
                        if(nounit)
                           t=t/a[i];
                        b[i]=t;
                     }
                     b+=ldB;
                  }
               }
            }
         }
         else
         {
            if(trans=='N')
            {
               if(uplo=='U')
               {
                  b=B;
                  a=A;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     bt=B;
                     for(k=0;k<j;k++)
                     {
                        for(i=0;i<m;i++)
                           b[i]-=a[k]*bt[i];
                        bt+=ldB;
                     }
                     if(nounit)
                     {
                        t=one/a[j];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     a+=ldA;
                     b+=ldB;
                  }
               }
               else
               {
                  b=B+n*ldB;
                  a=A+n*ldA;
                  for(j=n-1;j>=0;j--)
                  {
                     a-=ldA;
                     b-=ldB;
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     bt=B+(j+1)*ldB;
                     for(k=j+1;k<n;k++)
                     {
                        for(i=0;i<m;i++)
                           b[i]-=a[k]*bt[i];
                        bt+=ldB;
                     }
                     if(nounit)
                     {
                        t=one/a[j];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                  }
               }
            }
            else
            {
               if(uplo=='U')
               {
                  b=B+n*ldB;
                  a=A+n*ldA;
                  for(k=n-1;k>=0;k--)
                  {
                     a-=ldA;
                     b-=ldB;
                     if(nounit)
                     {
                        t=one/a[k];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B;
                     for(j=0;j<k;j++)
                     {
                        t=a[j];
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                  }
               }
               else
               {
                  a=A;
                  b=B;
                  for(k=0;k<n;k++)
                  {
                     if(nounit)
                     {
                        t=one/a[k];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B+(k+1)*ldB;
                     for(j=k+1;j<n;j++)
                     {
                        t=a[j];
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a+=ldA;
                     b+=ldB;
                  }
               }
            }
         }
      }
      return 0;
   }

   /// @brief Solves a complex triangular matrix equation.
   ///
   /// For a complex upper or lower triangular matrix A, complex rectangular matrix B and complex scalar alpha,
   ///
   ///        A*X=alpha*B or A'*X=alpha*B or A.'*X=alpha*B or X*A=alpha*B or X*A'=alpha*B or X*A.'=alpha*B
   /// is solved for the matrix X.
   /// @return 0 if success.
   /// @return -i if the ith argument is invalid.
   /// @tparam real_t Floating point type.
   /// @param side Specifies whether the matrix A appears on the left or right side as follows:
   ///
   ///        if side = 'L' or 'l' then A*X=alpha*B  or  A'*X=alpha*B is solved
   ///        if side = 'R' or 'r' then X*A=alpha*B  or  X*A'=alpha*B is solved
   /// @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 trans Specifies whether the transpose of A is to be used or not:
   ///
   ///        if trans = 'N' or 'n' then A*X=alpha*B or X*A=alpha*B is solved
   ///        if trans = 'T' or 't' then A.'*X=alpha*B or X*A.'=alpha*B is solved
   ///        if trans = 'C' or 'c' then A'*X=alpha*B or X*A'=alpha*B is solved
   /// @param diag specifies whether or not A is unit triangular as follows:
   ///
   ///        if diag = 'U' or 'u' then A is assumed to be unit triangular
   ///        if diag = 'N' or 'n' then A is not assumed to be unit triangular.
   /// If A is unit triangular, the diagonal elements are assumed to be unity and are not referenced.
   /// @param m Specifies the number of rows of the matrix B.  m>=0
   /// @param n Specifies the number of columns of the matrix B.  n>=0
   /// @param alpha Complex scalar.
   /// @param A Pointer to complex triangular matrix A.  The order of A is n if side = 'L' or 'l';
   /// the order of A is m is side = 'R' or 'r'.
   /// If uplo = 'U' or 'u, A is upper triangular and the lower triangular part is not referenced.
   /// If uplo = 'L' or 'l A is lower triangular and the upper triangular part is not referenced.
   /// @param ldA Column length of the matrix A.  If side='L' or 'l' then ldA>=n; if side='R' or 'r' then ldA>=m.
   /// @param B Pointer to complex m-by-n matrix B.  On entry, B contains the right hand side matrix B.
   /// On exit, B is overwritten with the solution matrix X.
   /// @param ldB Column length of the matrix B.  ldB>=m.
   /// @ingroup BLAS

   template <typename real_t>
   int TRSM(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<real_t> alpha, complex<real_t> *A, int_t ldA, complex<real_t> *B, int_t ldB)
   {

      using std::conj;
      using std::toupper;
      const complex<real_t> zero(0.0,0.0);
      const complex<real_t> one(1.0,0.0);
      int_t i,j,k;
      complex<real_t> *a,*b,*bt;
      complex<real_t> t;
      bool nounit;

      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);

      nounit=(diag=='N')?1:0;

      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      if(alpha==zero)
      {
         b=B;
         for(j=0;j<n;j++)
         {
            for(i=0;i<m;i++)
               b[i]=zero;
            b+=ldB;
         }
      }
      else
      {
         if(side=='L')
         {
            if(trans=='N')
            {
               if(uplo=='U')
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a=A+m*ldA;
                     for(k=m-1;k>=0;k--)
                     {
                        a-=ldA;
                        if(nounit)
                           b[k]=b[k]/a[k];
                        for(i=0;i<k;i++)
                           b[i]-=b[k]*a[i];
                     }
                     b+=ldB;
                  }
               }
               else
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a=A;
                     for(k=0;k<m;k++)
                     {
                        if(nounit)
                           b[k]=b[k]/a[k];
                        for(i=k+1;i<m;i++)
                           b[i]-=b[k]*a[i];
                        a+=ldA;
                     }
                     b+=ldB;
                  }
               }
            }
            else if(trans=='T')
            {
               if(uplo=='U')
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A;
                     for(i=0;i<m;i++)
                     {
                        t=alpha*b[i];
                        for(k=0;k<i;k++)
                           t-=a[k]*b[k];
                        if(nounit)
                           t=t/a[i];
                        b[i]=t;
                        a+=ldA;
                     }
                     b+=ldB;
                  }
               }
               else
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A+m*ldA;
                     for(i=m-1;i>=0;i--)
                     {
                        a-=ldA;
                        t=alpha*b[i];
                        for(k=i+1;k<m;k++)
                           t-=a[k]*b[k];
                        if(nounit)
                           t=t/a[i];
                        b[i]=t;
                     }
                     b+=ldB;
                  }
               }
            }
            else if(trans=='C')
            {
               if(uplo=='U')
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A;
                     for(i=0;i<m;i++)
                     {
                        t=alpha*b[i];
                        for(k=0;k<i;k++)
                           t-=conj(a[k])*b[k];
                        if(nounit)
                           t=t/conj(a[i]);
                        b[i]=t;
                        a+=ldA;
                     }
                     b+=ldB;
                  }
               }
               else
               {
                  b=B;
                  for(j=0;j<n;j++)
                  {
                     a=A+m*ldA;
                     for(i=m-1;i>=0;i--)
                     {
                        a-=ldA;
                        t=alpha*b[i];
                        for(k=i+1;k<m;k++)
                           t-=conj(a[k])*b[k];
                        if(nounit)
                           t=t/conj(a[i]);
                        b[i]=t;
                     }
                     b+=ldB;
                  }
               }
            }
         }
         else
         {
            if(trans=='N')
            {
               if(uplo=='U')
               {
                  b=B;
                  a=A;
                  for(j=0;j<n;j++)
                  {
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     bt=B;
                     for(k=0;k<j;k++)
                     {
                        for(i=0;i<m;i++)
                           b[i]-=a[k]*bt[i];
                        bt+=ldB;
                     }
                     if(nounit)
                     {
                        t=one/a[j];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     a+=ldA;
                     b+=ldB;
                  }
               }
               else
               {
                  b=B+n*ldB;
                  a=A+n*ldA;
                  for(j=n-1;j>=0;j--)
                  {
                     a-=ldA;
                     b-=ldB;
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     bt=B+(j+1)*ldB;
                     for(k=j+1;k<n;k++)
                     {
                        for(i=0;i<m;i++)
                           b[i]-=a[k]*bt[i];
                        bt+=ldB;
                     }
                     if(nounit)
                     {
                        t=one/a[j];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                  }
               }
            }
            else if(trans=='T')
            {
               if(uplo=='U')
               {
                  b=B+n*ldB;
                  a=A+n*ldA;
                  for(k=n-1;k>=0;k--)
                  {
                     a-=ldA;
                     b-=ldB;
                     if(nounit)
                     {
                        t=one/a[k];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B;
                     for(j=0;j<k;j++)
                     {
                        t=a[j];
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                  }
               }
               else
               {
                  a=A;
                  b=B;
                  for(k=0;k<n;k++)
                  {
                     if(nounit)
                     {
                        t=one/a[k];
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B+(k+1)*ldB;
                     for(j=k+1;j<n;j++)
                     {
                        t=a[j];
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a+=ldA;
                     b+=ldB;
                  }
               }
            }
            else if(trans=='C')
            {
               if(uplo=='U')
               {
                  b=B+n*ldB;
                  a=A+n*ldA;
                  for(k=n-1;k>=0;k--)
                  {
                     a-=ldA;
                     b-=ldB;
                     if(nounit)
                     {
                        t=one/conj(a[k]);
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B;
                     for(j=0;j<k;j++)
                     {
                        t=conj(a[j]);
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                  }
               }
               else
               {
                  a=A;
                  b=B;
                  for(k=0;k<n;k++)
                  {
                     if(nounit)
                     {
                        t=one/conj(a[k]);
                        for(i=0;i<m;i++)
                           b[i]*=t;
                     }
                     bt=B+(k+1)*ldB;
                     for(j=k+1;j<n;j++)
                     {
                        t=conj(a[j]);
                        for(i=0;i<m;i++)
                           bt[i]-=t*b[i];
                        bt+=ldB;
                     }
                     for(i=0;i<m;i++)
                        b[i]*=alpha;
                     a+=ldA;
                     b+=ldB;
                  }
               }
            }
         }
      }
      return 0;
   }

#ifdef __latl_cblas
#include <cblas.h>

   template <> int TRSM<float>(char side, char uplo, char trans, char diag, int_t m, int_t n, float alpha, float *A, int_t ldA, float *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;

      cblas_strsm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,alpha,A,ldA,B,ldB);

      return 0;
   }

   template <> int TRSM<double>(char side, char uplo, char trans, char diag, int_t m, int_t n, double alpha, double *A, int_t ldA, double *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;

      cblas_dtrsm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,alpha,A,ldA,B,ldB);

      return 0;
   }
   template <> int TRSM<float>(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<float> alpha, complex<float> *A, int_t ldA, complex<float> *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;

      cblas_ctrsm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,&alpha,A,ldA,B,ldB);

      return 0;
   }

   template <> int TRSM<double>(char side, char uplo, char trans, char diag, int_t m, int_t n, complex<double> alpha, complex<double> *A, int_t ldA, complex<double> *B, int_t ldB)
   {
      using std::toupper;
      side=toupper(side);
      uplo=toupper(uplo);
      trans=toupper(trans);
      diag=toupper(diag);
      if((side!='L')&&(side!='R'))
         return -1;
      else if((uplo!='U')&&(uplo!='L'))
         return -2;
      else if((trans!='N')&&(trans!='T')&&(trans!='C'))
         return -3;
      else if((diag!='U')&&(diag!='N'))
         return -4;
      else if(m<0)
         return -5;
      else if(n<0)
         return -6;
      else if(ldA<((side=='L')?m:n))
         return -9;
      else if(ldB<m)
         return -11;
      else if((m==0)||(n==0))
         return 0;

      const CBLAS_SIDE Side=(side=='L')?CblasLeft:CblasRight;
      const CBLAS_UPLO Uplo=(uplo=='U')?CblasUpper:CblasLower;
      const CBLAS_TRANSPOSE Trans=(trans=='N')?CblasNoTrans:((trans=='T')?CblasTrans:CblasConjTrans);
      const CBLAS_DIAG Diag=(diag=='N')?CblasNonUnit:CblasUnit;

      cblas_ztrsm(CblasColMajor,Side,Uplo,Trans,Diag,m,n,&alpha,A,ldA,B,ldB);

      return 0;
   }
#endif
}
#endif

max_stars_repo_name: langou/latl

max_issues_repo_name: langou/latl

max_forks_repo_name: langou/latl