text:
// Standard libraries.
#include <math.h>
#include <signal.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// Libraries from packages outside ASF.
#include <glib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_statistics_double.h>

// Libraries developed at ASF.
#include <asf.h>
#include <asf_meta.h>
#include <asf_raster.h>
#include "float_image.h"
#include "uint8_image.h"
#include <libasf_proj.h>
#include <spheroids.h>
#include <asf_contact.h>

void usage()
{
    printf("Usage: flip [v|h|vh] <input basename> <output basename>\n");
    exit(1);
}

int main(int argc, char *argv[])
{
  handle_common_asf_args(&argc, &argv, "flip");

  if (argc != 4) {
    usage();
  }

  if ( ! (strcmp(argv[1], "v") == 0 ||
	  strcmp(argv[1], "h") == 0 ||
	  strcmp(argv[1], "hv") == 0 ||
	  strcmp(argv[1], "vh") == 0))
  {
    usage();
  }

  int vert = strchr(argv[1], 'v') != NULL;
  int horz = strchr(argv[1], 'h') != NULL;

  if (!vert && !horz) {
    usage();
  }

  char *input_meta_name = appendExt(argv[2], ".meta");
  char *input_data_name = appendExt(argv[2], ".img");

  char *output_meta_name = appendExt(argv[3], ".meta");
  char *output_data_name = appendExt(argv[3], ".img");

  printf("Flipping image %s.\n",
	 vert && horz ? "vertically and horizontally" :
	 (vert ? "vertically" : "horizontally"));

  printf("Input data file: %s\n", input_data_name);
  printf("Output data file: %s\n", output_data_name);

  meta_parameters *imd = meta_read(input_meta_name);
  meta_write(imd, output_meta_name);

  FloatImage *finput = NULL;
  UInt8Image *binput = NULL;
  if (imd->optical || imd->general->data_type == BYTE) {
    binput = uint8_image_new_from_file(imd->general->sample_count,
                                       imd->general->line_count,
                                       input_data_name, 0);
  }
  else {
    finput = float_image_new_from_file(imd->general->sample_count,
                                       imd->general->line_count,
                                       input_data_name, 0,
                                       FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);
  }

  if (imd->optical || imd->general->data_type == BYTE) {
    if (horz)
      uint8_image_flip_x(binput);

    if (vert)
      uint8_image_flip_y(binput);

    uint8_image_store(binput, output_data_name);
    uint8_image_free(binput);
  }
  else {
    if (horz)
      float_image_flip_x(finput);

    if (vert)
      float_image_flip_y(finput);

    float_image_store(finput, output_data_name,
                      FLOAT_IMAGE_BYTE_ORDER_BIG_ENDIAN);
    float_image_free(finput);
  }

  meta_free(imd);

  free(input_data_name);
  free(input_meta_name);
  free(output_data_name);
  free(output_meta_name);

  return EXIT_SUCCESS;
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/* Cholesky Decomposition
 *
 * Copyright (C) 2000 Thomas Walter
 * Copyright (C) 2000, 2001, 2002, 2003, 2005, 2007 Brian Gough, Gerard Jungman
 * Copyright (C) 2016, 2019 Patrick Alken
 *
 * This 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, or (at your option) any
 * later version.
 *
 * This source 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.
 *
 * 03 May 2000: Modified for GSL by Brian Gough
 * 29 Jul 2005: Additions by Gerard Jungman
 * 04 Mar 2016: Change Cholesky algorithm to gaxpy version by Patrick Alken
 * 26 May 2019: implement recursive Cholesky with Level 3 BLAS by Patrick Alken
 */

/*
 * Cholesky decomposition of a symmetric positive definite matrix.
 *
 * This algorithm does:
 *   A = L * L'
 * with
 *   L  := lower left triangle matrix
 *   L' := the transposed form of L.
 *
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

#include "recurse.h"

static double cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work);
static int cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
static int cholesky_decomp_L2 (gsl_matrix * A);
static int cholesky_decomp_L3 (gsl_matrix * A);

/*
In GSL 2.2, we decided to modify the behavior of the Cholesky decomposition
to store the Cholesky factor in the lower triangle, and store the original
matrix in the upper triangle. Previous versions stored the Cholesky factor in
both places. The routine gsl_linalg_cholesky_decomp1 was added for the new
behavior, and gsl_linalg_cholesky_decomp is maintained for backward compatibility.
It will be removed in a future release.
*/

int
gsl_linalg_cholesky_decomp (gsl_matrix * A)
{
  int status;

  status = gsl_linalg_cholesky_decomp1(A);
  if (status == GSL_SUCCESS)
    {
      gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, A, A);
    }

  return status;
}

/*
gsl_linalg_cholesky_decomp1()
  Perform Cholesky decomposition of a symmetric positive
definite matrix using lower triangle using Level 3 BLAS algorithm.

Inputs: A - (input) symmetric, positive definite matrix
            (output) lower triangle contains Cholesky factor

Return: success/error

Notes:
1) original matrix is saved in upper triangle on output
*/

int
gsl_linalg_cholesky_decomp1 (gsl_matrix * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      /* save original matrix in upper triangle for later rcond calculation */
      gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, A, A);

      return cholesky_decomp_L3(A);
    }
}

int
gsl_linalg_cholesky_solve (const gsl_matrix * LLT,
                           const gsl_vector * b,
                           gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* copy x <- b */
      gsl_vector_memcpy (x, b);

      status = gsl_linalg_cholesky_svx(LLT, x);

      return status;
    }
}

int
gsl_linalg_cholesky_svx (const gsl_matrix * LLT,
                         gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for c using forward-substitution, L c = b */
      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* perform back-substitution, L^T x = c */
      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_solve_mat (const gsl_matrix * LLT,
                               const gsl_matrix * B,
                               gsl_matrix * X)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != B->size1)
    {
      GSL_ERROR ("matrix size must match B size", GSL_EBADLEN);
    }
  else if (LLT->size2 != X->size1)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* copy X <- B */
      gsl_matrix_memcpy (X, B);

      status = gsl_linalg_cholesky_svx_mat(LLT, X);

      return status;
    }
}

int
gsl_linalg_cholesky_svx_mat (const gsl_matrix * LLT,
                             gsl_matrix * X)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size2 != X->size1)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for C using forward-substitution, L C = B */
      gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, 1.0,
                      LLT, X);

      /* perform back-substitution, L^T X = C */
      gsl_blas_dtrsm (CblasLeft, CblasLower, CblasTrans, CblasNonUnit, 1.0,
                      LLT, X);

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_cholesky_invert()
  Compute the inverse of a symmetric positive definite matrix in
Cholesky form.

Inputs: LLT - matrix in cholesky form on input
              A^{-1} = L^{-t} L^{-1} on output

Return: success or error
*/

int
gsl_linalg_cholesky_invert(gsl_matrix * LLT)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else
    {
      int status;

      /* invert the lower triangle of LLT */
      status = gsl_linalg_tri_invert(CblasLower, CblasNonUnit, LLT);
      if (status)
        return status;

      /* compute A^{-1} = L^{-T} L^{-1} */
      status = gsl_linalg_tri_LTL(LLT);
      if (status)
        return status;

      /* copy lower triangle to upper */
      gsl_matrix_transpose_tricpy(CblasLower, CblasUnit, LLT, LLT);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_decomp_unit(gsl_matrix * A, gsl_vector * D)
{
  const size_t N = A->size1;
  size_t i, j;

  /* initial Cholesky */
  int stat_chol = gsl_linalg_cholesky_decomp1(A);

  if(stat_chol == GSL_SUCCESS)
  {
    /* calculate D from diagonal part of initial Cholesky */
    for(i = 0; i < N; ++i)
    {
      const double C_ii = gsl_matrix_get(A, i, i);
      gsl_vector_set(D, i, C_ii*C_ii);
    }

    /* multiply initial Cholesky by 1/sqrt(D) on the right */
    for(i = 0; i < N; ++i)
    {
      for(j = 0; j < N; ++j)
      {
        gsl_matrix_set(A, i, j, gsl_matrix_get(A, i, j) / sqrt(gsl_vector_get(D, j)));
      }
    }

    /* Because the initial Cholesky contained both L and transpose(L),
       the result of the multiplication is not symmetric anymore;
       but the lower triangle _is_ correct. Therefore we reflect
       it to the upper triangle and declare victory.
       */
    for(i = 0; i < N; ++i)
      for(j = i + 1; j < N; ++j)
        gsl_matrix_set(A, i, j, gsl_matrix_get(A, j, i));
  }

  return stat_chol;
}

/*
gsl_linalg_cholesky_scale()
  This function computes scale factors diag(S), such that

diag(S) A diag(S)

has a condition number within a factor N of the matrix
with the smallest condition number over all possible
diagonal scalings. See Corollary 7.6 of:

N. J. Higham, Accuracy and Stability of Numerical Algorithms (2nd Edition),
SIAM, 2002.

Inputs: A - symmetric positive definite matrix
        S - (output) scale factors, S_i = 1 / sqrt(A_ii)
*/

int
gsl_linalg_cholesky_scale(const gsl_matrix * A, gsl_vector * S)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
    }
  else if (N != S->size)
    {
      GSL_ERROR("S must have length N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute S_i = 1/sqrt(A_{ii}) */
      for (i = 0; i < N; ++i)
        {
          double Aii = gsl_matrix_get(A, i, i);

          if (Aii <= 0.0)
            gsl_vector_set(S, i, 1.0); /* matrix not positive definite */
          else
            gsl_vector_set(S, i, 1.0 / sqrt(Aii));
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_cholesky_scale_apply()
  This function applies scale transformation to A:

A <- diag(S) A diag(S)

Inputs: A     - (input/output)
                on input, symmetric positive definite matrix
                on output, diag(S) * A * diag(S) in lower triangle
        S     - (input) scale factors
*/

int
gsl_linalg_cholesky_scale_apply(gsl_matrix * A, const gsl_vector * S)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("A is not a square matrix", GSL_ENOTSQR);
    }
  else if (N != S->size)
    {
      GSL_ERROR("S must have length N", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      /* compute: A <- diag(S) A diag(S) using lower triangle */
      for (j = 0; j < N; ++j)
        {
          double sj = gsl_vector_get(S, j);

          for (i = j; i < N; ++i)
            {
              double si = gsl_vector_get(S, i);
              double *Aij = gsl_matrix_ptr(A, i, j);
              *Aij *= si * sj;
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_decomp2(gsl_matrix * A, gsl_vector * S)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR("cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (N != S->size)
    {
      GSL_ERROR("S must have length N", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* compute scaling factors to reduce cond(A) */
      status = gsl_linalg_cholesky_scale(A, S);
      if (status)
        return status;

      /* apply scaling factors */
      status = gsl_linalg_cholesky_scale_apply(A, S);
      if (status)
        return status;

      /* compute Cholesky decomposition of diag(S) A diag(S) */
      status = gsl_linalg_cholesky_decomp1(A);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_svx2 (const gsl_matrix * LLT,
                          const gsl_vector * S,
                          gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size2 != S->size)
    {
      GSL_ERROR ("matrix size must match S", GSL_EBADLEN);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* b~ = diag(S) b */
      gsl_vector_mul(x, S);

      /* Solve for c using forward-substitution, L c = b~ */
      gsl_blas_dtrsv (CblasLower, CblasNoTrans, CblasNonUnit, LLT, x);

      /* Perform back-substitution, L^T x~ = c */
      gsl_blas_dtrsv (CblasLower, CblasTrans, CblasNonUnit, LLT, x);

      /* compute original solution vector x = S x~ */
      gsl_vector_mul(x, S);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_solve2 (const gsl_matrix * LLT,
                            const gsl_vector * S,
                            const gsl_vector * b,
                            gsl_vector * x)
{
  if (LLT->size1 != LLT->size2)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != S->size)
    {
      GSL_ERROR ("matrix size must match S size", GSL_EBADLEN);
    }
  else if (LLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LLT->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* Copy x <- b */
      gsl_vector_memcpy (x, b);

      status = gsl_linalg_cholesky_svx2(LLT, S, x);

      return status;
    }
}

int
gsl_linalg_cholesky_rcond (const gsl_matrix * LLT, double * rcond,
                           gsl_vector * work)
{
  const size_t M = LLT->size1;
  const size_t N = LLT->size2;

  if (M != N)
    {
      GSL_ERROR ("cholesky matrix must be square", GSL_ENOTSQR);
    }
  else if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      int status;
      double Anorm = cholesky_norm1(LLT, work); /* ||A||_1 */
      double Ainvnorm;                          /* ||A^{-1}||_1 */

      *rcond = 0.0;

      /* don't continue if matrix is singular */
      if (Anorm == 0.0)
        return GSL_SUCCESS;

      /* estimate ||A^{-1}||_1 */
      status = gsl_linalg_invnorm1(N, cholesky_Ainv, (void *) LLT, &Ainvnorm, work);

      if (status)
        return status;

      if (Ainvnorm != 0.0)
        *rcond = (1.0 / Anorm) / Ainvnorm;

      return GSL_SUCCESS;
    }
}

/* compute 1-norm of original matrix, stored in upper triangle of LLT;
 * diagonal entries have to be reconstructed */
static double
cholesky_norm1(const gsl_matrix * LLT, gsl_vector * work)
{
  const size_t N = LLT->size1;
  double max = 0.0;
  size_t i, j;

  for (j = 0; j < N; ++j)
    {
      double sum = 0.0;
      gsl_vector_const_view lj = gsl_matrix_const_subrow(LLT, j, 0, j + 1);
      double Ajj;

      /* compute diagonal (j,j) entry of A */
      gsl_blas_ddot(&lj.vector, &lj.vector, &Ajj);

      for (i = 0; i < j; ++i)
        {
          double *wi = gsl_vector_ptr(work, i);
          double Aij = gsl_matrix_get(LLT, i, j);
          double absAij = fabs(Aij);

          sum += absAij;
          *wi += absAij;
        }

      gsl_vector_set(work, j, sum + fabs(Ajj));
    }

  for (i = 0; i < N; ++i)
    {
      double wi = gsl_vector_get(work, i);
      max = GSL_MAX(max, wi);
    }

  return max;
}

/* x := A^{-1} x = A^{-t} x, A = L L^T */
static int
cholesky_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  int status;
  gsl_matrix * A = (gsl_matrix * ) params;

  (void) TransA; /* unused parameter warning */

  /* compute L^{-1} x */
  status = gsl_blas_dtrsv(CblasLower, CblasNoTrans, CblasNonUnit, A, x);
  if (status)
    return status;

  /* compute L^{-t} x */
  status = gsl_blas_dtrsv(CblasLower, CblasTrans, CblasNonUnit, A, x);
  if (status)
    return status;

  return GSL_SUCCESS;
}

/*
cholesky_decomp_L2()
  Perform Cholesky decomposition of a symmetric positive
definite matrix using lower triangle

Inputs: A - (input) symmetric, positive definite matrix
            (output) lower triangle contains Cholesky factor

Return: success/error

Notes:
1) Based on algorithm 4.2.1 (Gaxpy Cholesky) of Golub and
Van Loan, Matrix Computations (4th ed), using Level 2 BLAS.
*/

static int
cholesky_decomp_L2 (gsl_matrix * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else
    {
      size_t j;

      for (j = 0; j < N; ++j)
        {
          double ajj;
          gsl_vector_view v = gsl_matrix_subcolumn(A, j, j, N - j); /* A(j:n,j) */

          if (j > 0)
            {
              gsl_vector_view w = gsl_matrix_subrow(A, j, 0, j);           /* A(j,1:j-1)^T */
              gsl_matrix_view m = gsl_matrix_submatrix(A, j, 0, N - j, j); /* A(j:n,1:j-1) */

              gsl_blas_dgemv(CblasNoTrans, -1.0, &m.matrix, &w.vector, 1.0, &v.vector);
            }

          ajj = gsl_matrix_get(A, j, j);

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          gsl_vector_scale(&v.vector, 1.0 / ajj);
        }

      return GSL_SUCCESS;
    }
}

/*
cholesky_decomp_L3()
  Perform Cholesky decomposition of a symmetric positive
definite matrix using Level 3 BLAS.

Inputs: A - (input) symmetric, positive definite matrix in lower triangle
            (output) lower triangle contains Cholesky factor

Return: success/error

Notes:
1) Based on ReLAPACK recursive block Cholesky algorithm using Level 3 BLAS

2) 28 May 2019: performed several benchmark tests of this recursive variant
against the right-looking block variant from LAPACK. This recursive variant
performed faster in all cases, so it is now the default algorithm.
*/

static int
cholesky_decomp_L3 (gsl_matrix * A)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("Cholesky decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (N <= CROSSOVER_CHOLESKY)
    {
      /* use unblocked Level 2 algorithm */
      return cholesky_decomp_L2(A);
    }
  else
    {
      /*
       * partition matrix:
       *
       * A11 A12
       * A21 A22
       *
       * where A11 is N1-by-N1
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT(N);
      const size_t N2 = N - N1;
      gsl_matrix_view A11 = gsl_matrix_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_view A21 = gsl_matrix_submatrix(A, N1, 0, N2, N1);
      gsl_matrix_view A22 = gsl_matrix_submatrix(A, N1, N1, N2, N2);

      /* recursion on A11 */
      status = cholesky_decomp_L3(&A11.matrix);
      if (status)
        return status;

      /* A21 = A21 * L11^{-T} */
      gsl_blas_dtrsm(CblasRight, CblasLower, CblasTrans, CblasNonUnit, 1.0, &A11.matrix, &A21.matrix);

      /* A22 -= L21 L21^T */
      gsl_blas_dsyrk(CblasLower, CblasNoTrans, -1.0, &A21.matrix, 1.0, &A22.matrix);

      /* recursion on A22 */
      status = cholesky_decomp_L3(&A22.matrix);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include "ll.h"
#include "gslsupp.h"

#include <math.h>
#include <gsl/gsl_sf.h>

double log_fact(size_t n)
{
    gsl_sf_result res;
    if (n <= UINT_MAX)
    {
        if (gsl_sf_lnfact_e((unsigned) n, &res) == GSL_SUCCESS) return res.val;
    }
    else
    {
        if (gsl_sf_lngamma_e((double) n + 1., &res) == GSL_SUCCESS) return res.val;
    }
    return nan(__func__);
}

double log_choose(size_t n, size_t m)
{
    if (m > n) return nan(__func__);
    if (m == n || !m) return 0.;
    return m > (n >> 1) ? log_fact(n) - log_fact(n - m) - log_fact(m) : log_fact(n) - log_fact(m) - log_fact(n - m); // Addition here may be not associative
}

double pdf_hypergeom(size_t k, size_t n1, size_t n2, size_t t)
{
    size_t car, n12 = size_add(&car, n1, n2);
    if (car) return nan(__func__);
    if (t > n12) t = n12;
    if (k > n1 || k > t || (t > n2 && k < t - n2)) return 0.;
    return exp(log_choose(n1, k) + log_choose(n2, t - k) - log_choose(n12, t));
}

double gamma_inc_P(double a, double x)
{
    gsl_sf_result res;
    if (gsl_sf_gamma_inc_P_e(a, x, &res) == GSL_SUCCESS) return res.val;
    return nan(__func__);
}

double gamma_inc_Q(double a, double x)
{
    gsl_sf_result res;
    if (gsl_sf_gamma_inc_Q_e(a, x, &res) == GSL_SUCCESS) return res.val;
    return nan(__func__);
}

double cdf_gamma_Q(double x, double a, double b)
{
    if (x <= 0.) return 1.;
    double y = x / b;
    return y < a ? 1. - gamma_inc_P(a, y) : gamma_inc_Q(a, y);
}

double cdf_chisq_Q(double x, double df)
{
    return cdf_gamma_Q(x, df / 2., 2.);
}


text:
#ifndef KGS_GSL_HELPERS_H
#define KGS_GSL_HELPERS_H

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

/** Print matrix to cout */
void gsl_matrix_cout (      const gsl_matrix *m);

/** Print matrix to file specified by `filename` */
void gsl_matrix_outtofile ( const gsl_matrix *m, const std::string& filename    );

/** Print vector to file specified by `filename` */
void gsl_vector_outtofile ( const gsl_vector *v, const std::string& filename    );

/** Print vector to output-stream */
void gsl_vector_out (const gsl_vector *v, std::ostream& os);

/** Print vector to cout */
void gsl_vector_cout (const gsl_vector *v);

/** Return length of vector */
double gsl_vector_length(   const gsl_vector *v );

/** Scale vector to length 1.0 */
void gsl_vector_normalize ( gsl_vector *v );

/** Scale vector to specified length */
void gsl_vector_scale_to_length(gsl_vector* ret, double length);

/** Scale vector so largest absolute value of any component is at most `maxComponent`.
 * If the values of all components are less than `maxComponent` no scaling is performed. */
void gsl_vector_scale_max_component(gsl_vector* v, double maxComponent);

/** Make a copy of the vector */
gsl_vector* gsl_vector_copy(gsl_vector*);

gsl_matrix* gsl_matrix_trans(gsl_matrix* A);

gsl_matrix* gsl_matrix_mul(gsl_matrix* A, gsl_matrix* B);

gsl_vector* gsl_matrix_vector_mul(gsl_matrix* A, gsl_vector* v);

gsl_matrix* pca (gsl_matrix* sample_matrix);

gsl_vector* RandomUnitVector (int size);

double frobenius_norm (const gsl_matrix *m);

/** Compute Shannon entropy (as in JCIM paper: fraction of significant contributors) for given input vector**/
double fractionOfSignificantContributors(const gsl_vector *v);

/** Compute Shannon entropy unnormalized: number of significant contributors of given input vector**/
double significantContributors(const gsl_vector *v);

/** Compute Shannon entropy from information theory in bits for given input vector**/
double shannonEntropyInBits(const gsl_vector* v);

/** Compute Shannon entropy from information theory in bits for given input vector**/
double shannonEntropyUnnormalizedInBits(const gsl_vector* v);

#endif //KGS_GSL_HELPERS_H

ext: h


text:
/* tz_fcmp.c
 *
 * 19-Nov-2007 Initial write: Ting Zhao
 */

#include <gsl/gsl_math.h>
#include "tz_mexutils.h"

void check_argin(int nrhs, const mxArray *prhs[])
{
  if (nrhs > 3  || nrhs < 2)
    mxErrMsgTxt("TZ_FCMP takes 2~3 arguments.");

  if (!(mxIsDouble(prhs[0]) && mxIsDouble(prhs[1])))
    mxErrMsgTxt("The first two arguments must be both double arrays.");

  if (!tz_mxSameSize(prhs[0], prhs[1]) && !tz_mxIsScalar(prhs[0]) &&
      !tz_mxIsScalar(prhs[1]))
    mxErrMsgTxt("The first two arguments must have the same size if neither of them is a scalar.");

  if (nrhs == 3) {
    if (!tz_mxIsDoubleScalar(prhs[2]))
      mxErrMsgTxt("The third arguments must be a double scalar.");
  }
}

/* tz_fcmp(x1, x2) compare two double float numbers.
 * tz_fcmp(x1, x2, ep) compare the double float numbers at the accuracy ep,
 * which is 1e-5 by default.
 */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  check_argin(nrhs, prhs);

  double epsilon = 1e-5;

  if (nrhs == 3) {
    epsilon = *mxGetPr(prhs[2]);
  }

  mwSize length;
  mwSize length1 = tz_mxGetL(prhs[0]);
  mwSize length2 = tz_mxGetL(prhs[1]);
  const mxArray *array = NULL;

  if (mxIsScalar(prhs[0]) && mxIsScalar(prhs[1])) {
    plhs[0] = mxCreateNumericMatrix(1, 1, mxINT8_CLASS, mxREAL);
    length = 1;
  } else if (!mxIsScalar(prhs[0]) && !mxIsScalar(prhs[1])) {
    plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(prhs[0]),
				   mxGetDimensions(prhs[0]), mxINT8_CLASS,
				   mxREAL);
    length = length1;
  } else {
    if (length1 == 1) {
      array = prhs[1];
      length = length2;
    } else {
      array = prhs[0];
      length = length1;
    }

    plhs[0] = mxCreateNumericArray(mxGetNumberOfDimensions(array),
				   mxGetDimensions(array), mxINT8_CLASS,
				   mxREAL);
  }

  INT8_T *y = (INT8_T *) mxGetPr(plhs[0]);
  double *x1 = mxGetPr(prhs[0]);
  double *x2 = mxGetPr(prhs[1]);

  mwSize i;
  if (length1 == length2) {
    for (i = 0; i < length; i++) {
      y[i] = (INT8_T) gsl_fcmp(x1[i], x2[i], epsilon);
    }
  } else if (length1 > length2) {
    for (i = 0; i < length; i++) {
      y[i] = (INT8_T) gsl_fcmp(x1[i], x2[0], epsilon);
    }
  } else {
    for (i = 0; i < length; i++) {
      y[i] = (INT8_T) gsl_fcmp(x1[0], x2[i], epsilon);
    }
  }
}


text:
#pragma once

#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <cctype>
#include <cerrno>
#include <chrono>
#include <climits>
#include <cmath>
#include <csignal>
#include <cstdarg>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <memory>
#include <numeric>
#include <sstream>
#include <stdexcept>
#include <string>
#include <thread>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>

#ifdef __linux__
#define LINUX
#endif

#ifdef _MSC_VER
#pragma warning(disable : 4267)
#endif

#include "Base64.h"
#include <GLTFAccessor.h>
#include <GLTFAsset.h>
#include <GLTFBuffer.h>
#include <GLTFBufferView.h>
#include <GLTFMesh.h>
#include <GLTFPrimitive.h>
#include <GLTFScene.h>
#include <GLTFTargetNames.h>

#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable : 4996)
#pragma warning(default : 4267)
#endif
#include "rapidjson/document.h"
#ifdef _MSC_VER
#pragma warning(pop)
#endif
#include "rapidjson/prettywriter.h"
#include "rapidjson/stringbuffer.h"
#include "rapidjson/writer.h"

#include <gsl/span>

#include <coveo/enumerable.h>
#include <coveo/linq.h>

#include <maya/M3dView.h>
#include <maya/MAnimControl.h>
#include <maya/MAnimUtil.h>
#include <maya/MArgDatabase.h>
#include <maya/MArgList.h>
#include <maya/MDagModifier.h>
#include <maya/MDagPath.h>
#include <maya/MDagPathArray.h>
#include <maya/MFileIO.h>
#include <maya/MFileObject.h>
#include <maya/MFloatMatrix.h>
#include <maya/MFloatPointArray.h>
#include <maya/MFloatVectorArray.h>
#include <maya/MFnAttribute.h>
#include <maya/MFnBlendShapeDeformer.h>
#include <maya/MFnBlinnShader.h>
#include <maya/MFnCamera.h>
#include <maya/MFnComponentListData.h>
#include <maya/MFnLambertShader.h>
#include <maya/MFnMatrixData.h>
#include <maya/MFnMesh.h>
#include <maya/MFnMessageAttribute.h>
#include <maya/MFnNumericAttribute.h>
#include <maya/MFnPhongShader.h>
#include <maya/MFnSet.h>
#include <maya/MFnSingleIndexedComponent.h>
#include <maya/MFnSkinCluster.h>
#include <maya/MFnStringArrayData.h>
#include <maya/MFnTransform.h>
#include <maya/MFnTypedAttribute.h>
#include <maya/MGlobal.h>
#include <maya/MIOStream.h>
#include <maya/MImage.h>
#include <maya/MItDependencyGraph.h>
#include <maya/MItDependencyNodes.h>
#include <maya/MItGeometry.h>
#include <maya/MItMeshFaceVertex.h>
#include <maya/MItMeshPolygon.h>
#include <maya/MMatrix.h>
#include <maya/MPointArray.h>
#include <maya/MPxCommand.h>
#include <maya/MQuaternion.h>
#include <maya/MRenderSetup.h>
#include <maya/MSelectionList.h>
#include <maya/MStreamUtils.h>
#include <maya/MSyntax.h>
#include <maya/MTime.h>
#include <maya/MUuid.h>

#ifdef isnan
#   undef isnan
#endif


text:
#pragma once
#include "Sampler.h"
#include "GradUtil.h"
#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#else
#include "FakeGSL.h"
#endif
#include "ActualEvaluators.h"

class BasicSampler : public Sampler {
	Interface* interf;
	BooleanDAG* dag;
	ActualEvaluators* actualEval;
	doublereal* xlow;
    doublereal* xupp;
    int ncontrols;

	gsl_vector* tmp;
	int RANDOM_SEARCH = 10;

    set<int> assertConstraints;
	int minimizeNode;

public:

    BasicSampler(BooleanDAG* dag_, Interface* interf_, ActualEvaluators* actualEval_, int ncontrols_,  doublereal* xlow_, doublereal* xupp_): dag(dag_), interf(interf_), actualEval(actualEval_), ncontrols(ncontrols_), xlow(xlow_), xupp(xupp_) {

    	minimizeNode = -1;

	    for (int i = 0; i < dag->size(); i++) { // TODO: this should also be set by caller class
	    	bool_node* n = (*dag)[i];
	    	if (n->type == bool_node::ASSERT && ((ASSERT_node*)n)->isHard()) {
	    		minimizeNode = i;
	    	} else if (n->type == bool_node::ASSERT || Util::isSqrt(n)) {
	    		assertConstraints.insert(i);
	    	} else if (n->type == bool_node::CTRL && n->getOtype() == OutType::BOOL) {
	    		assertConstraints.insert(i);
	    	}
	    }

	    tmp = gsl_vector_alloc(ncontrols);
    }

	virtual void sampleState(gsl_vector* state) {
		double best = GradUtil::MAXVAL;
        bool foundValid = false;
        for (int i = 0; i < RANDOM_SEARCH; i++) {
            randomize(tmp);
            cout << "Trying: ";
            for (int j = 0; j < tmp->size; j++) {
                cout << gsl_vector_get(tmp, j) << ", ";
            }
            actualEval->run(tmp);
            double error = getError();
            cout << "Error: " << error << endl;

            if (error < best) {
                best = error;
                gsl_vector_memcpy(state, tmp);
            }
        }
    }

    void randomize(gsl_vector* state) {
        for (int i = 0; i < state->size; i++) {
            double low = xlow[i];
            double high = xupp[i];
            double r = low + (rand() % (int)((high - low) * 10.0))/10.0;
            gsl_vector_set(state, i, r);
        }
    }

    double getError() {
        double error = 0.0;
        double e;
        if (minimizeNode >= 0) {
            e = actualEval->getErrorOnConstraint(minimizeNode);
            if (e > 0.0) {
                error += e; // TODO : change this
            }
        }
        for (auto it = assertConstraints.begin(); it != assertConstraints.end(); it++) {
            e = actualEval->getErrorOnConstraint(*it);
            if (e < 0.0) {
                error += -e;
            }
        }
        return error;
    }
};

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
#include <math.h>
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "matrix_operations.c"
#include "ellipsoid/parametric_function_roots.c"
#include "polynomial.c"

void characteristic_ellipsoid_matrix(double *X, double *R, double phi, double *rotax, double exponent)
{

    double w, x, y, z;
    w = cos(phi / 2.0);
    x = sin(phi / 2.0) * rotax[0];
    y = sin(phi / 2.0) * rotax[1];
    z = sin(phi / 2.0) * rotax[2];

    // Rotation matrix
    double Q[3][3] = {{ 0 }};

    Q[0][0] = 1.0 - 2.0 * (pow(y, 2) + pow(z, 2));
    Q[0][1] = 2.0 * (x * y + w * z);
    Q[0][2] = 2.0 * (x * z - w * y);

    Q[1][0] = 2.0 * (x * y - w * z);
    Q[1][1] = 1.0 - 2.0 * (pow(x, 2) + pow(z, 2));
    Q[1][2] = 2.0 * (y * z + w * x);

    Q[2][0] = 2.0 * (x * z + w * y);
    Q[2][1] = 2.0 * (y * z - w * x);
    Q[2][2] = 1.0 - 2.0 * (pow(x, 2) + pow(y, 2));

    //Rotation matrix transpose
    double Qt[3][3] = {{ 0 }};
    matrix_transpose(&Qt[0][0], &Q[0][0], 3, 3);

    // radii matrix
    double O[3][3] = {{ 0 }};
    double diag_vals[3];
    diag_vals[0] = pow(R[0], exponent * -2.0);
    diag_vals[1] = pow(R[1], exponent * -2.0);
    diag_vals[2] = pow(R[2], exponent * -2.0);
    set_diagonal(&O[0][0], diag_vals, 3, 3);

    // characteristic ellipse matrix
    double X_temp[3][3] = {{ 0 }};
    matrix_multiply(&X_temp[0][0], &O[0][0], &Q[0][0], 3, 3, 3);
    matrix_multiply(X, &Qt[0][0], &X_temp[0][0], 3, 3, 3);

}

double ellipsoid_overlap(double *rA, double *radiiA, double phiA, double *rotaxA,
    double *rB, double *radiiB, double phiB, double *rotaxB)
{

    // find XA^(-1) and XB^(1/2)
    double XA[3][3] = {{ 0 }};
    double XB[3][3] = {{ 0 }};
    characteristic_ellipsoid_matrix(&XA[0][0], &radiiA[0], phiA, &rotaxA[0], -1.0);
    characteristic_ellipsoid_matrix(&XB[0][0], &radiiB[0], phiB, &rotaxB[0], 0.5);

    // find A_AB
    double A_AB[3][3] = {{ 0 }};
    double A_temp[3][3] = {{ 0 }};
    matrix_multiply(&A_temp[0][0], &XA[0][0], &XB[0][0], 3, 3, 3);
    matrix_multiply(&A_AB[0][0], &XB[0][0], &A_temp[0][0], 3, 3, 3);

    // find r_AB
    double rAB[3];
    rAB[0] = rB[0] - rA[0];
    rAB[1] = rB[1] - rA[1];
    rAB[2] = rB[2] - rA[2];

    // find a_AB
    double a_AB[3];
    matrix_multiply(&a_AB[0], &XB[0][0], &rAB[0], 3, 3, 1);

    // extract elements of the matrix A_AB and a_AB
    double a11, a12, a13, a21, a22, a23, a31, a32, a33;
    double b1, b2, b3;
    a11 = A_AB[0][0];
    a12 = A_AB[0][1];
    a13 = A_AB[0][2];
    a21 = A_AB[1][0];
    a22 = A_AB[1][1];
    a23 = A_AB[1][2];
    a31 = A_AB[2][0];
    a32 = A_AB[2][1];
    a33 = A_AB[2][2];
    b1 = a_AB[0];
    b2 = a_AB[1];
    b3 = a_AB[2];

    // find coefficients for the parametric polynomial derivative used to find max
    double h[7];
    double z[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
    size_t n = 7;
    h[0] = h0(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[1] = h1(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[2] = h2(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[3] = h3(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[4] = h4(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[5] = h5(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);
    h[6] = h6(a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);

    // find roots
    size_t m;
    m = find_roots(&z[0], &h[0], n);

    double F;
    int i;

    if (m > 1)
    {
        if (f(0, a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3) > f(1, a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3))
        {

           F = f(0, a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);

        }
        else
        {

           F = f(1, a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);

        }

        for(i=0; i<=m-2; i++)
        {

           if (( z[2 * i + 1] == 0 ) & (z[2 * i] > 0) & (z[2 * i] < 1))
           {

              F = f(z[2 * i], a11, a12, a13, a21, a22, a23, a31, a32, a33, b1, b2, b3);

           }

        }
    }
    else
    {
        F = 0.;
    }

    return F;

}

double container_cube_overlap_potential(double *rA, double *radiiA, double phiA, double *rotaxA)
{

    double rB[3];
    double radiiB[3];
    double phiB;
    double rotaxB[3];

    double top, bottom, left, right, front, back;

    // top
    rB[0] = 0.5;
    rB[1] = 0.5;
    rB[2] = 2;

    radiiB[0] = INFINITY;
    radiiB[1] = INFINITY;
    radiiB[2] = 1;

    rotaxB[0] = 1;
    rotaxB[1] = 0;
    rotaxB[2] = 0;

    phiB = 0.0;

    top = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    // bottom
    rB[0] = 0.5;
    rB[1] = 0.5;
    rB[2] = -1;

    radiiB[0] = INFINITY;
    radiiB[1] = INFINITY;
    radiiB[2] = 1.0;

    rotaxB[0] = 1;
    rotaxB[1] = 0;
    rotaxB[2] = 0;

    phiB = 0;

    bottom = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    // left
    rB[0] = 0.5;
    rB[1] = -1.0;
    rB[2] = 0.5;

    radiiB[0] = INFINITY;
    radiiB[1] = 1;
    radiiB[2] = INFINITY;

    rotaxB[0] = 1;
    rotaxB[1] = 0;
    rotaxB[2] = 0;

    phiB = 0;

    left = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    // right
    rB[0] = 0.5;
    rB[1] = 2;
    rB[2] = 0.5;

    radiiB[0] = INFINITY;
    radiiB[1] = 1;
    radiiB[2] = INFINITY;

    rotaxB[0] = 1;
    rotaxB[1] = 0;
    rotaxB[2] = 0;

    phiB = 0;

    right = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    // front
    rB[0] = -1.0;
    rB[1] = 0.5;
    rB[2] = 0.5;

    radiiB[0] = 1.0;
    radiiB[1] = INFINITY;
    radiiB[2] = INFINITY;

    rotaxB[0] = 1.0;
    rotaxB[1] = 0.0;
    rotaxB[2] = 0.0;

    phiB = 0.0;

    front = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    // back
    rB[0] = 2.0;
    rB[1] = 0.5;
    rB[2] = 0.5;

    radiiB[0] = 1;
    radiiB[1] = INFINITY;
    radiiB[2] = INFINITY;

    rotaxB[0] = 1.0;
    rotaxB[1] = 0.0;
    rotaxB[2] = 0.0;

    phiB = 0;

    back = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0], &rB[0], &radiiB[0], phiB, &rotaxB[0]);

    return fminf(top, fminf(bottom, fminf(left, fminf(right, fminf(front, back)))));

}

void random_point_on_sphere(double *rotax, double theta, double z)
{

    double x, y;

    x = sqrt(1 - pow(z, 2)) * cos(theta);
    y = sqrt(1 - pow(z, 2)) * sin(theta);

    rotax[0] = x;
    rotax[1] = y;
    rotax[2] = z;

}

size_t rsa_align_cube(double *x, double *y, double *z,
    size_t npoints, double *radius, double *rotax, double phi,
    int step_limit, unsigned long randSeed)
{

    // Setup GSL random number generator
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    // Set the seed
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    double rA[3];
    double radiiA[3];
    double phiA;
    double rotaxA[3];

    double rB[3];
    double radiiB[3];
    double phiB;
    double rotaxB[3];

    size_t valid_pts;
    double F, G;
    int k, flag, step;

    step = 0;
    valid_pts = 1;

    // Get new ellipsoid position and orientation
    G = 0;
    while (G < 1)
    {

        rB[0] = gsl_rng_uniform (r);
        rB[1] = gsl_rng_uniform (r);
        rB[2] = gsl_rng_uniform (r);
        radiiB[0] = radius[0];
        radiiB[1] = radius[1];
        radiiB[2] = radius[2];
        phiB = phi;
        rotaxB[0] = rotax[0];
        rotaxB[1] = rotax[1];
        rotaxB[2] = rotax[2];

        G = container_cube_overlap_potential(&rB[0], &radiiB[0], phiB, &rotaxB[0]);

    }

    // store ellipsoid position
    x[0] = rB[0];
    y[0] = rB[1];
    z[0] = rB[2];

    // Generate ellipsoid positions
    while ((valid_pts < npoints) & (step < step_limit))
    {

        // Get new ellipsoid position and orientation
        G = 0;
        while (G < 1)
        {

            rB[0] = gsl_rng_uniform (r);
            rB[1] = gsl_rng_uniform (r);
            rB[2] = gsl_rng_uniform (r);
            radiiB[0] = radius[0];
            radiiB[1] = radius[1];
            radiiB[2] = radius[2];
            phiB = phi;
            rotaxB[0] = rotax[0];
            rotaxB[1] = rotax[1];
            rotaxB[2] = rotax[2];

            G = container_cube_overlap_potential(&rB[0], &radiiB[0], phiB, &rotaxB[0]);

        }

        flag = 1;
        for (k = 0; k < valid_pts; k++)
        {

            rA[0] = x[k];
            rA[1] = y[k];
            rA[2] = z[k];
            radiiA[0] = radius[0];
            radiiA[1] = radius[1];
            radiiA[2] = radius[2];
            phiA = phi;
            rotaxA[0] = rotax[0];
            rotaxA[1] = rotax[1];
            rotaxA[2] = rotax[2];

            F = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0],
                &rB[0], &radiiB[0], phiB, &rotaxB[0]);

            if (F < 1.0)
            {

                flag = 0;
                break;

            }
        }
        if (flag == 1)
        {

           x[valid_pts] = rB[0];
           y[valid_pts] = rB[1];
           z[valid_pts] = rB[2];
           valid_pts += 1;

        }

        step += 1;

    }

    gsl_rng_free (r);

    return valid_pts;

}

size_t rsa_cube(double *x, double *y, double *z,
    size_t npoints, double *radius,
    double *rtx, double *rty, double *rtz, double *phi,
    int step_limit, unsigned long randSeed)
{

    // Setup GSL random number generator
    const gsl_rng_type * T;
    gsl_rng * r;
    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    // Set the seed
    // srand ( time(NULL) );
    // unsigned long randSeed = rand();
    gsl_rng_set(r, randSeed);

    double rA[3];
    double radiiA[3];
    double phiA;
    double rotaxA[3];

    double rB[3];
    double radiiB[3];
    double phiB;
    double rotaxB[3];

    double theta, sz;

    size_t valid_pts;
    double F, G;
    int k, flag, step;

    step = 0;
    valid_pts = 1;

    // Get new ellipsoid position and orientation
    G = 0;
    while (G < 1)
    {

        rB[0] = gsl_rng_uniform (r);
        rB[1] = gsl_rng_uniform (r);
        rB[2] = gsl_rng_uniform (r);
        radiiB[0] = radius[0];
        radiiB[1] = radius[1];
        radiiB[2] = radius[2];
        phiB = 2 * M_PI * gsl_rng_uniform (r);

        // random rotation axis
        theta = 2 * M_PI * gsl_rng_uniform (r);
        sz = 2 * gsl_rng_uniform (r) - 1;
        random_point_on_sphere(&rotaxB[0], theta, sz);

        G = container_cube_overlap_potential(&rB[0], &radiiB[0], phiB, &rotaxB[0]);

    }

    // store ellipsoid position
    x[0] = rB[0];
    y[0] = rB[1];
    z[0] = rB[2];
    rtx[0] = rotaxB[0];
    rty[0] = rotaxB[1];
    rtz[0] = rotaxB[2];
    phi[0] = phiB;

    // Generate ellipsoid positions
    while ((valid_pts < npoints) & (step < step_limit))
    {

        // Get new ellipsoid position and orientation
        G = 0;
        while (G < 1)
        {

            rB[0] = gsl_rng_uniform (r);
            rB[1] = gsl_rng_uniform (r);
            rB[2] = gsl_rng_uniform (r);
            radiiB[0] = radius[0];
            radiiB[1] = radius[1];
            radiiB[2] = radius[2];
            phiB = 2 * M_PI * gsl_rng_uniform (r);

            // random rotation axis
            theta = 2 * M_PI * gsl_rng_uniform (r);
            sz = 2 * gsl_rng_uniform (r) - 1;
            random_point_on_sphere(&rotaxB[0], theta, sz);

            G = container_cube_overlap_potential(&rB[0], &radiiB[0], phiB, &rotaxB[0]);

        }

        flag = 1;
        for (k = 0; k < valid_pts; k++)
        {

            rA[0] = x[k];
            rA[1] = y[k];
            rA[2] = z[k];
            radiiA[0] = radius[0];
            radiiA[1] = radius[1];
            radiiA[2] = radius[2];
            phiA = phi[k];
            rotaxA[0] = rtx[k];
            rotaxA[1] = rty[k];
            rotaxA[2] = rtz[k];

            F = ellipsoid_overlap(&rA[0], &radiiA[0], phiA, &rotaxA[0],
                &rB[0], &radiiB[0], phiB, &rotaxB[0]);

            if (F < 1.0)
            {

                flag = 0;
                break;

            }
        }
        if (flag == 1)
        {

            x[valid_pts] = rB[0];
            y[valid_pts] = rB[1];
            z[valid_pts] = rB[2];
            rtx[valid_pts] = rotaxB[0];
            rty[valid_pts] = rotaxB[1];
            rtz[valid_pts] = rotaxB[2];
            phi[valid_pts] = phiB;

            valid_pts += 1;

        }

        step += 1;

    }

    gsl_rng_free (r);

    return valid_pts;

}


text:
#pragma once
#include <gsl/span>
#include "halley/data_structures/vector.h"
#include "component_schema.h"
#include "custom_type_schema.h"
#include "halley/data_structures/hash_map.h"
#include "halley/text/halleystring.h"
#include "message_schema.h"
#include "system_message_schema.h"
#include "system_schema.h"

namespace YAML
{
	class Node;
}

namespace Halley {
	struct CodegenSourceInfo {
		String filename;
		gsl::span<const gsl::byte> data;
		bool generate = false;
	};

	class Codegen;

    class ECSData {
    public:
		void loadSources(Vector<CodegenSourceInfo> files);

		const HashMap<String, ComponentSchema>& getComponents() const;
		const HashMap<String, SystemSchema>& getSystems() const;
		const HashMap<String, MessageSchema>& getMessages() const;
		const HashMap<String, SystemMessageSchema>& getSystemMessages() const;
		const HashMap<String, CustomTypeSchema>& getCustomTypes() const;

    	void clear();
		int getRevision() const;

	private:
    	void addSource(CodegenSourceInfo sourceInfo);
		void addComponent(YAML::Node rootNode, bool generate);
		void addSystem(YAML::Node rootNode, bool generate);
		void addMessage(YAML::Node rootNode, bool generate);
    	void addSystemMessage(YAML::Node rootNode, bool generate);
		void addType(YAML::Node rootNode);
		String getInclude(String typeName) const;

    	void validate();
		void process();

    	HashMap<String, ComponentSchema> components;
		HashMap<String, SystemSchema> systems;
		HashMap<String, MessageSchema> messages;
    	HashMap<String, SystemMessageSchema> systemMessages;
		HashMap<String, CustomTypeSchema> types;
		int revision = 0;
    };
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
#ifndef __INTEGRAL_H__
#define __INTEGRAL_H__

#include "basis.h"
#include "gslextra.h"

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf.h>

typedef struct gaussian_chain{
    double R[3];
    int a[3];

    double exponent;

    double coefficient;

    gaussian_chain * NEXT;

}gaussian_chain;

gaussian_chain * gaussian_chain_calloc();

void gaussian_chain_free(gaussian_chain * HEAD);

double gaussian_chain_get(gaussian_chain * HEAD,double x, double y, double z);
double orbital_get(orbital * orbital,double x, double y, double z);

double SIntegral(double ra[3], double rb[3], int ax, int ay, int az, int bx, int by, int bz, double alpha,double beta);

double gaussian_chain_SIntegral(gaussian_chain * a, gaussian_chain * b);
double gaussian_chain_full_SIntegral(gaussian_chain * a_HEAD, gaussian_chain * b_HEAD);
void single_electron_transform(gaussian_chain * HEAD, orbital * a);

double orbital_SIntegral(orbital * a,orbital * b);

void orbital_S_matrix(gsl_matrix * dest, orbital * HEAD, int length);

#endif

max_stars_repo_name: Walter-Feng/Overlap

max_issues_repo_name: Walter-Feng/Overlap

max_forks_repo_name: Walter-Feng/Overlap


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

#pragma once

#include <cblas.h>
#include <cuda_runtime_api.h>

#include <boost/shared_ptr.hpp>
#include <fstream>
#include <iostream>
#include <map>
#include <string>
#include <vector>

#include "modules/perception/base/blob.h"
#include "modules/perception/base/image.h"

namespace apollo {
namespace perception {
namespace inference {
class GPUL2Norm {
 public:
  void L2Norm(base::Blob<float> *input_data);

 private:
  base::Blob<float> scale_;
  base::Blob<float> ones_;
  base::Blob<float> square_;
};

void GPUGemmFloat(const CBLAS_TRANSPOSE TransA, const CBLAS_TRANSPOSE TransB,
                  const int M, const int N, const int K, const float alpha,
                  const float *A, const float *B, const float beta, float *C);
void GPUMultiFloat(const int n, const float *a, const float *b, float *result);
void GPUMSetFloat(const int n, const float alpha, float *result);

}  // namespace inference
}  // namespace perception
}  // namespace apollo

max_stars_repo_name: ghdawn/apollo


text:
/* specfunc/hyperg.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

/* Miscellaneous implementations of use
 * for evaluation of hypergeometric functions.
 */
#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_exp.h"
#include "gsl_sf_gamma.h"

#include "error.h"
#include "hyperg.h"

#define SUM_LARGE  (1.0e-5*GSL_DBL_MAX)

int
gsl_sf_hyperg_1F1_series_e(const double a, const double b, const double x,
                              gsl_sf_result * result
                              )
{
  double an  = a;
  double bn  = b;
  double n   = 1.0;
  double del = 1.0;
  double abs_del = 1.0;
  double max_abs_del = 1.0;
  double sum_val = 1.0;
  double sum_err = 0.0;

  while(abs_del/fabs(sum_val) > GSL_DBL_EPSILON) {
    double u, abs_u;

    if(bn == 0.0) {
      DOMAIN_ERROR(result);
    }
    if(an == 0.0 || n > 1000.0) {
      result->val  = sum_val;
      result->err  = sum_err;
      result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);
      return GSL_SUCCESS;
    }

    u = x * (an/(bn*n));
    abs_u = fabs(u);
    if(abs_u > 1.0 && max_abs_del > GSL_DBL_MAX/abs_u) {
      result->val = sum_val;
      result->err = fabs(sum_val);
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }
    del *= u;
    sum_val += del;
    if(fabs(sum_val) > SUM_LARGE) {
      result->val = sum_val;
      result->err = fabs(sum_val);
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }

    abs_del = fabs(del);
    max_abs_del = GSL_MAX_DBL(abs_del, max_abs_del);
    sum_err += 2.0*GSL_DBL_EPSILON*abs_del;

    an += 1.0;
    bn += 1.0;
    n  += 1.0;
  }

  result->val  = sum_val;
  result->err  = sum_err;
  result->err += abs_del;
  result->err += 2.0 * GSL_DBL_EPSILON * n * fabs(sum_val);

  return GSL_SUCCESS;
}

int
gsl_sf_hyperg_1F1_large_b_e(const double a, const double b, const double x, gsl_sf_result * result)
{
  if(fabs(x/b) < 1.0) {
    const double u = x/b;
    const double v = 1.0/(1.0-u);
    const double pre = pow(v,a);
    const double uv  = u*v;
    const double uv2 = uv*uv;
    const double t1  = a*(a+1.0)/(2.0*b)*uv2;
    const double t2a = a*(a+1.0)/(24.0*b*b)*uv2;
    const double t2b = 12.0 + 16.0*(a+2.0)*uv + 3.0*(a+2.0)*(a+3.0)*uv2;
    const double t2  = t2a*t2b;
    result->val  = pre * (1.0 - t1 + t2);
    result->err  = pre * GSL_DBL_EPSILON * (1.0 + fabs(t1) + fabs(t2));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    DOMAIN_ERROR(result);
  }
}

int
gsl_sf_hyperg_U_large_b_e(const double a, const double b, const double x,
                             gsl_sf_result * result,
			     double * ln_multiplier
			     )
{
  double N   = floor(b);  /* b = N + eps */
  double eps = b - N;

  if(fabs(eps) < GSL_SQRT_DBL_EPSILON) {
    double lnpre_val;
    double lnpre_err;
    gsl_sf_result M;
    if(b > 1.0) {
      double tmp = (1.0-b)*log(x);
      gsl_sf_result lg_bm1;
      gsl_sf_result lg_a;
      gsl_sf_lngamma_e(b-1.0, &lg_bm1);
      gsl_sf_lngamma_e(a, &lg_a);
      lnpre_val = tmp + x + lg_bm1.val - lg_a.val;
      lnpre_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(x) + fabs(tmp));
      gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, -x, &M);
    }
    else {
      gsl_sf_result lg_1mb;
      gsl_sf_result lg_1pamb;
      gsl_sf_lngamma_e(1.0-b, &lg_1mb);
      gsl_sf_lngamma_e(1.0+a-b, &lg_1pamb);
      lnpre_val = lg_1mb.val - lg_1pamb.val;
      lnpre_err = lg_1mb.err + lg_1pamb.err;
      gsl_sf_hyperg_1F1_large_b_e(a, b, x, &M);
    }

    if(lnpre_val > GSL_LOG_DBL_MAX-10.0) {
      result->val  = M.val;
      result->err  = M.err;
      *ln_multiplier = lnpre_val;
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }
    else {
      gsl_sf_result epre;
      int stat_e = gsl_sf_exp_err_e(lnpre_val, lnpre_err, &epre);
      result->val  = epre.val * M.val;
      result->err  = epre.val * M.err + epre.err * fabs(M.val);
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      *ln_multiplier = 0.0;
      return stat_e;
    }
  }
  else {
    double omb_lnx = (1.0-b)*log(x);
    gsl_sf_result lg_1mb;    double sgn_1mb;
    gsl_sf_result lg_1pamb;  double sgn_1pamb;
    gsl_sf_result lg_bm1;    double sgn_bm1;
    gsl_sf_result lg_a;      double sgn_a;
    gsl_sf_result M1, M2;
    double lnpre1_val, lnpre2_val;
    double lnpre1_err, lnpre2_err;
    double sgpre1, sgpre2;
    gsl_sf_hyperg_1F1_large_b_e(    a,     b, x, &M1);
    gsl_sf_hyperg_1F1_large_b_e(1.0-a, 2.0-b, x, &M2);

    gsl_sf_lngamma_sgn_e(1.0-b,   &lg_1mb,   &sgn_1mb);
    gsl_sf_lngamma_sgn_e(1.0+a-b, &lg_1pamb, &sgn_1pamb);

    gsl_sf_lngamma_sgn_e(b-1.0, &lg_bm1, &sgn_bm1);
    gsl_sf_lngamma_sgn_e(a,     &lg_a,   &sgn_a);

    lnpre1_val = lg_1mb.val - lg_1pamb.val;
    lnpre1_err = lg_1mb.err + lg_1pamb.err;
    lnpre2_val = lg_bm1.val - lg_a.val - omb_lnx - x;
    lnpre2_err = lg_bm1.err + lg_a.err + GSL_DBL_EPSILON * (fabs(omb_lnx)+fabs(x));
    sgpre1 = sgn_1mb * sgn_1pamb;
    sgpre2 = sgn_bm1 * sgn_a;

    if(lnpre1_val > GSL_LOG_DBL_MAX-10.0 || lnpre2_val > GSL_LOG_DBL_MAX-10.0) {
      double max_lnpre_val = GSL_MAX(lnpre1_val,lnpre2_val);
      double max_lnpre_err = GSL_MAX(lnpre1_err,lnpre2_err);
      double lp1 = lnpre1_val - max_lnpre_val;
      double lp2 = lnpre2_val - max_lnpre_val;
      double t1  = sgpre1*exp(lp1);
      double t2  = sgpre2*exp(lp2);
      result->val  = t1*M1.val + t2*M2.val;
      result->err  = fabs(t1)*M1.err + fabs(t2)*M2.err;
      result->err += GSL_DBL_EPSILON * exp(max_lnpre_err) * (fabs(t1*M1.val) + fabs(t2*M2.val));
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      *ln_multiplier = max_lnpre_val;
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }
    else {
      double t1 = sgpre1*exp(lnpre1_val);
      double t2 = sgpre2*exp(lnpre2_val);
      result->val  = t1*M1.val + t2*M2.val;
      result->err  = fabs(t1) * M1.err + fabs(t2)*M2.err;
      result->err += GSL_DBL_EPSILON * (exp(lnpre1_err)*fabs(t1*M1.val) + exp(lnpre2_err)*fabs(t2*M2.val));
      result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
      *ln_multiplier = 0.0;
      return GSL_SUCCESS;
    }
  }
}

/* [Carlson, p.109] says the error in truncating this asymptotic series
 * is less than the absolute value of the first neglected term.
 *
 * A termination argument is provided, so that the series will
 * be summed at most up to n=n_trunc. If n_trunc is set negative,
 * then the series is summed until it appears to start diverging.
 */
int
gsl_sf_hyperg_2F0_series_e(const double a, const double b, const double x,
                              int n_trunc,
                              gsl_sf_result * result
                              )
{
  const int maxiter = 2000;
  double an = a;
  double bn = b;
  double n   = 1.0;
  double sum = 1.0;
  double del = 1.0;
  double abs_del = 1.0;
  double max_abs_del = 1.0;
  double last_abs_del = 1.0;

  while(abs_del/fabs(sum) > GSL_DBL_EPSILON && n < maxiter) {

    double u = an * (bn/n * x);
    double abs_u = fabs(u);

    if(abs_u > 1.0 && (max_abs_del > GSL_DBL_MAX/abs_u)) {
      result->val = sum;
      result->err = fabs(sum);
      GSL_ERROR ("overflow", GSL_EOVRFLW);
    }

    del *= u;
    sum += del;

    abs_del = fabs(del);

    if(abs_del > last_abs_del) break; /* series is probably starting to grow */

    last_abs_del = abs_del;
    max_abs_del  = GSL_MAX(abs_del, max_abs_del);

    an += 1.0;
    bn += 1.0;
    n  += 1.0;

    if(an == 0.0 || bn == 0.0) break;        /* series terminated */

    if(n_trunc >= 0 && n >= n_trunc) break;  /* reached requested timeout */
  }

  result->val = sum;
  result->err = GSL_DBL_EPSILON * n + abs_del;
  if(n >= maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}


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

#ifndef FFTW_MPI_H
#define FFTW_MPI_H

#include <fftw.h>
#include <mpi.h> /* need access to the MPI type definitions */

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

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

typedef fftw_real TRANSPOSE_EL_TYPE;

typedef struct {
     int block_num, dest_pe, send_size, recv_size;
} transpose_mpi_exchange;

typedef struct {
     MPI_Comm comm;
     int n_pes, my_pe;

     int nx,ny,local_nx,local_ny;

     transpose_mpi_exchange *exchange;
     int num_steps, send_block_size, recv_block_size;

     MPI_Datatype el_type;

     MPI_Request request[2];

     int *perm_block_dest;
     int num_perm_blocks, perm_block_size;

     int all_blocks_equal;
     int *send_block_sizes, *send_block_offsets;
     int *recv_block_sizes, *recv_block_offsets;

     char *move;
     int move_size;
} transpose_mpi_plan_struct;

typedef transpose_mpi_plan_struct *transpose_mpi_plan;

extern void transpose_mpi_get_local_size(int n, int my_pe, int n_pes,
					 int *local_n, int *local_start);
extern int transpose_mpi_get_local_storage_size(int nx, int ny,
						int my_pe, int n_pes);

extern transpose_mpi_plan transpose_mpi_create_plan(int nx, int ny,
						    MPI_Comm comm);
extern void transpose_mpi_destroy_plan(transpose_mpi_plan p);

extern void transpose_mpi(transpose_mpi_plan p, int el_size,
			  TRANSPOSE_EL_TYPE *local_data,
			  TRANSPOSE_EL_TYPE *work);

typedef enum { BEFORE_TRANSPOSE, AFTER_TRANSPOSE } transpose_in_place_which;

typedef enum { TRANSPOSE_SYNC, TRANSPOSE_ASYNC } transpose_sync_type;

extern void transpose_in_place_local(transpose_mpi_plan p,
                              int el_size, TRANSPOSE_EL_TYPE *local_data,
                              transpose_in_place_which which);

extern TRANSPOSE_EL_TYPE *transpose_allocate_send_buf(transpose_mpi_plan p,
						      int el_size);
extern void transpose_get_send_block(transpose_mpi_plan p, int step,
				     int *block_y_start, int *block_ny);
extern void transpose_start_exchange_step(transpose_mpi_plan p,
					  int el_size,
					  TRANSPOSE_EL_TYPE *local_data,
					  TRANSPOSE_EL_TYPE *send_buf,
					  int step,
					  transpose_sync_type sync_type);
extern void transpose_finish_exchange_step(transpose_mpi_plan p, int step);

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

typedef struct {
     fftw_plan p_fft_x;  /* plan for first dimension */
     fftwnd_plan p_fft;  /* plan for subsequent dimensions */
     transpose_mpi_plan p_transpose, p_transpose_inv;
     fftw_complex *work; /* extra workspace, if needed */
} fftwnd_mpi_plan_data;

typedef fftwnd_mpi_plan_data *fftwnd_mpi_plan;

typedef enum {
    FFTW_NORMAL_ORDER,
    FFTW_TRANSPOSED_ORDER
} fftwnd_mpi_output_order;

extern fftwnd_mpi_plan fftwnd_mpi_create_plan(MPI_Comm comm,
					      int rank, const int *n,
					      fftw_direction dir,
					      int flags);
extern fftwnd_mpi_plan fftw2d_mpi_create_plan(MPI_Comm comm,
					      int nx, int ny,
					  fftw_direction dir, int flags);
extern fftwnd_mpi_plan fftw3d_mpi_create_plan(MPI_Comm comm,
					      int nx, int ny, int nz,
					  fftw_direction dir, int flags);

extern void fftwnd_mpi_destroy_plan(fftwnd_mpi_plan p);

extern void fftwnd_mpi_local_sizes(fftwnd_mpi_plan p,
				   int *local_nx,
				   int *local_x_start,
				   int *local_ny_after_transpose,
				   int *local_y_start_after_transpose,
				   int *total_local_size);

extern void fftwnd_mpi(fftwnd_mpi_plan p,
		       int n_fields,
		       fftw_complex *local_data, fftw_complex *work,
		       fftwnd_mpi_output_order output_order);

extern void fftw_mpi_die(const char *error_string);

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

typedef struct fftw_mpi_twiddle_struct {
     int rows, rowstart, cols, n;
     fftw_complex *W;
     int refcount;
     struct fftw_mpi_twiddle_struct *next;
} fftw_mpi_twiddle;

typedef struct fftw_mpi_plan_struct {
     int n, m, r, local_m, local_m_start, local_r, local_r_start;
     fftw_complex *fft_work;
     fftw_mpi_twiddle *tw;
     transpose_mpi_plan p_transpose, p_transpose_inv;
     fftw_plan pm, pr;
     int flags;
     fftw_direction dir;
} *fftw_mpi_plan;

/* new flags for the MPI planner: */
#define FFTW_SCRAMBLED_INPUT (8192)
#define FFTW_SCRAMBLED_OUTPUT (16384)

extern void fftw_mpi_local_sizes(fftw_mpi_plan p,
				 int *local_n,
				 int *local_start,
				 int *local_n_after_transform,
				 int *local_start_after_transform,
				 int *total_local_size);

extern fftw_mpi_plan fftw_mpi_create_plan(MPI_Comm comm,
					  int n,
					  fftw_direction dir, int flags);

extern void fftw_mpi_destroy_plan(fftw_mpi_plan p);

extern void fftw_mpi(fftw_mpi_plan p, int n_fields,
		     fftw_complex *local_data, fftw_complex *work);

extern void fftw_mpi_print_plan(fftw_mpi_plan p);

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

#ifdef __cplusplus
} /* extern "C" */
#endif /* __cplusplus */

#endif /* FFTW_MPI_H */

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_syrk (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = -1.0f;
   float beta = 0.1f;
   float A[] = { 0.412f, -0.229f };
   int lda = 1;
   float C[] = { 0.628f, -0.664f, -0.268f, 0.096f };
   int ldc = 2;
   float C_expected[] = { -0.106944f, 0.027948f, -0.268f, -0.042841f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1566)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = -1.0f;
   float beta = 0.1f;
   float A[] = { 0.101f, -0.653f };
   int lda = 2;
   float C[] = { 0.432f, 0.107f, -0.952f, -0.532f };
   int ldc = 2;
   float C_expected[] = { 0.032999f, 0.107f, -0.029247f, -0.479609f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1567)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = 1.0f;
   float beta = 0.1f;
   float A[] = { 0.79f, 0.595f };
   int lda = 2;
   float C[] = { 0.257f, 0.183f, -0.021f, -0.053f };
   int ldc = 2;
   float C_expected[] = { 0.6498f, 0.48835f, -0.021f, 0.348725f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1568)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = 1.0f;
   float beta = 0.1f;
   float A[] = { -0.181f, -0.654f };
   int lda = 1;
   float C[] = { -0.4f, 0.615f, 0.147f, -0.163f };
   int ldc = 2;
   float C_expected[] = { -0.007239f, 0.615f, 0.133074f, 0.411416f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1569)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = 0.0f;
   float beta = -1.0f;
   float A[] = { -0.191f, 0.584f };
   int lda = 1;
   float C[] = { -0.719f, -0.681f, -0.003f, 0.544f };
   int ldc = 2;
   float C_expected[] = { 0.719f, -0.681f, 0.003f, -0.544f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1570)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha = 0.0f;
   float beta = -1.0f;
   float A[] = { 0.788f, 0.041f };
   int lda = 2;
   float C[] = { 0.029f, 0.365f, 0.739f, -0.769f };
   int ldc = 2;
   float C_expected[] = { -0.029f, -0.365f, 0.739f, 0.769f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1571)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = -0.3f;
   float beta = -1.0f;
   float A[] = { 0.733f, 0.678f };
   int lda = 2;
   float C[] = { -0.941f, 0.96f, 0.07f, -0.295f };
   int ldc = 2;
   float C_expected[] = { 0.779813f, 0.96f, -0.219092f, 0.157095f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1572)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha = -0.3f;
   float beta = -1.0f;
   float A[] = { -0.87f, 0.675f };
   int lda = 1;
   float C[] = { -0.602f, -0.432f, -0.984f, 0.384f };
   int ldc = 2;
   float C_expected[] = { 0.37493f, 0.608175f, -0.984f, -0.520687f };
   cblas_ssyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyrk(case 1573)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 0.1;
   double beta = -0.3;
   double A[] = { 0.169, -0.875 };
   int lda = 1;
   double C[] = { 0.159, 0.277, 0.865, 0.346 };
   int ldc = 2;
   double C_expected[] = { -0.0448439, -0.0978875, 0.865, -0.0272375 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1574)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 0.1;
   double beta = -0.3;
   double A[] = { 0.536, -0.725 };
   int lda = 2;
   double C[] = { 0.154, -0.445, -0.841, -0.91 };
   int ldc = 2;
   double C_expected[] = { -0.0174704, -0.445, 0.21344, 0.3255625 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1575)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = 0;
   double beta = -1;
   double A[] = { -0.07, 0.8 };
   int lda = 2;
   double C[] = { 0.823, -0.88, -0.136, 0.793 };
   int ldc = 2;
   double C_expected[] = { -0.823, 0.88, -0.136, -0.793 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1576)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = 0;
   double beta = -1;
   double A[] = { -0.058, 0.649 };
   int lda = 1;
   double C[] = { -0.187, 0.294, -0.004, -0.933 };
   int ldc = 2;
   double C_expected[] = { 0.187, 0.294, 0.004, 0.933 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1577)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 1;
   double beta = -1;
   double A[] = { 0.263, -0.289 };
   int lda = 1;
   double C[] = { 0.554, -0.679, 0.993, 0.758 };
   int ldc = 2;
   double C_expected[] = { -0.484831, -0.679, -1.069007, -0.674479 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1578)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha = 1;
   double beta = -1;
   double A[] = { -0.265, -0.837 };
   int lda = 2;
   double C[] = { -0.994, 0.967, -0.34, -0.069 };
   int ldc = 2;
   double C_expected[] = { 1.064225, -0.745195, -0.34, 0.769569 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1579)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { -0.464, 0.394 };
   int lda = 2;
   double C[] = { -0.45, -0.447, 0.649, 0.055 };
   int ldc = 2;
   double C_expected[] = { -0.5145888, -0.447, 0.7038448, 0.0084292 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1580)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { 0.815, 0.168 };
   int lda = 1;
   double C[] = { 0.817, -0.957, -0.395, -0.382 };
   int ldc = 2;
   double C_expected[] = { 0.6177325, -0.998076, -0.395, -0.3904672 };
   cblas_dsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyrk(case 1581)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.447f, -0.507f, -0.425f, 0.701f };
   int lda = 1;
   float C[] = { 0.16f, -0.245f, 0.922f, -0.437f, 0.24f, 0.008f, -0.095f, 0.749f };
   int ldc = 2;
   float C_expected[] = { -0.0235f, 0.0895f, -0.2329f, 0.2233f, 0.24f, 0.008f, -0.0464f, -0.2342f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1582) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1582) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { -0.421f, -0.435f, -0.914f, -0.493f };
   int lda = 2;
   float C[] = { -0.761f, -0.38f, 0.043f, -0.999f, 0.779f, 0.238f, 0.082f, 0.394f };
   int ldc = 2;
   float C_expected[] = { 0.2663f, 0.0379f, 0.043f, -0.999f, -0.2575f, 0.0065f, -0.064f, -0.11f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1583) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1583) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.827f, -0.896f, 0.417f, 0.865f };
   int lda = 2;
   float C[] = { -0.349f, -0.31f, 0.972f, 0.794f, -0.906f, -0.595f, -0.089f, -0.333f };
   int ldc = 2;
   float C_expected[] = { 0.254587f, 1.54008f, -1.4909f, -0.482723f, -0.906f, -0.595f, 0.634336f, -0.63041f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1584) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1584) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.607f, 0.747f, -0.889f, 0.333f };
   int lda = 1;
   float C[] = { 0.244f, 0.564f, 0.009f, 0.578f, -0.827f, 0.558f, -0.337f, 0.731f };
   int ldc = 2;
   float C_expected[] = { 0.05996f, -1.05166f, 0.009f, 0.578f, 0.980674f, 0.211852f, -0.651432f, 0.339074f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1585) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1585) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 113;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.784f, -0.281f, -0.88f, 0.479f };
   int lda = 2;
   float C[] = { 0.491f, 0.531f, 0.805f, -0.097f, 0.728f, 0.674f, -0.705f, -0.754f };
   int ldc = 2;
   float C_expected[] = { 0.004695f, 0.050392f, -0.458321f, 1.42782f, 0.728f, 0.674f, 1.29896f, -1.54804f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1586) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1586) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 113;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.272f, -0.146f, 0.155f, 0.038f };
   int lda = 1;
   float C[] = { 0.533f, -0.41f, -0.904f, 0.301f, -0.836f, 0.57f, -0.374f, -0.293f };
   int ldc = 2;
   float C_expected[] = { 0.462668f, 0.453576f, -0.904f, 0.301f, -0.522292f, -0.848294f, 0.315581f, -0.36222f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1587) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1587) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { -0.055f, -0.127f, -0.896f, -0.625f };
   int lda = 1;
   float C[] = { -0.619f, 0.511f, -0.877f, 0.557f, -0.801f, -0.437f, -0.922f, 0.332f };
   int ldc = 2;
   float C_expected[] = { 0.60503f, -0.524104f, -0.877f, 0.557f, 0.652833f, 0.406905f, -0.198f, 0.080191f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1588) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1588) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { -0.528f, 0.759f, -0.079f, 0.952f };
   int lda = 2;
   float C[] = { 0.775f, 0.855f, 0.786f, 0.525f, 0.85f, 0.044f, 0.658f, 0.947f };
   int ldc = 2;
   float C_expected[] = { 0.026504f, -1.1523f, -0.223383f, -1.20586f, 0.85f, 0.044f, -0.507584f, -1.84706f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1589) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1589) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { -0.049f, -0.687f, -0.434f, 0.294f };
   int lda = 2;
   float C[] = { 0.937f, -0.113f, 0.796f, 0.293f, 0.876f, -0.199f, -0.757f, -0.103f };
   int ldc = 2;
   float C_expected[] = { 0.467432f, -0.045674f, 0.796f, 0.293f, 1.09924f, 0.084752f, -0.65508f, -0.358192f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1590) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1590) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { 0.359f, -0.364f, 0.926f, -0.69f };
   int lda = 1;
   float C[] = { 0.306f, 0.249f, 0.28f, 0.229f, 0.866f, 0.092f, 0.886f, -0.283f };
   int ldc = 2;
   float C_expected[] = { 0.302385f, -0.012352f, 0.361274f, -0.355774f, 0.866f, 0.092f, 1.26738f, -1.56088f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1591) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1591) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 113;
   int N = 2;
   int K = 1;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { 0.607f, 0.555f, -0.85f, 0.831f };
   int lda = 2;
   float C[] = { 0.069f, 0.368f, 0.551f, -0.912f, -0.243f, -0.063f, -0.924f, 0.192f };
   int ldc = 2;
   float C_expected[] = { -0.0855042f, -0.196089f, 0.551f, -0.912f, 0.28988f, -0.107516f, 0.131688f, 0.427004f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1592) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1592) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 113;
   int N = 2;
   int K = 1;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { 0.427f, 0.86f, -0.136f, 0.002f };
   int lda = 1;
   float C[] = { 0.398f, -0.47f, 0.011f, -0.547f, -0.106f, 0.016f, 0.681f, 0.246f };
   int ldc = 2;
   float C_expected[] = { 0.0937373f, -0.276059f, 0.0295482f, 0.0288526f, -0.106f, 0.016f, -0.0054932f, 0.0020124f };
   cblas_csyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyrk(case 1593) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyrk(case 1593) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.718, 0.023, 0.355, -0.492 };
   int lda = 1;
   double C[] = { -0.637, -0.727, -0.475, -0.776, 0.802, -0.55, -0.837, 0.222 };
   int ldc = 2;
   double C_expected[] = { -0.7948013, -0.6854089, -0.5203527, -0.6458521, 0.802, -0.55, -0.7672563, 0.3151921 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1594) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1594) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.209, 0.139, -0.202, -0.223 };
   int lda = 2;
   double C[] = { -0.695, 0.524, 0.212, -0.88, -0.752, 0.291, 0.684, -0.124 };
   int ldc = 2;
   double C_expected[] = { -0.7081182, 0.5090054, 0.212, -0.88, -0.7411652, 0.3122834, 0.6776683, -0.1519201 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1595) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1595) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.365, -0.624, 0.632, 0.348 };
   int lda = 2;
   double C[] = { 0.877, 0.927, -0.377, 0.967, 0.008, 0.292, -0.779, 0.794 };
   int ldc = 2;
   double C_expected[] = { 0.9082933, 0.7647289, -0.3208028, 1.1220636, 0.008, 0.292, -0.9064832, 0.6898704 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1596) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1596) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.067, -0.586, 0.208, 0.331 };
   int lda = 1;
   double C[] = { 0.584, -0.454, 0.93, 0.782, 0.489, -0.278, 0.081, -0.919 };
   int ldc = 2;
   double C_expected[] = { 0.6778197, -0.5114479, 0.93, 0.782, 0.4493975, -0.2167775, 0.0871195, -0.9669385 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1597) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1597) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 113;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { -0.617, 0.179, -0.626, 0.334 };
   int lda = 2;
   double C[] = { 0.346, -0.903, 0.022, -0.839, -0.715, 0.049, -0.338, 0.149 };
   int ldc = 2;
   double C_expected[] = { 0.1123886, 0.0694648, 0.1157132, 0.0348456, -0.715, 0.049, 0.0269168, -0.005768 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1598) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1598) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 113;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { -0.356, -0.308, 0.493, -0.351 };
   int lda = 1;
   double C[] = { -0.898, -0.905, 0.002, -0.219, 0.881, 0.879, 0.275, -0.351 };
   int ldc = 2;
   double C_expected[] = { 0.0685704, -0.0866128, 0.002, -0.219, -0.0852112, 0.0597384, 0.0697086, 0.0394848 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1599) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1599) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.103, -0.951, -0.601, -0.041 };
   int lda = 1;
   double C[] = { -0.918, -0.018, 0.991, -0.789, -0.698, -0.067, 0.956, -0.599 };
   int ldc = 2;
   double C_expected[] = { -0.9375906, -0.1073792, 0.991, -0.789, -0.7555774, -0.0647088, 0.9510718, -0.563048 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1600) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1600) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { -0.237, 0.925, -0.904, -0.091 };
   int lda = 2;
   double C[] = { -0.572, 0.915, 0.398, 0.222, 0.016, 0.288, -0.078, -0.507 };
   int ldc = 2;
   double C_expected[] = { -0.528155, 0.8350544, 0.4794633, 0.2518423, 0.016, 0.288, -0.0944528, -0.4261065 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1601) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1601) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0};
   double A[] = { 0.963, -0.23, -0.435, 0.289 };
   int lda = 2;
   double C[] = { 0.282, -0.272, -0.516, -0.594, -0.001, 0.155, -0.39, -0.354 };
   int ldc = 2;
   double C_expected[] = { -0.2180427, 0.2203409, -0.516, -0.594, 0.0678948, -0.1487506, -0.0065682, 0.0859994 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1602) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1602) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 2;
   int K = 1;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0};
   double A[] = { 0.674, 0.1, -0.098, 0.552 };
   int lda = 1;
   double C[] = { 0.089, -0.523, -0.551, 0.618, 0.67, 0.247, 0.975, -0.714 };
   int ldc = 2;
   double C_expected[] = { -0.1467628, 0.0039876, 0.0001508, -0.1207996, 0.67, 0.247, 0.0993492, 0.0029476 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1603) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1603) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 113;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { 0.033, -0.864, 0.168, 0.524 };
   int lda = 2;
   double C[] = { 0.788, 0.016, -0.436, 0.749, -0.89, -0.87, 0.421, -0.203 };
   int ldc = 2;
   double C_expected[] = { 0.845024, -0.729407, -0.436, 0.749, -0.76214, -0.41172, 0.244936, -0.449352 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1604) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1604) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 113;
   int N = 2;
   int K = 1;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { 0.957, -0.079, 0.935, 0.232 };
   int lda = 1;
   double C[] = { -0.744, -0.061, 0.195, -0.574, 0.551, 0.478, -0.337, 0.1 };
   int ldc = 2;
   double C_expected[] = { -0.592794, 0.848608, 0.046841, 0.339123, 0.551, 0.478, -0.77084, 0.920401 };
   cblas_zsyrk(order, uplo, trans, N, K, alpha, A, lda, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 4; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyrk(case 1605) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyrk(case 1605) imag");
     };
   };
  };

}


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

#ifndef KATLA_SOCKET_H
#define KATLA_SOCKET_H

#include "posix-errors.h"

#include "katla/core/error.h"
#include "katla/core/core.h"

#include <gsl/span>

#include <net/ethernet.h>

#include <memory>
#include <optional>
#include <system_error>

#include <chrono>
#include <string>

namespace outcome = OUTCOME_V2_NAMESPACE;

namespace katla {

class PosixSocket {
public:
    enum class ProtocolDomain { Unix, IPv4, IPv6, Packet, Can, Bluetooth, VSock};
    enum class Type { Stream, Datagram, SequencedPacket, Raw};
    enum class FrameType : uint16_t { All = ETH_P_ALL, EtherCat = 0x88a4 };

    struct SocketOptions {
        ProtocolDomain domain {ProtocolDomain::IPv4};
        Type type {Type::Stream};
        FrameType frameType {FrameType::All};
        bool nonBlocking {false};
        bool reuseAddress {false};
    };

    static constexpr SocketOptions IPv4SocketOptions {ProtocolDomain::IPv4, Type::Stream, FrameType::All, false, false};

    struct WaitResult {
        bool dataToRead {0};
        bool urgentDataToRead {0};
        bool writingWillNotBlock {0};
        bool readHangup {0};
        bool writeHangup {0};
        bool error {0};
        bool invalid {0};
        bool wakeup {0};
    };

    PosixSocket();
    PosixSocket(ProtocolDomain protocolDomain, Type type, FrameType frameType, bool nonBlocking);
    PosixSocket(const PosixSocket&) = delete;
    ~PosixSocket();

    static outcome::result<std::array<std::shared_ptr<PosixSocket>, 2>, Error> createUnnamedPair(ProtocolDomain protocolDomain, Type type, FrameType frameType, bool nonBlocking);

    outcome::result<void, Error> bind(std::string url);
    outcome::result<void, Error> bindIPv4(std::string ip, int port, SocketOptions options = IPv4SocketOptions);

    outcome::result<void, Error> listen();
    outcome::result<std::unique_ptr<PosixSocket>, Error> accept();
    std::optional<Error> error();

    outcome::result<void, Error> connect(std::string url, SocketOptions options);
    outcome::result<void, Error> connectIPv4(std::string ip, int port, SocketOptions options = IPv4SocketOptions);

    // TODO add wakeup to interrupt wait -> multithreading??
    outcome::result<WaitResult, Error> poll(std::chrono::milliseconds timeout, bool writePending = false);

    outcome::result<ssize_t, Error> read(const gsl::span<std::byte>& buffer);
    outcome::result<ssize_t, Error> write(const gsl::span<std::byte>& buffer);

    outcome::result<ssize_t, Error> receiveFrom(const gsl::span<std::byte>& buffer);
    outcome::result<ssize_t, Error> sendTo(std::string url, const gsl::span<std::byte>& buffer);

    outcome::result<void, Error> wakeup();

    outcome::result<void, Error> close();
private:
    outcome::result<void, Error> create();

    PosixSocket & operator=(const PosixSocket&) = delete;

    PosixSocket(ProtocolDomain protocolDomain, Type type, FrameType frameType, bool nonBlocking, int fd, int wakeupFd);

    static int mapProtocolDomain (ProtocolDomain protocolDomain);
    static int mapType (Type type);

    int _fd {-1};
    int _wakeupFd {-1};

    ProtocolDomain _protocolDomain {ProtocolDomain::IPv4};
    Type _type {Type::Stream};
    FrameType _frameType {FrameType::All};

    std::string _url;

    bool _nonBlocking {false};
};

}

#endif // KATLA_SOCKET_H

max_stars_repo_name: plok/katla

max_issues_repo_name: plok/katla

max_forks_repo_name: plok/katla


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

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>

#include "find.c"

int
gsl_histogram_increment (gsl_histogram * h, double x)
{
  int status = gsl_histogram_accumulate (h, x, 1.0);
  return status;
}

int
gsl_histogram_accumulate (gsl_histogram * h, double x, double weight)
{
  const size_t n = h->n;
  size_t index = 0;

  int status = find (h->n, h->range, x, &index);

  if (status)
    {
      return GSL_EDOM;
    }

  if (index >= n)
    {
      GSL_ERROR ("index lies outside valid range of 0 .. n - 1",
                 GSL_ESANITY);
    }

  h->bin[index] += weight;

  return GSL_SUCCESS;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE