text:
#pragma once
#include <gsl/span>

/**
 * \brief Class for managing Bluetooth state.
 */
class Bluetooth
{
public:
	/**
	 * \brief Disconnects a device with matching MAC address from
	 * the first Bluetooth radio it is connected to.
	 * \param macAddress The MAC address to search for.
	 * \return \c true on success.
	 */
	static bool disconnectDevice(const gsl::span<uint8_t>& macAddress);
};


text:
/*
MIT License

Copyright (c) 2020 Huy Vo

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#ifndef PACMENSL_UTIL_H
#define PACMENSL_UTIL_H

#define ARMA_DONT_PRINT_ERRORS
#define ARMA_DONT_USE_WRAPPER
#include <armadillo>
#include <petscmat.h>
#include <petscvec.h>
#include <petscao.h>
#include <petscis.h>
#include <petscistypes.h>
#include <petscoptions.h>
#include <petsc.h>
#include <petscsys.h>
#include <petscconf.h>
#include <cassert>
#include <memory>
#include <mpi.h>
#include <zoltan.h>
#include <parmetis.h>
#include <string>
#include <sstream>
#include "ErrorHandling.h"

#define NOT_COPYABLE_NOT_MOVABLE(object)\
            object( const object & ) = delete;\
            object &operator=( const object & ) = delete;\

namespace pacmensl {
#define SQR1 sqrt(0.1e0)

/*! Round to 2 significant digits
 */
double round2digit(double x);

/*! Initialize and finalize Parallel context
 *
 */
int PACMENSLInit(int *argc, char ***argv, const char *help);

int PACMENSLFinalize();

void sequential_action(MPI_Comm comm, std::function<void(void *)> action, void *data);

class Environment {
 public:
  Environment();

  Environment(int *argc, char ***argv, const char *help);

  ~Environment();

 private:
  bool initialized = false;
  bool init_petsc = false; // If PETSc or MPI were already set, we do not meddle with them
  bool init_mpi = false;
};
}

#endif

ext: h

max_stars_repo_name: voduchuy/pacmensl

max_issues_repo_name: voduchuy/pacmensl

max_forks_repo_name: voduchuy/pacmensl


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

#include "energy.h"
#include <stdlib.h>
#include <assert.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>

/// \brief Allocate memory for energy levels structure
void AllocateEnergyLevels(const int m_max, const int num, en_levels_t *levels)
{
	levels->m_max = m_max;
	levels->num   = num;

	levels->en = malloc((m_max+1)*num * sizeof(double));
}

/// \brief Free memory of energy levels structure
void FreeEnergyLevels(en_levels_t *levels)
{
	free(levels->en);

	levels->num   = 0;
	levels->m_max = 0;
}

//_______________________________________________________________________________________________________________________
//

/// \brief Comparison function used for sorting
static int en_occupancy_cmp(const void *e1, const void *e2)
{
	const en_occupancy_t *a = (en_occupancy_t *)e1;
	const en_occupancy_t *b = (en_occupancy_t *)e2;

	// sort according to energy
	if (a->en < b->en)
	{
		return -1;
	}
	else if (a->en > b->en)
	{
		return 1;
	}

	// should actually never reach this point

	// compare 'm' quantum numbers
	if (a->m < b->m)
	{
		return -1;
	}
	else if (a->m > b->m)
	{
		return 1;
	}

	// compare energy quantum numbers
	if (a->n < b->n)
	{
		return -1;
	}
	else if (a->n > b->n)
	{
		return 1;
	}

	return 0;
}

/// \brief Get the maximum occupancy of an energy level, depending on 'm' quantum number and whether state is spin-polarized
static inline int GetMaxOccupancy(const int m, const bool spinpol)
{
	int occ = 1;

	if (m != 0)
	{
		occ *= 2;
	}
	if (!spinpol)
	{
		occ *= 2;
	}

	return occ;
}

//_______________________________________________________________________________________________________________________
//

/// \brief Fermi-Dirac distribution function
static double FermiDirac(const double kBT, const double mu, const double en)
{
	return 1 / (exp((en - mu)/kBT) + 1);
}

/// \brief GSL parameters for NumParticleDifference()
typedef struct
{
	int nelec;			//!< physical number of electrons
	int nlevels;		//!< number of energy levels
	double *en;			//!< energies
	int *occ;			//!< maximum occupancy of each energy level
	double kBT;			//!< k_B T (Boltzmann constant times temperature)
}
gsl_FD_params_t;

/// \brief Calculate difference in number of particles given a Fermi-Dirac distribution with chemical potential 'mu'
static double NumParticleDifference(double mu, void *p)
{
	int i;

	const gsl_FD_params_t *params = (gsl_FD_params_t *)p;

	double N = 0;
	for (i = 0; i < params->nlevels; i++)
	{
		N += params->occ[i] * FermiDirac(params->kBT, mu, params->en[i]);
	}

	return N - params->nelec;
}

//_______________________________________________________________________________________________________________________
///
/// \brief Calculate chemical potential 'mu' such that FermiDirac distribution gives exactly 'nelec' electrons
///
static double CalculateChemicalPotential(const en_levels_t *levels, const int nelec, const bool spinpol, const double kBT)
{
	int m, n, i;

	const int maxiter = 128;

	gsl_FD_params_t params;
	params.nelec = nelec;
	params.nlevels = (levels->m_max+1)*levels->num;
	params.en  = malloc(params.nlevels * sizeof(double));
	params.occ = malloc(params.nlevels * sizeof(int));
	params.kBT = kBT;
	// fill energy levels and occupancies
	for (m = 0; m <= levels->m_max; m++)
	{
		for (n = 0; n < levels->num; n++)
		{
			params.en [n + m*levels->num] = levels->en[n + m*levels->num];
			params.occ[n + m*levels->num] = GetMaxOccupancy(m, spinpol);
		}
	}

	// first interval used in interval bisection method
	double mu_lo = 0.0;
	double mu_hi = levels->en[(levels->m_max+1)*levels->num-1];		// use "last" energy level as upper bound

	gsl_function F;
	F.function = &NumParticleDifference;
	F.params = &params;

	const gsl_root_fsolver_type *T = gsl_root_fsolver_brent;
	gsl_root_fsolver *s = gsl_root_fsolver_alloc(T);
	gsl_root_fsolver_set(s, &F, mu_lo, mu_hi);

	double mu = 0;
	int status = GSL_CONTINUE;
	for (i = 0; i < maxiter && status == GSL_CONTINUE; i++)
	{
		status = gsl_root_fsolver_iterate(s);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "CalculateOccupancy() warning: gsl_root_fsolver_iterate() returned with status code %i\n", status);
		}

		mu = gsl_root_fsolver_root(s);
		mu_lo  = gsl_root_fsolver_x_lower(s);
		mu_hi  = gsl_root_fsolver_x_upper(s);
		status = gsl_root_test_interval(mu_lo, mu_hi, 0, 1e-8);
	}

	if (status != GSL_SUCCESS) {
		fprintf(stderr, "CalculateOccupancy() warning: not converged after %i iterations, status: %i\n", maxiter, status);
	}

	// clean up
	free(params.occ);
	free(params.en);
	gsl_root_fsolver_free(s);

	return mu;
}

//_______________________________________________________________________________________________________________________
///
/// \brief Calculate occupancy of energy levels
///
void CalculateOccupancy(const en_levels_t *levels, const int nelec, const bool spinpol, const double kBT, en_occlist_t *occlist)
{
	int i, m, n;

	// accumulate all energies into one list
	occlist->occupancy = malloc((levels->m_max+1)*levels->num * sizeof(en_occupancy_t));

	if (kBT == 0)	// zero temperature
	{
		for (m = 0; m <= levels->m_max; m++)
		{
			for (n = 0; n < levels->num; n++)
			{
				en_occupancy_t *encur = &occlist->occupancy[n + m*levels->num];
				encur->en  = levels->en[n + m*levels->num];
				encur->m   = m;
				encur->n   = n;
				encur->occ = 0;		// initially set to zero
			}
		}

		// sort energies
		qsort(occlist->occupancy, (levels->m_max+1)*levels->num, sizeof(en_occupancy_t), en_occupancy_cmp);

		// fill up energy levels
		en_occupancy_t *encur = occlist->occupancy;
		occlist->length = 1;
		for (i = 0; i < nelec; i++)
		{
			if (encur->occ >= GetMaxOccupancy(encur->m, spinpol)) {
				encur++;	// use next level
				occlist->length++;
			}

			encur->occ++;
		}
	}
	else	// k_B T > 0
	{
		assert(kBT > 0);

		const double mu = CalculateChemicalPotential(levels, nelec, spinpol, kBT);

		// all energy levels have nonzero occupancy since Fermi-Dirac function is strictly positive
		occlist->length = (levels->m_max+1)*levels->num;
		for (m = 0; m <= levels->m_max; m++)
		{
			for (n = 0; n < levels->num; n++)
			{
				en_occupancy_t *encur = &occlist->occupancy[n + m*levels->num];
				encur->en  = levels->en[n + m*levels->num];
				encur->m   = m;
				encur->n   = n;
				encur->occ = GetMaxOccupancy(m, spinpol) * FermiDirac(kBT, mu, encur->en);
			}
		}

		// sort energies
		qsort(occlist->occupancy, (levels->m_max+1)*levels->num, sizeof(en_occupancy_t), en_occupancy_cmp);
	}
}

/// \brief Free memory of energy level occupancy list
void FreeOccupancy(en_occlist_t *occlist)
{
	free(occlist->occupancy);
	occlist->length = 0;
}

ext: c

max_stars_repo_name: Arfaouim/quantum-dot

max_issues_repo_name: Arfaouim/quantum-dot

max_forks_repo_name: Arfaouim/quantum-dot


text:
/** @file */
#ifndef __CCL_HALOFIT_H_INCLUDED__
#define __CCL_HALOFIT_H_INCLUDED__

#include <gsl/gsl_spline.h>
#include <gsl/gsl_spline.h>

CCL_BEGIN_DECLS

typedef struct halofit_struct {
  gsl_spline *rsigma;
  gsl_spline *sigma2;
  gsl_spline *n_eff;
  gsl_spline *C;
  gsl_spline *weff;
  gsl_spline *omeff;
  gsl_spline *deeff;
} halofit_struct;

/*
 * Allocate a new struct for storing halofit data
 * @param cosmo Cosmological data
 * @return int, status of computations
 */
halofit_struct* ccl_halofit_struct_new(ccl_cosmology *cosmo, int *status);

/*
 * Free a halofit struct
 * @param hf, pointer to halofit struct to free
 */
void ccl_halofit_struct_free(halofit_struct *hf);

/**
 * Computes the halofit non-linear power spectrum
 * @param cosmo: cosmology object containing parameters
 * @param k: wavenumber in units of Mpc^{-1}
 * @param a: scale factor normalised to a=1 today
 * @param status: Status flag: 0 if there are no errors, non-zero otherwise
 * @param hf: halofit splines for evaluating the power spectrum
 * @return halofit_matter_power: halofit power spectrum, P(k), units of Mpc^{3}
 */
double ccl_halofit_power(ccl_cosmology *cosmo, double k, double a, halofit_struct *hf, int *status);

CCL_END_DECLS

#endif

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


text:
#pragma once

#include "HttpHeaders.h"
#include "Library.h"

#include <gsl/span>

#include <cstddef>
#include <map>
#include <string>

namespace CesiumAsync {

/**
 * @brief A completed response for a 3D Tiles asset.
 */
class CESIUMASYNC_API IAssetResponse {
public:
  /**
   * @brief Default destructor
   */
  virtual ~IAssetResponse() = default;

  /**
   * @brief Returns the HTTP response code.
   */
  virtual uint16_t statusCode() const = 0;

  /**
   * @brief Returns the HTTP content type
   */
  virtual std::string contentType() const = 0;

  /**
   * @brief Returns the HTTP headers of the response
   */
  virtual const HttpHeaders& headers() const = 0;

  /**
   * @brief Returns the data of this response
   */
  virtual gsl::span<const std::byte> data() const = 0;
};

} // namespace CesiumAsync

max_stars_repo_name: yieryi/cesium-native

max_issues_repo_name: yieryi/cesium-native

max_forks_repo_name: yieryi/cesium-native


text:
//! \file openmc_driver.h
//! Driver to initialize and run OpenMC in stages
#ifndef ENRICO_OPENMC_DRIVER_H
#define ENRICO_OPENMC_DRIVER_H

#include "enrico/cell_instance.h"
#include "enrico/geom.h"
#include "enrico/neutronics_driver.h"

#include "openmc/cell.h"
#include "openmc/tallies/filter_cell_instance.h"
#include "openmc/tallies/tally.h"
#include <gsl/gsl>
#include <mpi.h>

#include <vector>

namespace enrico {

//! Driver to initialize and run OpenMC in stages
class OpenmcDriver : public NeutronicsDriver {
public:
  //! One-time initalization of OpenMC and member variables
  //! \param comm An existing MPI communicator used to inialize OpenMC
  explicit OpenmcDriver(MPI_Comm comm);

  //! One-time finalization of OpenMC
  ~OpenmcDriver();

  //////////////////////////////////////////////////////////////////////////////
  // NeutronicsDriver interface

  //! Find cells corresponding to a vector of positions
  //! \param positions (x,y,z) coordinates to search for
  //! \return Handles to cells
  std::vector<CellHandle> find(const std::vector<Position>& position) override;

  //! Set the density of the material in a cell
  //! \param cell Handle to a cell
  //! \param rho Density in [g/cm^3]
  void set_density(CellHandle cell, double rho) const override;

  //! Set the temperature of a cell
  //! \param cell Handle to a cell
  //! \param T Temperature in [K]
  void set_temperature(CellHandle cell, double T) const override;

  //! Get the density of a cell
  //! \param cell Handle to a cell
  //! \return Cell density in [g/cm^3]
  double get_density(CellHandle cell) const override;

  //! Get the temperature of a cell
  //! \param cell Handle to a cell
  //! \return Temperature in [K]
  double get_temperature(CellHandle cell) const override;

  //! Get the volume of a cell
  //! \param cell Handle to a cell
  //! \return Volume in [cm^3]
  double get_volume(CellHandle cell) const override;

  //! Detemrine whether a cell contains fissionable nuclides
  //! \param cell Handle to a cell
  //! \return Whether the cell contains fissionable nuclides
  bool is_fissionable(CellHandle cell) const override;

  std::size_t n_cells() const override { return cells_.size(); }

  //! Create energy production tallies
  void create_tallies() override;

  //! Determine number of cells participating in coupling
  //! \return Number of cells
  xt::xtensor<double, 1> heat_source(double power) const final;

  std::string cell_label(CellHandle cell) const;

  //////////////////////////////////////////////////////////////////////////////
  // Driver interface

  //! Initialization required in each Picard iteration
  void init_step() final;

  //! Runs OpenMC for one Picard iteration
  void solve_step() final;

  //! Writes OpenMC output for given timestep and iteration
  //! \param timestep timestep index
  //! \param iteration iteration index
  void write_step(int timestep, int iteration) final;

  //! Finalization required in each Picard iteration
  void finalize_step() final;

private:
  // Data members
  openmc::Tally* tally_;               //!< Fission energy deposition tally
  openmc::CellInstanceFilter* filter_; //!< Cell instance filter
  std::vector<CellInstance> cells_;    //!< Array of cell instances
  int n_fissionable_cells_;            //!< Number of fissionable cells in model
};

} // namespace enrico

#endif // ENRICO_OPENMC_DRIVER_H

max_stars_repo_name: pshriwise/enrico

max_issues_repo_name: pshriwise/enrico

max_forks_repo_name: lebuller/enrico


text:
/**
 *
 * @file testing_zgels.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_zmain.h"

#undef REAL
#define COMPLEX

enum blas_order_type {
            blas_rowmajor = 101,
            blas_colmajor = 102 };

enum blas_uplo_type  {
            blas_upper = 121,
            blas_lower = 122 };

enum blas_cmach_type {
            blas_base      = 151,
            blas_t         = 152,
            blas_rnd       = 153,
            blas_ieee      = 154,
            blas_emin      = 155,
            blas_emax      = 156,
            blas_eps       = 157,
            blas_prec      = 158,
            blas_underflow = 159,
            blas_overflow  = 160,
            blas_sfmin     = 161};

enum blas_norm_type {
            blas_one_norm       = 171,
            blas_real_one_norm  = 172,
            blas_two_norm       = 173,
            blas_frobenius_norm = 174,
            blas_inf_norm       = 175,
            blas_real_inf_norm  = 176,
            blas_max_norm       = 177,
            blas_real_max_norm  = 178 };

static void
BLAS_error(char *rname, int err, int val, int x) {
  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
  abort();
}

static
void
BLAS_zsy_norm(enum blas_order_type order, enum blas_norm_type norm,
  enum blas_uplo_type uplo, int n, const PLASMA_Complex64_t *a, int lda, double *res) {
  int i, j; double anorm, v;
  char rname[] = "BLAS_zsy_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_inf_norm) {
    anorm = 0.0;
    if (blas_upper == uplo) {
      for (i = 0; i < n; ++i) {
        v = 0.0;
        for (j = 0; j < i; ++j) {
          v += cabs( a[j + i * lda] );
        }
        for (j = i; j < n; ++j) {
          v += cabs( a[i + j * lda] );
        }
        if (v > anorm)
          anorm = v;
      }
    } else {
      BLAS_error( rname, -3, norm, 0 );
      return;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
void
BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm,
  int m, int n, const PLASMA_Complex64_t *a, int lda, double *res) {
  int i, j; float anorm, v;
  char rname[] = "BLAS_zge_norm";

  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );

  if (norm == blas_frobenius_norm) {
    anorm = 0.0f;
    for (j = n; j; --j) {
      for (i = m; i; --i) {
        v = a[0];
        anorm += v * v;
        a++;
      }
      a += lda - m;
    }
    anorm = sqrt( anorm );
  } else if (norm == blas_inf_norm) {
    anorm = 0.0f;
    for (i = 0; i < m; ++i) {
      v = 0.0f;
      for (j = 0; j < n; ++j) {
        v += cabs( a[i + j * lda] );
      }
      if (v > anorm)
        anorm = v;
    }
  } else {
    BLAS_error( rname, -2, norm, 0 );
    return;
  }

  if (res) *res = anorm;
}

static
double
BLAS_dpow_di(double x, int n) {
  double rv = 1.0;

  if (n < 0) {
    n = -n;
    x = 1.0 / x;
  }

  for (; n; n >>= 1, x *= x) {
    if (n & 1)
      rv *= x;
  }

  return rv;
}

static
double
BLAS_dfpinfo(enum blas_cmach_type cmach) {
  double eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
  int t = 53, l = 1024, m = -1021;
  char rname[] = "BLAS_dfpinfo";

  if ((sizeof eps) == sizeof(float)) {
    t = 24;
    l = 128;
    m = -125;
  } else {
    t = 53;
    l = 1024;
    m = -1021;
  }

  /* for (i = 0; i < t; ++i) eps *= half; */
  eps = BLAS_dpow_di( b, -t );
  /* for (i = 0; i >= m; --i) r *= half; */
  r = BLAS_dpow_di( b, m-1 );

  o -= eps;
  /* for (i = 0; i < l; ++i) o *= b; */
  o = (o * BLAS_dpow_di( b, l-1 )) * b;

  switch (cmach) {
    case blas_eps: return eps;
    case blas_sfmin: return r;
    default:
      BLAS_error( rname, -1, cmach, 0 );
      break;
  }
  return 0.0;
}

static int check_orthogonality(int, int, int, PLASMA_Complex64_t*, double);
static int check_factorization(int, int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, PLASMA_Complex64_t*, double);
static int check_solution(int, int, int, PLASMA_Complex64_t*, int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, double);

int testing_zgels(int argc, char **argv)
{
    int mode = 0;

    if ( argc < 1 ){
        goto usage;
    } else {
        mode = atoi(argv[0]);
    }

    /* Check for number of arguments*/
    if ( ((mode == 0) && (argc != 6)) ||
         ((mode != 0) && (argc != 7)) ){
      usage:
        USAGE("GELS", "MODE M N LDA NRHS LDB [RH]",
              "   - MODE : 0: flat, 1: tree (RH needed)\n"
              "   - M    : number of rows of the matrix A\n"
              "   - N    : number of columns of the matrix A\n"
              "   - LDA  : leading dimension of the matrix A\n"
              "   - NRHS : number of RHS\n"
              "   - LDB  : leading dimension of the matrix B\n"
              "   - RH   : Size of each subdomains\n");
        return -1;
    }

    int M     = atoi(argv[1]);
    int N     = atoi(argv[2]);
    int LDA   = atoi(argv[3]);
    int NRHS  = atoi(argv[4]);
    int LDB   = atoi(argv[5]);
    int rh;

    int K = min(M, N);
    double eps;
    int info_ortho, info_solution, info_factorization;
    int i,j;
    int LDAxN = LDA*N;
    int LDBxNRHS = LDB*NRHS;

    PLASMA_Complex64_t *A1 = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *A2 = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B1 = (PLASMA_Complex64_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B2 = (PLASMA_Complex64_t *)malloc(LDB*NRHS*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Q  = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!B1)||(!B2)||(!Q)){
        printf("Out of Memory \n ");
        return -2;
    }

    if ( mode ) {
        rh = atoi(argv[6]);

        PLASMA_Set(PLASMA_HOUSEHOLDER_MODE, PLASMA_TREE_HOUSEHOLDER);
        PLASMA_Set(PLASMA_HOUSEHOLDER_SIZE, rh);
    }

    PLASMA_Alloc_Workspace_zgels(M, N, &T);
    eps = BLAS_dfpinfo( blas_eps );

    /*----------------------------------------------------------
    *  TESTING ZGELS
    */

    /* Initialize A1 and A2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i] ;

    /* Initialize B1 and B2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxNRHS, B1);
    for (i = 0; i < M; i++)
        for (j = 0; j < NRHS; j++)
            B2[LDB*j+i] = B1[LDB*j+i] ;

    memset((void*)Q, 0, LDA*N*sizeof(PLASMA_Complex64_t));
    for (i = 0; i < K; i++)
        Q[LDA*i+i] = 1.0;

    /* PLASMA ZGELS */
    PLASMA_zgels(PlasmaNoTrans, M, N, NRHS, A2, LDA, T, B2, LDB);

    /* PLASMA ZGELS */
    if (M >= N)
       /* Building the economy-size Q */
       PLASMA_zungqr(M, N, K, A2, LDA, T, Q, LDA);
    else
       /* Building the economy-size Q */
       PLASMA_zunglq(M, N, K, A2, LDA, T, Q, LDA);

    printf("\n");
    printf("------ TESTS FOR PLASMA ZGELS ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 60.\n");

    /* Check the orthogonality, factorization and the solution */
    info_ortho = check_orthogonality(M, N, LDA, Q, eps);
    info_factorization = check_factorization(M, N, A1, A2, LDA, Q, eps);
    info_solution = check_solution(M, N, NRHS, A1, LDA, B1, B2, LDB, eps);

    if ((info_solution == 0)&(info_factorization == 0)&(info_ortho == 0)) {
        printf("***************************************************\n");
        printf(" ---- TESTING ZGELS ...................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("************************************************\n");
        printf(" - TESTING ZGELS ... FAILED !\n");
        printf("************************************************\n");
    }

    /*-------------------------------------------------------------
    *  TESTING ZGEQRF + ZGEQRS or ZGELQF + ZGELQS
    */

    /* Initialize A1 and A2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i];

    /* Initialize B1 and B2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxNRHS, B1);
    for (i = 0; i < M; i++)
        for (j = 0; j < NRHS; j++)
             B2[LDB*j+i] = B1[LDB*j+i];

    memset((void*)Q, 0, LDA*N*sizeof(PLASMA_Complex64_t));
    for (i = 0; i < K; i++)
        Q[LDA*i+i] = 1.0;

    if (M >= N) {
        printf("\n");
        printf("------ TESTS FOR PLASMA ZGEQRF + ZGEQRS ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", M, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Plasma routines */
        PLASMA_zgeqrf(M, N, A2, LDA, T);
        PLASMA_zungqr(M, N, K, A2, LDA, T, Q, LDA);
        PLASMA_zgeqrs(M, N, NRHS, A2, LDA, T, B2, LDB);

        /* Check the orthogonality, factorization and the solution */
        info_ortho = check_orthogonality(M, N, LDA, Q, eps);
        info_factorization = check_factorization(M, N, A1, A2, LDA, Q, eps);
        info_solution = check_solution(M, N, NRHS, A1, LDA, B1, B2, LDB, eps);

        if ((info_solution == 0)&(info_factorization == 0)&(info_ortho == 0)) {
            printf("***************************************************\n");
            printf(" ---- TESTING ZGEQRF + ZGEQRS ............ PASSED !\n");
            printf("***************************************************\n");
        }
        else{
            printf("***************************************************\n");
            printf(" - TESTING ZGEQRF + ZGEQRS ... FAILED !\n");
            printf("***************************************************\n");
        }
    }
    else  {
        printf("\n");
        printf("------ TESTS FOR PLASMA ZGELQF + ZGELQS ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", M, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n", eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        /* Plasma routines */
        PLASMA_zgelqf(M, N, A2, LDA, T);
        PLASMA_zunglq(M, N, K, A2, LDA, T, Q, LDA);
        PLASMA_zgelqs(M, N, NRHS, A2, LDA, T, B2, LDB);

       /* Check the orthogonality, factorization and the solution */
       info_ortho = check_orthogonality(M, N, LDA, Q, eps);
       info_factorization = check_factorization(M, N, A1, A2, LDA, Q, eps);
       info_solution = check_solution(M, N, NRHS, A1, LDA, B1, B2, LDB, eps);

       if ( (info_solution == 0) & (info_factorization == 0) & (info_ortho == 0) ) {
          printf("***************************************************\n");
          printf(" ---- TESTING ZGELQF + ZGELQS ............ PASSED !\n");
          printf("***************************************************\n");
       }
       else {
          printf("***************************************************\n");
          printf(" - TESTING ZGELQF + ZGELQS ... FAILED !\n");
          printf("***************************************************\n");
        }
    }

    /*----------------------------------------------------------
    *  TESTING ZGEQRF + ZORMQR + ZTRSM
    */

    /* Initialize A1 and A2 */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i];

    /* Initialize B1 and B2 */
    memset(B2, 0, LDB*NRHS*sizeof(PLASMA_Complex64_t));
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxNRHS, B1);
    for (i = 0; i < M; i++)
        for (j = 0; j < NRHS; j++)
            B2[LDB*j+i] = B1[LDB*j+i];

    /* PLASMA ZGEQRF+ ZUNMQR + ZTRSM */
    memset((void*)Q, 0, LDA*N*sizeof(PLASMA_Complex64_t));
    for (i = 0; i < K; i++)
        Q[LDA*i+i] = 1.0;

    if (M >= N) {
        printf("\n");
        printf("------ TESTS FOR PLASMA ZGEQRF + ZUNMQR + ZTRSM  ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", M, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n",eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        PLASMA_zgeqrf(M, N, A2, LDA, T);
        PLASMA_zungqr(M, N, K, A2, LDA, T, Q, LDA);
        PLASMA_zunmqr(PlasmaLeft, PlasmaConjTrans, M, NRHS, N, A2, LDA, T, B2, LDB);
        PLASMA_ztrsm(PlasmaLeft, PlasmaUpper, PlasmaNoTrans, PlasmaNonUnit, N, NRHS, 1.0, A2, LDA, B2, LDB);
    }
    else {
        printf("\n");
        printf("------ TESTS FOR PLASMA ZGELQF + ZUNMLQ + ZTRSM  ROUTINE -------  \n");
        printf("            Size of the Matrix %d by %d\n", M, N);
        printf("\n");
        printf(" The matrix A is randomly generated for each test.\n");
        printf("============\n");
        printf(" The relative machine precision (eps) is to be %e \n",eps);
        printf(" Computational tests pass if scaled residuals are less than 60.\n");

        PLASMA_zgelqf(M, N, A2, LDA, T);
        PLASMA_ztrsm(PlasmaLeft, PlasmaLower, PlasmaNoTrans, PlasmaNonUnit, M, NRHS, 1.0, A2, LDA, B2, LDB);
        PLASMA_zunglq(M, N, K, A2, LDA, T, Q, LDA);
        PLASMA_zunmlq(PlasmaLeft, PlasmaConjTrans, N, NRHS, M, A2, LDA, T, B2, LDB);
    }

    /* Check the orthogonality, factorization and the solution */
    info_ortho = check_orthogonality(M, N, LDA, Q, eps);
    info_factorization = check_factorization(M, N, A1, A2, LDA, Q, eps);
    info_solution = check_solution(M, N, NRHS, A1, LDA, B1, B2, LDB, eps);

    if ( (info_solution == 0) & (info_factorization == 0) & (info_ortho == 0) ) {
        if (M >= N) {
            printf("***************************************************\n");
            printf(" ---- TESTING ZGEQRF + ZUNMQR + ZTRSM .... PASSED !\n");
            printf("***************************************************\n");
        }
        else {
            printf("***************************************************\n");
            printf(" ---- TESTING ZGELQF + ZTRSM + ZUNMLQ .... PASSED !\n");
            printf("***************************************************\n");
        }
    }
    else {
        if (M >= N) {
            printf("***************************************************\n");
            printf(" - TESTING ZGEQRF + ZUNMQR + ZTRSM ... FAILED !\n");
            printf("***************************************************\n");
        }
        else {
            printf("***************************************************\n");
            printf(" - TESTING ZGELQF + ZTRSM + ZUNMLQ ... FAILED !\n");
            printf("***************************************************\n");
        }
    }

    free(A1); free(A2); free(B1); free(B2); free(Q);
    PLASMA_Dealloc_Handle_Tile( &T );

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

static int check_orthogonality(int M, int N, int LDQ, PLASMA_Complex64_t *Q, double eps)
{
    double alpha, beta;
    double normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);

    double *work = (double *)malloc(minMN*sizeof(double));

    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    PLASMA_Complex64_t *Id = (PLASMA_Complex64_t *) malloc(minMN*minMN*sizeof(PLASMA_Complex64_t));
    memset((void*)Id, 0, minMN*minMN*sizeof(PLASMA_Complex64_t));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (PLASMA_Complex64_t)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_zherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    BLAS_zsy_norm( blas_colmajor, blas_inf_norm, blas_upper, minMN, Id, minMN, &normQ );

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n", normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || isinf(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 60.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

static int check_factorization(int M, int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, PLASMA_Complex64_t *Q, double eps )
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha, beta;
    int info_factorization;
    int i,j;

    PLASMA_Complex64_t *Ql       = (PLASMA_Complex64_t *)malloc(M*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(M*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(max(M,N)*sizeof(double));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        PLASMA_Complex64_t *R = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
        memset((void*)R, 0, N*N*sizeof(PLASMA_Complex64_t));
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex64_t));
        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, CBLAS_SADDR(alpha), Q, LDA, R, N, CBLAS_SADDR(beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        PLASMA_Complex64_t *L = (PLASMA_Complex64_t *)malloc(M*M*sizeof(PLASMA_Complex64_t));
        memset((void*)L, 0, M*M*sizeof(PLASMA_Complex64_t));
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(PLASMA_Complex64_t));
        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, CBLAS_SADDR(alpha), L, M, Q, LDA, CBLAS_SADDR(beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, Residual, M, &Rnorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, A2, LDA, &Anorm );

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || isinf(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 60.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}

/*--------------------------------------------------------------
 * Check the solution
 */

static int check_solution(int M, int N, int NRHS, PLASMA_Complex64_t *A, int LDA, PLASMA_Complex64_t *B, PLASMA_Complex64_t *X, int LDB, double eps)
{
    int info_solution;
    double Rnorm, Anorm, Xnorm, Bnorm;
    PLASMA_Complex64_t zone, mzone, zzero;
    double result;
    double *work = (double *)malloc(max(M, N)* sizeof(double));

    zone  =  1.0;
    mzone = -1.0;
    zzero =  0.0;

    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N,    A, LDA, &Anorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, NRHS, B, LDB, &Bnorm );
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, X, LDB, &Xnorm );

    /* Compute Ax - b */
    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, CBLAS_SADDR(zone), A, LDA, X, LDB, CBLAS_SADDR(mzone), B, LDB);

    /* Compute A' * (Ax - b) */
    cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, NRHS, M, CBLAS_SADDR(zone), A, LDA, B, LDB, CBLAS_SADDR(zzero), X, LDB);
    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, X, LDB, &Rnorm );

    if (getenv("PLASMA_TESTING_VERBOSE"))
       printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
     }
    else{
        printf("-- The solution is CORRECT ! \n");
        info_solution = 0;
    }
    free(work);
    return info_solution;
}


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

#pragma once

#include <gsl/gsl>

namespace _winml {

void LoadSpanFromDisjointBuffers(
    size_t num_buffers,
    std::function<gsl::span<byte>(size_t)> get_buffer,
    gsl::span<byte>& buffer_span);

void StoreSpanIntoDisjointBuffers(
    size_t num_buffers,
    std::function<gsl::span<byte>(size_t)> get_buffer,
    gsl::span<byte>& buffer_span);

} // namespace _winml

max_stars_repo_name: dennyac/onnxruntime

max_issues_repo_name: dennyac/onnxruntime

max_forks_repo_name: dennyac/onnxruntime


text:
/* Standalone c program to use viscdisk to run the Pringle (1981) ring
   problem, with a complex equation of state */

#include <math.h>
#include <stdio.h>
#include "init.h"
#include "vader.h"
#include "vader_common.h"
#include <gsl/gsl_sf_bessel.h>

/* Physical constants */
#define KB              1.3806488e-16
#define MH              1.6737236e-24
#define MU              0.61
#define A               7.5657314e-15

/* Problem parameters */

/* Set this to a non-empty value to restart from a checkpoint */
#define CHECKPOINT      "ringrad_00001.vader"

/* Grid parmeters */
#define NR              4096
#define RMIN            1.5e10
#define RMAX            1.5e12
#define LINEAR          1

/* EOS parameters */
#define EOS_FUNC        1
#define GAMMA           1.66666666667
#define DELTA           0.0

/* Viscosity parameters */
#define ALPHA_FUNC      1
#define ALPHA           0.0

/* Inner boundary condition */
#define IBC_PRES_TYPE   FIXED_TORQUE
#define IBC_ENTH_TYPE   FIXED_ENTHALPY_VALUE
#define IBC_FUNC        1
#define IBC_PRES_VAL    0.0
#define IBC_ENTH_VAL    0.0

/* Outer boundary condition */
#define OBC_PRES_TYPE   FIXED_TORQUE
#define OBC_ENTH_TYPE   FIXED_ENTHALPY_VALUE
#define OBC_FUNC        1
#define OBC_PRES_VAL    0.0
#define OBC_ENTH_VAL    0.0

/* Source functions */
#define MASS_SRC_FUNC   0
#define MASS_SRC_VAL    0.0
#define INT_EN_SRC_FUNC 0
#define INT_EN_SRC_VAL  0.0

/* Control parameters */
#define ERR_TOL         1.0e-6
#define DT_TOL          0.1
#define MAX_ITER        40
#define INTERP_ORDER    1
#define MAX_DT_INCREASE 1.5
#define DT_MIN          1.0e-20
#define MAX_STEP        100000
#define USE_BE          0
#define PRE_TIMESTEP    0
#define POST_TIMESTEP   0
#define VERBOSITY       1

/* Output parameters */
#define NSAVE            65
#define OUTFILE         "ringrad.out"
#define NUSEROUT        0
#define CHECKNAME       "ringrad"
#define USERREADCHK     false
#define USERWRITECHK    false

/* Time parameters */
#define DT_START        -1
#define END_TIME        0.128

/* Problem-specific parameters*/
#define MSTAR           (1.99e33*1.4)
#define RING_MASS       1.99e27
#define INIT_TEMP       1.0e4
#define COL_RATIO       1.0e10
#define RING_LOC        7.5e11
#define NU              5.93e8
#define FZ0             7.5e9
#define GAMMA_GAS       (5.0/3.0)

/* Main */
int main () {

  grid *grd;
  double col0, t0, colAnalyt, x, tau;
  double *col, *pres, *eInt;
  double *colOut, *presOut, *mBndOut, *eBndOut, *eIntOut, *tOut;
  double param[5];
  double tSave[NSAVE];
  bool writeCheckpoint[NSAVE];
  char *checkpoint = CHECKPOINT;
  int i, j, idx;
  FILE *fp;
  unsigned long nStep, nIter, nFail, nOut;

  /* Allocate grid and workspace */
  grd = gridInitKeplerian(NR, RMIN, RMAX, MSTAR, LINEAR);

  /* Allocate data */
  col = (double *) calloc(NR, sizeof(double));
  pres = (double *) calloc(NR, sizeof(double));
  eInt = (double *) calloc(NR, sizeof(double));

  /* Initialize column density and pressure arrays, and associated variables */
  col0 = RING_MASS / (M_PI*SQR(RING_LOC));
  t0 = SQR(RING_LOC)/(12.0*NU);
  for (i=0, idx=-1; i<NR; i++)
    if ((grd->r_h[i] < RING_LOC) && (RING_LOC <= grd->r_h[i+1])) idx = i;
  for (i=0; i<NR; i++)
    col[i] = RING_MASS / grd->area[idx] / COL_RATIO;
  col[idx] = RING_MASS / grd->area[idx];
  for (i=0; i<NR; i++) {
    pres[i] = col[i]*INIT_TEMP*KB/(MU*MH) + (1./3.)*FZ0*A*SQR(SQR(INIT_TEMP));
    eInt[i] = col[i]*INIT_TEMP*KB/(MU*MH*(GAMMA_GAS-1)) +
      FZ0*A*SQR(SQR(INIT_TEMP));
  }
  param[0] = NU;
  param[1] = RING_LOC;
  param[2] = RING_MASS;
  param[3] = RING_MASS / grd->area[idx] / COL_RATIO;
  param[4] = INIT_TEMP*KB/(MU*MH);

  /* Set output times */
  for (i=0; i<NSAVE; i++) {
    tSave[i] = t0 * i*END_TIME/((float) NSAVE-1);
    writeCheckpoint[i] = true;
  }

  /* Run the simulation */
  vader(
	/* Restart file name; if this is empty or NULL, calculation is
	   run from the start */
	checkpoint,
	/* Start and end time */
	0.0, END_TIME*t0,
	/* Equation of state parameters */
	EOS_FUNC, GAMMA, DELTA,
	/* Viscosity parameters */
	ALPHA_FUNC, ALPHA,
	/* Inner boundary condition */
	IBC_PRES_TYPE, IBC_ENTH_TYPE,
	IBC_FUNC, IBC_PRES_VAL, IBC_ENTH_VAL,
	/* Outer boundary condition */
	OBC_PRES_TYPE, OBC_ENTH_TYPE,
	OBC_FUNC, OBC_PRES_VAL, OBC_ENTH_VAL,
	/* Source functions */
	MASS_SRC_FUNC, MASS_SRC_VAL,
	INT_EN_SRC_FUNC, INT_EN_SRC_VAL,
	/* Control and method parameters */
	DT_START, DT_MIN, DT_TOL, ERR_TOL, MAX_DT_INCREASE,
	MAX_ITER, INTERP_ORDER, MAX_STEP, USE_BE,
	PRE_TIMESTEP, POST_TIMESTEP, VERBOSITY,
	/* Output control parameters */
	NSAVE, tSave, NUSEROUT, NULL,
	writeCheckpoint, CHECKNAME,
	USERREADCHK, USERWRITECHK,
	/* Computational grid */
	&grd,
	/* Initial conditions, holders for final conditions */
	&col, &pres, &eInt,
	/* User-defined extra parameters */
	&param,
	/* Diagnostic outputs */
	&nStep, &nIter, &nFail, &nOut,
	/* Simulation outputs */
	&tOut, &colOut, &presOut, &eIntOut, &mBndOut, &eBndOut,
	NULL, NULL, NULL
	);

  /* Print diagnostic output */
  printf("Total iterations = %ld, failed convergences = %ld\n",
	 nIter, nFail);

  /* Write output to file */
  if (!(fp=fopen(OUTFILE, "w"))) {
    fprintf(stderr, "Unable to open output file %s!\n", OUTFILE);
    exit(1);
  }
  fprintf(fp, "time/ts        x              col            colExact       pres           eInt\n");
  for (i=0; i<nOut; i++) {
    for (j=0; j<NR; j++) {
      x = grd->r_g[j+1]/RING_LOC;
      tau = tOut[i]/t0;
      colAnalyt = col0/ (pow(x, 0.25)*tau) * exp(-SQR(x-1.0)/tau) *
	gsl_sf_bessel_Inu_scaled(0.25, 2*x/tau);
      if (!isfinite(colAnalyt)) colAnalyt=0.0;
      if (colAnalyt < RING_MASS / grd->area[idx] / COL_RATIO)
	colAnalyt = RING_MASS / grd->area[idx] / COL_RATIO;
      fprintf(fp, "%e   %e   %e   %e   %e   %e\n",
	      tau, x,
	      colOut[grd->nr*i+j],
	      colAnalyt,
	      presOut[grd->nr*i+j],
	      eIntOut[grd->nr*i+j]);
    }
  }
  fclose(fp);

  /* Return successful exit */
  return(0);
}

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
/*
    This code is written by <albanese@fbk.it>.
    (C) 2008 mlpy Developers.

    See DWT in the GSL Library.

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

#include <Python.h>
#include <numpy/arrayobject.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_math.h>

static PyObject *dwt_dwt(PyObject *self, PyObject *args, PyObject *keywds)
{
  PyObject *x = NULL; PyObject *xcopy = NULL;

  char wf;
  int k, n;
  double *_xcopy;
  PyObject *centered = Py_False;
  gsl_wavelet *w;
  gsl_wavelet_workspace *work;

  /* Parse Tuple*/
  static char *kwlist[] = {"x", "wf", "k", "centered", NULL};
  if (!PyArg_ParseTupleAndKeywords(args, keywds, "Oci|O", kwlist,
				   &x, &wf, &k, &centered))
    return NULL;

  /* Build xcopy */
  xcopy = PyArray_FROM_OTF(x, NPY_DOUBLE, NPY_OUT_ARRAY | NPY_ENSURECOPY);
  if (xcopy == NULL) return NULL;

  n = (int) PyArray_DIM(xcopy, 0);
  _xcopy = (double *) PyArray_DATA(xcopy);

  switch (wf)
    {
    case 'd':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_daubechies, k);
      break;

    case 'h':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_haar_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_haar, k);
      break;

    case 'b':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_bspline_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_bspline, k);
      break;

    default:
      PyErr_SetString(PyExc_ValueError, "wavelet family is not valid");
      return NULL;
    }

  work = gsl_wavelet_workspace_alloc (n);

  gsl_wavelet_transform_forward (w, _xcopy, 1, n, work);

  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);

  return Py_BuildValue("N", xcopy);
}

static PyObject *dwt_idwt(PyObject *self, PyObject *args, PyObject *keywds)
{
  PyObject *x = NULL; PyObject *xcopy = NULL;

  char wf;
  int k, n;
  double *_xcopy;
  PyObject *centered = Py_False;
  gsl_wavelet *w;
  gsl_wavelet_workspace *work;

  /* Parse Tuple*/
  static char *kwlist[] = {"X", "wf", "k", "centered", NULL};
  if (!PyArg_ParseTupleAndKeywords(args, keywds, "Oci|O", kwlist,
				   &x, &wf, &k, &centered))
    return NULL;

  /* Build xcopy */
  xcopy = PyArray_FROM_OTF(x, NPY_DOUBLE, NPY_OUT_ARRAY | NPY_ENSURECOPY);
  if (xcopy == NULL) return NULL;

  n = (int) PyArray_DIM(xcopy, 0);
  _xcopy = (double *) PyArray_DATA(xcopy);

  switch (wf)
    {
    case 'd':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_daubechies_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_daubechies, k);
      break;

    case 'h':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_haar_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_haar, k);
      break;

    case 'b':
      if (centered == Py_True)
	w = gsl_wavelet_alloc (gsl_wavelet_bspline_centered, k);
      else
	w = gsl_wavelet_alloc (gsl_wavelet_bspline, k);
      break;

    default:
      PyErr_SetString(PyExc_ValueError, "wavelet family is not valid");
      return NULL;
    }

  work = gsl_wavelet_workspace_alloc (n);

  gsl_wavelet_transform_inverse (w, _xcopy, 1, n, work);

  gsl_wavelet_free (w);
  gsl_wavelet_workspace_free (work);

  return Py_BuildValue("N", xcopy);
}

/* Doc strings: */
static char module_doc[]  = "Discrete Wavelet Transform Module from GSL";

static char dwt_dwt_doc[] =
  "Discrete Wavelet Tranform\n\n"
  ":Parameters:\n"
  "   x : 1d array_like object (the length is restricted to powers of two)\n"
  "      data\n"
  "   wf : string ('d': daubechies, 'h': haar, 'b': bspline)\n"
  "      wavelet family\n"
  "   k : integer\n"
  "      member of the wavelet family\n\n"
  "         * daubechies : k = 4, 6, ..., 20 with k even\n"
  "         * haar : the only valid choice of k is k = 2\n"
  "         * bspline : k = 103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309\n"
  "   centered : bool\n"
  "      align the coefficients of the various sub-bands on edges.\n"
  "      Thus the resulting visualization of the coefficients of the\n"
  "      wavelet transform in the phase plane is easier to understand.\n\n"
  ":Returns:\n"
  "  X : 1d numpy array\n"
  "    discrete wavelet transformed data\n\n"
  "Example\n\n"
  ">>> import numpy as np\n"
  ">>> import mlpy.wavelet as wave\n"
  ">>> x = np.array([1,2,3,4,3,2,1,0])\n"
  ">>> wave.dwt(x=x, wf='d', k=6)\n"
  "array([ 5.65685425,  3.41458985,  0.29185347, -0.29185347, -0.28310081,\n"
  "       -0.07045258,  0.28310081,  0.07045258])\n";

static char dwt_idwt_doc[] =
  "Inverse Discrete Wavelet Tranform\n\n"
  ":Parameters:\n"
  "   X : 1d array_like object\n"
  "      discrete wavelet transformed data\n"
  "   wf : string ('d': daubechies, 'h': haar, 'b': bspline)\n"
  "      wavelet type\n"
  "   k : integer\n"
  "      member of the wavelet family\n\n"
  "         * daubechies : k = 4, 6, ..., 20 with k even\n"
  "         * haar : the only valid choice of k is k = 2\n"
  "         * bspline : k = 103, 105, 202, 204, 206, 208, 301, 303, 305 307, 309\n\n"
  "   centered : bool\n"
  "      if the coefficients are aligned\n\n"
  ":Returns:\n"
  "   x : 1d numpy array\n"
  "      data\n\n"
  "Example:\n\n"
  ">>> import numpy as np\n"
  ">>> import mlpy.wavelet as wave\n"
  ">>> X = np.array([ 5.65685425,  3.41458985,  0.29185347, -0.29185347, -0.28310081,\n"
  "...               -0.07045258,  0.28310081,  0.07045258])\n"
  ">>> wave.idwt(X=X, wf='d', k=6)\n"
  "array([  1.00000000e+00,   2.00000000e+00,   3.00000000e+00,\n"
  "         4.00000000e+00,   3.00000000e+00,   2.00000000e+00,\n"
  "         1.00000000e+00,  -3.53954610e-09])\n";

/* Method table */
static PyMethodDef dwt_methods[] = {
  {"dwt",
   (PyCFunction)dwt_dwt,
   METH_VARARGS | METH_KEYWORDS,
   dwt_dwt_doc},
  {"idwt",
   (PyCFunction)dwt_idwt,
   METH_VARARGS | METH_KEYWORDS,
   dwt_idwt_doc},
  {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3

static struct PyModuleDef moduledef = {
  PyModuleDef_HEAD_INIT,
  "_dwt",
  module_doc,
  -1,
  dwt_methods,
  NULL, NULL, NULL, NULL
};

PyObject *PyInit__dwt(void)
{
  PyObject *m;
  m = PyModule_Create(&moduledef);
  if (!m) {
    return NULL;
  }

  import_array();

  return m;
}

#else

PyMODINIT_FUNC init_dwt(void)
{
  PyObject *m;

  m = Py_InitModule3("_dwt", dwt_methods, module_doc);
  if (m == NULL) {
    return;
  }

  import_array();
}

#endif


text:
/*
  NAME:
     read_IC
  PURPOSE:
     read the initial conditions file
  CALLING SEQUENCE:
     read_IC(char ICfilename[])
  INPUT:
     ICfilename   - initial conditions filename
  OUTPUT:
     sets the options and the initial conditions
  REVISION HISTORY:
     2008-09-21 - Written Bovy
*/
#include <stdbool.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <proj_gauss_mixtures.h>
#include <proj_gauss_main.h>

bool read_IC(char ICfilename[]){

  FILE *ICfile;
  if ( (ICfile= fopen(ICfilename,"r")) == NULL){
    printf ("Opening the initial conditions file failed...\n");
    return false;
  }
  char line[100];
  printf("Reading the options section of the initial conditions file...\n");
  while (fgets(line,100,ICfile) != NULL){
    if (line[0] != '#'){
      if (line[0] == '\n') break;
      if (parse_option(line) == false){
	printf("One of the lines in the options section of the initial conditions file is corrupted\n");
	printf("Please check the initial conditions file and try again\n");
	return false;
      }
    }
  }

  printf("Successfully read the options\n");

  printf("Reading the initial model parameters...\n");

  //Read first block, establish the dimension d of the modeled quantities
  int countd=0;
  double *mmtemp = (double *) malloc (1000000 * sizeof (double));
  while (fgets(line,100,ICfile) != NULL){
    if (line[0] != '#'){
      if (line[0] == '\n') break;
      *(mmtemp++) = atof(line);
      ++countd;
    }
  }
  //now determine d
  d = (int) (-3 + sqrt(9 + 8 * (countd-1)))/2 ;
  dV = (int) (d*(d+1)/2);

  //allocate the alpha, mm and VV matrices
  int kk;
  for (kk=0; kk != K; ++kk){
    gaussians->mm = gsl_vector_alloc (d);
    gaussians->VV = gsl_matrix_alloc (d,d);;
    ++gaussians;
  }
  gaussians -= K;

  //first map the mmtemp values on the right alpha, mm, VV
  mmtemp -= countd;
  (*gaussians).alpha = *(mmtemp++);
  int dd;
  for (dd=0; dd != d; ++dd)
    gsl_vector_set(gaussians->mm,dd,*(mmtemp++));
  int dd1,dd2;
  for (dd1=0; dd1 != d; ++dd1)
    gsl_matrix_set(gaussians->VV,dd1,dd1,*(mmtemp++));
  for (dd1=0; dd1 != d-1; ++dd1)
    for (dd2=dd1+1; dd2 != d; ++dd2){
      gsl_matrix_set(gaussians->VV,dd1,dd2,*mmtemp);
      gsl_matrix_set(gaussians->VV,dd2,dd1,*mmtemp);
      mmtemp++;
    }

  ++gaussians;

  //reallocate mmtemp
  mmtemp -= countd;
  mmtemp = (double *) realloc (mmtemp,countd * sizeof (double) );
  if (mmtemp == NULL){
    printf("Error reallocating memory\n");
    printf("Returning\n");
    return false;
  }

  //Then read the rest of the Gaussians.
  for (kk=1; kk != K; ++kk){
    while (fgets(line,100,ICfile) != NULL){
      if (line[0] != '#'){
	if (line[0] == '\n') break;
	*(mmtemp++) = atof(line);
      }
    }
    mmtemp -=countd;
    (*gaussians).alpha = *(mmtemp++);
    for (dd=0; dd != d; ++dd)
      gsl_vector_set(gaussians->mm,dd,*(mmtemp++));
    for (dd1=0; dd1 != d; ++dd1)
      gsl_matrix_set(gaussians->VV,dd1,dd1,*(mmtemp++));
    for (dd1=0; dd1 != d-1; ++dd1)
      for (dd2=dd1+1; dd2 != d; ++dd2){
	gsl_matrix_set(gaussians->VV,dd1,dd2,*mmtemp);
	gsl_matrix_set(gaussians->VV,dd2,dd1,*mmtemp);
	mmtemp++;
      }
    ++gaussians;
    mmtemp -= countd;
  }
  gaussians -= K;

  free(mmtemp);

  fclose(ICfile);

  printf("Successfully read initial model parameters from the initial conditions file\n");

  //Print options
  printf("\nThe options are set to:\n");
  printf("K\t\t=\t");
  printf("%i",K);
  printf("\n");
  printf("maxiter\t\t=\t");
  printf("%lli",maxiter);
  printf("\n");
  printf("tol\t\t=\t");
  printf("%f",tol);
  printf("\n");
  printf("splitnmerge\t=\t");
  printf("%i",splitnmerge);
  printf("\n");
  printf("likeonly\t=\t");
  printf("%i",likeonly);
  printf("\n");
  printf("w\t\t=\t");
  printf("%f",w);
  printf("\n");
  printf("fixamp\t\t=\t");
  int ii;
  for (ii=0; ii != K; ++ii){
    printf("%i",*fixampP);
    if (ii < K-1) printf("\t");
    fixampP++;
  }
  fixampP -= K;
  printf("\n");
  printf("fixmean\t\t=\t");
  for (ii=0; ii != K; ++ii){
    printf("%i",*fixmeanP);
    if (ii < K-1) printf("\t");
    fixmeanP++;
  }
  fixmeanP -= K;
  printf("\n");
  printf("fixcovar\t=\t");
  for (ii=0; ii != K; ++ii){
    printf("%i",*fixcovarP);
    if (ii < K-1) printf("\t");
    fixcovarP++;
  }
  fixcovarP -= K;
  printf("\n");

  //Print the initial model parameters
  printf("\nInitial model parameters used:\n\n");
  for (kk=0; kk != K; ++kk){
    printf("Gaussian ");
    printf("%i",kk);
    printf("\n");
    printf("amp\t=\t");
    printf("%f",(*gaussians).alpha);
    printf("\n");
    printf("mean\t=\t");
    for (dd=0; dd != d; ++dd){
      printf("%f",gsl_vector_get(gaussians->mm,dd));
      if (dd < d-1) printf("\t");
    }
    printf("\n");
    printf("covar\t=\t");
    for (dd1=0; dd1 != d; ++dd1)
      printf("%f\t",gsl_matrix_get(gaussians->VV,dd1,dd1));
    for (dd1=0; dd1 != d-1; ++dd1)
      for (dd2=dd1+1; dd2 != d; ++dd2){
	printf("%f\t",gsl_matrix_get(gaussians->VV,dd1,dd2));
      }
    ++gaussians;
    printf("\n\n");
  }
  gaussians -= K;

  return true;

}