text:
// Copyright (C) 2015 Vincent Lejeune
// For conditions of distribution and use, see copyright notice in License.txt
#pragma once

#include <vector>
#include <tuple>
#include <array>
#include <memory>
#include <gsl/gsl>
#include <variant>
#include <optional>
#include "..\Core\SColor.h"

namespace irr
{
	namespace video
	{
		enum class E_INDEX_TYPE
		{
			EIT_16BIT,
			EIT_32BIT
		};

		enum class E_POLYGON_MODE
		{
			EPM_FILL,
			EPM_LINE,
			EPM_POINT
		};

		enum class E_CULL_MODE
		{
			ECM_NONE,
			ECM_FRONT,
			ECM_BACK,
			ECM_FRONT_AND_BACK
		};

		enum class E_FRONT_FACE
		{
			EFF_CCW,
			EFF_CW
		};

		enum class E_SAMPLE_COUNT
		{
			ESC_1,
			ESC_2,
			ESC_4,
			ESC_8
		};

		enum class E_PRIMITIVE_TYPE
		{
			EPT_POINTS,
			EPT_LINE_STRIP,
			EPT_LINES,
			EPT_TRIANGLE_STRIP,
			EPT_TRIANGLES,
		};

		enum class E_COMPARE_FUNCTION
		{
			ECF_LESS,
			ECF_LEQUAL,
			ECF_NEVER,
		};

		enum class E_STENCIL_OP
		{
			ESO_KEEP,
		};

		enum class E_ASPECT
		{
			EA_COLOR,
			EA_DEPTH,
			EA_STENCIL,
			EA_DEPTH_STENCIL
		};

		enum class E_MEMORY_POOL
		{
			EMP_CPU_WRITEABLE,
			EMP_GPU_LOCAL,
			EMP_CPU_READABLE,
		};

		enum class E_TEXTURE_TYPE
		{
			ETT_2D,
			ETT_CUBE,
		};
	}
}

enum class RESOURCE_USAGE
{
	PRESENT,
	COPY_DEST,
	COPY_SRC,
	RENDER_TARGET,
	READ_GENERIC,
	DEPTH_WRITE,
	undefined,
	uav,
};

enum image_flags
{
	usage_transfer_src = 0x1,
	usage_transfer_dst = 0x2,
	usage_sampled = 0x4,
	usage_render_target = 0x8,
	usage_depth_stencil = 0x10,
	usage_input_attachment = 0x20,
	usage_cube = 0x40,
};

enum buffer_flags
{
	none = 0,
	usage_uav = 0x1,
	usage_texel_buffer = 0x2,
    usage_buffer_transfer_src = 0x4,
	usage_uniform = 0x8,
	usage_index = 0x10,
	usage_vertex = 0x20,
};

enum class SAMPLER_TYPE
{
	NEAREST,
	BILINEAR_CLAMPED,
	TRILINEAR,
	ANISOTROPIC,
};

struct MipLevelData
{
	uint64_t Offset;
	uint32_t Width;
	uint32_t Height;
	uint32_t RowPitch;
};

enum class RESOURCE_VIEW
{
	CONSTANTS_BUFFER,
	INPUT_ATTACHMENT,
	SHADER_RESOURCE,
	TEXEL_BUFFER,
	SAMPLER,
	UAV_BUFFER,
	UAV_IMAGE,
};

enum class shader_stage
{
	vertex_shader,
	fragment_shader,
	all,
};

struct range_of_descriptors
{
	const RESOURCE_VIEW range_type;
	const uint32_t bind_point;
	const uint32_t count;

	constexpr range_of_descriptors(const RESOURCE_VIEW rt, const uint32_t bindpoint, const uint32_t range_size)
		: range_type(rt), bind_point(bindpoint), count(range_size)
	{}
};

struct descriptor_set
{
	const std::vector<range_of_descriptors> descriptors_ranges;
	const shader_stage stage;

	descriptor_set(const std::initializer_list<range_of_descriptors>& ranges, const shader_stage& st)
		: descriptors_ranges(ranges), stage(st)
	{ }
};

struct pipeline_vertex_attributes
{
	uint32_t location;
	irr::video::ECOLOR_FORMAT format;
	uint32_t binding;
	uint32_t stride;
	uint32_t offset;
};
/*	VkPipelineColorBlendAttachmentState blend_attachment_state{ true, VK_BLEND_FACTOR_ONE, VK_BLEND_FACTOR_ONE , VK_BLEND_OP_ADD, VK_BLEND_FACTOR_ONE , VK_BLEND_FACTOR_ONE , VK_BLEND_OP_ADD, -1 };
VkPipelineColorBlendStateCreateInfo blend_state{ VK_STRUCTURE_TYPE_PIPELINE_COLOR_BLEND_STATE_CREATE_INFO, nullptr, 0, false, VK_LOGIC_OP_NO_OP, 1, &blend_attachment_state };*/

enum class blend_factor
{
	zero,
	one,
};

enum class blend_op
{
	add
};

struct color_output
{
	bool blend;
	blend_op op;
	blend_factor src_color;
	blend_factor src_alpha;
	blend_factor dst_color;
	blend_factor dst_alpha;
};

struct compute_pipeline_state_description
{
	std::vector<uint32_t> compute_binary;

	compute_pipeline_state_description set_compute_shader(gsl::span<const uint32_t> code)
	{
		std::copy(code.begin(), code.end(), std::back_inserter(compute_binary));
		return *this;
	}
};

struct graphic_pipeline_state_description
{
	std::vector<uint32_t> vertex_binary;
	std::vector<uint32_t> fragment_binary;

	std::vector<pipeline_vertex_attributes> attributes;
	std::vector<color_output> color_outputs;

	bool rasterization_depth_clamp_enable;
	bool rasterization_discard_enable;
	irr::video::E_POLYGON_MODE rasterization_polygon_mode;
	irr::video::E_CULL_MODE rasterization_cull_mode;
	irr::video::E_FRONT_FACE rasterization_front_face;
	bool rasterization_depth_bias_enable;
	float rasterization_depth_bias_constant_factor;
	float rasterization_depth_bias_clamp;
	float rasterization_depth_bias_slope_factor;
	float rasterization_line_width;
	bool rasterization_conservative_enable;

	bool multisample_multisample_enable;
	irr::video::E_SAMPLE_COUNT multisample_sample_count;
	float multisample_min_sample_shading;
	// sample mask ?
	bool multisample_alpha_to_coverage;
	bool multisample_alpha_to_one;

	irr::video::E_PRIMITIVE_TYPE input_assembly_topology;
	bool input_assembly_primitive_restart;

	bool depth_stencil_depth_test;
	bool depth_stencil_depth_write;
	irr::video::E_COMPARE_FUNCTION depth_stencil_depth_compare_op;
	bool depth_stencil_depth_clip_enable;
	bool depth_stencil_stencil_test;
	irr::video::E_STENCIL_OP depth_stencil_front_stencil_fail_op;
	irr::video::E_STENCIL_OP depth_stencil_front_stencil_depth_fail_op;
	irr::video::E_STENCIL_OP depth_stencil_front_stencil_pass_op;
	irr::video::E_COMPARE_FUNCTION depth_stencil_front_stencil_compare_op;
	irr::video::E_STENCIL_OP depth_stencil_back_stencil_fail_op;
	irr::video::E_STENCIL_OP depth_stencil_back_stencil_depth_fail_op;
	irr::video::E_STENCIL_OP depth_stencil_back_stencil_pass_op;
	irr::video::E_COMPARE_FUNCTION depth_stencil_back_stencil_compare_op;
	float depth_stencil_min_depth_clip;
	float depth_stencil_max_depth_clip;

	static graphic_pipeline_state_description get()
	{
		return graphic_pipeline_state_description(false, false, irr::video::E_POLYGON_MODE::EPM_FILL, irr::video::E_CULL_MODE::ECM_BACK, irr::video::E_FRONT_FACE::EFF_CW, false, 0.f, 0.f, 0.f, 1.f, false,
			false, irr::video::E_SAMPLE_COUNT::ESC_1, 0.f, false, false, irr::video::E_PRIMITIVE_TYPE::EPT_TRIANGLES, false, true, true, irr::video::E_COMPARE_FUNCTION::ECF_LESS, false, false,
			irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_COMPARE_FUNCTION::ECF_NEVER,
			irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_STENCIL_OP::ESO_KEEP, irr::video::E_COMPARE_FUNCTION::ECF_NEVER,
			0.f, 1.f);
	}

	graphic_pipeline_state_description set_depth_compare_function(irr::video::E_COMPARE_FUNCTION depth_compare) const
	{
		return graphic_pipeline_state_description(rasterization_depth_clamp_enable,
			rasterization_discard_enable,
			rasterization_polygon_mode,
			rasterization_cull_mode,
			rasterization_front_face,
			rasterization_depth_bias_enable,
			rasterization_depth_bias_constant_factor,
			rasterization_depth_bias_clamp,
			rasterization_depth_bias_slope_factor,
			rasterization_line_width,
			rasterization_conservative_enable,
			multisample_multisample_enable,
			multisample_sample_count,
			multisample_min_sample_shading,
			multisample_alpha_to_coverage,
			multisample_alpha_to_one,
			input_assembly_topology,
			input_assembly_primitive_restart,
			depth_stencil_depth_test,
			depth_stencil_depth_write,
			depth_compare,
			depth_stencil_depth_clip_enable,
			depth_stencil_stencil_test,
			depth_stencil_front_stencil_fail_op,
			depth_stencil_front_stencil_depth_fail_op,
			depth_stencil_front_stencil_pass_op,
			depth_stencil_front_stencil_compare_op,
			depth_stencil_back_stencil_fail_op,
			depth_stencil_back_stencil_depth_fail_op,
			depth_stencil_back_stencil_pass_op,
			depth_stencil_back_stencil_compare_op,
			depth_stencil_min_depth_clip,
			depth_stencil_max_depth_clip);
	}

	graphic_pipeline_state_description set_depth_write(bool depthwrite) const
	{
		return graphic_pipeline_state_description(rasterization_depth_clamp_enable,
			rasterization_discard_enable,
			rasterization_polygon_mode,
			rasterization_cull_mode,
			rasterization_front_face,
			rasterization_depth_bias_enable,
			rasterization_depth_bias_constant_factor,
			rasterization_depth_bias_clamp,
			rasterization_depth_bias_slope_factor,
			rasterization_line_width,
			rasterization_conservative_enable,
			multisample_multisample_enable,
			multisample_sample_count,
			multisample_min_sample_shading,
			multisample_alpha_to_coverage,
			multisample_alpha_to_one,
			input_assembly_topology,
			input_assembly_primitive_restart,
			depth_stencil_depth_test,
			depthwrite,
			depth_stencil_depth_compare_op,
			depth_stencil_depth_clip_enable,
			depth_stencil_stencil_test,
			depth_stencil_front_stencil_fail_op,
			depth_stencil_front_stencil_depth_fail_op,
			depth_stencil_front_stencil_pass_op,
			depth_stencil_front_stencil_compare_op,
			depth_stencil_back_stencil_fail_op,
			depth_stencil_back_stencil_depth_fail_op,
			depth_stencil_back_stencil_pass_op,
			depth_stencil_back_stencil_compare_op,
			depth_stencil_min_depth_clip,
			depth_stencil_max_depth_clip);
	}

	graphic_pipeline_state_description set_depth_test(bool depth_test) const
	{
		return graphic_pipeline_state_description(rasterization_depth_clamp_enable,
			rasterization_discard_enable,
			rasterization_polygon_mode,
			rasterization_cull_mode,
			rasterization_front_face,
			rasterization_depth_bias_enable,
			rasterization_depth_bias_constant_factor,
			rasterization_depth_bias_clamp,
			rasterization_depth_bias_slope_factor,
			rasterization_line_width,
			rasterization_conservative_enable,
			multisample_multisample_enable,
			multisample_sample_count,
			multisample_min_sample_shading,
			multisample_alpha_to_coverage,
			multisample_alpha_to_one,
			input_assembly_topology,
			input_assembly_primitive_restart,
			depth_test,
			depth_stencil_depth_write,
			depth_stencil_depth_compare_op,
			depth_stencil_depth_clip_enable,
			depth_stencil_stencil_test,
			depth_stencil_front_stencil_fail_op,
			depth_stencil_front_stencil_depth_fail_op,
			depth_stencil_front_stencil_pass_op,
			depth_stencil_front_stencil_compare_op,
			depth_stencil_back_stencil_fail_op,
			depth_stencil_back_stencil_depth_fail_op,
			depth_stencil_back_stencil_pass_op,
			depth_stencil_back_stencil_compare_op,
			depth_stencil_min_depth_clip,
			depth_stencil_max_depth_clip);
	}

	graphic_pipeline_state_description set_vertex_shader(gsl::span<const uint32_t> binary)
	{
		std::copy(binary.begin(), binary.end(), std::back_inserter(vertex_binary));
		return *this;
	}

	graphic_pipeline_state_description set_fragment_shader(gsl::span<const uint32_t> binary)
	{
		std::copy(binary.begin(), binary.end(), std::back_inserter(fragment_binary));
		return *this;
	}

	graphic_pipeline_state_description set_vertex_attributes(const gsl::span<pipeline_vertex_attributes> &attributes_)
	{
		attributes = std::vector<pipeline_vertex_attributes>(attributes_.begin(), attributes_.end());
		return *this;
	}

	graphic_pipeline_state_description set_color_outputs(const gsl::span<color_output> &color_outputs_)
	{
		color_outputs = std::vector<color_output>(color_outputs_.begin(), color_outputs_.end());
		return *this;
	}

	graphic_pipeline_state_description()
	{

	}

private:
	graphic_pipeline_state_description(
		bool depth_clamp_enable,
		bool rasterizer_discard_enable,
		irr::video::E_POLYGON_MODE polygon_mode,
		irr::video::E_CULL_MODE cull_mode,
		irr::video::E_FRONT_FACE front_face,
		bool depth_bias_enable,
		float depth_bias_constant_factor,
		float depth_bias_clamp,
		float depth_bias_slope_factor,
		float line_width,
		bool conservative_enable,
		bool multisample_enable,
		irr::video::E_SAMPLE_COUNT sample_count,
		float min_sample_shading,
		bool alpha_to_coverage,
		bool alpha_to_one,
		irr::video::E_PRIMITIVE_TYPE topology,
		bool primitive_restart,
		bool depth_test,
		bool depth_write,
		irr::video::E_COMPARE_FUNCTION depth_compare_op,
		bool depth_clip_enable,
		bool stencil_test,
		irr::video::E_STENCIL_OP front_stencil_fail_op,
		irr::video::E_STENCIL_OP front_stencil_depth_fail_op,
		irr::video::E_STENCIL_OP front_stencil_pass_op,
		irr::video::E_COMPARE_FUNCTION front_stencil_compare_op,
		irr::video::E_STENCIL_OP back_stencil_fail_op,
		irr::video::E_STENCIL_OP back_stencil_depth_fail_op,
		irr::video::E_STENCIL_OP back_stencil_pass_op,
		irr::video::E_COMPARE_FUNCTION back_stencil_compare_op,
		float min_depth_clip,
		float max_depth_clip
	)
		: rasterization_depth_clamp_enable(depth_clamp_enable),
		rasterization_discard_enable(rasterizer_discard_enable),
		rasterization_polygon_mode(polygon_mode),
		rasterization_cull_mode(cull_mode),
		rasterization_front_face(front_face),
		rasterization_depth_bias_enable(depth_bias_enable),
		rasterization_depth_bias_constant_factor(depth_bias_constant_factor),
		rasterization_depth_bias_clamp(depth_bias_clamp),
		rasterization_depth_bias_slope_factor(depth_bias_slope_factor),
		rasterization_line_width(line_width),
		rasterization_conservative_enable(conservative_enable),
		multisample_multisample_enable(multisample_enable),
		multisample_sample_count(sample_count),
		multisample_min_sample_shading(min_sample_shading),
		multisample_alpha_to_coverage(alpha_to_coverage),
		multisample_alpha_to_one(alpha_to_one),
		input_assembly_topology(topology),
		input_assembly_primitive_restart(primitive_restart),
		depth_stencil_depth_test(depth_test),
		depth_stencil_depth_write(depth_write),
		depth_stencil_depth_compare_op(depth_compare_op),
		depth_stencil_depth_clip_enable(depth_clip_enable),
		depth_stencil_stencil_test(stencil_test),
		depth_stencil_front_stencil_fail_op(front_stencil_fail_op),
		depth_stencil_front_stencil_depth_fail_op(front_stencil_depth_fail_op),
		depth_stencil_front_stencil_pass_op(front_stencil_pass_op),
		depth_stencil_front_stencil_compare_op(front_stencil_compare_op),
		depth_stencil_back_stencil_fail_op(back_stencil_fail_op),
		depth_stencil_back_stencil_depth_fail_op(back_stencil_depth_fail_op),
		depth_stencil_back_stencil_pass_op(back_stencil_pass_op),
		depth_stencil_back_stencil_compare_op(back_stencil_compare_op),
		depth_stencil_min_depth_clip(min_depth_clip),
		depth_stencil_max_depth_clip(max_depth_clip)
	{ }
};

struct framebuffer_t {
	virtual ~framebuffer_t() {}
};

struct pipeline_state_t {
	virtual ~pipeline_state_t() {}
};

struct compute_pipeline_state_t {
	virtual ~compute_pipeline_state_t() {}
};

struct pipeline_layout_t {
	virtual ~pipeline_layout_t() {}
};

struct render_pass_t {
	virtual ~render_pass_t() {}
};

using clear_value_t = std::variant<std::array<float, 4>, std::tuple<float, uint8_t> >;

struct image_view_t {
	virtual ~image_view_t() {}
};

struct sampler_t {
	virtual ~sampler_t() {}
};

struct buffer_view_t {
	virtual ~buffer_view_t() {}
};

struct allocated_descriptor_set {
};

struct descriptor_set_layout {
	virtual ~descriptor_set_layout() {}
};

struct buffer_t {
	virtual void* map_buffer() = 0;
	virtual void unmap_buffer() = 0;

	virtual ~buffer_t() {}
};

struct image_t {
	virtual ~image_t() {}
};

struct descriptor_storage_t {
	virtual std::unique_ptr<allocated_descriptor_set> allocate_descriptor_set_from_cbv_srv_uav_heap(uint32_t starting_index, const std::vector<descriptor_set_layout*> layouts, uint32_t descriptors_count) = 0;
	virtual std::unique_ptr<allocated_descriptor_set> allocate_descriptor_set_from_sampler_heap(uint32_t starting_index, const std::vector<descriptor_set_layout*> layouts, uint32_t descriptors_count) = 0;
	virtual ~descriptor_storage_t() {};
};

struct command_list_t {
	virtual void bind_graphic_descriptor(uint32_t bindpoint, const allocated_descriptor_set& descriptor_set, pipeline_layout_t& sig) = 0;
	virtual void bind_compute_descriptor(uint32_t bindpoint, const allocated_descriptor_set& descriptor_set, pipeline_layout_t& sig) = 0;
	virtual void copy_buffer_to_image_subresource(image_t& destination_image, uint32_t destination_subresource, buffer_t& source, uint64_t offset_in_buffer,
		uint32_t width, uint32_t height, uint32_t row_pitch, irr::video::ECOLOR_FORMAT format) = 0;
	virtual void set_pipeline_barrier(image_t& resource, RESOURCE_USAGE before, RESOURCE_USAGE after, uint32_t subresource, irr::video::E_ASPECT) = 0;
	virtual void set_uav_flush(image_t& resource) = 0;

	virtual void set_viewport(float x, float width, float y, float height, float min_depth, float max_depth) = 0;
	virtual void set_scissor(uint32_t left, uint32_t right, uint32_t top, uint32_t bottom) = 0;
	virtual void set_graphic_pipeline(pipeline_state_t& pipeline) = 0;
	virtual void set_graphic_pipeline_layout(pipeline_layout_t& sig) = 0;
	virtual void set_compute_pipeline(compute_pipeline_state_t& pipeline) = 0;
	virtual void set_compute_pipeline_layout(pipeline_layout_t& sig) = 0;
	virtual void set_descriptor_storage_referenced(descriptor_storage_t& main_heap, descriptor_storage_t* sampler_heap = nullptr) = 0;
	virtual void bind_index_buffer(buffer_t& buffer, uint64_t offset, uint32_t size, irr::video::E_INDEX_TYPE type) = 0;
	virtual void bind_vertex_buffers(uint32_t first_bind, const std::vector<std::tuple<buffer_t&, uint64_t, uint32_t, uint32_t> > &buffer_offset_stride_size) = 0;
	virtual void draw_indexed(uint32_t index_count, uint32_t instance_count, uint32_t base_index, int32_t base_vertex, uint32_t base_instance) = 0;
	virtual void draw_non_indexed(uint32_t vertex_count, uint32_t instance_count, int32_t base_vertex, uint32_t base_instance) = 0;
	virtual void dispatch(uint32_t x, uint32_t y, uint32_t z) = 0;
	virtual void copy_buffer(buffer_t& src, uint64_t src_offset, buffer_t& dst, uint64_t dst_offset, uint64_t size) = 0;

	virtual void clear_depth_stencil(image_t &img, float depth) = 0;
	virtual void clear_depth_stencil(image_t &img, uint8_t stencil) = 0;
	virtual void clear_depth_stencil(image_t &img, float depth, uint8_t stencil) = 0;
	virtual void clear_color(image_t &img, const std::array<float, 4> &clear_colors) = 0;

	virtual void begin_renderpass(render_pass_t& rp, framebuffer_t& fbo,
		gsl::span<clear_value_t> clear_values,
		uint32_t width, uint32_t height) = 0;
	virtual void next_subpass() = 0;
	virtual void end_renderpass() = 0;

	virtual void make_command_list_executable() = 0;
	virtual void start_command_list_recording(struct command_list_storage_t& storage) = 0;
};

struct semaphore_t {
	virtual ~semaphore_t() {};
};

struct fence_t {
	virtual ~fence_t() {};
};

struct command_list_storage_t {
	virtual std::unique_ptr<command_list_t> create_command_list() = 0;
	virtual void reset_command_list_storage() = 0;
	virtual ~command_list_storage_t() {}
};

struct command_queue_t {
	virtual void submit_executable_command_list(command_list_t& command_list, semaphore_t* wait_sem) = 0;
	virtual void wait_for_command_queue_idle() = 0;
};

struct swap_chain_t {
	virtual ~swap_chain_t() {}
	virtual uint32_t get_next_backbuffer_id(semaphore_t& semaphore) = 0;
	virtual std::vector<std::unique_ptr<image_t>> get_image_view_from_swap_chain() = 0;
	virtual void present(command_queue_t& cmdqueue, uint32_t backbuffer_index) = 0;
};

struct device_t {
	virtual std::unique_ptr<command_list_storage_t> create_command_storage() = 0;
	virtual std::unique_ptr<buffer_t> create_buffer(size_t size, irr::video::E_MEMORY_POOL memory_pool, uint32_t flags) = 0;
	virtual std::unique_ptr<buffer_view_t> create_buffer_view(buffer_t&, irr::video::ECOLOR_FORMAT, uint64_t offset, uint32_t size) = 0;
	virtual void set_constant_buffer_view(const allocated_descriptor_set& descriptor_set, uint32_t offset_in_set, uint32_t binding_location, buffer_t& buffer, uint32_t buffer_size, uint64_t offset_in_buffer = 0) = 0;
	virtual void set_uniform_texel_buffer_view(const allocated_descriptor_set& descriptor_set, uint32_t offset_in_set, uint32_t binding_location, buffer_view_t& buffer_view) = 0;
	virtual void set_uav_buffer_view(const allocated_descriptor_set& descriptor_set, uint32_t offset_in_set, uint32_t binding_location, buffer_t& buffer, uint64_t offset, uint32_t size) = 0;
	virtual std::unique_ptr<image_t> create_image(irr::video::ECOLOR_FORMAT format, uint32_t width, uint32_t height, uint16_t mipmap, uint32_t layers, uint32_t flags, clear_value_t *clear_value) = 0;
	virtual std::unique_ptr<image_view_t> create_image_view(image_t& img, irr::video::ECOLOR_FORMAT fmt, uint16_t base_mipmap, uint16_t mipmap_count, uint16_t base_layer, uint16_t layer_count, irr::video::E_TEXTURE_TYPE texture_type, irr::video::E_ASPECT aspect = irr::video::E_ASPECT::EA_COLOR) = 0;
	virtual void set_image_view(const allocated_descriptor_set& descriptor_set, uint32_t offset, uint32_t binding_location, image_view_t& img_view) = 0;
	virtual void set_input_attachment(const allocated_descriptor_set& descriptor_set, uint32_t offset, uint32_t binding_location, image_view_t& img_view) = 0;
	virtual void set_uav_image_view(const allocated_descriptor_set& descriptor_set, uint32_t offset, uint32_t binding_location, image_view_t& img_view) = 0;
	virtual void set_sampler(const allocated_descriptor_set& descriptor_set, uint32_t offset, uint32_t binding_location, sampler_t& sampler) = 0;
	virtual std::unique_ptr<sampler_t> create_sampler(SAMPLER_TYPE sampler_type) = 0;
	virtual std::unique_ptr<descriptor_storage_t> create_descriptor_storage(uint32_t num_sets, const std::vector<std::tuple<RESOURCE_VIEW, uint32_t> > &num_descriptors) = 0;
	virtual std::unique_ptr<framebuffer_t> create_frame_buffer(gsl::span<const image_view_t*> render_targets, uint32_t width, uint32_t height, render_pass_t* render_pass) = 0;
	virtual std::unique_ptr<framebuffer_t> create_frame_buffer(gsl::span<const image_view_t*> render_targets, const image_view_t& depth_stencil_texture, uint32_t width, uint32_t height, render_pass_t* render_pass) = 0;
	virtual std::unique_ptr<descriptor_set_layout> get_object_descriptor_set(const descriptor_set &ds) = 0;
	virtual std::unique_ptr<pipeline_state_t> create_graphic_pso(const graphic_pipeline_state_description&, const render_pass_t&, const pipeline_layout_t&, const uint32_t& subpass) = 0;
	virtual std::unique_ptr<compute_pipeline_state_t> create_compute_pso(const compute_pipeline_state_description&, const pipeline_layout_t&) = 0;
	virtual std::unique_ptr<pipeline_layout_t> create_pipeline_layout(gsl::span<const descriptor_set_layout *>) = 0;
	virtual std::unique_ptr<fence_t> create_fence() = 0;
	virtual std::unique_ptr<semaphore_t> create_semaphore() = 0;

	virtual std::unique_ptr<render_pass_t> create_ibl_sky_pass(const irr::video::ECOLOR_FORMAT&) = 0;
	virtual std::unique_ptr<render_pass_t> create_object_sunlight_pass(const irr::video::ECOLOR_FORMAT&) = 0;
	virtual std::unique_ptr<render_pass_t> create_ssao_pass() = 0;
	virtual std::unique_ptr<render_pass_t> create_blit_pass(const irr::video::ECOLOR_FORMAT& color_format) = 0;

	virtual ~device_t() {};
};

clear_value_t get_clear_value(irr::video::ECOLOR_FORMAT format, float depth, uint8_t stencil);
clear_value_t get_clear_value(irr::video::ECOLOR_FORMAT format, const std::array<float,4> &color);

max_stars_repo_name: vlj/YAGF

max_issues_repo_name: vlj/YAGF

max_forks_repo_name: vlj/YAGF


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

#include <config.h>
#include <stdlib.h>

#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_errno.h>

/*
gsl_spmatrix_equal()
  Return 1 if a = b, 0 otherwise
*/

int
gsl_spmatrix_equal(const gsl_spmatrix *a, const gsl_spmatrix *b)
{
  const size_t M = a->size1;
  const size_t N = a->size2;

  if (b->size1 != M || b->size2 != N)
    {
      GSL_ERROR_VAL("matrices must have same dimensions", GSL_EBADLEN, 0);
    }
  else if (a->sptype != b->sptype)
    {
      GSL_ERROR_VAL("trying to compare different sparse matrix types", GSL_EINVAL, 0);
    }
  else
    {
      const size_t nz = a->nz;
      size_t n;

      if (nz != b->nz)
        return 0; /* different number of non-zero elements */

      if (GSL_SPMATRIX_ISTRIPLET(a))
        {
          /*
           * triplet formats could be out of order but identical, so use
           * gsl_spmatrix_get() on b for each aij
           */
          for (n = 0; n < nz; ++n)
            {
              double bij = gsl_spmatrix_get(b, a->i[n], a->p[n]);

              if (a->data[n] != bij)
                return 0;
            }
        }
      else if (GSL_SPMATRIX_ISCCS(a))
        {
          /*
           * for CCS, both matrices should have everything
           * in the same order
           */

          /* check row indices and data */
          for (n = 0; n < nz; ++n)
            {
              if ((a->i[n] != b->i[n]) || (a->data[n] != b->data[n]))
                return 0;
            }

          /* check column pointers */
          for (n = 0; n < a->size2 + 1; ++n)
            {
              if (a->p[n] != b->p[n])
                return 0;
            }
        }
      else if (GSL_SPMATRIX_ISCRS(a))
        {
          /*
           * for CRS, both matrices should have everything
           * in the same order
           */

          /* check column indices and data */
          for (n = 0; n < nz; ++n)
            {
              if ((a->i[n] != b->i[n]) || (a->data[n] != b->data[n]))
                return 0;
            }

          /* check row pointers */
          for (n = 0; n < a->size1 + 1; ++n)
            {
              if (a->p[n] != b->p[n])
                return 0;
            }
        }
      else
        {
          GSL_ERROR_VAL("unknown sparse matrix type", GSL_EINVAL, 0);
        }

      return 1;
    }
} /* gsl_spmatrix_equal() */

ext: c

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#include <scraps/config.h>

#include <scraps/Temp.h>
#include <scraps/type-traits.h>
#include <scraps/utility.h>

SCRAPS_IGNORE_WARNINGS_PUSH
#include <gsl.h>
SCRAPS_IGNORE_WARNINGS_POP

#include <array>

namespace scraps {

/**
 * Convenience functions to create an std::array from a gsl::span
 */
template <typename T, std::ptrdiff_t ArraySize, template <typename, std::ptrdiff_t...> class Span>
constexpr std::array<std::remove_cv_t<T>, ArraySize> ToArray(Span<T, ArraySize> span) {
    static_assert(ArraySize != gsl::dynamic_range, "Dynamic range spans are not allowed.");
    static_assert(ArraySize > 0, "ArraySize must be greater than 0");

    std::array<std::remove_cv_t<T>, ArraySize> ret{};
    std::copy(span.begin(), span.end(), ret.begin());
    return ret;
}

template <typename T, std::ptrdiff_t ArraySize, template <typename, std::ptrdiff_t...> class Span>
constexpr std::array<std::remove_cv_t<T>, ArraySize> ToArray(Temp<Span<T, ArraySize>> span) {
    return ToArray(static_cast<Span<T, ArraySize>>(span));
}

template <typename T, size_t N>
constexpr std::array<std::remove_cv_t<T>, N> ToArray(T (&arr)[N]) {
    std::array<std::remove_cv_t<T>, N> ret{};
    std::copy(arr, arr + N, ret.begin());
    return ret;
}

/**
 * Convenience type for defining a c-style array using the template parameters
 * of a std::array.
 */
template <typename ArrayT>
using CArray = typename RemoveCVRType<ArrayT>::value_type[std::tuple_size<RemoveCVRType<ArrayT>>::value];

/**
 * Hashing struct suitable for using an std::array as the key in stl
 * associative containers such as unordered_map.
 */
struct ArrayHasher {
    template <typename T, size_t N>
    size_t operator()(const std::array<T, N>& arr) const {
        return HashRange(arr.begin(), arr.end());
    }
};

} // namespace scraps

namespace std {
    template<typename T, size_t N>
    struct hash<std::array<T, N>> {
        size_t operator()(const std::array<T, N>& array) const {
            return scraps::ArrayHasher{}(array);
        }
    };
} // namespace std

ext: h

max_stars_repo_name: carlbrown/scraps

max_issues_repo_name: carlbrown/scraps

max_forks_repo_name: carlbrown/scraps


text:
/* multimin/fdfminimizer.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 *
 * 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_multimin.h>

gsl_multimin_fdfminimizer *
gsl_multimin_fdfminimizer_alloc (const gsl_multimin_fdfminimizer_type * T,
				 size_t n)
{
  int status;

  gsl_multimin_fdfminimizer *s =
    (gsl_multimin_fdfminimizer *) malloc (sizeof (gsl_multimin_fdfminimizer));

  if (s == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for minimizer struct",
		     GSL_ENOMEM, 0);
    }

  s->type = T;

  s->x = gsl_vector_calloc (n);

  if (s->x == 0)
    {
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for x", GSL_ENOMEM, 0);
    }

  s->gradient = gsl_vector_calloc (n);

  if (s->gradient == 0)
    {
      gsl_vector_free (s->x);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for gradient", GSL_ENOMEM, 0);
    }

  s->dx = gsl_vector_calloc (n);

  if (s->dx == 0)
    {
      gsl_vector_free (s->x);
      gsl_vector_free (s->gradient);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for dx", GSL_ENOMEM, 0);
    }

  s->state = malloc (T->size);

  if (s->state == 0)
    {
      gsl_vector_free (s->x);
      gsl_vector_free (s->gradient);
      gsl_vector_free (s->dx);
      free (s);
      GSL_ERROR_VAL ("failed to allocate space for minimizer state",
		     GSL_ENOMEM, 0);
    }

  status = (T->alloc) (s->state, n);

  if (status != GSL_SUCCESS)
    {
      free (s->state);
      gsl_vector_free (s->x);
      gsl_vector_free (s->gradient);
      gsl_vector_free (s->dx);
      free (s);

      GSL_ERROR_VAL ("failed to initialize minimizer state", GSL_ENOMEM, 0);
    }

  return s;
}

int
gsl_multimin_fdfminimizer_set (gsl_multimin_fdfminimizer * s,
                               gsl_multimin_function_fdf * fdf,
                               const gsl_vector * x,
                               double step_size, double tol)
{
  if (s->x->size != fdf->n)
    {
      GSL_ERROR ("function incompatible with solver size", GSL_EBADLEN);
    }

  if (x->size != fdf->n)
    {
      GSL_ERROR ("vector length not compatible with function", GSL_EBADLEN);
    }

  s->fdf = fdf;

  gsl_vector_memcpy (s->x,x);
  gsl_vector_set_zero (s->dx);

  return (s->type->set) (s->state, s->fdf, s->x, &(s->f), s->gradient, step_size, tol);
}

void
gsl_multimin_fdfminimizer_free (gsl_multimin_fdfminimizer * s)
{
  (s->type->free) (s->state);
  free (s->state);
  gsl_vector_free (s->dx);
  gsl_vector_free (s->gradient);
  gsl_vector_free (s->x);
  free (s);
}

int
gsl_multimin_fdfminimizer_iterate (gsl_multimin_fdfminimizer * s)
{
  return (s->type->iterate) (s->state, s->fdf, s->x, &(s->f), s->gradient, s->dx);
}

int
gsl_multimin_fdfminimizer_restart (gsl_multimin_fdfminimizer * s)
{
  return (s->type->restart) (s->state);
}

const char *
gsl_multimin_fdfminimizer_name (const gsl_multimin_fdfminimizer * s)
{
  return s->type->name;
}

gsl_vector *
gsl_multimin_fdfminimizer_x (gsl_multimin_fdfminimizer * s)
{
  return s->x;
}

gsl_vector *
gsl_multimin_fdfminimizer_dx (gsl_multimin_fdfminimizer * s)
{
  return s->dx;
}

gsl_vector *
gsl_multimin_fdfminimizer_gradient (gsl_multimin_fdfminimizer * s)
{
  return s->gradient;
}

double
gsl_multimin_fdfminimizer_minimum (gsl_multimin_fdfminimizer * s)
{
  return s->f;
}


text:
#ifndef __EXECUTE_PLUGIN_H__
#define __EXECUTE_PLUGIN_H__

#include <memory>
#include <optional>
#include <stdexcept>
#include <string>
#include <vector>

#include <gsl/gsl>

#include "config/commandLineOptions.h"
#include "config/fleetingOptionsInterface.h"
#include "config/pattern.h"
#include "config/patternsHandler.h"
#include "config/settingsNode.h"

#include "plugin.h"

namespace execHelper::plugins {
using Plugins = std::map<std::string, std::shared_ptr<const Plugin>>;

/**
 * \brief Exception thrown when the requested plugin is invalid
 *
 * Exception thrown when the requested plugin is invalid e.g. due to the fact that it can not be found
 */
struct InvalidPlugin : public std::runtime_error {
  public:
    /**
     * Create an invalid plugin
     *
     * \param[in] msg   A message detailing the specifics of the exception
     */
    inline explicit InvalidPlugin(const std::string& msg)
        : std::runtime_error(msg) {}

    /*! @copydoc InvalidPlugin(const std::string&)
     */
    inline explicit InvalidPlugin(const char* msg) : std::runtime_error(msg) {}
};

/**
 * \brief Plugin for executing arbitrary configured commands and/or plugins
 *
 * The ExecutePlugin handles the context of and calls all plugins. It uses the prototype pattern for retrieving a map of plugins it can call.
 */
class ExecutePlugin : public Plugin {
  public:
    /**
     * Create an executePlugin instance
     *
     * \param[in] commandsToExecute  The commands to execute with this plugin
     * instance
     */
    explicit ExecutePlugin(
        const config::CommandCollection& commandsToExecute) noexcept;

    /**
     * \param[in] commandsToExecute The commands to execute
     * \param[in] initialCommand    The initial command that is being executed
     */
    ExecutePlugin(const config::CommandCollection& commandsToExecute,
                  const config::Command& initialCommand) noexcept;

    config::VariablesMap
    getVariablesMap(const config::FleetingOptionsInterface& fleetingOptions)
        const noexcept override;
    bool apply(core::Task task, const config::VariablesMap& variables,
               const config::Patterns& patterns) const noexcept override;

    std::string summary() const noexcept override;

    /**
     * Returns a list with the names of all known plugins
     *
     * @returns A list of plugin names
     */
    static auto getPluginNames() noexcept -> std::vector<std::string>;

    /**
     * Returns an instance of the plugin associated with the given name
     *
     * \param[in] pluginName    The plugin to get the associated instance from
     * \returns A pointer to the new instance
     * \throws  InvalidPlugin   When no plugin associated with the given pluginName is found
     */
    static std::shared_ptr<const Plugin>
    getPlugin(const std::string& pluginName);

    /**
     * Push the given fleeting options on the stack
     *
     * \param[in] fleetingOptions   The fleeting options to use. The last option
     * on the stack will be used for calling the commands. \returns True    if
     * the options were successfully pushed False   otherwise
     */
    static bool push(gsl::not_null<const config::FleetingOptionsInterface*>
                         fleetingOptions) noexcept;

    /**
     * Push the given settings node on the stack
     *
     * \param[in] settings   The settings node to use. The last settings node on
     * the stack will be used for calling the commands. \returns True    if the
     * settings node was successfully pushed False   otherwise
     */
    static bool push(config::SettingsNode&& settings) noexcept;

    /**
     * Push the given patterns on the stack
     *
     * \param[in] patterns   The patterns to use. The last patterns on the stack
     * will be used for calling the commands.
     * \returns True    if the patterns were successfully pushed
     *          False   otherwise
     */
    static bool push(config::Patterns&& patterns) noexcept;

    /**
     * Push the plugin prototypes to the stack
     *
     * \param[in] plugins   Mapping of discovered plugin prototypes
     */
    static void push(Plugins&& plugins) noexcept;

    /**
     * Pop the last fleeting options from the stack
     *
     */
    static void popFleetingOptions() noexcept;

    /**
     * Pop the last settings node from the stack
     */
    static void popSettingsNode() noexcept;

    /**
     * Pop the last patterns from the stack
     */
    static void popPatterns() noexcept;

    /**
     * Pop the last plugin prototypes from the stack
     */
    static void popPlugins() noexcept;

  private:
    static auto getNextStep(const config::Command& command,
                            const config::Command& originalCommand) noexcept
        -> std::shared_ptr<const Plugin>;
    static bool
    getVariablesMap(config::VariablesMap* variables,
                    const std::vector<config::SettingsKeys>& keys,
                    const config::SettingsNode& rootSettings) noexcept;
    static void index(config::VariablesMap* variables,
                      const config::SettingsNode& settings,
                      const config::SettingsKeys& key) noexcept;

    const config::CommandCollection m_commands;
    const config::CommandCollection m_initialCommands;

    static std::vector<gsl::not_null<const config::FleetingOptionsInterface*>>
        m_fleeting;
    static std::vector<config::SettingsNode> m_settings;
    static std::vector<config::PatternsHandler> m_patterns;
    static std::vector<Plugins> m_plugins;
};
} // namespace execHelper::plugins

#endif /* __EXECUTE_PLUGIN_H__ */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
/* roots/fdfsolver.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Reid Priedhorsky, 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 <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>

gsl_root_fdfsolver *
gsl_root_fdfsolver_alloc (const gsl_root_fdfsolver_type * T)
{

  gsl_root_fdfsolver * s = (gsl_root_fdfsolver *) malloc (sizeof (gsl_root_fdfsolver));

  if (s == 0)
    {
      GSL_ERROR_VAL ("failed to allocate space for root solver struct",
			GSL_ENOMEM, 0);
    };

  s->state = malloc (T->size);

  if (s->state == 0)
    {
      free (s);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for root solver state",
			GSL_ENOMEM, 0);
    };

  s->type = T ;
  s->fdf = NULL;

  return s;
}

int
gsl_root_fdfsolver_set (gsl_root_fdfsolver * s, gsl_function_fdf * f, double root)
{
  s->fdf = f;
  s->root = root;

  return (s->type->set) (s->state, s->fdf, &(s->root));
}

int
gsl_root_fdfsolver_iterate (gsl_root_fdfsolver * s)
{
  return (s->type->iterate) (s->state, s->fdf, &(s->root));
}

void
gsl_root_fdfsolver_free (gsl_root_fdfsolver * s)
{
  free (s->state);
  free (s);
}

const char *
gsl_root_fdfsolver_name (const gsl_root_fdfsolver * s)
{
  return s->type->name;
}

double
gsl_root_fdfsolver_root (const gsl_root_fdfsolver * s)
{
  return s->root;
}


text:
/*
* -----------------------------------------------------------------
*  Pairwise Mixing Stirred Reaction Library --- pmsr_lib.h
*  Version: 1.6180
*  Date: Nov 16, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for the PMSR model.
* -----------------------------------------------------------------
*/

#ifndef __PMSR_LIB_H__
#define __PMSR_LIB_H__

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

/*
* ------------------------------------------------------------
*  Types: struct pmsr, pmsr_wrk
* ------------------------------------------------------------
* This structure contains fields of PMSR model.
*
* last update: Nov 5, 2009
* ------------------------------------------------------------
*/

typedef struct pmsr
{
    int *pmsr_idx;         /* particles index vector */
    gsl_vector **pmsr_ps;  /* particle system        */
} pmsr_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void pmsr_title();

int pmsr_input(unsigned long int *seed,
               unsigned int *Np,
               unsigned int *K,
               double *t0,
               double *delta_t,
               double *tau_res,
               double *tau_mix,
	           double *tau_pair,
               double *Tmax,
               double *Tmin);

pmsr_wrk *pmsr_alloc();

int pmsr_init(int Np,
              int Neq,
              pmsr_wrk *pmsr);

void pmsr_free(int Np,
               void **pmsr_bl);

void pmsr_set_all(int Np,
                  gsl_vector *phi,
                  pmsr_wrk *pmsr);

int pmsr_Nexc(int Np,
              double delta_t,
              double tau_res);

int pmsr_Npair(int Np,
               double delta_t,
               double tau_pair);

void pmsr_meanvar(int Np,
                  int k,
                  gsl_vector **ps,
                  double *mean,
                  double *var);

void pmsr_mixture(int Np,
                  int Neq,
                  double delta_t,
                  double tau_mix,
                  int *idx,
                  gsl_vector **ps);

void pmsr_iops(int Np,
               int Nexc,
               int Npair,
               gsl_rng *r,
               gsl_vector *phi,
               pmsr_wrk *pmsr);

#endif /* __PMSR_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
/**
 *
 * @file example_dpotrf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

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

static void GENMAT_SYM_FULL(int m, double *A)
{
	srand48(time(NULL));

	int j;
	for (j = 0; j < m; ++j ) {
		int i;
		for( i = j; i < m; ++i ) {
			double dran = drand48();
			A[j*m+i] = A[i*m+j] = dran;
		}
  	}
	for(j = 0; j < m; ++j)
		A[j*m+j] += 10 * m;
}

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

    int N     = 4000;
	int LDA   = 4000;
    int info_factorization;

    double *A   = (double *)malloc(LDA*N*sizeof(double));
	double *X	 = (double *)malloc(LDA*N*sizeof(double));
	double *B	 = (double *)malloc(LDA*N*sizeof(double));

    /* Check if unable to allocate memory */
    if ((!A)||(!X)||(!B)){
        fprintf(stderr, "Out of Memory \n ");
        return 0;
    }

	PLASMA_Set( PLASMA_RUNTIME_MODE, PLASMA_OMPSS); //Toggle runtime mdoe to OmpSs
	PLASMA_Set( PLASMA_TILE_SIZE, 384); //Set up tile size

	GENMAT_SYM_FULL(N, A); //Initialize A with a SPD matrix

	LAPACKE_dlarnv(IONE, ISEED, LDA*N, X);

	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N,
			1.0, A, LDA, X, LDA, 0.0, B, LDA);

    PLASMA_dpotrf(PlasmaLower, N, A, LDA);
	PLASMA_dtrsm(PlasmaLeft, PlasmaLower, PlasmaNoTrans, PlasmaNonUnit, N, N, 1.0, A, LDA, B, LDA);
	PLASMA_dtrsm(PlasmaLeft, PlasmaLower, PlasmaTrans, PlasmaNonUnit, N, N, 1.0, A, LDA, B, LDA);

	char infnorm = 'F';
	double Normb = LAPACKE_dlange(LAPACK_COL_MAJOR, infnorm, N, N, B, LDA);
	double Normx = LAPACKE_dlange(LAPACK_COL_MAJOR, infnorm, N, N, X, LDA);
	printf("Residual= %e, norm(b): %e, norm(x): %e\n", fabs(Normb-Normx)/Normx, Normb, Normx);

    free(A); free(X); free(B);

    return 0;
}


text:
#pragma once

#include <cstddef>
#include <iostream>
#include <limits>

#include <cuda/define_specifiers.hpp>
#include <cuda/runtime_api.hpp>

#include <gsl-lite/gsl-lite.hpp>

#include <cub/cub.cuh>

#include <thrustshift/fill.h>
#include <thrustshift/math.h>
#include <thrustshift/not-a-vector.h>

namespace thrustshift {

namespace kernel {

template <typename IndexT>
__global__ void wrap_subgroups(gsl_lite::span<const IndexT> subgroup_ptrs,
                               gsl_lite::span<IndexT> group_ptrs,
                               IndexT mean_group_size) {
	const IndexT gtid = threadIdx.x + blockIdx.x * blockDim.x;
	const IndexT num_rows = subgroup_ptrs.size() - 1;
	const IndexT N = group_ptrs.size() - 1;
	if (gtid < num_rows) {
		const IndexT row_id = gtid;
		const IndexT curr_work = subgroup_ptrs[row_id + 1];
		const IndexT prev_work = subgroup_ptrs[row_id];
		const IndexT l = curr_work / mean_group_size;
		const IndexT m = prev_work / mean_group_size;
		if (prev_work != curr_work && l != 0 && (l != m) &&
		    l != N) { // only first thread is allowed to write to last element
			group_ptrs[l] = row_id + 1;
		}
	}
	if (gtid == 0) {
		group_ptrs[N] = num_rows;
		group_ptrs[0] = 0;
	}
}

} // namespace kernel

namespace {

template <typename StorageIndex>
struct unequal_to_functor {
	StorageIndex x;
	CUDA_FD bool operator()(StorageIndex y) const {
		return y != x;
	}
};

} // namespace

namespace async {

/*! \brief Wrap subsequent subgroups of varying size into groups.
 *
 *  No subgroup can be part of two groups. Each subgroup can
 *  contain a certain amount of elements (also no elements at all). This
 *  function tries to wrap the subgroups into bigger groups of size `mean_group_size`.
 *  This function is especially useful to distribute the work amongst CUDA thread blocks,
 *  which operate on rows of a sparse matrix. If a single subgroup is larger than `mean_group_size`
 *  the subgroup will form a group on its own, which is larger than `mean_group_size`.
 *  The last group is rather larger than `mean_group_size` instead of creating one
 *  group which is smaller than `mean_group_size`.
 *
 *  \param stream CUDA stream.
 *  \param subgroup_ptrs Array of length `num_subgroups + 1`. `subgroup_ptrs[i]` denotes
 *    the index of the first element of subgroup `i`. `subgroup_ptrs[num_subgroups]` must
 *    be equal to `num_elements`.
 *  \param num_elements The sum of the sizes of all subgroups.
 *  \param group_ptrs The result. `group_ptrs[k]` is the start subgroup ID of group `i`, whereas
 *    `group_ptrs[k+1]` is the end subgroup ID of group `i`. Thus group `i` contains all subgroups
 *    in range `[group_ptrs[k], group_ptrs[k+1])`. The last valid output
 *    is `group_ptrs[*group_ptrs_size - 1] = subgroup_ptrs.size() - 1`. The maximum number of groups
 *    is `num_elements / mean_group_size + 1`. The minimum number of groups is 2. Thus the buffer
 *    of `group_ptrs` must be of size `max(num_elements / mean_group_size + 1, 2)`.
 *  \param mean_group_size Intended group size, which is not always possible. E.g. if there
 *    is only one subgroup, there can also be only one group.
 *  \param group_ptrs_size Pointer to an element in GPU memory. After successful execution the
 *    total number of groups is written to this variable.
 */
template <typename IndexT, class MemoryResource>
void wrap_subgroups(cuda::stream_t& stream,
                    gsl_lite::span<const IndexT> subgroup_ptrs,
                    IndexT num_elements,
                    gsl_lite::span<IndexT> group_ptrs,
                    IndexT mean_group_size,
                    gsl_lite::not_null<IndexT*> group_ptrs_size,
                    MemoryResource& delayed_memory_resource) {

	constexpr cuda::grid::block_dimension_t block_dim = 128;

	if (subgroup_ptrs.size() < 2) {
		// There is not a single subgroup
		return;
	}

	if (gsl_lite::narrow<IndexT>(group_ptrs.size()) !=
	    std::max(num_elements / mean_group_size + 1, (IndexT) 2)) {
		std::cerr << "ERROR: group_ptrs buffer has wrong size! actual "
		             "size = "
		          << group_ptrs.size()
		          << ", required size = max(num_elements / "
		             "mean_group_size + 1, 2) = "
		          << std::max(num_elements / mean_group_size + 1, 2)
		          << std::endl;
		std::terminate();
	}

	unequal_to_functor<IndexT> f({std::numeric_limits<IndexT>::max()});

	gsl_lite::span<std::byte> tmp_storage;
	size_t tmp_storage_size = 0;

	auto enqueue_select_if = [&]() {
		cuda::throw_if_error(cub::DeviceSelect::If(tmp_storage.data(),
		                                           tmp_storage_size,
		                                           group_ptrs.data(),
		                                           group_ptrs.data(),
		                                           group_ptrs_size,
		                                           group_ptrs.size(),
		                                           f,
		                                           stream.handle()));
	};
	// set temporary memory size only
	enqueue_select_if();

	auto tmp =
	    make_not_a_vector<std::byte>(tmp_storage_size, delayed_memory_resource);
	tmp_storage = tmp.to_span();

	fill(stream, group_ptrs, std::numeric_limits<IndexT>::max());

	gsl_Expects(subgroup_ptrs.size() > 0);
	const IndexT num_rows = subgroup_ptrs.size() - 1;
	const IndexT grid_dim =
	    ceil_divide(num_rows, gsl_lite::narrow<IndexT>(block_dim));
	cuda::enqueue_launch(kernel::wrap_subgroups<IndexT>,
	                     stream,
	                     cuda::make_launch_config(grid_dim, block_dim),
	                     subgroup_ptrs,
	                     group_ptrs,
	                     mean_group_size);

	enqueue_select_if();
}

} // namespace async

} // namespace thrustshift

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
#pragma once
#include "utility/types.h"
#include "utility/c++std/optional.h"
#include <gsl/gsl>
#include <vector>
#include <string>

namespace cws80 {

enum class FileChooserMode {
    Open,
    Save,
};

struct FileChooserFilter {
    const char *name;
    std::vector<const char *> patterns;
};

class NativeUI {
public:
    virtual ~NativeUI() {}

    static NativeUI *create();

    virtual std::string choose_file(gsl::span<const FileChooserFilter> filters, const std::string &directory, const std::string &title, FileChooserMode mode) = 0;
    virtual cxx::optional<std::string> edit_line(const std::string &title, const std::string &initial_value) = 0;
    virtual uint select_by_menu(gsl::span<const std::string> choices, uint initial_selection) = 0;
};

}  // namespace cws80


text:
//====---- Sudoku/Location.h                                          ----====//
//
// class Location: represent locations within a board,
// and calculate location related properties.
//====--------------------------------------------------------------------====//
//
// Everything is `constexpr` and `noexcept`.
// All values are `const`.
//
// A separate Location_Block class offers block properties.
// There is no implicit conversion between these.
//
//====--------------------------------------------------------------------====//
#pragma once

#include "Size.h"

#include <gsl/gsl> // index

#include <type_traits> // is_signed

namespace Sudoku
{
template<int N>
class Location
{
	using Size  = ::Sudoku::Size<N>;
	using index = gsl::index;

	// prefer signed integers for calculations
	static_assert(std::is_signed_v<decltype(Size::base)>);
	static_assert(std::is_signed_v<decltype(Size::elem)>);

	static constexpr gsl::index location(index row, index col) noexcept
	{
		return row * Size::elem + col; // [0,full)
	}

public:
	// constructors
	constexpr Location() = default;
	explicit constexpr Location(index element) noexcept : id_(element) {}
	constexpr Location(index row, index col) noexcept : id_(location(row, col))
	{
	}

	// information
	[[nodiscard]] constexpr index element() const noexcept
	{
		return id_; // default [0,full)
	}
	[[nodiscard]] constexpr index row() const noexcept
	{
		return id_ / Size::elem; //[0,elem)
	}
	[[nodiscard]] constexpr index col() const noexcept
	{
		return id_ % Size::elem; //[0,elem)
	}
	[[nodiscard]] constexpr index block() const noexcept
	{
		return row() / Size::base * Size::base + col() / Size::base; // [0,elem)
	}
	[[nodiscard]] constexpr index block_row() const noexcept
	{
		return row() % Size::base; // [0,base)
	}
	[[nodiscard]] constexpr index block_col() const noexcept
	{
		return col() % Size::base; // [0,base)
	}
	[[nodiscard]] constexpr index block_elem() const noexcept
	{
		return block_row() * Size::base + block_col(); // [0,elem)
	}

	// comparison
	friend constexpr bool
		operator==(const Location& a, const Location& b) noexcept
	{ // Friend definitions: abseil.io/tips/99
		return a.id_ == b.id_;
	}

private:
	const index id_{};
};

//====--------------------------------------------------------------------====//

template<int N>
class Location_Block
{
	using Size     = ::Sudoku::Size<N>;
	using Location = ::Sudoku::Location<N>;

	// prefer signed integers for calculations
	static_assert(std::is_signed_v<decltype(Size::base)>);

	static constexpr gsl::index
		block_element(gsl::index row, gsl::index col) noexcept
	{
		return row * Size::base + col; // [0,elem)
	}
	static constexpr Location
		block_loc(gsl::index id, gsl::index element) noexcept
	{
		const gsl::index row{
			(id / Size::base) * Size::base + element / Size::base};
		const gsl::index col{
			(id % Size::base) * Size::base + element % Size::base};
		return Location(row, col);
	}
	static constexpr Location
		block_loc(gsl::index id, gsl::index row, gsl::index col) noexcept
	{
		return block_loc(id, block_element(row, col));
	}

public:
	constexpr Location_Block() = default;
	explicit constexpr Location_Block(Location loc) noexcept : loc_(loc)
	{ // empty constructor
	}
	constexpr Location_Block(gsl::index id, gsl::index element) noexcept
		: loc_(block_loc(id, element))
	{ // empty constructor
	}
	constexpr Location_Block(
		gsl::index id, gsl::index row, gsl::index col) noexcept
		: loc_(block_loc(id, row, col))
	{ // empty constructor
	}

	[[nodiscard]] constexpr gsl::index id() const noexcept
	{
		return loc_.block(); // [0,elem)
	}
	[[nodiscard]] constexpr gsl::index element() const noexcept
	{
		return loc_.block_elem(); // [0,elem)
	}
	[[nodiscard]] constexpr gsl::index row() const noexcept
	{
		return loc_.block_row(); // [0,base)
	}
	[[nodiscard]] constexpr gsl::index col() const noexcept
	{
		return loc_.block_col(); // [0,base)
	}

	// [[implicit]]
	// NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
	constexpr operator Location() const noexcept { return loc_; }

private:
	const Location loc_{};
};

//====--------------------------------------------------------------------====//
// free-function declarations
template<int N>
constexpr bool
	operator==(const Location_Block<N>&, const Location_Block<N>&) noexcept;
template<int N>
constexpr bool
	operator==(const Location<N>&, const Location_Block<N>&) noexcept;
template<int N>
constexpr bool
	operator==(const Location_Block<N>&, const Location<N>&) noexcept;

template<int N, typename typeB>
constexpr bool operator!=(const Location<N>&, const typeB&) noexcept;
template<int N, typename typeB>
constexpr bool operator!=(const Location_Block<N>&, const typeB&) noexcept;

template<int N>
constexpr bool operator<(const Location<N>&, const Location<N>&) noexcept;
template<int N>
constexpr bool
	operator<(const Location_Block<N>&, const Location_Block<N>&) noexcept;
template<int N>
constexpr bool operator<=(const Location<N>&, const Location<N>&) noexcept;
template<int N>
constexpr bool
	operator<=(const Location_Block<N>&, const Location_Block<N>&) noexcept;
template<int N>
constexpr bool operator>=(const Location<N>&, const Location<N>&) noexcept;
template<int N>
constexpr bool
	operator>=(const Location_Block<N>&, const Location_Block<N>&) noexcept;
template<int N>
constexpr bool operator>(const Location<N>&, const Location<N>&) noexcept;
template<int N>
constexpr bool
	operator>(const Location_Block<N>&, const Location_Block<N>&) noexcept;

//====--------------------------------------------------------------------====//
// definitions

template<int N>
constexpr bool operator==(
	const Location_Block<N>& left, const Location_Block<N>& right) noexcept
{
	return Location<N>{left} == Location<N>{right};
}

template<int N>
constexpr bool
	operator==(const Location_Block<N>& block, const Location<N>& loc) noexcept
{
	return Location<N>{block} == loc;
}

template<int N>
inline constexpr bool
	operator==(const Location<N>& loc, const Location_Block<N>& block) noexcept
{
	return block == loc;
}

template<int N, typename typeB>
inline constexpr bool
	operator!=(const Location<N>& left, const typeB& right) noexcept
{
	static_assert(
		std::is_same_v<Location<N>, typeB> ||
		std::is_same_v<Location_Block<N>, typeB>);
	return !(left == right);
}

template<int N, typename typeB>
inline constexpr bool
	operator!=(const Location_Block<N>& left, const typeB& right) noexcept
{
	static_assert(
		std::is_same_v<Location<N>, typeB> ||
		std::is_same_v<Location_Block<N>, typeB>);
	return !(left == right);
}

template<int N>
constexpr bool
	operator<(const Location<N>& left, const Location<N>& right) noexcept
{
	return left.element() < right.element();
}

template<int N>
constexpr bool
	operator<=(const Location<N>& left, const Location<N>& right) noexcept
{
	return left.element() <= right.element();
}

template<int N>
constexpr bool
	operator>=(const Location<N>& left, const Location<N>& right) noexcept
{
	return left.element() >= right.element();
}

template<int N>
constexpr bool
	operator>(const Location<N>& left, const Location<N>& right) noexcept
{
	return left.element() > right.element();
}

// Sorted by block first.
template<int N>
constexpr bool operator<(
	const Location_Block<N>& left, const Location_Block<N>& right) noexcept
{
	return (left.id() == right.id()) ? (left.element() < right.element())
									 : left.id() < right.id();
}

template<int N>
constexpr bool operator<=(
	const Location_Block<N>& left, const Location_Block<N>& right) noexcept
{
	return (left.id() == right.id()) ? (left.element() <= right.element())
									 : left.id() < right.id();
}

template<int N>
constexpr bool operator>=(
	const Location_Block<N>& left, const Location_Block<N>& right) noexcept
{
	return (left.id() == right.id()) ? (left.element() >= right.element())
									 : left.id() > right.id();
}

template<int N>
constexpr bool operator>(
	const Location_Block<N>& left, const Location_Block<N>& right) noexcept
{
	return (left.id() == right.id()) ? (left.element() > right.element())
									 : left.id() > right.id();
}

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
/* permutation/gsl_permutation.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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.
 */

#ifndef __GSL_PERMUTATION_H__
#define __GSL_PERMUTATION_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>

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

__BEGIN_DECLS

struct gsl_permutation_struct
{
  size_t size;
  size_t *data;
};

typedef struct gsl_permutation_struct gsl_permutation;

gsl_permutation *gsl_permutation_alloc (const size_t n);
gsl_permutation *gsl_permutation_calloc (const size_t n);
void gsl_permutation_init (gsl_permutation * p);
void gsl_permutation_free (gsl_permutation * p);

int gsl_permutation_fread (FILE * stream, gsl_permutation * p);
int gsl_permutation_fwrite (FILE * stream, const gsl_permutation * p);
int gsl_permutation_fscanf (FILE * stream, gsl_permutation * p);
int gsl_permutation_fprintf (FILE * stream, const gsl_permutation * p, const char *format);

size_t gsl_permutation_size (const gsl_permutation * p);
size_t * gsl_permutation_data (const gsl_permutation * p);

size_t gsl_permutation_get (const gsl_permutation * p, const size_t i);
int gsl_permutation_swap (gsl_permutation * p, const size_t i, const size_t j);

int gsl_permutation_valid (gsl_permutation * p);
void gsl_permutation_reverse (gsl_permutation * p);
int gsl_permutation_inverse (gsl_permutation * inv, const gsl_permutation * p);
int gsl_permutation_next (gsl_permutation * p);
int gsl_permutation_prev (gsl_permutation * p);

extern int gsl_check_range;

#ifdef HAVE_INLINE

extern inline
size_t
gsl_permutation_get (const gsl_permutation * p, const size_t i)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (i >= p->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return p->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_PERMUTATION_H__ */


text:
/*This program evaluates the differential equations for a photon's geodesic in a FRW perturbed spacetime in SPHERICAL coordinates and restricted to photon RADIAL MOTION only.*/

/*This program solves the particular case for flat FRW perturbed spacetime with metric: $g_{ab} = {[g]}_{ab} + h_{ab}$. Where ${[g]}_{ab}$ corresponds to the flat FRW metric and $h_{ab}$ corresponds to the perturbation in the conformal Newtonian gauge. A Plummer potential or a Hernquist potential with adequate parameters can be used to simulate the perturbation.
The equations are written in the form $\frac{d(x or p)^{\alpha}}{d\lambda}=f(x^{\alpha},p^{\alpha})$ and the indice $\alpha$ runs from 0 to 1 since the motion is only radial.
Where $p^{\alpha}={\dot{x}}^{\alpha}$.
The coordinates for the photon's geodesics are then: (ct,r) = (x0,x1).*/

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>          //GSL error management module
#include <gsl/gsl_spline.h>         //GSL interpolation module

/* PHYSICAL CONSTANTS AND PARAMETERS */
#define A 10.0     //Distance parameter of the perturbations
#define G 43007.01     //Gravitational constant
#define M 15000.0     //Mass of the perturbation
#define C 299792.458  //Speed of light
#define GAMMA 0.8    //Parameter for halo mass evolution

/* PROGRAM PARAMETERS */
#define NLINES 100000 //Number of lines in geodesic_solution.dat file
#define NSTEPS 600000000 //Number of steps for solving geodesics
#define NLINESFRW 10000 //Number of lines in frw.dat file
#define DLAMBDA 0.01   //Geodesics parameter step

typedef long double mydbl;

/*Interpolation of function inside *spline object evaluated at an abscisa 'x'.
Argument *spline is a pointer to a spline object which stores the type of interpolation to be made. x is the independent variable where the function is evaluated. *acc is a pointer to a lookup object for interpolations.*/
double interpolator(gsl_spline *spline, double x, gsl_interp_accel *acc)
{
  double a = gsl_spline_eval(spline, x, acc);  //Interpolates data to abcisa x using method in spline and acceleration object acc
  return a;  //Return value of interpolated function at x
}

/*Mass evolution. Function providing the mass at a given z. Mass at present time is obtained at z=0.*/
mydbl mass(mydbl a)
{
  mydbl z = 1.0/a - 1.0;
  return M*expl(-GAMMA*z);
}

/*Plummer model. Provided two variables pointers, this function stores the potential and the derivative of the potential for the Plummer model in these variables at a given radius.*/
void plummer_model(mydbl a, mydbl adot, mydbl r, mydbl *potential, mydbl *der_potential_radial, mydbl *der_potential_temporal)
{
  *potential = -G*mass(a)/(sqrtl(A*A + r*r));
  *der_potential_radial = G*mass(a)*r/(powl(A*A+r*r, 1.5));
  *der_potential_temporal = *potential*GAMMA*adot/(a*a);
}

/*Hernquist model. Provided two variables pointers, this function stores the potential and the derivative of the potential for the Hernquist model in these variables at a given radius.*/
void hernquist_model(mydbl a, mydbl adot, mydbl r, mydbl *potential, mydbl *der_potential_radial, mydbl *der_potential_temporal)
{
  double rsign = copysign(1.0,(double)(1.0*r));
  mydbl absr = fabsl(r);
  if(absr > 2000.0)
    {
      *potential = 0.0;
      *der_potential_radial = 0.0;
      *der_potential_temporal = 0.0;
    }
  else if(rsign == -1.0)
    {
      *potential = -G*mass(a)/(A + absr);
      *der_potential_radial = -G*mass(a)/powl(absr+A,2.0);
      *der_potential_temporal = *potential*GAMMA*adot/(a*a);
    }
  else if(rsign == 1.0)
    {
      *potential = -G*mass(a)/(A + absr);
      *der_potential_radial = G*mass(a)/powl(absr+A,2.0);
      *der_potential_temporal = *potential*GAMMA*adot/(a*a);
    }
}

/*Function of the 0th momentum component differential equation for the geodesics.
${p0}^{dot} = f0(x^{\alpha},p^{\alpha})$.*/
mydbl geodesic_equation_0(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl p0, mydbl pr, mydbl x0, mydbl r, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  double t = (double)(1.0*x0/C);
  mydbl a = (mydbl) 1.0*interpolator(spline1, t, acc1);
  mydbl adot = (mydbl) 1.0*interpolator(spline2, t, acc2);
  mydbl potential, der_potential_radial, der_potential_temporal;
  (*model)(a, adot, r, &potential, &der_potential_radial, &der_potential_temporal);
  mydbl f = -2.0*p0*(der_potential_temporal*p0/C + der_potential_radial*pr)/(C*C + 2.0*potential) + der_potential_temporal*(p0*p0 + a*a*pr*pr)/(C*C*C) - (1.0 - 2.0*potential/(C*C))*(a*adot)*(pr*pr)/(C + 2.0*potential/C);
  return f;
}

/*Function of the 1th (radial) momentum component differential equation for the geodesics.
${p1}^{dot} = f1(x^{\alpha},p^{\alpha})$.*/
mydbl geodesic_equation_r(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl p0, mydbl pr, mydbl x0, mydbl r, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  double t = (double)(1.0*x0/C);
  mydbl a = (mydbl) 1.0*interpolator(spline1, t, acc1);
  mydbl adot = (mydbl) 1.0*interpolator(spline2, t, acc2);
  mydbl potential, der_potential_radial, der_potential_temporal;
  (*model)(a, adot, r, &potential, &der_potential_radial, &der_potential_temporal);
  mydbl f = - (der_potential_radial*p0*p0)/(a*a*(C*C - 2.0*potential)) - (2.0*adot*p0*pr)/(C*a) + der_potential_radial*(pr*pr)/(C*C - 2.0*potential) + 2.0*der_potential_temporal*p0*pr/(a*a*C*(C*C - 2.0*potential));
  return f;
}

/*Function for solving the geodesics differential equations using 4th order Runge-Kutta method.
Arguments are pointer so variables in that memory addresses are changed every time this function is called.*/
void runge_kutta_4(gsl_spline *spline1, gsl_interp_accel *acc1, gsl_spline *spline2, gsl_interp_accel *acc2, mydbl *x0, mydbl *x1, mydbl *p0, mydbl *p1, mydbl *lambda, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  /*Increment in the variables of the differential equation we want to solve*/
  mydbl dx0, dx1, dp0, dp1;

  /*dxi = (k1,j + 2*k2,j + 2*k3,j + k4,j)/6. In this sections the ki,j are declared with i=1,2,3,4.*/
  mydbl k1x0, k1x1, k1p0, k1p1;
  mydbl k2x0, k2x1, k2p0, k2p1;
  mydbl k3x0, k3x1, k3p0, k3p1;
  mydbl k4x0, k4x1, k4p0, k4p1;

  /*This section calculates the k1 quantities*/
  k1x0 = *p0*DLAMBDA; k1x1 = *p1*DLAMBDA;
  k1p0 = DLAMBDA*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0, *p1, *x0, *x1, (*model));
  k1p1 = DLAMBDA*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0, *p1, *x0, *x1, (*model));

  /*This section calculates the k2 quantities*/
  k2x0 = DLAMBDA*(*p0 + 0.5*k1p0); k2x1 = DLAMBDA*(*p1 + 0.5*k1p1);
  k2p0 = DLAMBDA*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + 0.5*k1p0, *p1 + 0.5*k1p1, *x0 + 0.5*k1x0, *x1 + 0.5*k1x1, (*model));
  k2p1 = DLAMBDA*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + 0.5*k1p0, *p1 + 0.5*k1p1, *x0 + 0.5*k1x0, *x1 + 0.5*k1x1, (*model));

  /*This section calculates the k3 quantities*/
  k3x0 = DLAMBDA*(*p0 + 0.5*k2p0); k3x1 = DLAMBDA*(*p1 + 0.5*k2p1);
  k3p0 = DLAMBDA*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + 0.5*k2p0, *p1 + 0.5*k2p1, *x0 + 0.5*k2x0, *x1 + 0.5*k2x1, (*model));
  k3p1 = DLAMBDA*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + 0.5*k2p0, *p1 + 0.5*k2p1, *x0 + 0.5*k2x0, *x1 + 0.5*k2x1, (*model));

  /*This section calculates the k4 quantities*/
  k4x0 = DLAMBDA*(*p0 + k3p0); k4x1 = DLAMBDA*(*p1 + k3p1);
  k4p0 = DLAMBDA*geodesic_equation_0(spline1, acc1, spline2, acc2, *p0 + k3p0, *p1 + k3p1, *x0 + k3x0, *x1 + k3x1, (*model));
  k4p1 = DLAMBDA*geodesic_equation_r(spline1, acc1, spline2, acc2, *p0 + k3p0, *p1 + k3p1, *x0 + k3x0, *x1 + k3x1, (*model));

  /*Calculation of the increments*/
  dx0 = (k1x0 + 2.0*k2x0 + 2.0*k3x0 + k4x0)/6.0;
  dx1 = (k1x1 + 2.0*k2x1 + 2.0*k3x1 + k4x1)/6.0;
  dp0 = (k1p0 + 2.0*k2p0 + 2.0*k3p0 + k4p0)/6.0;
  dp1 = (k1p1 + 2.0*k2p1 + 2.0*k3p1 + k4p1)/6.0;

  /*New values of the variables of the differential equation. Since we are using pointers, when called the routine the value of variable change.*/
  *x0 = *x0 + dx0; *x1 = *x1 + dx1;
  *p0 = *p0 + dp0; *p1 = *p1 + dp1;

  /*Increment of parameter of geodesics*/
  *lambda = *lambda + DLAMBDA;
}

/*To set the initial value of pr, it must hold $g_{\mu\nu}p^{\mu}p^{\nu} = 0$.
This factor multiplies p0 to guarantee that p1 fulfill the null geodesic condition.*/
mydbl condition_factor(mydbl r, double a, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  mydbl potential, der_potential_radial, der_potential_temporal;
  (*model)((mydbl)(1.0*a), 0.0, r, &potential, &der_potential_radial, &der_potential_temporal);
  return (mydbl)(1.0/a)*sqrtl((C*C + 2.0*potential)/(C*C - 2.0*potential));
}

/*$cp^{0}$ multiplied by this factor allows to obtain the energy for a local inertial observer in this spacetime.*/
mydbl energy_factor(mydbl r, double a, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  mydbl potential, der_potential_radial, der_potential_temporal;
  (*model)((mydbl)(1.0*a), 0.0, r, &potential, &der_potential_radial, &der_potential_temporal);
  mydbl g = sqrtl(1.0 + 2.0*potential/(C*C));
  return g;
}

/*Violation of null geodesics condition $g_{\mu\nu}p^{\mu}p^{\nu} = 0$.*/
mydbl violation(mydbl r, mydbl p0, mydbl pr, double a, void (*model)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *))
{
  mydbl potential, der_potential_radial, der_potential_temporal;
  (*model)((mydbl)(1.0*a), 0.0, r, &potential, &der_potential_radial, &der_potential_temporal);
  mydbl f = -(1.0+2.0*potential/(C*C))*p0*p0 + (mydbl)(1.0*a*a)*(1.0-2.0*potential/(C*C))*pr*pr;
  return f;
}

int main(void)
{
  /***READ SCALE FACTOR DATA AND PREPARE OBJECTS FOR INTERPOLATION ***/
  int i;            //For array manipulation

  /*Pointer to scale_factor.data file*/
  FILE *frw;
  frw = fopen("../scale_factor.dat","r");

  /*Variables and arrays to read the data*/
  double cosmictime[NLINESFRW], conftime, scale_factor[NLINESFRW], der_scale_factor[NLINESFRW];

  /*Reading the data*/
  for(i=0; i<NLINESFRW; i++)
    {
      fscanf(frw,"%lf %lf %lf %lf", &cosmictime[i], &conftime, &scale_factor[i], &der_scale_factor[i]);
    }

  /*Free space in memory*/
  fclose(frw);

  /*** Initializes objects for interpolation. 1 is for interpolation of scale factor, 2 is for interpolation of derivative of scale factor ***/

  /*Allocate space in memory*/
  gsl_interp_accel *acc1 = gsl_interp_accel_alloc();  //Acceleration type object (for index lookup)
  gsl_interp_accel *acc2 = gsl_interp_accel_alloc();
  gsl_spline *spline1 = gsl_spline_alloc(gsl_interp_cspline, NLINESFRW);  //Spline type object (define interpolation type and space in memory, works for both)
  gsl_spline *spline2 = gsl_spline_alloc(gsl_interp_cspline, NLINESFRW);

  /*Initializes objects for interpolation*/
  gsl_spline_init(spline1, cosmictime, scale_factor, NLINESFRW);  //Initializes spline object for data cosmictime, scale_factor of size NLINES
  gsl_spline_init(spline2, cosmictime, der_scale_factor, NLINESFRW);  //Initializes spline object for data cosmictime, der_scale_factor of size NLINES

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

  /***SOLVES GEODESIC EQUATIONS FOR PERTURBED FRW UNIVERSE WITH A GIVEN POTENTIAL MODEL ***/

  /*Choose model to use for solving geodesics*/
  void (*model_in_use)(mydbl, mydbl, mydbl, mydbl *, mydbl *, mydbl *);
  model_in_use = hernquist_model;

  /*Initial conditions*/
  mydbl ti = 9.0 ,x0, r = -3000.0, p0 = 1.0e-3, pr, lambda = 0.0, energy1, energy, v, difft, difference, pfrw, isw;
  double difftfrw, aem, aobs;
  x0 = C*ti;
  aem = interpolator(spline1, (double)(1.0*ti), acc1);
  pr = condition_factor(r, aem, model_in_use)*p0;
  energy1 = C*energy_factor(r, aem, model_in_use)*p0;
  v = violation(r, p0, pr, aem, model_in_use);
  difft = (energy1 - energy1)/energy1;
  difftfrw = (aem/aem) - 1.0;
  difference = difft - (mydbl)(1.0*difftfrw);
  pfrw = ((mydbl)(1.0*aem)*1.0e-03/((mydbl)(1.0*aem)));
  isw = p0/pfrw - 1.0;

  /*Pointer to file where solution of differential equation will be saved.*/
  FILE *geodesic;
  geodesic = fopen("geodesic_solution.dat","w");

  /*Write line of initial values in file*/
  fprintf(geodesic, "%16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8e %16.8Le %16.8Le\n", lambda, x0, r, p0, pr, energy1, v, difft, difftfrw, difference, isw);

  long long int ii;

  /*Solution of the differential equation*/
  for(ii=0; ii< NSTEPS; ii++)
    {
      runge_kutta_4(spline1, acc1, spline2, acc2, &x0, &r, &p0, &pr, &lambda, model_in_use);
      if((ii%(NSTEPS/NLINES)) == 0)
	{
	  ti = x0/C;
	  aobs = interpolator(spline1, (double)(1.0*ti), acc1);
	  energy = C*energy_factor(r, aobs, model_in_use)*p0;
	  v = violation(r, p0, pr, aobs, model_in_use);
	  difft = (energy - energy1)/energy1;
	  difftfrw = (aem/aobs) - 1.0;
	  difference = difft - (mydbl)(1.0*difftfrw);
	  pfrw = (mydbl)(1.0*aem)*1.0e-03/((mydbl)(1.0*aobs));
	  isw = p0/pfrw - 1.0;
	  fprintf(geodesic, "%16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8Le %16.8e %16.8Le %16.8Le\n", lambda, x0, r, p0, pr, energy, v, difft, difftfrw, difference, isw);
	}
    }

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

  /*** Releasing all used space in memory ***/
  fclose(geodesic); //Close file storing the results
  gsl_spline_free(spline1);  //Free memory of spline object
  gsl_spline_free(spline2);
  gsl_interp_accel_free(acc1);  //Free memory of accel object
  gsl_interp_accel_free(acc2);
}


text:
#ifndef __REGIONAL_H_
#define __REGIONAL_H_

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "defs.h"

typedef struct stochasticity_grid {
    t_float *arr;
    int height, width;
    double weight_factor;
    double stddev;
    double mean;
    int depth;
    gsl_rng *rng;
    double *weights;
} t_stochasticity_grid;

t_stochasticity_grid *create_stochasticity_grid(t_float *arr, int height, int width, double weight_factor, double mean, double stddev, long int seed);
void free_stochasticity_grid(t_stochasticity_grid *g);
void generate_stochasticity(t_stochasticity_grid *g);

#endif /* __REGIONAL_H_ */

max_stars_repo_name: rybicki/corridor-spom

max_issues_repo_name: rybicki/corridor-spom

max_forks_repo_name: rybicki/corridor-spom


text:
#pragma once
#include <vector>
#include <algorithm>

#include <gsl/gsl.h>

#include "utils/linear_algebra.h"

namespace util{
namespace math{

template<typename T, int64_t dim_1>
struct VecLoop_void{
    template<typename OP, typename... Args>
    void operator()(OP const&fun, Args&&... args) const {
        #pragma clang loop vectorize(enable)
        for(int64_t i=0; i<dim_1; ++i){
            fun(i, std::forward<Args>(args)...);
        }
    }
};
template<typename T, int64_t dim_1>
struct VecLoop_vec{
    template<typename OP, typename... Args>
    auto operator()(OP const&fun, Args&&... args) const {
        util::math::Vector<T,dim_1> result{};
        #pragma clang loop vectorize(enable)
        for(int64_t i=0; i<dim_1; ++i){
            result.span[i]=fun(i, std::forward<Args>(args)...);
        }
        return result;
    }
};

template<typename T, int64_t dim_1, int64_t dim_2>
struct MatLoop_mat{
    template<typename OP, typename... Args>
    void operator()(OP const&fun, Args&&... args) const {
        util::math::Matrix<T,dim_1,dim_2> result{};
        for(int64_t i=0; i<dim_1; ++i){
            for(int64_t j=0; j<dim_2; ++j){
                result.span[i][j]=fun(i,j, std::forward<Args>(args)...);
            }
        }
    }
};
template<typename T, int64_t dim_1, int64_t dim_2>
struct MatLoop_void{
    template<typename OP, typename... Args>
    void operator()(OP const&fun, Args&&... args) const {
        for(int64_t i=0; i<dim_1; ++i){
            for(int64_t j=0; j<dim_2; ++j){
                fun(i,j, std::forward<Args>(args)...);
            }
        }
    }
};

}//namespace util::math
}//namespace util

namespace gsl{

auto add_assign_vec=[](int64_t i, auto &out, auto const & x) {
    out[i]+=x[i];
};
auto sub_assign_vec=[](int64_t i, auto &out, auto const & x) {
    out[i]-=x[i];
};
auto mul_assign_vec=[](int64_t i, auto &out, auto x) {
    out[i]*=x;
};

auto add_assign_mat=[](int64_t i,int64_t j, auto &out, auto const & x) {
    out[i][j]+=x[i][j];
};
auto sub_assign_mat=[](int64_t i,int64_t j, auto &out, auto const & x) {
    out[i][j]-=x[i][j];
};
auto mul_assign_mat=[](int64_t i,int64_t j, auto &out, auto x) {
    out[i][j]*=x;
};

template<typename T, int64_t M>
span<T,M>& operator +=(span<T,M>& out, const span<T,M>& x){
    auto vecloop_void=util::math::VecLoop_void<T,M>{};
    vecloop_void(add_assign_vec, out, x);
    return out;
}
template<typename T, int64_t M>
span<T,M>& operator -=(span<T,M>& out, const span<T,M>& x){
    auto vecloop_void=util::math::VecLoop_void<T,M>{};
    vecloop_void(sub_assign_vec, out, x);
    return out;
}
template<typename T, int64_t M>
span<T,M>& operator *=(span<T,M>& out, T x){
    auto vecloop_void=util::math::VecLoop_void<T,M>{};
    vecloop_void(mul_assign_vec, out, x);
    return out;
}

template<typename T, int64_t M, int64_t N>
span<T,M,N>& operator +=(span<T,M,N>& out, const span<T,M,N>& x){
    auto matloop_void=util::math::MatLoop_void<T,M,N>{};
    matloop_void(add_assign_mat, out, x);
    return out;
}
template<typename T, int64_t M, int64_t N>
span<T,M,N>& operator -=(span<T,M,N>& out, const span<T,M,N>& x){
    auto matloop_void=util::math::MatLoop_void<T,M,N>{};
    matloop_void(sub_assign_mat, out, x);
    return out;
}
template<typename T, int64_t M, int64_t N>
span<T,M,N>& operator *=(span<T,M,N>& out, T x){
    auto matloop_void=util::math::MatLoop_void<T,M,N>{};
    matloop_void(mul_assign_mat, out, x);
    return out;
}

auto grad_update_mat=[](int64_t i,int64_t j,
                        auto &param, auto const &grad, auto scale) {
    param[i][j] += scale*grad[i][j];
};
auto grad_update_vec=[](int64_t i,
                        auto &param, auto const &grad, auto scale) {
    param[i] += scale*grad[i];
};
template<typename T, int64_t M, int64_t N>
void grad_update(span<T,M,N>& param, span<T,M,N> const &grad, T scale){
    auto matloop_void=util::math::MatLoop_void<T,M,N>{};
    matloop_void(grad_update_mat, param, grad, scale);
}
template<typename T, int64_t M>
void grad_update(span<T,M>& param, span<T,M> const &grad, T scale){
    auto vecloop_void=util::math::VecLoop_void<T,M>{};
    vecloop_void(grad_update_vec, param, grad, scale);
}

}//namespace gsl

max_stars_repo_name: uphere-co/nlp-prototype

max_issues_repo_name: uphere-co/nlp-prototype

max_forks_repo_name: uphere-co/nlp-prototype


text:
//
// Created by daeyun on 4/11/17.
//

#pragma once

#ifndef NDEBUG
#define USE_OMP false
#else
#define USE_OMP true
#endif

#include <iostream>
#include <limits>
#include <cstdlib>
#include <vector>
#include <array>
#include <string>
#include <map>
#include <set>
#include <memory>

#include <gsl/gsl_assert>
#include <Eigen/Dense>
#include "spdlog/spdlog.h"
#include "spdlog/sinks/stdout_color_sinks.h"

using Eigen::Matrix;
using Eigen::MatrixXd;
using Eigen::Dynamic;

using Vec = Eigen::VectorXd;
using Vec2 = Eigen::Vector2d;
using Vec2i = Eigen::Matrix<int, 2, 1>;
using Vec3 = Eigen::Vector3d;
using Vec4 = Eigen::Vector4d;
using Mat44 = Eigen::Matrix<double, 4, 4>;
using Mat34 = Eigen::Matrix<double, 3, 4>;
using Mat33 = Eigen::Matrix<double, 3, 3>;
using MatrixX = Eigen::MatrixXd;

using Points4d = Eigen::Matrix<double, 4, Eigen::Dynamic>;
using Points3d = Eigen::Matrix<double, 3, Eigen::Dynamic>;
using Points2d = Eigen::Matrix<double, 2, Eigen::Dynamic>;
using Points2i = Eigen::Matrix<int, 2, Eigen::Dynamic>;
using Points1d = Eigen::Matrix<double, 1, Eigen::Dynamic>;

using std::vector;
using std::array;
using std::string;
using std::unique_ptr;
using std::shared_ptr;
using std::make_unique;
using std::make_shared;
using std::move;
using std::tuple;
using std::pair;
using std::map;
using std::set;
using std::get;

constexpr double kInfinity = std::numeric_limits<double>::infinity();
constexpr double kOneOverTwoPi = M_1_PI * 0.5;
constexpr double kEpsilon = 1e-9;
constexpr double kDeg2Rad = M_PI / 180.0;

// Make sure to have the following in the executable: spdlog::stdout_color_mt("console");
#define LOGGER spdlog::get("console")


text:
/*
CUDA speed test without the memory overhead

gcc-mp-4.8 -O3 cudadgemmtest.c common.c -o cudadgemmtest -I../../../../netlib/CBLAS -I/usr/local/cuda/include/ -L/usr/local/cuda/lib -lcublas
export DYLD_LIBRARY_PATH=/usr/local/cuda/lib
./cudadgemmtest  > ../../../results/mac_os_x-x86_64-dgemm-cuda_nooh.csv

*/

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

// does AB == C ? If not, complain on stderr
void test(int m, double* a, double *b, double *c) {
	int i, j, k, exact = 0, wrong = 0;
	double diff;
	double* d = calloc(m * m, sizeof(double));
	for (i = 0 ; i < m ; i++) {
		for (j = 0 ; j < m ; j++) {
			for (k = 0 ; k < m ; k++) {
				d[i + j * m] += a[i + k * m] * b[j * m + k];
			}
		}
	}
	for (i = 0 ; i < m ; i++) {
		for (j = 0 ; j < m ; j++) {
			diff = c[i * m + j] - d[i * m + j];
			if (diff != 0.0) {
				exact++;
			}
			if (abs(diff) > 0.000001) {
				wrong++;
			}
		}
	}
	free(d);
	if (wrong > 0) {
		fprintf(stderr, "not exact = %d, wrong = %d\n", exact, wrong);
	}
}

void checkStatus(char* message, cublasStatus status) {
    if (status != CUBLAS_STATUS_SUCCESS) {
    	fprintf (stderr, "!!!! %s fail %d\n", message, status);
    	exit(EXIT_FAILURE);
    }
}

long benchmark(int size) {
    int m = sqrt(size);
	long requestStart, requestEnd;

    double* a = random_array(m * m);
    double* b = random_array(m * m);
    double* c = calloc(m * m, sizeof(double));

	double *cuA, *cuB, *cuC;
	cublasStatus status;

	status = cublasAlloc(m * m, sizeof(double),(void**)&cuA);
	checkStatus("A", status);
	status = cublasAlloc(m * m, sizeof(double),(void**)&cuB);
	checkStatus("B", status);
	status = cublasAlloc(m * m, sizeof(double),(void**)&cuC);
	checkStatus("C", status);

	status = cublasSetMatrix(m, m, sizeof(double), a, m, cuA, m);
	checkStatus("setA", status);

	status = cublasSetMatrix(m, m, sizeof(double), b, m, cuB, m);
	checkStatus("setB", status);

	requestStart = currentTimeNanos();

	cublasDgemm('N', 'N', m, m, m, 1, cuA, m, cuB, m, 0, cuC, m);

	requestEnd = currentTimeNanos();

	status = cublasGetMatrix(m, m, sizeof(double), cuC, m, c, m);
	checkStatus("setB", status);

	status = cublasFree(cuA);
	checkStatus("freeA", status);
	status = cublasFree(cuB);
	checkStatus("freeB", status);
	status = cublasFree(cuC);
	checkStatus("freeC", status);

#ifdef __TEST__
	test(m, a, b, c);
#endif

    free(a);
    free(b);
    free(c);

    return (requestEnd - requestStart);
  }

main()
{
	cublasStatus status;

	srand(time(NULL));

	status = cublasInit();
	checkStatus("init",  status);

    double factor = 6.0 / 100.0;
    int i, j;
    for (i = 0 ; i < 10 ; i++) {
        for (j = 1 ; j <= 100 ; j++) {
            int size = (int) pow(10.0, factor * j);
            if (size < 10) continue;
            long took = benchmark(size);
            printf("\"%d\",\"%lu\"\n", size, took);
			fflush(stdout);
        }
    }
}


text:
static char help[] = "Solves the Liouville-Bratu reaction-diffusion problem in 1D.  Option prefix lb_.  Solves\n"
"    - u'' - lambda e^u = 0\n"
"on [0,1] subject to homogeneous Dirichlet boundary conditions.  Optionally uses manufactured solution to problem with f(x) on right-hand-side.\n\n";

#include <petsc.h>

typedef struct {
    PetscBool  manufactured;
    PetscReal  lambda;
} AppCtx;

extern PetscErrorCode Exact(DM, DMDALocalInfo*, Vec, AppCtx*);
extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*, PetscReal*,
                                 PetscReal*, AppCtx*);
extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*, PetscReal*,
                                 Mat, Mat, AppCtx*);

int main(int argc,char **args) {
  PetscErrorCode ierr;
  DM                  da;
  SNES                snes;
  AppCtx              user;
  Vec                 u, uexact;
  PetscReal           unorm, errnorm;
  DMDALocalInfo       info;

  PetscInitialize(&argc,&args,(char*)0,help);
  user.lambda = 1.0;
  user.manufactured = PETSC_FALSE;
  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,
                           "lb_","options for bratu1D",""); CHKERRQ(ierr);
  ierr = PetscOptionsReal("-lambda","coefficient of nonlinear zeroth-order term",
                          "bratu1D.c",user.lambda,&(user.lambda),NULL); CHKERRQ(ierr);
  ierr = PetscOptionsBool("-manu","if set, use a manufactured solution",
                          "bratu1D.c",user.manufactured,&(user.manufactured),NULL); CHKERRQ(ierr);
  ierr = PetscOptionsEnd(); CHKERRQ(ierr);

  ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,9,1,1,NULL,&da); CHKERRQ(ierr);
  ierr = DMSetFromOptions(da); CHKERRQ(ierr);
  ierr = DMSetUp(da); CHKERRQ(ierr);
  ierr = DMDASetUniformCoordinates(da,0.0,1.0,-1.0,-1.0,-1.0,-1.0); CHKERRQ(ierr);
  ierr = DMSetApplicationContext(da,&user); CHKERRQ(ierr);

  ierr = SNESCreate(PETSC_COMM_WORLD,&snes); CHKERRQ(ierr);
  ierr = SNESSetDM(snes,da); CHKERRQ(ierr);
  ierr = DMDASNESSetFunctionLocal(da,INSERT_VALUES,
             (DMDASNESFunction)FormFunctionLocal,&user); CHKERRQ(ierr);
  ierr = DMDASNESSetJacobianLocal(da,
             (DMDASNESJacobian)FormJacobianLocal,&user); CHKERRQ(ierr);
  ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);

  ierr = DMCreateGlobalVector(da,&u); CHKERRQ(ierr);
  ierr = VecSet(u,0.0); CHKERRQ(ierr);
  ierr = SNESSolve(snes,NULL,u); CHKERRQ(ierr);

  ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
  if (user.manufactured) {
      ierr = VecDuplicate(u,&uexact); CHKERRQ(ierr);
      ierr = Exact(da,&info,uexact,&user); CHKERRQ(ierr);
      ierr = VecNorm(u,NORM_INFINITY,&unorm); CHKERRQ(ierr);
      ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uxact
      ierr = VecNorm(u,NORM_INFINITY,&errnorm); CHKERRQ(ierr);
      ierr = PetscPrintf(PETSC_COMM_WORLD,
              "on %d point grid:  |u-u_exact|_inf/|u|_inf = %g\n",
              info.mx,errnorm/unorm); CHKERRQ(ierr);
      VecDestroy(&uexact);
  } else {
      ierr = PetscPrintf(PETSC_COMM_WORLD,"done on %d point grid\n",info.mx); CHKERRQ(ierr);
  }

  VecDestroy(&u);  SNESDestroy(&snes);  DMDestroy(&da);
  return PetscFinalize();
}

PetscErrorCode Exact(DM da, DMDALocalInfo *info, Vec uex, AppCtx *user) {
    PetscErrorCode ierr;
    PetscInt   i;
    PetscReal  h = 1.0 / (info->mx-1), x, *auex;
    ierr = DMDAVecGetArray(da,uex,&auex); CHKERRQ(ierr);
    for (i=info->xs; i<info->xs+info->xm; i++) {
        x = i * h;
        auex[i] = sin(PETSC_PI * x);
    }
    ierr = DMDAVecRestoreArray(da,uex,&auex); CHKERRQ(ierr);
    return 0;
}

PetscErrorCode FormFunctionLocal(DMDALocalInfo *info, PetscReal *u,
                                 PetscReal *f, AppCtx *user) {
    PetscInt   i;
    PetscReal  h = 1.0 / (info->mx-1), x, R, compf, uex;
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if ((i == 0) || (i == info->mx-1)) {
            f[i] = u[i];
        } else {  // interior location
            R = user->lambda * PetscExpReal(u[i]);
            f[i] = - u[i+1] + 2.0 * u[i] - u[i-1] - h*h * R;
            if (user->manufactured) {
                x = i * h;
                uex = sin(PETSC_PI * x);
                compf = PETSC_PI * PETSC_PI * uex - user->lambda * PetscExpReal(uex);
                f[i] -= h * h * compf;
            }
        }
    }
    return 0;
}

PetscErrorCode FormJacobianLocal(DMDALocalInfo *info, PetscReal *u,
                                 Mat J, Mat P, AppCtx *user) {
    PetscErrorCode ierr;
    PetscInt   i, col[3];
    PetscReal  h = 1.0 / (info->mx-1), dRdu, v[3];
    for (i=info->xs; i<info->xs+info->xm; i++) {
        if ((i == 0) | (i == info->mx-1)) {
            v[0] = 1.0;
            ierr = MatSetValues(P,1,&i,1,&i,v,INSERT_VALUES); CHKERRQ(ierr);
        } else {
            dRdu = user->lambda * PetscExpReal(u[i]);
            col[0] = i-1;  v[0] = - 1.0;
            col[1] = i;    v[1] = 2.0 - h*h * dRdu;
            col[2] = i+1;  v[2] = - 1.0;
            ierr = MatSetValues(P,1,&i,3,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }
    ierr = MatAssemblyBegin(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    ierr = MatAssemblyEnd(P,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
    if (J != P) {
        ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    }
    return 0;
}


text:
/// @file mlr.h
/// @brief multiple linear regression support
/// @author Jeff Perry <jeffsp@gmail.com>
/// @version 1.0
/// @date 2013-06-13

#ifndef MLR_H
#define MLR_H

#include "raster_utils.h"
//#include <lapacke.h>

namespace horny_toad
{

    /// @brief helper function
    template<typename T>
    T matrix_multiply (const T &a, const T &b)
    {
        assert (a.cols () == b.rows ());
        const size_t N = a.rows ();
        const size_t M = a.cols ();
        const size_t P = b.cols ();
        T y (a.rows (), b.cols ());
        for (size_t i = 0; i < N; ++i)
            for (size_t j = 0; j < P; ++j)
                for (size_t k = 0; k < M; ++k)
                    y (i, j) += a (i, k) * b (k, j);
        return y;
    }

    /// @brief multiple linear regression
    ///
    /// @tparam T matrix types
    /// @param y responses
    /// @param x predictors
    ///
    /// @return linear estimates of y=b*x
    template<typename T>
    T mlr_inverse (const T &y, const T &x)
    {
        // add a column of 1's to x on the left
        T xx (x.rows (), x.cols () + 1);
        for (size_t i = 0; i < x.rows (); ++i)
        {
            xx (i, 0) = 1;
            for (size_t j = 0; j < x.cols (); ++j)
            {
                xx (i, j + 1) = x (i, j);
            }
        }
        // b = (x^T * x)^-1 * x^T * y
        T xt = transpose (xx);
        T tmp = matrix_multiply (xt, xx);
        tmp = invert (tmp);
        tmp = matrix_multiply (tmp, xt);
        tmp = matrix_multiply (tmp, y);
        return tmp;
    }

    /// @brief multiple linear regression
    ///
    /// @tparam T matrix types
    /// @param y responses
    /// @param x predictors
    ///
    /// @return linear estimates of y=b*x
    /*
    template<typename T>
    T mlr_lapack (const T &y, const T &x)
    {
        assert (x.cols () <= x.rows ());
        assert (y.rows () == x.rows ());
        // result gets stored in b
        T b (y);
        // add a column of 1's to x on the left
        T xx (x.rows (), x.cols () + 1);
        for (size_t i = 0; i < x.rows (); ++i)
        {
            xx (i, 0) = 1;
            for (size_t j = 0; j < x.cols (); ++j)
            {
                xx (i, j + 1) = x (i, j);
            }
        }
        const size_t M = xx.rows ();
        const size_t N = xx.cols ();
        // solve it
        if (LAPACKE_dgels (LAPACK_ROW_MAJOR, 'N', M, N, 1, &xx[0], N, &b[0], 1) != 0)
            throw std::runtime_error ("LAPACK error: can't solve mlr");
        // copy result into matrix with correct dimensions
        T z (xx.cols (), 1);
        for (size_t i = 0; i < z.rows (); ++i)
            z[i] = b[i];
        return z;
    }
    */

    /// @brief multiple linear regression
    ///
    /// @tparam T matrix types
    /// @param y responses
    /// @param x predictors
    ///
    /// @return linear estimates of y=b*x
    template<typename T>
    T mlr (const T &y, const T &x)
    {
        return mlr_inverse (y, x);
        //return mlr_lapack (y, x);
    }
} // namespace horny_toad

#endif // MLR_H


text:
/****************************************************************************
 *                                                                          *
 *  Author : lukasz.iwaszkiewicz@gmail.com                                  *
 *  ~~~~~~~~                                                                *
 *  License : see COPYING file for details.                                 *
 *  ~~~~~~~~~                                                               *
 ****************************************************************************/

#pragma once

/**
 * This file is for maintaining compatibility with systems where C++ standard library is
 * not available i.e. for Arduino. I was totally unaware, that there is no C++ standard
 * library implemented for this toy platform.
 */

#if defined(__GNUC__) && (defined(AVR) || defined(ARDUINO) || defined (ARDUINO_ARCH_AVR))
#include <assert.h>
#include <stddef.h>
#include <stdint.h>

#include <etl/array.h>
#include <etl/map.h>
#include <etl/optional.h>

#if !defined Expects
#define Expects(x) assert (x)
#endif

namespace gsl {

template <class T, std::size_t N> constexpr T &at (T (&arr)[N], const int i)
{
        Expects (i >= 0 && i < static_cast<int> (N));
        return arr[static_cast<size_t> (i)];
}

template <class Cont> constexpr auto at (Cont &cont, const int i) -> decltype (cont[cont.size ()])
{
        Expects (i >= 0 && i < static_cast<int> (cont.size ()));
        using size_type = decltype (cont.size ());
        return cont[static_cast<size_type> (i)];
}
} // namespace gsl

namespace std {
template <typename _Tp, typename _Up = _Tp &&> _Up __declval (int);

template <typename _Tp> _Tp __declval (long);

template <typename _Tp> auto declval () noexcept -> decltype (__declval<_Tp> (0));

template <typename _Tp> struct __declval_protector {
        static const bool __stop = false;
};

template <typename _Tp> auto declval () noexcept -> decltype (__declval<_Tp> (0))
{
        static_assert (__declval_protector<_Tp>::__stop, "declval() must not be used!");
        return __declval<_Tp> (0);
}
} // namespace std
#else
#include <cstdint>
#include <gsl/gsl>
#endif

#include <etl/map.h>
#include <etl/optional.h>

max_stars_repo_name: semasquare/cpp-can-isotp

max_issues_repo_name: semasquare/cpp-can-isotp

max_forks_repo_name: semasquare/cpp-can-isotp