text:
/*! \file rho.h
    \brief 電子密度ρ(r)を求めるクラスの宣言

    Copyright © 2015 @dc1394 All Rights Reserved.
    This software is released under the BSD 2-Clause License.
*/

#include "data.h"
#include "property.h"
#include <memory>           // for std::unique_ptr
#include <vector>           // for std::vector
#include <gsl/gsl_spline.h> // for gsl_interp_accel, gsl_interp_accel_free, gsl_spline, gsl_spline_free

namespace schrac {
    //! A class.
    /*!
        電子密度ρ(r)を求めるクラス
    */
    class Rho final {
        // #region コンストラクタ・デストラクタ

    public:
        //! A constructor.
        /*!
            唯一のコンストラクタ
            \param pdiffdata 微分方程式のデータオブジェクト
        */
        Rho(std::shared_ptr<DiffData> const & pdiffdata);

        //! A destructor.
        /*!
            デフォルトデストラクタ
        */
        ~Rho() = default;

        // #region メンバ関数

        //!  A public member function.
        /*!
            現在の電子密度ρ(r)で初期化する
        */
        void init();

        //!  A public member function (const).
        /*!
            電子密度ρ(r)を返す
            \param r rの値
            \return ρ(r)の値
        */
        double operator()(double r) const;

        //!  A public member function (const).
        /*!
            新しい電子密度ρnew(r)と、電子密度ρ(r)を一次混合する
            \param rhonew 新しい電子密度ρnew(r)
        */
        void rhomix(dvector const & rhonew);

        // #endregion メンバ関数

        // #region プロパティ

    public:
        //! A property.
        /*!
            電子密度が格納された可変長配列へのプロパティ
        */
        Property< std::vector<double> > const PRho;

        // #endregion プロパティ

        // #region メンバ変数

    private:
        //! A private member variable.
        /*!
            gsl_interp_accelへのスマートポインタ
        */
        std::unique_ptr<gsl_interp_accel, decltype(&gsl_interp_accel_free)> const acc_;

        //!  A private member variable.
        /*!
            データオブジェクト
        */
        std::shared_ptr<DiffData> pdiffdata_;

        //!  A private member variable.
        /*!
            密度ρ(r)
        */
        std::vector<double> rho_;

        //! A private member variable.
        /*!
            gsl_splineへのスマートポインタ
        */
        std::unique_ptr<gsl_spline, decltype(&gsl_spline_free)> const spline_;

        // #endregion メンバ変数

        // #region 禁止されたコンストラクタ・メンバ関数

        //! A private constructor (deleted).
        /*!
            デフォルトコンストラクタ（禁止）
        */
        Rho() = delete;

        //! A private copy constructor (deleted).
        /*!
            コピーコンストラクタ（禁止）
        */
        Rho(Rho const &) = delete;

        //! A private member function (deleted).
        /*!
            operator=()の宣言（禁止）
            \param コピー元のオブジェクト（未使用）
            \return コピー元のオブジェクト
        */
        Rho & operator=(Rho const &) = delete;

        // #endregion 禁止されたコンストラクタ・メンバ関数
    };
}


text:
#ifndef _LOGIT_H
#define _LOGIT_H

#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

class logit
{
  private:
    gsl_vector *y; // response vector
    gsl_matrix *X; // observation matrix

    size_t n; // sample size;
    size_t p; // number of parameters (including intercept)

    double psi; // dispersion

    bool free_data; // depend on different constructors

    gsl_vector *beta;
    gsl_vector *Jbeta;
    gsl_vector *pvalue;

  public:
    logit(gsl_vector *yv, gsl_matrix *Xv);
    logit(const std::vector<double> & yv, const std::vector<std::vector<double> > Xv);
    ~logit();
    void fit();
    double calculate_pi(gsl_vector* xi) const;
    double calculate_err(const gsl_vector* beta2) const;
    double calculate_err() const;
    int calculate_J(gsl_matrix* J) const;
    int calculate_U(gsl_vector* U) const;
    void display(gsl_matrix* m) const;
    void displayv(gsl_vector* v) const;
}; //logit

#endif // _LOGIT_H

ext: h


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include <gsl/gsl>

#include "core/optimizer/graph_transformer.h"
#include "orttraining/core/optimizer/graph_transformer_config.h"
#include "orttraining/core/session/training_session.h"

namespace onnxruntime {
struct FreeDimensionOverride;

namespace training {
namespace transformer_utils {

/** Generates all pre-training transformers for this level. */
std::vector<std::unique_ptr<GraphTransformer>> GeneratePreTrainingTransformers(
    TransformerLevel level,
    const std::unordered_set<std::string>& weights_to_train,
    const TrainingGraphTransformerConfiguration& config,
    const IExecutionProvider& execution_provider,  // required for constant folding
    const std::unordered_set<std::string>& rules_and_transformers_to_disable = {});

/** Generates all predefined (both rule-based and non-rule-based) transformers for this level.
    If transformers_and_rules_to_enable is not empty, it returns the intersection between the predefined transformers/rules
    and the transformers_and_rules_to_enable. */
InlinedVector<std::unique_ptr<GraphTransformer>> GenerateTransformers(
    TransformerLevel level,
    const std::unordered_set<std::string>& weights_to_train,
    gsl::span<const FreeDimensionOverride> free_dimension_overrides,
    const InlinedHashSet<std::string>& rules_and_transformers_to_disable = {});

}  // namespace transformer_utils
}  // namespace training
}  // namespace onnxruntime

max_stars_repo_name: SiriusKY/onnxruntime

max_issues_repo_name: SiriusKY/onnxruntime

max_forks_repo_name: SiriusKY/onnxruntime


text:
/*
Copyright (c) 2003-2008 Rudi Cilibrasi, Rulers of the RHouse
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
2. Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
3. Neither the name of the University nor the names of its contributors
   may be used to endorse or promote products derived from this software
   without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE RULERS AND CONTRIBUTORS ``AS IS'' AND
ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED.  IN NO EVENT SHALL THE RULERS 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.
*/
#include <stdio.h>
#include <string.h>
#if OPENMP_ENABLED
#include <omp.h>
#endif
#include <assert.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <stdlib.h>
#include <glib.h>
#include <glib/gstdio.h>
//#include <gsl/gsl_cblas.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <libintl.h>
#include <unistd.h>

#include "complearn/util.h"

#define _(O) gettext(O)

void complearn_mrbreaker(void)
{
  printf("break!\n");
}

int complearn_make_directory_if_necessary(const char *name)
{
  return g_mkdir_with_parents(name, 0700);
}

char *complearn_get_hostname()
{
  static char *answer;
  if (answer == NULL) {
    static char buf[1024];
#ifdef __MINGW32__
	strcpy(buf, "localhost");
#else
    gethostname(buf, sizeof(buf));
#endif
    buf[sizeof(buf)-1] = '\0';
    answer = buf;
  }
  return g_strdup(answer);
}

int complearn_get_pid(void)
{
  return getpid();
}

GSList *complearn_read_directory_of_files(const char *dirname)
{
  GError *err = NULL;
  GSList *result = NULL;
  GDir *d = g_dir_open(dirname, 0, &err);
  if (err)
    g_error(_("Cannot open directory %s"), dirname);
  const char *n;
  while ( (n = g_dir_read_name(d)) ) {
    char *cur = g_build_filename(dirname, n, NULL);
    if (g_file_test(cur, G_FILE_TEST_IS_REGULAR)) {
      result = g_slist_append(result, complearn_new_arg(complearn_read_whole_file(cur),cur));
    }
    g_free(cur);
  }
  g_dir_close(d);
  return result;
}

static GSList *complearn_read_list_of_files_or_strings(const char *filename, int isfiles)
{
  GIOChannel *in;
  GError *err = NULL;
  GSList *result = NULL;
  in = g_io_channel_new_file(filename, "r", &err);
  if (err)
    g_error(_("Cannot open file %s"), filename);
  for (;;) {
    GIOStatus s;
    gsize tpos;
    err = NULL;
    GString *block = g_string_new("");
    s = g_io_channel_read_line_string(in, block, &tpos, &err);
    if (s == G_IO_STATUS_EOF)
      break;
    if (err != NULL)
      g_error(_("Cannot read line from file %s: %s"), filename, err->message);
    g_string_truncate(block, tpos);
    GString *tokeep;
    tokeep = (isfiles != 0) ? complearn_read_whole_file(block->str) : g_string_new_len(block->str, block->len);
    result = g_slist_append(result, complearn_new_arg(tokeep, block->str));
  }
  g_io_channel_close(in);
  return result;
}

GSList *complearn_read_list_of_files(const char *filename)
{
  return complearn_read_list_of_files_or_strings(filename, 1);
}

GSList *complearn_read_list_of_strings(const char *filename)
{
  return complearn_read_list_of_files_or_strings(filename, 0);
}

GString *complearn_read_whole_file_ptr(FILE *fp)
{
  GString *result = g_string_new("");
  char buf[512];
  int rl;
  while ( (rl = fread(buf, 1, 512, fp)) == 512)
    g_string_append_len(result, buf, rl);
  if (rl > 0)
    g_string_append_len(result, buf, rl);
  return result;
}

char *complearn_chdir(const char *newdir)
{
  char *retval = g_get_current_dir();
  int rv;
  rv = g_chdir(newdir);
  if (rv != 0) {
    g_error(_("Cannot change directory to %s"), newdir);
  }
  return retval;
}

int complearn_remove_directory_recursively(const char *dirname, int f_err_out)
{
  GError *err = NULL;
  int result = 0;
  GDir *d = g_dir_open(dirname, 0, &err);
  if ( (err) ) {
    if (f_err_out)
      g_error(_("Cannot open directory %s. Error: %s"), dirname, err->message);
    else
      return 1;
  }
  const char *n;
  while ( (n = g_dir_read_name(d)) ) {
    char *cur = g_build_filename(dirname, n, NULL);
    if (g_file_test(cur, G_FILE_TEST_IS_REGULAR))
      g_unlink(cur);
    g_free(cur);
  }
  g_dir_close(d);
  result = g_rmdir(dirname);
  return result;
}

int complearn_write_file(const char *fname, const GString *f)
{
  g_assert(fname != NULL && strlen(fname) > 0);
  FILE *fp = fopen(fname, "wb");
  if (fp == NULL) {
    g_error(_("Cannot open file %s for writing."), fname);
  }
  fwrite(f->str, 1, f->len, fp);
  fclose(fp);
  return 0;
}

