text:
/**
 * @file xdegnome.c
 * @author Daniel R. Tabin
 * @brief Unit tests for Degnome
 */

#include "degnome.h"
#include <stdio.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <time.h>
#include <unistd.h>
#include <limits.h>

unsigned long rngseed=0;

int main(int argc, char **argv) {
	int verbose = 0;
	int seeded =0;

	if (argc == 2) {
		if (strncmp(argv[1], "-v", 2) != 0) {
			fprintf(stderr, "usage: xdegnome [-v] [-s] [0000000000]\n");
			exit(EXIT_FAILURE);
		}
		verbose = 1;
	}
	else if (argc == 3) {
		if (strncmp(argv[1], "-s", 2) != 0) {
			fprintf(stderr, "usage: xdegnome [-v] [-s] [0000000000]\n");
			exit(EXIT_FAILURE);
		}
		rngseed = (unsigned long) atoi(argv[2]);
	}
	else if (argc == 4) {
		if (strncmp(argv[1], "-v", 2) != 0) {
			fprintf(stderr, "usage: xdegnome [-v] [-s] [0000000000]\n");
			exit(EXIT_FAILURE);
		}
		else if (strncmp(argv[2], "-s", 2) != 0) {
			fprintf(stderr, "usage: xdegnome [-v] [-s]  [0000000000]\n");
			exit(EXIT_FAILURE);
		}
		verbose = 1;
		seeded = 1;
		rngseed = (unsigned long) atoi(argv[3]);

	}
	else if (argc != 1) {
		fprintf(stderr, "usage: xdegnome [-v] [-s] [seed]\n");
		exit(EXIT_FAILURE);
	}

	chrom_size = 15;

	if (verbose) {
		printf("Chromosome is length %u\n", chrom_size);
	}

	Degnome* bom_mom = Degnome_new(); // good parent
	Degnome* bad_dad = Degnome_new();	// bad parent
	Degnome* tst_bby = Degnome_new(); // their child

	// ============ change this to seed rn
	if (seeded == 0) {
		time_t currtime = time(NULL);                  // time
		unsigned long pid = (unsigned long) getpid();  // process id
		rngseed = currtime ^ pid;					   // random seed
	}
	gsl_rng* rng = gsl_rng_alloc(gsl_rng_taus);    // rand generator

	gsl_rng_set(rng, rngseed);
	rngseed = (rngseed == ULONG_MAX ? 0 : rngseed + 1);

	// ============ change this to seed rn

	for (int i = 0; i < chrom_size; i++) {
		bom_mom->dna_array[i] = 2*i;
		bad_dad->dna_array[i] = 1*i;

		bom_mom->hat_size += bom_mom->dna_array[i];
		bad_dad->hat_size += bad_dad->dna_array[i];
	}

	if (verbose) {
		printf("pre-mating values:\n");
		for (int i = 0; i < chrom_size; i++) {
			printf("Mom: %lf\t Dad: %f\n", bom_mom->dna_array[i], bad_dad->dna_array[i]);
		}
		printf("Mom hat_size: %lf\t Dad hat_size: %f\n", bom_mom->hat_size, bad_dad->hat_size);
	}

	Degnome_mate(tst_bby, bom_mom, bad_dad, rng, 0, 0,2);

	if (verbose) {
		printf("post-mating values:\n");
		for (int i = 0; i < chrom_size; i++) {
			printf("Mom: %lf\t Dad: %f\t Kid: %f\n", bom_mom->dna_array[i], bad_dad->dna_array[i], tst_bby->dna_array[i]);
		}
		printf("Mom hat_size: %lf\t Dad hat_size: %f\t Kid hat_size: %f\n", bom_mom->hat_size, bad_dad->hat_size, tst_bby->hat_size);
	}

	Degnome_free(bom_mom);
	Degnome_free(bad_dad);
	Degnome_free(tst_bby);

	gsl_rng_free(rng);

	printf("All tests for xdegnome completed\n");
}

max_stars_repo_name: masalemi/PopGenSim

max_issues_repo_name: masalemi/PopGenSim

max_forks_repo_name: masalemi/PopGenSim


text:
/*
 * qubit.h
 *
 *  Created on: 30 Jan 2019
 *      Author: jobinrjohnson
 */

#ifndef IMPLEMENTATION_QUBIT_H_
#define IMPLEMENTATION_QUBIT_H_

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_rng.h>
#include <iostream>

class Qubit
{

  protected:
	gsl_vector_complex *v_state;

  public:
	Qubit();
	Qubit(double, double, double, double);
	virtual ~Qubit();
	void print_state();
	gsl_vector_complex *get_state();
	bool set_state(gsl_vector_complex *);
};

#endif /* IMPLEMENTATION_QUBIT_H_ */

ext: h

max_stars_repo_name: jobinrjohnson/QSim

max_issues_repo_name: jobinrjohnson/QSim

max_forks_repo_name: jobinrjohnson/QSim


text:
/* Copyright (c) 2014, Giuseppe Argentieri <giuseppe.argentieri@ts.infn.it>

 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright notice,
 * this list of conditions and the following disclaimer in the documentation
 * and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 *
 */
/*
 *
 *
 *       Filename:  mat_file.c
 *
 *    Description:  Write and read the generator matrices on files
 *
 *        Version:  1.0
 *        Created:  03/05/2014 13:52:15
 *       Revision:  none
 *        License:  BSD
 *
 *         Author:  Giuseppe Argentieri (ga), giuseppe.argentieri@ts.infn.it
 *   Organization:  Università degli Studi di Trieste
 *
 *
 */

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_errno.h>
#include <string.h>
#include "funcs.h"

/*
 *      FUNCTION
 *         Name:  mat_write
 *  Description:  Save the generator matrix in a file.
 *
 */
int mat_write ( gsl_matrix* mat, char* name )
{
	FILE* f = fopen ( name, "w" ) ;
	if ( f == NULL )
		printf("Error: %s.\nFailed to open %s.\n", strerror(errno),
				name) ;

	int status = gsl_matrix_fprintf ( f, mat, "%.9f" ) ;

	fclose (f) ;

	return status;
}		/* -----  end of function mat_write  ----- */

/*
 *      FUNCTION
 *         Name:  mat_read
 *  Description:  Read the Redfield matrix from a file
 *
 */
int mat_read ( gsl_matrix* mat, char* name )
{
	FILE* f = fopen ( name, "r" ) ;
	if ( f == NULL )
		printf( "Error: %s.\nFailed to open %s.\n", strerror(errno),
				name ) ;

	int status = gsl_matrix_fscanf ( f, mat ) ;

	fclose (f) ;

	return status;
}		/* -----  end of function mat_read  ----- */

ext: c


text:
///
/// @file
///
/// @brief This file contains code that is used to validate the output of the
/// example codes.
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @section LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cblas.h>

static inline double squ(double x)
{
    return x*x;
}

void check_orthogonality(int n, size_t ldQ, double const *Q)
{
    size_t ldT = ((n/8)+1)*8;
    double *T = malloc(n*ldT*sizeof(double));

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0,
        Q, ldQ, Q, ldQ, 0.0, T, ldT);

    double dot = 0.0;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            dot += squ(T[i*ldT+j] - (i == j ? 1.0 : 0.0));

    double norm = ((long long)1<<52) * sqrt(dot)/sqrt(n);

    if (norm < 1000) {
        printf("The matrix is orthogonal.\n");
    }
    else {
        fprintf(stderr, "Matrix is not orthogonal.\n");
        exit(EXIT_FAILURE);
    }

    free(T);
}

void check_residual(
    int n, size_t ldQ, size_t ldA, size_t ldZ, size_t ldC,
    double const *Q, double const *A, double const *Z, double const *C)
{
    size_t ldT = ((n/8)+1)*8;
    double *T = malloc(n*ldT*sizeof(double));

    size_t ldY = ((n/8)+1)*8;
    double *Y = malloc(n*ldT*sizeof(double));

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0,
        Q, ldQ, A, ldA, 0.0, T, ldT);
    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0,
        T, ldT, Z, ldZ, 0.0, Y, ldY);

    double dot = 0.0, a_dot = 0.0;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            dot += squ(Y[i*ldY+j] - C[i*ldC+j]);
            a_dot += squ(C[i*ldC+j]);
        }
    }

    double norm = ((long long)1<<52) * sqrt(dot)/sqrt(a_dot);

    if (norm < 1000) {
        printf("The residual is small enough.\n");
    }
    else {
        fprintf(stderr, "The residual is too large.\n");
        exit(EXIT_FAILURE);
    }

    free(T);
    free(Y);
}

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

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

  gsl_matrix_complex * qr = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * a  = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * q  = gsl_matrix_complex_alloc(M, M);
  gsl_matrix_complex * r  = gsl_matrix_complex_alloc(M, N);
  gsl_vector_complex * d = gsl_vector_complex_alloc(N);

  gsl_matrix_complex_memcpy(qr, m);

  s += gsl_linalg_complex_QR_decomp(qr, d);
  s += gsl_linalg_complex_QR_unpack(qr, d, q, r);

  /* compute a = q r */
  gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, q, r, GSL_COMPLEX_ZERO, a);

  for (i = 0; i < M; i++) {
    for (j = 0; j < N; j++) {
      gsl_complex aij = gsl_matrix_complex_get(a, i, j);
      gsl_complex mij = gsl_matrix_complex_get(m, i, j);
      int foo_r = check(GSL_REAL(aij), GSL_REAL(mij), eps);
      int foo_i = check(GSL_IMAG(aij), GSL_IMAG(mij), eps);
      if (foo_r || foo_i)
        {
          printf("%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", desc, M, N, i, j, GSL_REAL(aij), GSL_REAL(mij));
          printf("%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", desc, M, N, i, j, GSL_IMAG(aij), GSL_IMAG(mij));
        }
      s += foo_r + foo_i;
    }
  }

  gsl_vector_complex_free(d);
  gsl_matrix_complex_free(qr);
  gsl_matrix_complex_free(a);
  gsl_matrix_complex_free(q);
  gsl_matrix_complex_free(r);

  return s;
}

static int
test_QRc_decomp(gsl_rng * r)
{
  int s = 0;
  size_t N, M;

  /* test random matrices */
  for (N = 1; N <= 20; ++N)
    {
      for (M = 1; M <= 20; ++M)
        {
          gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, M);

          create_random_complex_matrix(A, r);
          test_QRc_decomp_eps(A, 1.0e6 * GSL_MAX(N, M) * GSL_DBL_EPSILON, "complex_QR_decomp random");

          gsl_matrix_complex_free(A);
        }
    }

  return s;
}

static int
test_QRc_solve_eps(const gsl_matrix_complex * A, const gsl_vector_complex * rhs, const gsl_vector_complex * sol,
                   const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix_complex * QR = gsl_matrix_complex_alloc(M, N);
  gsl_vector_complex * tau = gsl_vector_complex_alloc(N);
  gsl_vector_complex * x  = gsl_vector_complex_alloc(N);

  gsl_matrix_complex_memcpy(QR, A);

  s += gsl_linalg_complex_QR_decomp(QR, tau);
  s += gsl_linalg_complex_QR_solve(QR, tau, rhs, x);

  for (i = 0; i < M; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

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

  gsl_matrix_complex_free(QR);
  gsl_vector_complex_free(tau);
  gsl_vector_complex_free(x);

  return s;
}

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

  for (N = 1; N <= 50; ++N)
    {
      gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
      gsl_vector_complex * sol = gsl_vector_complex_alloc(N);
      gsl_vector_complex * rhs = gsl_vector_complex_alloc(N);

      create_random_complex_matrix(A, r);
      create_random_complex_vector(sol, r);
      gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, A, sol, GSL_COMPLEX_ZERO, rhs);

      s += test_QRc_solve_eps(A, rhs, sol, 1.0e5 * N * GSL_DBL_EPSILON, "complex QR_solve_r random");

      gsl_matrix_complex_free(A);
      gsl_vector_complex_free(sol);
      gsl_vector_complex_free(rhs);
    }

  return s;
}

