text:
/* specfunc/bessel_K1.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * Copyright (C) 2016 Pavel Holoborodko, 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/*
 Minimax rational approximation for [0,1), peak relative error = 1.83*GSL_DBL_EPSILON.
 Source: http://www.advanpix.com/?p=3987
*/
static double k1_poly[9] = {
  -3.0796575782920622440538935e-01,
  -8.5370719728650778045782736e-02,
  -4.6421827664715603298154971e-03,
  -1.1253607036630425931072996e-04,
  -1.5592887702110907110292728e-06,
  -1.4030163679125934402498239e-08,
  -8.8718998640336832196558868e-11,
  -4.1614323580221539328960335e-13,
  -1.5261293392975541707230366e-15
};

static double i1_poly[7] = {
  8.3333333333333325191635191e-02,
  6.9444444444467956461838830e-03,
  3.4722222211230452695165215e-04,
  1.1574075952009842696580084e-05,
  2.7555870002088181016676934e-07,
  4.9724386164128529514040614e-09
};

/*
 Chebyshev expansion for [1,8], peak relative error = 1.28*GSL_DBL_EPSILON.
 Source: Pavel Holoborodko.
*/
static double ak1_data[25] = {
  +2.07996868001418246e-01,
  +1.62581565017881476e-01,
  -5.87070423518863640e-03,
  +4.95021520115789501e-04,
  -5.78958347598556986e-05,
  +8.18614610209334726e-06,
  -1.31604832009487277e-06,
  +2.32546031520101213e-07,
  -4.42206518311557987e-08,
  +8.92163994883100361e-09,
  -1.89046270526983427e-09,
  +4.17568808108504702e-10,
  -9.55912361791375794e-11,
  +2.25769353153867758e-11,
  -5.48128000211158482e-12,
  +1.36386122546441926e-12,
  -3.46936690565986409e-13,
  +9.00354564415705942e-14,
  -2.37950577776254432e-14,
  +6.39447503964025336e-15,
  -1.74498363492322044e-15,
  +4.82994547989290473e-16,
  -1.35460927805445606e-16,
  +3.84604274446777234e-17,
  -1.10456856122581316e-17
};

static cheb_series ak1_cs = {
  ak1_data,
  24,
  -1, 1,
  9
};

/*
 Chebyshev expansion for [8,inf), peak relative error = 1.25*GSL_DBL_EPSILON.
 Source: SLATEC/dbsk1e.f
*/
static double ak12_data[14] = {
  +.637930834373900104E-1,
  +.283288781304972094E-1,
  -.247537067390525035E-3,
  +.577197245160724882E-5,
  -.206893921953654830E-6,
  +.973998344138180418E-8,
  -.558533614038062498E-9,
  +.373299663404618524E-10,
  -.282505196102322545E-11,
  +.237201900248414417E-12,
  -.217667738799175398E-13,
  +.215791416161603245E-14,
  -.229019693071826928E-15,
  +.258288572982327496E-16
};

static cheb_series ak12_cs = {
  ak12_data,
  13,
  -1, 1,
  7
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*GSL_DBL_MIN) {
    OVERFLOW_ERROR(result);
  }
  else if(x < 1.0) {
    const double lx = log(x);
    const double ex = exp(x);
    const double x2 = x*x;
    const double t  = 0.25*x2;
    const double i1 = 0.5 * x * (1.0 + t * (0.5 + t * gsl_poly_eval(i1_poly,6,t)));
    result->val  = ex * (x2 * gsl_poly_eval(k1_poly,9,x2) + x * lx * i1 + 1) / x;
    result->err  = ex * (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 8.0) {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak1_cs, (16.0/x-9.0)/7.0, &c);
    result->val  = (1.375 + c.val) / sx; /* 1.375 = 11/8 */
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak12_cs, 16.0/x-1.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*GSL_DBL_MIN) {
    OVERFLOW_ERROR(result);
  }
  else if(x < 1.0) {
    const double lx = log(x);
    const double x2 = x*x;
    const double t  = 0.25*x2;
    const double i1 = 0.5 * x * (1.0 + t * (0.5 + t * gsl_poly_eval(i1_poly,6,t)));
    result->val  = (x2 * gsl_poly_eval(k1_poly,9,x2) + x * lx * i1 + 1) / x;
    result->err  = (1.6+fabs(lx)*0.6) * GSL_DBL_EPSILON;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    gsl_sf_result K1_scaled;
    int stat_K1 = gsl_sf_bessel_K1_scaled_e(x, &K1_scaled);
    int stat_e  = gsl_sf_exp_mult_err_e(-x, 0.0,
                                           K1_scaled.val, K1_scaled.err,
                                           result);
    result->err = fabs(result->val) * (GSL_DBL_EPSILON*fabs(x) + K1_scaled.err/K1_scaled.val);
    return GSL_ERROR_SELECT_2(stat_e, stat_K1);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_K1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K1_scaled_e(x, &result));
}

double gsl_sf_bessel_K1(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K1_e(x, &result));
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#pragma once

#include <unordered_map>
#include <vector>

#include <gsl/span>
#include <json.hpp>
#include <optional.hpp>

#include "Quiver/Animation/AnimationData.h"
#include "Quiver/Animation/AnimationId.h"
#include "Quiver/Animation/Rect.h"
#include "Quiver/Animation/TimeUnit.h"

namespace qvr
{

// Tells us about where to find an Animation on disk.
struct AnimationSourceInfo {
	std::string name;
	std::string filename;
};

inline bool operator==(const AnimationSourceInfo& a, const AnimationSourceInfo& b) {
	return a.name == b.name && a.filename == b.filename;
}

class AnimationLibrary
{
public:
	AnimationId Add(const AnimationData& anim);
	AnimationId Add(
		const AnimationData& anim,
		const AnimationSourceInfo& sourceInfo);

	bool Remove(const AnimationId anim);

	bool Contains(const AnimationId anim) const;

	auto GetCount() const -> int;

	auto GetAnimation(const AnimationSourceInfo& sourceInfo)
		const -> AnimationId;

	auto GetSourceInfo(const AnimationId anim)
		const -> std::experimental::optional<AnimationSourceInfo>;

	auto GetFrameCount(const AnimationId anim) const -> int;
	auto GetViewCount(const AnimationId anim) const -> int;
	auto HasAltViews(const AnimationId anim) const -> bool;

	auto GetRect(
		const AnimationId anim,
		const int frameIndex,
		const int viewIndex = 0)
			const -> Animation::Rect;

	auto GetRects(
		const AnimationId anim,
		const int frameIndex)
			const -> gsl::span<const Animation::Rect>;

	auto GetTime(
		const AnimationId anim,
		const int frameIndex)
			const -> Animation::TimeUnit;

	auto GetIds() const -> std::vector<AnimationId>;

	friend void to_json(nlohmann::json& j, const AnimationLibrary& animations);

private:
	struct AnimationInfo {
		AnimationInfo(
			const unsigned indexOfFirstRect,
			const unsigned indexOfFirstTime,
			const unsigned numRects,
			const unsigned numRectsPerFrame)
			: mIndexOfFirstRect(indexOfFirstRect)
			, mIndexOfFirstTime(indexOfFirstTime)
			, mNumRects(numRects)
			, mNumRectsPerFrame(numRectsPerFrame)
		{}

		AnimationInfo(const AnimationInfo& other) = default;

		AnimationInfo() = default;

		unsigned NumFrames() const { return mNumRects / mNumRectsPerFrame; }

		std::experimental::optional<AnimationSourceInfo> mSourceInfo;

		unsigned mIndexOfFirstRect = 0;
		unsigned mIndexOfFirstTime = 0;
		unsigned mNumRects = 0;
		unsigned mNumRectsPerFrame = 0;
	};

	std::unordered_map<AnimationId, AnimationInfo> infos;

	// Time values for each frame in every animation.
	std::vector<Animation::TimeUnit> allFrameTimes;

	// Rects for each frame in every animation, including alt view rects.
	std::vector<Animation::Rect> allFrameRects;
};

void to_json(nlohmann::json& j, const AnimationLibrary& animations);
void to_json(nlohmann::json& j, const AnimationSourceInfo& sourceInfo);

void from_json(const nlohmann::json& j, AnimationLibrary& animations);
void from_json(const nlohmann::json& j, AnimationSourceInfo& animationSource);

}

ext: h

max_stars_repo_name: rachelnertia/Quarrel

max_issues_repo_name: rachelnertia/Quarrel

max_forks_repo_name: rachelnertia/Quiver


text:
#include <lapacke.h>
#include <mpi.h>
#include <stdio.h>

#include "../include/paralleltt.h"

#include <unistd.h>

// X is r1 x n x r2, transpose is r2 x n x r1
double* train_transpose(const double* X, int r1, int n, int r2){
    double* X_T = (double*) calloc(r1*n*r2, sizeof(double));

    for (int ii = 0; ii < r1; ++ii){
        for (int jj = 0; jj < n; ++jj){
            for (int kk = 0; kk < r2; ++kk){
                X_T[kk + jj*r2 + ii*n*r2] = X[ii + jj*r1 + kk*n*r1];
            }
        }
    }

    return X_T;
}

void two_sketches_to_train(sketch* s1, sketch* s2, double** train_ptr)
{
    double* train = *train_ptr;
    int head = 0;
    if ((s1 == NULL) || (s2 == NULL)){
        sketch* s = (s1 == NULL) ? s2 : s1;
        MPI_tensor* ten = s->ten;
        MPI_Comm comm = ten->comm;
        int iscol = s->iscol;
        int flattening = s->flattening;

        flattening_info* fi = flattening_info_init(ten, flattening, iscol, 0);

        int d = ten->d;
        int n = ten->n[(iscol) ? 0 : d-1];
        int r = s->r;

        matrix_tt* train_mat = (iscol) ? matrix_tt_wrap(n, r, train) : matrix_tt_wrap(r, n, train);
        matrix_tt_fill_zeros(train_mat);

        int* op = s->owner_partition;
        int rank = ten->rank;

        for (int ii = op[rank]; ii < op[rank+1]; ++ii){
            // X submatrix
            matrix_tt* X_submat = own_submatrix(s, ii, 0);

            // train submatrix
            flattening_info_f_update(fi, ten, ii);
            int i0 = fi->f_t_index[0];
            int i1 = i0 + fi->f_t_sizes[0];

            matrix_tt* train_submat = (iscol) ? submatrix(train_mat, i0, i1, 0, r) : submatrix(train_mat, 0, r, i0, i1);
            train_submat->transpose = (iscol) ? 0 : 1;
            matrix_tt_copy_data(train_submat, X_submat);

            free(train_submat);
            free(X_submat);
        }
        matrix_tt_reduce(comm, rank, train_mat, NULL, head);

        flattening_info_free(fi);
        free(train_mat);
        return;
    }

    int iscol = s1->iscol ? 1 : 0;
    int iscol2 = s2->iscol ? 1 : 0;

    if (iscol != iscol2){
        printf("two_sketches_to_train: Cannot convert a column and a row to a train\n");
        return;
    }

    int flattening1 = s1->flattening;
    int flattening2 = s2->flattening;

    // We want s2 to be the taller sketch, so that we only communicate s1 if necessary (which should be cheaper!)
    if (((iscol) && (flattening1 > flattening2)) || (!(iscol) && (flattening2 > flattening1))){
        sketch* tmp = s1;
        s1 = s2;
        s2 = tmp;

        flattening1 = s1->flattening;
        flattening2 = s2->flattening;
    }

    // We want the flattening dimensions to be one off of each other
    MPI_tensor* ten = s1->ten;
    flattening_info* fi1 = flattening_info_init(ten, flattening1, iscol, 0);
    flattening_info* fi2 = flattening_info_init(ten, flattening2, iscol, 0);
    if (fi2->f_d - fi1->f_d != 1){
        printf("two_sketches_to_train: The difference between the two flattening dimensions is not correct!\n");
        return;
    }

    MPI_Comm comm = ten->comm;
    int size = ten->comm_size;
    int rank = ten->rank;

    int r1 = s1->r;
    int r2 = s2->r;
    int n_index = (iscol) ? flattening1 : flattening2;
    int n = ten->n[n_index];
    matrix_tt* train_mat = matrix_tt_wrap(r1, n*r2, train);
    matrix_tt_fill_zeros(train_mat);

    // Loop over all blocks of fi2
    matrix_tt* s_mat1 = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    for (int block2 = 0; block2 < fi2->f_Nblocks; ++block2){
        flattening_info_f_update(fi2, ten, block2);
        int owner2 = s_get_owner(s2, block2);

        // Get the corresponding block1
        int* f_t_block1 = fi2->f_t_block + ((iscol) ? 0 : 1);
        int* f_nps1 = fi2->f_nps + ((iscol) ? 0 : 1);
        int block1 = 0;
        for (int ii = (fi1->f_d) - 1; ii >= 0; --ii){
            block1 = block1*f_nps1[ii] + f_t_block1[ii];
        }

        flattening_info_f_update(fi1, ten, block1);
        int owner1 = s_get_owner(s1, block1);
        sendrecv_sketch_block(s_mat1, s1, fi1, owner2, 0);

        if (rank == owner2){

            // Get submatrix of s2
            matrix_tt* s_mat2 = own_submatrix(s2, block2 , 0);

            // Multiply
            matrix_tt* train_submat = submatrix(train_mat, 0, 0, 0, 0);
            matrix_tt* s_submat2 = submatrix(s_mat2, 0, 0, 0, 0);
            s_submat2->transpose = (iscol) ? 0 : 1;
            for (int kk = 0; kk < r2; ++kk){
                int jj0 = (fi2->f_t_index[(iscol) ? (fi2->f_d) - 1 : 0]) + kk * n;
                int sz = fi2->f_t_sizes[(iscol) ? (fi2->f_d) - 1 : 0];
                submatrix_update(train_submat, 0, r1, jj0, jj0+sz);

                s_submat2->X = s_mat2->X + s_mat2->offset + kk*(s_mat2->lda);
                if (iscol){
                    s_submat2->m = s_mat1->m;
                    s_submat2->lda = s_mat1->m;
                    s_submat2->n = sz;
                }
                else{
                    s_submat2->m = sz;
                    s_submat2->lda = sz;
                    s_submat2->n = s_mat1->m;
                }

                s_mat1->transpose = 1;
                matrix_tt_dgemm(s_mat1, s_submat2, train_submat, 1.0, 1.0);
            }
            free(train_submat);
            free(s_mat2);
            free(s_submat2);
        }
    }
    free(s_mat1);

    matrix_tt_reduce(comm, rank, train_mat, NULL, head);

    // If it was a row sketch, we need to transpose
    if ((!iscol) && (rank == head)){
        double* train_cp = train_transpose(train, r1, n, r2);
        *train_ptr = train_cp;

        matrix_tt_free(train_mat); train_mat = NULL;
    }
    else{
        free(train_mat); train_mat = NULL;
    }

    flattening_info_free(fi1);
    flattening_info_free(fi2);
}

void PSTT2_final_train(tensor_train* tt, sketch** sketches, int mid)
{
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int head = 0;

    sketch* s1 = sketches[mid-1];
    int r1 = s1->r;
    sketch* s2 = sketches[mid];
    int r2 = s2->r;
    int n = tt->n[mid];

    MPI_tensor* ten = s1->ten;
    MPI_Comm comm = ten->comm;
    int size = ten->comm_size;
    int* schedule_rank = ten->schedule[rank];

    flattening_info* fi1 = flattening_info_init(ten, s1->flattening, s1->iscol, 0);
    flattening_info* fi2 = flattening_info_init(ten, s2->flattening, s2->iscol, 0);

    matrix_tt* train_mat = matrix_tt_wrap(r1*n, r2, tt->trains[mid]);
    matrix_tt_fill_zeros(train_mat);

    int tmp_n = 1;
    int* mid_partition = ten->partitions[mid];
    for (int ii = 0; ii < ten->nps[mid]; ++ii){
        int candidate = mid_partition[ii+1] - mid_partition[ii];
        tmp_n = (candidate > tmp_n) ? candidate : tmp_n;
    }
    matrix_tt* tmp = matrix_tt_init(r1, tmp_n);

    matrix_tt* recv_mat = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* s_mat1 = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* s_mat2 = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* s_submat2 = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* subtensor_mat = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    // Loop over schedule and stream
    for (int ii = 0; ii < ten->n_schedule; ++ii){
        int block = schedule_rank[ii];
        stream(ten, block);

        // Get the correct sketch matrices for multiplication
        for (int jj = 0; jj < size; ++jj){
            int* schedule_jj = ten->schedule[jj];
            int block_jj = schedule_jj[ii];

            if (block_jj != -1){
                // NOTE: Valgrind says I have a memory leak with these allocations.
                //       I think it's just confused by MPI, but maybe it's worth checking out
                flattening_info_update(fi1, ten, block_jj);
                sendrecv_sketch_block(recv_mat, s1, fi1, jj, 0);
                if (recv_mat->X != NULL){
                    matrix_tt* swtch = recv_mat;
                    recv_mat = s_mat1;
                    s_mat1 = swtch;
                }

                flattening_info_update(fi2, ten, block_jj);
                sendrecv_sketch_block(recv_mat, s2, fi2, jj, 0);
                if (recv_mat->X != NULL){
                    matrix_tt* swtch = recv_mat;
                    recv_mat = s_mat2;
                    s_mat2 = swtch;
                }
            }
        }

        flattening_info_update(fi1, ten, block);
        flattening_info_update(fi2, ten, block);

        // Multiply
        if (block != -1){
            int jj0 = fi1->t_t_index[mid];
            int sz = fi1->t_t_sizes[mid];
            submatrix_update(tmp, 0, r1, 0, sz);
            matrix_tt_wrap_update(subtensor_mat, fi1->f_N, fi1->s_N, get_X(ten));
            s_mat1->transpose = 1;
            submatrix_update(train_mat, jj0*r1, (jj0 + sz)*r1, 0, r2);
            submatrix_update(subtensor_mat, 0, 0, 0, 0);

            s_submat2->transpose = s_mat2->transpose;
            s_submat2->offset = 0;
            s_submat2->lda = s_mat2->lda;
            s_submat2->X_size = s_mat2->X_size;
            s_submat2->X = s_mat2->X + s_mat2->offset;

            for (int jj = 0; jj < fi2->f_N; ++jj){
                matrix_tt_reshape(r1, sz, tmp);
                submatrix_update(subtensor_mat, 0, fi1->f_N, jj*sz, (jj+1)*sz);

                matrix_tt_dgemm(s_mat1, subtensor_mat, tmp, 1.0, 0.0);
                matrix_tt_reshape(r1*sz, 1, tmp);

                submatrix_update(s_submat2, jj, jj+1, 0, r2);
                matrix_tt_dgemm(tmp, s_submat2, train_mat, 1.0, 1.0);
            }
        }

    }
    free(recv_mat);
    free(s_mat1);
    free(s_mat2);
    free(s_submat2);
    free(subtensor_mat);

    flattening_info_free(fi1);
    flattening_info_free(fi2);

    submatrix_update(train_mat, 0, n*r1, 0, r2);
    matrix_tt_reduce(comm, rank, train_mat, NULL, head);
    matrix_tt_free(tmp);

    free(train_mat);
}

int PSTT2_get_mid(MPI_tensor* ten, int mid)
{
    int d = ten->d;
    if ((mid > 0) && (mid < d)){
        return mid;
    }

    int* n = ten->n;

    long n_left = 1;
    long n_right = 1;
    for (int ii = 0; ii < d; ++ii){
        n_right = (long) n_right * n[ii];
    }

    for (int ii = 0; ii < d; ++ii){
        n_left = (long) n_left * n[ii];
        n_right = (long) n_right / n[ii];
        if ( (mid == -1) && (n_right < n_left) ){
            mid = ii;
        }
    }

    // Just in case, but this should never actually happen
    if (mid == -1){
        printf("You entered a weird tensor (probably n_i = 1 identically), setting mid = d-1\n");
        mid = d-1;
    }

    return mid;
}

void PSTT2(tensor_train* tt, MPI_tensor* ten, int mid)
{
    mid = PSTT2_get_mid(ten, mid);
    int rank = ten->rank;
    int BUF = 2;
    int* r = tt->r;
    int d = ten->d;

    // Create sketches
    sketch** sketches = (sketch**) calloc(d-1, sizeof(sketch*));
    for (int ii = 0; ii < d-1; ++ii){
        int iscol = (ii < mid) ? 1 : 0; // Column sketch below mid, row sketch above it
        sketches[ii] = sketch_init(ten, ii+1, r[ii+1], BUF, iscol);
    }

    // Sketch the tensor
    multi_perform_sketch(sketches, d-1);

    for (int ii = 0; ii < d-1; ++ii){
        sketch_qr(sketches[ii]);
    }

    // multiply out the sketches
    for (int ii = 0; ii < mid; ++ii){
        sketch* s1 = (ii == 0) ? NULL : sketches[ii-1];
        sketch* s2 = sketches[ii];
        two_sketches_to_train(s1, s2, tt->trains + ii);
    }

    for (int ii = mid; ii < d-1; ++ii){
        sketch* s1 = (ii == d-2) ? NULL : sketches[ii+1];
        sketch* s2 = sketches[ii];
        two_sketches_to_train(s1, s2, tt->trains + ii + 1);
    }

    // Get the middle train
    PSTT2_final_train(tt, sketches, mid);

    // free
    for (int ii = 0; ii < d-1; ++ii){
        sketch_free(sketches[ii]); sketches[ii] = NULL;
    }
    free(sketches);
}

void PSTT2_onepass_final_train(tensor_train* tt, sketch** sketches, int mid)
{
    sketch* sketch_mid = sketches[mid];
    sketch* Q_right = sketches[mid+1];
    sketch** Q_left = sketches + mid - 1; //Note that we are using the pointer here to feed into sketch_to_tensor

    // Multiply the middle not-QRed sketch against the sketch to its right
    int p = sketch_mid->r; // this is r_left + buf
    int r_left = tt->r[mid];
    int n_mid = tt->n[mid];
    int r_right = tt->r[mid+1];
    double* b = (double*) malloc(p*n_mid*r_right*sizeof(double)); // Named for the Ax=b solve that we will do later
    two_sketches_to_train(sketch_mid, Q_right, &b);

    MPI_tensor* ten_left = sketch_to_tensor(Q_left);

    sketch* sketch_left = sketch_init_with_Omega(ten_left, ten_left->d - 1, p, 0, 0, sketch_mid->Omegas); // buf = 0, is_col = 0
    perform_sketch(sketch_left);
    double* A = (double*) calloc(p * r_left, sizeof(double));

    two_sketches_to_train(sketch_left, NULL, &A); // This reduces the matrix A

    // Solve least squares
    matrix_tt* A_mat = matrix_tt_wrap(p, r_left, A);
    matrix_tt* b_mat = matrix_tt_wrap(p, n_mid*r_right, b);
    matrix_tt* x_mat = matrix_tt_wrap(r_left, n_mid*r_right, tt->trains[mid]);

    int head = 0;
    if (ten_left->rank == head){
        matrix_tt_dgels(x_mat, A_mat, b_mat);
    }

    // Free
    sketch_free(sketch_left); sketch_left = NULL;
    MPI_tensor_free(ten_left); ten_left = NULL;
    matrix_tt_free(A_mat); A_mat = NULL;
    matrix_tt_free(b_mat); b_mat = NULL;
    free(x_mat);
}

void PSTT2_onepass(tensor_train* tt, MPI_tensor* ten, int mid)
{
    // Number of times in PSTT2: 5
    // Number of times in multi_perform_sketch: 3
    mid = PSTT2_get_mid(ten, mid);
    int rank = ten->rank;
    int BUF = 2;
    int* r = tt->r;
    int d = ten->d;

    // Create sketches
    sketch** sketches = (sketch**) calloc(d, sizeof(sketch*));
    for (int ii = 0; ii < mid; ++ii){
        int iscol = 1; // Column sketch below mid, row sketch above it
        sketches[ii] = sketch_init(ten, ii+1, r[ii+1], BUF, iscol);
    }

    // putting the BUF into the rank here, so that two_sketches_to_train and perform_sketch don't get confused in the future
    sketches[mid] = sketch_init(ten, mid, r[mid] + BUF, 0, 0);

    for (int ii = mid + 1; ii < d; ++ii){
        int iscol = 0;
        sketches[ii] = sketch_init(ten, ii, r[ii], BUF, iscol);
    }

    // Sketch the tensor
    multi_perform_sketch(sketches, d);

    // Find orthogonal basis
    for (int ii = 0; ii < mid; ++ii){
        sketch_qr(sketches[ii]);
    }
    for (int ii = mid+1; ii < d; ++ii){
        sketch_qr(sketches[ii]);
    }

    // multiply out the sketches
    for (int ii = 0; ii < mid; ++ii){
        sketch* s1 = (ii == 0) ? NULL : sketches[ii-1];
        sketch* s2 = sketches[ii];
        two_sketches_to_train(s1, s2, tt->trains + ii);
    }

    for (int ii = mid+1; ii < d; ++ii){
        sketch* s1 = (ii == d-1) ? NULL : sketches[ii+1];
        sketch* s2 = sketches[ii];
        two_sketches_to_train(s1, s2, tt->trains + ii);
    }

    // Get the middle train
    PSTT2_onepass_final_train(tt, sketches, mid);

    // free
    for (int ii = 0; ii < d; ++ii){
        if(sketches[ii] != NULL){
            sketch_free(sketches[ii]); sketches[ii] = NULL;
        }
    }
    free(sketches);
}

ext: c


text:
#pragma once

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

#include <boost/test/data/test_case.hpp>
#include <boost/test/unit_test.hpp>

#include <cuda/runtime_api.hpp>

#include <thrustshift/managed-vector.h>
#include <thrustshift/copy.h>

namespace thrustshift {

template<class MemoryResource>
void touch_all_memory_resource_pages(MemoryResource& memory_resource) {
	auto device = cuda::device::current::get();
	auto stream = device.default_stream();

	for (const auto [k, v] : memory_resource.get_book()) {
		for (const auto& page : v) {
			BOOST_TEST(!page.allocated);
			using T = std::byte;
			const size_t N = k.bytes / sizeof(T);
			managed_vector<T> dst(N);
			gsl_lite::span<T> src(
			    reinterpret_cast<T*>(page.ptr), N);
			async::copy(stream, src, dst);
			stream.synchronize();
			for (size_t i = 0; i < N; ++i) {
				dst[i] = src[i];
			}
		}
	}
}

}

max_stars_repo_name: pauleonix/thrustshift

max_issues_repo_name: pauleonix/thrustshift

max_forks_repo_name: pauleonix/thrustshift


text:
/* vector/gsl_vector_long_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_VECTOR_LONG_DOUBLE_H__
#define __GSL_VECTOR_LONG_DOUBLE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_long_double.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  long double *data;
  gsl_block_long_double *block;
  int owner;
}
gsl_vector_long_double;

typedef struct
{
  gsl_vector_long_double vector;
} _gsl_vector_long_double_view;

typedef _gsl_vector_long_double_view gsl_vector_long_double_view;

typedef struct
{
  gsl_vector_long_double vector;
} _gsl_vector_long_double_const_view;

typedef const _gsl_vector_long_double_const_view gsl_vector_long_double_const_view;

/* Allocation */

GSL_FUN gsl_vector_long_double *gsl_vector_long_double_alloc (const size_t n);
GSL_FUN gsl_vector_long_double *gsl_vector_long_double_calloc (const size_t n);

GSL_FUN gsl_vector_long_double *gsl_vector_long_double_alloc_from_block (gsl_block_long_double * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_long_double *gsl_vector_long_double_alloc_from_vector (gsl_vector_long_double * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_long_double_free (gsl_vector_long_double * v);

/* Views */

GSL_FUN _gsl_vector_long_double_view
gsl_vector_long_double_view_array (long double *v, size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_vector_long_double_view_array_with_stride (long double *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_long_double_const_view_array (const long double *v, size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_long_double_const_view_array_with_stride (const long double *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_vector_long_double_subvector (gsl_vector_long_double *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_long_double_view
gsl_vector_long_double_subvector_with_stride (gsl_vector_long_double *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_long_double_const_subvector (const gsl_vector_long_double *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_long_double_const_view
gsl_vector_long_double_const_subvector_with_stride (const gsl_vector_long_double *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_long_double_set_zero (gsl_vector_long_double * v);
GSL_FUN void gsl_vector_long_double_set_all (gsl_vector_long_double * v, long double x);
GSL_FUN int gsl_vector_long_double_set_basis (gsl_vector_long_double * v, size_t i);

GSL_FUN int gsl_vector_long_double_fread (FILE * stream, gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_fwrite (FILE * stream, const gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_fscanf (FILE * stream, gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_fprintf (FILE * stream, const gsl_vector_long_double * v,
                              const char *format);

GSL_FUN int gsl_vector_long_double_memcpy (gsl_vector_long_double * dest, const gsl_vector_long_double * src);

GSL_FUN int gsl_vector_long_double_reverse (gsl_vector_long_double * v);

GSL_FUN int gsl_vector_long_double_swap (gsl_vector_long_double * v, gsl_vector_long_double * w);
GSL_FUN int gsl_vector_long_double_swap_elements (gsl_vector_long_double * v, const size_t i, const size_t j);

GSL_FUN long double gsl_vector_long_double_max (const gsl_vector_long_double * v);
GSL_FUN long double gsl_vector_long_double_min (const gsl_vector_long_double * v);
GSL_FUN void gsl_vector_long_double_minmax (const gsl_vector_long_double * v, long double * min_out, long double * max_out);

GSL_FUN size_t gsl_vector_long_double_max_index (const gsl_vector_long_double * v);
GSL_FUN size_t gsl_vector_long_double_min_index (const gsl_vector_long_double * v);
GSL_FUN void gsl_vector_long_double_minmax_index (const gsl_vector_long_double * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_long_double_add (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_FUN int gsl_vector_long_double_sub (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_FUN int gsl_vector_long_double_mul (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_FUN int gsl_vector_long_double_div (gsl_vector_long_double * a, const gsl_vector_long_double * b);
GSL_FUN int gsl_vector_long_double_scale (gsl_vector_long_double * a, const double x);
GSL_FUN int gsl_vector_long_double_add_constant (gsl_vector_long_double * a, const double x);

GSL_FUN int gsl_vector_long_double_isnull (const gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_ispos (const gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_isneg (const gsl_vector_long_double * v);
GSL_FUN int gsl_vector_long_double_isnonneg (const gsl_vector_long_double * v);

GSL_FUN INLINE_DECL long double gsl_vector_long_double_get (const gsl_vector_long_double * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_long_double_set (gsl_vector_long_double * v, const size_t i, long double x);
GSL_FUN INLINE_DECL long double * gsl_vector_long_double_ptr (gsl_vector_long_double * v, const size_t i);
GSL_FUN INLINE_DECL const long double * gsl_vector_long_double_const_ptr (const gsl_vector_long_double * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
long double
gsl_vector_long_double_get (const gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_long_double_set (gsl_vector_long_double * v, const size_t i, long double x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
long double *
gsl_vector_long_double_ptr (gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (long double *) (v->data + i * v->stride);
}

INLINE_FUN
const long double *
gsl_vector_long_double_const_ptr (const gsl_vector_long_double * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const long double *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_LONG_DOUBLE_H__ */

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

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

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

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file population.c
 * \brief Header file to define the population functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <glib.h>
#include "bits.h"
#include "entity.h"
#include "population.h"

/**
 * Function to init a population.
 *
 * \return 1 on succes, 0 on error.
 */
int
population_new (Population * population,        ///< Population.
                GeneticVariable * variable,     ///< Variables data.
                unsigned int nvariables,        ///< Number of variables.
                unsigned int genome_nbits,
                ///< Number of bits of each entity genome.
                unsigned int nentities, ///< Number of entities.
                double mutation_ratio,  ///< Mutation ratio.
                double reproduction_ratio,      ///< Reproduction ratio.
                double adaptation_ratio,        ///< Adaptation ratio.
                double threshold)
                ///< Threshold to finish the simulations.
{
  unsigned int i, nmutations, nreproductions, nadaptations;
  nmutations = mutation_ratio * nentities;
  nreproductions = reproduction_ratio * nentities;
  nadaptations = adaptation_ratio * nentities;
  i = nmutations + nreproductions + nadaptations;
  if (!i)
    {
      fprintf (stderr, "ERROR: no evolution\n");
      return 0;
    }
  if (i >= nentities)
    {
      fprintf (stderr, "ERROR: no survival of entities\n");
      return 0;
    }
  population->nsurvival = nentities - i;
  if (population->nsurvival < 2)
    {
      fprintf (stderr, "ERROR: unable to reproduce the entities\n");
      return 0;
    }
  population->variable = variable;
  population->nvariables = nvariables;
  population->nentities = nentities;
  population->mutation_max = nentities;
  population->mutation_min = population->reproduction_max
    = nentities - nmutations;
  population->reproduction_min = population->adaptation_max
    = population->reproduction_max - nreproductions;
  population->adaptation_min = population->nsurvival;
  population->genome_nbits = genome_nbits;
  population->genome_nbytes = bit_sizeof (genome_nbits);
  population->objective
		= (double *) g_slice_alloc (nentities * sizeof (double));
  population->entity = (Entity *) g_slice_alloc (nentities * sizeof (Entity));
  for (i = 0; i < population->nentities; ++i)
    {
      entity_new (population->entity + i, population->genome_nbytes, i);
      population->objective[i] = G_MAXDOUBLE;
    }
  population->threshold = threshold;
  population->stop = 0;
  return 1;
}

/**
 * Function to free the memory allocated in a population.
 */
void
population_init_genomes (Population * population,       ///< Population.
                         gsl_rng * rng) ///< GSL random numbers generator.
{
  unsigned int i;
  for (i = 0; i < population->nentities; ++i)
    entity_init (population->entity + i, rng);
}

/**
 * Function to free the memory allocated in a population.
 */
void
population_free (Population * population)       ///< Population.
{
  unsigned int i, nentities;
	nentities = population->nentities;
  for (i = 0; i < nentities; ++i)
    entity_free (population->entity + i);
  g_slice_free1 (nentities * sizeof (Entity), population->entity);
  g_slice_free1 (nentities * sizeof (double), population->objective);
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


text:
/* matrix/gsl_matrix_int.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_MATRIX_INT_H__
#define __GSL_MATRIX_INT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_int.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  int * data;
  gsl_block_int * block;
  int owner;
} gsl_matrix_int;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_view;

typedef _gsl_matrix_int_view gsl_matrix_int_view;

typedef struct
{
  gsl_matrix_int matrix;
} _gsl_matrix_int_const_view;

typedef const _gsl_matrix_int_const_view gsl_matrix_int_const_view;

/* Allocation */

GSL_FUN gsl_matrix_int *
gsl_matrix_int_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_int *
gsl_matrix_int_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_int *
gsl_matrix_int_alloc_from_block (gsl_block_int * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_int *
gsl_matrix_int_alloc_from_matrix (gsl_matrix_int * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_int *
gsl_vector_int_alloc_row_from_matrix (gsl_matrix_int * m,
                                        const size_t i);

GSL_FUN gsl_vector_int *
gsl_vector_int_alloc_col_from_matrix (gsl_matrix_int * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_int_free (gsl_matrix_int * m);

/* Views */

GSL_FUN _gsl_matrix_int_view
gsl_matrix_int_submatrix (gsl_matrix_int * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_row (gsl_matrix_int * m, const size_t i);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_column (gsl_matrix_int * m, const size_t j);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_diagonal (gsl_matrix_int * m);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_subdiagonal (gsl_matrix_int * m, const size_t k);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_superdiagonal (gsl_matrix_int * m, const size_t k);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_subrow (gsl_matrix_int * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_int_view
gsl_matrix_int_subcolumn (gsl_matrix_int * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_int_view
gsl_matrix_int_view_array (int * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_int_view
gsl_matrix_int_view_array_with_tda (int * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_int_view
gsl_matrix_int_view_vector (gsl_vector_int * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_int_view
gsl_matrix_int_view_vector_with_tda (gsl_vector_int * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_int_const_view
gsl_matrix_int_const_submatrix (const gsl_matrix_int * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_row (const gsl_matrix_int * m,
                            const size_t i);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_column (const gsl_matrix_int * m,
                               const size_t j);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_diagonal (const gsl_matrix_int * m);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_subdiagonal (const gsl_matrix_int * m,
                                    const size_t k);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_superdiagonal (const gsl_matrix_int * m,
                                      const size_t k);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_subrow (const gsl_matrix_int * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_int_const_view
gsl_matrix_int_const_subcolumn (const gsl_matrix_int * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_int_const_view
gsl_matrix_int_const_view_array (const int * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_int_const_view
gsl_matrix_int_const_view_array_with_tda (const int * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector (const gsl_vector_int * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_int_const_view
gsl_matrix_int_const_view_vector_with_tda (const gsl_vector_int * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_int_set_zero (gsl_matrix_int * m);
GSL_FUN void gsl_matrix_int_set_identity (gsl_matrix_int * m);
GSL_FUN void gsl_matrix_int_set_all (gsl_matrix_int * m, int x);

GSL_FUN int gsl_matrix_int_fread (FILE * stream, gsl_matrix_int * m) ;
GSL_FUN int gsl_matrix_int_fwrite (FILE * stream, const gsl_matrix_int * m) ;
GSL_FUN int gsl_matrix_int_fscanf (FILE * stream, gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_fprintf (FILE * stream, const gsl_matrix_int * m, const char * format);

GSL_FUN int gsl_matrix_int_memcpy(gsl_matrix_int * dest, const gsl_matrix_int * src);
GSL_FUN int gsl_matrix_int_swap(gsl_matrix_int * m1, gsl_matrix_int * m2);

GSL_FUN int gsl_matrix_int_swap_rows(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_int_swap_columns(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_int_swap_rowcol(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_int_transpose (gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_transpose_memcpy (gsl_matrix_int * dest, const gsl_matrix_int * src);

GSL_FUN int gsl_matrix_int_max (const gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_min (const gsl_matrix_int * m);
GSL_FUN void gsl_matrix_int_minmax (const gsl_matrix_int * m, int * min_out, int * max_out);

GSL_FUN void gsl_matrix_int_max_index (const gsl_matrix_int * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_int_min_index (const gsl_matrix_int * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_int_minmax_index (const gsl_matrix_int * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_int_equal (const gsl_matrix_int * a, const gsl_matrix_int * b);

GSL_FUN int gsl_matrix_int_isnull (const gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_ispos (const gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_isneg (const gsl_matrix_int * m);
GSL_FUN int gsl_matrix_int_isnonneg (const gsl_matrix_int * m);

GSL_FUN int gsl_matrix_int_add (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_FUN int gsl_matrix_int_sub (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_FUN int gsl_matrix_int_mul_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_FUN int gsl_matrix_int_div_elements (gsl_matrix_int * a, const gsl_matrix_int * b);
GSL_FUN int gsl_matrix_int_scale (gsl_matrix_int * a, const double x);
GSL_FUN int gsl_matrix_int_add_constant (gsl_matrix_int * a, const double x);
GSL_FUN int gsl_matrix_int_add_diagonal (gsl_matrix_int * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_int_get_row(gsl_vector_int * v, const gsl_matrix_int * m, const size_t i);
GSL_FUN int gsl_matrix_int_get_col(gsl_vector_int * v, const gsl_matrix_int * m, const size_t j);
GSL_FUN int gsl_matrix_int_set_row(gsl_matrix_int * m, const size_t i, const gsl_vector_int * v);
GSL_FUN int gsl_matrix_int_set_col(gsl_matrix_int * m, const size_t j, const gsl_vector_int * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL int   gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x);
GSL_FUN INLINE_DECL int * gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const int * gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
int
gsl_matrix_int_get(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_int_set(gsl_matrix_int * m, const size_t i, const size_t j, const int x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
int *
gsl_matrix_int_ptr(gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (int *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const int *
gsl_matrix_int_const_ptr(const gsl_matrix_int * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const int *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_INT_H__ */


text:
/* err/gsl_message.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_MESSAGE_H__
#define __GSL_MESSAGE_H__
#include <gsl/gsl_types.h>

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

__BEGIN_DECLS

/* Provide a general messaging service for client use.  Messages can
 * be selectively turned off at compile time by defining an
 * appropriate message mask. Client code which uses the GSL_MESSAGE()
 * macro must provide a mask which is or'ed with the GSL_MESSAGE_MASK.
 *
 * The messaging service can be completely turned off
 * by defining GSL_MESSAGING_OFF.  */

void gsl_message(const char * message, const char * file, int line,
                 unsigned int mask);

#ifndef GSL_MESSAGE_MASK
#define GSL_MESSAGE_MASK 0xffffffffu /* default all messages allowed */
#endif

GSL_VAR unsigned int gsl_message_mask ;

/* Provide some symolic masks for client ease of use. */

enum {
  GSL_MESSAGE_MASK_A = 1,
  GSL_MESSAGE_MASK_B = 2,
  GSL_MESSAGE_MASK_C = 4,
  GSL_MESSAGE_MASK_D = 8,
  GSL_MESSAGE_MASK_E = 16,
  GSL_MESSAGE_MASK_F = 32,
  GSL_MESSAGE_MASK_G = 64,
  GSL_MESSAGE_MASK_H = 128
} ;

#ifdef GSL_MESSAGING_OFF        /* throw away messages */
#define GSL_MESSAGE(message, mask) do { } while(0)
#else                           /* output all messages */
#define GSL_MESSAGE(message, mask) \
       do { \
       if (mask & GSL_MESSAGE_MASK) \
         gsl_message (message, __FILE__, __LINE__, mask) ; \
       } while (0)
#endif

__END_DECLS

#endif /* __GSL_MESSAGE_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
/* W. H. Bell
** A program to generate random numbers obeying a Landau distribution
** using the GNU scientific libraries.
*/

#include <stdio.h>

/* Header files needed for the Landau distribution generator */
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

int main() {
  int i;
  double ran_num;

  const gsl_rng_type * T;
  gsl_rng * r;

  /* Initialise the random number generator */
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  printf("Loop Num | Value\n");
  printf("----------------\n");
  for (i=0;i<2500;i++) {

    /* Generate a random number obeying the Landau distribution */
    ran_num=gsl_ran_landau (r);

    printf("%.4d     | %lf\n", i, ran_num);
  }

  return 0;
}

max_stars_repo_name: williamhbell/PhysCIntro

max_issues_repo_name: williamhbell/PhysCIntro

max_forks_repo_name: williamhbell/PhysCIntro