GString *complearn_read_whole_file(const char *fname)
{
  FILE *fp = fopen(fname, "rb");
  if (fp == NULL)
    g_error(_("Cannot open file %s for reading"), fname);
  fseek(fp, 0, SEEK_END);
  int sz = ftell(fp);
  fseek(fp, 0, SEEK_SET);
  char *b = calloc(sz, 1);
  GString *result;
  fread(b, 1, sz, fp);
  fclose(fp);
  result = g_string_new_len(b, sz);
  g_free(b);
  return result;
}

char *complearn_make_temp_dir(void)
{
  char *result = complearn_make_temp_dir_name();
#ifdef __MINGW32__
  mkdir(result);
#else
  mkdir(result, 0700);
#endif
  return result;
}

/* different each time */
char *complearn_make_temp_file_name(void)
{
  char *tnam, *result;
  static int i = 0;
  char pidbuf[32];
  sprintf(pidbuf, "%d-%d", ++i, getpid());
  tnam = getenv("TMPDIR");
  if (tnam == NULL)
    tnam = "/tmp";
  result = g_strconcat(tnam, "/clfile", pidbuf, NULL);
  return result;
}

char *complearn_make_temp_dir_name(void)
{
  int tid = 0;
#if OPENMP_ENABLED
    if (!g_thread_supported ()) g_thread_init (NULL);
    tid = omp_get_thread_num();
#endif

  char *result = NULL, *tnam;
  int incr = 0;
  char pidbuf[64];
  for (;;) {
    sprintf(pidbuf, "%dx%dx%d", getpid(),incr,tid);
    tnam = getenv("TMPDIR");
    if (tnam == NULL)
      tnam = "/tmp";
    result = tnam;
    result = g_strconcat(tnam, "/complearn-", pidbuf, NULL);
    if (!g_file_test(result, G_FILE_TEST_EXISTS) && !g_file_test(result, G_FILE_TEST_IS_DIR))
      break;
    incr += 1;
  }
  return result;
}

char ** complearn_dupe_strings(const char * const * str)
{
  int size = complearn_count_strings(str);
  char **result;
  int i;
  result = calloc(size+1, sizeof(gpointer));
  for (i = 0; i < size; i += 1)
    result[i] = g_strdup(str[i]);
  return result;
}

int complearn_count_strings(const char * const * str)
{
  int acc = 0;
  const char * const *cur = str;
  while (cur && *cur) {
    acc += 1;
    cur += 1;
  }
  return acc;
}

CLArgumentEntry *complearn_new_arg(GString *block, const char *label) {
  CLArgumentEntry *cl = calloc(1, sizeof(*cl));
  cl->block = block;
  cl->label = g_strdup(label);
  return cl;
}

/* This function makes the off-diagonal elements symmetric by averaging.  */
gsl_matrix *complearn_average_matrix(gsl_matrix *a)
{
  int i,j;
  gsl_matrix *res;
  res = gsl_matrix_alloc(a->size1, a->size2);

  for (i = 0;i<a->size1;i+=1)
    for (j = 0; j <a->size2; j+=1) {
      double v;
      if (i == j)
        gsl_matrix_set(res,i,j,gsl_matrix_get(a,i,j));
      if (i <= j)
        continue;
      v = (gsl_matrix_get(a,i,j) + gsl_matrix_get(a,j,i))/2.0;
      gsl_matrix_set(res,i,j,v);
      gsl_matrix_set(res,j,i,v);
    }
  return res;
}

gsl_matrix *complearn_svd_project(gsl_matrix *a)
{
  int retval;
  gsl_matrix *res;
  gsl_matrix *u, *v;
  gsl_vector *s;
  u = gsl_matrix_alloc(a->size1, a->size2);
  gsl_matrix_memcpy(u, a);
  v = gsl_matrix_alloc(a->size2, a->size2);
  s = gsl_vector_alloc(a->size1);
  retval = gsl_linalg_SV_decomp_jacobi(u, v, s);
//  g_assert(retval == GSL_OK);
  res = gsl_matrix_alloc(a->size1, a->size2);
  gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, a, u, 0.0, res);
  return res;
}

static gboolean is_last_common_char(char **inp, int which)
{
  int i;
  if (strlen(inp[0]) <= which)
    return FALSE;
  for (i = 1; inp[i]; i += 1)
    if (strlen(inp[i]) <= which || inp[i][strlen(inp[i])-which-1] != inp[0][strlen(inp[0])-which-1])
      return FALSE;
  return TRUE;
}

static gboolean is_common_char(char **inp, int which)
{
  int i;
  if (strlen(inp[0]) <= which)
    return FALSE;
  for (i = 1; inp[i]; i += 1)
    if (strlen(inp[i]) <= which || inp[i][which] != inp[0][which])
      return FALSE;
  return TRUE;
}

static char **chop_strings(char **inp, int prefixc, int suffixc)
{
  int sc = complearn_count_strings((const char * const *) inp);
  char **result = calloc(sc+1,sizeof(gpointer));
  int i;
  for (i = 0; i < sc; i += 1) {
    char *p = inp[i]+prefixc;
    char *res;
    res = g_strdup(p);
    if (suffixc && strlen(p) > suffixc)
      res[strlen(res)-suffixc] = '\0';
    result[i] = res;
  }
  return result;
}

char ** complearn_fix_labels(char **inp)
{
  int prefixc = 0, suffixc = 0;
  if (inp[1] != NULL) {
    while (is_common_char(inp, prefixc))
      prefixc += 1;
    while ( is_last_common_char(inp, suffixc))
      suffixc += 1;
  }
  return chop_strings(inp, prefixc, suffixc);
}

ext: c

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
/* specfunc/bessel_Y1.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_trig.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "bessel.h"
#include "bessel_amp_phase.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC besy1, 1977 version, w. fullerton */

/* chebyshev expansions

 series for by1        on the interval  0.          to  1.60000d+01
                                        with weighted error   1.87e-18
                                         log weighted error  17.73
                               significant figures required  17.83
                                    decimal places required  18.30
*/

static double by1_data[14] = {
  0.03208047100611908629,
  1.262707897433500450,
  0.00649996189992317500,
 -0.08936164528860504117,
  0.01325088122175709545,
 -0.00089790591196483523,
  0.00003647361487958306,
 -0.00000100137438166600,
  0.00000001994539657390,
 -0.00000000030230656018,
  0.00000000000360987815,
 -0.00000000000003487488,
  0.00000000000000027838,
 -0.00000000000000000186
};
static cheb_series by1_cs = {
  by1_data,
  13,
  -1, 1,
  10
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_Y1_e(const double x, gsl_sf_result * result)
{
  const double two_over_pi = 2.0/M_PI;
  const double xmin = 1.571*GSL_DBL_MIN; /*exp ( amax1(alog(r1mach(1)), -alog(r1mach(2)))+.01) */
  const double x_small = 2.0 * GSL_SQRT_DBL_EPSILON;
  const double xmax    = 1.0/GSL_DBL_EPSILON;

  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < xmin) {
    OVERFLOW_ERROR(result);
  }
  else if(x < x_small) {
    const double lnterm = log(0.5*x);
    gsl_sf_result J1;
    gsl_sf_result c;
    int status = gsl_sf_bessel_J1_e(x, &J1);
    cheb_eval_e(&by1_cs, -1.0, &c);
    result->val = two_over_pi * lnterm * J1.val + (0.5 + c.val)/x;
    result->err = fabs(lnterm) * (fabs(GSL_DBL_EPSILON * J1.val) + J1.err) + c.err/x;
    return status;
  }
  else if(x < 4.0) {
    const double lnterm = log(0.5*x);
    int status;
    gsl_sf_result J1;
    gsl_sf_result c;
    cheb_eval_e(&by1_cs, 0.125*x*x-1.0, &c);
    status = gsl_sf_bessel_J1_e(x, &J1);
    result->val = two_over_pi * lnterm * J1.val + (0.5 + c.val)/x;
    result->err = fabs(lnterm) * (fabs(GSL_DBL_EPSILON * J1.val) + J1.err) + c.err/x;
    return status;
  }
  else if(x < xmax) {
    const double z = 32.0/(x*x) - 1.0;
    gsl_sf_result ca;
    gsl_sf_result ct;
    gsl_sf_result cp;
    const int stat_ca = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bm1_cs,  z, &ca);
    const int stat_ct = cheb_eval_e(&_gsl_sf_bessel_amp_phase_bth1_cs, z, &ct);
    const int stat_cp = gsl_sf_bessel_cos_pi4_e(x, ct.val/x, &cp);
    const double sqrtx = sqrt(x);
    const double ampl  = (0.75 + ca.val) / sqrtx;
    result->val  = -ampl * cp.val;
    result->err  = fabs(cp.val) * ca.err/sqrtx + fabs(ampl) * cp.err;
    result->err += GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_3(stat_ca, stat_ct, stat_cp);
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_Y1(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_Y1_e(x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#ifndef _BCQL_H_
#define _BCQL_H_

#include <cblas.h>
#include <lapacke/lapacke.h>
#include <math.h>

//TODO Only one order is currently supported...
enum BCSK_Q_ORDER { BcskWXYZ = 130 };
enum BCSK_V_ORDER { BcskXYZ  = 140 };
enum BCSK_M_ORDER { BcskColMajor = 102 }; //TODO Add Row Major

//copies first "len" entries of src to corresponding entries in dst
int dDub ( double *dst , double *src , unsigned int len ) ;

// Sets first "len" entries of array to 0;
int dZeros ( double *array , unsigned int len ) ;

// puts in dst the element-by-element sum of src1 and src2
int dSum2  ( double *dst , double alpha , double *src1 , double beta , double *src2 , unsigned int len ) ;

// sum up 3 arrays
int dSum3 ( double *dst ,
            double a1 , double *src1 ,
            double a2 , double *src2 ,
            double a3 , double *src3 ,
            unsigned int len ) ;

// sum up 4 arrays
int dSum4 ( double *dst ,
            double a1 , double *src1 ,
            double a2 , double *src2 ,
            double a3 , double *src3 ,
            double a4 , double *src4 ,
            unsigned int len ) ;

// sum up 5 arrays
int dSum5 ( double *dst ,
            double a1 , double *src1 ,
            double a2 , double *src2 ,
            double a3 , double *src3 ,
            double a4 , double *src4 ,
            double a5 , double *src5 ,
            unsigned int len ) ;

// Compute cross product between 2 3D vecs
// c = a cross b ;
int dCross ( double *c , double *a , double *b ) ;

// Set the square matrix mat to the identity matrix
int dEye ( unsigned int order , double scale , double *mat ) ;

// Transpose square matrix
int dSqTr ( int n , double *src , double *dst ) ;

// Stores in quaternion c the hamilton product between quaternions a and b
int dHamilton (const enum BCSK_Q_ORDER Order , double *c , double *a , double *b ) ;

// Store in b the inverse of quaternion a
int dQInv( const enum BCSK_Q_ORDER Order , double *b , double *a) ;

//TODO vec arma_quat_between_vecs(const vec& a,const vec& b);

/***********************************************
* Compute angular velocity given the orientation quaternion
* and its derivative.
***********************************************/
//int dQuatInv( const enum LIBQUAT_ORDER BOrder , double *v , const enum LIBQUAT_ORDER AOrder , double *a) ;
int dQD2Vel ( const enum BCSK_V_ORDER v_order    ,
              const enum BCSK_Q_ORDER q_order    ,
              double *v                          ,
              double *q                          ,
              double *dqdt                       );

/*TODO**********************************************
* Compute angular accelleration given the orientation quaternion
* and its derivatives.
***********************************************/
/*vec arma_quat_d_to_acc(const vec& q0,const vec& q1,const vec& q2){
	vec w1(4);
	w1 = 2 * arma_quat_hamilton( q2-arma_quat_hamilton(arma_quat_hamilton(q1,arma_quat_inv(q0)),q1) , arma_quat_inv(q0) );
	return(w1.subvec(1,3));
}
*/

/***********************************************
* Rotate vector v by rotation q.
* This formula should be faster than:
* v_new = q * v * q^-1
***********************************************/
//TODO vec arma_quat_rot ( const vec& q,const vec& v ) ;

/***********************************************
* Return matrix M from vector v s.t. for any x
* M*x = cross(v,x)
***********************************************/
int dCrossMat ( const enum BCSK_M_ORDER MOrder , double *M , const double *v ) ;

/***********************************************
* Applyes Rodriguez formula to convert quaternion q
* to rotation matrix R
***********************************************/
int dQ2R ( const enum BCSK_M_ORDER Order , const enum BCSK_Q_ORDER , double *R , double *q) ;

/***********************************************
* Conversion from rotation matrix R to quaternion q
***********************************************/
int dR2Q ( const enum BCSK_M_ORDER ROrder , const enum BCSK_Q_ORDER QOrder , double *R , double *q) ;

/* Saturate each elements of an array between -abs(sat_level) and +abs(sat_level)
 *
 *                _______  _  +abs(sat_level)
 *             | /
 *      _______|/_______y  _ 0
 *             /
 *      ______/|           _ -abs(sat_level)
 *             x
 *
 *     x: input value
 *     y: output value
 *
 * N.B.
 */

int dSat( double *src , double *dst , double sat_level , int len ) ;

// normalize vector X
int dNormalise ( const int N , double *X  );

// compute rotation quaternion between v1 and v2
int dQbetVs ( double *Q , double *v1 , double *v2);

#endif


text:
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <libiomp/omp.h>
#include <string.h>
#include <stdbool.h>

#ifdef __APPLE__

#include <cblas.h>

#define set_num_threads(x) openblas_set_num_threads(x)
#define get_num_threads() openblas_get_num_threads()
#else
#include <mkl.h>
#define set_num_threads(x) mkl_set_num_threads(x)
#define get_num_threads() mkl_get_num_threads()
#endif

int main(int argc, char **argv) {

    char *filename = NULL;
    int option = 0, omp_threads = 1;
    while ((option = getopt(argc, argv, "f:t:")) != -1) {
        switch (option) {
            case 'f':
                filename = optarg;
                break;
            case 't':
                omp_threads = atoi(optarg);
                break;
            default:
                printf("Usage: mpi_cg -f string -t num_threads \n");
                return 0;
        }
    }
    if (filename == NULL) {
        printf("Usage: mpi_cg -f string -t num_threads \n");
        return 0;
    }

    set_num_threads(omp_threads);
    omp_set_num_threads(omp_threads);

    int numTasks, rank;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &numTasks);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int qProcs = (int) sqrt(numTasks);
    if (qProcs - sqrt(numTasks) != 0) {
        if (rank == 0)
            printf("np must be a square number\n");
        MPI_Finalize();
        return 0;
    }

    int n, dims[2] = {qProcs, qProcs}, periods[2] = {0, 0}, reorder = 0, coords[2], rc;
    MPI_Comm cartComm;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, reorder, &cartComm);
    MPI_Comm_rank(cartComm, &rank);
    MPI_Cart_coords(cartComm, rank, 2, coords);

    MPI_File matFile;
    rc = MPI_File_open(cartComm, filename, MPI_MODE_RDONLY, MPI_INFO_NULL, &matFile);
    if (rc && (rank == 0)) {
        printf("Unable to open file %s\n", filename);
        fflush(stdout);
    }

    MPI_Status status;
    MPI_File_read(matFile, &n, 1, MPI_INT, &status);
    int blockSize = (n + qProcs - 1) / qProcs;  /* number of rows in _block_ */
    int lastBlock = n - (qProcs - 1) * blockSize;

    double *LocalA = (double *) malloc(blockSize * blockSize * sizeof(double));

    MPI_Datatype readFileType;
    MPI_Type_vector(blockSize, blockSize, n, MPI_DOUBLE, &readFileType);
    MPI_Type_commit(&readFileType);
    MPI_Offset offset = (MPI_Offset) (1 * sizeof(int) +
                                      sizeof(double) * (blockSize * coords[1] + blockSize * n * coords[0]));
    MPI_File_set_view(matFile, offset, MPI_DOUBLE, readFileType,
                      "native", MPI_INFO_NULL);

    MPI_File_read_at(matFile, 0, LocalA, blockSize * blockSize, MPI_DOUBLE, &status);
    MPI_Type_free(&readFileType);
    if (coords[0] == qProcs - 1) {
        memset(LocalA + lastBlock * blockSize, 0, (blockSize - lastBlock) * blockSize * sizeof(double));
    }
    if (coords[1] == qProcs - 1) {
        int kk = 0;
        for (kk = 0; kk < blockSize; ++kk) {
            memset(LocalA + kk * blockSize + lastBlock, 0, (blockSize - lastBlock) * sizeof(double));
        }
    }

    double *LocalB = (double *) malloc(blockSize * sizeof(double));
    memset(LocalB, 0, blockSize * sizeof(double));
    offset = (MPI_Offset) (1 * sizeof(int) + sizeof(double) * n * n);
    MPI_File_set_view(matFile, offset, MPI_DOUBLE, MPI_DOUBLE,
                      "native", MPI_INFO_NULL);
    offset = coords[1] * blockSize;
    MPI_File_read_at_all(matFile, offset, LocalB, coords[1] == qProcs - 1 ? lastBlock : blockSize, MPI_DOUBLE, &status);

    MPI_File_close(&matFile);
    double mpi_start = MPI_Wtime();
    MPI_Comm rowComm, colComm;
    MPI_Comm_split(MPI_COMM_WORLD, coords[0], rank, &rowComm);
    MPI_Comm_split(MPI_COMM_WORLD, coords[1], rank, &colComm);

    double tolerance = 1e-8;
    double alpha, beta, rho_new = 0.0, rho_old = 0.0;
    int iteration = 0;
    int world_root = numTasks - 1;
    double *loc_Vr = (double *) malloc(blockSize * sizeof(double));
    double *loc_Vp = (double *) malloc(blockSize * sizeof(double));
    double *loc_Vw = (double *) malloc(blockSize * sizeof(double));
    double *loc_Vx = (double *) malloc(blockSize * sizeof(double));

    int row_rank, col_rank;
    MPI_Comm_rank(colComm, &col_rank);
    MPI_Comm_rank(rowComm, &row_rank);

    //localr <- localB
    memset(loc_Vx, 0, blockSize * sizeof(double));
    cblas_dcopy(blockSize, LocalB, 1, loc_Vr, 1);
    double b_norm2 = 0;
    if (col_rank == qProcs - 1) {
        double loc_b_dot = cblas_ddot((row_rank == qProcs - 1 ? lastBlock : blockSize), loc_Vr, 1, loc_Vr, 1);
        MPI_Reduce(&loc_b_dot, &rho_new, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
        b_norm2 = sqrt(rho_new);
    }

    int loopBreak = 0;
    while (true) {
        if (rank == world_root) {
            loopBreak = (sqrt(rho_new) / b_norm2 < tolerance) ? 1 : 0;
        }
        MPI_Bcast(&loopBreak, 1, MPI_INT, world_root, MPI_COMM_WORLD);
        if (loopBreak > 0)
            break;

        iteration++;
        if (iteration == 1) {
            cblas_dcopy(blockSize, loc_Vr, 1, loc_Vp, 1);
        } else {
            if (rank == world_root) {
                beta = rho_new / rho_old;
            }
            MPI_Bcast(&beta, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);
            cblas_dscal((row_rank == qProcs - 1 ? lastBlock : blockSize), beta, loc_Vp, 1);
            cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), 1.0, loc_Vr, 1, loc_Vp, 1);
        }

        cblas_dgemv(CblasRowMajor, CblasNoTrans, (col_rank == qProcs - 1 ? lastBlock : blockSize),
                    (row_rank == qProcs - 1 ? lastBlock : blockSize),
                    1.0, LocalA, blockSize, loc_Vp, 1, 0.0, loc_Vw, 1);

        if (col_rank == row_rank) {
            MPI_Reduce(MPI_IN_PLACE, loc_Vw, blockSize, MPI_DOUBLE, MPI_SUM, col_rank, rowComm);
        } else {
            MPI_Reduce(loc_Vw, loc_Vw, blockSize, MPI_DOUBLE, MPI_SUM, col_rank, rowComm);
        }

        MPI_Bcast(loc_Vw, (row_rank == qProcs - 1 ? lastBlock : blockSize), MPI_DOUBLE, row_rank, colComm);

        if (col_rank == qProcs - 1) {
            double loc_pw_dot = cblas_ddot((row_rank == qProcs - 1 ? lastBlock : blockSize), loc_Vp, 1, loc_Vw,
                                           1), glo_pw_dot;
            MPI_Reduce(&loc_pw_dot, &glo_pw_dot, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
            if (rank == world_root) {
                alpha = rho_new / glo_pw_dot;
            }
        }
        MPI_Bcast(&alpha, 1, MPI_DOUBLE, world_root, MPI_COMM_WORLD);

        cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), -alpha, loc_Vw, 1, loc_Vr, 1);

        if (coords[0] == qProcs - 1) {
            rho_old = rho_new;
            cblas_daxpy((row_rank == qProcs - 1 ? lastBlock : blockSize), alpha, loc_Vp, 1, loc_Vx, 1);
            double loc_r_dot = cblas_ddot(blockSize, loc_Vr, 1, loc_Vr, 1);

            MPI_Reduce(&loc_r_dot, &rho_new, 1, MPI_DOUBLE, MPI_SUM, qProcs - 1, rowComm);
        }
    }

    if (rank == world_root) {
        printf("mpi_cg_it: %d\t", iteration);
        printf("mpi_cg_time: %f\n", MPI_Wtime() - mpi_start);

    }
    MPI_Comm_free(&cartComm);
    MPI_Comm_free(&colComm);
    MPI_Comm_free(&rowComm);

    if (col_rank == qProcs - 1) {
        for (rc = 0; rc < (row_rank == qProcs - 1 ? lastBlock : blockSize); ++rc) {
            if (fabs(loc_Vx[rc] - 1) > 1e-3) {
                printf("MPI_CG_ERR %f \n", loc_Vx[rc]);
                break;
            }
        }
    }
    MPI_Finalize();
}

max_stars_repo_name: baishuai/MEXP

max_issues_repo_name: baishuai/MEXP

max_forks_repo_name: baishuai/MEXP


text:
#pragma once
#include <halley/maths/vector2.h>
#include <halley/text/halleystring.h>
#include <gsl/gsl>

namespace Halley
{
	class FontFacePimpl;
	class Image;

	struct FontMetrics
	{
		Vector2f advance;
		Vector2f bearingHorizontal;
		Vector2f bearingVertical;
	};

	struct KerningPair
	{
		int left;
		int right;
		Vector2f kerning;

		KerningPair(int left, int right, Vector2f kerning) : left(left), right(right), kerning(kerning) {}
	};

	class FontFace
	{
	public:
		explicit FontFace(String filename);
		explicit FontFace(gsl::span<const gsl::byte> data);
		~FontFace();

		void setSize(float size);
		String getName() const;
		float getSize() const;
		float getHeight() const;
		float getAscender() const;

		Vector<int> getCharCodes() const;
		Vector2i getGlyphSize(int charCode) const;

		void drawGlyph(Image& image, int charcode, Vector2i pos) const;
		FontMetrics getMetrics(int charcode, float scale = 1.0f) const;

		Vector<KerningPair> getKerning(const Vector<int>& codes) const;