static int
test_QRc_lssolve_eps(const gsl_matrix_complex * A, const gsl_vector_complex * rhs, const gsl_vector_complex * sol,
                     const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix_complex * QR = gsl_matrix_complex_alloc(M, N);
  gsl_vector_complex * tau = gsl_vector_complex_alloc(N);
  gsl_vector_complex * x  = gsl_vector_complex_alloc(N);
  gsl_vector_complex * r  = gsl_vector_complex_alloc(M);
  gsl_vector_complex * r2  = gsl_vector_complex_alloc(M);

  gsl_matrix_complex_memcpy(QR, A);

  s += gsl_linalg_complex_QR_decomp(QR, tau);
  s += gsl_linalg_complex_QR_lssolve(QR, tau, rhs, x, r);

  for (i = 0; i < N; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

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

  /* compute residual and check */
  gsl_vector_complex_memcpy(r2, rhs);
  gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_NEGONE, A, sol, GSL_COMPLEX_ONE, r2);

  for (i = 0; i < M; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(r, i);
      gsl_complex yi = gsl_vector_complex_get(r2, i);

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

  gsl_matrix_complex_free(QR);
  gsl_vector_complex_free(tau);
  gsl_vector_complex_free(x);
  gsl_vector_complex_free(r);
  gsl_vector_complex_free(r2);

  return s;
}

static int
test_QRc_lssolve()
{
  int s = 0;
  const double tol = 1.0e-10;

  {
    const double A_data[] = { 6.70178096995592e-01, 2.51765675689489e-01, 7.18387884271362e-01, 3.61106008613644e-01,
                              3.36342428526987e-02, 5.25755806084206e-01, 5.87596968507183e-01, 5.27700895952418e-01,
                              6.85871226549483e-01, 3.71216051930735e-01, 4.40045953339814e-01, 2.08875308141557e-01,
                              8.62676230425306e-01, 5.79995188556082e-01, 4.86983443637474e-02, 2.82915411954634e-01,
                              2.92843706013013e-01, 5.61536446182319e-01, 6.85137614758495e-02, 8.90853425372411e-01,
                              4.38527971314087e-01, 4.78136089625096e-01, 1.57942824868494e-01, 8.38451530279972e-01,
                              6.36273325487226e-01, 7.74039464290391e-02, 5.45646256301364e-01, 7.80075219450629e-01,
                              9.27400956530167e-01, 7.01700239834713e-02, 1.09682812589509e-01, 5.64047584357803e-01,
                              4.46922541620762e-02, 3.03438549353353e-01, 8.09200219159660e-01, 1.44245237525133e-01 };
    const double rhs_data[] = { 5.82246876454474e-01, 1.42259458199622e-02, 6.25588177982770e-01, 3.79195077388159e-01,
                                8.45448918385455e-02, 3.20808711881935e-01, 8.02701461544476e-01, 5.65141425118050e-01,
                                8.34227120735637e-01, 4.69005326248388e-01, 9.73712086117648e-01, 3.47197650692321e-01 };
    const double sol_data[] = { -2.07103028285037e-01, -3.77392587461962e-01,
                                 7.46590544020302e-01, -1.10976592416587e-01,
                                 6.07653916072011e-01,  2.05471935567110e-01 };
    gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array(A_data, 6, 3);
    gsl_vector_complex_const_view rhs = gsl_vector_complex_const_view_array(rhs_data, 6);
    gsl_vector_complex_const_view sol = gsl_vector_complex_const_view_array(sol_data, 3);

    test_QRc_lssolve_eps(&A.matrix, &rhs.vector, &sol.vector, tol, "complex QR_lssolve test1");
  }

  return s;
}

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

  gsl_matrix_complex * QR = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * T = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * A  = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * R  = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * Q  = gsl_matrix_complex_alloc(M, M);
  gsl_matrix_complex_view Q1 = gsl_matrix_complex_submatrix(Q, 0, 0, M, N);
  gsl_vector_complex_const_view tau = gsl_matrix_complex_const_diagonal(T);

  gsl_matrix_complex_memcpy(QR, m);

  s += gsl_linalg_complex_QR_decomp_r(QR, T);
  s += gsl_linalg_complex_QR_unpack_r(QR, T, Q, R);

  /* compute A = Q R */
  gsl_matrix_complex_memcpy(A, &Q1.matrix);
  gsl_blas_ztrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, GSL_COMPLEX_ONE, R, A);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          gsl_complex aij = gsl_matrix_complex_get(A, i, j);
          gsl_complex mij = gsl_matrix_complex_get(m, i, j);

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

  if (M > N)
    {
      gsl_matrix_complex * R_alt  = gsl_matrix_complex_alloc(M, N);
      gsl_matrix_complex * Q_alt  = gsl_matrix_complex_alloc(M, M);

      /* test that Q2 was computed correctly by comparing with Level 2 algorithm */
      gsl_linalg_complex_QR_unpack(QR, &tau.vector, Q_alt, R_alt);

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < M; j++)
            {
              gsl_complex aij = gsl_matrix_complex_get(Q, i, j);
              gsl_complex bij = gsl_matrix_complex_get(Q_alt, i, j);

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

      gsl_matrix_complex_free(R_alt);
      gsl_matrix_complex_free(Q_alt);
    }

  gsl_matrix_complex_free(QR);
  gsl_matrix_complex_free(T);
  gsl_matrix_complex_free(A);
  gsl_matrix_complex_free(Q);
  gsl_matrix_complex_free(R);

  return s;
}

static int
test_QRc_decomp_r(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

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

          create_random_complex_matrix(A, r);
          s += test_QRc_decomp_r_eps(A, 1.0e6 * M * GSL_DBL_EPSILON, "complex_QR_decomp_r random");

          gsl_matrix_complex_free(A);
        }
    }

  return s;
}

static int
test_QRc_solve_r_eps(const gsl_matrix_complex * A, const gsl_vector_complex * rhs, const gsl_vector_complex * sol,
                     const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix_complex * QR = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * T = gsl_matrix_complex_alloc(N, N);
  gsl_vector_complex * x  = gsl_vector_complex_alloc(N);

  gsl_matrix_complex_memcpy(QR, A);

  s += gsl_linalg_complex_QR_decomp_r(QR, T);
  s += gsl_linalg_complex_QR_solve_r(QR, T, rhs, x);

  for (i = 0; i < M; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

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

  gsl_matrix_complex_free(QR);
  gsl_matrix_complex_free(T);
  gsl_vector_complex_free(x);

  return s;
}

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

  for (N = 1; N <= 50; ++N)
    {
      gsl_matrix_complex * A = gsl_matrix_complex_alloc(N, N);
      gsl_vector_complex * sol = gsl_vector_complex_alloc(N);
      gsl_vector_complex * rhs = gsl_vector_complex_alloc(N);

      create_random_complex_matrix(A, r);
      create_random_complex_vector(sol, r);
      gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, A, sol, GSL_COMPLEX_ZERO, rhs);

      s += test_QRc_solve_r_eps(A, rhs, sol, 1.0e6 * N * GSL_DBL_EPSILON, "complex_QR_solve_r random");

      gsl_matrix_complex_free(A);
      gsl_vector_complex_free(sol);
      gsl_vector_complex_free(rhs);
    }

  return s;
}

static int
test_QRc_lssolve_r_eps(const gsl_matrix_complex * A, const gsl_vector_complex * rhs, const gsl_vector_complex * sol,
                       const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix_complex * QR = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * T = gsl_matrix_complex_alloc(N, N);
  gsl_vector_complex * x  = gsl_vector_complex_alloc(M);
  gsl_vector_complex * r  = gsl_vector_complex_alloc(M);
  gsl_vector_complex * work = gsl_vector_complex_alloc(N);
  gsl_vector_complex_const_view x1 = gsl_vector_complex_const_subvector(x, N, M - N);
  double rnorm_expected, rnorm;

  gsl_matrix_complex_memcpy(QR, A);

  s += gsl_linalg_complex_QR_decomp_r(QR, T);
  s += gsl_linalg_complex_QR_lssolve_r(QR, T, rhs, x, work);

  for (i = 0; i < N; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

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

  /* compute residual and check */
  gsl_vector_complex_memcpy(r, rhs);
  gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_NEGONE, A, sol, GSL_COMPLEX_ONE, r);

  rnorm_expected = gsl_blas_dznrm2(r);
  rnorm = gsl_blas_dznrm2(&x1.vector);
  gsl_test_rel(rnorm, rnorm_expected, eps, "%s rnorm (%3lu): %22.18g   %22.18g\n",
               desc, N, rnorm, rnorm_expected);

  gsl_matrix_complex_free(QR);
  gsl_matrix_complex_free(T);
  gsl_vector_complex_free(x);
  gsl_vector_complex_free(r);
  gsl_vector_complex_free(work);

  return s;
}

static int
test_QRc_lssolve_r()
{
  int s = 0;
  const double tol = 1.0e-10;

  {
    const double A_data[] = { 6.70178096995592e-01, 2.51765675689489e-01, 7.18387884271362e-01, 3.61106008613644e-01,
                              3.36342428526987e-02, 5.25755806084206e-01, 5.87596968507183e-01, 5.27700895952418e-01,
                              6.85871226549483e-01, 3.71216051930735e-01, 4.40045953339814e-01, 2.08875308141557e-01,
                              8.62676230425306e-01, 5.79995188556082e-01, 4.86983443637474e-02, 2.82915411954634e-01,
                              2.92843706013013e-01, 5.61536446182319e-01, 6.85137614758495e-02, 8.90853425372411e-01,
                              4.38527971314087e-01, 4.78136089625096e-01, 1.57942824868494e-01, 8.38451530279972e-01,
                              6.36273325487226e-01, 7.74039464290391e-02, 5.45646256301364e-01, 7.80075219450629e-01,
                              9.27400956530167e-01, 7.01700239834713e-02, 1.09682812589509e-01, 5.64047584357803e-01,
                              4.46922541620762e-02, 3.03438549353353e-01, 8.09200219159660e-01, 1.44245237525133e-01 };
    const double rhs_data[] = { 5.82246876454474e-01, 1.42259458199622e-02, 6.25588177982770e-01, 3.79195077388159e-01,
                                8.45448918385455e-02, 3.20808711881935e-01, 8.02701461544476e-01, 5.65141425118050e-01,
                                8.34227120735637e-01, 4.69005326248388e-01, 9.73712086117648e-01, 3.47197650692321e-01 };
    const double sol_data[] = { -2.07103028285037e-01, -3.77392587461962e-01,
                                 7.46590544020302e-01, -1.10976592416587e-01,
                                 6.07653916072011e-01,  2.05471935567110e-01 };
    gsl_matrix_complex_const_view A = gsl_matrix_complex_const_view_array(A_data, 6, 3);
    gsl_vector_complex_const_view rhs = gsl_vector_complex_const_view_array(rhs_data, 6);
    gsl_vector_complex_const_view sol = gsl_vector_complex_const_view_array(sol_data, 3);

    test_QRc_lssolve_r_eps(&A.matrix, &rhs.vector, &sol.vector, tol, "complex_QR_lssolve_r test1");
  }

  return s;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* ============================================================ *
 * decomp_eb.c							*
 * Martin Kilbinger, Liping Fu 2008, 2009			*
 * ============================================================ */
#include "decomp_eb.h"
#include <gsl/gsl_sf_gamma.h>

/* === CFHTLS Wide 3rd data release, Fu&Kilbinger (2010) === */

/* S/N, Psi=19' eta=1/50 */
const double a_FK10_SN[N_FK10]        = {0.1197730890, -0.3881211865, 0.5212557875, -0.3440507036, 0.2761305382, -0.07286690971};

/* FoM, Psi=222', eta=1/10 */
const double a_FK10_FoM_eta10[N_FK10] = {0.009877788826, 0.1061397843, -0.4300211814, 0.5451016406, -0.3372272549, 0.1716983151};

/* FoM, Psi=222', eta=1/50 */
const double a_FK10_FoM_eta50[N_FK10] = {0.1239456383, -0.3881431858, 0.5579593467, -0.3679282338, 0.1540941993, 0.01293361618};

/* First-kind Chebyshev polynomial of order n */
#define EPS 1.0e-10
double Cheby(double x, int n, error **err)
{
   double Cn;

   testErrorRetVA(x<-1-EPS || x>1+EPS, mr_range, "x = %g out of range", *err, __LINE__, -1, x);

   if (x<-1) x = -1.0;
   if (x>1) x = 1.0;

   Cn = cos (n * acos (x));
   return Cn;
}
#undef EPS

/* Second-kind Chebyshev polynomial of order n */
#define EPS 1.0e-10
double Cheby2(double x, int n, error **err)
{
   double Un;

   testErrorRetVA(x<-1-EPS || x>1+EPS, mr_range, "x = %g out of range", *err, __LINE__, -1, x);

   if (x<-1) x = -1.0;
   if (x>1) x = 1.0;

   if (x == 1) {
      Un = n + 1;
   } else if (x == -1) {
      Un = pow(-1.0, n) * (n + 1.0);
   } else {
      Un = sin((n+1.0)*acos(x))/sin(acos (x));
   }

   return Un;
}
#undef EPS

/* Legendre polynomial of order n */
double Legen(double x, int n)
{

   switch (n) {
      case 0:
	 return 1.0;
      case 1:
	 return x;
      default:
	 return (2.0*n-1.0)/(double)n*x*Legen(x,n-1)-(n-1.0)/(double)n*Legen(x,n-2);
   }
}

/* General basis function of order n */
#define EPS 1.0e-10
double C(double x, int n, poly_t poly, error **err)
{
   double c;

   testErrorRetVA(x<-1-EPS || x>1+EPS, mr_range, "x = %g out of range", *err, __LINE__, -1, x);

   if (x<-1) x = -1.0;
   if (x>1) x = 1.0;

   switch (poly) {
      case cheby :
	 c = Cheby(x, n, err);
	 forwardError(*err, __LINE__, -1.0);
	 break;
      case cheby2 :
	 c = Cheby2(x, n, err);
	 forwardError(*err, __LINE__, -1.0);
	 break;
      case legen :
	 c = Legen(x, n);
	 break;
      default :
	 *err = addErrorVA(mr_poly, "Unknown polynomial type %d", *err, __LINE__, poly);
	 return -1.0;
   }

   return c;
}
#undef EPS

/* ============================================================ *
 * The filter function for the generalised ring statistics.	*
 * FK09 (11).							*
 * ============================================================ */
#define	 EPS 1.0e-6
double Tp(double x, const double *a, int N, poly_t poly, error **err)
{
   int n;
   double res, Cn, Cnm1=0.0, Cnm2;

   if (x<-1-EPS || x>+1+EPS) return 0;

   /* NEW! */
   //if (x<-1) x = -1;
   //if (x>+1) x = +1;

   testErrorRetVA(N<=0, mr_range, "N has to be larger than zero but is %d", *err, __LINE__, 0.0, N);

   testErrorRetVA(x>+1+EPS, mr_range, "x=%.10f out of range", *err, __LINE__, 0.0, x);
   testErrorRetVA(x<-1-EPS, mr_range, "x=%.10f out of range", *err, __LINE__, 0.0, x);

   if (poly==cheby || poly==cheby2) {

      Cnm2 = 1.0;
      if (poly==cheby)       Cnm1 = x;
      else if (poly==cheby2) Cnm1 = 2.0*x;

      res  = 0.0;
      res += a[0] * Cnm2;
      if (N==1) return res;
      res += a[1] * Cnm1;
      if (N==2) return res;

      for (n=2; n<N; n++) {
	 /* Chebyshev (both T,U) recurrence relation *
	  * C_n(x) = 2 x C_{n-1}(x) - C_{n-2}(x)     */
	 Cn   = 2.0*x*Cnm1 - Cnm2;
	 res += a[n] * Cn;
	 Cnm2 = Cnm1;
	 Cnm1 = Cn;
      }

   } else {

      for (n=0,res=0.0; n<N; n++) {
	 Cn   = C(x, n, poly, err);
	 forwardError(*err, __LINE__, -1);
	 res += a[n] * Cn;
      }

   }

   return res;
}
#undef EPS

/* ============================================================ *
 * FK09 (24).							*
 * ============================================================ */
double Fn0(double x, int n, poly_t poly, error **err)
{
   double fn0, denom=0.0;
   int nn;

   if (n<0) nn = -n;
   else nn = n;

   if (poly==cheby) {
      if (nn==0) return x + 1.0;
      if (nn==1) return 0.5*(x*x - 1.0);
   } else if (poly==cheby2) {
      if (n==-1) return 0.0;  /* n, not nn: no symmetry for cheby2! */
   } else {
      *err = addErrorVA(mr_poly, "Tm for polynomial type %d not defined", *err, __LINE__, poly);
      return 0.0;
   }

   fn0 = 0.0;

   if (nn%2==0) fn0 += 1.0;   /* n even */
   else fn0 -= 1.0;           /* n odd  */

   //printf("# ** 1 %g\n", fn0);

   fn0 += x*Cheby(x, nn, err);
   forwardError(*err, __LINE__, 0.0);

   //printf("# ** 2 %g\n", fn0);

   if (poly==cheby) fn0 += n*(1 - x*x)*Cheby2(x, n-1, err);
   else if (poly==cheby2) fn0 -= (1 - x*x)*Cheby2(x, n-1, err);
   forwardError(*err, __LINE__, 0.0);

   //printf("# ** 3 %g\n", fn0);

   if (poly==cheby) denom = 1.0 - (double)(n*n);
   else if (poly==cheby2) denom = 1.0+(double)n;         /* n, not nn! */

   fn0 /= denom;

   return fn0;
}

/* ============================================================ *
 * FK09 (26)							*
 * ============================================================ */
void Fnnu(double x, int n, poly_t poly, double Fn[], error **err)
{
   double f0, fp1, fp2, fp3, fm1, fm2, fm3;

   f0  = Fn0(x, n, poly, err);      forwardError(*err, __LINE__,);
   fp1 = Fn0(x, n+1, poly, err);    forwardError(*err, __LINE__,);
   fp2 = Fn0(x, n+2, poly, err);    forwardError(*err, __LINE__,);
   fp3 = Fn0(x, n+3, poly, err);    forwardError(*err, __LINE__,);
   fm1 = Fn0(x, n-1, poly, err);    forwardError(*err, __LINE__,);
   fm2 = Fn0(x, n-2, poly, err);    forwardError(*err, __LINE__,);
   fm3 = Fn0(x, n-3, poly, err);    forwardError(*err, __LINE__,);

   Fn[0] = f0;
   Fn[1] = 0.5*(fp1 + fm1);
   Fn[2] = 0.25*(fp2 + 2.0*f0 + fm2);
   Fn[3] = 0.125*(fp3 + 3.0*fp1 + 3.0*fm1 + fm3);
}

/* ============================================================ *
 * FK09 (23)							*
 * ============================================================ */
#define EPS 1.0e-10
double alpha_explicit(double x, int n, double R, poly_t poly, error **err)
{
   double r, Fn[4], alpha;

   testErrorRet(poly!=cheby && poly!=cheby2, mr_poly,
		"Explicit expression only obtained so far for cheby1 and cheby2",
		*err, __LINE__, 0.0);

   if (x<-1.0+EPS) return 0.0;
   assert(x>-1.0);

   r = 3.0/dsqr(x + R);
   Fnnu(x, n, poly, Fn, err);
   forwardError(*err, __LINE__, 0.0);

   alpha  = R*(1.0-r*R*R)*Fn[0];
   alpha += (1.0 - 3.0*r*R*R)*Fn[1];
   alpha += -3.0*r*R*Fn[2];
   alpha += -r*Fn[3];
   alpha *= 4.0/3.0*r;

   return alpha;
}

/* ============================================================ *
 * The filter function T_- obtained from T_+, FK09 (20).	*
 * ============================================================ */
#define EPS 1.0e-10
double Tm(double x, const double *a, int N, poly_t poly, double R, error **err)
{
   double tm, c, alph;
   int n;

   testErrorRetVA(x<-1-EPS, mr_range, "x=%.20f smaller than -1", *err, __LINE__, 0.0, x);

   for (n=0,tm=0.0; n<N; n++) {
      c    = C(x, n, poly, err);                forwardError(*err, __LINE__, 0.0);
      alph = alpha_explicit(x, n, R, poly, err);
      forwardError(*err, __LINE__, 0.0);

      tm  += a[n]*(c + alph);
   }
   return tm;
}
#undef EPS

/* ================================================================== *
 * Returns the generalised ring statistic using a as the coefficients *
 * of the T_+ decomposition, see FK09 (4, 8, 11).		      *
 * The correlation functions xi+ and xi- are read from the            *
 * Nxi-dimensional arrays xip and xim, on angular scales th.          *
 * N is the size of the coefficient-array a.                          *
 * ================================================================== */
double RR_data(const double *xip, const double *xim, const double *th, const int Nxi,
	       double THETA_MIN, double THETA_MAX, const double *a, int N,
	       poly_t poly,  int pm, error **err)
{
   double theta, dtheta, res, A, B, x, summand;
   double dt1, dt2;
   int i;

   testErrorRetVA(abs(pm)!=1, mr_incompatible, "pm=%d not valid, has to be +1 or -1",
		  *err, __LINE__, 0.0, pm);

   testErrorRetVA(THETA_MIN<th[0], mr_range,
		  "THETA_MIN=%g' is smaller than minimum angular scale for xi+-, %g'",
		  *err, __LINE__, 0.0, THETA_MIN/arcmin, th[0]/arcmin);

   testErrorRetVA(THETA_MAX>th[Nxi-1], mr_range,
		  "THETA_MAX=%g' is larger than maximum angular scale for xi+-, %g'",
		  *err, __LINE__, 0.0, THETA_MAX/arcmin, th[Nxi-1]/arcmin);

   if (a==NULL) {
      *err = addError(mr_null, "Coefficients a=NULL", *err, __LINE__);
      return 0.0;
   }

   /* FK09 (8) */
   A = (THETA_MAX - THETA_MIN)/2.0;
   B = (THETA_MAX + THETA_MIN)/2.0;

   for (i=0,res=0.0; i<Nxi; i++) {
      theta  = th[i];
      if (theta < THETA_MIN) continue;
      else if (theta > THETA_MAX) break;

      /* theta[i] is the bin center */

      if (i==0) {
	 dt1 = (th[i+1] - th[i])/2.;
	 dtheta = 2*dt1;
      } else if (i==Nxi-1) {
	 dt2 = (th[i] - th[i-1])/2.;
	 dtheta = 2*dt2;
      } else {
	 dt1 = (th[i+1] - th[i])/2.;
	 dt2 = (th[i] - th[i-1])/2.;
	 dtheta = dt2 + dt1;
      }
      x = (theta-B)/A;

      if (pm==+1)  {
	 summand  = xip[i]*theta/dsqr(THETA_MAX);
	 summand *= Tp(x, a, N, poly, err);
	 forwardError(*err, __LINE__, -1.0);
      } else if (pm==-1) {
	 summand  = xim[i]*theta/dsqr(THETA_MAX);
	 summand *= Tm(x, a, N, poly, B/A, err);
	 forwardError(*err, __LINE__, -1.0);
      } else {
	 summand = 0.0;
      }

      res += summand*dtheta;
   }

   testErrorRet(!finite(res), math_infnan, "R is not finite", *err, __LINE__,
		  0.0);

   return res;
}

/* If cov_mode=outer, n, m are not used. For cov_mode=inner, a, N are not used.  *
 * If cov_mode=fromZ, the covariance using Z+ is returned and a, N are not used. */
#define EPS 1.0e-5
double cov_RR(const double *THETA_MIN, const double *THETA_MAX, const double *a, int N, poly_t poly,
	      const double *theta, const double *cov_xi, int Ntheta,
	      cov_mode_t cov_mode, int n, int m, double fac, error **err)
{
   double sum, dlogtheta, tmp, A[2], B[2], xi, xj, yi, yj, thetai, thetaj;
   int i, j;
   static double  **Cov_xi = NULL;

   /* Testing for a==NULL not possible because cov_RR is also called from get_cov_RR_inner_array */
   if (cov_mode==fromZ) {
#ifdef __MRING_H
      sum = cov_RR_Z(THETA_MIN, THETA_MAX, theta, cov_xi, Ntheta, err);
      forwardError(*err, __LINE__, 0.0);
      return sum;
#else
      *err = addError(mr_type, "cov_mode 'fromZ' not supported here", *err, __LINE__);
      return 0.0;
#endif
   }

   for (i=0; i<2; i++) {
      A[i] = (THETA_MAX[i] - THETA_MIN[i])/2.0;
      B[i] = (THETA_MAX[i] + THETA_MIN[i])/2.0;
   }

   dlogtheta = log(theta[1]) - log(theta[0]);

   for (i=0; i<2; i++) {
      testErrorRetVA(theta[0]>THETA_MIN[i], mr_range,
		     "Minumum of xi-covariance (%g=%g') larger than THETA_MIN (%g=%g')",
		     *err, __LINE__, -1, theta[0], theta[0]/arcmin, THETA_MIN[i], THETA_MIN[i]/arcmin);
      testErrorRetVA(theta[Ntheta-1]<THETA_MAX[i], mr_range,
		     "Maximum of xi-covariance scale (%g=%g') smaller than THETA_MAX (%g=%g')",
		     *err, __LINE__, -1, theta[Ntheta-1], theta[Ntheta-1]/arcmin, THETA_MAX[i],
		     THETA_MAX[i]/arcmin);
   }

   if (Cov_xi==NULL) {
      Cov_xi = sm2_matrix(0, Ntheta-1, 0, Ntheta-1, err);
      forwardError(*err, __LINE__, -1.0);
      for (i=0; i<Ntheta; i++) {
	 for (j=0; j<Ntheta; j++) {

	    Cov_xi[i][j] = cov_xi[i*Ntheta+j];

	 }

      }

      //if (fac<1.12)
	//fprintf(stderr, "WARNING: Cov-xi-factor = %g, too small for stable inversion of covariance!!!\n", fac);
   }

   sum = 0.0;
   //for (i=0; i<Ntheta; i++) {
   //thetai = theta[i];
   for (thetai=theta[0]; thetai<=theta[Ntheta-1]; thetai*=fac) {
      if (thetai<THETA_MIN[0] || thetai>THETA_MAX[0]) continue;
      xi = (thetai-B[0])/A[0];
      yi = thetai/THETA_MAX[0];

      //for (j=0; j<Ntheta; j++) {
      //thetaj = theta[j];
      for (thetaj=theta[0]; thetaj<=theta[Ntheta-1]; thetaj*=fac) {
	 if (thetaj<THETA_MIN[1] || thetaj>THETA_MAX[1]) continue;
	 xj = (thetaj-B[1])/A[1];
	 yj = thetaj/THETA_MAX[1];

	 tmp  = dsqr(dlogtheta);

	 if (cov_mode==outer) {
	    tmp *= Tp(xi, a, N, poly, err);
	    forwardError(*err, __LINE__, -1);
	    tmp *= Tp(xj, a, N, poly, err);
	    forwardError(*err, __LINE__, -1);
	 } else {
	    tmp *= C(xi, n, poly, err);  	    forwardError(*err, __LINE__, -1);
	    tmp *= C(xj, m, poly, err);  	    forwardError(*err, __LINE__, -1);
	 }

	 tmp *= dsqr(yi*yj);

	 //tmp *= cov_xi[i*Ntheta+j];
	 tmp *= sm2_interpol2d(Cov_xi, Ntheta, log(theta[0]), log(theta[Ntheta-1]), dlogtheta, log(thetai),
	 		       Ntheta, log(theta[0]), log(theta[Ntheta-1]), dlogtheta, log(thetaj), 0.0, 0.0, err);
	 forwardError(*err, __LINE__, -1.0);

	 sum += tmp;

      }

  }

   /* Comment the following two lines for faster code, but be aware that no check is done
    * whether the input covariance changes! */
   sm2_free_matrix(Cov_xi, 0, Ntheta-1, 0, Ntheta-1);
   Cov_xi = NULL;

   return sum;
}
#undef EPS

/* RR-Covariance using a diagonal xi-covariance.				 *
 * If cov_mode=outer, n, m are not used. For cov_mode=inner, a, N are not used.  *
 * If cov_mode=fromZ, the covariance using Z+ is returned and a, N are not used. */
#define EPS 1.0e-5
double cov_RR_diag_xi(const double *THETA_MIN, const double *THETA_MAX, const double *a, int N, poly_t poly,
		      const double *theta, const double *var_xi, int Ntheta,
		      int islog, error **err)
{
   double sum, dlogtheta, dtheta, tmp, A[2], B[2], xi, xj, yi, yj, thetai, thetaj;
   int i;

   for (i=0; i<2; i++) {
      A[i] = (THETA_MAX[i] - THETA_MIN[i])/2.0;
      B[i] = (THETA_MAX[i] + THETA_MIN[i])/2.0;
   }

   if (islog==0) {
      dtheta = theta[1] - theta[0];
      dlogtheta = 0.0;
   } else if (islog==1) {
      dlogtheta = log(theta[1]) - log(theta[0]);
      dtheta = 0.0;
   } else {
      *err = addErrorVA(mr_type, "Invalid flag islog = %d\n", *err, __LINE__, islog);
      return 0.0;
   }

   for (i=0; i<2; i++) {
      testErrorRetVA(theta[0]>THETA_MIN[i], mr_range,
		     "Minumum of xi-covariance (%g=%g') larger than THETA_MIN (%g=%g')",
		     *err, __LINE__, -1, theta[0], theta[0]/arcmin, THETA_MIN[i], THETA_MIN[i]/arcmin);
      testErrorRetVA(theta[Ntheta-1]<THETA_MAX[i], mr_range,
		     "Maximum of xi-covariance scale (%g=%g') smaller than THETA_MAX (%g=%g')",
		     *err, __LINE__, -1, theta[Ntheta-1], theta[Ntheta-1]/arcmin, THETA_MAX[i],
		     THETA_MAX[i]/arcmin);
   }

   sum = 0.0;
   for (i=0; i<Ntheta; i++) {
      thetai = theta[i];
      thetaj = thetai;

      if (thetai<THETA_MIN[0] || thetai>THETA_MAX[0]) continue;
      xi = (thetai-B[0])/A[0];
      yi = thetai/THETA_MAX[0];

      xj = (thetaj-B[1])/A[1];
      yj = thetaj/THETA_MAX[1];

      if (islog==0) {
	 //	 tmp = dsqr(dtheta) /thetai/thetai;
	 tmp = dsqr(dtheta)*yi*yj/THETA_MAX[0]/THETA_MAX[1];

      } else {

	 tmp  = dsqr(dlogtheta);
	 tmp *= dsqr(yi*yj);
      }

      tmp *= Tp(xi, a, N, poly, err);
      forwardError(*err, __LINE__, -1);
      tmp *= Tp(xj, a, N, poly, err);
      forwardError(*err, __LINE__, -1);

      tmp *= var_xi[i];
      //      tmp *= dsqr(yi*yj);

      sum += tmp;
  }

   return sum;
}
#undef EPS

/* ============================================================ *
 * Returns the chi^2 for a null test (RR=0).			*
 * ============================================================ */
double chi2_RB_null(const double *RB, const double *covRB, int NRB)
{
   int i, j;
   double c2;

   for (i=0,c2=0.0; i<NRB; i++) {
      for (j=0; j<NRB; j++) {
	 c2 += RB[i] * covRB[i*NRB+j] * RB[j];
      }
   }

   return c2;
}

/* ============================================================ *
 * Reads COSEBIs zeros and normalisation coefficients from file *
 * and returns polynomial coefficients.	File name is automatic, *
 * scale checks are done.					*
 * ============================================================ */
double *read_zeros_norm_cosebi_auto_check(double Psimin, double Psimax, const char *path, error **err)
{
   double *c_cosebi, psimin, psimax;
   char rname[1024];

   sprintf(rname, "%s/cosebi_tplog_rN_%d_%.3f_%.3f",
	   path == NULL ? "." : path, NMAX_COSEBI, Psimin/arcmin, Psimax/arcmin);
   c_cosebi = read_zeros_norm_cosebi(rname, &psimin, &psimax, err);
   forwardError(*err, __LINE__, NULL);

   testErrorRetVA(fabs(psimin-Psimin/arcmin)>EPSILON, mr_file,
		  "Inconsistent file %s, psimin=%g should be %g according to the file name",
		  *err, __LINE__, NULL, rname, psimin, Psimin/arcmin);
   testErrorRetVA(fabs(psimax-Psimax/arcmin)>EPSILON, mr_file,
		  "Inconsistent file %s, psimax=%g should be %g according to the file name",
		  *err, __LINE__, NULL, rname, psimax, Psimax/arcmin);

   return c_cosebi;
}

/* ============================================================ *
 * Reads COSEBIs zeros and normalisation coefficients from file *
 * and sets corresponding min and max scales.			*
 * Calculates polynomial coefficients c and returns them.       *
 * ============================================================ */
double *read_zeros_norm_cosebi(const char *rname, double *psimin, double *psimax, error **err)
{
   FILE *F;
   int Nzeros, Ncoeff, k, n, off_c, off_R, nmax;
   ssize_t nread;
   double *Rn, *Norm, *c;
   char *str, *line=NULL;

   F = fopen_err(rname, "r", err);                  forwardError(*err, __LINE__, NULL);

   /* Read two header lines */
   line = malloc_err(1024*sizeof(char), err);       forwardError(*err, __LINE__, NULL);
   str = fgets(line, 1024, F);
   str = fgets(line, 1024, F);
   free(line);

   nread = fscanf(F, "%d  %lg %lg\n", &nmax, psimin, psimax);
   testErrorRet(nread != 3, mr_file, "File has wrong format.", *err, __LINE__, NULL);

   testErrorRetVA(nmax > NMAX_COSEBI, mr_range,
		  "COSEBI number of modes n=%d read from file %s cannot be larger than NMAX_COSEBI=%d",
		  *err, __LINE__, NULL, nmax, rname, NMAX_COSEBI);

   /* Number of zeros for nmax polynomials,               *
    * sum_{i=1}^{nmax+1} (i+1) = sum_{i=0}^{nmax+1} i - 1 */
   Nzeros = (nmax+1) * (nmax+2) / 2 - 1;

   //fprintf(stderr, "nmax = %d , Psi = [%g, %g]\n", nmax, *psimin, *psimax);
   //fprintf(stderr, "Reading %d Rn, %d Norm, from %s\n", Nzeros, nmax, rname);

   Rn   = malloc_err(sizeof(double)*Nzeros, err);   forwardError(*err, __LINE__, NULL);
   Norm = malloc_err(sizeof(double)*nmax, err);     forwardError(*err, __LINE__, NULL);

   /* Read zeros */
   for (k=0; k<Nzeros; k++) {
      fscanf(F, "%lg ", Rn+k);
   }
   /* Read normalisations */
   for (k=0; k<nmax; k++) {
      fscanf(F, "%lg ", Norm+k);
   }

   fclose(F);

   /* Calculate polynomial coefficients */
   Ncoeff = NMAX_COSEBI * (NMAX_COSEBI + 5) / 2;
   c =  malloc_err(sizeof(double)*Ncoeff, err);        forwardError(*err, __LINE__, NULL);

   /* Polynomial T_i of order i+1 */
   for (n=1,off_R=0; n<=nmax; n++) {

      off_c = n*(n + 3)/2 - 2;   // Offset for coefficients for polynomial #n

      for (k=0; k<n+2; k++) {
	 c[k + off_c] = sum_combinations(n+1-k, n+1, Rn+off_R, err);
	 forwardError(*err, __LINE__, NULL);
	 c[k + off_c] *= Norm[n-1];
	 if (k % 2 == 1) c[k + off_c] = -c[k + off_c];
	 if (n % 2 == 0) c[k + off_c] = -c[k + off_c];
      }

      off_R = (n+1)*(n+2)/2 - 1; // Ncoeff for index i, degree i+1
   }

   free(Rn);
   free(Norm);

   return c;
}

/* ============================================================ *
 * COSEBI logarithmis filter function, SEK10 (28). The          *
 * coefficients have to be precalculated from the zeros, output *
 * by Mathematica.						*
 * Not normalised. 						*
 * ============================================================ */
double Tplog_c(double z, const double *c, int n, error **err)
{
   int i, off_c;
   double res, zpowerm;
   const double *cn;

   off_c = n*(n + 3)/2 - 2; // Offset for coefficients for polynomial #n
   cn    = c + off_c;

   for (i=0,res=0.0,zpowerm=1.0; i<=n+1; i++) {
      res     += cn[i] * zpowerm;
      zpowerm *= z;
   }

   testErrorRet(!finite(res), math_infnan, "Tplog is not finite", *err, __LINE__, 0.0);

   return res;
}

/* ============================================================ *
 * Filter function, SEK10 (38, 39). Normalised.                 *
 * ============================================================ */
double Tmlog(double z, const double *c, int n, error **err)
{
   int m, off_c;
   double res, zpowerm;
   const double *cn;

   off_c = n*(n + 3)/2 - 2; // Offset for coefficients for polynomial #n
   cn    = c + off_c;

   for (m=0,res=0.0,zpowerm=1; m<=n; m++) { // n+1 MKDEBUG ??
      res += dnm(n, m ,cn) * zpowerm;
      zpowerm *= z;
   }

   res += an2(n, cn) * exp(-2.0*z) - an4(n, cn) * exp(-4.0*z);
   testErrorRet(!finite(res), math_infnan, "Tmlog is not finite", *err, __LINE__, 0.0);

   return res;
}

double an2(int n, const double *c)
{
   int j;
   double res, jfac, m2powerjp1;

   for (j=0,res=0.0,jfac=1.0,m2powerjp1=-2.0; j<=n+1; j++) {
      res        += c[j]*jfac / m2powerjp1;
      jfac       *= j + 1.0;      /* factorial(j) */
      m2powerjp1 *= -2.0;        /* (-2)^(j+1)   */
   }

   res *= 4.0;

   return res;
}

double an4(int n, const double *c)
{
   int j;
   double res, jfac, m4powerjp1;

   for (j=0,res=0.0,jfac=1.0,m4powerjp1=-4.0; j<=n+1; j++) {
      res        += c[j]*jfac / m4powerjp1;
      jfac       *= j + 1.0;
      m4powerjp1 *= -4.0;        /* (-4)^(j+1) */
   }

   res *= 12.0;

   return res;
}

/* dnm(n+1, n, c) = 0 */
double dnm(int n, const int m, const double *c)
{
   int j;
   double res, jfac, m2powermjm1, p2powermjm1, mfac;

   p2powermjm1 = 0.5;
   m2powermjm1 = -0.5;
   mfac = gsl_sf_fact(m);
   for (j=m,res=0.0,jfac=mfac; j<=n+1; j++) {
      res         += c[j] * jfac * m2powermjm1 * (3.0 * p2powermjm1 - 1.0);
      jfac        *= j + 1.0;        /* j! */
      m2powermjm1 /= -2.0;    /* (-2)^(m-j-1) */
      p2powermjm1 /= 2.0;     /* 2^(m-k-1)    */
   }

   res *= 4.0/mfac;
   res += c[m];

   return res;
}

double sum_combinations(int j, int n, const double *r, error **err)
{
   gsl_combination *c;
   size_t i, c_i;
   double sum, prod;

   if (n==0) return 1.0;

   c = gsl_combination_calloc(n, j);

   sum = 0.0;
   do {
      for (i=0,prod=1.0; i<j; i++) {
	 c_i =  gsl_combination_get(c, i);
	 prod *= r[c_i];
	 //printf(" %zu", c_i);
      }
      sum += prod;
      //printf("{"); gsl_combination_fprintf(stdout, c, " %u"); printf (" }\n");
   } while (gsl_combination_next(c) == GSL_SUCCESS);

   gsl_combination_free(c);

   return sum;
}

/* ============================================================ *
 * E-/B-mode correlation functions from COSEBIs, SEK10 (40).    *
 * The results are stored in xi_pm_EB = (E+, E-, B+, B-).       *
 * ============================================================ */
void xipmEB(double theta, double THETA_MIN, double THETA_MAX, const double *c,
	    const double *E, const double *B, int N, double xi_pm_EB[4], error **err)
{
   double Tp, Tm, z;
   int n;

   for (n=0; n<4; n++) {
      xi_pm_EB[n] = 0.0;
   }

   for (n=1; n<=N; n++) {
      z  = log(theta/THETA_MIN);
      Tp = Tplog_c(z, c, n, err);
      Tm = Tmlog(z, c, n, err);

      xi_pm_EB[0] += E[n] * Tp;   // E+
      xi_pm_EB[1] += E[n] * Tm;   // E-
      xi_pm_EB[2] += B[n] * Tp;   // B+
      xi_pm_EB[3] += B[n] * Tm;   // B-
   }

   for (n=0; n<4; n++) {
      xi_pm_EB[n] *= 2.0 / theta / (THETA_MAX - THETA_MIN);
   }
}

max_stars_repo_name: danielgruen/ccv

max_issues_repo_name: danielgruen/ccv

max_forks_repo_name: danielgruen/ccv


text:
#include <errno.h>
#include <getopt.h>
#include <gsl/gsl_spmatrix.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/stat.h>

static int clock = 1;
static int results = 0;
static int verbose = 0;
static int help = 0;

int test (size_t m,
          size_t n,
          size_t nnz,
          const size_t *ptr,
          const size_t *ind,
          size_t B,
          double epsilon,
          double delta,
          int trials,
          int clock,
          int results,
          int verbose);

char *name ();

static void usage () {
  fprintf(stderr,"usage: %s [options] <input>\n"
  "  <input>                    MatrixMarket file (estimate fill of this matrix)\n"
  "  -B, --max-block-size <arg> Maximum block dimension for fill estimates\n"
  "  -e, --epsilon <arg>        Be accurate to relative error epsilon\n"
  "  -d, --delta <arg>          With probability (1 - delta)\n"
  "  -t, --trials <arg>         Number of trials to run\n"
  "  -c, --clock                Display timing information\n"
  "  -C, --noclock              Do not display timing information\n"
  "  -r, --results              Display fill estimates for all trials\n"
  "  -R, --noresults            Do not display fill estimates\n"
  "  -v, --verbose              Verbose mode\n"
  "  -q, --quiet                Quiet mode\n"
  "  -h, --help                 Display help message\n", name());
}

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

  size_t B = 12;
  double epsilon = 0.1;
  double delta = 0.01;
  int trials = 1;

  /* Beware. Option parsing below. */
  long longarg;
  double doublearg;
  while (1) {
    static char *options = "B:e:d:t:cCrRvqh";
    static struct option long_options[] = {
        {"max-block-size", required_argument, 0, 'B'},
        {"trials",         required_argument, 0, 't'},
        {"epsilon", required_argument, 0, 'e'},
        {"delta", required_argument, 0, 'd'},
        {"clock",     no_argument, &clock,   1},
        {"noclock",   no_argument, &clock,   0},
        {"results",   no_argument, &results, 1},
        {"noresults", no_argument, &results, 0},
        {"verbose",   no_argument, &verbose, 1},
        {"quiet",     no_argument, &verbose, 0},
        {"help",      no_argument, &help,    1},
        {0, 0, 0, 0}
      };

    /* getopt_long stores the option index here. */
    int option_index = 0;

    int c = getopt_long (argc, argv, options,
                     long_options, &option_index);

    /* Detect the end of the options. */
    if (c == -1)
      break;

    if (c == 0 && long_options[option_index].flag == 0)
      c = long_options[option_index].val;

    switch (c) {
      case 0:
        /* If this option set a flag, do nothing else now. */
        break;

      case 'B':
        errno = 0;
        longarg = strtol(optarg, 0, 10);
        if (errno != 0 || longarg < 1) {
          printf("option -B takes an integer maximum block size >= 1\n");
          usage();
          return 1;
        }
        B = longarg;
        break;

      case 'e':
        errno = 0;
        doublearg = strtod(optarg, 0);
        if (errno != 0 || doublearg < 0.0) {
          printf("option -e takes a desired relative error >= 0.0\n");
          usage();
          return 1;
        }
        epsilon = doublearg;
        break;

      case 'd':
        errno = 0;
        doublearg = strtod(optarg, 0);
        if (errno != 0 || doublearg < 0.0 || doublearg > 1.0) {
          printf("option -d takes a desired probability >= 0.0 and <= 1.0\n");
          usage();
          return 1;
        }
        delta = doublearg;
        break;

      case 't':
        errno = 0;
        longarg = strtol(optarg, 0, 10);
        if (errno != 0 || longarg < 1) {
          printf("option -t takes an integer number of trials >= 1\n");
          usage();
          return 1;
        }
        trials = longarg;
        break;

      case 'c':
        clock = 1;
        break;

      case 'C':
        clock = 0;
        break;

      case 'r':
        results = 1;
        break;

      case 'R':
        results = 0;
        break;

      case 'v':
        verbose = 1;
        break;

      case 'q':
        verbose = 0;
        break;

      case 'h':
        help = 1;
        break;

      case '?':
        usage();
        return 1;

      default:
        abort();
    }
  }

  if (help) {
    printf("Run a fill estimation algorithm!\n");
    usage();
    return 0;
  }

  if (argc - optind > 1) {
    printf("<input> cannot be more than one file\n");
    usage();
    return 1;
  }

  if (argc - optind < 1) {
    printf("<input> not specified\n");
    usage();
    return 1;
  }

  struct stat statthing;
  if (stat(argv[optind], &statthing) < 0 || !S_ISREG(statthing.st_mode)){
    printf("<input> must be filename of MatrixMarket matrix\n");
    usage();
    return 1;
  }

  FILE *f = fopen(argv[optind], "r");
  gsl_spmatrix *triples = gsl_spmatrix_fscanf(f);
  fclose(f);
  if (triples == 0) {
    printf("<input> must be filename of MatrixMarket matrix\n");
    usage();
    return 1;
  }

  gsl_spmatrix *csr = gsl_spmatrix_crs(triples);
  gsl_spmatrix_free(triples);

  //gsl_spmatrix_fprintf(stdout, csr, "%g");

  int ret = test(csr->size1, csr->size2, csr->nz, csr->p, csr->i, B, epsilon, delta, trials, clock, results, verbose);

  gsl_spmatrix_free(csr);

  return ret;
}


text:
/*
 * Copyright (C) 2015  University of Oregon
 *
 * You may distribute under the terms of either the GNU General Public
 * License or the Apache License, as specified in the LICENSE file.
 *
 * For more information, see the LICENSE file.
 */
/*---------------------------------------------------------------------------*/
/* 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 of the License, or         */
/* (at your option) any later version.                                       */
/*                                                                           */
/* This 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         */
/* If not, see <http://www.gnu.org/licenses/>.                               */
/*---------------------------------------------------------------------------*/
/**/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <sys/param.h>
#include <sys/stat.h>
#include <stddef.h>
#include <sys/time.h>

#define TRUE       1
#define FALSE      0

#define MAXSTR 1024
#define MAXELEM 1024

/* Some data types */
#define INT16      1  /* 16 bit integer */
#define INT32      2  /* 32 bit integer */
#define FLT32      3  /* 32 bit float */
#define DBL64      4  /* 64 bit double */

/* Include GSL libraries */
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>

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

/* Include local functions */
#include "gsl_functions.c"  /* local GSL functions */
#include "utils.c"
#include "getpars.c"
#include "fdf.c"
#include "image_functions.c"

/* Include Analyze libraries */
#include "AVW.c"  /* local AVW functions */

#ifdef NOTLIKELY
#include <AVW.h>
#include <AVW_ImageFile.h>
#include <AVW.c>
#endif
//#include "avw.c"

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
/*
 *  SPDX-FileCopyrightText: Copyright 2021, Siavash Ameli <sameli@berkeley.edu>
 *  SPDX-License-Identifier: BSD-3-Clause
 *  SPDX-FileType: SOURCE
 *
 *  This program is free software: you can redistribute it and/or modify it
 *  under the terms of the license found in the LICENSE.txt file in the root
 *  directory of this source tree.
 */

#ifndef _C_BASIC_ALGEBRA_CBLAS_INTERFACE_H_
#define _C_BASIC_ALGEBRA_CBLAS_INTERFACE_H_

// =======
// Headers
// =======

#include "../_definitions/definitions.h"  // USE_CBLAS

#if (USE_CBLAS == 1)

#include <cblas.h>  // CBLAS_LAYOUT, CBLAS_TRANSPOSE, cblas_sgemv, cblas_dgemv,
                    // cblas_scopy, cblas_dcopy, cblas_saxpy, cblas_daxpy,
                    // cblas_snrm2, cblas_dnrm2, cblas_sscal, cblas_dscal

typedef CBLAS_ORDER CBLAS_LAYOUT;  // backward compatibility with CBLAS_LAYOUT

// ===============
// cblas interface
// ===============

/// \namespace cblas_interface
///
/// \brief     A collection of templates to wrapper cblas functions.

namespace cblas_interface
{
    // cblas xgemv
    template <typename DataType>
    void xgemv(
            CBLAS_LAYOUT layout,
            CBLAS_TRANSPOSE TransA,
            const int M,
            const int N,
            const DataType alpha,
            const DataType* A,
            const int lda,
            const DataType* X,
            const int incX,
            const DataType beta,
            DataType* Y,
            const int incY);

    // cblas xcopy
    template <typename DataType>
    void xcopy(
            const int N,
            const DataType* X,
            const int incX,
            DataType* Y,
            const int incY);

    // cblas xaxpy
    template <typename DataType>
    void xaxpy(
            const int N,
            const DataType alpha,
            const DataType* X,
            const int incX,
            DataType* Y,
            const int incY);

    // cblas xdot
    template <typename DataType>
    DataType xdot(
            const int N,
            const DataType* X,
            const int incX,
            const DataType* Y,
            const int incY);

    // cblas xnrm2
    template <typename DataType>
    DataType xnrm2(
            const int N,
            const DataType* X,
            const int incX);

    // cblas xscal
    template <typename DataType>
    void xscal(
            const int N,
            const DataType alpha,
            DataType* X,
            const int incX);

}  // namespace cblas_interface

#endif  // USE_CBLAS
#endif  // _C_BASIC_ALGEBRA_CBLAS_INTERFACE_H_

max_stars_repo_name: ameli/TraceInv

max_issues_repo_name: ameli/TraceInv

max_forks_repo_name: ameli/TraceInv


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

/* Author:  G. Jungman */

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

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

#define ROOT_EIGHT (2.0*M_SQRT2)

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

/* based on SLATEC besi1(), besi1e() */

/* chebyshev expansions

 series for bi1        on the interval  0.          to  9.00000d+00
                                        with weighted error   2.40e-17
                                         log weighted error  16.62
                               significant figures required  16.23
                                    decimal places required  17.14

 series for ai1        on the interval  1.25000d-01 to  3.33333d-01
                                        with weighted error   6.98e-17
                                         log weighted error  16.16
                               significant figures required  14.53
                                    decimal places required  16.82

 series for ai12       on the interval  0.          to  1.25000d-01
                                       with weighted error   3.55e-17
                                        log weighted error  16.45
                              significant figures required  14.69
                                   decimal places required  17.12
*/

static double bi1_data[11] = {
  -0.001971713261099859,
   0.407348876675464810,
   0.034838994299959456,
   0.001545394556300123,
   0.000041888521098377,
   0.000000764902676483,
   0.000000010042493924,
   0.000000000099322077,
   0.000000000000766380,
   0.000000000000004741,
   0.000000000000000024
};
static cheb_series bi1_cs = {
  bi1_data,
  10,
  -1, 1,
  10
};

static double ai1_data[21] = {
  -0.02846744181881479,
  -0.01922953231443221,
  -0.00061151858579437,
  -0.00002069971253350,
   0.00000858561914581,
   0.00000104949824671,
  -0.00000029183389184,
  -0.00000001559378146,
   0.00000001318012367,
  -0.00000000144842341,
  -0.00000000029085122,
   0.00000000012663889,
  -0.00000000001664947,
  -0.00000000000166665,
   0.00000000000124260,
  -0.00000000000027315,
   0.00000000000002023,
   0.00000000000000730,
  -0.00000000000000333,
   0.00000000000000071,
  -0.00000000000000006
};
static cheb_series ai1_cs = {
  ai1_data,
  20,
  -1, 1,
  11
};

static double ai12_data[22] = {
   0.02857623501828014,
  -0.00976109749136147,
  -0.00011058893876263,
  -0.00000388256480887,
  -0.00000025122362377,
  -0.00000002631468847,
  -0.00000000383538039,
  -0.00000000055897433,
  -0.00000000001897495,
   0.00000000003252602,
   0.00000000001412580,
   0.00000000000203564,
  -0.00000000000071985,
  -0.00000000000040836,
  -0.00000000000002101,
   0.00000000000004273,
   0.00000000000001041,
  -0.00000000000000382,
  -0.00000000000000186,
   0.00000000000000033,
   0.00000000000000028,
  -0.00000000000000003
};
static cheb_series ai12_cs = {
  ai12_data,
  21,
  -1, 1,
  9
};

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

int gsl_sf_bessel_I1_scaled_e(const double x, gsl_sf_result * result)
{
  const double xmin    = 2.0 * GSL_DBL_MIN;
  const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
  const double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y < xmin) {
    UNDERFLOW_ERROR(result);
  }
  else if(y < x_small) {
    result->val = 0.5*x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y <= 3.0) {
    const double ey = exp(-y);
    gsl_sf_result c;
    cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
    result->val  = x * ey * (0.875 + c.val);
    result->err  = ey * c.err + y * GSL_DBL_EPSILON * fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(y <= 8.0) {
    const double sy = sqrt(y);
    gsl_sf_result c;
    double b;
    double s;
    cheb_eval_e(&ai1_cs, (48.0/y-11.0)/5.0, &c);
    b = (0.375 + c.val) / sy;
    s = (x > 0.0 ? 1.0 : -1.0);
    result->val  = s * b;
    result->err  = c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sy = sqrt(y);
    gsl_sf_result c;
    double b;
    double s;
    cheb_eval_e(&ai12_cs, 16.0/y-1.0, &c);
    b = (0.375 + c.val) / sy;
    s = (x > 0.0 ? 1.0 : -1.0);
    result->val  = s * b;
    result->err  = c.err / sy;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_I1_e(const double x, gsl_sf_result * result)
{
  const double xmin    = 2.0 * GSL_DBL_MIN;
  const double x_small = ROOT_EIGHT * GSL_SQRT_DBL_EPSILON;
  const double y = fabs(x);

  /* CHECK_POINTER(result) */

  if(y == 0.0) {
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y < xmin) {
    UNDERFLOW_ERROR(result);
  }
  else if(y < x_small) {
    result->val = 0.5*x;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
  else if(y <= 3.0) {
    gsl_sf_result c;
    cheb_eval_e(&bi1_cs, y*y/4.5-1.0, &c);
    result->val  = x * (0.875 + c.val);
    result->err  = y * c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(y < GSL_LOG_DBL_MAX) {
    const double ey = exp(y);
    gsl_sf_result I1_scaled;
    gsl_sf_bessel_I1_scaled_e(x, &I1_scaled);
    result->val  = ey * I1_scaled.val;
    result->err  = ey * I1_scaled.err + y * GSL_DBL_EPSILON * fabs(result->val);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    OVERFLOW_ERROR(result);
  }
}

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

#include "eval.h"

double gsl_sf_bessel_I1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_I1_scaled_e(x, &result));
}

double gsl_sf_bessel_I1(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_I1_e(x, &result));
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// -*- C++ -*-
//
// michael a.g. aïvázis <michael.aivazis@para-sim.com>
//
// (c) 2013-2020 parasim inc
// (c) 2010-2020 california institute of technology
// all rights reserved
//

// code guard
#if !defined(altar_bayesian_CoolingStep_h)
#define altar_bayesian_CoolingStep_h

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

// place everything in the local namespace
namespace altar {
    namespace bayesian {

        // forward declarations
        class CoolingStep;

    } // of namespace bayesian
} // of namespace altar

// declaration
class altar::bayesian::CoolingStep
{
    // types
public:
    typedef gsl_vector vector_t;
    typedef gsl_matrix matrix_t;

    // interface
public:
    inline auto samples() const;
    inline auto parameters() const;

    inline auto beta() const;
    inline void beta(double value);

    inline auto theta() const;
    inline auto prior() const;
    inline auto data() const;
    inline auto posterior() const;
    inline auto sigma() const;

    // meta-methods
public:
    virtual ~CoolingStep();
    inline CoolingStep(size_t samples, size_t parameters, double beta=0.0);

    // data
private:
    const size_t _samples; // the number of samples
    const size_t _parameters; // the number of model parameters

    double _beta; // the current temperature
    matrix_t * const _theta; // the sample set, a (samples x parameters) matrix

    vector_t * const _prior; // the log likelihoods of the samples in the model prior
    vector_t * const _data; // the log likelihoods of the samples given the data
    vector_t * const _posterior; // the posterior log likelihood

    matrix_t * const _sigma; // the parameter covariance

    // disallow
private:
    CoolingStep(const CoolingStep &) = delete;
    const CoolingStep & operator=(const CoolingStep &) = delete;
};

// get the inline definitions
#define altar_bayesian_CoolingStep_icc
#include "CoolingStep.icc"
#undef altar_bayesian_CoolingStep_icc

# endif
// end of file

max_stars_repo_name: AlTarFramework/altar

max_issues_repo_name: AlTarFramework/altar

max_forks_repo_name: AlTarFramework/altar


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

#pragma once

#include <memory>
#include <numeric>
#include <algorithm>
#include <vector>
#include <string>
#include <limits>
#include <memory>
#include <optional>
#include <list>
#include <map>
#include <deque>
#include <chrono>
#include <variant>
#include <cassert>

#include <wrl/client.h>
#include <wrl/implements.h>

#include <wil/wrl.h>
#include <wil/result.h>

#include <gsl/gsl>

#include <d3d12.h>
#include <d3d12sdklayers.h>
#include "External/D3DX12/d3dx12.h"

#include <DirectML.h>
#include "core/common/common.h"
#include "ErrorHandling.h"

// DirectML helper libraries
#include "External/DirectMLHelpers/ApiTraits.h"
#include "External/DirectMLHelpers/ApiHelpers.h"
#include "External/DirectMLHelpers/DirectMLSchema.h"
#include "External/DirectMLHelpers/AbstractOperatorDesc.h"
#include "External/DirectMLHelpers/GeneratedSchemaTypes.h"
#include "External/DirectMLHelpers/SchemaHelpers.h"
#include "External/DirectMLHelpers/GeneratedSchemaHelpers.h"

using Microsoft::WRL::ComPtr;

// Windows pollutes the macro space, causing a build break in schema.h.
#undef OPTIONAL

#include "core/providers/dml/DmlExecutionProvider/inc/DmlExecutionProvider.h"
#include "core/providers/dml/OperatorAuthorHelper/MLOperatorAuthorHelper.h"
#include "core/providers/dml/OperatorAuthorHelper/Common.h"

#include "DmlCommon.h"
#include "TensorDesc.h"
#include "DescriptorPool.h"
#include "IExecutionProvider.h"


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

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

/* The t-distribution has the form

   p(x) dx = (Gamma((nu + 1)/2)/(sqrt(pi nu) Gamma(nu/2))
   * (1 + (x^2)/nu)^-((nu + 1)/2) dx

   The method used here is the one described in Knuth */

double
gsl_ran_tdist (const gsl_rng * r, const double nu)
{
  if (nu <= 2)
    {
      double Y1 = gsl_ran_ugaussian (r);
      double Y2 = gsl_ran_chisq (r, nu);

      double t = Y1 / sqrt (Y2 / nu);

      return t;
    }
  else
    {
      double Y1, Y2, Z, t;
      do
	{
	  Y1 = gsl_ran_ugaussian (r);
	  Y2 = gsl_ran_exponential (r, 1 / (nu/2 - 1));

	  Z = Y1 * Y1 / (nu - 2);
	}
      while (1 - Z < 0 || exp (-Y2 - Z) > (1 - Z));

      /* Note that there is a typo in Knuth's formula, the line below
	 is taken from the original paper of Marsaglia, Mathematics of
	 Computation, 34 (1980), p 234-256 */

      t = Y1 / sqrt ((1 - 2 / nu) * (1 - Z));
      return t;
    }
}

double
gsl_ran_tdist_pdf (const double x, const double nu)
{
  double p;

  double lg1 = gsl_sf_lngamma (nu / 2);
  double lg2 = gsl_sf_lngamma ((nu + 1) / 2);

  p = ((exp (lg2 - lg1) / sqrt (M_PI * nu))
       * pow ((1 + x * x / nu), -(nu + 1) / 2));
  return p;
}


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

#pragma once

#ifndef type_f33c5e03_70f2_4879_9504_332323832a40_h
#define type_f33c5e03_70f2_4879_9504_332323832a40_h

#include <ariel/config.h>
#include <gslib/type.h>
#include <gslib/string.h>
#include <gslib/math.h>
#include <gslib/std.h>

__ariel_begin__

struct ariel_export color
{
    union
    {
        struct { byte red, green, blue, alpha; };
        uint _data;
    };

public:
    color() { data() = 0; }
    color(int r, int g, int b) { set_color(r, g, b); }
    color(int r, int g, int b, int a) { set_color(r, g, b, a); }
    uint& data() { return _data; }
    uint data() const { return _data; }
    void set_color(int r, int g, int b) { red = r, green = g, blue = b, alpha = 255; }
    void set_color(int r, int g, int b, int a) { red = r, green = g, blue = b, alpha = a; }
    bool operator !=(const color& cr) const { return blue != cr.blue || green != cr.green || red != cr.red || alpha != cr.alpha; }
    bool operator ==(const color& cr) const { return blue == cr.blue && green == cr.green && red == cr.red && alpha == cr.alpha; }
    color& lerp(const color& c1, const color& c2, float s)
    {
        assert(s >= 0.f && s <= 1.f);
        float t = 1.f - s;
        red = (int)(s * c1.red + t * c2.red);
        green = (int)(s * c1.green + t * c2.green);
        blue = (int)(s * c1.blue + t * c2.blue);
        alpha = (int)(s * c1.alpha + t * c2.alpha);
        return *this;
    }
};

struct ariel_export font
{
    enum
    {
        declare_mask(ftm_italic,    0),
        declare_mask(ftm_underline, 1),
        declare_mask(ftm_strikeout, 2),
    };

public:
    string      name;
    int         size;
    int         escape;
    int         orient;
    int         weight; /* 0-9 */
    uint        mask;

private:
    friend class fsys_win32;
    friend class fsys_dwrite;
    mutable uint sysfont;

public:
    font()
    {
        size = 0;
        escape = 0;
        orient = 0;
        weight = 3;
        mask = 0;
        sysfont = 0;
    }
    font(const font& that)
    {
        name    = that.name;
        size    = that.size;
        escape  = that.escape;
        orient  = that.orient;
        weight  = that.weight;
        mask    = that.mask;
        sysfont = that.sysfont;
    }
    font(const gchar* n, int sz)
    {
        name.assign(n);
        size = sz;
        escape = 0;
        orient = 0;
        weight = 3;
        mask = 0;
        sysfont = 0;
    }
    font& operator = (const font& that)
    {
        name    = that.name;
        size  = that.size;
        escape  = that.escape;
        orient  = that.orient;
        weight  = that.weight;
        mask    = that.mask;
        return *this;
    }
    bool operator == (const font& that) const
    {
        if(name != that.name)
            return false;
        if(size != that.size)
            return false;
        if(escape != that.escape)
            return false;
        if(orient != that.orient)
            return false;
        if(weight != that.weight)
            return false;
        if(mask != that.mask)
            return false;
        return true;
    }
    bool operator != (const font& that) const
    {
        if(name != that.name)
            return true;
        if(size != that.size)
            return true;
        if(escape != that.escape)
            return true;
        if(orient != that.orient)
            return true;
        if(weight != that.weight)
            return true;
        if(mask != that.mask)
            return true;
        return false;
    }
    size_t hash_value() const
    {
        hasher h;
        h.add_bytes((const byte*)name.c_str(), name.length() * sizeof(gchar));
        h.add_bytes((const byte*)&size, sizeof(size));
        h.add_bytes((const byte*)&escape, sizeof(escape));
        h.add_bytes((const byte*)&orient, sizeof(orient));
        h.add_bytes((const byte*)&weight, sizeof(weight));
        return h.add_bytes((const byte*)&mask, sizeof(mask));
    }
};

struct ariel_export viewport
{
    float           left;
    float           top;
    float           width;
    float           height;
    float           min_depth;
    float           max_depth;
};

struct ariel_export axis_aligned_bound_box
{
    float           left = FLT_MAX;
    float           right = -FLT_MAX;
    float           top = FLT_MAX;
    float           bottom = -FLT_MAX;
    float           front = FLT_MAX;
    float           back = -FLT_MAX;

public:
    void reset()
    {
        left = top = front = FLT_MAX;
        right = bottom = back = -FLT_MAX;
    }
    float width() const { return right - left; }
    float height() const { return bottom - top; }
    float depth() const { return back - front; }
};

struct ariel_export origin_bound_sphere
{
    float           radius = 0.f;
};

struct ariel_export bound_sphere
{
    vec3            origin;
    float           radius = 0.f;
};

enum res_type
{
    res_mesh,
};

class __gs_novtable ariel_export res_node abstract
{
public:
    virtual ~res_node() {}
    virtual res_type get_type() const = 0;
    virtual const string& get_name() const { return _name; }
    virtual bool has_name() const { return !_name.empty(); }

protected:
    string          _name;

public:
    void set_name(const string& name) { _name = name; }
};

__ariel_end__

namespace std {

template<>
class hash<gs::ariel::font>
#if defined(_MSC_VER) && (_MSC_VER < 1914)
    : public unary_function<gs::ariel::font, size_t>
#endif
{
public:
    size_t operator()(const gs::ariel::font& ft) const { return ft.hash_value(); }
};

};

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#ifndef ___BEST_SCORE_SO_FAR_FOR_NLOPT_H____
#define ___BEST_SCORE_SO_FAR_FOR_NLOPT_H____

#include "OptimizedNonlinearController.h"
#include <mutex>
#include <nlopt.h>
#include <thread>

class BestNLOptScoreSoFar
{
public:
	std::mutex my_mutex;
	std::thread* thethread;

	bool force_stop;
	nlopt_opt opt;

	MyNLOPTdata my_nloptdata;
	int num_iterations_so_far;

	double* my_best_params;
	double my_best_score;

	double original_starting_score;
	int which_algorithm;

	BestNLOptScoreSoFar() : my_best_params(nullptr), thethread(nullptr), opt(nullptr), force_stop(false) {}
	BestNLOptScoreSoFar(const BestNLOptScoreSoFar & other) : my_nloptdata(other.my_nloptdata), thethread(other.thethread), my_best_params(other.my_best_params), force_stop(other.force_stop), opt(other.opt) {}

	void Reset() {
		num_iterations_so_far = 0;
		if(my_best_params != nullptr) {free(my_best_params); my_best_params = nullptr;}
		my_best_score = 1e12;
		original_starting_score = 1e12;
		which_algorithm = 0;
		force_stop = false;
		if(opt != nullptr) {nlopt_destroy(opt);}
		opt = nullptr;
	}
};

#endif

max_stars_repo_name: jasonbunk/InvertedPendulumVisualControl

max_issues_repo_name: jasonbunk/InvertedPendulumVisualControl

max_forks_repo_name: jasonbunk/InvertedPendulumVisualControl


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cblas.h>
#include "linear_algebra.h"

Matrix *new_matrix(int row_num, int col_num)
{
    Matrix *matrix = malloc(sizeof(Matrix));
    if (matrix == NULL) return NULL;
    matrix->data = malloc((unsigned long) (row_num * col_num) * sizeof(double));

    if (matrix->data == NULL)
    {
        perror("Error allocating memory for a matrix");
        exit(EXIT_FAILURE);
    }

    matrix->row_num = row_num;
    matrix->col_num = col_num;

    return matrix;
}

Matrix *new_matrix_from_array(double *array, int row_num, int col_num)
{
   Matrix *matrix = new_matrix(row_num, col_num);
   memcpy(matrix->data, array, (unsigned long)(row_num * col_num) * sizeof(double));
   return matrix;
}

void free_matrix(Matrix *matrix)
{
    free(matrix->data);
    free(matrix);
}

Matrix *add_matrices(Matrix *matrix1, Matrix *matrix2)
{
    if (matrix1->row_num != matrix2->row_num ||
        matrix1->col_num != matrix2->col_num)
    {
        perror("Dimensions of the two matrices do not match");
        exit(EXIT_FAILURE);
    }

    Matrix *result_matrix = new_matrix(matrix1->row_num, matrix1->col_num);

    int i;
    int elements = matrix1->row_num * matrix1->col_num;

    for (i = 0; i < elements; i++)
    {
        result_matrix->data[i] = matrix1->data[i] + matrix2->data[i];
    }

    return result_matrix;
}

Matrix *multiply_matrices(Matrix *matrix1, Matrix *matrix2)
{
    if (matrix1->col_num != matrix2->row_num)
    {
        perror("Incompatible matrix dimensions");
        exit(EXIT_FAILURE);
    }

    Matrix *product = new_matrix(matrix1->row_num, matrix2->col_num);

    int m = matrix1->row_num;
    int n = matrix2->col_num;
    int k = matrix1->col_num;

    // Use high-performace matrix multiplication from BLAS
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1,
        matrix1->data, k, matrix2->data, n, 0, product->data, n);

    return product;
}

void multiply_vector_by_vector_transpose(Matrix *vector1, Matrix *vector2, double scalar, Matrix *matrix)
{
    if (vector1->col_num != 1 || vector2->col_num != 1)
    {
        perror("Not a vector.");
        exit(EXIT_FAILURE);
    }

    int m = vector1->row_num;
    int n = vector2->row_num;

    cblas_dger(CblasRowMajor, m, n, scalar, vector1->data, 1, vector2->data, 1, matrix->data, n);
}

Matrix *multiply_matrix_with_vector(Matrix *matrix, Matrix *vector, double scalar)
{
    if (matrix->col_num != vector->row_num)
    {
        perror("Incompatible matrix dimensions");
        exit(EXIT_FAILURE);
    }

    if (vector->col_num != 1)
    {
        perror("Not a vector.");
        exit(EXIT_FAILURE);
    }

    Matrix *product = new_matrix(matrix->row_num, 1);

    int m = matrix->row_num;
    int n = matrix->col_num;

    cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, scalar, matrix->data, n,
                vector->data, 1,  0, product->data, 1);

    return product;
}

Matrix *multiply_upper_symmetric_matrix_with_vector(Matrix *matrix, Matrix *vector)
{
    if (matrix->col_num != matrix->row_num)
    {
        perror("Matrix is not square");
        exit(EXIT_FAILURE);
    }

    if (matrix->col_num != vector->row_num)
    {
        perror("Incompatible matrix dimensions");
        exit(EXIT_FAILURE);
    }

    if (vector->col_num != 1)
    {
        perror("Not a vector.");
        exit(EXIT_FAILURE);
    }

    Matrix *product = new_matrix(matrix->row_num, 1);

    int n = matrix->col_num;

    cblas_dsymv(CblasRowMajor, CblasUpper, n, 1, matrix->data,
                n, vector->data, 1, 0, product->data, 1);

    return product;
}

void multiply_matrix_with_a_number(Matrix *matrix, double number)
{
    int i;
    int n = matrix->col_num;
    int row_num = matrix->row_num;

    for (i = 0; i < row_num; i++)
    {
        cblas_dscal(n, number, matrix->data + i * n, 1);
    }
}

Matrix *transpose_matrix(Matrix *matrix)
{
    Matrix *product = new_matrix(matrix->col_num, matrix->row_num);

    int i, j;

    for (i = 0; i < matrix->col_num; i++)
    {
        for (j = 0; j < matrix->row_num; j++)
        {
            product->data[i * matrix->row_num + j] = matrix->data[j * matrix->col_num + i];
        }
    }

    return product;
}

double dot_product(Matrix *matrix1, Matrix *matrix2)
{
    if (matrix1->row_num != matrix2->row_num ||
        matrix1->col_num != 1 || matrix2->col_num != 1)
    {
        perror("Matrices are not n by 1.");
        exit(EXIT_FAILURE);
    }

    int i;
    double result = 0;

    for (i = 0; i < matrix1->row_num; i++)
    {
        result += matrix1->data[i] * matrix2->data[i];
    }

    return result;
}

double norm(Matrix *matrix)
{
    return sqrt(dot_product(matrix, matrix));
}

Matrix *gramian(Matrix *matrix)
{
    Matrix *product = new_matrix(matrix->col_num, matrix->col_num);

    int n = matrix->col_num;
    int k = matrix->row_num;

    cblas_dsyrk(CblasRowMajor, CblasUpper, CblasTrans,
		 n, k, 1, matrix->data, n, 0, product->data, n);

    return product;
}


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

#ifndef __GSL_VECTOR_FLOAT_H__
#define __GSL_VECTOR_FLOAT_H__

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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  float *data;
  gsl_block_float *block;
  int owner;
}
gsl_vector_float;

typedef struct
{
  gsl_vector_float vector;
} _gsl_vector_float_view;

typedef _gsl_vector_float_view gsl_vector_float_view;

typedef struct
{
  gsl_vector_float vector;
} _gsl_vector_float_const_view;

typedef const _gsl_vector_float_const_view gsl_vector_float_const_view;

/* Allocation */

GSL_FUN gsl_vector_float *gsl_vector_float_alloc (const size_t n);
GSL_FUN gsl_vector_float *gsl_vector_float_calloc (const size_t n);

GSL_FUN gsl_vector_float *gsl_vector_float_alloc_from_block (gsl_block_float * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_float *gsl_vector_float_alloc_from_vector (gsl_vector_float * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_float_free (gsl_vector_float * v);

/* Views */

GSL_FUN _gsl_vector_float_view
gsl_vector_float_view_array (float *v, size_t n);

GSL_FUN _gsl_vector_float_view
gsl_vector_float_view_array_with_stride (float *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_float_const_view_array (const float *v, size_t n);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_float_const_view_array_with_stride (const float *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_float_view
gsl_vector_float_subvector (gsl_vector_float *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_float_view
gsl_vector_float_subvector_with_stride (gsl_vector_float *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_float_const_subvector (const gsl_vector_float *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_float_const_subvector_with_stride (const gsl_vector_float *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_float_set_zero (gsl_vector_float * v);
GSL_FUN void gsl_vector_float_set_all (gsl_vector_float * v, float x);
GSL_FUN int gsl_vector_float_set_basis (gsl_vector_float * v, size_t i);

GSL_FUN int gsl_vector_float_fread (FILE * stream, gsl_vector_float * v);
GSL_FUN int gsl_vector_float_fwrite (FILE * stream, const gsl_vector_float * v);
GSL_FUN int gsl_vector_float_fscanf (FILE * stream, gsl_vector_float * v);
GSL_FUN int gsl_vector_float_fprintf (FILE * stream, const gsl_vector_float * v,
                              const char *format);

GSL_FUN int gsl_vector_float_memcpy (gsl_vector_float * dest, const gsl_vector_float * src);

GSL_FUN int gsl_vector_float_reverse (gsl_vector_float * v);

GSL_FUN int gsl_vector_float_swap (gsl_vector_float * v, gsl_vector_float * w);
GSL_FUN int gsl_vector_float_swap_elements (gsl_vector_float * v, const size_t i, const size_t j);

GSL_FUN float gsl_vector_float_max (const gsl_vector_float * v);
GSL_FUN float gsl_vector_float_min (const gsl_vector_float * v);
GSL_FUN void gsl_vector_float_minmax (const gsl_vector_float * v, float * min_out, float * max_out);

GSL_FUN size_t gsl_vector_float_max_index (const gsl_vector_float * v);
GSL_FUN size_t gsl_vector_float_min_index (const gsl_vector_float * v);
GSL_FUN void gsl_vector_float_minmax_index (const gsl_vector_float * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_float_add (gsl_vector_float * a, const gsl_vector_float * b);
GSL_FUN int gsl_vector_float_sub (gsl_vector_float * a, const gsl_vector_float * b);
GSL_FUN int gsl_vector_float_mul (gsl_vector_float * a, const gsl_vector_float * b);
GSL_FUN int gsl_vector_float_div (gsl_vector_float * a, const gsl_vector_float * b);
GSL_FUN int gsl_vector_float_scale (gsl_vector_float * a, const float x);
GSL_FUN int gsl_vector_float_add_constant (gsl_vector_float * a, const double x);
GSL_FUN int gsl_vector_float_axpby (const float alpha, const gsl_vector_float * x, const float beta, gsl_vector_float * y);
GSL_FUN float gsl_vector_float_sum (const gsl_vector_float * a);

GSL_FUN int gsl_vector_float_equal (const gsl_vector_float * u,
                            const gsl_vector_float * v);

GSL_FUN int gsl_vector_float_isnull (const gsl_vector_float * v);
GSL_FUN int gsl_vector_float_ispos (const gsl_vector_float * v);
GSL_FUN int gsl_vector_float_isneg (const gsl_vector_float * v);
GSL_FUN int gsl_vector_float_isnonneg (const gsl_vector_float * v);

GSL_FUN INLINE_DECL float gsl_vector_float_get (const gsl_vector_float * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_float_set (gsl_vector_float * v, const size_t i, float x);
GSL_FUN INLINE_DECL float * gsl_vector_float_ptr (gsl_vector_float * v, const size_t i);
GSL_FUN INLINE_DECL const float * gsl_vector_float_const_ptr (const gsl_vector_float * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
float
gsl_vector_float_get (const gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_float_set (gsl_vector_float * v, const size_t i, float x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
float *
gsl_vector_float_ptr (gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (float *) (v->data + i * v->stride);
}

INLINE_FUN
const float *
gsl_vector_float_const_ptr (const gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const float *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_FLOAT_H__ */


text:
#ifndef structs_h
#define structs_h

#include <ceed.h>
#include <petsc.h>

// -----------------------------------------------------------------------------
// PETSc Operator Structs
// -----------------------------------------------------------------------------

// Data for PETSc Matshell
typedef struct UserO_ *UserO;
struct UserO_ {
  MPI_Comm comm;
  DM dm;
  Vec X_loc, Y_loc, diag;
  CeedVector x_ceed, y_ceed;
  CeedOperator op;
  Ceed ceed;
};

// Data for PETSc Prolong/Restrict Matshells
typedef struct UserProlongRestr_ *UserProlongRestr;
struct UserProlongRestr_ {
  MPI_Comm comm;
  DM dmc, dmf;
  Vec loc_vec_c, loc_vec_f, mult_vec;
  CeedVector ceed_vec_c, ceed_vec_f;
  CeedOperator op_prolong, op_restrict;
  Ceed ceed;
};

// -----------------------------------------------------------------------------
// libCEED Data Structs
// -----------------------------------------------------------------------------

// libCEED data struct for level
typedef struct CeedData_ *CeedData;
struct CeedData_ {
  Ceed ceed;
  CeedBasis basis_x, basis_u, basis_c_to_f;
  CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_u_i, elem_restr_qd_i;
  CeedQFunction qf_apply;
  CeedOperator op_apply, op_restrict, op_prolong;
  CeedVector q_data, x_ceed, y_ceed;
};

// BP specific data
typedef struct {
  CeedInt num_comp_x, num_comp_u, topo_dim, q_data_size, q_extra;
  CeedQFunctionUser setup_geo, setup_rhs, apply, error;
  const char *setup_geo_loc, *setup_rhs_loc, *apply_loc, *error_loc;
  CeedEvalMode in_mode, out_mode;
  CeedQuadMode q_mode;
  PetscBool enforce_bc;
  PetscErrorCode (*bc_func)(PetscInt, PetscReal, const PetscReal *,
                            PetscInt, PetscScalar *, void *);
} BPData;

#endif // structs_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED