text:
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <assert.h>
#include <math.h>

#if  __INTEL_COMPILER
#include "mkl.h"
#else
#include <gsl/gsl_cblas.h>
#endif

#include "defs.h"
#include "utils.h"
#include "progressbar.h"

long naive(const double * restrict x1, const double * restrict y1, const int N, double * restrict d)
{
	long numcomputed=0;

	for(int i=0;i<N;i++) {
		for(int j=0;j<N;j++) {
			double *dist = d + i*N + j;
			const double *y0 = (const double *) &y1[j*NDIM];
			const double *x0 = (const double *) &x1[i*NDIM];
			const double dx = x0[0] - y0[0];
			const double dy = x0[1] - y0[1];
			const double dz = x0[2] - y0[2];
#ifndef SQRT_DIST
			*dist = dx*dx + dy*dy + dz*dz;
#else
			*dist = sqrt(dx*dx + dy*dy + dz*dz);
#endif
		}
 		numcomputed+=N;
	}

	return numcomputed;
}

long chunked(const double * restrict x1, const double * restrict y1, const int N, double * restrict d)
{
	long numcomputed = 0;
	const int block_size = 16;
	for(int i=0;i<N;i+=block_size) {
		const int block_size1 = (N-i) > block_size ? block_size:(N-i);
		for(int j=0;j<N;j+=block_size) {
			const int block_size2 = (N-j) > block_size ? block_size:(N-j);
			for(int ii=0;ii<block_size1;ii++) {
				const long index = (i+ii)*(long) N + j;
				double *dist = (double *) &d[index];
				for(int jj=0;jj<block_size2;jj++) {
					const double *y0 = (const double *) (&y1[(j+jj)*NDIM]);
					const double *x0 = (const double *) (&x1[(i+ii)*NDIM]);
					const double dx = x0[0] - y0[0];
					const double dy = x0[1] - y0[1];
					const double dz = x0[2] - y0[2];
					const double sqr_ds = dx*dx + dy*dy + dz*dz;
#ifndef SQRT_DIST
					dist[jj] = sqr_ds;
#else
					dist[jj] = sqrt(sqr_ds);
#endif
				}
				numcomputed += block_size2;
			}
		}
	}

	return numcomputed;
}

long compiler_vectorized_chunked(const double * restrict x1, const double * restrict y1, const int N, double * restrict d)
{
	long numcomputed = 0;
	const int block_size = 4;
	for(int i=0;i<N;i++) {
		int j;
		for(j=0;j<=(N-block_size);j+=block_size){
			double *dist = &d[i*N + j];
			for(int jj=0;jj<block_size;jj++) {
				const double *y0 = (const double *) (&y1[(j+jj)*NDIM]);
				const double *x0 = (const double *) (&x1[i*NDIM]);
				const double dx = x0[0] - y0[0];
				const double dy = x0[1] - y0[1];
				const double dz = x0[2] - y0[2];
				dist[jj] = dx*dx + dy*dy + dz*dz;
			}
			numcomputed+=block_size;
		}
		double *dist = (double *) &d[i*N + j];
		for(int jj=0;j<N;jj++,j++) {
			const double *y0 = (const double *) (&y1[j*NDIM]);
			const double *x0 = (const double *) (&x1[i*NDIM]);
			const double dx = x0[0] - y0[0];
			const double dy = x0[1] - y0[1];
			const double dz = x0[2] - y0[2];
			dist[jj] = dx*dx + dy*dy + dz*dz;
			numcomputed++;
		}
	}

#ifdef SQRT_DIST
	for(long index=0;index<numcomputed;index++){
		d[index] = sqrt(d[index]);
	}
#endif

	return numcomputed;
}

long blas_computed(const double * restrict x1, const double * restrict y1, const int N, double * restrict d)
{
	long numcomputed = (long) N * (long) N;
	const double alpha = -2.0;
	/* struct timeval t0,t1; */
	/* gettimeofday(&t0,NULL); */
	cblas_dgemm(CblasRowMajor,CblasNoTrans, CblasTrans,
							N, N, NDIM,
							alpha, x1, NDIM , y1, NDIM, 0.0, d, N);
	/* gettimeofday(&t1,NULL); */
	/* printf(" \t\t cblas_time = %6.2lf ms\n",1e3*(t1.tv_sec-t0.tv_sec) + 1e-3*(t1.tv_usec - t0.tv_usec)); */
	double *x0 = (double *) x1;
	for(int i=0;i<N;i++){
		const double sqr_x0_norm = x0[0]*x0[0] + x0[1]*x0[1] + x0[2]*x0[2];
		for(int j=0;j<N;j++) {
			const double *y0 = (const double *) &y1[j*NDIM];
			const double sqr_y0_norm = y0[0]*y0[0] + y0[1]*y0[1] + y0[2]*y0[2];
			double *dist = (double *) &d[i*N + j];
			*dist += (sqr_x0_norm + sqr_y0_norm);

#ifdef SQRT_DIST
			*dist = sqrt(*dist);
#endif
		}
		x0 += NDIM;
	}
	return numcomputed;
}

long check_result(double *dist, const double *x, const double *y, const int N)
{
	int bad=0;
	int numbadprinted=0;
	for(int i=0;i<N;i++) {
		for(int j=0;j<N;j++) {
			const double *tmp_x = (const double *) &(x[i*NDIM]);
			const double *tmp_y = (const double *) &(y[j*NDIM]);
			double sqr_ds = 0.0;
			for(int k=0;k<NDIM;k++) {
				const double dx = tmp_x[k] - tmp_y[k];
				sqr_ds += dx*dx;
			}
			const long index = i*(long) N + j;
			const double this_dist = dist[index];

#ifdef SQRT_DIST
			const double d = sqrt(sqr_ds);
#else
			const double d = sqr_ds;
#endif

			if(fabs(d - this_dist) > 1e-8) {
				if(numbadprinted < 20) {
					fprintf(stderr,"i = %d j = %d d = %0.8lf this_dist = %0.8lf \n",i,j,d,this_dist);
					numbadprinted++;
				}
				bad++;
			}
		}
	}

	return bad;
}

int main(int argc, char **argv)
{
	int numpart = 0;
	if(argc > 1 ) {
		numpart = atoi(argv[1]);
	} else {
		numpart = NELEMENTS;
	}

	const size_t numbytes = NDIM*numpart;
	double *x = calloc(sizeof(*x), numbytes);
	double *y = calloc(sizeof(*y), numbytes);
	const long totnpairs = (long) numpart * (long) numpart;
	double *dist = calloc(sizeof(*dist), totnpairs);
	assert(x != NULL && y != NULL && dist != NULL && "memory allocation failed");
  const char allfunction_names[][MAXLEN] = {"naive","chunked","compiler_vectorized_chunked","blas_computed"};
	const int ntests = sizeof(allfunction_names)/(sizeof(char)*MAXLEN);
	long (*allfunctions[]) (const double * restrict, const double * restrict, const int, double * restrict)      = {naive, chunked,compiler_vectorized_chunked,blas_computed};

	srand(seed);

	double function_best_mean_time[ntests],function_sigma_time[ntests],function_best_time_in_ms[ntests],function_best_mcycles[ntests];
	int function_niterations[ntests];
	for(int i=0;i<ntests;i++) {
		function_best_mean_time[i] = 1e16;
		function_sigma_time[i] = 0.0;
		function_best_time_in_ms[i] = 1e16;
		function_best_mcycles[i] = 1e16;
	}
#ifndef SQRT_DIST
	const long totflop = (long) numpart * (long) numpart * (8);
#else
	const long totflop = (long) numpart * (long) numpart * (8 + 10); //some hand-wavy thing saying that sqrt is 10 flop
#endif

	if(clustered_data == 0) {
		fill_array(x, numpart);
		fill_array(y, numpart);
	} else {
		assert(numpart <= max_galaxy_in_source_file && "Clustered data does not contain enough galaxies..please reduce numpart in defs.h");
		assert(NDIM == 3 && "Clustered galaxy data contains *exactly* 3 spatial dimensions");
		read_ascii(x, numpart, source_galaxy_file);
		read_ascii(y, numpart, source_galaxy_file);
	}

	const int64_t totniterations = repeat*ntests*(int64_t) max_niterations;
	int64_t numdone = 0;
	int interrupted=0;

	fprintf(stderr,"# Running benchmarks with N = %05d particles\n",numpart);
	init_my_progressbar(totniterations, &interrupted);
	for(int irep=0;irep<repeat;irep++) {
		for(int i=0;i<ntests;i++) {
			struct timeval t0,t1;
			double sum_x=0.0, sum_sqr_x=0.0;

			//warm-up
			gettimeofday(&t0,NULL);
			long ncomputed = (allfunctions[i]) (x, y, numpart, dist);
			gettimeofday(&t1,NULL);

			//check-result
			long numbad = check_result(dist, x, y, numpart);
			if(numbad != 0 || ncomputed != totnpairs) {
				fprintf(stderr,"ERROR: In function `%s' Number of incorrectly calculated distances = %ld out of a total of (%ld) possible pairs. Npairs computed = %ld\n", allfunction_names[i],numbad, totnpairs, ncomputed);
				goto cleanup;
			}

			double best_time_in_ms=1e16, best_time_in_megacycles=1e16;
			uint64_t start_cycles, end_cycles;
			const int64_t numdone_before_iter_loop = numdone;
			for(int iter=0;iter<max_niterations;iter++) {
				gettimeofday(&t0,NULL);
				start_cycles = rdtsc();
				(allfunctions[i]) (x, y, numpart, dist);
				end_cycles = rdtsc();
				gettimeofday(&t1,NULL);

				numdone++;
				my_progressbar(numdone,&interrupted);

				const double this_time_in_ms = 1e3*(t1.tv_sec-t0.tv_sec) +  1e-3*(t1.tv_usec - t0.tv_usec);
				if(this_time_in_ms < best_time_in_ms) best_time_in_ms = this_time_in_ms;
				const double this_time_in_megacycles = (end_cycles - start_cycles)/(1024.*1024.);
				if(this_time_in_megacycles < best_time_in_megacycles) best_time_in_megacycles = this_time_in_megacycles;

				sum_sqr_x += this_time_in_ms*this_time_in_ms;
				sum_x     += this_time_in_ms;

				if(max_niterations <= 10) {
					printf("     %-35s  %0.2lf \n",allfunction_names[i], this_time_in_ms);
				}

				if(best_time_in_ms < function_best_time_in_ms[i]) {
					function_best_time_in_ms[i] = best_time_in_ms;
				}

				if(best_time_in_megacycles < function_best_mcycles[i]) {
					function_best_mcycles[i] = best_time_in_megacycles;
				}

				const double mean_time  = sum_x/(iter+1.0);
				const double sigma_time = sqrt(sum_sqr_x/(iter+1.0) - mean_time*mean_time);
				if(mean_time < function_best_mean_time[i]) {
					function_best_mean_time[i] = mean_time;
				}
				function_niterations[i] = iter + 1;
				function_sigma_time[i] = sigma_time;

        //If the std.dev is small compared to typical runtime and
				//the code has run for more than XXX milli-seconds, then break
				if(sigma_time/mean_time < convergence_ratio && sum_x >= 300 && iter >= 5) {
					numdone = numdone_before_iter_loop + max_niterations;
					break;
				}
			}
		}//i loop over ntests
	}//irep loop over nrepeats
  finish_myprogressbar(&interrupted);

	printf("#######################################################\n");
	printf("##  Function                            Time (ms)      \n");
	printf("#######################################################\n");
	for(int i=0;i<ntests;i++) {
		const double mean_time = function_best_mean_time[i];
		const double sigma_time = function_sigma_time[i];
		const double gflops = (double) totflop/(1e-3*mean_time)/1e9;
		const double best_time_in_ms = function_best_time_in_ms[i];
		const double best_time_in_megacycles = function_best_mcycles[i];
		const int actual_niterations = function_niterations[i];
		printf(ANSI_COLOR_RED "# %-35s  %6.2lf +- %5.2lf " ANSI_COLOR_GREEN " (best -- %6.2lf ms, %6.2lf Mcycles) " ANSI_COLOR_RESET "," ANSI_COLOR_BLUE " %5.2lf GFlops [%04d iterations]" ANSI_COLOR_RESET "\n",
					 allfunction_names[i], mean_time, sigma_time, best_time_in_ms, best_time_in_megacycles, gflops, actual_niterations);

	}

cleanup:
	{
		free((void *) x);free((void *) y);free((void *) dist);
		exit(EXIT_FAILURE);
	}

	exit(EXIT_SUCCESS);
}

max_stars_repo_name: manodeep/pairwise

max_issues_repo_name: manodeep/pairwise

max_forks_repo_name: manodeep/pairwise


text:
// Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
// LICENCE available at https://www.gnu.org/software/gsl/doc/html/gpl.html
#ifndef TRIDIAG_H_
#define TRIDIAG_H_

#include <gsl/gsl_errno.h>

#include <cstdlib>
#include <iostream>
#include <vector>

namespace GSL {

int solve_tridiag_nonsym(const std::vector<double>& diag, const std::vector<double>& abovediag,
                         const std::vector<double>& belowdiag, const std::vector<double>& rhs, std::vector<double>& x,
                         size_t N);

int gsl_linalg_solve_tridiag(const std::vector<double>& diag, const std::vector<double>& abovediag,
                             const std::vector<double>& belowdiag, const std::vector<double>& rhs,
                             std::vector<double>& solution);

}  // namespace GSL

#endif  // TRIDIAG_H_

ext: h

max_stars_repo_name: carmeloevoli/SelfTevHalos

max_issues_repo_name: carmeloevoli/SelfTevHalos

max_forks_repo_name: carmeloevoli/SelfTevHalos


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  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.
 *
 * @brief scale decoder
 * @file ScaleDecoderStream.cpp
 */
#pragma once
#include "../../libutilities/Common.h"
#include "../../libutilities/DataConvertUtility.h"
#include "../../libutilities/FixedBytes.h"
#include "Common.h"
#include "FixedWidthIntegerCodec.h"
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/optional.hpp>
#include <boost/variant.hpp>
#include <array>
#include <gsl/span>

namespace bcos
{
namespace codec
{
namespace scale
{
class ScaleDecoderStream
{
public:
    // special tag to differentiate decoding streams from others
    static constexpr auto is_decoder_stream = true;
    explicit ScaleDecoderStream(gsl::span<byte const> span);

    /**
     * @brief scale-decodes pair of values
     * @tparam F first value type
     * @tparam S second value type
     * @param p pair of values to decode
     * @return reference to stream
     */
    template <class F, class S>
    ScaleDecoderStream& operator>>(std::pair<F, S>& p)
    {
        static_assert(!std::is_reference_v<F> && !std::is_reference_v<S>);
        return *this >> const_cast<std::remove_const_t<F>&>(p.first)  // NOLINT
               >> const_cast<std::remove_const_t<S>&>(p.second);      // NOLINT
    }

    /**
     * @brief scale-decoding of tuple
     * @tparam T enumeration of tuples types
     * @param v reference to tuple
     * @return reference to stream
     */
    template <class... T>
    ScaleDecoderStream& operator>>(std::tuple<T...>& v)
    {
        if constexpr (sizeof...(T) > 0)
        {
            decodeElementOfTuple<0>(v);
        }
        return *this;
    }

    /**
     * @brief scale-decoding of variant
     * @tparam T enumeration of various types
     * @param v reference to variant
     * @return reference to stream
     */
    template <class... Ts>
    ScaleDecoderStream& operator>>(boost::variant<Ts...>& v)
    {
        // first byte means type index
        uint8_t type_index = 0u;
        *this >> type_index;  // decode type index

        // ensure that index is in [0, types_count)
        if (type_index >= sizeof...(Ts))
        {
            BOOST_THROW_EXCEPTION(
                ScaleDecodeException() << errinfo_comment("exception for WRONG_TYPE_INDEX"));
        }

        tryDecodeAsOneOfVariant<0>(v, type_index);
        return *this;
    }

    /**
     * @brief scale-decodes shared_ptr value
     * @tparam T value type
     * @param v value to decode
     * @return reference to stream
     */
    template <class T>
    ScaleDecoderStream& operator>>(std::shared_ptr<T>& v)
    {
        using mutableT = std::remove_const_t<T>;

        static_assert(std::is_default_constructible_v<mutableT>);

        v = std::make_shared<mutableT>();
        return *this >> const_cast<mutableT&>(*v);  // NOLINT
    }

    /**
     * @brief scale-decodes unique_ptr value
     * @tparam T value type
     * @param v value to decode
     * @return reference to stream
     */
    template <class T>
    ScaleDecoderStream& operator>>(std::unique_ptr<T>& v)
    {
        using mutableT = std::remove_const_t<T>;

        static_assert(std::is_default_constructible_v<mutableT>);

        v = std::make_unique<mutableT>();
        return *this >> const_cast<mutableT&>(*v);  // NOLINT
    }

    /**
     * @brief scale-encodes any integral type including bool
     * @tparam T integral type
     * @param v value of integral type
     * @return reference to stream
     */
    template <typename T, typename I = std::decay_t<T>,
        typename = std::enable_if_t<std::is_integral<I>::value>>
    ScaleDecoderStream& operator>>(T& v)
    {
        // check bool
        if constexpr (std::is_same<I, bool>::value)
        {
            v = decodeBool();
            return *this;
        }
        // check byte
        if constexpr (sizeof(T) == 1u)
        {
            v = nextByte();
            return *this;
        }
        // decode any other integer
        v = decodeInteger<I>(*this);
        return *this;
    }

    /**
     * @brief scale-decodes any optional value
     * @tparam T type of optional value
     * @param v optional value reference
     * @return reference to stream
     */
    template <class T>
    ScaleDecoderStream& operator>>(boost::optional<T>& v)
    {
        using mutableT = std::remove_const_t<T>;

        static_assert(std::is_default_constructible_v<mutableT>);

        // optional bool is special case of optional values
        // it is encoded as one byte instead of two
        // as described in specification
        if constexpr (std::is_same<mutableT, bool>::value)
        {
            v = decodeOptionalBool();
            return *this;
        }
        // detect if optional has value
        bool has_value = false;
        *this >> has_value;
        if (!has_value)
        {
            v.reset();
            return *this;
        }
        // decode value
        v.emplace();
        return *this >> const_cast<mutableT&>(*v);  // NOLINT
    }

    ScaleDecoderStream& operator>>(u256& v);
    /**
     * @brief scale-decodes compact integer value
     * @param v compact integer reference
     * @return
     */
    ScaleDecoderStream& operator>>(CompactInteger& v);
    ScaleDecoderStream& operator>>(s256& v)
    {
        u256 unsignedValue;
        *this >> unsignedValue;
        v = u2s(unsignedValue);
        return *this;
    }

    template <unsigned N>
    ScaleDecoderStream& operator>>(FixedBytes<N>& fixedData)
    {
        bytes decodedData;
        *this >> decodedData;
        if (decodedData.size() < FixedBytes<N>::size)
        {
            BOOST_THROW_EXCEPTION(ScaleDecodeException() << errinfo_comment(
                                      "exception for invalid FixedBytes, expected size:" +
                                      std::to_string(FixedBytes<N>::size) +
                                      ", decoded data size:" + std::to_string(decodedData.size())));
        }
        fixedData = FixedBytes<N>(decodedData.data(), FixedBytes<N>::ConstructorType::FromPointer);
        return *this;
    }
    /**
     * @brief decodes vector of items
     * @tparam T item type
     * @param v reference to vector
     * @return reference to stream
     */
    template <class T>
    ScaleDecoderStream& operator>>(std::vector<T>& v)
    {
        using mutableT = std::remove_const_t<T>;
        using size_type = typename std::list<T>::size_type;

        static_assert(std::is_default_constructible_v<mutableT>);

        CompactInteger size{0u};
        *this >> size;

        auto item_count = size.convert_to<size_type>();
        std::vector<mutableT> vec;
        try
        {
            vec.resize(item_count);
        }
        catch (const std::bad_alloc&)
        {
            BOOST_THROW_EXCEPTION(
                ScaleDecodeException()
                << errinfo_comment("exception for TOO_MANY_ITEMS: " + std::to_string(item_count)));
        }
        if constexpr (sizeof(T) == 1u)
        {
            vec.assign(m_currentIterator, m_currentIterator + item_count);
            m_currentIterator += item_count;
            m_currentIndex += item_count;
        }
        else
        {
            for (size_type i = 0u; i < item_count; ++i)
            {
                *this >> vec[i];
            }
        }
        v = std::move(vec);
        return *this;
    }

    /**
     * @brief decodes map of pairs
     * @tparam T item type
     * @tparam F item type
     * @param m reference to map
     * @return reference to stream
     */
    template <class T, class F>
    ScaleDecoderStream& operator>>(std::map<T, F>& m)
    {
        using mutableT = std::remove_const_t<T>;
        static_assert(std::is_default_constructible_v<mutableT>);
        using mutableF = std::remove_const_t<F>;
        static_assert(std::is_default_constructible_v<mutableF>);

        using size_type = typename std::map<T, F>::size_type;

        CompactInteger size{0u};
        *this >> size;

        auto item_count = size.convert_to<size_type>();
        std::map<mutableT, mutableF> map;
        for (size_type i = 0u; i < item_count; ++i)
        {
            std::pair<mutableT, mutableF> p;
            *this >> p;
            map.emplace(std::move(p));
        }
        m = std::move(map);
        return *this;
    }

    /**
     * @brief decodes collection of items
     * @tparam T item type
     * @param v reference to collection
     * @return reference to stream
     */
    template <class T>
    ScaleDecoderStream& operator>>(std::list<T>& v)
    {
        using mutableT = std::remove_const_t<T>;
        using size_type = typename std::list<T>::size_type;

        static_assert(std::is_default_constructible_v<mutableT>);

        CompactInteger size{0u};
        *this >> size;

        auto item_count = size.convert_to<size_type>();

        std::list<T> lst;
        try
        {
            lst.reserve(item_count);
        }
        catch (const std::bad_alloc&)
        {
            BOOST_THROW_EXCEPTION(ScaleDecodeException() << errinfo_comment(
                                      "exception for TOO_MANY_ITEMS" + std::to_string(item_count)));
        }

        for (size_type i = 0u; i < item_count; ++i)
        {
            lst.emplace_back();
            *this >> lst.back();
        }
        v = std::move(lst);
        return *this;
    }

    /**
     * @brief decodes array of items
     * @tparam T item type
     * @tparam size of the array
     * @param a reference to the array
     * @return reference to stream
     */
    template <class T, size_t size>
    ScaleDecoderStream& operator>>(std::array<T, size>& a)
    {
        using mutableT = std::remove_const_t<T>;
        for (size_t i = 0u; i < size; ++i)
        {
            *this >> const_cast<mutableT&>(a[i]);  // NOLINT
        }
        return *this;
    }

    /**
     * @brief decodes string from stream
     * @param v value to decode
     * @return reference to stream
     */
    ScaleDecoderStream& operator>>(std::string& v);

    /**
     * @brief hasMore Checks whether n more bytes are available
     * @param n Number of bytes to check
     * @return True if n more bytes are available and false otherwise
     */
    bool hasMore(uint64_t n) const;

    /**
     * @brief takes one byte from stream and
     * advances current byte iterator by one
     * @return current byte
     */
    uint8_t nextByte()
    {
        if (!hasMore(1))
        {
            BOOST_THROW_EXCEPTION(ScaleDecodeException()
                                  << errinfo_comment("nextByte exception for NOT_ENOUGH_DATA"));
        }
        ++m_currentIndex;
        return *m_currentIterator++;
    }
    using SizeType = gsl::span<const byte>::size_type;

    gsl::span<byte const> span() const { return m_span; }
    SizeType currentIndex() const { return m_currentIndex; }

private:
    bool decodeBool();
    /**
     * @brief special case of optional values as described in specification
     * @return boost::optional<bool> value
     */
    boost::optional<bool> decodeOptionalBool();

    template <size_t I, class... Ts>
    void decodeElementOfTuple(std::tuple<Ts...>& v)
    {
        using T = std::remove_const_t<std::tuple_element_t<I, std::tuple<Ts...>>>;
        *this >> const_cast<T&>(std::get<I>(v));  // NOLINT
        if constexpr (sizeof...(Ts) > I + 1)
        {
            decodeElementOfTuple<I + 1>(v);
        }
    }

    template <size_t I, class... Ts>
    void tryDecodeAsOneOfVariant(boost::variant<Ts...>& v, size_t i)
    {
        using T = std::remove_const_t<std::tuple_element_t<I, std::tuple<Ts...>>>;
        static_assert(std::is_default_constructible_v<T>);
        if (I == i)
        {
            T val;
            *this >> val;
            v = std::forward<T>(val);
            return;
        }
        if constexpr (sizeof...(Ts) > I + 1)
        {
            tryDecodeAsOneOfVariant<I + 1>(v, i);
        }
    }

private:
    gsl::span<byte const> m_span;
    gsl::span<byte const>::const_iterator m_currentIterator;
    SizeType m_currentIndex;
};
}  // namespace scale
}  // namespace codec
}  // namespace bcos


text:
/* tt.h:
 /   This file contains everything related to the tensor_train struct
*/

#ifndef TENSOR_TRAIN_H
#define TENSOR_TRAIN_H

#include <lapacke.h>
#include <mpi.h>
#include "tensor.h"

typedef struct tensor_train {
    int d;           // dimension of tensor
    int* n;          // tensor size of each dimension
    int* r;          // tensor-train ranks

    double** trains; // array of pointers storing the address of the trains
} tensor_train;

tensor_train* TT_init(const int d, const int* restrict n);
tensor_train* TT_init_rank(const int d, const int* restrict n, const int* restrict r);
tensor_train* TT_copy(tensor_train* X);
void TT_free(tensor_train* sim);
void TT_print(tensor_train* tt);

// Broadcast a tensor train to all nodes
void tt_broadcast(MPI_Comm comm, tensor_train* tt);

// Find how much the tensor_train has compressed the tensor
double get_compression(const tensor_train* tt);

#endif


text:
#ifndef __GSL_MATRIX_H__
#define __GSL_MATRIX_H__

#include <gsl/matrix/gsl_matrix_double.h>

#endif /* __GSL_MATRIX_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*! Implementation of the class `Solver`.
 * \file solver.h
 */

#if !defined(SOLVER_H)
#define SOLVER_H

#include <petsc.h>

/**
 * \class Solver
 * \brief Super class for an iterative solver.
 */
class Solver
{
public:
  virtual ~Solver(){ }
  virtual PetscErrorCode create(const Mat &A) = 0;
  virtual PetscErrorCode solve(Vec &x, Vec &b) = 0;
  virtual PetscErrorCode getIters(PetscInt &iters) = 0;

}; // Solver

#endif

max_stars_repo_name: hietwll/PetIBM

max_issues_repo_name: hietwll/PetIBM

max_forks_repo_name: hietwll/PetIBM


text:
/* multifit/fdfridge.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_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_blas.h>

static int fdfridge_f(const gsl_vector * x, void * params, gsl_vector * f);
static int fdfridge_df(const gsl_vector * x, void * params, gsl_matrix * J);

gsl_multifit_fdfridge *
gsl_multifit_fdfridge_alloc (const gsl_multifit_fdfsolver_type * T,
                             const size_t n, const size_t p)
{
  gsl_multifit_fdfridge * work;

  work = calloc(1, sizeof(gsl_multifit_fdfridge));
  if (work == NULL)
    {
      GSL_ERROR_VAL("failed to allocate workspace",
                    GSL_ENOMEM, 0);
    }

  work->s = gsl_multifit_fdfsolver_alloc (T, n + p, p);
  if (work->s == NULL)
    {
      gsl_multifit_fdfridge_free(work);
      GSL_ERROR_VAL("failed to allocate space for fdfsolver",
                    GSL_ENOMEM, 0);
    }

  work->wts = gsl_vector_alloc(n + p);
  if (work->wts == NULL)
    {
      gsl_multifit_fdfridge_free(work);
      GSL_ERROR_VAL("failed to allocate space for weight vector",
                    GSL_ENOMEM, 0);
    }

  work->f = gsl_vector_alloc(n);
  if (work->f == NULL)
    {
      gsl_multifit_fdfridge_free(work);
      GSL_ERROR_VAL("failed to allocate space for f vector",
                    GSL_ENOMEM, 0);
    }

  work->n = n;
  work->p = p;
  work->lambda = 0.0;

  /* initialize weights to 1 (for augmented portion of vector) */
  gsl_vector_set_all(work->wts, 1.0);

  return work;
} /* gsl_multifit_fdfridge_alloc() */

void
gsl_multifit_fdfridge_free(gsl_multifit_fdfridge *work)
{
  if (work->s)
    gsl_multifit_fdfsolver_free(work->s);

  if (work->wts)
    gsl_vector_free(work->wts);

  if (work->f)
    gsl_vector_free(work->f);

  free(work);
}

const char *
gsl_multifit_fdfridge_name(const gsl_multifit_fdfridge * w)
{
  return gsl_multifit_fdfsolver_name(w->s);
}

gsl_vector *
gsl_multifit_fdfridge_position (const gsl_multifit_fdfridge * w)
{
  return gsl_multifit_fdfsolver_position(w->s);
}

gsl_vector *
gsl_multifit_fdfridge_residual (const gsl_multifit_fdfridge * w)
{
  return gsl_multifit_fdfsolver_residual(w->s);
}

size_t
gsl_multifit_fdfridge_niter (const gsl_multifit_fdfridge * w)
{
  return w->s->niter;
}

int
gsl_multifit_fdfridge_set (gsl_multifit_fdfridge * w,
                           gsl_multifit_function_fdf * f,
                           const gsl_vector * x,
                           const double lambda)
{
  return gsl_multifit_fdfridge_wset(w, f, x, lambda, NULL);
} /* gsl_multifit_fdfridge_set() */

int
gsl_multifit_fdfridge_wset (gsl_multifit_fdfridge * w,
                            gsl_multifit_function_fdf * f,
                            const gsl_vector * x,
                            const double lambda,
                            const gsl_vector * wts)
{
  const size_t n = w->n;
  const size_t p = w->p;

  if (n != f->n || p != f->p)
    {
      GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
    }
  else if (p != x->size)
    {
      GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
    }
  else if (wts != NULL && n != wts->size)
    {
      GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
    }
  else
    {
      int status;
      gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);

      /* save user defined fdf */
      w->fdf = f;

      /* build modified fdf for Tikhonov terms */
      w->fdftik.f = &fdfridge_f;
      w->fdftik.df = &fdfridge_df;
      w->fdftik.n = n + p; /* add p for Tikhonov terms */
      w->fdftik.p = p;
      w->fdftik.params = (void *) w;

      /* store damping parameter */
      w->lambda = lambda;
      w->L = NULL;

      if (wts)
        {
          /* copy weight vector into user portion of w->wts */
          gsl_vector_memcpy(&wv.vector, wts);
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
        }
      else
        {
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
        }

      /* update function/Jacobian evaluations */
      f->nevalf = w->fdftik.nevalf;
      f->nevaldf = w->fdftik.nevaldf;

      return status;
    }
} /* gsl_multifit_fdfridge_wset() */

int
gsl_multifit_fdfridge_set2 (gsl_multifit_fdfridge * w,
                            gsl_multifit_function_fdf * f,
                            const gsl_vector * x,
                            const gsl_vector * lambda)
{
  return gsl_multifit_fdfridge_wset2(w, f, x, lambda, NULL);
} /* gsl_multifit_fdfridge_set2() */

int
gsl_multifit_fdfridge_wset2 (gsl_multifit_fdfridge * w,
                             gsl_multifit_function_fdf * f,
                             const gsl_vector * x,
                             const gsl_vector * lambda,
                             const gsl_vector * wts)
{
  const size_t n = w->n;
  const size_t p = w->p;

  if (n != f->n || p != f->p)
    {
      GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
    }
  else if (p != x->size)
    {
      GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
    }
  else if (lambda->size != p)
    {
      GSL_ERROR ("lambda vector size does not match solver", GSL_EBADLEN);
    }
  else if (wts != NULL && n != wts->size)
    {
      GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
    }
  else
    {
      int status;
      gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);

      /* save user defined fdf */
      w->fdf = f;
      w->fdf->nevalf = 0;
      w->fdf->nevaldf = 0;

      /* build modified fdf for Tikhonov terms */
      w->fdftik.f = &fdfridge_f;
      w->fdftik.df = &fdfridge_df;
      w->fdftik.n = n + p; /* add p for Tikhonov terms */
      w->fdftik.p = p;
      w->fdftik.params = (void *) w;

      /* store damping matrix */
      w->lambda = 0.0;
      w->L_diag = lambda;
      w->L = NULL;

      if (wts)
        {
          /* copy weight vector into user portion */
          gsl_vector_memcpy(&wv.vector, wts);
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
        }
      else
        {
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
        }

      /* update function/Jacobian evaluations */
      f->nevalf = w->fdftik.nevalf;
      f->nevaldf = w->fdftik.nevaldf;

      return status;
    }
} /* gsl_multifit_fdfridge_wset2() */

int
gsl_multifit_fdfridge_set3 (gsl_multifit_fdfridge * w,
                            gsl_multifit_function_fdf * f,
                            const gsl_vector * x,
                            const gsl_matrix * L)
{
  return gsl_multifit_fdfridge_wset3(w, f, x, L, NULL);
} /* gsl_multifit_fdfridge_set3() */

int
gsl_multifit_fdfridge_wset3 (gsl_multifit_fdfridge * w,
                             gsl_multifit_function_fdf * f,
                             const gsl_vector * x,
                             const gsl_matrix * L,
                             const gsl_vector * wts)
{
  const size_t n = w->n;
  const size_t p = w->p;

  if (n != f->n || p != f->p)
    {
      GSL_ERROR ("function size does not match solver", GSL_EBADLEN);
    }
  else if (p != x->size)
    {
      GSL_ERROR ("vector length does not match solver", GSL_EBADLEN);
    }
  else if (L->size2 != p)
    {
      GSL_ERROR ("L matrix size2 does not match solver", GSL_EBADLEN);
    }
  else if (wts != NULL && n != wts->size)
    {
      GSL_ERROR ("weight vector length does not match solver", GSL_EBADLEN);
    }
  else
    {
      int status;
      gsl_vector_view wv = gsl_vector_subvector(w->wts, 0, n);

      /* save user defined fdf */
      w->fdf = f;
      w->fdf->nevalf = 0;
      w->fdf->nevaldf = 0;

      /* build modified fdf for Tikhonov terms */
      w->fdftik.f = &fdfridge_f;
      w->fdftik.df = &fdfridge_df;
      w->fdftik.n = n + p; /* add p for Tikhonov terms */
      w->fdftik.p = p;
      w->fdftik.params = (void *) w;

      /* store damping matrix */
      w->lambda = 0.0;
      w->L_diag = NULL;
      w->L = L;

      if (wts)
        {
          /* copy weight vector into user portion */
          gsl_vector_memcpy(&wv.vector, wts);
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, w->wts);
        }
      else
        {
          status = gsl_multifit_fdfsolver_wset(w->s, &(w->fdftik), x, NULL);
        }

      /* update function/Jacobian evaluations */
      f->nevalf = w->fdftik.nevalf;
      f->nevaldf = w->fdftik.nevaldf;

      return status;
    }
} /* gsl_multifit_fdfridge_wset3() */

int
gsl_multifit_fdfridge_iterate (gsl_multifit_fdfridge * w)
{
  int status = gsl_multifit_fdfsolver_iterate(w->s);

  /* update function/Jacobian evaluations */
  w->fdf->nevalf = w->fdftik.nevalf;
  w->fdf->nevaldf = w->fdftik.nevaldf;

  return status;
}

int
gsl_multifit_fdfridge_driver (gsl_multifit_fdfridge * w,
                              const size_t maxiter,
                              const double xtol,
                              const double gtol,
                              const double ftol,
                              int *info)
{
  int status = gsl_multifit_fdfsolver_driver(w->s, maxiter, xtol,
                                             gtol, ftol, info);
  return status;
} /* gsl_multifit_fdfridge_driver() */

/*
fdfridge_f()
  Callback function to provide residuals, including extra p
Tikhonov terms. The residual vector will have the form:

f~ = [     f     ]
     [ \lambda x ]

where f is the user supplied residuals, x are the model
parameters, and \lambda is the Tikhonov damping parameter

Inputs: x      - model parameters (size p)
        params - pointer to fdfridge workspace
        f      - (output) (n+p) vector to store f~
*/

static int
fdfridge_f(const gsl_vector * x, void * params, gsl_vector * f)
{
  int status;
  gsl_multifit_fdfridge *w = (gsl_multifit_fdfridge *) params;
  const size_t n = w->n;
  const size_t p = w->p;
  gsl_vector_view f_user = gsl_vector_subvector(f, 0, n);
  gsl_vector_view f_tik = gsl_vector_subvector(f, n, p);

  /* call user callback function to get residual vector f */
  status = gsl_multifit_eval_wf(w->fdf, x, NULL, &f_user.vector);
  if (status)
    return status;

  if (w->L_diag)
    {
      /* store diag(L_diag) x in Tikhonov portion of f~ */
      gsl_vector_memcpy(&f_tik.vector, x);
      gsl_vector_mul(&f_tik.vector, w->L_diag);
    }
  else if (w->L)
    {
      /* store Lx in Tikhonov portion of f~ */
      gsl_blas_dgemv(CblasNoTrans, 1.0, w->L, x, 0.0, &f_tik.vector);
    }
  else
    {
      /* store \lambda x in Tikhonov portion of f~ */
      gsl_vector_memcpy(&f_tik.vector, x);
      gsl_vector_scale(&f_tik.vector, w->lambda);
    }

  return GSL_SUCCESS;
} /* fdfridge_f() */

static int
fdfridge_df(const gsl_vector * x, void * params, gsl_matrix * J)
{
  int status;
  gsl_multifit_fdfridge *w = (gsl_multifit_fdfridge *) params;
  const size_t n = w->n;
  const size_t p = w->p;
  gsl_matrix_view J_user = gsl_matrix_submatrix(J, 0, 0, n, p);
  gsl_matrix_view J_tik = gsl_matrix_submatrix(J, n, 0, p, p);
  gsl_vector_view diag = gsl_matrix_diagonal(&J_tik.matrix);

  /* compute Jacobian */
  if (w->fdf->df)
    status = gsl_multifit_eval_wdf(w->fdf, x, NULL, &J_user.matrix);
  else
    {
      /* compute f(x) and then finite difference Jacobian */
      status = gsl_multifit_eval_wf(w->fdf, x, NULL, w->f);
      status += gsl_multifit_fdfsolver_dif_df(x, NULL, w->fdf, w->f,
                                              &J_user.matrix);
    }

  if (status)
    return status;

  if (w->L_diag)
    {
      /* store diag(L_diag) in Tikhonov portion of J */
      gsl_matrix_set_zero(&J_tik.matrix);
      gsl_vector_memcpy(&diag.vector, w->L_diag);
    }
  else if (w->L)
    {
      /* store L in Tikhonov portion of J */
      gsl_matrix_memcpy(&J_tik.matrix, w->L);
    }
  else
    {
      /* store \lambda I_p in Tikhonov portion of J */
      gsl_matrix_set_zero(&J_tik.matrix);
      gsl_vector_set_all(&diag.vector, w->lambda);
    }

  return GSL_SUCCESS;
} /* fdfridge_df() */

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * 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.
 */

#ifndef LIB_JSON_LOAD_JSON_ARRAY_LOADER_H_
#define LIB_JSON_LOAD_JSON_ARRAY_LOADER_H_

#include <glog/logging.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <jansson.h>
#include <stdint.h>

#include <array>
#include <string>
#include <vector>

#include "common/c_math/vec3.h"
#include "lib/json_load/json_load_or_die.h"

// TODO: Put this under the lib namespace.
namespace json_load {

// Utility class for loading JSON arrays into various vector and array
// structures such as Vec3, gsl_vector, and gsl_matrix.
class JsonArrayLoader {
 public:
  JsonArrayLoader(int32_t num_rows, int32_t num_cols, json_t *root)
      : root_(CHECK_NOTNULL(root)),
        num_rows_(num_rows),
        num_cols_(num_cols),
        matrix_buffer_(num_rows),
        buffer_(0) {
    CHECK_LT(0, num_rows);
    CHECK_LT(0, num_cols);
  }

  virtual ~JsonArrayLoader() {}

  template <int32_t rank>
  void LoadVector(const std::string &field_name,
                  const std::array<int32_t, rank> &dims, gsl_vector *dest) {
    CHECK_NOTNULL(dest);
    int32_t size = 1;
    for (int32_t i = 0; i < rank; ++i) size *= dims[i];
    CHECK_EQ(size, dest->size);

    LoadBuffer(field_name, size, &buffer_);
    CopyVector(buffer_, dest);
  }

  void LoadVector(const std::string &field_name, gsl_vector *dest) {
    LoadVector<1>(field_name, {{static_cast<int32_t>(dest->size)}}, dest);
  }

  void LoadVec3(const std::string &field_name, Vec3 *dest) {
    LoadBuffer(field_name, 3, &buffer_);
    dest->x = buffer_[0];
    dest->y = buffer_[1];
    dest->z = buffer_[2];
  }

  void LoadMatrix(const std::string &field_name, gsl_matrix *dest) {
    CHECK_NOTNULL(dest);

    json_t *field = json_load::LoadFieldOrDie(root_, field_name);
    CHECK(json_is_array(field));

    LoadMatrixFromField(field, dest);
  }

 protected:
  void LoadMatrixFromField(json_t *field, gsl_matrix *dest) {
    for (int32_t i = 0; i < num_rows_; ++i) {
      json_t *slice = json_array_get(field, i);
      CHECK(json_is_array(slice));
      CHECK_EQ(num_cols_, json_array_size(slice));
      json_load::LoadArray1D_DoubleOrDie(slice, num_cols_, &matrix_buffer_[i]);

      for (int32_t j = 0; j < num_cols_; ++j) {
        gsl_matrix_set(dest, i, j, matrix_buffer_[i][j]);
      }
    }
  }

  json_t *root_;

 private:
  void LoadBuffer(const std::string &field_name, int32_t size,
                  std::vector<double> *buffer) {
    CHECK_NOTNULL(buffer);
    json_t *array = json_load::LoadFieldOrDie(root_, field_name);
    CHECK(json_is_array(array));
    CHECK_EQ(size, json_array_size(array));
    json_load::LoadArray1D_DoubleOrDie(array, size, buffer);
  }

  void CopyVector(const std::vector<double> &buffer, gsl_vector *dest) {
    CHECK_NOTNULL(dest);
    CHECK_EQ(buffer.size(), dest->size);
    for (uint32_t i = 0U; i < buffer.size(); ++i) {
      gsl_vector_set(dest, i, buffer[i]);
    }
  }

  const int32_t num_rows_;
  const int32_t num_cols_;

  std::vector<std::vector<double>> matrix_buffer_;
  std::vector<double> buffer_;

  DISALLOW_COPY_AND_ASSIGN(JsonArrayLoader);
};

}  // namespace json_load

#endif  // LIB_JSON_LOAD_JSON_ARRAY_LOADER_H_


text:
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, 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.
//
// Eigen 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 Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN_GSL_HELPER
#define EIGEN_GSL_HELPER

#include <Eigen/Core>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_poly.h>

namespace Eigen {

template<typename Scalar, bool IsComplex = NumTraits<Scalar>::IsComplex> struct GslTraits
{
  typedef gsl_matrix* Matrix;
  typedef gsl_vector* Vector;
  static Matrix createMatrix(int rows, int cols) { return gsl_matrix_alloc(rows,cols); }
  static Vector createVector(int size) { return gsl_vector_alloc(size); }
  static void free(Matrix& m) { gsl_matrix_free(m); m=0; }
  static void free(Vector& m) { gsl_vector_free(m); m=0; }
  static void prod(const Matrix& m, const Vector& v, Vector& x) { gsl_blas_dgemv(CblasNoTrans,1,m,v,0,x); }
  static void cholesky(Matrix& m) { gsl_linalg_cholesky_decomp(m); }
  static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_cholesky_solve(m,b,x); }
  static void eigen_symm(const Matrix& m, Vector& eval, Matrix& evec)
  {
    gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc(m->size1);
    Matrix a = createMatrix(m->size1, m->size2);
    gsl_matrix_memcpy(a, m);
    gsl_eigen_symmv(a,eval,evec,w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
    gsl_eigen_symmv_free(w);
    free(a);
  }
  static void eigen_symm_gen(const Matrix& m, const Matrix& _b, Vector& eval, Matrix& evec)
  {
    gsl_eigen_gensymmv_workspace * w = gsl_eigen_gensymmv_alloc(m->size1);
    Matrix a = createMatrix(m->size1, m->size2);
    Matrix b = createMatrix(_b->size1, _b->size2);
    gsl_matrix_memcpy(a, m);
    gsl_matrix_memcpy(b, _b);
    gsl_eigen_gensymmv(a,b,eval,evec,w);
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
    gsl_eigen_gensymmv_free(w);
    free(a);
  }

  template<class EIGEN_VECTOR, class EIGEN_ROOTS>
  static void eigen_poly_solve(const EIGEN_VECTOR& poly, EIGEN_ROOTS& roots )
  {
    const int deg = poly.size()-1;
    double *z = new double[2*deg];
    double *a = new double[poly.size()];
    for( int i=0; i<poly.size(); ++i ){ a[i] = poly[i]; }

    gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc (poly.size());
    gsl_poly_complex_solve(a, poly.size(), w, z);
    gsl_poly_complex_workspace_free (w);
    for( int i=0; i<deg; ++i )
    {
      roots[i].real() = z[2*i];
      roots[i].imag() = z[2*i+1];
    }

    delete[] a;
    delete[] z;
  }
};

template<typename Scalar> struct GslTraits<Scalar,true>
{
  typedef gsl_matrix_complex* Matrix;
  typedef gsl_vector_complex* Vector;
  static Matrix createMatrix(int rows, int cols) { return gsl_matrix_complex_alloc(rows,cols); }
  static Vector createVector(int size) { return gsl_vector_complex_alloc(size); }
  static void free(Matrix& m) { gsl_matrix_complex_free(m); m=0; }
  static void free(Vector& m) { gsl_vector_complex_free(m); m=0; }
  static void cholesky(Matrix& m) { gsl_linalg_complex_cholesky_decomp(m); }
  static void cholesky_solve(const Matrix& m, const Vector& b, Vector& x) { gsl_linalg_complex_cholesky_solve(m,b,x); }
  static void prod(const Matrix& m, const Vector& v, Vector& x)
  { gsl_blas_zgemv(CblasNoTrans,gsl_complex_rect(1,0),m,v,gsl_complex_rect(0,0),x); }
  static void eigen_symm(const Matrix& m, gsl_vector* &eval, Matrix& evec)
  {
    gsl_eigen_hermv_workspace * w = gsl_eigen_hermv_alloc(m->size1);
    Matrix a = createMatrix(m->size1, m->size2);
    gsl_matrix_complex_memcpy(a, m);
    gsl_eigen_hermv(a,eval,evec,w);
    gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
    gsl_eigen_hermv_free(w);
    free(a);
  }
  static void eigen_symm_gen(const Matrix& m, const Matrix& _b, gsl_vector* &eval, Matrix& evec)
  {
    gsl_eigen_genhermv_workspace * w = gsl_eigen_genhermv_alloc(m->size1);
    Matrix a = createMatrix(m->size1, m->size2);
    Matrix b = createMatrix(_b->size1, _b->size2);
    gsl_matrix_complex_memcpy(a, m);
    gsl_matrix_complex_memcpy(b, _b);
    gsl_eigen_genhermv(a,b,eval,evec,w);
    gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);
    gsl_eigen_genhermv_free(w);
    free(a);
  }
};

template<typename MatrixType>
void convert(const MatrixType& m, gsl_matrix* &res)
{
//   if (res)
//     gsl_matrix_free(res);
  res = gsl_matrix_alloc(m.rows(), m.cols());
  for (int i=0 ; i<m.rows() ; ++i)
    for (int j=0 ; j<m.cols(); ++j)
      gsl_matrix_set(res, i, j, m(i,j));
}

template<typename MatrixType>
void convert(const gsl_matrix* m, MatrixType& res)
{
  res.resize(int(m->size1), int(m->size2));
  for (int i=0 ; i<res.rows() ; ++i)
    for (int j=0 ; j<res.cols(); ++j)
      res(i,j) = gsl_matrix_get(m,i,j);
}

template<typename VectorType>
void convert(const VectorType& m, gsl_vector* &res)
{
  if (res) gsl_vector_free(res);
  res = gsl_vector_alloc(m.size());
  for (int i=0 ; i<m.size() ; ++i)
      gsl_vector_set(res, i, m[i]);
}

template<typename VectorType>
void convert(const gsl_vector* m, VectorType& res)
{
  res.resize (m->size);
  for (int i=0 ; i<res.rows() ; ++i)
    res[i] = gsl_vector_get(m, i);
}

template<typename MatrixType>
void convert(const MatrixType& m, gsl_matrix_complex* &res)
{
  res = gsl_matrix_complex_alloc(m.rows(), m.cols());
  for (int i=0 ; i<m.rows() ; ++i)
    for (int j=0 ; j<m.cols(); ++j)
    {
      gsl_matrix_complex_set(res, i, j,
        gsl_complex_rect(m(i,j).real(), m(i,j).imag()));
    }
}

template<typename MatrixType>
void convert(const gsl_matrix_complex* m, MatrixType& res)
{
  res.resize(int(m->size1), int(m->size2));
  for (int i=0 ; i<res.rows() ; ++i)
    for (int j=0 ; j<res.cols(); ++j)
      res(i,j) = typename MatrixType::Scalar(
        GSL_REAL(gsl_matrix_complex_get(m,i,j)),
        GSL_IMAG(gsl_matrix_complex_get(m,i,j)));
}

template<typename VectorType>
void convert(const VectorType& m, gsl_vector_complex* &res)
{
  res = gsl_vector_complex_alloc(m.size());
  for (int i=0 ; i<m.size() ; ++i)
      gsl_vector_complex_set(res, i, gsl_complex_rect(m[i].real(), m[i].imag()));
}

template<typename VectorType>
void convert(const gsl_vector_complex* m, VectorType& res)
{
  res.resize(m->size);
  for (int i=0 ; i<res.rows() ; ++i)
    res[i] = typename VectorType::Scalar(
        GSL_REAL(gsl_vector_complex_get(m, i)),
        GSL_IMAG(gsl_vector_complex_get(m, i)));
}

}

#endif // EIGEN_GSL_HELPER

max_stars_repo_name: mathstuf/ParaView

max_issues_repo_name: mathstuf/ParaView

max_forks_repo_name: mathstuf/ParaView


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

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm(Quark *quark, Quark_Task_Flags *task_flags,
                      PLASMA_enum transA, int transB,
                      int m, int n, int k, int nb,
                      PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                                                const PLASMA_Complex64_t *B, int ldb,
                      PLASMA_Complex64_t beta, PLASMA_Complex64_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    C,                 INOUT,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm2( Quark *quark, Quark_Task_Flags *task_flags,
                        PLASMA_enum transA, int transB,
                        int m, int n, int k, int nb,
                        PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                        const PLASMA_Complex64_t *B, int ldb,
                        PLASMA_Complex64_t beta, PLASMA_Complex64_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    C,                 INOUT | LOCALITY | GATHERV,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_quark = PCORE_zgemm_quark
#define CORE_zgemm_quark PCORE_zgemm_quark
#endif
void CORE_zgemm_quark(Quark *quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int m;
    int n;
    int k;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t *A;
    int lda;
    PLASMA_Complex64_t *B;
    int ldb;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t *C;
    int ldc;

    quark_unpack_args_13(quark, transA, transB, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        m, n, k,
        CBLAS_SADDR(alpha), A, lda,
        B, ldb,
        CBLAS_SADDR(beta), C, ldc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm_f2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                                                   const PLASMA_Complex64_t *B, int ldb,
                         PLASMA_Complex64_t beta, PLASMA_Complex64_t *C, int ldc,
                         PLASMA_Complex64_t *fake1, int szefake1, int flag1,
                         PLASMA_Complex64_t *fake2, int szefake2, int flag2)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_f2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t)*nb*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(PLASMA_Complex64_t)*szefake1, fake1,             flag1,
        sizeof(PLASMA_Complex64_t)*szefake2, fake2,             flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_f2_quark = PCORE_zgemm_f2_quark
#define CORE_zgemm_f2_quark PCORE_zgemm_f2_quark
#endif
void CORE_zgemm_f2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t *B;
    int LDB;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t *C;
    int LDC;
    void *fake1, *fake2;

    quark_unpack_args_15(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1, fake2);
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm_p2(Quark *quark, Quark_Task_Flags *task_flags,
                         PLASMA_enum transA, int transB,
                         int m, int n, int k, int nb,
                         PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                         const PLASMA_Complex64_t **B, int ldb,
                         PLASMA_Complex64_t beta, PLASMA_Complex64_t *C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_p2_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_p2_quark = PCORE_zgemm_p2_quark
#define CORE_zgemm_p2_quark PCORE_zgemm_p2_quark
#endif
void CORE_zgemm_p2_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t **B;
    int LDB;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t *C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        *B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm_p3(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                           const PLASMA_Complex64_t *B, int ldb,
                           PLASMA_Complex64_t beta, PLASMA_Complex64_t **C, int ldc)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_p3_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t)*ldb*nb,   B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t*),         C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_p3_quark = PCORE_zgemm_p3_quark
#define CORE_zgemm_p3_quark PCORE_zgemm_p3_quark
#endif
void CORE_zgemm_p3_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t *B;
    int LDB;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t **C;
    int LDC;

    quark_unpack_args_13(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC);
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        B, LDB,
        CBLAS_SADDR(beta), *C, LDC);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zgemm_p2f1(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum transA, int transB,
                           int m, int n, int k, int nb,
                           PLASMA_Complex64_t alpha, const PLASMA_Complex64_t *A, int lda,
                           const PLASMA_Complex64_t **B, int ldb,
                           PLASMA_Complex64_t beta, PLASMA_Complex64_t *C, int ldc,
                           PLASMA_Complex64_t *fake1, int szefake1, int flag1)
{
    DAG_CORE_GEMM;
    QUARK_Insert_Task(quark, CORE_zgemm_p2f1_quark, task_flags,
        sizeof(PLASMA_enum),                &transA,    VALUE,
        sizeof(PLASMA_enum),                &transB,    VALUE,
        sizeof(int),                        &m,         VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(int),                        &k,         VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha,     VALUE,
        sizeof(PLASMA_Complex64_t)*lda*nb,   A,                 INPUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_Complex64_t*),         B,                 INPUT,
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,      VALUE,
        sizeof(PLASMA_Complex64_t)*ldc*nb,    C,                 INOUT | LOCALITY,
        sizeof(int),                        &ldc,       VALUE,
        sizeof(PLASMA_Complex64_t)*szefake1, fake1,             flag1,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zgemm_p2f1_quark = PCORE_zgemm_p2f1_quark
#define CORE_zgemm_p2f1_quark PCORE_zgemm_p2f1_quark
#endif
void CORE_zgemm_p2f1_quark(Quark* quark)
{
    PLASMA_enum transA;
    PLASMA_enum transB;
    int M;
    int N;
    int K;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t *A;
    int LDA;
    PLASMA_Complex64_t **B;
    int LDB;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t *C;
    int LDC;
    void *fake1;

    quark_unpack_args_14(quark, transA, transB, M, N, K, alpha,
                         A, LDA, B, LDB, beta, C, LDC, fake1);
    cblas_zgemm(
        CblasColMajor,
        (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB,
        M, N, K,
        CBLAS_SADDR(alpha), A, LDA,
        *B, LDB,
        CBLAS_SADDR(beta), C, LDC);
}


text:
#pragma once
#include "halley/utils/utils.h"
#include "halley/text/halleystring.h"
#include <memory>
#include <gsl/span>
#include "halley/resources/resource_data.h"

namespace Halley {
	enum class AssetType;
	class Deserializer;
	class Serializer;
	class AssetDatabase;
	class ResourceData;
	class ResourceDataReader;

	struct AssetPackHeader {
		std::array<char, 8> identifier;
		std::array<char, 16> iv;
		uint64_t assetDbStartPos;
		uint64_t dataStartPos;

		void init(size_t assetDbSize);
	};

    class AssetPack {
    public:
		AssetPack();
		AssetPack(const AssetPack& other) = delete;
		AssetPack(AssetPack&& other) noexcept;
		AssetPack(std::unique_ptr<ResourceDataReader> reader, const String& encryptionKey = "", bool preLoad = false);
		~AssetPack();

		AssetPack& operator=(const AssetPack& other) = delete;
		AssetPack& operator=(AssetPack&& other) noexcept;

		AssetDatabase& getAssetDatabase();
		const AssetDatabase& getAssetDatabase() const;
		Bytes& getData();
		const Bytes& getData() const;

		Bytes writeOut() const;

		std::unique_ptr<ResourceData> getData(const String& asset, AssetType type, bool stream);

		void readToMemory();
		void encrypt(const String& key);
		void decrypt(const String& key);

    	void readData(size_t pos, gsl::span<gsl::byte> dst);

		std::unique_ptr<ResourceDataReader> extractReader();

    private:
		std::unique_ptr<AssetDatabase> assetDb;
		std::unique_ptr<ResourceDataReader> reader;
		std::atomic<bool> hasReader;
		std::mutex readerMutex;
		size_t dataOffset = 0;
		Bytes data;
		std::array<char, 16> iv;
    };

	class PackDataReader final : public ResourceDataReader {
	public:
		PackDataReader(AssetPack& pack, size_t startPos, size_t fileSize);

		size_t size() const override;
		int read(gsl::span<gsl::byte> dst) override;
		void seek(int64_t pos, int whence) override;
		size_t tell() const override;
		void close() override;

	private:
		AssetPack& pack;
		const size_t startPos;
		const size_t fileSize;
		size_t curPos = 0;
		mutable std::mutex mutex;
	};
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley


text:
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"
#include "cblas.h"

void
cblas_ssbmv (const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
	     const int N, const int K, const float alpha, const float *A,
	     const int lda, const float *X, const int incX, const float beta,
	     float *Y, const int incY)
{
#define BASE float
#include "source_sbmv.h"
#undef BASE
}


text:
#if !defined(INDENTCLASSIFIER_H_INCLUDED)
#define INDENTCLASSIFIER_H_INCLUDED

#include "FileEnumerator.h"
#include "Utils.h"

#include <filesystem>
#include <gsl/gsl>
#include <iosfwd>
#include <map>
#include <string>

/// \brief IndentType contains the enum constants that indicate the type of
/// indentation for a line of text or a text file.
enum class IndentType
{
	space, ///< Indents consist entirely of spaces
	tab, ///< Indents consist entirely of tabs
	javadocTab,	///< Indents are entirely tabs or tabs followed by exactly one space and asterisk, as in a tab-indented JavaDoc comment
	javadocLeft,	///< Indented exactly one space followed by an asterisk, as in an unindented JavaDoc comment (single-line only)
	mixed, ///< Indents are mixed tabs and spaces, or some are tabs and others are spaces
	indeterminate ///< Refers to lines with no whitespace indent
};

using LineTypeCounts = ::std::map<IndentType, size_t>;

/// \brief The same as lineTypeCounts.at(indentType) except that it returns zero
/// if the map does not contain an entry for indentType.
size_t get(const LineTypeCounts& lineTypeCounts, IndentType indentType);

class IndentClassifier
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName, const char* pMsg);

	IndentClassifier(::gsl::span<const char*const> args);
	int run() const;

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

PRIVATE_EXCEPT_IN_TEST:
	using Path = ::std::filesystem::path;

	void processFile(const Path& p) const;
	static char indicatorLetter(IndentType iType);
	static ::std::string displayPath(const Path& p);

	/// \brief Scans the input stream and counts lines of the various indent types.
	static LineTypeCounts scanFile(const Path& p);
	static LineTypeCounts scanFile(::std::istream& in, bool isJavaFile);
	static IndentType classifyLine(const ::std::string& line, bool isJavaFile);
	static IndentType classifyFile(const LineTypeCounts& lineTypeCounts);

	FileEnumerator m_fileEnumerator;
};

#endif // INDENTCLASSIFIER_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


text:
#ifndef IBS_HPP
#define IBS_HPP

#include <array>
#include <assert.h>
#include <memory>
#include <vector>

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

#include "jspec2/force.h"
#include "jspec2/rate.h"

class Lattice;
class IonBeam;

enum class IBSModel {MARTINI, BM, BMC, BMZ};

class IBSSolver {
protected:
    double log_c_ = 0.0;     //Coulomb logarithm.
    double k_ = 0.0;          //Coupling rate in transverse directions.
    bool cache_invalid = true;
    bool ibs_by_element = false; //Calculate and output the ibs rate contribution element by element.

    void ibs_coupling(double &rx, double &ry, double k, double emit_x, double emit_y);
    void ibs_by_element_sddshead(std::ofstream& outfile, int n_element);
public:
    double log_c() const { return log_c_; }
    double k() const { return k_; }
    void set_k(double x) { k_ = x; }
    void set_log_c(double x) { log_c_ = x; }
    void set_ibs_by_element(bool b) {ibs_by_element = b;}
    void invalidate_cache() { cache_invalid = true; }

    IBSSolver(double log_c, double k);

    virtual rate3d rate(const Lattice &lattice, const IonBeam &beam) = 0;
};

class IBSSolver_Martini : public IBSSolver {
private:
    struct TrigonometryStorageUV {
        double sin_u2_cos_v2;
        double g1;
        double g2_1;
        double g2_2;
    };
    struct TrigonometryStorageV {
        double sin_v;
        double cos_v;
    };
    struct TrigonometryStorageU {
        double sin_u;
        double sin_u2;
        double cos_u2;
        double g3;
        std::vector<TrigonometryStorageUV> uv;
    };
    struct OpticalStorage {
        double a;
        double b2;
        double c2;
        double d2;
        double dtld;
        double k1;
        double k2;
        double k3;
    };

    int nu_ = 0;                //Grid number in u direction.
    int nv_ = 0;                //Grid number in v direction.
    int nz_ = 0;                //Grid number in z direction.

    // Scratch variables for IBS calculation (Martini model)
    std::vector<double> sigma_xbet, sigma_xbetp, sigma_y, sigma_yp;
    std::vector<TrigonometryStorageU> storage_u;
    std::vector<TrigonometryStorageV> storage_v;
    std::vector<OpticalStorage> storage_opt;
    std::vector<double> f1, f2, f3;

    void bunch_size(const Lattice &lattice, const IonBeam &beam);
    void abcdk(const Lattice &lattice, const IonBeam &beam);
    void coef_f();
    void f();
    double coef_a(const Lattice &lattice, const IonBeam &beam) const;
public:
    int nu() const { return nu_; }
    int nv() const { return nv_; }
    int nz() const { return nz_; }
    void set_nu(int nu) { assert(nu>0&&"Wrong value of nu in IBS parameters!"); nu_ = nu; invalidate_cache(); }
    void set_nv(int nv) { assert(nv>0&&"Wrong value of nv in IBS parameters!"); nv_ = nv; invalidate_cache(); }
    void set_nz(int nz) { assert(nz>0&&"Wrong value of nz in IBS parameters!"); nz_ = nz; invalidate_cache(); }
    IBSSolver_Martini(int nu, int nv, int nz, double log_c, double k);
    virtual rate3d rate(const Lattice &lattice, const IonBeam &beam) override;
};

class IBSSolver_BM : public IBSSolver {
 private:
     struct OpticalStorage { //variables only depends on the TWISS parameters and the energy.
         double phi;
         double dx2; //D_x * D_x
         double dx_betax_phi_2; // D_x * D_x / (beta_x * beta_x) + phi * phi
         double sqrt_betay; // sqrt(beta_y)
         double gamma_phi_2; // gamma * gamma * phi * phi
     };
     struct Kernels {
         double  psi;
         double sx;
         double sp;
         double sxp;
         double inv_sigma;
     };

     // Scratch variables for IBS calculation (Bjorken-Mtingwa model using Sergei Nagitsev's formula)
     std::vector<OpticalStorage> optical_strage;
     std::vector<Kernels> kernels;
     void init_fixed_var(const Lattice &lattice, const IonBeam &beam);
     void calc_kernels(const Lattice &lattice, const IonBeam &beam);
     double coef_bm(const Lattice &lattice, const IonBeam &beam) const;
 public:
     IBSSolver_BM(double log_c, double k);
     virtual rate3d rate(const Lattice &lattice, const IonBeam &beam) override;

};

class IBSSolver_BMZ : public IBSSolver {
private:
    int nt_;     //Number of steps for integration.
    struct optcl{
        double phi_x2;
        double phi_y2;
        double hx;
        double hy;
        double dx_2_over_beta_x;
        double dy_2_over_beta_y;
        double beta_x_over_hx;
        double hy_beta_x_over_hx;
        double beta_phi_x2;
        double beta_phi_y2;
        double hy_over_beta_y;
    };
    std::vector<optcl> optical;
    double factor = 3;
    void init_optical(const Lattice &lattice);
    double calc_abc(const Lattice &lattice, const IonBeam& beam, int i, double& a, double& b, double& c,
                               double& ax, double& bx, double& ay, double& by,double& al, double& bl);
    double coef(const Lattice &lattice, const IonBeam &beam) const;
    void calc_integral(double a, double b, double c, double ax, double bx, double ay, double by, double al,
                                  double bl, double& ix, double& iy, double& is, int nt, double u);
public:
    IBSSolver_BMZ(int nt, double log_c, double k);
    void set_nt(int n){assert(n>0&&"Wrong value of nt in IBS parameters!"); nt_ = n; invalidate_cache();}
    virtual rate3d rate(const Lattice &lattice, const IonBeam &beam) override;
    void set_factor(double x){factor = x;}
};

class IBSSolver_BM_Complete : public IBSSolver {
private:
    int nt_;
    double inv_ex;
    double inv_ey;
    double inv_dp2;
    double gamma;
    double gamma2;
    double factor = 3;
    struct Itgrl{
        double lambda;
        double lambda_sqrt;
        double ct;   // ct = 1/(1-t)^2
    };
    struct Optc {
        double phix;
        double phiy;
        double hx;
        double hy;
    };

    std::vector<Optc> optc;
    std::vector<Itgrl> itgrl;
    void init_optc(const Lattice &lattice);
    double det(std::array<std::array<double, 3>,3>& l);
    double inv(std::array<std::array<double, 3>,3>& l, std::array<std::array<double, 3>,3>& v);
    double trace(std::array<std::array<double, 3>,3>& l){return l[0][0]+l[1][1]+l[2][2];}
    void calc_l(const Lattice& lattice, int i, std::array<std::array<double, 3>,3>& lh,
                std::array<std::array<double, 3>,3>& lv, std::array<std::array<double, 3>,3>& ls);
    void calc_itgl(int i, std::array<std::array<double, 3>,3>& ii, std::array<std::array<double, 3>,3>& l,
                  std::array<std::array<double, 3>,3>& ll, std::array<std::array<double, 3>,3>& lh,
                  std::array<std::array<double, 3>,3>& lv, std::array<std::array<double, 3>,3>& ls);
    void calc_beam_const(const IonBeam& beam);
    double coef(const Lattice &lattice, const IonBeam &beam) const;
public:
     IBSSolver_BM_Complete(int nt, double log_c, double k);
     void set_factor(double x){factor = x;}
     virtual rate3d rate(const Lattice &lattice, const IonBeam &beam) override;

};
//
//class IBSSolver_BMZ2 : public IBSSolver {
//private:
//    int nt_;     //Number of steps for integration.
//    struct itgrl{
//        double lambda;
////        double lambda_2;
////        double lambda_3;
//        double lambda_sqrt;
//        double ct;   // ct = 1/(1-t)^2
//    };
//    double gamma_2;
//    double gamma_2_inv;
//    double gamma_4;
//    double emit_x_inv;
//    double emit_y_inv;
//    double emit_x2_inv;
//    double emit_y2_inv;
//
//    struct optcl{
//        double phi_x;
//        double phi_y;
//        double dx_2_over_beta_x;
//        double dy_2_over_beta_y;
//    };
//    std::vector<itgrl> integral;
//    std::vector<optcl> optical;
//    struct debug {
//        double  a, b, c, ax, bx, ay, by, al, bl, ix, iy, is;
//    };
//    std::vector<debug> my_debug;
//
//    gsl_integration_workspace *gw;
//
//    size_t limit = 100;
//    double espabs = 1e-12;
//    double esprel = 1e-3;
//    struct P{
//        double a;
//        double b;
//        double c;
//        double ai;
//        double bi;
//    }p;
////    double core(double q, void* params);
//
////    void init_integral(int n);
//    void init_optical(const Lattice &lattice);
//    double calc_abc(const Lattice &lattice, const Beam& beam, int i, double& a, double& b, double& c,
//                               double& ax, double& bx, double& ay, double& by,double& al, double& bl, double& scale);
//    double coef(const Lattice &lattice, const Beam &beam) const;
//    void calc_integral(double a, double b, double c, double ax, double bx, double ay, double by, double al,
//                                  double bl, double& ix, double& iy, double& is, int nt, std::vector<itgrl>& g);
//public:
//    IBSSolver_BMZ(int nt, double log_c, double k);
//    set_nt(int n){assert(n>0&&"Wrong value of nt in IBS parameters!"); nt_ = n; invalidate_cache();}
//    virtual void rate(const Lattice &lattice, const Beam &beam, double &rx, double &ry, double &rs);
//
//    ~IBSSolver_BMZ(){gsl_integration_workspace_free(gw);};
//
//};
//
//
//
//};
#endif


text:
/*
  C code for actionAngle calculations
*/
#ifndef __GALPY_ACTIONANGLE_H__
#define __GALPY_ACTIONANGLE_H__
#ifdef __cplusplus
extern "C" {
#endif
#ifdef _WIN32
#include <Python.h>
#endif
#include <stdbool.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
#include "interp_2d.h"
/*
  Macro for dealing with potentially unused variables due to OpenMP
 */
/* If we're not using GNU C, elide __attribute__ if it doesn't exist*/
#ifndef __has_attribute      // Compatibility with non-clang compilers.
#define __has_attribute(x) 0
#endif
#if defined(__GNUC__) || __has_attribute(unused)
#  define UNUSED __attribute__((unused))
#else
#  define UNUSED /*NOTHING*/
#endif
/*
  Structure declarations
*/
struct pragmasolver{
  gsl_root_fsolver *s;
};
#ifdef __cplusplus
}
#endif
#endif /* actionAngle.h */

max_stars_repo_name: gusbeane/galpy

max_issues_repo_name: gusbeane/galpy

max_forks_repo_name: gusbeane/galpy


text:
/*
 * Test the effective difference of a single bit changing on the input.
 *
 * This test uses random data as input. A hash is calculates.  A bit is
 * changed on the input at one location. Another hash is calculated.
 * Both hashes are XOR'ed and sent to output as a binary stream. This
 * is repeated until the output length is reached.
 */

#include "hash_bbt.h"

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

#include <stdint.h>
#include <stdio.h>
#include <string.h>

#define HASH_PARAMS bbt_table_1
extern struct bbt_hash_patterns HASH_PARAMS;

#define HASH_INPUT_SZ 32

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

	if (argc < 4) {
		printf("Usage: %s <outfile> <bit position> <output size> <numeric random seed>\n", argv[0]);
		return 1;
	}

	unsigned argPos = 1;
	FILE *outfile = fopen(argv[argPos++], "w");
	unsigned long int bitPos = strtoul(argv[argPos++], NULL, 10);
	if (bitPos > ((8*HASH_INPUT_SZ)-1)) {
		fprintf(stderr, "ERROR: bitPos too large. Beyond end of buffer length.\n");
		return 1;
	}
	unsigned long int outputLength = strtoul(argv[argPos++], NULL, 10);
	unsigned long int seed = strtoul(argv[argPos++], NULL, 10);

	unsigned char hashInput[HASH_INPUT_SZ];
	gsl_rng *rng = gsl_rng_alloc(gsl_rng_mt19937); /* Mersenne Twister */
	gsl_rng_set(rng, seed);

	bbt_hash_ctxt hc;
	bbt_hash_init(&hc, &(HASH_PARAMS));

	unsigned bitIndex = bitPos/8;
	unsigned bitShift = 7 - (bitPos%8);
	bbt_hash_t bitMask = ((bbt_hash_t)1) << bitShift;

	printf("Bit change at index %u and shift %u.\n", bitIndex, bitShift);

	for (unsigned outputWritten = 0; outputWritten < outputLength;) {

		unsigned char *p = hashInput;
		for (unsigned i = 0; i < HASH_INPUT_SZ; i++) {
			*p = gsl_rng_get(rng) & 0xFF;
			p++;
		}

		bbt_hash_t hashA;
		bbt_hash_t hashB;

		bbt_hash_calc(&hc, hashInput, HASH_INPUT_SZ);
		hashA = bbt_hash_getHash(&hc);
		bbt_hash_reset(&hc);

		hashInput[bitIndex] = hashInput[bitIndex] ^ bitMask;

		bbt_hash_calc(&hc, hashInput, HASH_INPUT_SZ);
		hashB = bbt_hash_getHash(&hc);
		bbt_hash_reset(&hc);

		bbt_hash_t hashCombined = hashA ^ hashB;
		for (unsigned i = 0; i < sizeof(bbt_hash_t); i++) {
			fputc(hashCombined & 0xFF, outfile);
			hashCombined = hashCombined >> 8;
			outputWritten++;
		}
	}

	fclose(outfile);
	gsl_rng_free(rng);
	return 0;
}

max_stars_repo_name: jonkanderson/BitBalancedTableHash

max_issues_repo_name: jonkanderson/BitBalancedTableHash

max_forks_repo_name: jonkanderson/BitBalancedTableHash


text:
/*
 *  maths_base.h
 *  ecosstat_project
 *
 *  Created by Karim Benabed on 23/06/09.
 *  Copyright 2009 Institut d'Astrophysique de Paris. All rights reserved.
 *
 */

#ifndef __MATHS_BASE_H
#define __MATHS_BASE_H

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#if 0
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_int.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#endif

#ifdef __PLANCK__
#include "HL2_likely/tools/errorlist.h"
#include "HL2_likely/tools/io.h"
#else
#include "errorlist.h"
#include "io.h"
#endif

#define math_base                -300
#define math_negative            -1 + math_base
#define math_singularValue       -2 + math_base
#define math_tooManySteps        -3 + math_base
#define math_underflow	         -4 + math_base
#define math_infnan	         -5 + math_base
#define math_wrongValue	         -6 + math_base
#define math_alloc	         -7 + math_base
#define math_interpoloutofrange	 -8 + math_base
#define math_interpol2small      -9 + math_base
#define math_interpol2big       -10 + math_base
#define math_stackTooSmall      -11 + math_base
#define math_overflow	        -12 + math_base
#define math_unknown            -13 + math_base

typedef double my_complex[2];

/* Mathematical constants */
#define pi     3.14159265358979323846
#define pi_sqr 9.86960440108935861883
#define twopi  6.28318530717958647693
#define ln2    0.69314718
#define ln2pi  1.837877066409
#define arcmin 2.90888208665721580e-4
#define arcsec 4.84813681e-6

/* Confidence levels, = erf({1,2,3}/sqrt(2)) */
#define conf_68 0.6827
#define conf_90 0.9000
#define conf_95 0.9545
#define conf_99 0.9973

/* Small numbers */
#define EPSILON  1.0e-5
#define EPSILON1 1.0e-8
#define EPSILON2 1.0e-18

/* Square of the absolute value of x (=|x|^2) */
#define ABSSQR(x) ((x)[0]*(x)[0] + (x)[1]*(x)[1])

/* Scalar product (x,y) */
#define SP(x,y) ((x)[0]*(y)[0]+(x)[1]*(y)[1])

#endif


text:
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <galpy_potentials.h>
// MovingObjectPotential
// 3 arguments: amp, t0, tf
void constrain_range(double * d) {
  // Constrains index to be within interpolation range
  if (*d < 0) *d = 0.0;
  if (*d > 1) *d = 1.0;
}
double MovingObjectPotentialRforce(double R,double z, double phi,
				   double t,
				   struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z, Rdist,RF;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate R force
  RF= calcRforce(Rdist,(obj_z-z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg);
  return -amp*RF*(cos(phi)*(obj_x-x)+sin(phi)*(obj_y-y))/Rdist;
}

double MovingObjectPotentialzforce(double R,double z,double phi,
				      double t,
				      struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z, Rdist;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate z force
  return -amp * calczforce(Rdist,(obj_z-z),phi,t,potentialArgs->nwrapped,
			   potentialArgs->wrappedPotentialArg);
}

double MovingObjectPotentialphiforce(double R,double z,double phi,
					double t,
					struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,obj_z, Rdist,RF;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range(&d_ind);
  // Interpolate x, y, z
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate phiforce
  RF= calcRforce(Rdist,(obj_z-z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg);
  return -amp*RF*R*(cos(phi)*(obj_y-y)-sin(phi)*(obj_x-x))/Rdist;
}

double MovingObjectPotentialPlanarRforce(double R, double phi,
				      double t,
				      struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,Rdist,RF;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range(&d_ind);
  // Interpolate x, y
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate R force
  RF= calcPlanarRforce(Rdist, phi, t, potentialArgs->nwrapped,
		       potentialArgs->wrappedPotentialArg);
  return -amp*RF*(cos(phi)*(obj_x-x)+sin(phi)*(obj_y-y))/Rdist;
}

double MovingObjectPotentialPlanarphiforce(double R, double phi,
					double t,
					struct potentialArg * potentialArgs){
  double amp,t0,tf,d_ind,x,y,obj_x,obj_y,Rdist,RF;
  double * args= potentialArgs->args;
  // Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  x= R*cos(phi);
  y= R*sin(phi);
  constrain_range(&d_ind);
  // Interpolate x, y
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  // Calculate phiforce
  RF= calcPlanarRforce(Rdist, phi, t, potentialArgs->nwrapped,
		       potentialArgs->wrappedPotentialArg);
  return -amp*RF*R*(cos(phi)*(obj_y-y)-sin(phi)*(obj_x-x))/Rdist;
}


text:
//IN method.
//Linear transformation (weights only, no biases) of Ni inputs to No outputs.
//This version uses CBLAS

//Input X has Ni neurons and output Y has No neurons.

//The vecs of length Ni are always contiguous in memory, such that:

//If col-major: Y[:,l] = W' * X[:,l]
//where:
//X has size Ni x L
//Y has size No x L
//W has size Ni x No

//If row-major: Y[l,:] = X[l,:] * W'
//X has size L x Ni
//Y has size L x No
//W has size No x Ni

//For a different set-up that allows linear transformation of vecs in
//any orientation, use the linear function from math.

//I retain the for loop through L for compatibility with real-time streaming.

#include <cblas.h>

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int linear_cblas_s (float *Y, const float *X, const float *W, const size_t Ni, const size_t No, const size_t L);
int linear_cblas_d (double *Y, const double *X, const double *W, const size_t Ni, const size_t No, const size_t L);
int linear_cblas_c (float *Y, const float *X, const float *W, const size_t Ni, const size_t No, const size_t L);
int linear_cblas_z (double *Y, const double *X, const double *W, const size_t Ni, const size_t No, const size_t L);

int linear_cblas_s (float *Y, const float *X, const float *W, const size_t Ni, const size_t No, const size_t L)
{
    for (size_t l=L; l>0u; --l, X+=Ni, Y+=No)
    {
        cblas_sgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,1.0f,W,(int)Ni,X,1,0.0f,Y,1);
    }

    return 0;
}

int linear_cblas_d (double *Y, const double *X, const double *W, const size_t Ni, const size_t No, const size_t L)
{
    for (size_t l=L; l>0u; --l, X+=Ni, Y+=No)
    {
        cblas_dgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,1.0,W,(int)Ni,X,1,0.0,Y,1);
    }

    return 0;
}

int linear_cblas_c (float *Y, const float *X, const float *W, const size_t Ni, const size_t No, const size_t L)
{
    const float z[2] = {0.0f,0.0f}, o[2] = {1.0f,0.0f};

    for (size_t l=L; l>0u; --l, X+=2u*Ni, Y+=2u*No)
    {
        cblas_cgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,o,W,(int)Ni,X,1,z,Y,1);
    }

    return 0;
}

int linear_cblas_z (double *Y, const double *X, const double *W, const size_t Ni, const size_t No, const size_t L)
{
    const double z[2] = {0.0,0.0}, o[2] = {1.0,0.0};

    for (size_t l=L; l>0u; --l, X+=2u*Ni, Y+=2u*No)
    {
        cblas_zgemv(CblasRowMajor,CblasNoTrans,(int)No,(int)Ni,o,W,(int)Ni,X,1,z,Y,1);
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif


text:
/*----------------------------------------------------------------------------
/
/ Filename: fits_wavelet.c
/ Author: Jay Billings
/ Author's email: jayjaybillings@gmail.com
/ Description: This program will perform a Debauchies wavelet transform on
/              an input image and turn off the wavelet coefficients up to
/              a certain percentage of the total size.
/
/              Usage:
/                 ./fits_wavelet <image> <limit>
/
/                 <image> should be an image in FITS format and <limit>
/                 should be a number between 0.0 and 1.0.
/
/              Output:
/                 filtered_<image>
/
/                 The output image is also in the FITS format and is readily
/                 viewable in FV, available from the NASA website, or the
/                 Sky Image Processor, useable as a java applet from Dr.
/                 John Simonetti, Virginia Tech,
/
/                 http://www.phys.vt.edu/~jhs/SIP.
/
/              This code requires the CFITSIO library and the GNU Scientific
/              Library. Both are freely available on the internet. This code
/              compiles fine with gcc-4.1.2 using
/
/                 gcc -o fits_wavelet fits_wavelet.c -lm -lgsl -lgslcblas
/                  -lcfitiso
/
/ Copyright (c) 2008       Jay Jay Billings
/
/ 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/
/ The license is also available at:
/
/		http://www.gnu.org/copyleft/gpl.html .
/
/ Date: 2008/03/24
/
*///--------------------------------------------------------------------------

#include <string.h>
#include <stdio.h>
#include <fitsio.h>
#include <gsl/gsl_wavelet.h>
#include <gsl/gsl_sort.h>
#include <math.h>

void read_image(double *,char *,int []);
void write_image(double *, char *, int []);

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

    int i, j, size[3], a;
    double f, limit, *pixels, *pixels_copy;
    gsl_wavelet *w;
    gsl_wavelet_workspace *work;
    size_t *p;

    // Check to make sure that the inputs are OK.
    if (argc < 3) {
       printf("You need to specify a file! \n\nUsage: ");
       printf("\t fits_wavelet <input_file> <limit>\n\n");
       exit(0);
    }
    f = atof(argv[2]);
    if (f >= 1.0) {
       printf("The limit must be between 0.0 and 1.0!\n\n");
       printf("Try something else! (Like 0.85...)\n");
       exit(0);
    }

    // Allocate some storage arrays. The maximum size is currently
    // set at 1024*1024, or 2^20, initilized to zero.
    if (!(pixels = calloc(1024*1024,sizeof(double)))) {
         printf("Memory allocation error\n");
         exit(0);
    }
    if (!(pixels_copy = calloc(1024*1024,sizeof(double)))) {
         printf("Memory allocation error\n");
         exit(0);
    }
    if (!(p = malloc(1024*1024*sizeof(size_t)))) {
         printf("Memory allocation error\n");
         exit(0);
    }

    // Read the image into the pixels array.
    read_image(pixels,argv[1],size);

    // Allocate the wavelet workspace.
    w = gsl_wavelet_alloc(gsl_wavelet_daubechies,20);
    work = gsl_wavelet_workspace_alloc(1024*1024);

    // Perform the wavelet transfer.
    gsl_wavelet_transform_forward(w,pixels,1,1024,work);

    // Make a copy of the pixels array and sort it.
    for (i = 0; i < size[2]; ++i) {
        pixels_copy[i] = fabs(pixels[i]);

    }
    gsl_sort_index(p,pixels_copy,1,size[0]*size[1]);

    // Loop over the wavelet coefficients and turn them off
    // up to a certain limit.
    limit = (double) size[2] * f;
    for (i = 0; i < (int) limit; ++i) {
        pixels[p[i]] = 0.0;
    }

    // Perform the inverse transform.
    gsl_wavelet_transform_inverse(w,pixels,1,1024,work);

    // Write the image.
    write_image(pixels,argv[1],size);

    // Free all the memory and exit.
    gsl_wavelet_free(w);
    gsl_wavelet_workspace_free(work);
    free(pixels);free(pixels_copy);free(p);
    return 0;
}

// This subroutine was adapted from documentation provided on the FITSIO website.
void read_image(double *pixels_vector, char *filename, int size[3])
{
    fitsfile *fptr;   /* FITS file pointer, defined in fitsio.h */
    int status = 0;   /* CFITSIO status value MUST be initialized to zero! */
    int bitpix, naxis, ii, anynul,a;
    long naxes[2] = {1,1}, fpixel[2] = {1,1};
    double *pixels;
    char format[20], hdformat[20];

    /*if (argc != 2) {
      printf("Usage:  imlist filename[ext][section filter] \n");
      printf("\n");
      printf("List the the pixel values in a FITS image \n");
      printf("\n");
      printf("Example: \n");
      printf("  imlist image.fits                    - list the whole image\n");
      printf("  imlist image.fits[100:110,400:410]   - list a section\n");
      printf("  imlist table.fits[2][bin (x,y) = 32] - list the pixels in\n");
      printf("         an image constructed from a 2D histogram of X and Y\n");
      printf("         columns in a table with a binning factor = 32\n");
      return(0);
    }*/

    if (!fits_open_file(&fptr, filename, READONLY, &status))
    {
        if (!fits_get_img_param(fptr, 2, &bitpix, &naxis, naxes, &status) )
        {
          if (naxis > 2 || naxis == 0)
             printf("Error: only 1D or 2D images are supported\n");
          else
          {
            /* get memory for 1 row */
            if (!(pixels = malloc(naxes[0] * sizeof(double)))) {
                printf("Memory allocation error\n");
                exit(0);
            }

            if (bitpix > 0) {  /* set the default output format string */
               strcpy(hdformat, " %7d");
               strcpy(format,   " %7.0f");
            } else {
               strcpy(hdformat, " %15d");
               strcpy(format,   " %15.5f");
            }

            //printf("\n      ");          /* print column header */
            //for (ii = 1; ii <= naxes[0]; ii++)
            //   printf(hdformat, ii);
            //printf("\n");                /* terminate header line */

            /* loop over all the rows in the image, top to bottom */
            a = 0;
            for (fpixel[1] = naxes[1]; fpixel[1] >= 1; fpixel[1]--)
            {
               if (fits_read_pix(fptr, TDOUBLE, fpixel, naxes[0], NULL,
                    pixels, NULL, &status) )  /* read row of pixels */
                  break;  /* jump out of loop on error */

            //   printf(" %4d \n",fpixel[1]);  /* print row number */
               for (ii = 0; ii < naxes[0]; ii++) {
                  pixels_vector[a] = pixels[ii];
//                  printf(format, pixels_vector[a]);   /* print each value  */
 //                 printf(" %d\n",a);                    /* terminate line */
                  ++a;
               }
            }
            size[0] = naxes[0];
            size[1] = naxes[1];
            size[2] = naxes[0]*naxes[1];
            free(pixels);
          }
        }
        fits_close_file(fptr, &status);
    }

    if (status) {
       fits_report_error(stderr, status); /* print any error message */
       exit(0);
    }
    //return(status);
    return;
}

void write_image(double *pixels_vector, char *f1, int size[3]) {

    fitsfile *infptr, *outfptr;   /* FITS file pointer, defined in fitsio.h */
    int status = 0;   /* CFITSIO status value MUST be initialized to zero! */
    int ii = 1;
    char f2[200];
    long size_l[2], fpixel[2] = {1,1};
    void *pix_ptr = &pixels_vector;

    size_l[0] = (long) size[0];
    size_l[1] = (long) size[1];

    sprintf(f2,"filtered_%s",f1);

    /* Create the output file */
    if (!fits_create_file(&outfptr, f2, &status) ) {
       if (!fits_create_img(outfptr,32,2,size_l,&status)) {
          if (!fits_write_pix(outfptr,TDOUBLE,fpixel,(long) size[2],pixels_vector \
           ,&status)) {
          if (status == END_OF_FILE) status = 0;
          }
          //printf("Write status: %d \n",status);
       }
       fits_close_file(outfptr,  &status);
    }

    return;
}