text:
#ifndef DEINOS_ALGORITHM_H
#define DEINOS_ALGORITHM_H
#include "chess.h"
#include <vector>
#include <array>
#include <gsl/pointers>
#include <mutex>
#include <variant>
#include <functional>
#include <string>
#include <thread>
#include <future>
#include <condition_variable>

//memory consmption ideas:
//can halve Move size by removing piece type tracking and putting boold inside promotion piece
//can reduce AnalysedPosition size by doubling up control inside one byte
//perhaps dynamically delete parts of the tree e.g. only positions and regenerate if necessary
//deletion thread?

//Other
//use gsl::index in for loops?

namespace algorithm {
	class AnalysedPosition{
	public:
		constexpr AnalysedPosition() = default;
		explicit AnalysedPosition(const chess::Position&); //generate from scratch

		void advance_by(chess::MoveRecord);
		struct occlusion_info {
			std::array<std::array<chess::Square, 16>, 2> squares;
			std::array<int, 2> counts = {0};
		};
		AnalysedPosition::occlusion_info get_occlusion(chess::MoveRecord);

		inline const chess::Position& pos() const {return m_position;}
		inline uint8_t ctrl(chess::Almnt a, chess::Square s) const {return m_control[chess::as_index(a)].get(s.file(), s.rank());}
		inline const std::vector<chess::MoveRecord>& moves(chess::Almnt a) const {return m_moves[chess::as_index(a)];}
		inline const std::vector<chess::MoveRecord>& moves() const {return moves(pos().to_move());}
		inline chess::Move get_move(int index) const {return chess::Move(pos(), moves().at(index));}
		inline const chess::Square king_sq(chess::Almnt a) const {return m_king_sq[chess::as_index(a)];}
		inline bool in_check(chess::Almnt a) const {return ctrl(!a, king_sq(a)) > 0;}
		inline bool legal_check() const {return in_check(pos().to_move());}
		inline bool illegal_check() const {return in_check(!pos().to_move());}
		std::optional<chess::Move> find_record(const std::string& name) const;
		friend std::ostream& operator<<(std::ostream& os, const AnalysedPosition& ap);

	private:
		void append_calculation(chess::Square start); //calculate data associated with this square and append to state (for initialisation)
		void append_castling();
		std::array<chess::Square, 2> m_king_sq;
		chess::Position m_position;
		std::array<chess::HalfByteBoard, 2> m_control;
		std::array<std::vector<chess::MoveRecord>, 2> m_moves = {};
	};
	std::ostream& operator<<(std::ostream& os, const AnalysedPosition& ap);

	class Node;

	struct Edge {
		//std::mutex ptr_mutex;
		float total_value;
		int visits = 0;
		//bool legal = true;
		inline bool legal() const {return visits != -1;}
		std::unique_ptr<Node> node = nullptr;
	};

	struct EdgeDatum {
		float total_value = 0.0;
		int visits = 0;
		bool legal = true;
	};

	struct EdgeData {
		EdgeData(int moves_num) : edges(std::vector<EdgeDatum>(moves_num)) {}
		int total_n = 1;
		std::vector<EdgeDatum> edges;
	};

	class Node { //TODO
	public:
		explicit Node(std::unique_ptr<const AnalysedPosition> t_apos);
		explicit Node(const AnalysedPosition& t_apos);
		int preferred_index();
		std::unique_ptr<Node>* find_child(const std::string& fen);
		std::unique_ptr<Node>* find_child(const chess::Move& mv);
		inline Node* child(int index) const {return edges[index].node.get();}
		const chess::Move best_move() {return chess::Move(apos->pos(), apos->moves()[preferred_index()]);}
		inline std::optional<chess::GameResult> result() const {return m_result;}
		//inline int res_dist() const {return m_res_dist;} //TODO
		std::string display() const;

		inline int total_n() const {return m_total_n;} //synchronise?

		const std::unique_ptr<const AnalysedPosition> apos;
		const bool white_to_play;

	private:
		void update(int index, float t_value);
		void increment_n();
		void set_illegal(int index);

		std::optional<chess::GameResult> m_result = std::nullopt;
		//int m_res_dist = 0; //TODO
		//std::mutex data_mutex2;
		//EdgeData data;
		int m_total_n = 1;
		std::vector<Edge> edges;
		//std::vector<Edge> edges2;
		std::mutex data_mutex; //40B

		std::array<uint8_t, 3> node_prefetch = {0};

		friend class Tree;
	};

	class Tree {
	public:
		Tree(const AnalysedPosition& base_apos, std::function<float(const AnalysedPosition&)> t_value_fn,
			std::function<float(const AnalysedPosition&, const chess::Move&)> t_prior_fn, float t_expl_c = 0.2)
			: base(std::make_unique<Node>(base_apos)), value_fn(t_value_fn), prior_fn(t_prior_fn), expl_c(t_expl_c) {}
		void search(bool update_prefetch = false);
		std::unique_ptr<Node> base;
		std::function<float(const AnalysedPosition&)> value_fn;
		std::function<float(const AnalysedPosition&, const chess::Move&)> prior_fn;
		float expl_c = 0.2; //exploration coefficient

	private:
		float evaluate_node(Node& node); //used in search() for recursion
		void update_node(int index, float t_value);
		std::optional<int> edge_to_search(Node& node);
		void update_prefetch(Node& node);
	};

	class TreeEngine {
	public:
		TreeEngine(
			const AnalysedPosition& initial_position,
			std::function<float(const AnalysedPosition&)> t_value_fn,
			std::function<float(const AnalysedPosition&, const chess::Move&)> t_prior_fn,
			float exploration_coefficient);

		~TreeEngine() {
			halt_promise.set_value();
			for (auto& t : m_threads) t.join();
		}

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

		//void start();
		const chess::Move choose_move(); //maybe add exploration?
		bool advance_to(const std::string& fen);
		bool advance_by(const chess::Move& mv); //TODO
		std::string display() const; //TODO

		inline int total_n() const {return m_tree.base->total_n();};

	private:
		bool should_pause();
		void pause(); //blocks until all threads are paused
		void resume(); //blocks until all threads are resumed

		Tree m_tree;
		std::promise<void> halt_promise;
		std::mutex pause_mx;
		std::condition_variable pause_cv;
		bool pause_bool = false;
		int pause_count = 0;
		std::vector<std::thread> m_threads;
	};
}
#endif


text:
/*
 * Gsl_vector: test gsl_vector
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_vector.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_vector.c $(pkg-config --libs gsl)
 */

#include <stdio.h>
#include <gsl/gsl_vector.h>

static double v2;

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

    gsl_vector * v = gsl_vector_calloc(3);

    gsl_vector_set_all(v, 84.0);

    gsl_vector_fprintf(stdout, v, "%f");

    gsl_vector_free(v);

    return 0;
}

ext: c

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
//==============================================
//---------- by Dmitry R. Gulevich -------------
//--------- drgulevich@corp.ifmo.ru ------------
//--- ITMO University, St Petersburg, Russia ---
//==============================================
// Translation to C++ and augmentation
// S. R. Whiteley, wrcad.com,  Synopsys, Inc
//==============================================
//
// License:  GNU General Public License Version 3, 29 June 2007.

#ifndef MMJCO_H
#define MMJCO_H

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <complex>
#include <cmath>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

// Uncomment this to use the cminpack library for the least-squares
// fitting, otherwise we will use the cmpfit library provided.
//#define WITH_CMINPACK

#define MM_ECHG 1.602176634e-19     // Electron charge (SI units).
#define MM_BOLTZ 1.3806226e-23      // Boltzmann constant (SI units).

using namespace std;

class mmjco
{
public:
    // Constructor/destructor.
    mmjco(double, double, double, double);
    ~mmjco();

    // Functions return values for arbitrary x (values x<0 obtained by
    // symmetry), not smoothed.

    complex<double> Jpair(double x)
        {
            double absx = maximum(fabs(x), 1e-5);
            return (Repair(absx) + sign(x)*1j*Impair(absx));
        }

    complex<double> Jqp(double x)
        {
            double absx = maximum(fabs(x), 1e-5);
            return (Reqp(absx) + sign(x)*1j*Imqp(absx));
        }

    // As above, but smoothed.

    complex<double> Jpair_smooth(double x)
        {
            return (Jpair(x) + Jpair_correction(x));
        }

    complex<double> Jqp_smooth(double x)
        {
            return (Jqp(x) + Jqp_correction(x));
        }

    // Utilities.

    static double sign(double x)
        {
            if (x > 0.0)
                return (1.0);
            if (x == 0.0)
                return (0.0);
            return (-1.0);
        }

    static double minimum(double x, double y) { return (x < y ? x : y); }
    static double maximum(double x, double y) { return (x > y ? x : y); }

    static void save_data(const char*, const double*, const complex<double>*,
        const complex<double>*, int);

    static const char *version()    { return (mm_version); }

private:
    double Repair(double);
    double Impair(double);
    double Reqp(double);
    double Imqp(double);

    double Repair_integrand_part1(double y, double x)
        {
            return (tanh(mm_b*fabs(y))/
                ( sqrt(mm_dd1-(y-x)*(y-x)) * sqrt(y*y-mm_dd2) ));
        }

    double Repair_integrand_part2(double y, double x)
        {
            return (tanh(mm_b*fabs(y))/
                ( sqrt(y*y-mm_dd1) * sqrt(mm_dd2-(y+x)*(y+x)) ));
        }

    double Impair_integrand(double y, double x)
        {
            return (( tanh(mm_b*(y+x)) - tanh(mm_b*y) )*sign(y)*sign(y+x)/(
                sqrt(y*y-mm_dd1) * sqrt((y+x)*(y+x)-mm_dd2) ));
        }

    double Reqp_integrand_part1(double y, double x)
        {
            return (fabs(y)*tanh(mm_b*y)*(y-x)/( sqrt(y*y-mm_dd1) *
                sqrt(mm_dd2-(y-x)*(y-x)) ));
        }

    double Reqp_integrand_part1b_d1_xd2(double y, double x)
        {
            return (fabs(y)*tanh(mm_b*y)*(y-x)/sqrt((y+mm_d1)*(mm_d2+y-x)));
        }

    double Reqp_integrand_part1b_xd2_xd2(double y, double x)
        {
            return (fabs(y)*tanh(mm_b*y)*(y-x)/sqrt(y*y-mm_dd1));
        }

    double Reqp_integrand_part2_xd1_d2(double y, double x)
        {
            return (fabs(y)*tanh(mm_b*y)*(y+x)/sqrt((mm_d1-y-x)*(mm_d2-y)));
        }

    double Reqp_integrand_part2_xd1_xd1(double y, double x)
        {
            return (fabs(y)*tanh(mm_b*y)*(y+x)/sqrt(y*y-mm_dd2));
        }

    double Imqp_integrand(double y, double x)
        {
            return (( tanh(mm_b*(y+x))-tanh(mm_b*y) )*fabs(y)*fabs(y+x)/
                ( sqrt((y+x)*(y+x)-mm_dd1) * sqrt(y*y-mm_dd2) ));
        }

    // Passed to integration functions.

    static double re_p_i1(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Repair_integrand_part1(x, m->mm_arg));
        }

    static double re_p_i2(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Repair_integrand_part2(x, m->mm_arg));
        }

    static double im_p_i(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Impair_integrand(x, m->mm_arg));
        }

    static double re_q_i1(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Reqp_integrand_part1(x, m->mm_arg));
        }

    static double re_q_i1b1x2(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Reqp_integrand_part1b_d1_xd2(x, m->mm_arg));
        }

    static double re_q_i1bx2x2(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Reqp_integrand_part1b_xd2_xd2(x, m->mm_arg));
        }

    static double re_q_i2x12(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Reqp_integrand_part2_xd1_d2(x, m->mm_arg));
        }

    static double re_q_i2x1x1(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Reqp_integrand_part2_xd1_xd1(x, m->mm_arg));
        }

    static double im_q_i(double x, void *ptr)
        {
            mmjco *m = (mmjco*)ptr;
            return (m->Imqp_integrand(x, m->mm_arg));
        }

    // Smoothing.

    complex<double> Jpair_correction(double x)
        {
            if (mm_symj) {
                double absx = maximum(fabs(x),1.e-5);
                return (dRe(absx) + sign(x)*1j*(dIm(absx) + dIm_at_0(absx)));
            }
            else {
                double absx = maximum(fabs(x),1.e-5);
                return (dRe(absx) + dRe_minus(absx) + sign(x)*1j*(dIm(absx) +
                    dIm_minus(absx)));
            }
        }

    complex<double> Jqp_correction(double x)
        {

            if (mm_symj) {
                double absx = maximum(fabs(x),1.e-5);
                return (dRe(absx) + sign(x)*1j*(-dIm(absx) + dIm_at_0(absx)));
            }
            else {
                double absx = maximum(fabs(x),1.e-5);
                return (dRe(absx) - dRe_minus(absx) + sign(x)*1j*(-dIm(absx) +
                    dIm_minus(absx)));
            }
        }

    // Smoothing for Repair, Reqp.
    double dRe(double x)
        {
            double sqpos = (x-1.0)*(x-1.0);
            double sqneg = (x+1.0)*(x+1.0);
            return (-x*(1.0/M_PI)*mm_ip0*0.5*log(
                ((sqpos+mm_dsm*mm_dsm)/sqpos)*(sqneg/(sqneg+mm_dsm*mm_dsm)) ));
        }

    // Smoothing for Repair, -Reqp at x=mm_d2-mm_d1.
    double dRe_minus(double x)
        {
            return (M_PI*x*sqrt(mm_d1*mm_d2)*
                (tanh(mm_b*mm_d2)-tanh(mm_b*mm_d1)) *
                0.5*( (2.0/M_PI)*atan((x-mm_d21)/mm_dsm)
                - sign(x-mm_d21) + (2.0/M_PI)*atan((x+mm_d21)/mm_dsm) -
                sign(x+mm_d21) ) / (4.0*mm_d21));
        }

    // Smoothing for Impair, Imqp at x=mm_d2-mm_d1.
    double dIm_minus(double x)
        {
            double square1 = (x-mm_d21)*(x-mm_d21);
            double square2 = (x+mm_d21)*(x+mm_d21);
            return (-x*sqrt(mm_d1*mm_d2)*(tanh(mm_b*mm_d2)-tanh(mm_b*mm_d1))*
                0.5*log((square1+mm_dsm*mm_dsm)*(square2+mm_dsm*mm_dsm)/
                (square1*square2)) / (4.0*mm_d21));
        }

    // Smoothing for Impair, -Imqp.
    double dIm(double x)
        {
            return (x*0.5*mm_ip0*((2.0/M_PI)*atan((1.0-x)/mm_dsm) -
                sign(1.0-x) + (2.0/M_PI)*atan((1.0+x)/mm_dsm) - sign(1.0+x)));
        }

    // Smoothing for Impair, Imqp at mm_d2-mm_d1=0.
    double dIm_at_0(double x)
        {
            double x2=x*x;
            return (-mm_b*x*mm_expb*0.5*log(
                (x2+mm_dsm*mm_dsm)/x2)/((mm_expb+1.0)*(mm_expb+1.0)));
        }

    static const char *mm_version;

    double mm_d1, mm_d2;
    double mm_d21;
    double mm_dd1, mm_dd2;
    double mm_b;
    double mm_arg;
    int mm_limit;
    int mm_key;
    double mm_epsabs;
    double mm_epsrel;
    gsl_integration_workspace *mm_ws;
    gsl_integration_qaws_table *mm_tbl;

    // Smoothing.
    double mm_expb;
    double mm_ip0;
    double mm_dsm;
    bool mm_symj;;
};

//
// Class to compute optimized TCA fitting parameters.
//

// Largest acceptable fitting table.
#define MAX_NTERMS 20

class mmjco_fit
{
public:
    struct mm_adata
    {
        const double *x;
        double thr;
        const complex<double> *Jpair_data;
        const complex<double> *Jqp_data;
    };

    mmjco_fit()
        {
            mmf_pAB = 0;
            mmf_nterms = 0;
        }

    ~mmjco_fit()
        {
            delete [] mmf_pAB;
        }

    void new_fit_parameters(const double*, const complex<double>*,
        const complex<double>*, int, int, double);
    void save_fit_parameters(const char*);
    void load_fit_parameters(const char*);
    void tca_fit(const double*, int, complex<double>**, complex<double>**);
    double residual(const complex<double>*, const complex<double>*,
        const complex<double>*, const complex<double>*, int, double);

    const complex<double> *params()     const { return (mmf_pAB); }
    int numterms()                      const { return (mmf_nterms); }

private:
    static complex<double> *modelJpair(const complex<double>*, int,
        const double*, int);
    static complex<double> *modelJqp(const complex<double>*, int,
        const double*, int);

    static double rep(double zeta)     { return (-fabs(zeta)); }
    static double invrep(double xi)    { return (-xi); }

    // Relative difference with threshold thr.
    static double Drel(double X, double Xref, double thr)
        {
            double axref = fabs(Xref);
            if (thr > axref)
                axref = thr;
            return (fabs(X-Xref)/axref);
        }

    // Flatten complex array.
    static double *realimag(complex<double> *carray, int sz)
        {
            double *param_list = new double[2*sz];
            for (int i = 0; i < sz; i++) {
                param_list[i*2] = carray[i].real();
                param_list[i*2+1] = carray[i].imag();
            }
            return (param_list);
        }

    // Form complex array.
    static complex<double> *ccombine(const double *param_list, int plsz)
        {
            complex<double> *cpars = new complex<double>[plsz/2];
            for (int i = 0; i < plsz; i++) {
                int n = i/2;
                if (i & 1)
                    cpars[n].imag(param_list[i]);
                else
                    cpars[n].real(param_list[i]);
            }
            return (cpars);
        }

#ifdef WITH_CMINPACK
    static int func(void*, int, int, const double*, double*, int);
#else
    static int func(int, int, double*, double*, double**, void*);
#endif

    complex<double> *mmf_pAB;
    int mmf_nterms;
};

#endif // MMJCO_H

ext: h

max_stars_repo_name: wrcad/xictools

max_issues_repo_name: markvolkmann/xictools

max_forks_repo_name: markvolkmann/xictools


text:
/**
 *
 * @file qwrapper_zlansy.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlansy(Quark *quark, Quark_Task_Flags *task_flags,
                       int norm, PLASMA_enum uplo, int N,
                       const PLASMA_Complex64_t *A, int LDA, int szeA,
                       int szeW, double *result)
{
    szeW = max(1, szeW);
    DAG_CORE_LANSY;
    QUARK_Insert_Task(quark, CORE_zlansy_quark, task_flags,
        sizeof(PLASMA_enum),                &norm,  VALUE,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
        sizeof(int),                        &LDA,   VALUE,
        sizeof(double)*szeW,                 NULL,          SCRATCH,
        sizeof(double),                      result,        OUTPUT,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlansy_f1(Quark *quark, Quark_Task_Flags *task_flags,
                          PLASMA_enum norm, PLASMA_enum uplo, int N,
                          const PLASMA_Complex64_t *A, int LDA, int szeA,
                          int szeW, double *result,
                          double *fake, int szeF)
{
    szeW = max(1, szeW);
    DAG_CORE_LANSY;

    if ( result == fake ) {
        QUARK_Insert_Task(quark, CORE_zlansy_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double)*szeF,                 result,        OUTPUT | GATHERV,
            0);
    } else {
        QUARK_Insert_Task(quark, CORE_zlansy_f1_quark, task_flags,
            sizeof(PLASMA_enum),                &norm,  VALUE,
            sizeof(PLASMA_enum),                &uplo,  VALUE,
            sizeof(int),                        &N,     VALUE,
            sizeof(PLASMA_Complex64_t)*szeA,     A,             INPUT,
            sizeof(int),                        &LDA,   VALUE,
            sizeof(double)*szeW,                 NULL,          SCRATCH,
            sizeof(double),                      result,        OUTPUT,
            sizeof(double)*szeF,                 fake,          OUTPUT | GATHERV,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlansy_quark = PCORE_zlansy_quark
#define CORE_zlansy_quark PCORE_zlansy_quark
#endif
void CORE_zlansy_quark(Quark *quark)
{
    double *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;

    quark_unpack_args_7(quark, norm, uplo, N, A, LDA, work, normA);
    *normA = LAPACKE_zlansy_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm), lapack_const(uplo),
        N, A, LDA, work);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlansy_f1_quark = PCORE_zlansy_f1_quark
#define CORE_zlansy_f1_quark PCORE_zlansy_f1_quark
#endif
void CORE_zlansy_f1_quark(Quark *quark)
{
    double *normA;
    int norm;
    PLASMA_enum uplo;
    int N;
    PLASMA_Complex64_t *A;
    int LDA;
    double *work;
    double *fake;

    quark_unpack_args_8(quark, norm, uplo, N, A, LDA, work, normA, fake);
    *normA = LAPACKE_zlansy_work(
        LAPACK_COL_MAJOR,
        lapack_const(norm),
        lapack_const(uplo),
        N, A, LDA, work);
}


text:
#pragma once

#include "gl_helpers.h"
#include "GL_Loader.h"
#include "MappedGL.h"
#include "Range.h"

#include <handy/Guard.h>

#include <gsl/span>

#include <type_traits>
#include <vector>

namespace ad {
namespace graphics {

struct [[nodiscard]] VertexArrayObject : public ResourceGuard<GLuint>
{
    /// \note Deleting a bound VAO reverts the binding to zero
    VertexArrayObject() :
        ResourceGuard<GLuint>{reserve(glGenVertexArrays),
                              [](GLuint aIndex){glDeleteVertexArrays(1, &aIndex);}}
    {}
};

inline void bind(const VertexArrayObject & aVertexArray)
{
    glBindVertexArray(aVertexArray);
}

// TODO Ad 2022/02/02: Is it a good idea to "expect" the object to unbind
// when the underlying unbinding mechanism does not use it (just reset a default)?
inline void unbind(const VertexArrayObject & aVertexArray)
{
    glBindVertexArray(0);
}

/// \TODO understand when glDisableVertexAttribArray should actually be called
///       (likely not specially before destruction, but more when rendering other objects
///        since it is a current(i.e. global) VAO state)
/// \Note Well note even that: Activated vertex attribute array are per VAO, so changing VAO
///       Already correctly handles that.
struct [[nodiscard]] VertexBufferObject : public ResourceGuard<GLuint>
{
    VertexBufferObject() :
        ResourceGuard<GLuint>{reserve(glGenBuffers),
                              [](GLuint aIndex){glDeleteBuffers(1, &aIndex);}}
    {}
};

inline void bind(const VertexBufferObject & aVertexBuffer)
{
    glBindBuffer(GL_ARRAY_BUFFER, aVertexBuffer);
}

inline void unbind(const VertexBufferObject &)
{
    glBindBuffer(GL_ARRAY_BUFFER, 0);
}

struct [[nodiscard]] IndexBufferObject : public ResourceGuard<GLuint>
{
    IndexBufferObject() :
        ResourceGuard<GLuint>{reserve(glGenBuffers),
                              [](GLuint aIndex){glDeleteBuffers(1, &aIndex);}}
    {}
};

inline void bind(const IndexBufferObject & aIndexBuffer)
{
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, aIndexBuffer);
}

inline void unbind(const IndexBufferObject &)
{
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, 0);
}

/// \brief A VertexArray with a vector of VertexBuffers
struct [[nodiscard]] VertexSpecification
{
    VertexSpecification(VertexArrayObject aVertexArray={},
                        std::vector<VertexBufferObject> aVertexBuffers={}) :
        mVertexArray{std::move(aVertexArray)},
        mVertexBuffers{std::move(aVertexBuffers)}
    {}

    VertexArrayObject mVertexArray;
    std::vector<VertexBufferObject> mVertexBuffers;
};

inline void bind(const VertexSpecification & aVertexSpecification)
{
    bind(aVertexSpecification.mVertexArray);
}

inline void unbind(const VertexSpecification & aVertexSpecification)
{
    unbind(aVertexSpecification.mVertexArray);
}

/// \brief Describes the attribute access from the shader (layout id, value type, normalization)
struct Attribute
{
    enum class Access
    {
        Float,
        Integer,
    };

    constexpr Attribute(GLuint aValue) :
        mIndex(aValue)
    {}

    constexpr Attribute(GLuint aValue, Access aAccess, bool aNormalize=false) :
        mIndex(aValue),
        mTypeInShader(aAccess),
        mNormalize(aNormalize)
    {}

    GLuint mIndex; // index to match in vertex shader.
    Access mTypeInShader{Access::Float}; // destination data type
    bool mNormalize{false}; // if destination is float and source is integral, should it be normalized (value/type_max_value)
};

/// \brief The complete description of an attribute expected by OpenGL
struct AttributeDescription : public Attribute
{
    GLuint mDimension;  // from 1 to 4 (explicit distinct attributes must be used for matrix data)
    size_t mOffset;     // offset for the attribute within the vertex data structure (interleaved)
    GLenum mDataType;   // attribute source data type
};

std::ostream & operator<<(std::ostream &aOut, const AttributeDescription & aDescription);

typedef std::initializer_list<AttributeDescription> AttributeDescriptionList;

/***
 *
 * Vertex Buffer
 *
 ***/

/// \brief Attach an existing VertexBuffer to an exisiting VertexArray,
/// without providing initial data.
void attachVertexBuffer(const VertexBufferObject & aVertexBuffer,
                        const VertexArrayObject & aVertexArray,
                        AttributeDescriptionList aAttributes,
                        GLsizei aStride,
                        GLuint aAttributeDivisor = 0);

/// \brief This overload deduces the stride from T_vertex.
template <class T_vertex>
void attachVertexBuffer(const VertexBufferObject & aVertexBuffer,
                        const VertexArrayObject & aVertexArray,
                        AttributeDescriptionList aAttributes,
                        GLuint aAttributeDivisor = 0)
{
    return attachVertexBuffer(aVertexBuffer, aVertexArray, aAttributes, sizeof(T_vertex), aAttributeDivisor);
}

/// \brief Intialize a VertexBufferObject, without providing initial data.
///
/// This is an extension of `attachVertexBuffer()`, which constructs the vertex buffer it attaches,
/// instead of expecting it as argument.
VertexBufferObject initVertexBuffer(const VertexArrayObject & aVertexArray,
                                    AttributeDescriptionList aAttributes,
                                    GLsizei aStride,
                                    GLuint aAttributeDivisor = 0);

/// \brief This overload deduces the stride from T_vertex.
template <class T_vertex>
VertexBufferObject initVertexBuffer(const VertexArrayObject & aVertexArray,
                                    AttributeDescriptionList aAttributes,
                                    GLuint aAttributeDivisor = 0)
{
    return initVertexBuffer(aVertexArray, aAttributes, sizeof(T_vertex), aAttributeDivisor);
}

/// \brief Create a VertexBufferObject with provided attributes, load it with data,
/// and associate the data to attributes of `aVertexArray`.
///
/// This is an extension of `initVertexBuffer()`, which loads data into the initialized vertex buffer.
///
/// \param aAttribute describes the associaiton.
///
/// \note This is the lowest level overload, with explicit attribute description and raw data pointer.
/// Other overloads end-up calling it.
VertexBufferObject loadVertexBuffer(const VertexArrayObject & aVertexArray,
                                    AttributeDescriptionList aAttributes,
                                    GLsizei aStride,
                                    size_t aSize,
                                    const GLvoid * aData,
                                    GLuint aAttributeDivisor = 0);

/// \brief This overload deduces the stride and size from T_vertex,
/// which could itself be deduced from the provided span.
template <class T_vertex>
VertexBufferObject loadVertexBuffer(const VertexArrayObject & aVertexArray,
                                    const AttributeDescriptionList & aAttributes,
                                    const gsl::span<T_vertex> aVertices,
                                    GLuint aAttributeDivisor = 0)
{
    return loadVertexBuffer(aVertexArray,
                            aAttributes,
                            sizeof(T_vertex),
                            aVertices.size_bytes(),
                            aVertices.data(),
                            aAttributeDivisor);
}

/// \brief Create a VertexBufferObject and load it with provided data.
/// But do **not** attach it to a VertexArrayObject / do **not** associate the vertex data to attributes.
///
/// \note Attachment to a VertexArrayObject as well as attributes association might be done later with
/// `attachVertexBuffer()`.
template <class T_vertex>
VertexBufferObject loadUnattachedVertexBuffer(gsl::span<const T_vertex> aVertices,
                                              GLenum aHint = GL_STATIC_DRAW)
{
    VertexBufferObject vbo;
    glBindBuffer(GL_ARRAY_BUFFER, vbo);
    glBufferData(GL_ARRAY_BUFFER, aVertices.size_bytes(), aVertices.data(), aHint);
    return vbo;
}

/// \brief High-level function directly appending a loaded VertexBuffer to a VertexSpecification.
///
/// This is an extension to `loadVertexBuffer()`, which appends the loaded vertex buffer to `aSpecification`.
template <class T_vertex>
void appendToVertexSpecification(VertexSpecification & aSpecification,
                                 const AttributeDescriptionList & aAttributes,
                                 const gsl::span<T_vertex> aVertices,
                                 GLuint aAttributeDivisor = 0)
{
    aSpecification.mVertexBuffers.push_back(
            loadVertexBuffer(aSpecification.mVertexArray,
                             aAttributes,
                             std::move(aVertices),
                             aAttributeDivisor));
}

/***
 *
 * Index Buffer
 *
 ***/

/// \brief Attach an existing IndexBuffer to an exisiting VertexArray,
/// without providing initial data.
inline void attachIndexBuffer(const IndexBufferObject & aIndexBuffer,
                              const VertexArrayObject & aVertexArray)
{
    glBindVertexArray(aVertexArray);
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, aIndexBuffer);
}

/// \brief Intialize and attach an IndexBufferObject, without providing initial data.
///
/// This is an extension of `attachIndexBuffer()`, which constructs the index buffer it attaches,
/// instead of expecting it as argument.
inline IndexBufferObject initIndexBuffer(const VertexArrayObject & aVertexArray)
{
    IndexBufferObject ibo;
    attachIndexBuffer(ibo, aVertexArray);
    return ibo;
}

/// \brief Initialize, attach and load data into an IndexBufferObject.
///
/// This is an extension of `initIndexBuffer()`, which loads data into the initialized vertex buffer.
template <class T_index>
IndexBufferObject loadIndexBuffer(const VertexArrayObject & aVertexArray,
                                  const gsl::span<T_index> aIndices,
                                  const BufferHint aHint)
{
    IndexBufferObject ibo = initIndexBuffer(aVertexArray);
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, aIndices.size_bytes(), aIndices.data(), getGLBufferHint(aHint));
    return ibo;
}

/***
 * Buffer re-specification
 *
 * see: https://www.khronos.org/opengl/wiki/Buffer_Object_Streaming#Buffer_re-specification
 ***/

/// \brief Respecify content of a vertex buffer.
inline void respecifyBuffer(const VertexBufferObject & aVBO, const GLvoid * aData, GLsizei aSize)
{
    glBindBuffer(GL_ARRAY_BUFFER, aVBO);

    // Orphan the previous buffer
    glBufferData(GL_ARRAY_BUFFER, aSize, NULL, GL_STATIC_DRAW);

    // Copy value to new buffer
    glBufferSubData(GL_ARRAY_BUFFER, 0, aSize, aData);
}

/// \brief Respecify content of an index buffer.
inline void respecifyBuffer(const IndexBufferObject & aIBO, const GLvoid * aData, GLsizei aSize)
{
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, aIBO);

    // Orphan the previous buffer
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, aSize, NULL, GL_STATIC_DRAW);

    // Copy value to new buffer
    glBufferSubData(GL_ELEMENT_ARRAY_BUFFER, 0, aSize, aData);
}

/// \brief Overload accepting a span of generic values, instead of low-level void pointer.
/// It works with both vertex and index buffers.
template <class T_values, class T_buffer>
void respecifyBuffer(const T_buffer & aBufferObject, const gsl::span<T_values> aValues)
{
    respecifyBuffer(aBufferObject,
                    aValues.data(),
                    static_cast<GLsizei>(aValues.size_bytes()));
}

/// \brief Respecify a vertex buffer with the exactly same size (allowing potential optimizations).
///
/// \attention This is undefined behaviour is aData does not point to at least the same amount
/// of data that was present before in the re-specified vertex buffer.
inline void respecifyBufferSameSize(const VertexBufferObject & aVBO, const GLvoid * aData)
{
    GLint size;
    glGetBufferParameteriv(GL_ARRAY_BUFFER, GL_BUFFER_SIZE, &size);

    respecifyBuffer(aVBO, aData, size);
}

} // namespace graphics
} // namespace ad

max_stars_repo_name: Adnn/Graphics

max_issues_repo_name: Adnn/Graphics

max_forks_repo_name: Adnn/Graphics


text:
#ifndef __INC_MODELSTRUCT__
#define __INC_MODELSTRUCT__

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

/**
 * @file modelstruct.h
 * \brief defines modelstruct which contains the raw data to be worked on by emulator and estimator
 */
struct optstruct;
//! holds the main data structures for a single model emulator run

/**
 * \struct modelstruct
 * this structure is everything you need to do run the estimation and the
 * emulation
 *
 * it contains:
 * training vector
 * model points
 * theta vector (hyperparameters)
 * the options used to define this model
 * \note: the options used to be kept separately, they are now always contained in the modelstruct
 */
typedef struct modelstruct{
	/**
	 * \brief x_values (nmodel_points * nparams)
	 *
	 * the points at which the model (that which is to be emulated) was evaluated,
	 * each row is a new model point, the columns are the dimensions
	 */
	gsl_matrix* xmodel;
	/**
	 * \brief y_values (nmodel_points)
	 *
	 * the values of the model evaluated at each model point
	 * this is what we aim to interpolate
	 */
	gsl_vector* training_vector;
	/**
	 * \brief hyperparameters (nthetas)
	 *
	 * working and eventually best values of the covariance fns (hyper)parameters which
	 * should force the emulated model to pass through all of the training points
	 */
	gsl_vector* thetas;
	/**
   * the average distances between samples in the design in each dimension
	 * this sets soft minima on the correlation function length
	 */
	gsl_vector* sample_scales;
	/**
	 * the options struct used to init this structure
	 */
	struct optstruct* options;

	/**
	 * fnptrs which define various aspects of the regression etc
	 * should these be in optstruct?
	 */

	/**
	 * setup the regression H vector, depends upon the regression order
	 *
	 * (needed in regression.h)
	 */
	void (*makeHVector)(gsl_vector *h_vector, gsl_vector *x_location, int nparams);

	/**
	 * points to the covariance function,
	 *
	 * (needed in emulator.h)
	 */
	double (*covariance_fn)(gsl_vector*, gsl_vector*, gsl_vector*, int, int);

	/**
	 * compute the gradient matrix for the length setting theta values
	 * dC/dTheta = (C-nugget) * (1/2)*(x_i - x_j)^(alpha) * alpha / (thetaLength)
	 *
	 * this is a fn-ptr which will be set when the cov function is setup
	 * the different target fns are in emulator.c called derivative_l_<covfnname>
	 *
	 * @param covsub the covariance matrix with the nugget term subtracted
	 * @param thetaLength the current value of the length scale we're differentiating wrt
	 * @param index, the direction we're looking in
	 * @return dCdTheta the matrix of the derivative of C wrt index
	 *
	 * (needed in maxmultimin.h)
	 */
	void (*makeGradMatLength)(gsl_matrix *dCdTheta,  gsl_matrix* xmodel,
														double thetaLength, int index, int nmodel_points, int nparams);

} modelstruct;

#include "optstruct.h"

void alloc_modelstruct(modelstruct* the_model, optstruct* options);
void free_modelstruct(modelstruct* the_model);
void fill_modelstruct(modelstruct* the_model, optstruct* options, char** input_data);

void load_modelstruct(FILE* fptr, modelstruct* the_model, optstruct* opts);
void dump_modelstruct(FILE* fptr, modelstruct* the_model, optstruct* opts);

void copy_modelstruct(modelstruct* dst, modelstruct* src);

modelstruct * alloc_modelstruct_2(gsl_matrix* xmodel, gsl_vector* training_vector, int cov_fn_index, int regression_order);
void free_modelstruct_2(modelstruct * model);
void dump_modelstruct_2(FILE *fptr, modelstruct* the_model);
modelstruct* load_modelstruct_2(FILE *fptr);

void set_global_ptrs(modelstruct * model);
gsl_vector * fill_sample_scales_vec(gsl_matrix* xmodel);

#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void free_MCMC_covariance(MCMC_info *MCMC) {
    int           i_DS;
    MCMC_DS_info *current_DS;
    MCMC_DS_info *next_DS;

    if(MCMC->n_M != NULL) {
        SID_log("Freeing MCMC covariance matrix...", SID_LOG_OPEN);

        SID_free(SID_FARG MCMC->V);
        if(MCMC->m != NULL) {
            gsl_matrix_free(MCMC->m);
            MCMC->m = NULL;
        }
        if(MCMC->b != NULL) {
            gsl_vector_free(MCMC->b);
            MCMC->b = NULL;
        }

        SID_log("Done.", SID_LOG_CLOSE);
    }
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#ifndef _GDRDGEMM_H
#define _GDRDGEMM_H

#ifndef MKL
#include <cblas.h>
#else
#include <mkl_cblas.h>
#endif

#define NMAT 2048
void gdr_check_and_restart(double a[][NMAT],
			   double b[][NMAT],
			   double c[][NMAT]);

void gdrblas_dgemm
(
#ifndef MKL
   const enum CBLAS_ORDER             ORDER,
   const enum CBLAS_TRANSPOSE             TRANSA,
   const enum CBLAS_TRANSPOSE             TRANSB,
#else
   const CBLAS_ORDER             ORDER,
   const CBLAS_TRANSPOSE             TRANSA,
   const CBLAS_TRANSPOSE             TRANSB,
#endif
   const int                        M,
   const int                        N,
   const int                        K,
   const double                     ALPHA,
   const double *                   A,
   const int                        LDA,
   const double *                   B,
   const int                        LDB,
   const double                     BETA,
   double *                         C,
   const int                        LDC
   );

void mygdrdgemm(int m,
		int n,
		int k,
		double alpha,
		double * a,
		int na,
		double * b,
		int nb,
		double beta,
		double * c,
		int nc);

#endif


text:
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <unistd.h>
#include <dirent.h>
#include "hdf5.h"
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_bessel.h>
#include "mclib_3d.h"
//#include "mclib.h"
#include <omp.h>

#define R_DIM 1260
#define THETA_DIM 280
#define PHI_DIM 280

void read_hydro(char hydro_prefix[200], int frame, double r_inj, double **x, double **y,  double **z, double **szx, double **szy, double **r, double **theta, double **phi,\
  double **velx, double **vely, double **velz, double **dens, double **pres, double **gamma, double **dens_lab, double **temp, int *number,  int ph_inj, double min_r, double max_r, double fps, FILE *fPtr)
{
    FILE *hydroPtr=NULL;
    char hydrofile[200]="", file_num[200]="", full_file[200]=""  ;
    char buf[10]="";
    int i=0, j=0, k=0, elem=0;
    float buffer=0;
    float *dens_unprc=malloc(sizeof(float)*R_DIM*THETA_DIM*PHI_DIM);
    float *vel_r_unprc=malloc(sizeof(float)*R_DIM*THETA_DIM*PHI_DIM);
    float *vel_theta_unprc=malloc(sizeof(float)*R_DIM*THETA_DIM*PHI_DIM);
    float *vel_phi_unprc=malloc(sizeof(float)*R_DIM*THETA_DIM*PHI_DIM);
    float *pres_unprc=malloc(sizeof(float)*R_DIM*THETA_DIM*PHI_DIM);
    double ph_rmin=0, ph_rmax=0;
    double r_in=1e10, r_ref=2e13;
    double *r_edge=malloc(sizeof(double)*(R_DIM+1));
    double *dr=malloc(sizeof(double)*(R_DIM));
    double *r_unprc=malloc(sizeof(double)*R_DIM);
    double *theta_unprc=malloc(sizeof(double)*THETA_DIM);
    double *phi_unprc=malloc(sizeof(double)*PHI_DIM);

    if (ph_inj==0)
    {
        ph_rmin=min_r;
        ph_rmax=max_r;
    }

    //density
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 1,"-" );
    modifyFlashName(file_num, hydrofile, frame, 1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, ".data");

    fprintf(fPtr,"Reading Density: %s\n", full_file);
    fflush(fPtr);

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(dens_unprc, sizeof(float)*R_DIM*THETA_DIM*PHI_DIM,R_DIM*THETA_DIM*PHI_DIM, hydroPtr); //data
    fclose(hydroPtr);

    for (i=0;i<R_DIM*THETA_DIM*PHI_DIM;i++)
    {
        if ((i>98784000-5) || (i<5))
        {
            fprintf(fPtr,"Density %d: %0.7e\n", i, *(dens_unprc+i));
            fflush(fPtr);
        }
    }

    //velocities divided by c
    //v_r
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 2,"-" );
    modifyFlashName(file_num, hydrofile, frame, 1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, ".data");

    fprintf(fPtr,"Reading v_r: %s\n", full_file);
    fflush(fPtr);

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(vel_r_unprc, sizeof(float)*R_DIM*THETA_DIM*PHI_DIM,R_DIM*THETA_DIM*PHI_DIM, hydroPtr);
    fclose(hydroPtr);

    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_r %d: %e\n", i, *(vel_r_unprc+i));
        fflush(fPtr);
    }

    //v_theta
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 3,"-" );
    modifyFlashName(file_num, hydrofile, frame, 1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, ".data");

    fprintf(fPtr,"Reading v_theta: %s\n", full_file);
    fflush(fPtr);

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(vel_theta_unprc, sizeof(float)*R_DIM*THETA_DIM*PHI_DIM,R_DIM*THETA_DIM*PHI_DIM, hydroPtr);
    fclose(hydroPtr);

    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_theta %d: %e\n", i, *(vel_theta_unprc+i));
        fflush(fPtr);
    }

    //v_phi
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 4,"-" );
    modifyFlashName(file_num, hydrofile, frame, 1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, ".data");

    fprintf(fPtr,"Reading v_phi: %s\n", full_file);
    fflush(fPtr);

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(vel_phi_unprc, sizeof(float)*R_DIM*THETA_DIM*PHI_DIM,R_DIM*THETA_DIM*PHI_DIM, hydroPtr);
    fclose(hydroPtr);

    for (i=0;i<5;i++)
    {
        fprintf(fPtr,"V_phi %d: %e\n", i, *(vel_phi_unprc+i));
        fflush(fPtr);
    }

    //pressure (divided by c^2)
    snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"u0", 8,"-" );
    modifyFlashName(file_num, hydrofile, frame, 1);
    snprintf(full_file, sizeof(full_file), "%s%s", file_num, ".data");

    fprintf(fPtr,"Reading pres: %s\n", full_file);
    fflush(fPtr);

    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    fread(pres_unprc, sizeof(float)*R_DIM*THETA_DIM*PHI_DIM,R_DIM*THETA_DIM*PHI_DIM, hydroPtr);
    fclose(hydroPtr);

    for (i=PHI_DIM-1;i<PHI_DIM;i++)
    {
        for (j=THETA_DIM-1;j<THETA_DIM;j++)
        {
            for (k=R_DIM-5;k<R_DIM;k++)
            {

            fprintf(fPtr,"Pres %d: %e\n", (i*R_DIM*THETA_DIM + j*R_DIM + k  ), *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )));
            fflush(fPtr);

            }
        }
    }

    // see how many elements there are to test if reading correctly
    /*
    hydroPtr=fopen(full_file, "rb");
    fread(&buffer, sizeof(float), 1,hydroPtr); //random stuff about the file from fortran
    while (1 == fread(&buffer,sizeof(float),1,hydroPtr))
    {
        elem++;
    }
    //fread(pres_unprc, sizeof(double)*HYDRO_DIM,HYDRO_DIM, hydroPtr);
    fclose(hydroPtr);

    fprintf(fPtr,"Elem %d\n", elem);
    */

    //R
    if (frame<=1300)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 0,"-x1.data" );
    }
    else if (frame<=2000)
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 1,"-x1.data" );
    }
    else
    {
        snprintf(hydrofile,sizeof(hydrofile),"%s%s%d%s",hydro_prefix,"grid0", 2,"-x1.data" );
    }

    fprintf(fPtr,"Reading Radius: %s\n", hydrofile);
    fflush(fPtr);

    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<R_DIM)
    {
        fscanf(hydroPtr, "%lf", (r_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma

        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(r_unprc+i));
            fflush(fPtr);
        }

        i++;
    }

    fclose(hydroPtr);

    //calculate radial grid edges
    *(r_edge+0)=r_in;
    i=0;
    for (i=1;i<R_DIM;i++)
    {
        *(r_edge+i)=(*(r_edge+i-1))+((*(r_edge+i-1))*(M_PI/560)/(1+((*(r_edge+i-1))/r_ref))); //r_i = r_(i-1) + Dq r_(i-1) [1 + r_(i-1)/r0]-1
        *(dr+i-1)=(*(r_edge+i))-(*(r_edge+i-1));

        if (i<5)
        {
            fprintf(fPtr,"R Edge: %d: %e Dr: %e\n", i, *(r_edge+i), *(dr+i-1));
            fflush(fPtr);
        }

    }
    free(r_edge);

    //Theta
    snprintf(hydrofile,sizeof(hydrofile),"%s%s",hydro_prefix,"grid-x2.data" );
    fprintf(fPtr,"Reading Theta: %s\n", hydrofile);
    fflush(fPtr);

    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<THETA_DIM)
    {
        fscanf(hydroPtr, "%lf", (theta_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma

        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(theta_unprc+i));
            fflush(fPtr);
        }

        i++;
    }

    fclose(hydroPtr);

    //Phi
    snprintf(hydrofile,sizeof(hydrofile),"%s%s",hydro_prefix,"grid-x3.data" );
    fprintf(fPtr,"Reading Phi: %s\n", hydrofile);
    fflush(fPtr);

    hydroPtr=fopen(hydrofile, "r");

    i=0;
    while (i<PHI_DIM)
    {
        fscanf(hydroPtr, "%lf", (phi_unprc+i));  //read value
        fgets(buf, 3,hydroPtr); //read comma

        if (i<5)
        {
            fprintf(fPtr,"R %d: %e\n", i, *(phi_unprc+i));
            fflush(fPtr);
        }

        i++;
    }

    fclose(hydroPtr);

    //limit number of array elements
    elem=0;
    for (i=0;i<PHI_DIM;i++)
    {
        for (j=0;j<THETA_DIM;j++)
        {
            for (k=0;k<R_DIM;k++)
            {
                //if I have photons do selection differently than if injecting photons
                if (ph_inj==0)
                {
                    //if calling this function when propagating photons, choose blocks based on where the photons are
                    if (((ph_rmin - C_LIGHT/fps)<(*(r_unprc+k))) && (*(r_unprc+k)  < (ph_rmax + C_LIGHT/fps) ))
                    {
                        // *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )
                        elem++;
                    }
                }
                else
                {
                    //if calling this function to inject photons choose blocks based on injection parameters, r_inj, which is sufficient
                    if (((r_inj - C_LIGHT/fps)<(*(r_unprc+k))) && (*(r_unprc+k)  < (r_inj + C_LIGHT/fps) ))
                    {
                        // *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )
                        elem++;
                    }

                }

            }
        }
    }

    fprintf(fPtr,"Number of post restricted Elems: %d %e\n", elem, r_inj);

    //allocate space for new set of data
    (*pres)=malloc (elem * sizeof (double ));
    (*velx)=malloc (elem * sizeof (double ));
    (*vely)=malloc (elem * sizeof (double ));
    (*velz)=malloc (elem * sizeof (double ));
    (*dens)=malloc (elem * sizeof (double ));
    (*x)=malloc (elem * sizeof (double ));
    (*y)=malloc (elem * sizeof (double ));
    (*z)=malloc (elem * sizeof (double ));
    (*r)=malloc (elem * sizeof (double ));
    (*theta)=malloc (elem * sizeof (double ));
    (*phi)=malloc (elem * sizeof (double ));
    (*gamma)=malloc (elem * sizeof (double ));
    (*dens_lab)=malloc (elem * sizeof (double ));
    (*szx)=malloc (elem * sizeof (double )); //theta and phi resolution
    (*szy)=malloc (elem * sizeof (double )); //r resolution
    (*temp)=malloc (elem * sizeof (double ));

    //limit number of array elements
    elem=0;
    for (i=0;i<PHI_DIM;i++)
    {
        for (j=0;j<THETA_DIM;j++)
        {
            for (k=0;k<R_DIM;k++)
            {
                //if I have photons do selection differently than if injecting photons
                if (ph_inj==0)
                {
                    //if calling this function when propagating photons, choose blocks based on where the photons are
                    if (((ph_rmin - C_LIGHT/fps)<(*(r_unprc+k))) && (*(r_unprc+k)  < (ph_rmax + C_LIGHT/fps) ))
                    {
                        (*pres)[elem] = *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  ));
                        (*dens)[elem] = *(dens_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  ));
                        (*temp)[elem] =  pow(3*(*(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
                        //fprintf(fPtr, "Selected Values: pres, %e, temp, %e\n",*(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )),pow(3*(*(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0) );

                        (*gamma)[elem] = pow(pow(1.0-(pow(*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+ pow(*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+pow(*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)),0.5),-1);
                        (*dens_lab)[elem] = (*(dens_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )))*pow(pow(1.0-(pow(*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+ pow(*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+pow(*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)),0.5),-1);
                        (*r)[elem] = *(r_unprc+k);
                        (*theta)[elem] = *(theta_unprc+j);
                        (*phi)[elem] = *(phi_unprc+i);
                        (*x)[elem] = (*(r_unprc+k))*sin(*(theta_unprc+j))*cos(*(phi_unprc+i));
                        (*y)[elem] = (*(r_unprc+k))*sin(*(theta_unprc+j))*sin(*(phi_unprc+i));
                        (*z)[elem] = (*(r_unprc+k))*cos(*(theta_unprc+j));
                        (*szx)[elem] = M_PI/560;
                        (*szy)[elem] =  *(dr+k);
                        (*velx)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j))*cos(*(phi_unprc+i))) + ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))*cos(*(phi_unprc+i))) - ((*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(phi_unprc+i)));
                        (*vely)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j))*sin(*(phi_unprc+i))) + ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))*sin(*(phi_unprc+i))) + ((*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(phi_unprc+i)));
                        (*velz)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))) - ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j)));

                        elem++;

                    }
                }
                else
                {
                    //if calling this function to inject photons choose blocks based on injection parameters, r_inj, which is sufficient
                    if (((r_inj - C_LIGHT/fps)<(*(r_unprc+k))) && (*(r_unprc+k)  < (r_inj + C_LIGHT/fps) ))
                    {
                        (*pres)[elem] = *(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  ));
                        (*dens)[elem] = *(dens_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  ));
                        (*temp)[elem] =  pow(3*(*(pres_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
                        (*gamma)[elem] = pow(pow(1.0-(pow(*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+ pow(*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+pow(*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)),0.5),-1);
                        (*dens_lab)[elem] = (*(dens_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k  )))*pow(pow(1.0-(pow(*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+ pow(*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)+pow(*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )),2)),0.5),-1);
                        (*r)[elem] = *(r_unprc+k);
                        (*theta)[elem] = *(theta_unprc+j);
                        (*phi)[elem] = *(phi_unprc+i);
                        (*x)[elem] = (*(r_unprc+k))*sin(*(theta_unprc+j))*cos(*(phi_unprc+i));
                        (*y)[elem] = (*(r_unprc+k))*sin(*(theta_unprc+j))*sin(*(phi_unprc+i));
                        (*z)[elem] = (*(r_unprc+k))*cos(*(theta_unprc+j));
                        (*szx)[elem] = M_PI/560;
                        (*szy)[elem] =  *(dr+k);
                        (*velx)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j))*cos(*(phi_unprc+i))) + ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))*cos(*(phi_unprc+i))) - ((*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(phi_unprc+i)));
                        (*vely)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j))*sin(*(phi_unprc+i))) + ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))*sin(*(phi_unprc+i))) + ((*(vel_phi_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(phi_unprc+i)));
                        (*velz)[elem]=((*(vel_r_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*cos(*(theta_unprc+j))) - ((*(vel_theta_unprc+(i*R_DIM*THETA_DIM + j*R_DIM + k )))*sin(*(theta_unprc+j)));

                        elem++;
                    }

                }

            }
        }
    }

    *number=elem;

    free(pres_unprc); free(dens_unprc); free(r_unprc); free(theta_unprc); free(phi_unprc);free(dr);free(vel_r_unprc); free(vel_theta_unprc); free(vel_phi_unprc);

}

void photonInjection3D( struct photon **ph, int *ph_num, double r_inj, double ph_weight, int min_photons, int max_photons, char spect, int array_length, double fps, double theta_min, double theta_max,\
 double *x, double *y, double *z, double *szx, double *szy, double *r, double *theta, double *phi, double *temps, double *vx, double *vy, double *vz, gsl_rng * rand)
{
    int i=0, block_cnt=0, *ph_dens=NULL, ph_tot=0, j=0,k=0;
    double ph_dens_calc=0.0, fr_dum=0.0, y_dum=0.0, yfr_dum=0.0, fr_max=0, bb_norm=0, position_phi, ph_weight_adjusted, theta_prime=0;
    double com_v_phi, com_v_theta, *p_comv=NULL, *boost=NULL; //comoving phi, theta, comoving 4 momentum for a photon, and boost for photon(to go to lab frame)
    double *l_boost=NULL; //pointer to hold array of lorentz boost, to lab frame, values
    float num_dens_coeff;

    if (spect=='w') //from MCRAT paper, w for wien spectrum
    {
        num_dens_coeff=8.44;
        printf("in wien spectrum\n");
    }
    else
    {
        num_dens_coeff=20.29; //this is for black body spectrum
        printf("in BB spectrum");
    }

    //find how many blocks are near the injection radius within the angles defined in mc.par, get temperatures and calculate number of photons to allocate memory for
    //and then rcord which blocks have to have "x" amount of photons injected there
    printf("%e, %e\n",*(phi+i), theta_max);
    for(i=0;i<array_length;i++)
    {
        //look at all boxes in width delta r=c/fps and within angles we are interested in NEED TO modify for RIKEN data- dont need r anymore, just theta and phi? (didnt work), just look at pojection on x-z plane
        theta_prime=acos(*(y+i)/(*(r+i))); //jet axis here is the y axis
            if ( (theta_prime< theta_max) && (theta_prime > theta_min) ) //(*(r+i) > (r_inj - C_LIGHT/fps))  &&   (*(r+i)  < (r_inj + C_LIGHT/fps)  ) &&
            {
                //printf("%e\n", theta_prime );
                block_cnt++;
            }
    }
    printf("Blocks: %d\n", block_cnt);

    ph_dens=malloc(block_cnt * sizeof(int));

    //calculate the photon density for each block and save it to the array
    j=0;
    ph_tot=0;
    ph_weight_adjusted=ph_weight;
    //printf("%d %d\n", max_photons, min_photons);
    while ((ph_tot>max_photons) || (ph_tot<min_photons) )
    {
        j=0;
        ph_tot=0;
        //allocate memory to record density of photons for each block
        //ph_dens=malloc(block_cnt * sizeof(int));

        for (i=0;i<array_length;i++)
        {
            //printf("%d\n",i);
            //printf("%e, %e, %e, %e, %e, %e\n", *(r+i),(r_inj - C_LIGHT/fps), (r_inj + C_LIGHT/fps), *(theta+i) , theta_max, theta_min);
                //NEED TO modify for RIKEN data - modified
                theta_prime=acos(*(y+i)/(*(r+i)));
                if ( (theta_prime< theta_max) && (theta_prime > theta_min) )
                {
                    //NEED TO modify for RIKEN data - modified
                    ph_dens_calc=(num_dens_coeff*pow(*(temps+i),3.0)*pow(*(r+i),2)*sin(*(theta+i))* pow(*(szx+i),2.0)*(*(szy+i)) /(ph_weight_adjusted))*pow(pow(1.0-(pow(*(vx+i),2)+pow(*(vy+i),2)),0.5),-1) ; //a*T^3/(weight) dV, dV=2*PI*x*dx^2,

                     (*(ph_dens+j))=gsl_ran_poisson(rand,ph_dens_calc) ; //choose from poission distribution with mean of ph_dens_calc

                    //printf("%d, %lf \n",*(ph_dens+j), ph_dens_calc);

                     //sum up all the densities to get total number of photons
                     ph_tot+=(*(ph_dens+j));

                     j++;
                }
        }

        if (ph_tot>max_photons)
        {
            //if the number of photons is too big make ph_weight larger
            ph_weight_adjusted*=10;
            //free(ph_dens);
        }
        else if (ph_tot<min_photons)
        {
            ph_weight_adjusted*=0.5;
            //free(ph_dens);
        }
        //printf("dens: %d, photons: %d\n", *(ph_dens+(j-1)), ph_tot);

    }

    printf("%d\n", ph_tot);

    //allocate memory for that many photons and also allocate memory to hold comoving 4 momentum of each photon and the velocity of the fluid
    (*ph)=malloc (ph_tot * sizeof (struct photon ));

    p_comv=malloc(4*sizeof(double));
    boost=malloc(3*sizeof(double));
    l_boost=malloc(4*sizeof(double));

    //go through blocks and assign random energies/locations to proper number of photons
    ph_tot=0;
    k=0;
    for (i=0;i<array_length;i++)
    {
        theta_prime=acos(*(y+i)/(*(r+i)));
        if ( (theta_prime< theta_max) && (theta_prime > theta_min) )  //NEED TO modify for RIKEN data - modified
        {

            //*(temps+i)=0.76*(*(temps+i));
            for(j=0;j<( *(ph_dens+k) ); j++ )
            {
                    //have to get random frequency for the photon comoving frequency
                    y_dum=1; //initalize loop
                    yfr_dum=0;
                    while (y_dum>yfr_dum)
                    {
                        fr_dum=gsl_rng_uniform_pos(rand)*6.3e11*(*(temps+i)); //in Hz
                        //printf("%lf, %lf ",gsl_rng_uniform_pos(rand), (*(temps+i)));
                        y_dum=gsl_rng_uniform_pos(rand);
                        //printf("%lf ",fr_dum);

                        if (spect=='w')
                        {
                            yfr_dum=(1.0/(1.29e31))*pow((fr_dum/(*(temps+i))),3.0)/(exp((PL_CONST*fr_dum)/(K_B*(*(temps+i)) ))-1); //curve is normalized to maximum
                        }
                        else
                        {
                            fr_max=(5.88e10)*(*(temps+i));//(C_LIGHT*(*(temps+i)))/(0.29); //max frequency of bb
                            bb_norm=(PL_CONST*fr_max * pow((fr_max/C_LIGHT),2.0))/(exp(PL_CONST*fr_max/(K_B*(*(temps+i))))-1); //find value of bb at fr_max
                            yfr_dum=((1.0/bb_norm)*PL_CONST*fr_dum * pow((fr_dum/C_LIGHT),2.0))/(exp(PL_CONST*fr_dum/(K_B*(*(temps+i))))-1); //curve is normalized to vaue of bb @ max frequency
                        }
                        //printf("%lf, %lf,%lf,%e \n",(*(temps+i)),fr_dum, y_dum, yfr_dum);

                    }
                   //printf("%lf\n ",fr_dum);
                   //position_phi= gsl_rng_uniform(rand)*2*M_PI; //NEED TO modify for RIKEN data-modified, dont need anymore
                   com_v_phi=gsl_rng_uniform(rand)*2*M_PI;
                   com_v_theta=acos((gsl_rng_uniform(rand)*2)-1);
                   //printf("%lf, %lf, %lf\n", position_phi, com_v_phi, com_v_theta);

                   //populate 4 momentum comoving array
                   *(p_comv+0)=PL_CONST*fr_dum/C_LIGHT;
                   *(p_comv+1)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*cos(com_v_phi);
                   *(p_comv+2)=(PL_CONST*fr_dum/C_LIGHT)*sin(com_v_theta)*sin(com_v_phi);
                   *(p_comv+3)=(PL_CONST*fr_dum/C_LIGHT)*cos(com_v_theta);

                    //populate boost matrix, not sure why multiplying by -1, seems to give correct answer in old python code...
                    //NEED TO modify for RIKEN data - modified
                    *(boost+0)=-1*(*(vx+i));
                    *(boost+1)=-1*(*(vy+i));
                    *(boost+2)=-1*(*(vz+i));
                    //printf("%lf, %lf, %lf\n", *(boost+0), *(boost+1), *(boost+2));

                    //boost to lab frame
                    lorentzBoost(boost, p_comv, l_boost, 'p');
                    //printf("Assignemnt: %e, %e, %e, %e\n", *(l_boost+0), *(l_boost+1), *(l_boost+2),*(l_boost+3));

                (*ph)[ph_tot].p0=(*(l_boost+0));
                (*ph)[ph_tot].p1=(*(l_boost+1));
                (*ph)[ph_tot].p2=(*(l_boost+2));
                (*ph)[ph_tot].p3=(*(l_boost+3));
                //NEED TO modify for RIKEN data-modified
                (*ph)[ph_tot].r0= (*(x+i)); //put photons @ center of box that they are supposed to be in with random phi
                (*ph)[ph_tot].r1=(*(y+i)) ;
                (*ph)[ph_tot].r2=(*(z+i)); //y coordinate in flash becomes z coordinate in MCRaT
                (*ph)[ph_tot].num_scatt=0;
                (*ph)[ph_tot].weight=ph_weight_adjusted;
                //printf("%d\n",ph_tot);
                ph_tot++;
            }
            k++;
        }
    }

    *ph_num=ph_tot; //save number of photons
    //printf(" %d: %d\n", *(ph_dens+(k-1)), *ph_num);
    free(ph_dens); free(p_comv);free(boost); free(l_boost);

}

 void phMinMax(struct photon *ph, int ph_num, double *min, double *max)
 {
      double temp_r_max=0, temp_r_min=-1;
        int i=0;
      double ph_r=0;

      for (i=0;i<ph_num;i++)
    {
        ph_r=pow(pow( ((ph+i)->r0), 2.0) + pow(((ph+i)->r1),2.0 ) + pow(((ph+i)->r2) , 2.0),0.5);
        if (ph_r > temp_r_max )
        {
            temp_r_max=ph_r;
            printf("The new max is: %e\n", temp_r_max);
        }

        if ((i==0) || (ph_r<temp_r_min))
        {
            temp_r_min=ph_r;
            printf("The new min is: %e\n", temp_r_min);
        }
    }

    *max=temp_r_max;
    *min=temp_r_min;

 }

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


text:
#include "gemm.h"
#include "utils.h"
#include "cuda.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>

#include <omp.h>

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

#ifdef OPENGEMM
#include "gemm_kernels.h"
#endif

void gemm_bin(int M, int N, int K, real_t ALPHA, char *A, int lda, real_t *B,
		int ldb, real_t *C, int ldc) {
	int i, j, k;
	for (i = 0; i < M; ++i) {
		for (k = 0; k < K; ++k) {
			char A_PART = A[i * lda + k];
			if (A_PART) {
				for (j = 0; j < N; ++j) {
					C[i * ldc + j] += B[k * ldb + j];
				}
			} else {
				for (j = 0; j < N; ++j) {
					C[i * ldc + j] -= B[k * ldb + j];
				}
			}
		}
	}
}

real_t *random_matrix(int rows, int cols) {
	int i;
	real_t *m = calloc(rows * cols, sizeof(real_t));
	for (i = 0; i < rows * cols; ++i) {
		m[i] = (real_t) rand() / RAND_MAX;
	}
	return m;
}

void time_random_matrix(int TA, int TB, int m, int k, int n) {
	real_t *a;
	if (!TA)
		a = random_matrix(m, k);
	else
		a = random_matrix(k, m);
	int lda = (!TA) ? k : m;
	real_t *b;
	if (!TB)
		b = random_matrix(k, n);
	else
		b = random_matrix(n, k);
	int ldb = (!TB) ? n : k;

	real_t *c = random_matrix(m, n);
	int i;
	clock_t start = clock(), end;
	for (i = 0; i < 10; ++i) {
		gemm_cpu(TA, TB, m, n, k, 1, a, lda, b, ldb, 1, c, n);
	}
	end = clock();
	printf("Matrix Multiplication %dx%d * %dx%d, TA=%d, TB=%d: %lf ms\n", m, k,
			k, n, TA, TB, (real_t)(end - start) / CLOCKS_PER_SEC);
	free(a);
	free(b);
	free(c);
}

void gemm(int TA, int TB, int M, int N, int K, real_t ALPHA, real_t *A, int lda,
		real_t *B, int ldb, real_t BETA, real_t *C, int ldc) {
	gemm_cpu(TA, TB, M, N, K, ALPHA, A, lda, B, ldb, BETA, C, ldc);
}

void gemm_nn(int M, int N, int K, real_t ALPHA, real_t *A, int lda, real_t *B,
		int ldb, real_t *C, int ldc) {
	int i, j, k;
#pragma omp parallel for
	for (i = 0; i < M; ++i) {
		for (k = 0; k < K; ++k) {
			register real_t A_PART = ALPHA * A[i * lda + k];
			for (j = 0; j < N; ++j) {
				C[i * ldc + j] += A_PART * B[k * ldb + j];
			}
		}
	}
}

void gemm_nt(int M, int N, int K, real_t ALPHA, real_t *A, int lda, real_t *B,
		int ldb, real_t *C, int ldc) {
	int i, j, k;
#pragma omp parallel for
	for (i = 0; i < M; ++i) {
		for (j = 0; j < N; ++j) {
			register real_t sum = 0;
			for (k = 0; k < K; ++k) {
				sum += ALPHA * A[i * lda + k] * B[j * ldb + k];
			}
			C[i * ldc + j] += sum;
		}
	}
}

void gemm_tn(int M, int N, int K, real_t ALPHA, real_t *A, int lda, real_t *B,
		int ldb, real_t *C, int ldc) {
	int i, j, k;
#pragma omp parallel for
	for (i = 0; i < M; ++i) {
		for (k = 0; k < K; ++k) {
			register real_t A_PART = ALPHA * A[k * lda + i];
			for (j = 0; j < N; ++j) {
				C[i * ldc + j] += A_PART * B[k * ldb + j];
			}
		}
	}
}

void gemm_tt(int M, int N, int K, real_t ALPHA, real_t *A, int lda, real_t *B,
		int ldb, real_t *C, int ldc) {
	int i, j, k;
#pragma omp parallel for
	for (i = 0; i < M; ++i) {
		for (j = 0; j < N; ++j) {
			register real_t sum = 0;
			for (k = 0; k < K; ++k) {
				sum += ALPHA * A[i + k * lda] * B[k + j * ldb];
			}
			C[i * ldc + j] += sum;
		}
	}
}

void gemm_cpu(int TA, int TB, int M, int N, int K, real_t ALPHA, real_t *A,
		int lda, real_t *B, int ldb, real_t BETA, real_t *C, int ldc) {
	//printf("cpu: %d %d %d %d %d %f %d %d %f %d\n",TA, TB, M, N, K, ALPHA, lda, ldb, BETA, ldc);

#ifdef OPENBLAS
	/*
	void cblas_sgemm 	( 	const CBLAS_LAYOUT  	layout,
		const CBLAS_TRANSPOSE  	TransA,
		const CBLAS_TRANSPOSE  	TransB,
		const int  	M,
		const int  	N,
		const int  	K,
		const float  	alpha,
		const float *  	A,
		const int  	lda,
		const float *  	B,
		const int  	ldb,
		const float  	beta,
		float *  	C,
		const int  	ldc
	)
	 */
	CBLAS_TRANSPOSE transa, transb;

	if (!TA && !TB){
		transa = CblasNoTrans;
		transb = CblasNoTrans;
	}else if (TA && !TB){
		transa = CblasTrans;
		transb = CblasNoTrans;
	}else if (!TA && TB){
		transa = CblasNoTrans;
		transb = CblasTrans;
	}else{
		transa = CblasTrans;
		transb = CblasTrans;
	}
	//	              0          1       2     3  4  5    6    7   8   9   10  11     12  13
	cblas_sgemm(CblasRowMajor, transa, transb, M, N, K, ALPHA, A, lda, B, ldb, BETA,  C, ldc);
#else

	int i, j;
	for (i = 0; i < M; ++i) {
		for (j = 0; j < N; ++j) {
			C[i * ldc + j] *= BETA;
		}
	}

	if (!TA && !TB)
		gemm_nn(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
	else if (TA && !TB)
		gemm_tn(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
	else if (!TA && TB)
		gemm_nt(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
	else
		gemm_tt(M, N, K, ALPHA, A, lda, B, ldb, C, ldc);
#endif
}

#ifdef GPU

void gemm_gpu(int TA, int TB, int M, int N, int K, real_t ALPHA, real_t *A_gpu,
		int lda, real_t *B_gpu, int ldb, real_t BETA, real_t *C_gpu, int ldc,
		unsigned char use_tensor_cores, cudaStream_t st) {
	cublasHandle_t handle = blas_handle(use_tensor_cores);
	cublasSetStream(handle, st);
//	if(TB || TA)
//	{
//		printf("Matrix need to be transposed %d %d\n", TB, TA);
//	}
#ifndef OPENGEMM

#if REAL_TYPE == HALF
	//run_cuda_gemm_half(int TA, int TB, int M, int N, int K, real_t ALPHA, real_t *A_gpu,
//	int lda, real_t *B_gpu, int ldb, real_t BETA, real_t *C_gpu, int ldc)
	run_cuda_gemm_half(handle, TA, TB, M, N, K, ALPHA, A_gpu, lda, B_gpu, ldb, BETA, C_gpu, ldc, st);
#elif REAL_TYPE == FLOAT
	cudaError_t status = (cudaError_t) cublasSgemm(handle, (TB ? CUBLAS_OP_T : CUBLAS_OP_N),
			(TA ? CUBLAS_OP_T : CUBLAS_OP_N), N, M, K, &ALPHA, B_gpu, ldb,
			A_gpu, lda, &BETA, C_gpu, ldc);
	check_error(status);
#elif REAL_TYPE == DOUBLE
	cudaError_t status = (cudaError_t) cublasDgemm(handle, (TB ? CUBLAS_OP_T : CUBLAS_OP_N),
			(TA ? CUBLAS_OP_T : CUBLAS_OP_N), N, M, K, &ALPHA, B_gpu, ldb,
			A_gpu, lda, &BETA, C_gpu, ldc);
	check_error(status);
#endif

#else

#if REAL_TYPE == HALF
	run_cuda_gemm_half(handle, TA, TB, M, N, K, ALPHA, A_gpu, lda, B_gpu, ldb, BETA, C_gpu, ldc, st);
#elif REAL_TYPE == FLOAT
	sgemm((TB ? CUBLAS_OP_T : CUBLAS_OP_N),
			(TA ? CUBLAS_OP_T : CUBLAS_OP_N), N, M, K, &ALPHA, B_gpu, ldb,
			A_gpu, lda, &BETA, C_gpu, ldc);

#elif REAL_TYPE == DOUBLE
	dgemm((TB ? CUBLAS_OP_T : CUBLAS_OP_N),
			(TA ? CUBLAS_OP_T : CUBLAS_OP_N), N, M, K, &ALPHA, B_gpu, ldb,
			A_gpu, lda, &BETA, C_gpu, ldc);

#endif

#endif

//	cublasHandle_t handle = blas_handle();
//	cudaError_t status = cublasSgemm(handle, (TB ? CUBLAS_OP_T : CUBLAS_OP_N),
//			(TA ? CUBLAS_OP_T : CUBLAS_OP_N), N, M, K, &ALPHA, B_gpu, ldb,
//			A_gpu, lda, &BETA, C_gpu, ldc);

}

void time_gpu_random_matrix(int TA, int TB, int m, int k, int n) {
	real_t *a;
	if (!TA)
		a = random_matrix(m, k);
	else
		a = random_matrix(k, m);
	int lda = (!TA) ? k : m;
	real_t *b;
	if (!TB)
		b = random_matrix(k, n);
	else
		b = random_matrix(n, k);
	int ldb = (!TB) ? n : k;

	real_t *c = random_matrix(m, n);
	int i;
	clock_t start = clock(), end;
	for (i = 0; i < 32; ++i) {
		gemm_gpu(TA, TB, m, n, k, 1, a, lda, b, ldb, 1, c, n, 0, 0x0);
	}
	end = clock();
	printf("Matrix Multiplication %dx%d * %dx%d, TA=%d, TB=%d: %lf s\n", m, k,
			k, n, TA, TB, (real_t)(end - start) / CLOCKS_PER_SEC);
	free(a);
	free(b);
	free(c);
}

void time_gpu(int TA, int TB, int m, int k, int n) {
	int iter = 10;
	real_t *a = random_matrix(m, k);
	real_t *b = random_matrix(k, n);

	int lda = (!TA) ? k : m;
	int ldb = (!TB) ? n : k;

	real_t *c = random_matrix(m, n);

	real_t *a_cl = cuda_make_array(a, m * k);
	real_t *b_cl = cuda_make_array(b, k * n);
	real_t *c_cl = cuda_make_array(c, m * n);

	int i;
	clock_t start = clock(), end;
	for (i = 0; i < iter; ++i) {
		gemm_gpu(TA, TB, m, n, k, 1, a_cl, lda, b_cl, ldb, 1, c_cl, n, 0, 0x0);
		cudaDeviceSynchronize();
	}
	double flop = ((double) m) * n * (2. * k + 2.) * iter;
	double gflop = flop / pow(10., 9);
	end = clock();
	double seconds = sec(end - start);
	printf(
			"Matrix Multiplication %dx%d * %dx%d, TA=%d, TB=%d: %lf s, %lf GFLOPS\n",
			m, k, k, n, TA, TB, seconds, gflop / seconds);
	cuda_free(a_cl);
	cuda_free(b_cl);
	cuda_free(c_cl);
	free(a);
	free(b);
	free(c);
}

void test_gpu_accuracy(int TA, int TB, int m, int k, int n) {
	srand(0);
	real_t *a;
	if (!TA)
		a = random_matrix(m, k);
	else
		a = random_matrix(k, m);
	int lda = (!TA) ? k : m;
	real_t *b;
	if (!TB)
		b = random_matrix(k, n);
	else
		b = random_matrix(n, k);
	int ldb = (!TB) ? n : k;

	real_t *c = random_matrix(m, n);
	real_t *c_gpu = random_matrix(m, n);
	memset(c, 0, m * n * sizeof(real_t));
	memset(c_gpu, 0, m * n * sizeof(real_t));
	int i;
	//pm(m,k,b);
	gemm_gpu(TA, TB, m, n, k, 1, a, lda, b, ldb, 1, c_gpu, n, 0, 0x0);
	//printf("GPU\n");
	//pm(m, n, c_gpu);

	gemm_cpu(TA, TB, m, n, k, 1, a, lda, b, ldb, 1, c, n);
	//printf("\n\nCPU\n");
	//pm(m, n, c);
	double sse = 0;
	for (i = 0; i < m * n; ++i) {
		//printf("%f %f\n", c[i], c_gpu[i]);
		sse += pow(c[i] - c_gpu[i], 2);
	}
	printf("Matrix Multiplication %dx%d * %dx%d, TA=%d, TB=%d: %g SSE\n", m, k,
			k, n, TA, TB, sse / (m * n));
	free(a);
	free(b);
	free(c);
	free(c_gpu);
}

int test_gpu_blas() {
	/*
	 test_gpu_accuracy(0,0,10,576,75);

	 test_gpu_accuracy(0,0,17,10,10);
	 test_gpu_accuracy(1,0,17,10,10);
	 test_gpu_accuracy(0,1,17,10,10);
	 test_gpu_accuracy(1,1,17,10,10);

	 test_gpu_accuracy(0,0,1000,10,100);
	 test_gpu_accuracy(1,0,1000,10,100);
	 test_gpu_accuracy(0,1,1000,10,100);
	 test_gpu_accuracy(1,1,1000,10,100);

	 test_gpu_accuracy(0,0,10,10,10);

	 time_gpu(0,0,64,2916,363);(cudaError_t) cublasDgemm
	 time_gpu(0,0,64,2916,363);
	 time_gpu(0,0,64,2916,363);
	 time_gpu(0,0,192,729,1600);
	 time_gpu(0,0,384,196,1728);
	 time_gpu(0,0,256,196,3456);
	 time_gpu(0,0,256,196,2304);
	 time_gpu(0,0,128,4096,12544);
	 time_gpu(0,0,128,4096,4096);
	 */
	time_gpu(0, 0, 64, 75, 12544);
	time_gpu(0, 0, 64, 75, 12544);
	time_gpu(0, 0, 64, 75, 12544);
	time_gpu(0, 0, 64, 576, 12544);
	time_gpu(0, 0, 256, 2304, 784);
	time_gpu(1, 1, 2304, 256, 784);
	time_gpu(0, 0, 512, 4608, 196);
	time_gpu(1, 1, 4608, 512, 196);

	return 0;
}
#endif

max_stars_repo_name: divadnauj-GB/nvbitfi

max_issues_repo_name: divadnauj-GB/nvbitfi

max_forks_repo_name: divadnauj-GB/nvbitfi