	private:
		std::unique_ptr<FontFacePimpl> pimpl;
		float size = 0;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


text:
/*
    COPYRIGHT DISCLAIMER

    Vincent Le Guilloux, Peter Schmidtke and Pierre Tuffery, hereby
	disclaim all copyright interest in the program “fpocket” (which
	performs protein cavity detection) written by Vincent Le Guilloux and Peter
	Schmidtke.

    Vincent Le Guilloux  28 November 2008
    Peter Schmidtke      28 November 2008
    Pierre Tuffery       28 November 2008

    GNU GPL

    This file is part of the fpocket package.

    fpocket 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.

    fpocket 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 fpocket.  If not, see <http://www.gnu.org/licenses/>.

**/

#ifndef DH_UTILS
#define DH_UTILS

/* ------------------------------ INCLUDES ---------------------------------- */

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <ctype.h>
#include <limits.h>
#include <time.h>

#ifdef MD_USE_GSL 	/* GSL */

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

#endif 				/* /GSL */

#include "memhandler.h"

/* ------------------------------- PUBLIC MACROS ---------------------------- */

#define M_MAX_PDB_NAME_LEN 200  /**< maximum pdb filename length*/

#define M_SIGN 1
#define M_NO_SIGN 0

#ifdef MD_USE_GSL	/* GSL */

#define M_GEN_MTWISTER gsl_rng_mt19937
#define M_GEN_GFSR4 gsl_rng_gfsr4
#define M_GEN_TAUS gsl_rng_taus2
#define M_GEN_RANLXS0 gsl_rng_ranlxs0
#define M_GEN_RANLXS1 gsl_rng_ranlxs1

#endif				/* /GSL */

/* ------------------------------ PUBLIC STRUCTURES ------------------------- */

typedef struct tab_str
{
	char **t_str ;
	int nb_str ;

} tab_str ;

/* ------------------------------- PROTOTYPES ------------------------------- */

int str_is_number(const char *str, const int sign) ;
int str_is_float(const char *str, const int sign) ;
void str_trim(char *str) ;
tab_str* str_split(const char *str, const int sep) ;

tab_str* f_readl(const char *fpath, int nchar_max) ;
void free_tab_str(tab_str *tstr) ;
void print_tab_str(tab_str* strings) ;

int in_tab(int *tab, int size, int val) ;
int index_of(int *tab, int size, int val)  ;

void remove_ext(char *str) ;
void remove_path(char *str) ;
void extract_ext(char *str, char *dest) ;
void extract_path(char *str, char *dest) ;

void start_rand_generator(void) ;
float rand_uniform(float min, float max) ;

FILE* fopen_pdb_check_case(char *name, const char *mode)  ;

float float_get_min_in_2D_array(float **t,size_t n,int col);
float float_get_max_in_2D_array(float **t,size_t n,int col);

#endif


text:
/// \file angles.c
/// \brief Optimize polar angles with respect to Coulomb repulsion, given the radii.
//
//	Copyright (c) 2014, Christian B. Mendl
//	All rights reserved.
//	http://christian.mendl.net
//
//	This program is free software; you can redistribute it and/or
//	modify it under the terms of the Simplified BSD License
//	http://www.opensource.org/licenses/bsd-license.php
//
//	Reference:
//	  Christian B. Mendl, Francesc Malet, Paola Gori-Giorgi
//	  Wigner localization in quantum dots from Kohn-Sham density functional theory without symmetry breaking
//	  Physical Review B 89, 125106 (2014)
//	  (preprint http://arxiv.org/abs/1311.6011)
//_______________________________________________________________________________________________________________________
//

#include "angles.h"
#include "coulomb.h"
#include "util.h"
#define _USE_MATH_DEFINES
#include <math.h>
#include <stdlib.h>
#include <memory.h>
#include <float.h>
#include <assert.h>
#include <stdio.h>
// LAPACK functions
#include <acml.h>
// GSL library
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_siman.h>

static const int maxiter		= 64;			//!< maximum number of iterations for local angle optimization
static const double step_size	= 0.01;			//!< step size for angle optimization

//_______________________________________________________________________________________________________________________
//
// GSL wrappers of Vee() for local minimization

/// \brief GSL parameters for angle optimization
typedef struct
{
	double phi0;
	const double *r;
}
gsl_angle_params_t;

/// \brief GSL wrapper for Vee()
static double gslVee(const gsl_vector *phiRest, void *p)
{
	int i;
	const gsl_angle_params_t *params = (gsl_angle_params_t *)p;
	double v_val;

	// first angle is held fixed
	const int nelec = (int)phiRest->size + 1;

	double *phi = (double *)malloc(nelec*sizeof(double));
	phi[0] = params->phi0;
	for (i = 1; i < nelec; i++) {
		phi[i] = gsl_vector_get(phiRest, i-1);
	}

	v_val = Vee(nelec, params->r, phi);

	free(phi);

	return v_val;
}

/// \brief GSL wrapper for VeeGradPhi()
static void gslVeeGradPhi(const gsl_vector *phiRest, void *p, gsl_vector *df)
{
	int i;
	const gsl_angle_params_t *params = (gsl_angle_params_t *)p;
	double *grad;

	// first angle is held fixed
	const int nelec = (int)phiRest->size + 1;

	double *phi = (double *)malloc(nelec*sizeof(double));
	phi[0] = params->phi0;
	for (i = 1; i < nelec; i++) {
		phi[i] = gsl_vector_get(phiRest, i-1);
	}

	grad = (double *)malloc(nelec*sizeof(double));

	// calculate gradient
	VeeGradPhi(nelec, params->r, phi, grad);

	// ignore derivative with respect to first angle (held fixed)
	for (i = 1; i < nelec; i++) {
		gsl_vector_set(df, i-1, grad[i]);
	}

	free(grad);
	free(phi);
}

/// \brief GSL wrapper for both 'Vee()' and 'VeeGradPhi()' together
static void gslVeeValGradPhi(const gsl_vector *phiRest, void *p, double *f, gsl_vector *df)
{
	int i;
	const gsl_angle_params_t *params = (gsl_angle_params_t *)p;
	double *phi;
	double *grad;

	// first angle is held fixed
	const int nelec = (int)phiRest->size + 1;

	phi = (double *)malloc(nelec*sizeof(double));
	phi[0] = params->phi0;
	for (i = 1; i < nelec; i++) {
		phi[i] = gsl_vector_get(phiRest, i-1);
	}

	*f = Vee(nelec, params->r, phi);

	grad = (double *)malloc(nelec*sizeof(double));

	// calculate gradient
	VeeGradPhi(nelec, params->r, phi, grad);

	// ignore derivative with respect to first angle (which is held fixed)
	for (i = 1; i < nelec; i++) {
		gsl_vector_set(df, i-1, grad[i]);
	}

	free(grad);
	free(phi);
}

//_______________________________________________________________________________________________________________________
///
/// \brief Optimize polar angles with respect to Coulomb repulsion, given the radii
///
/// \param nelec number of electrons
/// \param r radii, array of length 'nelec'
/// \param phi serves as input (starting values) and output (optimized angles); array of length 'nelec';
///        phi[0] is not changed due to invariance under global rotation
///
/// \return detailed error description, if any
///
error_desc_t OptimizeAngles(const int nelec, const double *r, double *phi)
{
	int i, j, k;
	error_desc_t err = { 0 };
	int status;

	// tolerance for convergence test based on gradient
	const double tol = nelec < 7 ? 1e-12 : 1e-10;

	const gsl_multimin_fdfminimizer_type *T = gsl_multimin_fdfminimizer_vector_bfgs2;	// alternative: gsl_multimin_fdfminimizer_conjugate_fr
	gsl_multimin_fdfminimizer *s = gsl_multimin_fdfminimizer_alloc(T, nelec-1);

	gsl_angle_params_t params;
	params.phi0 = phi[0];
	params.r = r;

	gsl_multimin_function_fdf my_func;
	my_func.n = nelec-1;
	my_func.f = gslVee;
	my_func.df = gslVeeGradPhi;
	my_func.fdf = gslVeeValGradPhi;
	my_func.params = &params;

	// starting point
	// omit first angle phi[0]
	gsl_vector *x = gsl_vector_alloc(nelec-1);
	for (i = 1; i < nelec; i++) {
		gsl_vector_set(x, i-1, phi[i]);
	}

	gsl_multimin_fdfminimizer_set(s, &my_func, x, step_size, 0.1);

	for (k = 0; k < maxiter; k++)
	{
		status = gsl_multimin_fdfminimizer_iterate(s);

		// interval [-pi, pi]
		for (i = 0; i < nelec-1; i++)
		{
			double curphi = gsl_vector_get(s->x, i);
			assert(!isinf(curphi));

			if (curphi < -M_PI)		// using "if" instead of "while" to avoid infinite loop if curphi is +-1.#INF
			{
				curphi += 2*M_PI;
			}
			else if (curphi > M_PI)
			{
				curphi -= 2*M_PI;
			}

			gsl_vector_set(s->x, i, curphi);
		}

		// check status
		if (status != 0) {
			err.codes |= E_ITERATION;
			err.msg = gsl_strerror(status);
			fprintf(stderr, "Warning: 'gsl_multimin_fdfminimizer_iterate' failed at iteration %i, error: %s.\n", k, err.msg);
			break;
		}

		// convergence test
		status = gsl_multimin_test_gradient(s->gradient, tol);
		if (status == GSL_SUCCESS) {
			printf("'gsl_multimin_test_gradient()' successful.\n");
			break;
		}
	}

	// norm of gradient
	double delta = 0;
	for (i = 0; i < nelec-1; i++)
	{
		delta += square(gsl_vector_get(s->gradient, i));
	}
	// 2-norm
	delta = sqrt(delta);
	printf("final delta: %g\n", delta);

	if (k == maxiter)
	{
		fprintf(stderr, "Warning: 'OptimizeAngles' did not converge within %i iterations, final delta: %g.\n", maxiter, delta);
		err.codes |= E_CONVERGENCE;
		err.conv_delta = delta;
	}

	// copy values back
	for (i = 1; i < nelec; i++) {
		phi[i] = gsl_vector_get(s->x, i-1);
	}

	gsl_multimin_fdfminimizer_free(s);
	gsl_vector_free(x);

	// consistency check: all eigenvalues of Hessian matrix should be non-negative
	{
		double *hessfull = (double *)malloc(nelec*nelec*sizeof(double));
		double *hesssubm = (double *)malloc((nelec-1)*(nelec-1)*sizeof(double));
		double *eigenval = (double *)malloc((nelec-1)*sizeof(double));

		// calculate Hessian matrix
		VeeHessPhi(nelec, r, phi, hessfull);
		// omit first row and column in Hessian matrix since phi[0] is fixed
		for (j = 1; j < nelec; j++)
		{
			for (i = 1; i < nelec; i++)
			{
				hesssubm[(i-1)+(j-1)*(nelec-1)] = hessfull[i+j*nelec];
			}
		}
		// calculate eigenvalues
		dsyev('N', 'U', nelec-1, hesssubm, nelec-1, eigenval, &status);
		// check for errors
		if (status != 0) {
			fprintf(stderr, "Warning: 'dsyev' failed, info: %i\n", status);
			err.codes |= E_LAPACK;
		}
		if (eigenval[0] < 0) {
			fprintf(stderr, "Warning: Hessian has negative eigenvalue %g.\n", eigenval[0]);
			err.codes |= E_NEG_HESS;
			err.hess_eig = eigenval[0];
		}

		free(eigenval);
		free(hesssubm);
		free(hessfull);
	}

	return err;
}

//_______________________________________________________________________________________________________________________
///
/// \brief Perform a global optimization of the polar angles with respect to Coulomb repulsion, given the radii
///
/// \param nelec number of electrons
/// \param r radii, array of length 'nelec'
/// \param phi serves as input (starting values) and output (optimized angles); array of length 'nelec';
///        phi[0] is not changed due to invariance under global rotation
/// \param optVee optimized Coulomb repulsion (output)
/// \return detailed error description, if any
///
error_desc_t OptimizeAnglesGlobal(const int nelec, const double *r, double *phi, double *optVee)
{
	int i, k, l;

	error_desc_t err = OptimizeAngles(nelec, r, phi);

	// temporary storage for angles
	double *chi	= (double *)malloc(nelec*sizeof(double));

	// try permuting angles to find minimum value
	(*optVee) = Vee(nelec, r, phi);
	for (i = 0; i < 32; i++)		// iterate several times
	{
		// for all pairs of electrons...
		for (l = 1; l < nelec; l++)		// angle of first electron always zero
		{
			for (k = 1; k < l; k++)
			{
				double curV;

				memcpy(chi, phi, nelec*sizeof(double));
				// permute k <-> l
				chi[k] = phi[l];
				chi[l] = phi[k];
				error_desc_t cur_err = OptimizeAngles(nelec, r, chi);

				curV = Vee(nelec, r, chi);
				if (curV < (*optVee))
				{
					// only accumulate error if we really use the optimized angles
					AccumulateError(&err, &cur_err);

					// save values
					memcpy(phi, chi, nelec*sizeof(double));
					(*optVee) = curV;
				}
			}
		}
	}

	// clean up
	free(chi);

	return err;
}

/*
//_______________________________________________________________________________________________________________________
//
// GSL wrappers of Vee() for global minimization by simulated annealing

typedef struct
{
	int nelec;
	const double *r;
	double phi[10];
}
saConfig_t;

static double saVee(void *xp)
{
	const saConfig_t *x = (saConfig_t *)xp;

	return Vee(x->nelec, x->r, x->phi);
}

static void saStep(const gsl_rng *randgen, void *xp, double step_size)
{
	int i;
	saConfig_t *x = (saConfig_t *)xp;

	// first angle always fixed
	for (i = 1; i < x->nelec; i++)
	{
		double u = gsl_rng_uniform(randgen);
		x->phi[i] += step_size * (u - 0.5);		// divide by nelec ?

		// [-pi,pi] range
		if (x->phi[i] < -M_PI)
		{
			x->phi[i] += 2*M_PI;
		}
		else if (x->phi[i] > M_PI)
		{
			x->phi[i] -= 2*M_PI;
		}
	}
}

static double saDistance(void *xp, void *yp)
{
	int i;

	const saConfig_t *x = (saConfig_t *)xp;
	const saConfig_t *y = (saConfig_t *)yp;

	// distance between angles
	double dphi[10];
	for (i = 0; i < x->nelec; i++)
	{
		dphi[i] = x->phi[i] - y->phi[i];

		// [-pi,pi] range
		if (dphi[i] < -M_PI)
		{
			dphi[i] += 2*M_PI;
		}
		else if (dphi[i] > M_PI)
		{
			dphi[i] -= 2*M_PI;
		}
	}

	assert(x->nelec == y->nelec);

	return Norm(dphi + 1, x->nelec - 1);	// skip first angle (always held fixed)
}

static void saPrint(void *xp)
{
	int i;
	const saConfig_t *x = (saConfig_t *)xp;

	printf("   angles phi: ");
	for (i = 0; i < x->nelec; i++) {
		printf("%f ", x->phi[i]);
	}
}

//_______________________________________________________________________________________________________________________
///
/// \brief Globally optimize polar angles with respect to Coulomb repulsion, given the radii
///
/// \param nelec number of electrons
/// \param r radii, array of length 'nelec'
/// \param phi serves as input (starting values) and output (optimized angles); array of length 'nelec';
///        phi[0] is not changed due to invariance under global rotation
///
/// \return detailed error description, if any
///
error_desc_t OptimizeAnglesGlobal(const int nelec, const double *r, double *phi)
{
	const gsl_rng_type *T;
	gsl_rng *randgen;

	gsl_siman_params_t params;

	error_desc_t err = { 0 };

	// initial configuration
	saConfig_t config;
	config.nelec = nelec;
	config.r = r;
	memcpy(config.phi, phi, nelec*sizeof(double));

	gsl_rng_env_setup();

	// random number generator
	T = gsl_rng_default;
	randgen = gsl_rng_alloc(T);

	// parameters
	params.n_tries			= 200;				// how many points do we try before stepping
	params.iters_fixed_T	= 1000;				// how many iterations for each T?
	params.step_size		= 0.01 * M_PI;		// max step size in random walk
	params.k				= 1.0;				// virtual Boltzmann constant
	params.t_initial		= 0.001;			// initial temperature
	params.mu_t				= 1.004;			// damping factor for temperature
	params.t_min			= 2.0e-7;			// minimum temperature

	gsl_siman_solve(randgen, &config, saVee, saStep, saDistance, saPrint, NULL, NULL, NULL, sizeof(saConfig_t), params);

	// first angle must not change
	assert(config.phi[0] == phi[0]);

	// store optimized angles in 'phi'
	memcpy(phi, config.phi, nelec*sizeof(double));

	// clean up
	gsl_rng_free(randgen);

	// consistency check: all eigenvalues of Hessian matrix should be non-negative
	{
		int i, j;
		int status;
		double *hessfull = (double *)malloc(nelec*nelec*sizeof(double));
		double *hesssubm = (double *)malloc((nelec-1)*(nelec-1)*sizeof(double));
		double *eigenval = (double *)malloc((nelec-1)*sizeof(double));

		// calculate Hessian matrix
		VeeHessPhi(nelec, r, phi, hessfull);
		// omit first row and column in Hessian matrix since phi[0] is fixed
		for (j = 1; j < nelec; j++)
		{
			for (i = 1; i < nelec; i++)
			{
				hesssubm[(i-1)+(j-1)*(nelec-1)] = hessfull[i+j*nelec];
			}
		}
		// calculate eigenvalues
		dsyev('N', 'U', nelec-1, hesssubm, nelec-1, eigenval, &status);
		// check for errors
		if (status != 0) {
			fprintf(stderr, "Warning: 'dsyev' failed, info: %i\n", status);
			err.codes |= E_LAPACK;
		}
		if (eigenval[0] < 0) {
			fprintf(stderr, "Warning: Hessian has negative eigenvalue %g.\n", eigenval[0]);
			err.codes |= E_NEG_HESS;
			err.hess_eig = eigenval[0];
		}

		free(eigenval);
		free(hesssubm);
		free(hessfull);
	}

	return err;
}
*/

max_stars_repo_name: Arfaouim/quantum-dot

max_issues_repo_name: Arfaouim/quantum-dot

max_forks_repo_name: Arfaouim/quantum-dot


text:
#pragma once
#include <gsl/span>
#include "halley/core/api/audio_api.h"
#include "audio_buffer.h"

namespace Halley
{
	class AudioMixer
	{
	public:
		static void mixAudio(AudioSamplesConst src, AudioSamples dst, float gainStart, float gainEnd);
		static void mixAudio(AudioMultiChannelSamplesConst src, AudioMultiChannelSamples dst, float gainStart, float gainEnd);
		static void mixAudio(AudioMultiChannelSamples src, AudioMultiChannelSamples dst, float gainStart, float gainEnd);

		static void interleaveChannels(AudioSamples dst, gsl::span<AudioBuffer*> srcs);
		static void concatenateChannels(AudioSamples dst, gsl::span<AudioBuffer*> srcs);
		static void compressRange(AudioSamples buffer);

		static void zero(AudioSamples dst);
		static void zero(AudioMultiChannelSamples dst, size_t nChannels = 8);
		static void zeroRange(AudioMultiChannelSamples dst, size_t nChannels, size_t start, size_t len = std::numeric_limits<size_t>::max());
		static void copy(AudioMultiChannelSamples dst, AudioMultiChannelSamples src, size_t nChannels = 8);
		static void copy(AudioSamples dst, AudioSamples src);
		static void copy(AudioSamples dst, AudioSamples src, float gainStart, float gainEnd);
	};
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
/* blas/gsl_cblas.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* This is a copy of the CBLAS standard header.
 * We carry this around so we do not have to
 * break our model for flexible BLAS functionality.
 */

#ifndef __GSL_CBLAS_H__
#define __GSL_CBLAS_H__
#include <stddef.h>
#include <gsl/gsl_types.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

/*
 * Enumerated and derived types
 */
#define CBLAS_INDEX size_t  /* this may vary between platforms */

enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
enum CBLAS_UPLO {CblasUpper=121, CblasLower=122};
enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132};
enum CBLAS_SIDE {CblasLeft=141, CblasRight=142};

/*
 * ===========================================================================
 * Prototypes for level 1 BLAS functions (complex are recast as routines)
 * ===========================================================================
 */
GSL_EXPORT float  cblas_sdsdot(const int N, const float alpha, const float *X,
                               const int incX, const float *Y, const int incY);
GSL_EXPORT double cblas_dsdot(const int N, const float *X, const int incX, const float *Y,
                              const int incY);
GSL_EXPORT float  cblas_sdot(const int N, const float  *X, const int incX,
                             const float  *Y, const int incY);
GSL_EXPORT double cblas_ddot(const int N, const double *X, const int incX,
                             const double *Y, const int incY);

/*
 * Functions having prefixes Z and C only
 */
GSL_EXPORT void   cblas_cdotu_sub(const int N, const void *X, const int incX,
                                  const void *Y, const int incY, void *dotu);
GSL_EXPORT void   cblas_cdotc_sub(const int N, const void *X, const int incX,
                                  const void *Y, const int incY, void *dotc);

GSL_EXPORT void   cblas_zdotu_sub(const int N, const void *X, const int incX,
                                  const void *Y, const int incY, void *dotu);
GSL_EXPORT void   cblas_zdotc_sub(const int N, const void *X, const int incX,
                                  const void *Y, const int incY, void *dotc);

/*
 * Functions having prefixes S D SC DZ
 */
GSL_EXPORT float  cblas_snrm2(const int N, const float *X, const int incX);
GSL_EXPORT float  cblas_sasum(const int N, const float *X, const int incX);

GSL_EXPORT double cblas_dnrm2(const int N, const double *X, const int incX);
GSL_EXPORT double cblas_dasum(const int N, const double *X, const int incX);

GSL_EXPORT float  cblas_scnrm2(const int N, const void *X, const int incX);
GSL_EXPORT float  cblas_scasum(const int N, const void *X, const int incX);

GSL_EXPORT double cblas_dznrm2(const int N, const void *X, const int incX);
GSL_EXPORT double cblas_dzasum(const int N, const void *X, const int incX);

/*
 * Functions having standard 4 prefixes (S D C Z)
 */
GSL_EXPORT CBLAS_INDEX cblas_isamax(const int N, const float  *X, const int incX);
GSL_EXPORT CBLAS_INDEX cblas_idamax(const int N, const double *X, const int incX);
GSL_EXPORT CBLAS_INDEX cblas_icamax(const int N, const void   *X, const int incX);
GSL_EXPORT CBLAS_INDEX cblas_izamax(const int N, const void   *X, const int incX);

/*
 * ===========================================================================
 * Prototypes for level 1 BLAS routines
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (s, d, c, z)
 */
GSL_EXPORT void cblas_sswap(const int N, float *X, const int incX,
                            float *Y, const int incY);
GSL_EXPORT void cblas_scopy(const int N, const float *X, const int incX,
                            float *Y, const int incY);
GSL_EXPORT void cblas_saxpy(const int N, const float alpha, const float *X,
                            const int incX, float *Y, const int incY);

GSL_EXPORT void cblas_dswap(const int N, double *X, const int incX,
                            double *Y, const int incY);
GSL_EXPORT void cblas_dcopy(const int N, const double *X, const int incX,
                            double *Y, const int incY);
GSL_EXPORT void cblas_daxpy(const int N, const double alpha, const double *X,
                            const int incX, double *Y, const int incY);

GSL_EXPORT void cblas_cswap(const int N, void *X, const int incX,
                            void *Y, const int incY);
GSL_EXPORT void cblas_ccopy(const int N, const void *X, const int incX,
                            void *Y, const int incY);
GSL_EXPORT void cblas_caxpy(const int N, const void *alpha, const void *X,
                            const int incX, void *Y, const int incY);

GSL_EXPORT void cblas_zswap(const int N, void *X, const int incX,
                            void *Y, const int incY);
GSL_EXPORT void cblas_zcopy(const int N, const void *X, const int incX,
                            void *Y, const int incY);
GSL_EXPORT void cblas_zaxpy(const int N, const void *alpha, const void *X,
                            const int incX, void *Y, const int incY);

/*
 * Routines with S and D prefix only
 */
GSL_EXPORT void cblas_srotg(float *a, float *b, float *c, float *s);
GSL_EXPORT void cblas_srotmg(float *d1, float *d2, float *b1, const float b2, float *P);
GSL_EXPORT void cblas_srot(const int N, float *X, const int incX,
                           float *Y, const int incY, const float c, const float s);
GSL_EXPORT void cblas_srotm(const int N, float *X, const int incX,
                           float *Y, const int incY, const float *P);

GSL_EXPORT void cblas_drotg(double *a, double *b, double *c, double *s);
GSL_EXPORT void cblas_drotmg(double *d1, double *d2, double *b1, const double b2, double *P);
GSL_EXPORT void cblas_drot(const int N, double *X, const int incX,
                           double *Y, const int incY, const double c, const double  s);
GSL_EXPORT void cblas_drotm(const int N, double *X, const int incX,
                            double *Y, const int incY, const double *P);

/*
 * Routines with S D C Z CS and ZD prefixes
 */
GSL_EXPORT void cblas_sscal(const int N, const float alpha, float *X, const int incX);
GSL_EXPORT void cblas_dscal(const int N, const double alpha, double *X, const int incX);
GSL_EXPORT void cblas_cscal(const int N, const void *alpha, void *X, const int incX);
GSL_EXPORT void cblas_zscal(const int N, const void *alpha, void *X, const int incX);
GSL_EXPORT void cblas_csscal(const int N, const float alpha, void *X, const int incX);
GSL_EXPORT void cblas_zdscal(const int N, const double alpha, void *X, const int incX);

/*
 * ===========================================================================
 * Prototypes for level 2 BLAS
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_EXPORT void cblas_sgemv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const float alpha, const float *A, const int lda,
                            const float *X, const int incX, const float beta,
                            float *Y, const int incY);
GSL_EXPORT void cblas_sgbmv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const int KL, const int KU, const float alpha,
                            const float *A, const int lda, const float *X,
                            const int incX, const float beta, float *Y, const int incY);
GSL_EXPORT void cblas_strmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const float *A, const int lda,
                            float *X, const int incX);
GSL_EXPORT void cblas_stbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const float *A, const int lda,
                            float *X, const int incX);
GSL_EXPORT void cblas_stpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const float *Ap, float *X, const int incX);
GSL_EXPORT void cblas_strsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const float *A, const int lda, float *X,
                            const int incX);
GSL_EXPORT void cblas_stbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const float *A, const int lda,
                            float *X, const int incX);
GSL_EXPORT void cblas_stpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const float *Ap, float *X, const int incX);

GSL_EXPORT void cblas_dgemv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const double alpha, const double *A, const int lda,
                            const double *X, const int incX, const double beta,
                            double *Y, const int incY);
GSL_EXPORT void cblas_dgbmv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const int KL, const int KU, const double alpha,
                            const double *A, const int lda, const double *X,
                            const int incX, const double beta, double *Y, const int incY);
GSL_EXPORT void cblas_dtrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const double *A, const int lda,
                            double *X, const int incX);
GSL_EXPORT void cblas_dtbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const double *A, const int lda,
                            double *X, const int incX);
GSL_EXPORT void cblas_dtpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const double *Ap, double *X, const int incX);
GSL_EXPORT void cblas_dtrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const double *A, const int lda, double *X,
                            const int incX);
GSL_EXPORT void cblas_dtbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const double *A, const int lda,
                            double *X, const int incX);
GSL_EXPORT void cblas_dtpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const double *Ap, double *X, const int incX);

GSL_EXPORT void cblas_cgemv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *X, const int incX, const void *beta,
                            void *Y, const int incY);
GSL_EXPORT void cblas_cgbmv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const int KL, const int KU, const void *alpha,
                            const void *A, const int lda, const void *X,
                            const int incX, const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_ctrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ctbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ctpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *Ap, void *X, const int incX);
GSL_EXPORT void cblas_ctrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *A, const int lda, void *X,
                            const int incX);
GSL_EXPORT void cblas_ctbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ctpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *Ap, void *X, const int incX);

GSL_EXPORT void cblas_zgemv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *X, const int incX, const void *beta,
                            void *Y, const int incY);
GSL_EXPORT void cblas_zgbmv(const enum CBLAS_ORDER order,
                            const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
                            const int KL, const int KU, const void *alpha,
                            const void *A, const int lda, const void *X,
                            const int incX, const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_ztrmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ztbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ztpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *Ap, void *X, const int incX);
GSL_EXPORT void cblas_ztrsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *A, const int lda, void *X,
                            const int incX);
GSL_EXPORT void cblas_ztbsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const int K, const void *A, const int lda,
                            void *X, const int incX);
GSL_EXPORT void cblas_ztpsv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,
                            const int N, const void *Ap, void *X, const int incX);

/*
 * Routines with S and D prefixes only
 */
GSL_EXPORT void cblas_ssymv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const float alpha, const float *A,
                            const int lda, const float *X, const int incX,
                            const float beta, float *Y, const int incY);
GSL_EXPORT void cblas_ssbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const int K, const float alpha, const float *A,
                            const int lda, const float *X, const int incX,
                            const float beta, float *Y, const int incY);
GSL_EXPORT void cblas_sspmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const float alpha, const float *Ap,
                            const float *X, const int incX,
                            const float beta, float *Y, const int incY);
GSL_EXPORT void cblas_sger(const enum CBLAS_ORDER order, const int M, const int N,
                           const float alpha, const float *X, const int incX,
                           const float *Y, const int incY, float *A, const int lda);
GSL_EXPORT void cblas_ssyr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const float alpha, const float *X,
                           const int incX, float *A, const int lda);
GSL_EXPORT void cblas_sspr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const float alpha, const float *X,
                           const int incX, float *Ap);
GSL_EXPORT void cblas_ssyr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const float alpha, const float *X,
                            const int incX, const float *Y, const int incY, float *A,
                            const int lda);
GSL_EXPORT void cblas_sspr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const float alpha, const float *X,
                            const int incX, const float *Y, const int incY, float *A);

GSL_EXPORT void cblas_dsymv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const double alpha, const double *A,
                            const int lda, const double *X, const int incX,
                            const double beta, double *Y, const int incY);
GSL_EXPORT void cblas_dsbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const int K, const double alpha, const double *A,
                            const int lda, const double *X, const int incX,
                            const double beta, double *Y, const int incY);
GSL_EXPORT void cblas_dspmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const double alpha, const double *Ap,
                            const double *X, const int incX,
                            const double beta, double *Y, const int incY);
GSL_EXPORT void cblas_dger(const enum CBLAS_ORDER order, const int M, const int N,
                           const double alpha, const double *X, const int incX,
                           const double *Y, const int incY, double *A, const int lda);
GSL_EXPORT void cblas_dsyr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const double alpha, const double *X,
                           const int incX, double *A, const int lda);
GSL_EXPORT void cblas_dspr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const double alpha, const double *X,
                           const int incX, double *Ap);
GSL_EXPORT void cblas_dsyr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const double alpha, const double *X,
                            const int incX, const double *Y, const int incY, double *A,
                            const int lda);
GSL_EXPORT void cblas_dspr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const double alpha, const double *X,
                            const int incX, const double *Y, const int incY, double *A);

/*
 * Routines with C and Z prefixes only
 */
GSL_EXPORT void cblas_chemv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const void *alpha, const void *A,
                            const int lda, const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_chbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const int K, const void *alpha, const void *A,
                            const int lda, const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_chpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const void *alpha, const void *Ap,
                            const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_cgeru(const enum CBLAS_ORDER order, const int M, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_cgerc(const enum CBLAS_ORDER order, const int M, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_cher(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const float alpha, const void *X, const int incX,
                           void *A, const int lda);
GSL_EXPORT void cblas_chpr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const float alpha, const void *X,
                           const int incX, void *A);
GSL_EXPORT void cblas_cher2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_chpr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *Ap);

GSL_EXPORT void cblas_zhemv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const void *alpha, const void *A,
                            const int lda, const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_zhbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const int K, const void *alpha, const void *A,
                            const int lda, const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_zhpmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                            const int N, const void *alpha, const void *Ap,
                            const void *X, const int incX,
                            const void *beta, void *Y, const int incY);
GSL_EXPORT void cblas_zgeru(const enum CBLAS_ORDER order, const int M, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_zgerc(const enum CBLAS_ORDER order, const int M, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_zher(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const double alpha, const void *X, const int incX,
                           void *A, const int lda);
GSL_EXPORT void cblas_zhpr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                           const int N, const double alpha, const void *X,
                           const int incX, void *A);
GSL_EXPORT void cblas_zher2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *A, const int lda);
GSL_EXPORT void cblas_zhpr2(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo, const int N,
                            const void *alpha, const void *X, const int incX,
                            const void *Y, const int incY, void *Ap);

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_EXPORT void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                            const int K, const float alpha, const float *A,
                            const int lda, const float *B, const int ldb,
                            const float beta, float *C, const int ldc);
GSL_EXPORT void cblas_ssymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const float alpha, const float *A, const int lda,
                            const float *B, const int ldb, const float beta,
                            float *C, const int ldc);
GSL_EXPORT void cblas_ssyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const float alpha, const float *A, const int lda,
                            const float beta, float *C, const int ldc);
GSL_EXPORT void cblas_ssyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const float alpha, const float *A, const int lda,
                             const float *B, const int ldb, const float beta,
                             float *C, const int ldc);
GSL_EXPORT void cblas_strmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const float alpha, const float *A, const int lda,
                            float *B, const int ldb);
GSL_EXPORT void cblas_strsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const float alpha, const float *A, const int lda,
                            float *B, const int ldb);

GSL_EXPORT void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                            const int K, const double alpha, const double *A,
                            const int lda, const double *B, const int ldb,
                            const double beta, double *C, const int ldc);
GSL_EXPORT void cblas_dsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const double alpha, const double *A, const int lda,
                            const double *B, const int ldb, const double beta,
                            double *C, const int ldc);
GSL_EXPORT void cblas_dsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const double alpha, const double *A, const int lda,
                            const double beta, double *C, const int ldc);
GSL_EXPORT void cblas_dsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const double alpha, const double *A, const int lda,
                             const double *B, const int ldb, const double beta,
                             double *C, const int ldc);
GSL_EXPORT void cblas_dtrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const double alpha, const double *A, const int lda,
                            double *B, const int ldb);
GSL_EXPORT void cblas_dtrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const double alpha, const double *A, const int lda,
                            double *B, const int ldb);

GSL_EXPORT void cblas_cgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                            const int K, const void *alpha, const void *A,
                            const int lda, const void *B, const int ldb,
                            const void *beta, void *C, const int ldc);
GSL_EXPORT void cblas_csymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *B, const int ldb, const void *beta,
                            void *C, const int ldc);
GSL_EXPORT void cblas_csyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const void *alpha, const void *A, const int lda,
                            const void *beta, void *C, const int ldc);
GSL_EXPORT void cblas_csyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const void *alpha, const void *A, const int lda,
                             const void *B, const int ldb, const void *beta,
                             void *C, const int ldc);
GSL_EXPORT void cblas_ctrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            void *B, const int ldb);
GSL_EXPORT void cblas_ctrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            void *B, const int ldb);

GSL_EXPORT void cblas_zgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                            const int K, const void *alpha, const void *A,
                            const int lda, const void *B, const int ldb,
                            const void *beta, void *C, const int ldc);
GSL_EXPORT void cblas_zsymm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *B, const int ldb, const void *beta,
                            void *C, const int ldc);
GSL_EXPORT void cblas_zsyrk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const void *alpha, const void *A, const int lda,
                            const void *beta, void *C, const int ldc);
GSL_EXPORT void cblas_zsyr2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const void *alpha, const void *A, const int lda,
                             const void *B, const int ldb, const void *beta,
                             void *C, const int ldc);
GSL_EXPORT void cblas_ztrmm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            void *B, const int ldb);
GSL_EXPORT void cblas_ztrsm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const enum CBLAS_TRANSPOSE TransA,
                            const enum CBLAS_DIAG Diag, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            void *B, const int ldb);

/*
 * Routines with prefixes C and Z only
 */
GSL_EXPORT void cblas_chemm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *B, const int ldb, const void *beta,
                            void *C, const int ldc);
GSL_EXPORT void cblas_cherk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const float alpha, const void *A, const int lda,
                            const float beta, void *C, const int ldc);
GSL_EXPORT void cblas_cher2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const void *alpha, const void *A, const int lda,
                             const void *B, const int ldb, const float beta,
                             void *C, const int ldc);

GSL_EXPORT void cblas_zhemm(const enum CBLAS_ORDER Order, const enum CBLAS_SIDE Side,
                            const enum CBLAS_UPLO Uplo, const int M, const int N,
                            const void *alpha, const void *A, const int lda,
                            const void *B, const int ldb, const void *beta,
                            void *C, const int ldc);
GSL_EXPORT void cblas_zherk(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                            const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                            const double alpha, const void *A, const int lda,
                            const double beta, void *C, const int ldc);
GSL_EXPORT void cblas_zher2k(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
                             const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
                             const void *alpha, const void *A, const int lda,
                             const void *B, const int ldb, const double beta,
                             void *C, const int ldc);

GSL_EXPORT void cblas_xerbla(int p, const char *rout, const char *form, ...);

__END_DECLS

#endif /* __GSL_CBLAS_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#include <pygsl/utils.h>
#include <pygsl/block_helpers.h>
#include <pygsl/error_helpers.h>
#include <gsl/gsl_qrng.h>
#include <pygsl/error_helpers.h>
#include "qrng_module_defines.h"
static PyObject *module = NULL;

static void			/* generic instance destruction */
generic_dealloc (PyObject *self)
{
  DEBUG_MESS(1, " *** generic_dealloc %p\n", (void *) self);
  PyObject_Del(self);
}

/*---------------------------------------------------------------------------*/
/* encapsulation for the PyGSL_qrng_type */
typedef struct {
     PyObject_HEAD
     gsl_qrng_type * qrng_type;
     const char * py_name; /* can not use the name here... niederreiter-*/
} PyGSL_qrng_type;

/*
PyObject*
qrng_type_getattr(PyGSL_qrng_type *self, char *name)
{
     if(strcmp(name, "__doc__") == 0){
	  return PyGSL_get_docobject_for_object(self->py_name);
     }else{
	  Py_INCREF(Py_None);
	  return Py_None;
     }
}
*/
static PyTypeObject PyGSL_qrng_type_pytype = {
  PyObject_HEAD_INIT(NULL)	/* fix up the type slot in initcrng */
  0,				/* ob_size */
  "PyGSL_qrng_type",     	/* tp_name */
  sizeof(PyGSL_qrng_type),       /* tp_basicsize */
  0,				/* tp_itemsize */

  /* standard methods */
  (destructor)  generic_dealloc,   /* tp_dealloc  ref-count==0  */
  (printfunc)   0,		   /* tp_print    "print x"     */
  (getattrfunc) 0/*qrng_type_getattr*/, /* tp_getattr  "x.attr"      */
  (setattrfunc) 0,		   /* tp_setattr  "x.attr=v"    */
  (cmpfunc)     0,		   /* tp_compare  "x > y"       */
  (reprfunc)    0,                 /* tp_repr     `x`, print x  */

  /* type categories */
  0,				/* tp_as_number   +,-,*,/,%,&,>>,pow...*/
  0,				/* tp_as_sequence +,[i],[i:j],len, ...*/
  0,				/* tp_as_mapping  [key], len, ...*/

  /* more methods */
  (hashfunc)     0,		/* tp_hash    "dict[x]" */
  (ternaryfunc)  0,             /* tp_call    "x()"     */
  (reprfunc)     0,             /* tp_str     "str(x)"  */
  (getattrofunc) 0,		/* tp_getattro */
  (setattrofunc) 0,		/* tp_setattro */
  0,				/* tp_as_buffer */
  0L,				/* tp_flags */
  NULL/* PyGSL_qrng_type_pytype_doc*/       /* tp_doc */
};

void
create_qrng_types(PyObject *module)
{
     PyGSL_qrng_type *a_rng = NULL;
     const gsl_qrng_type* thisRNGType,
	  *types[3] = {gsl_qrng_niederreiter_2, gsl_qrng_sobol, NULL};
     const char* gsl_qrng_names[3] = {"_qrng.niederreiter_2", "_qrng.sobol", NULL};

     PyObject* module_dict=PyModule_GetDict(module);
     PyObject* item = NULL;
     int i;

     assert(module_dict);
     FUNC_MESS_BEGIN();
     /* provide other rng types as subclasses of rng  */
     for(i=0; types[i]!=NULL; i++){
	  thisRNGType = types[i];
	  a_rng =  PyObject_NEW(PyGSL_qrng_type, &PyGSL_qrng_type_pytype);
	  a_rng->qrng_type = (gsl_qrng_type *) thisRNGType;
	  item = PyString_FromString((thisRNGType)->name);
	  assert(item);
	  PyGSL_clear_name(PyString_AsString(item), PyString_Size(item));
	  assert(gsl_qrng_names[i]);
	  a_rng->py_name = gsl_qrng_names[i];
	  PyDict_SetItem(module_dict, item, (PyObject *) a_rng);
	  Py_DECREF(item);
	  item = NULL;
	  thisRNGType++;
     }
     FUNC_MESS_END();
}
#define PyGSL_QRngType_Check(v) ((v)->ob_type == &PyGSL_qrng_type_pytype)

/*---------------------------------------------------------------------------*/
/* qrng */

staticforward PyTypeObject PyGSL_qrng_pytype;
#if 1
typedef struct {
  PyObject_HEAD
  gsl_qrng * qrng;
} PyGSL_qrng;
#endif
static void qrng_delete(PyGSL_qrng *self);

static PyObject *		/* on "instance.attr" */
qrng_getattr (PyGSL_qrng *self, char *name);

static PyObject *
qrng_get(PyGSL_qrng *self, PyObject *args);

static PyTypeObject PyGSL_qrng_pytype = {
  PyObject_HEAD_INIT(NULL)	/* fix up the type slot in initcrng */
  0,				/* ob_size */
  "PyGSL_qrng",			/* tp_name */
  sizeof(PyGSL_qrng),	        /* tp_basicsize */
  0,				/* tp_itemsize */

  /* standard methods */
  (destructor)  qrng_delete,       /* tp_dealloc  ref-count==0  */
  (printfunc)   0,		   /* tp_print    "print x"     */
  (getattrfunc) qrng_getattr,       /* tp_getattr  "x.attr"      */
  (setattrfunc) 0,		   /* tp_setattr  "x.attr=v"    */
  (cmpfunc)     0,		   /* tp_compare  "x > y"       */
  (reprfunc)    0,                 /* tp_repr     `x`, print x  */

  /* type categories */
  0,				/* tp_as_number   +,-,*,/,%,&,>>,pow...*/
  0,				/* tp_as_sequence +,[i],[i:j],len, ...*/
  0,				/* tp_as_mapping  [key], len, ...*/

  /* more methods */
  (hashfunc)     0,		/* tp_hash    "dict[x]" */
  (ternaryfunc)  qrng_get,      /* tp_call    "x()"     */
  (reprfunc)     0,             /* tp_str     "str(x)"  */
  (getattrofunc) 0,		/* tp_getattro */
  (setattrofunc) 0,		/* tp_setattro */
  0,				/* tp_as_buffer */
  0L,				/* tp_flags */
  0		/* tp_doc */
};

#define PyGSLQRng_Check(v) ((v)->ob_type == &PyGSL_qrng_pytype)

PyObject *
qrng_init(PyObject *self, PyObject *args)
{
     PyObject *type = NULL;
     PyGSL_qrng *qrng = NULL;
     int dimension;
     assert(args);
     if (0 == PyArg_ParseTuple(args, "O!i:rng.__init__", &PyGSL_qrng_type_pytype,
			       &type, &dimension)){
	  PyGSL_add_traceback(module, __FILE__, "rng.__init__", __LINE__ - 2);
	  return NULL;
     }
     if (dimension <= 0){
	  PyErr_SetString(PyExc_ValueError, "The sample number must be positive!");
	  PyGSL_add_traceback(module, __FILE__, "qrng.__init__", __LINE__ - 2);
	  return NULL;
     }
     qrng =  (PyGSL_qrng *) PyObject_NEW(PyGSL_qrng, &PyGSL_qrng_pytype);
     qrng->qrng = gsl_qrng_alloc(((PyGSL_qrng_type *)type)->qrng_type, dimension);
     return (PyObject *) qrng;
}

static void
qrng_delete(PyGSL_qrng *self)
{
     assert(PyGSLQRng_Check(self));
     gsl_qrng_free(self->qrng);
     DEBUG_MESS(1, " self %p\n",(void *) self);
}

static PyObject *
qrng_reinit(PyGSL_qrng *self, PyObject *args)
{
     assert(PyGSLQRng_Check(self));
     gsl_qrng_init(self->qrng);
     Py_INCREF(Py_None);
     return Py_None;
}

static PyObject *
qrng_get(PyGSL_qrng *self, PyObject *args)
{
     int dimension = 1, lineno=0;
     PyArrayObject *a_array;
     PyGSL_array_index_t dims[2];
     double *data;
     int i;

     FUNC_MESS_BEGIN();
     assert(PyGSLQRng_Check(self));
     if(0 == PyArg_ParseTuple(args, "|i", &dimension)){
	  goto fail;
     }
     if(dimension <= 0){
	  lineno = __LINE__ - 1;
	  PyErr_SetString(PyExc_ValueError,
			  "The sample number must be positive!");
	  goto fail;
     }
     dims[0] = dimension;
     dims[1] = self->qrng->dimension;
     DEBUG_MESS(5, "Building return array with dimensions (%ld,%ld)", (long)dims[0], (long)dims[1]);
     a_array = (PyArrayObject *) PyGSL_New_Array(2, dims, PyArray_DOUBLE);
     if(a_array == NULL){lineno = __LINE__ - 1; goto fail;}
     DEBUG_MESS(5, "Its strides are (%d,%d)", a_array->strides[0], a_array->strides[1]);
     assert((a_array->strides[1] / sizeof(double)) == 1);

     for(i=0; i<dimension; i++){
	  DEBUG_MESS(6, "Setting slice %d", i);
	  data = (double *) (a_array->data + a_array->strides[0] * i);
	  DEBUG_MESS(6, "Data at %p", (void *) data);
	  gsl_qrng_get(self->qrng, data);
     }

     FUNC_MESS_END();
     return (PyObject *) a_array;

 fail:
     FUNC_MESS("In Fail!");
     PyGSL_add_traceback(module, __FILE__, "_qrng.__attr__",  lineno);
     return NULL;
}

static PyObject *
qrng_name(PyGSL_qrng *self, PyObject *args)
{
     assert(PyGSLQRng_Check(self));
     if(0 == PyArg_ParseTuple(args, ":name"))
	  return NULL;
     return PyString_FromString(gsl_qrng_name(self->qrng));
}

static PyObject *
qrng_clone(PyGSL_qrng *self, PyObject *args)
{
     PyGSL_qrng * qrng;
     assert(PyGSLQRng_Check(self));
     if(0 == PyArg_ParseTuple(args, ":clone"))
	  return NULL;
     qrng =  (PyGSL_qrng *) PyObject_NEW(PyGSL_qrng, &PyGSL_qrng_pytype);
     qrng->qrng = gsl_qrng_clone(self->qrng);
     return (PyObject *) qrng;
}

static struct PyMethodDef qrng_methods[] = {
     {"init", (PyCFunction) qrng_reinit, METH_VARARGS, QRNG_INIT},
     {"get", (PyCFunction) qrng_get, METH_VARARGS, QRNG_GET},
     {"name", (PyCFunction) qrng_name, METH_VARARGS, QRNG_NAME},
     {"clone", (PyCFunction) qrng_clone, METH_VARARGS, QRNG_CLONE},
     {"__copy__", (PyCFunction) qrng_clone, METH_VARARGS, QRNG_CLONE},
     {NULL, NULL}
};

static PyObject *		/* on "instance.attr" */
qrng_getattr (PyGSL_qrng *self, char *name){
     PyObject *tmp = NULL;
     assert(PyGSLQRng_Check(self));

     /*
     if(strcmp(name, "__doc__") == 0)
	  return PyGSL_get_docobject_for_object("qrng.qrng");
     */
     tmp = Py_FindMethod(qrng_methods, (PyObject *) self, name);
     if(NULL == tmp){
	  PyGSL_add_traceback(module, __FILE__, "qrng.__attr__", __LINE__ - 1);
	  return NULL;
     }
     return tmp;

}

static PyMethodDef PyGSL_qrng_module_functions[] = {
     {"qrng", qrng_init, METH_VARARGS},
     {0, 0}
};

void
init_qrng(void)
{
     PyObject *m;

     m = Py_InitModule("_qrng", PyGSL_qrng_module_functions);
     init_pygsl();

     assert(m);
     create_qrng_types(m);
     module = m;

     PyGSL_qrng_type_pytype.ob_type = &PyType_Type;
     PyGSL_qrng_pytype.ob_type = &PyType_Type;

     /*
       {
       PyObject *o, *s, *d;
       d = PyModule_GetDict(m);

       s = PyGSL_get_docobject_for_object("_qrng");
       if(s)
          PyDict_SetItemString(d, "__doc__", s);

     o = PyGSL_get_docobject_for_object("_qrng.PyGSL_qrng_type");
     if(o){
	  Py_INCREF(o);
	  PyGSL_qrng_type_pytype.tp_doc = PyString_AsString(o);
     }
     o = PyGSL_get_docobject_for_object("_qrng.PyGSL_qrng");
     if(o){
	  Py_INCREF(o);
	  PyGSL_qrng_type_pytype.tp_doc = PyString_AsString(o);
     }
     }
     */

}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod