text:
#include <pthread.h>
#include <stdlib.h>
#include <string.h>

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

// MAX and MIN aren't defined in the standard library. O:
#define MAX(x, y) (((x) > (y)) ? (x) : (y))
#define MIN(x, y) (((x) < (y)) ? (x) : (y))

typedef struct PopulationParametersStruct {
    int n; ///< Population size, used for determining the number of generations between events.
    float f; ///< Scaling factor for the total mutation rate. Smaller = shorter time to coalescence.
    float c; ///< Constant bias for the upward mutation rate.
    float d; ///< Linear bias for the downward mutation rate.
    int kappa; ///< Lower bound of repeat lengths.
    int omega; ///< Upper bound of repeat lengths.
} PopulationParameters;

typedef struct CoalescentStruct {
    PopulationParameters theta; ///< Mutation model parameters.
    int *coalescent_tree; ///< Pointer to our tree, stored as an array.
    int offset; ///< We generalize to include 1+ ancestors. Determine the offset for array representation of tree.
    gsl_rng *r; ///< Pointer to our RNG. This is preserved across the trace and evolve steps.
} PopulationTree;

int _triangle (int a) { return (int) (a * (a + 1) / 2.0); }
int _round_num (int a) { return (a < 0) ? (int) (a - 0.5) : (int) (a + 0.5); }
typedef int (*_mutate_f) (int, int, float, float, int, int, const gsl_rng *);

int _mutate_generation (int t, int ell, float c, float d, int kappa, int omega, const gsl_rng *r) {
    for (int k = 0; k < t; k++) {
        // Compute our upward mutation rate. We are bounded by omega.
        ell = (gsl_ran_flat(r, 0, 1) < c) ? MIN(omega, ell + 1) : ell;

        // Compute our downward mutation rate. We are bounded by kappa.
        ell = (gsl_ran_flat(r, 0, 1) < ell * d) ? MAX(kappa, ell - 1) : ell;
    }

    // Return the result of mutating over several generations.
    return ell;
}

int _mutate_draw (int t, int ell, float c, float d, int kappa, int omega, const gsl_rng *r) {
    return MIN((unsigned) omega, MAX((unsigned) kappa, ell + gsl_ran_poisson(r, t * c) -
                                                       gsl_ran_poisson(r, t * ell * d)));
}

void _trace_tree (int *coalescent_tree, int n, const gsl_rng *r) {
    for (int tau = 0; tau < 2 * n - 1; tau++) { // Diploid!
        int tau_0 = _triangle(tau + 0), tau_1 = _triangle(tau + 1); // Start at 2nd coalescence.

        // We save the indices of our ancestors to our descendants.
        for (int k = 0; k < tau_1 - tau_0; k++) {
        	coalescent_tree[tau_1 + k + 1] = tau_0 + k;
        }

        // Determine the individual whose frequency increases.
        gsl_ran_choose(r, &coalescent_tree[tau_1], 1, &coalescent_tree[tau_1 + 1], tau_1 - tau_0, sizeof(int));

        // Finally, shuffle the individuals in this generation.
        gsl_ran_shuffle(r, &coalescent_tree[tau_1], tau_1 - tau_0 + 1, sizeof(int));
    }
}

void _evolve_individual (int k, int *descendants, int t_coalescence, _mutate_f mutate, PopulationTree *p) {
    int *descendant_to_evolve = &(p->coalescent_tree[descendants[k]]);

    // Evolve each ancestor according to the average time to coalescence and the scaling factor f.
    *descendant_to_evolve = (*mutate)(t_coalescence, *descendant_to_evolve, p->theta.c, p->theta.d,
                                      p->theta.kappa, p->theta.omega, p->r);

    // Save our descendant state.
    descendants[k] = *descendant_to_evolve;
}

/**
 * We assumed our tree has been traced. Given the population tree structure and a pointer to a mutation function,
 * determine the repeat length of all individuals in our tree. We do so by iterating through each coalescent event.
 */
void _evolve_event (PopulationTree *p, _mutate_f mutate) {
    int descendants_size, expected_time, t_coalescence;
    int *descendants = NULL;

    for (int tau = p->offset; tau < 2 * p->theta.n - 1; tau++) {
        // We define our descendants for the tau'th coalescent.
        descendants = &(p->coalescent_tree[_triangle(tau + 1)]);
        descendants_size = _triangle(tau + 2) - _triangle(tau + 1);

        // Determine time to coalescence. This is exponentially distributed, but the mean stays the same. Scale by f.
        expected_time = (int) (p->theta.f * 2 * p->theta.n / (float) _triangle(tau + 1));
        t_coalescence = MAX(1, _round_num(gsl_ran_exponential(p->r, expected_time)));

        // Iterate through each of the descendants (currently indices) and determine each ancestor.
        for (int k = 0; k < descendants_size; k++) {
            _evolve_individual(k, descendants, t_coalescence, mutate, p);
        }
    }
}

void _evolve (int *i_0, int i_0_size, PopulationTree *p) {
    // Determine our offset, and seed our ancestors for the tree.
    p->offset = i_0_size - 1;
    for (int k = 0; k < i_0_size; k++) {
        p->coalescent_tree[k + _triangle(p->offset)] = i_0[k];
    }

    // From our common ancestors, descend forward in time and populate our tree with repeat lengths.
    _evolve_event(p, _mutate_draw);
}

void _cleanup (PopulationTree *p) {
    free(p->coalescent_tree);
    gsl_rng_free(p->r);
}

max_stars_repo_name: glennga/kumulaau

max_issues_repo_name: glennga/kumulaau

max_forks_repo_name: glennga/kumulaau


text:
#ifndef VECTORS_H
#define VECTORS_H

#include <cmath>
#include <gsl/gsl_blas.h>
#include "types.h"
#include "cuda.h"

__host__ __device__ void angle_ab(const double dir1[3], const double vel[3], double dir2[3]);
__host__ __device__ double doppler(const double dir1[3], const double vel[3]);
__host__ __device__ void scatter_dir(const double dir_in[3], double cos_theta, double dir_out[3]);
__host__ __device__ void get_rand_isotropic_unitvec(double vecout[3]);
__host__ __device__ void move_pkt(PKT *pkt_ptr, double distance, const double time);

__host__ __device__
inline double vec_len(const double x[3])
// return the the magnitude of a vector
{
  #ifdef __CUDA_ARCH__
    return sqrt((x[0] * x[0]) + (x[1] * x[1]) + (x[2] * x[2]));
  #else
    return cblas_dnrm2(3, x, 1);
  #endif
}

__host__ __device__
inline void vec_norm(const double vec_in[3], double vec_out[3])
// Routine for normalizing a vector.
{
  const double magnitude = vec_len(vec_in);

  vec_out[0] = vec_in[0] / magnitude;
  vec_out[1] = vec_in[1] / magnitude;
  vec_out[2] = vec_in[2] / magnitude;

  // vec_copy(vec_out, vec_in);
  // cblas_dscal(3, 1 / magnitude, vec_out, 1)
}

__host__ __device__
inline double dot(const double x[3], const double y[3])
// Routine for taking dot product.
{
  #ifdef __CUDA_ARCH__
    return (x[0] * y[0]) + (x[1] * y[1]) + (x[2] * y[2]);
  #else
    return cblas_ddot(3, x, 1, y, 1);
  #endif
}

__host__ __device__
inline void get_velocity(const double x[3], double y[3], const double t)
// Routine for getting velocity vector of the flow at a position with homologous expansion.
{
  y[0] = x[0] / t;
  y[1] = x[1] / t;
  y[2] = x[2] / t;
}

__host__ __device__
inline void cross_prod(const double vec1[3], const double vec2[3], double vecout[3])
{
  vecout[0] = (vec1[1] * vec2[2]) - (vec2[1] * vec1[2]);
  vecout[1] = (vec1[2] * vec2[0]) - (vec2[2] * vec1[0]);
  vecout[2] = (vec1[0] * vec2[1]) - (vec2[0] * vec1[1]);
}

__host__ __device__
inline void vec_scale(double vec[3], const double scalefactor)
{
  #ifdef __CUDA_ARCH__
    for (int d = 0; d < 3; d++)
    {
      vec[d] *= scalefactor;
    }
  #else
    cblas_dscal(3, scalefactor, vec, 1);
  #endif
}

__host__ __device__
inline void vec_copy(double destination[3], const double source[3])
{
  #ifdef __CUDA_ARCH__
    for (int d = 0; d < 3; d++)
    {
      destination[d] = source[d];
    }
  #else
    cblas_dcopy(3, source, 1, destination, 1);
  #endif
}

__host__ __device__
inline double doppler_packetpos(const PKT *const pkt_ptr)
{
  double vel_vec[3];
  get_velocity(pkt_ptr->pos, vel_vec, pkt_ptr->prop_time); // homologous flow velocity
  const double dopplerfactor = doppler(pkt_ptr->dir, vel_vec);
  return dopplerfactor;
}

#endif //VECTORS_H

ext: h

max_stars_repo_name: artis-mcrt/artis

max_issues_repo_name: artis-mcrt/artis

max_forks_repo_name: artis-mcrt/artis


text:
#ifndef UTIL_H
#define UTIL_H

#define PRECISION 2
#define MAX_LINE_LENGTH 9046
#define USE_CBLAS
#define USE_LAPACK

#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/time.h>
#include <time.h>

#ifdef USE_CBLAS
#include <cblas.h>
#endif

#ifdef USE_LAPACK
#include <lapacke.h>
#endif

/******************************************************************************
 *                                LOGGING
 ******************************************************************************/

/* DEBUG */
#ifdef NDEBUG
#define DEBUG(M, ...)
#else
#define DEBUG(M, ...) \
  fprintf(stderr, "[DEBUG] %s:%d: " M "\n", __func__, __LINE__, ##__VA_ARGS__)
#endif

/* LOG */
#define LOG_ERROR(M, ...) \
  fprintf(stderr, "[ERROR] [%s] " M "\n", __func__, ##__VA_ARGS__)
#define LOG_WARN(M, ...) fprintf(stderr, "[WARN] " M "\n", ##__VA_ARGS__)
#define LOG_INFO(M, ...) fprintf(stderr, "[INFO] " M "\n", ##__VA_ARGS__)

/* FATAL */
#define FATAL(M, ...)                                \
  fprintf(stderr, "[FATAL] " M "\n", ##__VA_ARGS__); \
  exit(-1);

/* CHECK */
#define CHECK(A, M, ...)       \
  if (!(A)) {                  \
    log_err(M, ##__VA_ARGS__); \
    goto error;                \
  }

/******************************************************************************
 *                                   DATA
 ******************************************************************************/

#if PRECISION == 1
typedef float real_t;
#elif PRECISION == 2
typedef double real_t;
#else
#error "Precision not defined!"
#endif

char *malloc_string(const char *s);
int dsv_rows(const char *fp);
int dsv_cols(const char *fp, const char delim);
char **dsv_fields(const char *fp, const char delim, int *nb_fields);
real_t **dsv_data(const char *fp, const char delim, int *nb_rows, int *nb_cols);
real_t **csv_data(const char *fp, int *nb_rows, int *nb_cols);
int **load_iarrays(const char *csv_path, int *nb_arrays);
real_t **load_darrays(const char *csv_path, int *nb_arrays);
real_t *load_vector(const char *file_path);

/******************************************************************************
 *                                 MATHS
 ******************************************************************************/

#define UNUSED(expr) \
  do {               \
    (void)(expr);    \
  } while (0)

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#define MIN(x, y) ((x) < (y) ? (x) : (y))
#define MAX(x, y) ((x) > (y) ? (x) : (y))
#define SIGN(a, b) ((b) >= 0.0 ? fabs(a) : -fabs(a))

float randf(float a, float b);
real_t deg2rad(const real_t d);
real_t rad2deg(const real_t r);
int fltcmp(const real_t x, const real_t y);
real_t pythag(const real_t a, const real_t b);
real_t lerp(const real_t a, const real_t b, const real_t t);
void lerp3(const real_t *a, const real_t *b, const real_t t, real_t *x);
real_t sinc(const real_t x);

/******************************************************************************
 *                              LINEAR ALGEBRA
 ******************************************************************************/

void print_matrix(const char *prefix, const real_t *data, const size_t m,
                  const size_t n);
void print_vector(const char *prefix, const real_t *data, const size_t length);

void eye(real_t *A, const size_t m, const size_t n);
void ones(real_t *A, const size_t m, const size_t n);
void zeros(real_t *A, const size_t m, const size_t n);

real_t *mat_new(const size_t m, const size_t n);
int mat_cmp(const real_t *A, const real_t *B, const size_t m, const size_t n);
int mat_equals(const real_t *A, const real_t *B, const size_t m, const size_t n,
               const real_t tol);
int mat_save(const char *save_path, const real_t *A, const int m, const int n);
real_t *mat_load(const char *save_path, int *nb_rows, int *nb_cols);
void mat_set(real_t *A, const size_t stride, const size_t i, const size_t j,
             const real_t val);
real_t mat_val(const real_t *A, const size_t stride, const size_t i,
               const size_t j);
void mat_copy(const real_t *src, const int m, const int n, real_t *dest);
void mat_block_get(const real_t *A, const size_t stride, const size_t rs,
                   const size_t cs, const size_t re, const size_t ce,
                   real_t *block);
void mat_block_set(real_t *A, const size_t stride, const size_t rs,
                   const size_t cs, const size_t re, const size_t ce,
                   const real_t *block);
void mat_diag_get(const real_t *A, const int m, const int n, real_t *d);
void mat_diag_set(real_t *A, const int m, const int n, const real_t *d);
void mat_triu(const real_t *A, const size_t n, real_t *U);
void mat_tril(const real_t *A, const size_t n, real_t *L);
real_t mat_trace(const real_t *A, const size_t m, const size_t n);
void mat_transpose(const real_t *A, size_t m, size_t n, real_t *A_t);
void mat_add(const real_t *A, const real_t *B, real_t *C, size_t m, size_t n);
void mat_sub(const real_t *A, const real_t *B, real_t *C, size_t m, size_t n);
void mat_scale(real_t *A, const size_t m, const size_t n, const real_t scale);

real_t *vec_new(const size_t length);
void vec_copy(const real_t *src, const size_t length, real_t *dest);
int vec_equals(const real_t *x, const real_t *y, const size_t length);
void vec_add(const real_t *x, const real_t *y, real_t *z, size_t length);
void vec_sub(const real_t *x, const real_t *y, real_t *z, size_t length);
void vec_scale(real_t *x, const size_t length, const real_t scale);
real_t vec_norm(const real_t *x, const size_t length);

void dot(const real_t *A, const size_t A_m, const size_t A_n, const real_t *B,
         const size_t B_m, const size_t B_n, real_t *C);
void skew(const real_t x[3], real_t A[3 * 3]);
void fwdsubs(const real_t *L, const real_t *b, real_t *y, const size_t n);
void bwdsubs(const real_t *U, const real_t *y, real_t *x, const size_t n);
int check_jacobian(const char *jac_name, const real_t *fdiff, const real_t *jac,
                   const size_t m, const size_t n, const real_t tol,
                   const int print);

#ifdef USE_CBLAS
void cblas_dot(const real_t *A, const size_t A_m, const size_t A_n,
               const real_t *B, const size_t B_m, const size_t B_n, real_t *C);
#endif

/******************************************************************************
 *                                   SVD
 ******************************************************************************/

int svd(real_t *A, int m, int n, real_t *U, real_t *s, real_t *V_t);
int svdcomp(real_t *A, int m, int n, real_t *w, real_t *V);
int pinv(real_t *A, const int m, const int n, real_t *A_inv);

#ifdef USE_LAPACK

#endif

/******************************************************************************
 *                                   CHOL
 ******************************************************************************/

void chol(const real_t *A, const size_t n, real_t *L);
void chol_solve(const real_t *A, const real_t *b, real_t *x, const size_t n);

#ifdef USE_LAPACK
void lapack_chol_solve(const real_t *A, const real_t *b, real_t *x,
                       const size_t n);
#endif

/******************************************************************************
 *                                   TIME
 ******************************************************************************/

typedef uint64_t timestamp_t;

struct timespec tic();
float toc(struct timespec *tic);
float mtoc(struct timespec *tic);
float time_now();

/******************************************************************************
 *                                TRANSFORMS
 ******************************************************************************/

void tf(const real_t params[7], real_t T[4 * 4]);
void tf_params(const real_t T[4 * 4], real_t params[7]);
void tf_rot_set(real_t T[4 * 4], const real_t C[3 * 3]);
void tf_trans_set(real_t T[4 * 4], const real_t r[3]);
void tf_trans_get(const real_t T[4 * 4], real_t r[3]);
void tf_rot_get(const real_t T[4 * 4], real_t C[3 * 3]);
void tf_quat_get(const real_t T[4 * 4], real_t q[4]);
void tf_inv(const real_t T[4 * 4], real_t T_inv[4 * 4]);
void tf_point(const real_t T[4 * 4], const real_t p[3], real_t retval[3]);
void tf_hpoint(const real_t T[4 * 4], const real_t p[4], real_t retval[4]);
void tf_perturb_rot(real_t T[4 * 4], const real_t step_size, const int i);
void tf_perturb_trans(real_t T[4 * 4], const real_t step_size, const int i);
void rvec2rot(const real_t *rvec, const real_t eps, real_t *R);
void euler321(const real_t euler[3], real_t C[3 * 3]);
void rot2quat(const real_t C[3 * 3], real_t q[4]);
void quat2euler(const real_t q[4], real_t euler[3]);
void quat2rot(const real_t q[4], real_t C[3 * 3]);
void quat_lmul(const real_t p[4], const real_t q[4], real_t r[4]);
void quat_rmul(const real_t p[4], const real_t q[4], real_t r[4]);
void quat_mul(const real_t p[4], const real_t q[4], real_t r[4]);
void quat_delta(const real_t dalpha[3], real_t dq[4]);

/******************************************************************************
 *                                 IMAGE
 ******************************************************************************/

typedef struct image_t {
  int width;
  int height;
  uint8_t *data;
} image_t;

void image_init(image_t *img, uint8_t *data, int width, int height);

/******************************************************************************
 *                                  CV
 ******************************************************************************/

/******************************** RADTAN **************************************/

void radtan4_distort(const real_t params[4], const real_t p[2], real_t p_d[2]);
void radtan4_point_jacobian(const real_t params[4], const real_t p[2],
                            real_t J_point[2 * 2]);
void radtan4_params_jacobian(const real_t params[4], const real_t p[2],
                             real_t J_param[2 * 4]);

/********************************* EQUI ***************************************/

void equi4_distort(const real_t params[4], const real_t p[2], real_t p_d[2]);
void equi4_point_jacobian(const real_t params[4], const real_t p[2],
                          real_t J_point[2 * 2]);
void equi4_params_jacobian(const real_t params[4], const real_t p[2],
                           real_t J_param[2 * 4]);

/******************************** PINHOLE *************************************/

real_t pinhole_focal(const int image_width, const real_t fov);
int pinhole_project(const real_t params[4], const real_t p_C[3], real_t x[2]);

void pinhole_point_jacobian(const real_t params[4], real_t J_point[2 * 3]);
void pinhole_params_jacobian(const real_t params[4], const real_t x[2],
                             real_t J[2 * 4]);

/**************************** PINHOLE-RADTAN4 *********************************/

void pinhole_radtan4_project(const real_t params[8], const real_t p_C[3],
                             real_t x[2]);
void pinhole_radtan4_project_jacobian(const real_t params[8],
                                      const real_t p_C[3], real_t J[2 * 3]);
void pinhole_radtan4_params_jacobian(const real_t params[8],
                                     const real_t p_C[3], real_t J[2 * 8]);

/***************************** PINHOLE-EQUI4 **********************************/

void pinhole_equi4_project(const real_t params[8], const real_t p_C[3],
                           real_t x[2]);
void pinhole_equi4_project_jacobian(const real_t params[8], const real_t p_C[3],
                                    real_t J[2 * 3]);
void pinhole_equi4_params_jacobian(const real_t params[8], const real_t p_C[3],
                                   real_t J[2 * 8]);

#endif  // ZERO_H

max_stars_repo_name: daoran/ba

max_issues_repo_name: daoran/ba

max_forks_repo_name: daoran/ba


text:
/* Copyright (C) 2015 Atsushi Togo */
/* All rights reserved. */

/* This file is part of phonopy. */

/* Redistribution and use in source and binary forms, with or without */
/* modification, are permitted provided that the following conditions */
/* are met: */

/* * Redistributions of source code must retain the above copyright */
/*   notice, this list of conditions and the following disclaimer. */

/* * 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. */

/* * Neither the name of the phonopy project 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. */

#ifndef __frequency_shift4_H__
#define __frequency_shift4_H__

#include <lapacke.h>
#include <phonoc_array.h>

void get_fc4_frequency_shifts(double *frequency_shifts,
			      const double *fc4_normal_real,
			      const double *frequencies,
			      const Iarray *grid_points1,
			      const Darray *temperatures,
			      const int *band_indicies,
			      const int num_band0,
			      const int num_band,
			      const double unit_conversion_factor);
void
get_fc4_normal_for_frequency_shift(double *fc4_normal_real,
				   const double *frequencies,
				   const lapack_complex_double *eigenvectors,
				   const int grid_point0,
				   const Iarray *grid_points1,
				   const int *grid_address,
				   const int *mesh,
				   const double *fc4,
				   const Darray *shortest_vectors,
				   const Iarray *multiplicity,
				   const double *masses,
				   const int *p2s_map,
				   const int *s2p_map,
				   const Iarray *band_indicies,
				   const double cutoff_frequency);
void reciprocal_to_normal4(lapack_complex_double *fc4_normal,
			   const lapack_complex_double *fc4_reciprocal,
			   const double *freqs0,
			   const double *freqs1,
			   const lapack_complex_double *eigvecs0,
			   const lapack_complex_double *eigvecs1,
			   const double *masses,
			   const int *band_indices,
			   const int num_band0,
			   const int num_band,
			   const double cutoff_frequency);
void set_phonons_for_frequency_shift(Darray *frequencies,
				     Carray *eigenvectors,
				     char *phonon_done,
				     const Iarray *grid_points,
				     const int *grid_address,
				     const int *mesh,
				     const Darray *fc2,
				     const Darray *svecs_fc2,
				     const Iarray *multi_fc2,
				     const double *masses_fc2,
				     const int *p2s_fc2,
				     const int *s2p_fc2,
				     const double unit_conversion_factor,
				     const double *born,
				     const double *dielectric,
				     const double *reciprocal_lattice,
				     const double *q_direction,
				     const double nac_factor,
				     const char uplo);

#endif

max_stars_repo_name: atztogo/forcefit

max_issues_repo_name: atztogo/forcefit

max_forks_repo_name: atztogo/forcefit


text:
/*
 obtainerr2_par.c
 written by Eunseok Lee

 function: solve a linear equation system and obtain the squared error

 v1: Feb 2, 2018
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//#include <cem.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

double obtainerr2_par(int row_ini, int row_end, double *matA, int n_row, int n_col, int *col_ids, int n_col_ids, double *y) {

    int i, j, k, tmp;
    double err[n_row];
    double sum;

    gsl_matrix * A = gsl_matrix_alloc (n_row-1,n_col_ids);
    gsl_matrix * V = gsl_matrix_alloc (n_col_ids,n_col_ids);
    gsl_vector * S = gsl_vector_alloc (n_col_ids);
    gsl_vector * work = gsl_vector_alloc (n_col_ids);
//    gsl_matrix * U = gsl_matrix_alloc (n_row-1,n_col_ids);
    gsl_vector * b = gsl_vector_alloc (n_row-1);
    gsl_vector * x = gsl_vector_alloc (n_col_ids);

    for (i=0;i<n_row;i++)
        err[i] = 0.0;

    // need to copy y to b!
    for (i=row_ini;i<row_end;i++) {
        for (j=0;j<i;j++) {
            for (k=0; k<n_col_ids; k++) {
                tmp = *(col_ids+k);
                gsl_matrix_set (A, j, k, *(matA+j*n_col+tmp));
            }
            gsl_vector_set (b, j, y[j]);
        }
        for (j=i+1;j<n_row;j++) {
            for (k=0; k<n_col_ids; k++) {
                tmp = *(col_ids+k);
                gsl_matrix_set (A, j-1, k, *(matA+j*n_col+tmp));
            }
            gsl_vector_set (b, j-1, y[j]);
        }

        gsl_linalg_SV_decomp(A, V, S, work);    // on output, A is replaced by U!
        gsl_linalg_SV_solve (A, V, S, b, x);    // So, A should be used instead of U, here.

        sum = 0.0;
        for (k=0;k<n_col_ids;k++) {
            sum += *(matA+i*n_col+col_ids[k]) * gsl_vector_get(x,k);
        }
        err[i] = y[i] - sum;
    }

    sum = 0.0;
    for (i=row_ini;i<row_end;i++)
        sum += err[i]*err[i];

    gsl_matrix_free(A);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(b);
    gsl_vector_free(x);

    return(sum);
}

ext: c


text:
/* -*- c++ -*- */
/*
 * Copyright (C) 2017 Free Software Foundation, Inc.
 *
 * This file is part of GNU Radio
 *
 * GNU Radio is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3, or (at your option)
 * any later version.
 *
 * GNU Radio 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 GNU Radio; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street,
 * Boston, MA 02110-1301, USA.
 */

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>

#define	NSTEPS		 10	// how many steps of mu are in the generated table
#define	MAX_NSTEPS	256
#define	NTAPS		  8	// # of taps in the interpolator
#define	MAX_NTAPS	128

extern void initpt (double x[], int ntaps);
extern double objective (double x[], int *ntaps);
extern double global_mu;
extern double global_B;
extern gsl_integration_workspace *global_gsl_int_workspace;

// fortran
extern double prax2_ (double (fct)(double x[], int *ntaps),
		      double initv[], int *ntaps, double result[]);

static void
usage (char *name)
{
  fprintf (stderr, "usage: %s [-v] [-n <nsteps>] [-t <ntaps>] [-B <bw>]\n", name);
  exit (1);
}

static void
printline (double x[], int ntaps, int imu, int nsteps)
{
  int	i;

  printf ("  { ");
  for (i = 0; i < ntaps; i++){
    printf ("%12.5e", x[i]);
    if (i != ntaps - 1)
      printf (", ");
    else
      printf (" }, // %3d/%d\n", imu, nsteps);
  }
}

int
main (int argc, char **argv)
{
  double	xx[MAX_NSTEPS+1][MAX_NTAPS];
  int	 	ntaps = NTAPS;
  int		nsteps = NSTEPS;
  int		i, j;
  double 	result;
  double	step_size;
  int		c;
  int		verbose = 0;

  global_B = 0.25;

  while ((c = getopt (argc, argv, "n:t:B:v")) != EOF){
    switch (c){
    case 'n':
      nsteps = strtol (optarg, 0, 0);
      break;

    case 't':
      ntaps = strtol (optarg, 0, 0);
      break;

    case 'B':
      global_B = strtod (optarg, 0);
      break;

    case 'v':
      verbose = 1;
      break;

    default:
      usage (argv[0]);
      break;
    }
  }

  if ((nsteps & 1) != 0){
    fprintf (stderr, "%s: nsteps must be even\n", argv[0]);
    exit (1);
  }

  if (nsteps > MAX_NSTEPS){
    fprintf (stderr, "%s: nsteps must be < %d\n", argv[0], MAX_NSTEPS);
    exit (1);
  }

  if ((ntaps & 1) != 0){
    fprintf (stderr, "%s: ntaps must be even\n", argv[0]);
    exit (1);
  }

  if (nsteps > MAX_NTAPS){
    fprintf (stderr, "%s: ntaps must be < %d\n", argv[0], MAX_NTAPS);
    exit (1);
  }

  if (global_B < 0 || global_B > 0.5){
    fprintf (stderr, "%s: bandwidth must be in the range (0, 0.5)\n", argv[0]);
    exit (1);
  }

  global_gsl_int_workspace = gsl_integration_workspace_alloc(4000);
  if (global_gsl_int_workspace == NULL) {
    fprintf (stderr, "%s: unable to allocate GSL integration work space\n",
             argv[0]);
    exit (1);
  }

  step_size = 1.0/nsteps;

  // compute optimal values for mu <= 1.0

  for (j = 0; j <= nsteps; j++){

    global_mu = j * step_size;		// this determines the MU for which we're computing the taps

    // initialize X to a reasonable starting value

    initpt (&xx[j][0], ntaps);

    // find the value of X that minimizes the value of OBJECTIVE

    result = prax2_ (objective, &xx[j][0], &ntaps, &xx[j][0]);

    if (verbose){
      fprintf (stderr, "Mu: %10.8f\t", global_mu);
      fprintf (stderr, "Objective: %g\n", result);
    }
  }

  gsl_integration_workspace_free(global_gsl_int_workspace);

  // now print out the table

  printf ("\
/*\n\
 * This file was machine generated by gen_interp_differentiator_taps.\n\
 * DO NOT EDIT BY HAND.\n\
 */\n\n");

  printf ("static const int	DNTAPS     = %4d;\n", ntaps);
  printf ("static const int	DNSTEPS    = %4d;\n", nsteps);
  printf ("static const double	DBANDWIDTH = %g;\n\n", global_B);

  printf ("static const float Dtaps[DNSTEPS+1][DNTAPS] = {\n");
  printf ("  //    -4            -3            -2            -1             0             1             2             3                mu\n");

  for (i = 0; i <= nsteps; i++)
    printline (xx[i], ntaps, i, nsteps);

  printf ("};\n\n");

  return 0;
}


text:
/* vector/vector.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 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 <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>

/* turn off range checking at runtime if zero */
int gsl_check_range = 1;

#define BASE_GSL_COMPLEX_LONG
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_LONG

#define BASE_GSL_COMPLEX
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX

#define BASE_GSL_COMPLEX_FLOAT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_GSL_COMPLEX_FLOAT

#define BASE_LONG_DOUBLE
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_LONG_DOUBLE

#define BASE_DOUBLE
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_DOUBLE

#define BASE_FLOAT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_FLOAT

#define BASE_ULONG
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_ULONG

#define BASE_LONG
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_LONG

#define BASE_UINT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_UINT

#define BASE_INT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_INT

#define BASE_USHORT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_USHORT

#define BASE_SHORT
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_SHORT

#define BASE_UCHAR
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_UCHAR

#define BASE_CHAR
#include "templates_on.h"
#include "vector_source.c"
#include "templates_off.h"
#undef  BASE_CHAR


text:
/********************************************************************
 * BenchIT - Performance Measurement for Scientific Applications
 * Contact: developer@benchit.org
 *
 * $Id: kernel_main.c 1 2009-09-11 12:26:19Z william $
 * $URL: svn+ssh://william@rupert.zih.tu-dresden.de/svn-base/benchit-root/BenchITv6/kernel/numerical/BLAS3/C/0/MKL/sgemm/kernel_main.c $
 * For license details see COPYING in the package base directory
 *******************************************************************/
/* Kernel: Measurment of sgemm performance
 *******************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "interface.h"
#include "aligned_memory.h"

#ifdef USE_MKL
  #include <mkl_cblas.h>
  #define BLASOP cblas_sgemm	/* blas operation */
#endif
#ifdef USE_ACML
  #include <acml.h>
  #define BLASOP sgemm		/* blas operation */
#endif
#ifdef USE_ATLAS
  #include <cblas.h>
  #define BLASOP cblas_sgemm	/* blas operation */
#endif

#define DT float		/* datatype */
#define ALIGN 16		/* alignment */
#define EPSILON 1.0e-6		/* epsilon for result validation, 1.0e-6 for single and 1.0e-15 for double */

#define KERNELDES "sgemm"
#define LEGENDY "FLOPS (sgemm)"

#ifndef __work_h
#define __work_h

#if (defined (_CRAYMPP) || \
     defined (_USE_SHORT_AS_INT))
typedef short myinttype;
#elif (defined (_USE_LONG_AS_INT))
typedef long myinttype;
#else
typedef int myinttype;
#endif

/** The data structure that holds all the data.
 *  Please use this construct instead of global variables.
 *  Global Variables seem to have large access times (?).
 */
typedef struct mydata
{
   myinttype maxsize;
   DT *A, *B, *C;
} mydata_t;

#endif

/**  The implementation of the bi_getinfo from the BenchIT interface.
 *   Here the infostruct is filled with information about the
 *   kernel.
 *   @param infostruct  a pointer to a structure filled with zeros
 */
void bi_getinfo(bi_info * infostruct)
{
   char * p = 0;
   mydata_t * penv;

   penv = (mydata_t *) malloc(sizeof(mydata_t));

   /* get environment variables for the kernel */
   /* parameter list */
   p = bi_getenv("BENCHIT_KERNEL_PROBLEMLIST", 0);
   bi_parselist(p);

   /* get environment variables for the kernel */
   infostruct->codesequence = bi_strdup("start kernel; matrix-matrix multiplication; ");
   infostruct->kerneldescription = bi_strdup(KERNELDES);
   infostruct->xaxistext = bi_strdup("Problem Size");
   infostruct->num_measurements = infostruct->listsize;
   infostruct->num_processes = 1;
   infostruct->num_threads_per_process = 0;
   infostruct->kernel_execs_mpi1 = 0;
   infostruct->kernel_execs_mpi2 = 0;
   infostruct->kernel_execs_pvm = 0;
   infostruct->kernel_execs_omp = 0;
   infostruct->kernel_execs_pthreads = 0;
   infostruct->numfunctions = 1;

   /* allocating memory for y axis texts and properties */
   allocYAxis(infostruct);
   /* setting up y axis texts and properties */
   infostruct->yaxistexts[0] = bi_strdup("flop/s");
   infostruct->selected_result[0] = SELECT_RESULT_HIGHEST;
   infostruct->base_yaxis[0] = 0; //logarythmic axis 10^x
   infostruct->legendtexts[0] = bi_strdup(LEGENDY);

   /* free all used space */
   if (penv) free(penv);
}

/** Implementation of the bi_init of the BenchIT interface.
 *  Here you have the chance to allocate the memory you need.
 *  It is also possible to allocate the memory at the beginning
 *  of every single measurement and to free the memory thereafter.
 *  But always making use of the same memory is faster.
 *  HAVE A LOOK INTO THE HOWTO !
 */
void* bi_init(int problemSizemax)
{
  mydata_t * pmydata;
  myinttype i, m, tmp;

  pmydata = (mydata_t*)malloc(sizeof(mydata_t));
  if (pmydata == 0)
  {
     fprintf(stderr, "Allocation of structure mydata_t failed\n"); fflush(stderr);
     exit(127);
  }

  m = (myinttype)bi_get_list_element(1);
  for(i=2; i<=problemSizemax; i++){
    tmp = (myinttype)bi_get_list_element(i);
    if(tmp>m) m = tmp;
  }
  pmydata->maxsize = m;

  IDL(3, printf("\nAllocate (%d,%d)-Matrix\n",m,m));

  pmydata->A = (DT*)_aligned_malloc(m*m * sizeof(DT), ALIGN);
  pmydata->B = (DT*)_aligned_malloc(m*m * sizeof(DT), ALIGN);
  pmydata->C = (DT*)_aligned_malloc(m*m * sizeof(DT), ALIGN);

  /* matrix A is stored column by column like in fortran
   * elseif you must call cblas_sgemm with "CblasRowMajor" and thru that the algo
   * first transpose the matrix -> less performance
   */
  for(i=0; i<m*m; i++){
    pmydata->A[i] = 0.5;
    pmydata->B[i] = 2.0;
  }

  return (void *)pmydata;
}

/** initialize array with 0
 */
void initzero(DT *array, int size)
{
  myinttype i;

  for(i=0; i<size; i++){
    array[i] = 0.0;
  }
}

/** The central function within each kernel. This function
 *  is called for each measurement step seperately.
 *  @param  mdpv         a pointer to the structure created in bi_init,
 *                       it is the pointer the bi_init returns
 *  @param  problemSize  the actual problemSize
 *  @param  results      a pointer to a field of doubles, the
 *                       size of the field depends on the number
 *                       of functions, there are #functions+1
 *                       doubles
 *  @return 0 if the measurement was sucessfull, something
 *          else in the case of an error
 */
int bi_entry(void * mdpv, int iproblemSize, double * dresults)
{
  /* dstart, dend: the start and end time of the measurement */
  /* dtime: the time for a single measurement in seconds */
  double dstart = 0.0, dend = 0.0, dtime = 0.0;
  /* flops stores the calculated FLOPS */
  double dres = 0.0, dm = 0.0;
  /* ii is used for loop iterations */
  myinttype imyproblemSize;
  /* cast void* pointer */
  mydata_t * pmydata = (mydata_t *) mdpv;
#ifdef USE_ACML
  char trans = 'N';
#endif

  DT alpha=1.0;
  myinttype m, incxy=1;

  /* get current problem size from problemlist */
  imyproblemSize = (myinttype)bi_get_list_element(iproblemSize);
  m = imyproblemSize;

  /* check wether the pointer to store the results in is valid or not */
  if (dresults == NULL) return 1;

  initzero(pmydata->C, pmydata->maxsize*pmydata->maxsize);

  /* get the actual time
   * do the measurement / your algorythm
   * get the actual time
   */
  dstart = bi_gettime();
#ifdef USE_ACML
  BLASOP(trans, trans, m, m, m, alpha, pmydata->A, m, pmydata->B, m, alpha, pmydata->C, m);
#else
  BLASOP(CblasColMajor, CblasNoTrans, CblasNoTrans, m, m, m, alpha, pmydata->A, m, pmydata->B, m, alpha, pmydata->C, m);
#endif
  dend = bi_gettime();

  if(abs((pmydata->C)[0]-(DT)m) > abs((DT)m)*EPSILON || abs((pmydata->C)[m*m-1]-(DT)m) > abs((DT)m)*EPSILON) {
    IDL(0, printf("ERROR: the result is not valid!\n"));
    IDL(0, printf("expected: C(1,1)=%e, C(m,m)=%e,   got: C(1,1)=%e, C(m,m)=%e\n",(DT)m,(DT)m,(pmydata->C)[0],(pmydata->C)[m*m-1]));
  }

  /* calculate the used time and FLOPS */
  dtime = dend - dstart;
  dtime -= dTimerOverhead;
  dm = (double)m;
  dres  = dm * dm * dm + (dm-1) * dm * dm;	/* m*m*m mult and m-1*m*m add for A*B */
  dres += dm * dm;				/* m*m mult for alpha */
  dres += 2.0 * dm * dm;			/* m*m mult for C and m*m add for final result */

  /* If the operation was too fast to be measured by the timer function,
   * mark the result as invalid
   */
  if(dtime < dTimerGranularity) dtime = INVALID_MEASUREMENT;

  /* store the results in results[1], results[2], ...
   * [1] for the first function, [2] for the second function
   * and so on ...
   * the index 0 always keeps the value for the x axis
   */
  dresults[0] = (double)imyproblemSize;
  dresults[1] = (dtime!=INVALID_MEASUREMENT) ? dres / dtime : INVALID_MEASUREMENT;

  return 0;
}

/** Clean up the memory
 */
void bi_cleanup(void* mdpv)
{
   mydata_t * pmydata = (mydata_t*)mdpv;

   _aligned_free(pmydata->A);
   _aligned_free(pmydata->B);
   _aligned_free(pmydata->C);

   if (pmydata) free(pmydata);
   return;
}


text:
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_math.h>

#include "allvars.h"
#include "proto.h"

#ifdef CR_DIFFUSION_GREEN
#include "cosmic_rays.h"

void compute_diff_weights(int target, int mode);
void scatter_diffusion(int target, int mode);

#define m_p (PROTONMASS * All.HubbleParam / All.UnitMass_in_g)

void greenf_diffusion(void)
{
  int *noffset, *nbuffer, *nsend, *nsend_local, *numlist, *ndonelist;
  int i, j, n;
  int ndone;
  long long ntot, ntotleft;
  int maxfill, source;
  int level, ngrp, sendTask, recvTask;
  int place, nexport;
  double cr_efac_i, kappa_egy;
  double kappa, a3inv, egysum, egytot, egytot_before;
  double meanKineticEnergy, qmeanKin;
  double CR_q_i;
  MPI_Status status;

  if(ThisTask == 0)
    {
      printf("Doing diffusion step with Green function method\n");
      fflush(stdout);
    }

  for(CRpop = 0; CRpop < NUMCRPOP; CRpop++)
    {

      noffset = mymalloc(sizeof(int) * NTask);	/* offsets of bunches in common list */
      nbuffer = mymalloc(sizeof(int) * NTask);
      nsend_local = mymalloc(sizeof(int) * NTask);
      nsend = mymalloc(sizeof(int) * NTask * NTask);
      ndonelist = mymalloc(sizeof(int) * NTask);

      if(All.ComovingIntegrationOn)
	a3inv = 1 / (All.Time * All.Time * All.Time);
      else
	a3inv = 1;

      egysum = 0;

      for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
	{
	  if(P[n].Type == 0)
	    {
	      NumSphUpdate++;

	      egysum += SphP[n].CR_E0[CRpop];

	      SphP[n].CR_DeltaE[CRpop] = 0;
	      SphP[n].CR_DeltaN[CRpop] = 0;

	      kappa = All.CR_DiffusionCoeff;

	      if(All.CR_DiffusionDensScaling != 0.0)
		kappa *= pow(SphP[n].d.Density * a3inv / All.CR_DiffusionDensZero, All.CR_DiffusionDensScaling);

	      if(All.CR_DiffusionEntropyScaling != 0.0)
		kappa *= pow(SphP[n].Entropy / All.CR_DiffusionEntropyZero, All.CR_DiffusionEntropyScaling);

	      if(SphP[n].CR_E0[CRpop] > 0)
		{
		  CR_q_i = SphP[n].CR_q0[CRpop] * pow(SphP[n].d.Density * a3inv, 0.33333);

		  cr_efac_i =
		    CR_Tab_MeanEnergy(CR_q_i, All.CR_Alpha[CRpop] - 0.3333, CRpop)
		    / CR_Tab_MeanEnergy(CR_q_i, All.CR_Alpha[CRpop], CRpop);

		  kappa *= (All.CR_Alpha[CRpop] - 1) / (All.CR_Alpha[CRpop] - 1.33333) * pow(CR_q_i, 0.3333);

		  kappa_egy = kappa * cr_efac_i;
		}
	      else
		kappa_egy = kappa;

	      SphP[n].CR_Kappa[CRpop] = kappa;
	      SphP[n].CR_Kappa_egy[CRpop] = kappa_egy;
	    }
	}

      MPI_Allreduce(&egysum, &egytot_before, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

      numlist = mymalloc(NTask * sizeof(int) * NTask);
      MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
      for(i = 0, ntot = 0; i < NTask; i++)
	ntot += numlist[i];
      myfree(numlist);

      /* first, sum the weights */

      i = 0;			/* beginn with this index */
      ntotleft = ntot;		/* particles left for all tasks together */

      while(ntotleft > 0)
	{
	  for(j = 0; j < NTask; j++)
	    nsend_local[j] = 0;

	  /* do local particles and prepare export list */

	  for(nexport = 0, ndone = 0; i < NumPart && nexport < All.BunchSizeDensity - NTask; i++)
	    if(P[i].Type == 0)
	      {
		ndone++;

		for(j = 0; j < NTask; j++)
		  Exportflag[j] = 0;

		compute_diff_weights(i, 0);

		for(j = 0; j < NTask; j++)
		  {
		    if(Exportflag[j])
		      {
			DensDataIn[nexport].Pos[0] = P[i].Pos[0];
			DensDataIn[nexport].Pos[1] = P[i].Pos[1];
			DensDataIn[nexport].Pos[2] = P[i].Pos[2];
			DensDataIn[nexport].Hsml = PPP[i].Hsml;
			DensDataIn[nexport].CR_Kappa[CRpop] = SphP[i].CR_Kappa[CRpop];
			DensDataIn[nexport].CR_Kappa_egy[CRpop] = SphP[i].CR_Kappa_egy[CRpop];

			DensDataIn[nexport].Index = i;
			DensDataIn[nexport].Task = j;
			nexport++;
			nsend_local[j]++;
		      }
		  }
	      }

	  qsort(DensDataIn, nexport, sizeof(struct densdata_in), dens_compare_key);

	  for(j = 1, noffset[0] = 0; j < NTask; j++)
	    noffset[j] = noffset[j - 1] + nsend_local[j - 1];

	  MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);

	  /* now do the particles that need to be exported */

	  for(level = 1; level < (1 << PTask); level++)
	    {
	      for(j = 0; j < NTask; j++)
		nbuffer[j] = 0;
	      for(ngrp = level; ngrp < (1 << PTask); ngrp++)
		{
		  maxfill = 0;
		  for(j = 0; j < NTask; j++)
		    {
		      if((j ^ ngrp) < NTask)
			if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
			  maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		    }
		  if(maxfill >= All.BunchSizeDensity)
		    break;

		  sendTask = ThisTask;
		  recvTask = ThisTask ^ ngrp;

		  if(recvTask < NTask)
		    {
		      if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
			{
			  /* get the particles */
			  MPI_Sendrecv(&DensDataIn[noffset[recvTask]],
				       nsend_local[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				       recvTask, TAG_CONDUCT_A,
				       &DensDataGet[nbuffer[ThisTask]],
				       nsend[recvTask * NTask + ThisTask] * sizeof(struct densdata_in),
				       MPI_BYTE, recvTask, TAG_CONDUCT_A, MPI_COMM_WORLD, &status);
			}
		    }

		  for(j = 0; j < NTask; j++)
		    if((j ^ ngrp) < NTask)
		      nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
		}

	      for(j = 0; j < nbuffer[ThisTask]; j++)
		compute_diff_weights(j, 1);

	      for(j = 0; j < NTask; j++)
		nbuffer[j] = 0;
	      for(ngrp = level; ngrp < (1 << PTask); ngrp++)
		{
		  maxfill = 0;
		  for(j = 0; j < NTask; j++)
		    {
		      if((j ^ ngrp) < NTask)
			if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
			  maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		    }
		  if(maxfill >= All.BunchSizeDensity)
		    break;

		  sendTask = ThisTask;
		  recvTask = ThisTask ^ ngrp;

		  if(recvTask < NTask)
		    {
		      if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
			{
			  /* send the results */
			  MPI_Sendrecv(&DensDataResult[nbuffer[ThisTask]],
				       nsend[recvTask * NTask + ThisTask] * sizeof(struct densdata_out),
				       MPI_BYTE, recvTask, TAG_CONDUCT_B,
				       &DensDataPartialResult[noffset[recvTask]],
				       nsend_local[recvTask] * sizeof(struct densdata_out),
				       MPI_BYTE, recvTask, TAG_CONDUCT_B, MPI_COMM_WORLD, &status);

			  /* add the result to the particles */
			  for(j = 0; j < nsend_local[recvTask]; j++)
			    {
			      source = j + noffset[recvTask];
			      place = DensDataIn[source].Index;

			      SphP[place].CR_WeightSum += DensDataPartialResult[source].CR_WeightSum;
			      SphP[place].CR_WeightSum_egy += DensDataPartialResult[source].CR_WeightSum_egy;
			    }
			}
		    }

		  for(j = 0; j < NTask; j++)
		    if((j ^ ngrp) < NTask)
		      nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
		}
	      level = ngrp - 1;
	    }

	  MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
	  for(j = 0; j < NTask; j++)
	    ntotleft -= ndonelist[j];
	}

  /*************** now do the diffusion step itself */

      numlist = mymalloc(NTask * sizeof(int) * NTask);
      MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
      for(i = 0, ntot = 0; i < NTask; i++)
	ntot += numlist[i];
      myfree(numlist);

      /* first, sum the weights */

      i = 0;			/* beginn with this index */
      ntotleft = ntot;		/* particles left for all tasks together */

      while(ntotleft > 0)
	{
	  for(j = 0; j < NTask; j++)
	    nsend_local[j] = 0;

	  /* do local particles and prepare export list */

	  for(nexport = 0, ndone = 0; i < NumPart && nexport < All.BunchSizeDensity - NTask; i++)
	    if(P[i].Type == 0)
	      {
		ndone++;

		for(j = 0; j < NTask; j++)
		  Exportflag[j] = 0;

		scatter_diffusion(i, 0);

		for(j = 0; j < NTask; j++)
		  {
		    if(Exportflag[j])
		      {
			DensDataIn[nexport].Pos[0] = P[i].Pos[0];
			DensDataIn[nexport].Pos[1] = P[i].Pos[1];
			DensDataIn[nexport].Pos[2] = P[i].Pos[2];

			DensDataIn[nexport].CR_Kappa[CRpop] = SphP[i].CR_Kappa[CRpop];
			DensDataIn[nexport].CR_Kappa_egy[CRpop] = SphP[i].CR_Kappa_egy[CRpop];
			DensDataIn[nexport].CR_WeightSum = SphP[i].CR_WeightSum;
			DensDataIn[nexport].CR_WeightSum_egy = SphP[i].CR_WeightSum_egy;
			DensDataIn[nexport].Hsml = PPP[i].Hsml;
			DensDataIn[nexport].CR_E0[CRpop] = SphP[i].CR_E0[CRpop];
			DensDataIn[nexport].CR_n0[CRpop] = SphP[i].CR_n0[CRpop];

			DensDataIn[nexport].Index = i;
			DensDataIn[nexport].Task = j;
			nexport++;
			nsend_local[j]++;
		      }
		  }
	      }

	  qsort(DensDataIn, nexport, sizeof(struct densdata_in), dens_compare_key);

	  for(j = 1, noffset[0] = 0; j < NTask; j++)
	    noffset[j] = noffset[j - 1] + nsend_local[j - 1];

	  MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);

	  /* now do the particles that need to be exported */

	  for(level = 1; level < (1 << PTask); level++)
	    {
	      for(j = 0; j < NTask; j++)
		nbuffer[j] = 0;
	      for(ngrp = level; ngrp < (1 << PTask); ngrp++)
		{
		  maxfill = 0;
		  for(j = 0; j < NTask; j++)
		    {
		      if((j ^ ngrp) < NTask)
			if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
			  maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
		    }
		  if(maxfill >= All.BunchSizeDensity)
		    break;

		  sendTask = ThisTask;
		  recvTask = ThisTask ^ ngrp;

		  if(recvTask < NTask)
		    {
		      if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
			{
			  /* get the particles */
			  MPI_Sendrecv(&DensDataIn[noffset[recvTask]],
				       nsend_local[recvTask] * sizeof(struct densdata_in), MPI_BYTE,
				       recvTask, TAG_CONDUCT_A,
				       &DensDataGet[nbuffer[ThisTask]],
				       nsend[recvTask * NTask + ThisTask] * sizeof(struct densdata_in),
				       MPI_BYTE, recvTask, TAG_CONDUCT_A, MPI_COMM_WORLD, &status);
			}
		    }

		  for(j = 0; j < NTask; j++)
		    if((j ^ ngrp) < NTask)
		      nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
		}

	      for(j = 0; j < nbuffer[ThisTask]; j++)
		scatter_diffusion(j, 1);

	      level = ngrp - 1;
	    }

	  MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
	  for(j = 0; j < NTask; j++)
	    ntotleft -= ndonelist[j];
	}

      /* now set the new cosmic ray prorperties */

      egysum = 0;

      for(n = 0; n < N_gas; n++)
	{
	  if(P[n].Type == 0)
	    {
	      SphP[n].CR_E0[CRpop] = SphP[n].CR_DeltaE[CRpop];
	      SphP[n].CR_n0[CRpop] = SphP[n].CR_DeltaN[CRpop];

	      egysum += SphP[n].CR_E0[CRpop];

	      SphP[n].CR_DeltaE[CRpop] = 0;
	      SphP[n].CR_DeltaN[CRpop] = 0;

	      if(SphP[n].CR_n0[CRpop] > 1.0e-12 && SphP[n].CR_E0[CRpop] > 0)
		{
		  meanKineticEnergy = SphP[n].CR_E0[CRpop] * m_p / SphP[n].CR_n0[CRpop];

		  qmeanKin = CR_q_from_mean_kinetic_energy(meanKineticEnergy[CRpop], CRpop);

		  SphP[n].CR_q0[CRpop] = qmeanKin * pow(SphP[n].d.Density * a3inv, -(1.0 / 3.0));
		  SphP[n].CR_C0[CRpop] = SphP[n].CR_n0[CRpop] * (All.CR_Alpha[CRpop] - 1.0) *
		    pow(SphP[n].CR_q0[CRpop], All.CR_Alpha[CRpop] - 1.0);
		}
	      else
		{
		  SphP[n].CR_E0[CRpop] = 0.0;
		  SphP[n].CR_n0[CRpop] = 0.0;

		  SphP[n].CR_q0[CRpop] = 1.0e10;
		  SphP[n].CR_C0[CRpop] = 0.0;
		}
	    }
	}

      MPI_Allreduce(&egysum, &egytot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
      if(ThisTask == 0)
	{
	  printf("Energy before/after= %g | %g\n", egytot_before, egytot);
	  fflush(stdout);
	}

      myfree(ndonelist);
      myfree(nsend);
      myfree(nsend_local);
      myfree(nbuffer);
      myfree(noffset);

      All.TimeOfLastDiffusion = All.Time;
    }
}

void compute_diff_weights(int target, int mode)
{
  int j, n;
  int startnode, numngb_inbox;
  double h;
  double weightsum, weightsum_egy, kappa, kappa_egy, kappaeff, kappaeff_egy;;
  double dx, dy, dz, r, r2;
  MyFloat *pos;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      kappa = SphP[target].CR_Kappa;
      kappa_egy = SphP[target].CR_Kappa_egy;
    }
  else
    {
      pos = DensDataGet[target].Pos;
      h = DensDataGet[target].Hsml;
      kappa = DensDataGet[target].CR_Kappa;
      kappa_egy = DensDataGet[target].CR_Kappa_egy;
    }

  weightsum = weightsum_egy = 0;

  startnode = All.MaxPart;

  kappaeff = 4 * kappa * (All.Time - All.TimeOfLastDiffusion);
  kappaeff_egy = 4 * kappa_egy * (All.Time - All.TimeOfLastDiffusion);

  if(kappaeff > 0)
    {
      do
	{
	  numngb_inbox = ngb_treefind_variable(&pos[0], h, &startnode);

	  for(n = 0; n < numngb_inbox; n++)
	    {
	      j = Ngblist[n];

	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;

	      if(r2 < h * h)
		{
		  r = sqrt(r2);

		  weightsum += P[j].Mass * exp(-r2 / kappaeff);
		  weightsum_egy += P[j].Mass * exp(-r2 / kappaeff_egy);
		}
	    }
	}
      while(startnode >= 0);
    }

  if(mode == 0)
    {
      SphP[target].CR_WeightSum = weightsum;
      SphP[target].CR_WeightSum_egy = weightsum_egy;
    }
  else
    {
      DensDataResult[target].CR_WeightSum = weightsum;
      DensDataResult[target].CR_WeightSum_egy = weightsum_egy;
    }

}

void scatter_diffusion(int target, int mode)
{
  int j, n;
  int startnode, numngb_inbox;
  double h, weight;
  double weightsum, kappa, kappaeff, CR_E0, CR_n0;
  double weightsum_egy, kappa_egy, kappaeff_egy;
  double dx, dy, dz, r, r2;
  MyFloat *pos;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = PPP[target].Hsml;
      kappa = SphP[target].CR_Kappa;
      kappa_egy = SphP[target].CR_Kappa_egy;
      weightsum = SphP[target].CR_WeightSum;
      weightsum_egy = SphP[target].CR_WeightSum_egy;
      CR_E0 = SphP[target].CR_E0;
      CR_n0 = SphP[target].CR_n0;
    }
  else
    {
      pos = DensDataGet[target].Pos;
      h = DensDataGet[target].Hsml;
      kappa = DensDataGet[target].CR_Kappa;
      kappa_egy = DensDataGet[target].CR_Kappa_egy;
      weightsum = DensDataGet[target].CR_WeightSum;
      weightsum_egy = DensDataGet[target].CR_WeightSum_egy;
      CR_E0 = DensDataGet[target].CR_E0;
      CR_n0 = DensDataGet[target].CR_n0;
    }

  startnode = All.MaxPart;

  kappaeff = 4 * kappa * (All.Time - All.TimeOfLastDiffusion);
  kappaeff_egy = 4 * kappa_egy * (All.Time - All.TimeOfLastDiffusion);

  if(kappaeff > 0)
    {
      do
	{
	  numngb_inbox = ngb_treefind_variable(&pos[0], h, &startnode);

	  for(n = 0; n < numngb_inbox; n++)
	    {
	      j = Ngblist[n];

	      dx = pos[0] - P[j].Pos[0];
	      dy = pos[1] - P[j].Pos[1];
	      dz = pos[2] - P[j].Pos[2];

#ifdef PERIODIC			/*  now find the closest image in the given box size  */
	      if(dx > boxHalf_X)
		dx -= boxSize_X;
	      if(dx < -boxHalf_X)
		dx += boxSize_X;
	      if(dy > boxHalf_Y)
		dy -= boxSize_Y;
	      if(dy < -boxHalf_Y)
		dy += boxSize_Y;
	      if(dz > boxHalf_Z)
		dz -= boxSize_Z;
	      if(dz < -boxHalf_Z)
		dz += boxSize_Z;
#endif
	      r2 = dx * dx + dy * dy + dz * dz;

	      if(r2 < h * h)
		{
		  r = sqrt(r2);

		  weight = P[j].Mass * exp(-r2 / kappaeff) / weightsum;
		  SphP[j].CR_DeltaN += CR_n0 * weight;

		  weight = P[j].Mass * exp(-r2 / kappaeff_egy) / weightsum_egy;
		  SphP[j].CR_DeltaE += CR_E0 * weight;

		}
	    }
	}
      while(startnode >= 0);
    }
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/**
 *
 * @file pzgbrdb.c
 *
 *  PLASMA auxiliary 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 Azzam Haidar
 * @date 2011-05-15
 * @precisions normal z -> c d s
 *
 **/
#include "common.h"
#include <lapacke.h>

#undef REAL
#define COMPLEX

#define DEP(m)  &(DEP[m])
#define A(_m, _n) (PLASMA_Complex64_t *)plasma_geteltaddr(&A, (_m), (_n), eltsize)
/***************************************************************************//**
 *  Parallel Reduction from BAND Bidiagonal to the final condensed form - dynamic scheduler
 **/
void plasma_pzgbrdb_quark(PLASMA_enum uplo,
                          PLASMA_desc A, double *D, double *E, PLASMA_desc T,
                          PLASMA_sequence *sequence, PLASMA_request *request)
{
    plasma_context_t *plasma;
    Quark_Task_Flags task_flags = Quark_Task_Flags_Initializer;

#ifdef COMPLEX
    static double dzero = (double) 0.0;
    double absztmp;
#endif

    static PLASMA_Complex64_t zone  = (PLASMA_Complex64_t) 1.0;
    static PLASMA_Complex64_t zzero = (PLASMA_Complex64_t) 0.0;
    PLASMA_Complex64_t *VQ, *TAUQ,*VP, *TAUP;
    PLASMA_Complex64_t ztmp, V, TAU;
    int M, N, NB, MINMN, INgrsiz, INthgrsiz, BAND;
    int myid, grsiz, shift=3, stt, st, ed, stind, edind;
    int blklastind, colpt, PCOL, ACOL, MCOL;
    int stepercol, mylastid, grnb, grid;
    int *DEP,*MAXID;
    int i, sweepid, m;
    int thgrsiz, thgrnb, thgrid, thed;
    size_t eltsize = plasma_element_size(A.dtyp);

    plasma = plasma_context_self();
    if (sequence->status != PLASMA_SUCCESS)
        return;

    QUARK_Task_Flag_Set(&task_flags, TASK_SEQUENCE, (intptr_t)sequence->quark_sequence);

    M     = A.m;
    N     = A.n;
    NB    = A.mb;
    MINMN = min(M,N);

    /* Quick return */
    if ( MINMN == 0 ){
        return;
    }
    if ( NB == 0 ) {
        memset(D, 0,  MINMN   *sizeof(double));
        memset(E, 0, (MINMN-1)*sizeof(double));
#ifdef COMPLEX
        for (i=0; i<MINMN; i++)
            D[i] = cabs(*A(i,i));
#else
        for (i=0; i<MINMN; i++)
            D[i] = *A(i,i);
#endif
        return;
    }

    /*
     * Barrier is used because the bulge have to wait until
     * the reduction to band has been finish.
     * otherwise, I can remove this BARRIER when I integrate
     * the function dependencies link inside the reduction to
     * band. Keep in mind the case when NB=1, where no bulge-chasing.
     */
    /***************************************************************/
    QUARK_Barrier(plasma->quark);
    /***************************************************************/

    /*
     * Case NB=1 ==> matrix is already Bidiagonal. no need to bulge.
     * Make diagonal and superdiagonal elements real, storing them in
     * D and E. if PlasmaLower, first transform lower bidiagonal form
     * to upper bidiagonal by applying plane rotations/ Householder
     * from the left, overwriting superdiagonal elements then make
     * elements real of the resulting upper Bidiagonal. if PlasmaUpper
     * then make its elements real.  For Q, PT: ZSCAL should be done
     * in case of WANTQ.
     */
    if ( NB == 1 ) {
        memset(D, 0, MINMN*sizeof(double));
        memset(E, 0, (MINMN-1)*sizeof(double));

        if(uplo==PlasmaLower){
            for (i=0; i<(MINMN-1); i++)
            {
                /* generate Householder to annihilate a(i+1,i) and create a(i,i+1) */
                V             = *A((i+1), i);
                *A((i+1),  i) = zzero;
                LAPACKE_zlarfg_work( 2, A(i, i), &V, 1, &TAU);
                /* apply Left*/
                TAU  = conj(TAU);
                ztmp = TAU*V;
                V    = conj(V);
                *A(i,   i+1) = - V * TAU * (*A(i+1, i+1));
                *A(i+1, i+1) = *(A(i+1, i+1)) * (zone - V * ztmp);
            }
        }
        /* PlasmaLower or PlasmaUpper, both are now upper */
        /* Make diagonal and superdiagonal elements real,
         * storing them in D and E
         */
#ifdef COMPLEX
        ztmp = zone;
        for (i=0; i<MINMN; i++)
        {
            ztmp     = *A(i, i) * conj(ztmp);
            absztmp  = cabs(ztmp);
            D[i]     = absztmp;               /* diag value */
            if(absztmp != dzero)
                ztmp = (PLASMA_Complex64_t) (ztmp / absztmp);
            else
                ztmp = zone;
            if(i<(MINMN-1)) {
                ztmp     = *A(i, (i+1)) * conj(ztmp);
                absztmp  = cabs(ztmp);
                E[i]     = absztmp;            /* upper off-diag value */
                if(absztmp != dzero)
                    ztmp = (PLASMA_Complex64_t) (ztmp / absztmp);
                else
                    ztmp = zone;
            }
        }
#else
        for (i=0; i < MINMN-1; i++) {
            D[i] = *A(i, i  );
            E[i] = *A(i, i+1);
        }
        D[i] = *A(i, i);
#endif
        return;
    }

    /*
     * Case MINMN<NB ==> matrix is very small and better to call lapack ZGETRD.
     *
     *    Use fact that one row of block is stored the same way than in LAPACK
     *    Doesn't work if M > NB because of tile storage
     */
    if ( MINMN <= 0 )
    {
        PLASMA_Complex64_t *work, *taup, *tauq;
        int info, ldwork = N*N;
        work = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, ldwork, PlasmaComplexDouble);
        taup = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, MINMN,  PlasmaComplexDouble);
        tauq = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, MINMN,  PlasmaComplexDouble);

        info = LAPACKE_zgebrd_work(LAPACK_COL_MAJOR, M, N,
                                   A(0,0), A.lm, D, E, taup, tauq, work, ldwork);
        plasma_shared_free(plasma, (void*) work);
        plasma_shared_free(plasma, (void*) taup);
        plasma_shared_free(plasma, (void*) tauq);

        if( info == 0 )
            sequence->status = PLASMA_SUCCESS;
        else
            plasma_sequence_flush(plasma->quark, sequence, request, info);
        return;
    }

    /* General case NB > 1 && N > NB */
    DEP   = (int *)                plasma_shared_alloc(plasma, MINMN+1, PlasmaInteger      );
    MAXID = (int *)                plasma_shared_alloc(plasma, MINMN+1, PlasmaInteger      );
    memset(MAXID,0,(MINMN+1)*sizeof(int));
    VQ    = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN,   PlasmaComplexDouble);
    TAUQ  = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN,   PlasmaComplexDouble);
    VP    = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN,   PlasmaComplexDouble);
    TAUP  = (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, 2*MINMN,   PlasmaComplexDouble);
    memset(VQ,0,2*MINMN*sizeof(PLASMA_Complex64_t));
    memset(TAUQ,0,2*MINMN*sizeof(PLASMA_Complex64_t));
    memset(VP,0,2*MINMN*sizeof(PLASMA_Complex64_t));
    memset(TAUP,0,2*MINMN*sizeof(PLASMA_Complex64_t));

    /***************************************************************************
     *                       START BULGE CHASING CODE
     **************************************************************************/
    printf("converting tile to lap A.lm %d   ln %d \n",A.lm,A.ln);
    PLASMA_Complex64_t *AB =   (PLASMA_Complex64_t *) plasma_shared_alloc(plasma, A.lm*A.ln,   PlasmaComplexDouble);
    memset(AB,0,A.lm*A.ln*sizeof(PLASMA_Complex64_t));
    plasma_zooptile2lap( A,   AB, NB, NB,  A.lm, A.ln,  sequence, request);
    int LDAB = A.lm;
/*
    plasma_dynamic_call_6( plasma_pzhbcpy_t2bl,
        PLASMA_enum, uplo,
        PLASMA_desc, A,
        PLASMA_Complex64_t*, &(AB[NB]),
        int, M,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);
*/

    QUARK_Barrier(plasma->quark);

    /*
     * Initialisation of local parameter. those parameter should be
     * input or tuned parameter.
     */
    INgrsiz = 1;
    if( NB > 160 ) {
        INgrsiz = 2;
    }
    else if( NB > 100 ) {
        if( MINMN < 5000 )
            INgrsiz = 2;
        else
            INgrsiz = 4;
    } else {
        INgrsiz = 6;
    }
    INthgrsiz = MINMN;
    BAND      = 0;

    grsiz   = INgrsiz;
    thgrsiz = INthgrsiz;
    if( grsiz   == 0 ) grsiz   = 6;
    if( thgrsiz == 0 ) thgrsiz = MINMN;

    i = shift/grsiz;
    stepercol =  i*grsiz == shift ? i:i+1;

    i       = (MINMN-1)/thgrsiz;
    thgrnb  = i*thgrsiz == (MINMN-1) ? i:i+1;

    printf("starting code grsiz %d thgr %d\n", grsiz, thgrsiz);
    for (thgrid = 1; thgrid<=thgrnb; thgrid++){
        stt  = (thgrid-1)*thgrsiz+1;
        thed = min( (stt + thgrsiz -1), (MINMN-1));
        for (i = stt; i <= MINMN-1; i++){
            ed=min(i,thed);
            if(stt>ed)break;
            for (m = 1; m <=stepercol; m++){
                st=stt;
                for (sweepid = st; sweepid <=ed; sweepid++){
                    /* PCOL:  dependency on the ID of the master of the group of the previous column.  (Previous Column:PCOL). */
                    /* ACOL:  dependency on the ID of the master of the previous group of my column.   (Acctual  Column:ACOL). (it is 0(NULL) for myid=1) */
                    /* MCOL:  OUTPUT dependency on the my ID, to be used by the next ID. (My Column: MCOL). I am the master of this group. */
                    myid     = (i-sweepid)*(stepercol*grsiz) +(m-1)*grsiz + 1;
                    mylastid = myid+grsiz-1;
                    PCOL     = mylastid+shift-1;  /* to know the dependent ID of the previous column. need to know the master of its group*/
                    MAXID[sweepid] = myid;
                    PCOL     = min(PCOL,MAXID[sweepid-1]); /* for the last columns, we might do only 1 or 2 kernel, so the PCOL will be wrong. this is to force it to the last ID of the previous col.*/
                    grnb     = PCOL/grsiz;
                    grid     = grnb*grsiz == PCOL ? grnb:grnb+1;
                    PCOL     = (grid-1)*grsiz +1; /* give me the ID of the master of the group of the previous column.*/
                    ACOL     = myid-grsiz;
                    if(myid==1)ACOL=0;
                    MCOL     = myid;

                    QUARK_CORE_zbrdalg2(
                        plasma->quark, &task_flags,
                        uplo, MINMN, NB,
                        AB, LDAB, VQ, TAUQ, VP, TAUP, i, sweepid, m, grsiz, BAND,
                        DEP(PCOL), DEP(ACOL), DEP(MCOL) );

                    if(mylastid%2 ==0){
                        blklastind      = (mylastid/2)*NB+1+sweepid-1;
                    }else{
                        colpt      = ((mylastid+1)/2)*NB + 1 +sweepid -1 ;
                        stind      = colpt-NB+1;
                        edind      = min(colpt,MINMN);
                        if( (stind>=edind-1) && (edind==MINMN) )
                            blklastind=MINMN;
                        else
                            blklastind=0;
                    }
                    if(blklastind >= (MINMN-1))  stt=stt+1;
                } /* END for sweepid=st:ed    */
            } /* END for m=1:stepercol */
        } /* END for i=1:MINMN-2      */
    } /* END for thgrid=1:thgrnb     */

     /*
     * Barrier used only for now, to be sure that everything
     * is done before copying the D and E and free workspace.
     * this will be removed later when D and E are directly filled
     * during the bulge process.
     */
    QUARK_Barrier(plasma->quark);
    printf("converting back to tile\n");

    plasma_zooplap2tile( A, AB, NB, NB, A.lm, A.ln, 0, 0, M, N, sequence, request, plasma_desc_mat_free(&A));
/*
    A = plasma_desc_init(PlasmaComplexDouble, NB,NB, NB*NB, M, N, 0, 0, M, N);
        if ( plasma_desc_mat_alloc( &(A) ) ) {
        plasma_error( __func__, "plasma_shared_alloc() failed");
        {free;};
        return PLASMA_ERR_OUT_OF_RESOURCES;
        }
        */
    /*
    plasma_parallel_call_5( plasma_pzlapack_to_tile,
        PLASMA_Complex64_t*, AB,
        int,                 MINMN,
        PLASMA_desc,         A,
        PLASMA_sequence*,    sequence,
        PLASMA_request*,     request);

*/

    QUARK_Barrier(plasma->quark);
    printf("converting done\n");

    plasma_shared_free(plasma, (void*) DEP);
    plasma_shared_free(plasma, (void*) MAXID);
    plasma_shared_free(plasma, (void*) VQ);
    plasma_shared_free(plasma, (void*) TAUQ);
    plasma_shared_free(plasma, (void*) VP);
    plasma_shared_free(plasma, (void*) TAUP);

    /*
     * STORE THE RESULTING diagonal/off-diagonal in D AND E
     */
    memset(D, 0, MINMN*sizeof(double));
    memset(E, 0, (MINMN-1)*sizeof(double));

    /* Make diagonal and superdiagonal elements real,
     * storing them in D and E
     */
    if(uplo==PlasmaUpper){
        for (i=0; i < MINMN-1; i++) {
            D[i] = *A(i, i  );
            E[i] = *A(i, i+1);
        }
        D[i] = *A(i, i);
    }
    else{
        for (i=0; i < MINMN-1; i++) {
            D[i] = *A(i, i  );
            E[i] = *A(i+1, i);
        }
        D[i] = *A(i, i);
    }

} /* END FUNCTION */


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C header defining useful physical constants (values taken from LAL).
 * Also defines boolean conventions.
 *
 */

#ifndef _CONSTANTS_H
#define _CONSTANTS_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>

#if defined(__cplusplus)
extern "C" {
#elif 0
} /* so that editors will match preceding brace */
#endif

/***************************************************/
/****** Boolean conventions for loading files ******/

#define SUCCESS 0
#define FAILURE 1

/*********************************************************/
/*************** Mathematical constants ******************/

#define PI 3.1415926535897932384626433832795029
#define GAMMA 0.577215664901532860606512090082402431
#define sqrt2 1.4142135623730950488
#define sqrt3 1.7320508075688772935
#define sqrt6 2.4494897427831780982
#define invsqrt2 0.70710678118654752440
#define invsqrt3 0.57735026918962576451
#define invsqrt6 0.40824829046386301637

/**********************************************************************/
/**************** Physical constants in SI units **********************/

#define C_SI 299792458.
#define G_SI 6.67259e-11
#define MSUN_SI 1.98892e30
#define MTSUN_SI 4.9254923218988636432342917247829673e-6
#define PC_SI 3.0856775807e16
#define AU_SI 1.4959787066e11
#define YRSID_SI 3.15581497635e7 /* Sideral year as found on http://hpiers.obspm.fr/eop-pc/models/constants.html */

//#define Omega_SI 1.99098659277e-6 /* Orbital pulsation: 20pi/year for "fast-orbit-LISA2017 "*/
//#define Omega_SI 6e-7 /* Orbital pulsation: ~6pi/year for "fast-orbit-LISA2017 "*/
//#define Omega_SI 1.99098659277e-7 /* Orbital pulsation: 2pi/year - use sidereal year as found on http://hpiers.obspm.fr/eop-pc/models/constants.html */
#define EarthOrbitOmega_SI 1.99098659277e-7 /* Orbital pulsation: 2pi/year - use sidereal year as found on http://hpiers.obspm.fr/eop-pc/models/constants.html */
//#define Omega_SI 1.99098659277e-9 /* Orbital pulsation: 2pi/century - alternative model to test the impact of a non-moving LISA*/
//#define Omega_SI 1.99098659277e-11 /* Orbital pulsation: 2pi/century - alternative model to test the impact of a non-moving LISA*/
//#define f0_SI 3.168753575e-8 /* Orbital frequency: 1/year */   //not used?
//#define L_SI 5.0e9 /* Arm length of the detector (in m): (Standard LISA) */
//#define L_SI 2.5e9 /* Arm length of the detector (in m): (L3 reference LISA) */
//#define R_SI 1.4959787066e11 /* Radius of the orbit around the sun: 1AU */
#define AU_SI 1.4959787066e11 /* Radius of the orbit around the sun: 1AU */
//#define R_SI 1.4959787066e8 /* Radius of the orbit around the sun: 0.001 AU for hacked LISA test*/

/**********************************************************/
/********** Constants used to relate time scales **********/

#define EPOCH_J2000_0_TAI_UTC 32           /* Leap seconds (TAI-UTC) on the J2000.0 epoch (2000 JAN 1 12h UTC) */
#define EPOCH_J2000_0_GPS 630763213        /* GPS seconds of the J2000.0 epoch (2000 JAN 1 12h UTC) */

/*******************************************/
/**************** NaN **********************/

#ifndef NAN
# define NAN (INFINITY-INFINITY)
#endif

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _CONSTANTS_H */

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


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

#pragma once

#include <codecvt>
#include <string>
#include <locale>

#include <gsl/gsl>

namespace mira
{
    inline std::string utf16_to_utf8(const wchar_t* input)
    {
        std::wstring_convert<std::codecvt_utf8_utf16<wchar_t>> converter;
        return converter.to_bytes(input);
    }

    inline std::wstring utf8_to_utf16(const char* input)
    {
        std::wstring_convert<std::codecvt_utf8<wchar_t>> converter;
        return converter.from_bytes(input);
    }

    inline std::string utf16_to_utf8(const std::wstring& input)
    {
        return utf16_to_utf8(input.data());
    }

    inline std::wstring utf8_to_utf16(const std::string& input)
    {
        return utf8_to_utf16(input.data());
    }

    struct string_compare
    {
        using is_transparent = std::true_type;

        bool operator()(gsl::cstring_span<> a, gsl::cstring_span<> b) const
        {
            return std::lexicographical_compare(a.begin(), a.end(), b.begin(), b.end());
        }

        bool operator()(gsl::czstring_span<> a, gsl::czstring_span<> b) const
        {
            return (*this)(a.as_string_span(), b.as_string_span());
        }

        bool operator()(const char* a, const char* b) const
        {
            return strcmp(a, b) < 0;
        }
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#ifndef CG_H_
#define CG_H_

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "util.h"
#include "parameters.h"
#include <cblas.h>

void cgsolver( double *A, double *b, double *x, int m, int n );
double * init_source_term(int n, double h);
void smvm(int m, const double* val, const int* col, const int* row, const double* x, double* y);
void cgsolver_sparse( double *Aval, int *I, int *J, double *b, double *x, int n);

#endif /*CG_H_*/


text:
///
/// @file
///
/// @brief This file contains functions and macros that are used in the sanity
/// checks.
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @internal 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.
///

#ifndef STARNEIG_COMMON_SANITY
#define STARNEIG_COMMON_SANITY

#include <starneig_config.h>
#include <starneig/configuration.h>

#ifdef STARNEIG_ENABLE_SANITY_CHECKS

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

///
/// @brief Reports a sanity check error. Aborts the program.
///
/// @param[in] message
///         Message.
///
#define STARNEIG_SANITY_REPORT(message) { \
    fprintf(stderr, "[starneig][sanity] %s:%d: %s\n", \
        __FILE__, __LINE__, message); \
    abort(); \
}

///
/// @brief Reports a sanity check error. Aborts the program.
///
/// @param[in] message
///         printf formatted message.
///
/// @param[in] ...
///         Additional printf compatible arguments.
///
#define STARNEIG_SANITY_REPORT_ARGS(message, ...) { \
    fprintf(stderr, "[starneig][sanity] %s:%d: %s\n", \
        __FILE__, __LINE__, message, __VA_ARGS__); \
    abort(); \
}

///
/// @brief Checks a conditional statement and reports a sanity check error if
/// the condition is not satisfied. Aborts the program.
///
/// @param[in] cond
///         The conditional statement.
///
/// @param[in] message
///         Message.
///
#define STARNEIG_SANITY_CHECK(cond, message) { \
    if (!(cond)) { \
        fprintf(stderr, "[starneig][sanity] %s:%d: %s\n", \
            __FILE__, __LINE__, message); \
        abort(); \
    } \
}

///
/// @brief Checks a conditional statement and reports a sanity check error if
/// the condition is not satisfied. Aborts the program.
///
/// @param[in] cond
///         The conditional statement.
///
/// @param[in] message
///         printf formatted message.
///
/// @param[in] ...
///         Additional printf compatible arguments.
///
#define STARNEIG_SANITY_CHECK_ARGS(cond, message, ...) { \
    if (!(cond)) { \
        fprintf(stderr, "[starneig][sanity] %s:%d: %s\n", \
            __FILE__, __LINE__, message, __VA_ARGS__); \
        abort(); \
    } \
}

static inline void starneig_sanity_check_inf(
    int rbegin, int rend, int cbegin, int cend, int ldA, double const *A,
    char const *mat, char const *file, int line)
{
    if (A == NULL)
        return;

    for (int i = cbegin; i < cend; i++) {
        for (int j = rbegin; j < rend; j++) {
            if (isinf(A[i*ldA+j])) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix %s has an infinite "
                    "element.\n", file, line, mat);
                abort();
            }
            if (isnan(A[i*ldA+j])) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix %s has a NaN element.\n",
                    file, line, mat);
                abort();
            }
        }
    }
}

#define STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, mat) \
    starneig_sanity_check_inf(rbegin, rend, cbegin, cend, ldA, A, mat, \
        __FILE__, __LINE__)

static inline void starneig_sanity_check_multiplicities(
    int n, double *real, double *imag, char const *file, int line)
{
    for (int i = 0; i < n; i++) {
        double lambda_re = real[i];
        double lambda_im = imag[i];

        if (lambda_im == 0.0) { // real eigenvalue
            for (int j = i+1; j < n; j++) {
                if (real[j] == lambda_re && imag[j] == 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Multiple eigenvalues at "
                        "S(%d,%d) and S(%d,%d).\n", file, line, i, i, j, j);
                    abort();
                }
            }
        }
        else { // complex eigenvalue
            for (int j = i+1; j < n; j++) {
                if (real[j] == lambda_re && imag[j] == lambda_im) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Multiple eigenvalues at "
                        "S(%d:%d,%d:%d) and S(%d:%d,%d:%d).\n",
                        file, line, i, i+1, i, i+1, j, j+1, j, j+1);
                    abort();
                }
            }
        }
    }
}

///
/// @brief A sanity check that checks for multiple eigenvalues.
///
/// @param[in] n
///         The length of the vectors real and imag.
///
/// @param[in] real
///         A vector that contains the real parts of the eigenvalues.
///
/// @param[in] imag
///         A vector that contains the imaginary parts of the eigenvalues.
///
#define STARNEIG_SANITY_CHECK_MULTIPLICITIES(n, real, imag) \
    starneig_sanity_check_multiplicities(n, real, imag, __FILE__, __LINE__)

static inline void starneig_sanity_check_orthogonality(
    int n, int ldQ, double const *Q, char const *mat, char const *file,
    int line)
{
    if (Q == NULL)
        return;

    starneig_sanity_check_inf(0, n, 0, n, ldQ, Q, mat, file, line);

    size_t ldT;
    double *T = starneig_alloc_matrix(n, n, sizeof(double), &ldT);

    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 (10000 < norm || isnan(norm)) {
        fprintf(stderr,
            "[starneig][sanity] %s:%d: Matrix %s is not orthogonal.\n",
            file, line, mat);
        fprintf(stderr,
            "[starneig][sanity] %s:%d: |%s %s^T - I| / |I| = %.0f u.\n",
            file, line, mat, mat, norm);
        abort();
    }

    starneig_free_matrix(T);
}

///
/// @brief A sanity check that makes sure a matrix Q is orthogonal.
///
/// @param[in] n
///         The order of the matrix Q.
///
/// @param[in] ldQ
///         The leading dimension of the matrix Q.
///
/// @param[in] Q
///         The matrix Q.
///
/// @param[in] mat
///         A string that identifies the matrix Q.
///
#define STARNEIG_SANITY_CHECK_ORTHOGONALITY(n, ldQ, Q, mat) \
    starneig_sanity_check_orthogonality(n, ldQ, Q, mat, __FILE__, __LINE__)

struct starneig_sanity_check_args {
    double *A; size_t ldA;
    double *B; size_t ldB;
};

static inline struct starneig_sanity_check_args *
starneig_sanity_check_residuals_begin(
    int n, int ldQ, int ldZ, int ldA, int ldB, double const *Q, double const *Z,
    double const *A, double const *B, char const *file, int line)
{
    starneig_sanity_check_inf(0, n, 0, n, ldA, A, "A", file, line);
    starneig_sanity_check_inf(0, n, 0, n, ldB, B, "B", file, line);
    starneig_sanity_check_orthogonality(n, ldQ, Q, "Q", file, line);
    starneig_sanity_check_orthogonality(n, ldZ, Z, "Z", file, line);

    if (Z == NULL) {
        Z = Q; ldZ = ldQ;
    }

    struct starneig_sanity_check_args *ret =
        malloc(sizeof(struct starneig_sanity_check_args));

    size_t ldT;
    double *T = starneig_alloc_matrix(n, n, sizeof(double), &ldT);

    ret->A = starneig_alloc_matrix(n, n, sizeof(double), &ret->ldA);
    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, ret->A, ret->ldA);

    ret->B = NULL; ret->ldB = 0;
    if (B != NULL) {
        ret->B = starneig_alloc_matrix(n, n, sizeof(double), &ret->ldB);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, n, n, n, 1.0,
            Q, ldQ, B, ldB, 0.0, T, ldT);
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0,
            T, ldT, Z, ldZ, 0.0, ret->B, ret->ldB);
    }

    starneig_free_matrix(T);

    starneig_sanity_check_inf(0, n, 0, n, ret->ldA, ret->A, "A", file, line);
    starneig_sanity_check_inf(0, n, 0, n, ret->ldB, ret->B, "B", file, line);

    return ret;
}

///
/// @brief The first part of a sanity check that makes sure a matrix pencil
/// Q (A,B) X^T and an updated matrix pencil ~Q (~A,~B) ~X^T are equivalent.
///
/// @param[in] name
///         An unique identifier.
///
/// @param[in] n
///         The order of the matrices A, B, Q and Z.
///
/// @param[in] ldQ
///         The leading dimension of the matrix Q.
///
/// @param[in] ldZ
///         The leading dimension of the matrix Z.
///
/// @param[in] ldA
///         The leading dimension of the matrix A.
///
/// @param[in] ldB
///         The leading dimension of the matrix B.
///
/// @param[in] Q
///         The matrix Q.
///
/// @param[in] Z
///         The matrix Z.
///
/// @param[in] A
///         The matrix A.
///
/// @param[in] B
///         The matrix B.
///
#define STARNEIG_SANITY_CHECK_RESIDUALS_BEGIN( \
    name, n, ldQ, ldZ, ldA, ldB, Q, Z, A, B) \
    struct starneig_sanity_check_args *name = \
    starneig_sanity_check_residuals_begin( \
        n, ldQ, ldZ, ldA, ldB, Q, Z, A, B, __FILE__, __LINE__)

static inline void starneig_sanity_check_residuals_end(
    int n, int ldQ, int ldZ, int ldA, int ldB, double const *Q, double const *Z,
    double const *A, double const *B,
    struct starneig_sanity_check_args const *args, char const *file, int line)
{
    struct starneig_sanity_check_args *ret =
        starneig_sanity_check_residuals_begin(
            n, ldQ, ldZ, ldA, ldB, Q, Z, A, B, file, line);

    int failure = 0;

    {
        double dot = 0.0;
        for(int i = 0; i < n; ++i )
            for(int j = 0; j < n; ++j )
                dot += squ(ret->A[ret->ldA*i+j] - args->A[args->ldA*i+j]);

        double a_dot = 0.0;
        for(int i = 0; i < n; ++i )
            for(int j = 0; j < n; ++j )
                a_dot += squ(args->A[args->ldA*i+j]);

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

        if (10000 < norm || isnan(norm)) {
            fprintf(stderr,
                "[starneig][sanity] %s:%d: Residual check failed for the "
                "matrix A.\n", file, line);
            fprintf(stderr,
                "[starneig][sanity] %s:%d: The norm was %.0f u.\n",
                file, line, norm);
            failure++;
        }
    }

    if (B != NULL) {
        double dot = 0.0;
        for(int i = 0; i < n; ++i )
            for(int j = 0; j < n; ++j )
                dot += squ(ret->B[ret->ldB*i+j] - args->B[args->ldB*i+j]);

        double a_dot = 0.0;
        for(int i = 0; i < n; ++i )
            for(int j = 0; j < n; ++j )
                a_dot += squ(args->B[args->ldB*i+j]);

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

        if (10000 < norm || isnan(norm)) {
            fprintf(stderr,
                "[starneig][sanity] %s:%d: Residual check failed for the "
                "matrix B.\n", file, line);
            fprintf(stderr,
                "[starneig][sanity] %s:%d: The norm was %.0f u.\n",
                file, line, norm);
            failure++;
        }
    }

    if (failure)
        abort();

    starneig_free_matrix(ret->A);
    starneig_free_matrix(ret->B);
    free(ret);
}

///
/// @brief The second part of a sanity check that makes sure a matrix pencil
/// Q (A,B) X^T and an updated matrix pencil ~Q (~A,~B) ~X^T are equivalent.
///
/// @param[in] name
///         An unique identifier.
///
/// @param[in] n
///         The order of the matrices A, B, Q and Z.
///
/// @param[in] ldQ
///         The leading dimension of the matrix Q.
///
/// @param[in] ldZ
///         The leading dimension of the matrix Z.
///
/// @param[in] ldA
///         The leading dimension of the matrix A.
///
/// @param[in] ldB
///         The leading dimension of the matrix B.
///
/// @param[in] Q
///         The matrix Q.
///
/// @param[in] Z
///         The matrix Z.
///
/// @param[in] A
///         The matrix A.
///
/// @param[in] B
///         The matrix B.
///
#define STARNEIG_SANITY_CHECK_RESIDUALS_END( \
    name, n, ldQ, ldZ, ldA, ldB, Q, Z, A, B) \
    starneig_sanity_check_residuals_end( \
        n, ldQ, ldZ, ldA, ldB, Q, Z, A, B, name, __FILE__, __LINE__); \
    if (name != NULL) { \
        starneig_free_matrix(name->A); \
        starneig_free_matrix(name->B); \
        free(name); \
        name = NULL; \
    }

///
/// @brief Skips the second part of a sanity check that makes sure a matrix
/// pencil Q (A,B) X^T and an updated matrix pencil ~Q (~A,~B) ~X^T are
/// equivalent.
///
/// @param[in] name
///         An unique identifier.
///
#define STARNEIG_SANITY_CHECK_RESIDUALS_SKIP(name) \
    if (name != NULL) { \
        starneig_free_matrix(name->A); \
        starneig_free_matrix(name->B); \
        free(name); \
        name = NULL; \
    }

static inline void starneig_sanity_check_bulges(
    int begin, int shifts, int n, int ldA, int ldB,
    double const *A, double const *B, char const *file, int line)
{
    starneig_sanity_check_inf(
        0, n, begin, begin+3*(shifts/2)+1, ldA, A, "A", file, line);
    starneig_sanity_check_inf(
        0, n, begin, begin+3*(shifts/2)+1, ldB, B, "B", file, line);

    if (n < begin+3*(shifts/2)+1) {
        fprintf(stderr,
            "[starneig][sanity] %s:%d: Matrix A has an invalid bulge.\n",
            file, line);
        abort();
    }

    double const *_A = A+begin*ldA+begin;
    double const *_B = B != NULL ? B+begin*ldB+begin : NULL;
    int _n = n - begin;
    for (int i = 0; i < shifts/2; i++) {
        for (int j = 3*i+4; j < _n; j++) {
            if (_A[3*i*ldA+j] != 0.0 || _A[(3*i+1)*ldA+j] != 0.0 ||
            _A[(3*i+2)*ldA+j] != 0.0) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix A has an invalid "
                    "bulge.\n", file, line);
                abort();
            }
        }
        for (int j = 3*i+3; _B != NULL && j < _n; j++) {
            if (_B[3*i*ldB+j] != 0.0 || _B[(3*i+1)*ldB+j] != 0.0 ||
            _B[(3*i+2)*ldB+j] != 0.0) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix B has invalid bulge .\n",
                    file, line);
                abort();
            }
        }
        if (_B != NULL &&
        (_B[3*i*ldB+3*i+1] != 0.0 || _B[3*i*ldB+3*i+2] != 0.0)) {
            fprintf(stderr,
                "[starneig][sanity] %s:%d: Matrix B has invalid bulge.\n",
                file, line);
            abort();
        }
    }
}

///
/// @brief A sanity check that makes sure that matrix pencil (A,B) contains the
/// correct number of bulges.
///
/// @param[in] begin
///         The column that should contain the first bulge.
///
/// @param[in] shifts
///         The number of shifts (shifts/2 bulges).
///
/// @param[in] n
///         The order of the matrices A and B.
///
/// @param[in] ldA
///         The leading dimension of the matrix A.
///
/// @param[in] ldB
///         The leading dimension of the matrix B.
///
/// @param[in] A
///         The matrix A.
///
/// @param[in] B
///         The matrix B.
///
#define STARNEIG_SANITY_CHECK_BULGES(begin, shifts, n, ldA, ldB, A, B) \
    starneig_sanity_check_bulges( \
        begin, shifts, n, ldA, ldB, A, B, __FILE__, __LINE__)

static inline void starneig_sanity_check_schur(
    int begin, int end, int n, int ldA, int ldB,
    double const *A, double const *B, char const *file, int line)
{
    starneig_sanity_check_inf(0, n, begin, end, ldA, A, "A", file, line);
    starneig_sanity_check_inf(0, n, begin, end, ldB, B, "B", file, line);

    const double safmin = dlamch("S");

    int two_by_two = 0;
    for (int i = begin; i < end; i++) {
        if (i+1 < end && A[i*ldA+i+1] != 0.0) {
            if (two_by_two) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix A is not in Schur "
                    "form.\n", file, line);
                abort();
            }

            if (B != NULL) {
                if (B[(i+1)*ldB+i] != 0.0 || B[i*ldB+i+1] != 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix pencil (A,B) "
                        "contains a non-normalized 2-by-2 block.\n",
                        file, line);
                    abort();
                }

                if (B[i*ldB+i] == 0.0 || B[(i+1)*ldB+i+1] == 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix B is singular.",
                        file, line);
                    abort();
                }

                double s1, s2, wr1, wr2, wi;

                extern void dlag2_(double const *, int const *, double const *,
                    int const *, double const *, double *, double *,
                    double *, double *, double *);

                dlag2_(&A[i*ldA+i], &ldA, &B[i*ldB+i], &ldB, &safmin,
                    &s1, &s2, &wr1, &wr2, &wi);

                if (wi == 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix pencil (A,B) "
                        "contains a fake 2-by-2 block.\n", file, line);
                    abort();
                }
            }
            else {
                if (A[i*ldA+i] != A[(i+1)*ldA+i+1]) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix A contains a "
                        "non-normalized 2-by-2 block.\n",
                        file, line);
                    abort();
                }

                double a[] = {
                    A[i*ldA+i], A[i*ldA+i+1], A[(i+1)*ldA+i], A[(i+1)*ldA+i+1]
                };

                double rt1r, rt1i, rt2r, rt2i, cs, ss;

                extern void dlanv2_(
                    double *, double *, double *, double *, double *,
                    double *, double *, double *, double *, double *);

                dlanv2_(&a[0], &a[2], &a[1], &a[3],
                    &rt1r, &rt1i, &rt2r, &rt2i, &cs, &ss);

                if (rt1i == 0.0 || rt2i == 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix A contains a fake "
                        "2-by-2 block.\n", file, line);
                    abort();
                }
            }

            two_by_two = 1;
        }
        else {
            two_by_two = 0;
        }

        for (int j = i+2; j < n; j++) {
            if (A[i*ldA+j] != 0.0) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix A is not in Schur "
                    "form.\n", file, line);
                abort();
            }
        }
    }

    if (B != NULL) {
        for (int i = begin; i < end; i++) {
            for (int j = i+2; j < n; j++) {
                if (B[i*ldB+j] != 0.0) {
                    fprintf(stderr,
                        "[starneig][sanity] %s:%d: Matrix B is not in upper "
                        "Hessenberg form.\n", file, line);
                    abort();
                }
            }
        }
    }
}

///
/// @brief A sanity check that makes sure that matrix pencil (A,B) is in Schur
/// form.
///
/// @param[in] begin
///         The first column to check.
///
/// @param[in] end
///         The last column to check + 1;
///
/// @param[in] n
///         The order of the matrices A and B.
///
/// @param[in] ldA
///         The leading dimension of the matrix A.
///
/// @param[in] ldB
///         The leading dimension of the matrix B.
///
/// @param[in] A
///         The matrix A.
///
/// @param[in] B
///         The matrix B.
///
#define STARNEIG_SANITY_CHECK_SCHUR(begin, end, n, ldA, ldB, A, B) \
    starneig_sanity_check_schur( \
        begin, end, n, ldA, ldB, A, B, __FILE__, __LINE__)

static inline void starneig_sanity_check_hessenberg(
    int begin, int end, int n, int ldA, int ldB,
    double const *A, double const *B, char const *file, int line)
{
    starneig_sanity_check_inf(0, n, begin, end, ldA, A, "A", file, line);
    starneig_sanity_check_inf(0, n, begin, end, ldB, B, "B", file, line);

    for (int i = begin; i < end; i++) {
        for (int j = i+2; j < n; j++) {
            if (A[i*ldA+j] != 0.0) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix A is not in upper "
                    "Hessenberg form.\n", file, line);
                abort();
            }
        }
    }
    for (int i = begin; B != NULL && i < end; i++) {
        for (int j = i+1; j < n; j++) {
            if (B[i*ldB+j] != 0.0) {
                fprintf(stderr,
                    "[starneig][sanity] %s:%d: Matrix B is not in upper "
                    "triangular form.\n", file, line);
                abort();
            }
        }
    }
}

///
/// @brief A sanity check that makes sure that matrix pencil (A,B) is in
/// Hessenberg-triangular form.
///
/// @param[in] begin
///         The first column to check.
///
/// @param[in] end
///         The last column to check + 1;
///
/// @param[in] n
///         The order of the matrices A and B.
///
/// @param[in] ldA
///         The leading dimension of the matrix A.
///
/// @param[in] ldB
///         The leading dimension of the matrix B.
///
/// @param[in] A
///         The matrix A.
///
/// @param[in] B
///         The matrix B.
///
#define STARNEIG_SANITY_CHECK_HESSENBERG(begin, end, n, ldA, ldB, A, B) \
    starneig_sanity_check_hessenberg( \
        begin, end, n, ldA, ldB, A, B, __FILE__, __LINE__)

#else

#define STARNEIG_SANITY_REPORT(message) {}
#define STARNEIG_SANITY_REPORT_ARGS(message, ...) {}
#define STARNEIG_SANITY_CHECK(cond, message) {}
#define STARNEIG_SANITY_CHECK_ARGS(cond, message, ...) {}
#define STARNEIG_SANITY_CHECK_INF(rbegin, rend, cbegin, cend, ldA, A, mat) {}
#define STARNEIG_SANITY_CHECK_MULTIPLICITIES(n, real, imag) {}
#define STARNEIG_SANITY_CHECK_ORTHOGONALITY(n, ldQ, Q, mat) {}
#define STARNEIG_SANITY_CHECK_RESIDUALS_BEGIN( \
    name, n, ldQ, ldZ, ldA, ldB, Q, Z, A, B) {}
#define STARNEIG_SANITY_CHECK_RESIDUALS_END( \
    name, n, ldQ, ldZ, ldA, ldB, Q, Z, A, B) {}
#define STARNEIG_SANITY_CHECK_RESIDUALS_SKIP(name) {}
#define STARNEIG_SANITY_CHECK_BULGES(begin, shifts, n, ldA, ldB, A, B) {}
#define STARNEIG_SANITY_CHECK_SCHUR(begin, end, n, ldA, ldB, A, B) {}
#define STARNEIG_SANITY_CHECK_HESSENBERG(begin, end, n, ldA, ldB, A, B) {}

#endif // STARNEIG_ENABLE_SANITY_CHECKS

#endif // STARNEIG_COMMON_SANITY

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig