text:
#pragma once

#include "Cesium3DTiles/BoundingVolume.h"
#include "Cesium3DTiles/Gltf.h"
#include "Cesium3DTiles/Library.h"
#include "Cesium3DTiles/TileContentLoadResult.h"
#include "Cesium3DTiles/TileContentLoader.h"
#include "Cesium3DTiles/TileID.h"
#include "Cesium3DTiles/TileRefine.h"
#include "CesiumGltf/GltfReader.h"
#include <cstddef>
#include <glm/mat4x4.hpp>
#include <gsl/span>
#include <spdlog/fwd.h>

namespace Cesium3DTiles {

class Tileset;

/**
 * @brief Creates {@link TileContentLoadResult} from glTF data.
 */
class CESIUM3DTILES_API GltfContent final : public TileContentLoader {
public:
  /**
   * @copydoc TileContentLoader::load
   *
   * The result will only contain the `model`. Other fields will be
   * empty or have default values.
   */
  std::unique_ptr<TileContentLoadResult>
  load(const TileContentLoadInput& input) override;

  /**
   * @brief Create a {@link TileContentLoadResult} from the given data.
   *
   * (Only public to be called from `Batched3DModelContent`)
   *
   * @param pLogger Only used for logging
   * @param url The URL, only used for logging
   * @param data The actual glTF data
   * @return The {@link TileContentLoadResult}
   */
  static std::unique_ptr<TileContentLoadResult> load(
      const std::shared_ptr<spdlog::logger>& pLogger,
      const std::string& url,
      const gsl::span<const std::byte>& data);

  /**
   * @brief Creates texture coordinates for raster tiles that are mapped to 3D
   * tiles.
   *
   * This is not supposed to be called by clients.
   *
   * It will be called for all {@link RasterMappedTo3DTile} objects of a
   * {@link Tile}, and extend the accessors of the given glTF model with
   * accessors that contain the texture coordinate sets for different
   * projections. Further details are not specified here.
   *
   * @param gltf The glTF model.
   * @param textureCoordinateID The texture coordinate ID.
   * @param projection The {@link CesiumGeospatial::Projection}.
   * @param rectangle The {@link CesiumGeometry::Rectangle}.
   * @return The bounding region.
   */
  static CesiumGeospatial::BoundingRegion createRasterOverlayTextureCoordinates(
      CesiumGltf::Model& gltf,
      uint32_t textureCoordinateID,
      const CesiumGeospatial::Projection& projection,
      const CesiumGeometry::Rectangle& rectangle);

private:
  static CesiumGltf::GltfReader _gltfReader;
};

} // namespace Cesium3DTiles

ext: h


text:
#ifndef IOSLICEMUT_H
#define IOSLICEMUT_H

#include <WinSock2.h>

#include <variant>
#include <gsl/span>

namespace laio {

    template<typename T>
    using Result = std::variant<T, std::exception>;

    namespace net {

        class IoSpanMut {
            WSABUF _raw_wsa_buffer;
        public:
            explicit constexpr IoSpanMut(const WSABUF& wsaBuffer) noexcept
                : _raw_wsa_buffer{wsaBuffer} {}

            explicit constexpr IoSpanMut(WSABUF&& wsaBuffer) noexcept
                : _raw_wsa_buffer{std::move(wsaBuffer)} {}  // NOLINT(hicpp-move-const-arg,performance-move-const-arg)

            static inline IoSpanMut create(unsigned char buf[]) noexcept {
                static_assert(sizeof(*buf) <= (std::numeric_limits<ULONG>::max)());

                // Seriously? Casting an unsigned char to a signed char?
                // Well behaved for all 128 ASCII characters - Blows up for any value lesser than 0 or greater than 127
                return IoSpanMut{std::move(WSABUF{  // NOLINT(hicpp-move-const-arg,performance-move-const-arg)
                    sizeof(*buf),
                    reinterpret_cast<CHAR*>(buf),
                })};
            }

            inline Result<std::monostate> advance(std::size_t n) noexcept {
                if (static_cast<std::size_t>(sizeof(_raw_wsa_buffer.len)) < n) {
                    return std::out_of_range("Advancing out of range");
                }
                _raw_wsa_buffer.len -= static_cast<ULONG>(n);
                _raw_wsa_buffer.buf += n;
                return std::monostate{};
            }

            inline gsl::span<const unsigned char> as_span() noexcept {
                return gsl::span<const unsigned char>{
                    reinterpret_cast<const unsigned char*>(_raw_wsa_buffer.buf),
                    static_cast<int>(_raw_wsa_buffer.len)
                };
            }

            inline gsl::span<unsigned char> as_mut_span() noexcept {
                return gsl::span<unsigned char>{
                    reinterpret_cast<unsigned char*>(_raw_wsa_buffer.buf),
                    static_cast<int>(_raw_wsa_buffer.len)
                };
            }
        };

    } // namespace net

} // namespace laio

#endif // IOSLICEMUT_H

max_stars_repo_name: mjptree/laio

max_issues_repo_name: mjptree/laio

max_forks_repo_name: mjptree/laio


text:
#define _XOPEN_SOURCE 600

#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#include <stdint.h>
#include <time.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <errno.h>
#include "fileio.h"
#include "ipf_balance.h"
#include "likelihood.h"
#include "master_tasks.h"
#include "mutation.h"
#include "randomisation.h"
#include "resampling.h"
#include "typedefs.h"

static const int ALG = 5;
static const int MASTER = 0;

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

  struct Process process;
  struct Estimates estimates;
  long * inds;
  long * BM;
  int do_resampling;
  double * w;
  double * x;
  double * x_resamp;
  double * proc_w;
  double * output;
  double my_normal_weight;
  double rsamp_time = 0;
  double rsamp_start = 0;
  double comm_time = 0;
  double comm_start;
  double comm_time2 = 0;
  double total_w = 1.0;
  double weight_time = 0.0;
  double mutation_time = 0.0;
  double init_time = 0.0;
  long n_xfer;
  unsigned int *cnts;
  long *cnts_l;
  const gsl_rng_type *T;
  double timer_start;

  // Parse command line arguments
  int run_idx = 0;
  if(argc > 1)
    run_idx = atoi(argv[1]);

  long N = 1000;
  if(argc > 2)
    N = atol(argv[2]);

  int state_dim;
  if ( argc > 3 )
    state_dim = atoi( argv[ 3 ] );
  else {
    perror("State dimesion must be provided");
    return -1;
  }

  double epc_threshold;
  if (argc > 4 )
    epc_threshold = atof( argv[ 4 ] );
  else {
    perror("Reampling threshold must be provided");
    return -1;
  }

  // Initialize the MPI environment
  MPI_Init(NULL, NULL);

  gsl_rng_env_setup();
  T = gsl_rng_default;
  gsl_rng *rgen;
  rgen = gsl_rng_alloc( T );

  // Get the number of processes
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);

  // Get the rank of the process
  int world_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

  // Start the timer
  struct timespec start_time;;
  clock_gettime(CLOCK_REALTIME, &start_time);

  MPI_Barrier( MPI_COMM_WORLD );
  timer_start = MPI_Wtime();

  long data_length;
  process = readProcessFile( state_dim );
  data_length = (int)process.n;

  //---------------
  // Initialisation
  //---------------
  x = malloc( sizeof( double ) * N * state_dim);
  w = malloc( sizeof( double ) * N );

  inds = malloc( sizeof( long ) * N );
  x_resamp = malloc( sizeof( double ) * N * state_dim);

  proc_w = malloc( sizeof( double ) * world_size );
  cnts = malloc( sizeof( unsigned int ) * world_size );
  cnts_l = malloc( sizeof( long ) * world_size );

  BM = malloc( sizeof( long ) * world_size * world_size );
  initialise_balancing( world_size );

  // Initial distribution setup
  setupInitialSeed( world_rank, rgen );

  // Create initial sample
  createInitialSample( x, N, rgen, state_dim );

  // Setup the randomisation
  set_resamplingseed( world_rank );

  int output_dim = 2 + state_dim + 2 * world_size + world_size * state_dim;
  if ( world_rank == MASTER )
    output = malloc( sizeof( double ) * output_dim * data_length );

  MPI_Barrier( MPI_COMM_WORLD );
  init_time += MPI_Wtime() - timer_start;

  my_normal_weight = 1.0 / world_size;

  // main loop
  for( long n = 0; n < data_length; n++ ) {

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    timer_start = MPI_Wtime();

    evaluateLikelihoods( x, process, state_dim, data_length, N, n, w,
			 my_normal_weight );

    MPI_Barrier( MPI_COMM_WORLD );
    weight_time += MPI_Wtime() - timer_start ;

    // compute the relevant estimates
    estimates = estimateOutput( w, x, world_rank, N, run_idx, state_dim );

    // Combine the results of all filters
    combineEstimates( world_rank, state_dim, estimates, output, n,
		      output_dim, world_size, &total_w, &my_normal_weight,
		      &do_resampling, epc_threshold );
    //    printf("My weight %i: %f\n", (int)world_rank,my_normal_weight);
    // -----------------------------
    // RESAMPLING
    // -----------------------------

    // Resampling within processors

    MPI_Barrier( MPI_COMM_WORLD );
    rsamp_start = MPI_Wtime();

    serial_multinomial_resample_sng_noutk( N, w, inds, N );

    for ( long i = 0; i < N; i++)
      for ( int d = 0; d < state_dim; d++ )
	x_resamp[ state_dim * i + d ] = x[ state_dim * inds[ i ] + d ];

    for ( long i = 0; i < N; i++)
      for ( int d = 0; d < state_dim ; d++ )
	x[ state_dim * i + d ] = x_resamp[ state_dim * i + d ];

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    rsamp_time += MPI_Wtime() - rsamp_start;
    comm_start = MPI_Wtime();

    if ( do_resampling ) {

      //my_weight = estimates.w; // sum of weights for the current process

      // --------------------------
      // START THE FULL INTERACTION
      // --------------------------

      // Processes send their weights to the MASTER process
      MPI_Gather( &my_normal_weight, 1, MPI_DOUBLE, proc_w, 1, MPI_DOUBLE, MASTER,
		  MPI_COMM_WORLD );
      MPI_Barrier( MPI_COMM_WORLD );

      // Sample how many particles need to be taken from each process
      if ( world_rank == MASTER ) {

	gsl_ran_multinomial( rgen, ( size_t ) world_size,
			     ( unsigned int ) ( world_size ), proc_w, cnts );

	for( int j = 0; j < world_size; j++ )
	  cnts_l[ j ] = ( long ) cnts[ j ];

      }

      // Send the island resampled index array to all processes
      MPI_Bcast( cnts_l, world_size, MPI_LONG, MASTER, MPI_COMM_WORLD );

      //if ( world_rank == MASTER )
      //  for ( int j = 0; j < world_size; j++ )
      //    output[ output_dim * n + count_offset + j ] = cnts_l[ j ];

      // -----------------------------------------------------------------------------
      MPI_Barrier( MPI_COMM_WORLD );
      comm_time2 += MPI_Wtime() - comm_start;
      comm_start = MPI_Wtime();

      // Calculate how the islands should be communicated
      calculate_balancing_matrix( cnts_l, world_size , BM , 1 );

      if ( cnts_l[ world_rank ] < 1 ) { // this process needs an island

	// iterate over all processes (other than the current) and see if that
	// process has surplus particles
	for ( int i = 0; i < world_size; i++ ) {

	  if ( i != world_rank ) {

	    // NB! BM[ world_size * i + world_rank ] should be 0 or 1
	    n_xfer = BM[ world_size * i + world_rank ] * N;

	    if ( n_xfer > 0 ) {

	      MPI_Recv( x, N * state_dim, MPI_DOUBLE, i, 2, MPI_COMM_WORLD,
			MPI_STATUS_IGNORE );

	      cnts_l[ world_rank ] += n_xfer;

	    }
	  }
	}

      } else if ( cnts_l[ world_rank ] > 1 ) { // this process has excess islands

	// iterate over all processes (other than the current) and see if
	// they need particles
	for (int i = 0; i < world_size; i++ ) {
	  if( world_rank != i ) {

	    // NB! BM[ world_size * i + world_rank ] should be 0 or 1
	    n_xfer = BM[ world_size * world_rank + i] * N;

	    if ( n_xfer > 0 ) {

	      MPI_Send( x, n_xfer * state_dim, MPI_DOUBLE, i, 2, MPI_COMM_WORLD );

	      cnts_l[ world_rank ] -= n_xfer;

	    }
	  }
	}

      }

      my_normal_weight = 1.0 / world_size;
    }
    //else {
    //      printf("Skip resample!\n");
    //    }

    //if ( world_rank == MASTER ) {
    //  printf("%i: %f\n", (int)n,my_normal_weight);
    //}

    // ----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    comm_time += MPI_Wtime() - comm_start;
    timer_start = MPI_Wtime();

    mutation(N, x, state_dim, rgen );

    // -----------------------------------------------------------------------------
    MPI_Barrier( MPI_COMM_WORLD );
    mutation_time += MPI_Wtime() - timer_start ;

  }

  struct timespec end_time;
  clock_gettime(CLOCK_REALTIME, &end_time);

  createOutput( ALG, world_rank, output, data_length, world_size,
		run_idx, N, state_dim, start_time, end_time, comm_time, rsamp_time,
		comm_time2, weight_time, mutation_time, init_time,
		(int)( atof( argv[ 4 ] ) * 100 ) );

  finalise_balancing();

  // Finalize the MPI environment.
  MPI_Finalize();

  return 0;
}

ext: c

max_stars_repo_name: heinekmp/AIRPF

max_issues_repo_name: heinekmp/AIRPF

max_forks_repo_name: heinekmp/AIRPF


text:
/* integration/workspace.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

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

gsl_integration_workspace *
gsl_integration_workspace_alloc (const size_t n)
{
  gsl_integration_workspace * w ;

  if (n == 0)
    {
      GSL_ERROR_VAL ("workspace length n must be positive integer",
			GSL_EDOM, 0);
    }

  w = (gsl_integration_workspace *)
    malloc (sizeof (gsl_integration_workspace));

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

  w->alist = (double *) malloc (n * sizeof (double));

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

      GSL_ERROR_VAL ("failed to allocate space for alist ranges",
			GSL_ENOMEM, 0);
    }

  w->blist = (double *) malloc (n * sizeof (double));

  if (w->blist == 0)
    {
      free (w->alist);
      free (w);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for blist ranges",
			GSL_ENOMEM, 0);
    }

  w->rlist = (double *) malloc (n * sizeof (double));

  if (w->rlist == 0)
    {
      free (w->blist);
      free (w->alist);
      free (w);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for rlist ranges",
			GSL_ENOMEM, 0);
    }

  w->elist = (double *) malloc (n * sizeof (double));

  if (w->elist == 0)
    {
      free (w->rlist);
      free (w->blist);
      free (w->alist);
      free (w);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for elist ranges",
			GSL_ENOMEM, 0);
    }

  w->order = (size_t *) malloc (n * sizeof (size_t));

  if (w->order == 0)
    {
      free (w->elist);
      free (w->rlist);
      free (w->blist);
      free (w->alist);
      free (w);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for order ranges",
			GSL_ENOMEM, 0);
    }

  w->level = (size_t *) malloc (n * sizeof (size_t));

  if (w->level == 0)
    {
      free (w->order);
      free (w->elist);
      free (w->rlist);
      free (w->blist);
      free (w->alist);
      free (w);		/* exception in constructor, avoid memory leak */

      GSL_ERROR_VAL ("failed to allocate space for order ranges",
			GSL_ENOMEM, 0);
    }

  w->size = 0 ;
  w->limit = n ;
  w->maximum_level = 0 ;

  return w ;
}

void
gsl_integration_workspace_free (gsl_integration_workspace * w)
{
  free (w->level) ;
  free (w->order) ;
  free (w->elist) ;
  free (w->rlist) ;
  free (w->blist) ;
  free (w->alist) ;
  free (w) ;
}

/*
size_t
gsl_integration_workspace_limit (gsl_integration_workspace * w)
{
  return w->limit ;
}

size_t
gsl_integration_workspace_size (gsl_integration_workspace * w)
{
  return w->size ;
}
*/


text:
/**
 * Interface of spce_pathlength.h
 */

#ifndef _SPCE_PATHLENGTH_H

#define _SPCE_PATHLENGTH_H

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include "spc_trace_functions.h"

extern int
abscissa_to_pathlength (const trace_func * const func, gsl_vector * const data);

#endif /* !_SPCE_PATHLENGTH_H */

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*                  This file is part of the program                         */
/*          GCG --- Generic Column Generation                                */
/*                  a Dantzig-Wolfe decomposition based extension            */
/*                  of the branch-cut-and-price framework                    */
/*         SCIP --- Solving Constraint Integer Programs                      */
/*                                                                           */
/* Copyright (C) 2010-2019 Operations Research, RWTH Aachen University       */
/*                         Zuse Institute Berlin (ZIB)                       */
/*                                                                           */
/* This program 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.                    */
/*                                                                           */
/* 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 Lesser General Public License for more details.                       */
/*                                                                           */
/* You should have received a copy of the GNU Lesser General Public License  */
/* along with this program; if not, write to the Free Software               */
/* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.*/
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

/**@file   graph_gcg.h
 * @brief  Implementation of the graph which supports both node and edge weights.
 * @author Igor Pesic
 */

/*---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/

#ifndef GCG_GRAPH_GCG_H_
#define GCG_GRAPH_GCG_H_

#include <map>
#include "bridge.h"

#ifdef WITH_GSL
   #include <gsl/gsl_spmatrix.h>
   #include <gsl/gsl_spblas.h>
#endif

namespace gcg {

class EdgeGCG
{
public:
    int src, dest;
    double weight;
    EdgeGCG() {}
    EdgeGCG(int s, int d, double w): src(s), dest(d), weight(w) {}
} ;

class GraphGCG: public gcg::Bridge
{
private:
   bool undirected;
   bool locked;           // true if we are not allowed to change the graph anymore
   bool initialized;      // true if at least 1 node
   std::vector<int> nodes;
#ifdef WITH_GSL
   gsl_spmatrix* adj_matrix_sparse;
   gsl_spmatrix* working_adj_matrix;      // this one is used ONLY during MCL algorithm!
#else
   std::vector<std::vector<double>> adj_matrix;      /** For undirected graphs, this matrix is symmetrical */
#endif
   std::vector<EdgeGCG*> edges;

public:

   GraphGCG();
   GraphGCG(int _n_nodes, bool _undirected);

   virtual ~GraphGCG();
   virtual SCIP_RETCODE addNNodes(int _n_nodes);
   virtual SCIP_RETCODE addNNodes(int _n_nodes, std::vector<int> weights);
   virtual int getNNodes();
   virtual int getNEdges();

#ifdef WITH_GSL
   virtual gsl_spmatrix* getAdjMatrix();
   virtual void expand(int factor);
   virtual void inflate(double factor);
   virtual void colL1Norm();
   virtual void prune();
   virtual bool stopMCL(int iter);
   virtual std::vector<int> getClustersMCL();
   virtual void initMCL();
   virtual void clearMCL();
#else
   virtual std::vector<std::vector<double>> getAdjMatrix();
#endif
   virtual SCIP_RETCODE getEdges(std::vector<void*>& edges);
   virtual SCIP_Bool isEdge(int node_i, int node_j);
   virtual int getNNeighbors(int node);
   virtual std::vector<int> getNeighbors(int node);
   virtual std::vector<std::pair<int, double> > getNeighborWeights(int node);
   virtual SCIP_RETCODE addNode(int node, int weight);
   virtual SCIP_RETCODE addNode();                    /** Sets node weight to 0 and the ID to the next available. */
   virtual SCIP_RETCODE deleteNode(int node);
   virtual SCIP_RETCODE addEdge(int i, int j);        /** Sets edge weight to 1. */
   virtual SCIP_RETCODE addEdge(int node_i, int node_j, double weight);
   virtual SCIP_RETCODE setEdge(int node_i, int node_j, double weight);
   virtual SCIP_RETCODE deleteEdge(int node_i, int node_j);
   virtual int graphGetWeights(int node);
   virtual double getEdgeWeight(int node_i, int node_j);

   virtual int edgeComp(const EdgeGCG* a, const EdgeGCG* b);

   virtual SCIP_RETCODE flush();       // lock the graph and compresses the adj matrix if we use GSL
   virtual SCIP_RETCODE normalize();
   virtual double getEdgeWeightPercentile(double q);
};

} /* namespace gcg */
#endif /* GCG_GRAPH_GCG_H_ */


text:
/*
MPCOTool:
The Multi-Purposes Calibration and Optimization Tool. A software to perform
calibrations or optimizations of empirical parameters.

AUTHORS: Javier Burguete and Borja Latorre.

Copyright 2012-2019, AUTHORS.

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 AUTHORS ``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 AUTHORS 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 main.c
 * \brief Main source file.
 * \authors Javier Burguete and Borja Latorre.
 * \copyright Copyright 2012-2019, all rights reserved.
 */
#define _GNU_SOURCE
#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <locale.h>
#include <gsl/gsl_rng.h>
#include <libxml/parser.h>
#include <libintl.h>
#include <glib.h>
#include <json-glib/json-glib.h>
#ifdef G_OS_WIN32
#include <windows.h>
#endif
#if HAVE_MPI
#include <mpi.h>
#endif
#if HAVE_GTK
#include <gio/gio.h>
#include <gtk/gtk.h>
#endif
#include "genetic/genetic.h"
#include "utils.h"
#include "experiment.h"
#include "variable.h"
#include "input.h"
#include "optimize.h"
#if HAVE_GTK
#include "interface.h"
#endif
#include "mpcotool.h"

int
main (int argn, char **argc)
{
#if HAVE_GTK
  show_pending = process_pending;
#endif
  return mpcotool (argn, argc);
}

max_stars_repo_name: jburguete/mpcotool

max_issues_repo_name: jburguete/mpcotool

max_forks_repo_name: jburguete/mpcotool


text:
/**
 * \file DelayInterface.h
 */

#ifndef ATK_DELAY_DELAYINTERFACE_H
#define ATK_DELAY_DELAYINTERFACE_H

#include <ATK/Delay/config.h>

#include <gsl/gsl>

#include <map>

namespace ATK
{
  /// Interface for a fixed filter
  class ATK_DELAY_EXPORT DelayInterface
  {
  public:
    virtual ~DelayInterface() = default;
    /// Sets the delay of the filter
  virtual void set_delay(gsl::index delay) = 0;
    /// Returns the delay
    virtual gsl::index get_delay() const = 0;
  };

  /// Interface for a universal filter
  template<class DataType>
  class ATK_DELAY_EXPORT UniversalDelayInterface
  {
  public:
    virtual ~UniversalDelayInterface() = default;

    /// Sets the blend of the filter
    virtual void set_blend(DataType blend) = 0;
    /// Returns the blend
    virtual DataType get_blend() const = 0;
    /// Sets the feedback of the filter
    virtual void set_feedback(DataType feedback) = 0;
    /// Returns the feedback
    virtual DataType get_feedback() const = 0;
    /// Sets the feedforward of the filter
    virtual void set_feedforward(DataType feedforward) = 0;
    /// Returns the feedforward
    virtual DataType get_feedforward() const = 0;
  };
}

#endif

max_stars_repo_name: AudioTK/AudioTK

max_issues_repo_name: AudioTK/AudioTK

max_forks_repo_name: AudioTK/AudioTK


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

#include <math.h>
#include <libgen.h>

#include <options/options.h>
#include "../csm/csm_all.h"

void purify(LDP ld, double threshold_min, double threshold_max);

struct ld_purify_params {
	double threshold_min, threshold_max;
	const char* file_input;
	const char* file_output;
};

int main(int argc, const char * argv[]) {
	sm_set_program_name(argv[0]);

	options_banner("ld_purify: Makes sure that the file format is valid. \n * Sets valid=0 if reading is outside interval ");

	struct ld_purify_params p;

	struct option* ops = options_allocate(20);
	options_double(ops, "threshold_min", &p.threshold_min, 0.01,
		"Sets valid=0 if readings are less than this threshold.");
	options_double(ops, "threshold_max", &p.threshold_max, 79.0,
		"Sets valid=0 if readings are more than this threshold.");

	options_string(ops, "in", &p.file_input, "stdin", "Input file ");
	options_string(ops, "out", &p.file_output, "stdout", "Output file ");

	if(!options_parse_args(ops, argc, argv)) {
		options_print_help(ops, stderr);
		return -1;
	}

	FILE * in = open_file_for_reading(p.file_input);
	if(!in) return -3;

	FILE * out = open_file_for_writing(p.file_output);
	if(!out) return -2;

	LDP ld; int count = -1;
	while( (ld = ld_from_json_stream(in))) {

		purify(ld, p.threshold_min, p.threshold_max);

		if(!ld_valid_fields(ld))  {
			sm_error("Wait, we didn't purify enough  (#%d in file)\n", count);
			continue;
		}

		ld_write_as_json(ld, out);
		ld_free(ld);
	}

	return 0;
}

void purify(LDP ld, double threshold_min, double threshold_max) {
	for(int i=0;i<ld->nrays;i++) {
		if(!ld->valid[i]) continue;

		double rho = ld->readings[i];
		if( is_nan(rho) | (rho < threshold_min) | (rho > threshold_max) ) {
			ld->readings[i] = GSL_NAN;
			ld->valid[i] = 0;
			ld->alpha[i] = GSL_NAN;
			ld->alpha_valid[i] = 0;
		}

	}

}

max_stars_repo_name: ozaslan/basics

max_issues_repo_name: ozaslan/basics

max_forks_repo_name: ozaslan/basics


text:
#include <algorithm>
#include <numeric>

#ifdef HAS_MKL
#include <mkl.h>
#else
#include <cblas.h>
#endif

#ifdef __SSE4_2__
#include <emmintrin.h>
#endif

typedef double Real;

inline void GetTotalTraction_(Real *TotalTraction, const Real *CauchyStressTensor,
    const Real *ElectricDisplacementx, int ndim) {
    if (ndim==3) {
        TotalTraction[0] = CauchyStressTensor[0];
        TotalTraction[1] = CauchyStressTensor[4];
        TotalTraction[2] = CauchyStressTensor[8];
        TotalTraction[3] = CauchyStressTensor[1];
        TotalTraction[4] = CauchyStressTensor[2];
        TotalTraction[5] = CauchyStressTensor[5];

        TotalTraction[6] = ElectricDisplacementx[0];
        TotalTraction[7] = ElectricDisplacementx[1];
        TotalTraction[8] = ElectricDisplacementx[2];
    }
    else if (ndim == 2) {
        TotalTraction[0] = CauchyStressTensor[0];
        TotalTraction[1] = CauchyStressTensor[3];
        TotalTraction[2] = CauchyStressTensor[1];

        TotalTraction[3] = ElectricDisplacementx[0];
        TotalTraction[4] = ElectricDisplacementx[1];
    }
}

inline void FillConstitutiveB_(Real *B, const Real* SpatialGradient,
                     int ndim, int nvar, int rows, int cols) {
    int i = 0;

    if (ndim == 3) {

        for (; i<rows; ++i) {

            // Store in registers
            const Real a0 = SpatialGradient[i*ndim];
            const Real a1 = SpatialGradient[i*ndim+1];
            const Real a2 = SpatialGradient[i*ndim+2];

            // MECHANICAL TERMS
            B[i*cols*nvar]              = a0;
            B[i*cols*nvar+cols+1]       = a1;
            B[i*cols*nvar+2*(cols+1)]   = a2;

            B[i*cols*nvar+cols+5]       = a2;
            B[i*cols*nvar+2*cols+5]     = a1;

            B[i*cols*nvar+4]            = a2;
            B[i*cols*nvar+2*cols+4]     = a0;

            B[i*cols*nvar+3]            = a1;
            B[i*cols*nvar+cols+3]       = a0;

            // ELECTROSTATIC TERMS
            B[i*cols*nvar+3*cols+6]     = a0;
            B[i*cols*nvar+3*cols+7]     = a1;
            B[i*cols*nvar+3*cols+8]     = a2;
        }
    }

    else if (ndim == 2) {

        for (; i<rows; ++i) {

            // Store in registers
            const Real a0 = SpatialGradient[i*ndim];
            const Real a1 = SpatialGradient[i*ndim+1];

            // MECHANICAL TERMS
            B[i*cols*nvar]              = a0;
            B[i*cols*nvar+cols+1]       = a1;

            B[i*cols*nvar+2]            = a1;
            B[i*cols*nvar+cols+2]       = a0;

            // ELECTROSTATIC TERMS
            B[i*cols*nvar+2*cols+3]     = a0;
            B[i*cols*nvar+2*cols+4]     = a1;
        }
    }
}
inline void _ConstitutiveStiffnessIntegrandDPF_Filler_(Real *stiffness, Real *traction,
    const Real* SpatialGradient,
    const Real* ElectricDisplacementx,
    const Real* CauchyStressTensor,
    const Real* H_Voigt,
    const Real* detJ,
    int ngauss,
    int noderpelem,
    int ndim,
    int nvar,
    int H_VoigtSize,
    int requires_geometry_update) {

    int local_size = nvar*noderpelem;

    Real *t;

#ifdef __SSE4_2__
    if (ndim==3) {
        t = (Real*)_mm_malloc(9*sizeof(Real),32);
    }
    else if (ndim==2) {
        t = (Real*)_mm_malloc(5*sizeof(Real),32);
    }

    Real *B = (Real*)_mm_malloc(H_VoigtSize*local_size*sizeof(Real),32);
    Real *HBT = (Real*)_mm_malloc(H_VoigtSize*local_size*sizeof(Real),32);
    Real *BDB_1 = (Real*)_mm_malloc(local_size*local_size*sizeof(Real),32);
#else
    if (ndim==3) {
        t = (Real*)malloc(9*sizeof(Real));
    }
    else if (ndim==2) {
        t = (Real*)malloc(5*sizeof(Real));
    }
    // Real *local_traction = (Real*)malloc(local_size*sizeof(Real));
    // std::fill(local_traction,local_traction+local_size,0.);

    Real *B = (Real*)malloc(H_VoigtSize*local_size*sizeof(Real));
    Real *HBT = (Real*)malloc(H_VoigtSize*local_size*sizeof(Real));
    Real *BDB_1 = (Real*)malloc(local_size*local_size*sizeof(Real));
#endif

    std::fill(B,B+H_VoigtSize*local_size,0.);
    // std::fill(HBT,HBT+H_VoigtSize*local_size,0.);

    for (int igauss = 0; igauss < ngauss; ++igauss) {

        FillConstitutiveB_(B,&SpatialGradient[igauss*ndim*noderpelem],ndim,nvar,noderpelem,H_VoigtSize);

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
            H_VoigtSize, local_size, H_VoigtSize, 1.0, &H_Voigt[igauss*H_VoigtSize*H_VoigtSize], H_VoigtSize, B, H_VoigtSize, 0.0, HBT, local_size);

        cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
            local_size, local_size, H_VoigtSize, 1.0, B, H_VoigtSize, HBT, local_size, 0.0, BDB_1, local_size);

        // Multiply stiffness with detJ
        const Real detJ_igauss = detJ[igauss];
        for (int i=0; i<local_size*local_size; ++i) {
            stiffness[i] += BDB_1[i]*detJ_igauss;
        }

        if (requires_geometry_update==1) {
            // Compute tractions
            GetTotalTraction_(t,
                &CauchyStressTensor[igauss*ndim*ndim],
                &ElectricDisplacementx[igauss*ndim], ndim);

            // Multiply B with traction - for loop is okay
            // std::fill(local_traction,local_traction+local_size,0.);
            for (int i=0; i<local_size; ++i) {
                Real tmp = 0;
                for (int j=0; j<H_VoigtSize; ++j) {
                    tmp += B[i*H_VoigtSize+j]*t[j];
                }
                // local_traction[i] = tmp;
                traction[i] += tmp*detJ_igauss;
            }

            // // Multiply traction with detJ
            // for (int i=0; i<local_size; ++i) {
            //     traction[i] += local_traction[i]*detJ_igauss;
            // }
        }
    }

#ifdef __SSE4_2__
    _mm_free(t);
    _mm_free(B);
    _mm_free(HBT);
    _mm_free(BDB_1);
#else
    free(t);
    // free(local_traction);
    free(B);
    free(HBT);
    free(BDB_1);
#endif

    // Real AA[12]; // 3x4
    // Real BB[12]; // 4x3
    // Real CC[12];
    // std::iota(AA,AA+12,1);
    // std::iota(BB,BB+12,2);

    // cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
    //                 3, 3, 4, 1.0, AA, 4, BB, 3, 0.0, CC, 3);

    // for (int i=0; i<9; ++i) {
    //     // std::cout << AA[i] << " " << BB[i] << '\n';
    //     std::cout << CC[i] << '\n';
    // }
    // std::cout << '\n';
}

max_stars_repo_name: jdlaubrie/florence

max_issues_repo_name: jdlaubrie/florence

max_forks_repo_name: jdlaubrie/florence


text:
#ifndef _GALAXY_H_
#define _GALAXY_H_

#include <vector>
#include <gsl/gsl_integration.h>
#include "cosmology.h"
#include <vector_types.h>

class galaxy{
    double ra, dec, red, w, nbar;
    double3 cart;
    bool cart_set;

    public:
        galaxy(double RA, double DEC, double RED, double NZ, double W);

        void bin(std::vector<double> &delta, int3 N, double3 L, double3 r_min,
                 cosmology &cosmo, double3 &pk_nbw, double3 &bk_nbw,
                 gsl_integration_workspace *ws);

        void set_cartesian(cosmology &cosmo, double3 r_min, gsl_integration_workspace *ws);

        double3 get_unshifted_cart(cosmology &cosmo, gsl_integration_workspace *ws);

        double3 get_cart();
};

#endif


text:
/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * 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 or any later version.
 *
 * 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 program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */

/*****************************************************************************
 *                                  dgt_usefftw.h
 *
 * Discrete Gabor Transform by Using FFTW.
 *
 * These routines are designed for calculating discrete Gabor transform and
 * its inversion of 1D signals. In order to eliminate the border effect, the
 * input signal("signal") is extended by three forms: zeros padded("zpd"),
 * periodized extension("ppd") and symetric extension("sym").
 *
 * The analysis/synthesis function is given by users, and it's daul
 * (synthesis/analysis) function can be computed by "daul" routine. The over
 * sampling rate is equal to N/dM, where N denotes frequency sampling numbers
 * and dM denotes the time sampling interval.
 *
 * N and dM should can be devided evenly by the window length "Lw". The
 * recovered signal just has the elements from 1 to dM*floor(Ls/dM) of the
 * original signal. So you'd better let dM can be deviede evenly by  the
 * original signal length "Ls".
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/

#ifndef DGT_USEFFTW_H
#define DGT_USEFFTW_H

#include <string>
#include <complex>
#include <fftw.h>
#include <matrix.h>
#include <linequs1.h>
//#include <linequs3.h>
#include <utilities.h>

namespace splab
{

    template<typename Type> Vector<Type>
    daulFFTW( const Vector<Type>&, int, int );

    template<typename Type>
    Matrix< complex<Type> > dgtFFTW( const Vector<Type>&,
                                     const Vector<Type>&,
                                     int, int,
                                     const string &mode = "zpd" );
    template<typename Type>
    Vector<Type> idgtFFTW( const Matrix< complex<Type> >&,
                           const Vector<Type>&,
                           int, int );

    #include <dgt_usefftw-impl.h>

}
// namespace splab

#endif
// DGT_USEFFTW_H

max_stars_repo_name: cianfrocco-lab/Motion-correction

max_issues_repo_name: wjiang/motioncorr

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: wjiang/motioncorr

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
/* bst/gsl_bst_rb.h
 *
 * Copyright (C) 2018 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.
 */

#ifndef __GSL_BST_RB_H__
#define __GSL_BST_RB_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_bst_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

#ifndef GSL_BST_RB_MAX_HEIGHT
#define GSL_BST_RB_MAX_HEIGHT 48
#endif

/* red-black node */
struct gsl_bst_rb_node
{
  struct gsl_bst_rb_node *rb_link[2]; /* subtrees */
  void *rb_data;                      /* pointer to data */
  unsigned char rb_color;             /* color */
};

/* red-black tree data structure */
typedef struct
{
  struct gsl_bst_rb_node *rb_root;   /* tree's root */
  gsl_bst_cmp_function *rb_compare;  /* comparison function */
  void *rb_param;                    /* extra argument to |rb_compare| */
  const gsl_bst_allocator *rb_alloc; /* memory allocator */
  size_t rb_count;                   /* number of items in tree */
  unsigned long rb_generation;       /* generation number */
} gsl_bst_rb_table;

/* red-black traverser structure */
typedef struct
{
  const gsl_bst_rb_table *rb_table;  /* tree being traversed */
  struct gsl_bst_rb_node *rb_node;   /* current node in tree */
  struct gsl_bst_rb_node *rb_stack[GSL_BST_RB_MAX_HEIGHT];
                                     /* all the nodes above |rb_node| */
  size_t rb_height;                  /* number of nodes in |rb_parent| */
  unsigned long rb_generation;       /* generation number */
} gsl_bst_rb_traverser;

__END_DECLS

#endif /* __GSL_BST_RB_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/*
** read nuisance covariates file into a matrix
**
** G.Lohmann, MPI-KYB, 2018
*/

#include <viaio/Vlib.h>
#include <viaio/mu.h>

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>

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

extern int VStringToken (char *,char *,int,int);
extern int line_empty(char *buf,int len);
extern int CheckBuffer(char *buf,int len);
extern int VistaFormat(char *buf,int len);
extern int test_ascii(int val);

gsl_matrix *VReadCovariates(VString cfile,VBoolean demean)
{
  float x=0;
  int i=0,j=0,nrows=0,ncols=0,len=10000,tlen=80;
  char *buf = (char *)VCalloc(len,sizeof(char));
  char *token = (char *)VCalloc(tlen,sizeof(char));

  FILE *fp = fopen(cfile,"r");
  if (fp == NULL) VError(" err opening covariates file %s",cfile);

  /* get matrix dimensions */
  while (!feof(fp)) {

    for (j=0; j<len; j++) buf[j] = '\0';
    if (fgets(buf,len,fp) == NULL) break;
    if (VistaFormat(buf,len) > 1) VError(" File %s must be a text file",cfile);
    if (CheckBuffer(buf,len) < 1) continue;
    if (! test_ascii((int)buf[0])) VError(" File %s: line %d begins with an illegal character (%c)",cfile,nrows,buf[0]);

    j = 0;
    while (VStringToken(buf,token,j,tlen)) {
      sscanf(token,"%f",&x);
      j++;
    }
    if (j > ncols) ncols = j;
    nrows++;
  }
  rewind(fp);
  fprintf(stderr," Nuisance covariates:  %d x %d\n",nrows,ncols);
  if (ncols < 1) VError(" Nuisance covariates: no columns in file %s",cfile);

  /* read into matrix */
  gsl_matrix *dest = gsl_matrix_calloc(nrows,ncols);
  i=0;
  while (!feof(fp)) {
    for (j=0; j<len; j++) buf[j] = '\0';
    if (fgets(buf,len,fp) == NULL) break;
    if (CheckBuffer(buf,len) < 1) continue;

    j = 0;
    while (VStringToken(buf,token,j,tlen)) {
      if (j >= ncols) break;
      sscanf(token,"%f",&x);
      gsl_matrix_set(dest,i,j,(double)x);
      j++;
    }
    i++;
    if (i >= nrows) break;
  }
  if (!demean) return dest;

  /* subtract mean */
  double u=0,sum1=0,sum2=0,nx=0,mean=0,sigma=0,tiny=1.0e-6;
  for (i=0; i<dest->size2; i++) {
    sum1 = sum2 = nx = 0;
    for (j=0; j<dest->size1; j++) {
      u = gsl_matrix_get(dest,j,i);
      sum1 += u;
      sum2 += u*u;
      nx++;
    }
    mean = sum1/nx;
    sigma = sqrt((sum2 - nx * mean * mean) / (nx - 1.0));
    if (sigma < tiny) VError(" Nuisance regressors:  column %d is constant",i);

    for (j=0; j<dest->size1; j++) {
      u = gsl_matrix_get(dest,j,i);
      u = (u-mean);
      gsl_matrix_set(dest,j,i,u);
    }
  }
  return dest;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_integration.h>

#include "testint.c"

gsl_integration_workspace *w;
int inverse;

double
cauchy (double x, void *p)
{
  return gsl_ran_cauchy_pdf (x, *(double*)p);
}

double
gaussian (double x, void *p)
{
  return gsl_ran_gaussian_pdf (x, *(double*)p);
}

double
laplace (double x, void *p)
{
  return gsl_ran_laplace_pdf (x, *(double*)p);
}

double
rayleigh (double x, void *p)
{
  return gsl_ran_rayleigh_pdf (x, *(double*)p);
}

double
flat (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_flat_pdf (x, c[0], c[1]);
}

double
lognormal (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_lognormal_pdf (x, c[0], c[1]);
}

double
gamma (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_gamma_pdf (x, c[0], c[1]);
}

double
chisq (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_chisq_pdf (x, c[0]);
}

double
fdist (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_fdist_pdf (x, c[0], c[1]);
}

double
tdist (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_tdist_pdf (x, c[0]);
}

double
beta (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_beta_pdf (x, c[0], c[1]);
}

double
gumbel1 (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_gumbel1_pdf (x, c[0], c[1]);
}

double
gumbel2 (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_gumbel2_pdf (x, c[0], c[1]);
}

double
weibull (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_weibull_pdf (x, c[0], c[1]);
}

double
pareto (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_pareto_pdf (x, c[0], c[1]);
}

double
logistic (double x, void *p)
{
  double * c = (double *)p;
  return gsl_ran_logistic_pdf (x, c[0]);
}

double
integrate (gsl_function * f, double a, double b)
{
  double res = 0, err = 0;
  double tol = 1e-12;
  gsl_integration_qag (f, a, b, 0, tol, 1000, 0,  w, &res, &err);
  return res;
}

double
integrate_lower (gsl_function * f, double b)
{
  double res = 0, err = 0;
  qagil (f, b, 0, 1e-12, 1000, w, &res, &err);
  return res;
}

double
integrate_upper (gsl_function * f, double a)
{
  double res = 0, err = 0;
  qagiu (f, a, 0, 1e-12, 1000, w, &res, &err);
  return res;
}

int
test (const char * name, gsl_function * f, double x[], int N, char *fmt)
{
  int i;
  double res, err, sum = 0, sumerr = 0;

  printf ("void test_auto_%s (void);\n\n", name);

  printf ("void\ntest_auto_%s (void)\n{\n", name);

  /* gsl_set_error_handler_off(); */

  w = gsl_integration_workspace_alloc (1000);

  for (i = 0; i < N; i++)
    {
      res = 0;
      err = 0;

      if (x[0] < -1000)
        {
          if (x[i] < 0)
            {
              res = integrate_lower (f, x[i]);
            }
          else
            {
              res = 1 - integrate_upper (f, x[i]);
            }
          sum = res;
        }
      else
        {
          if (i == 0)
            sum += 0;
          else
            sum += integrate(f, x[i-1], x[i]);
        }

      if (res < 0)
        continue;

      printf (fmt, "_P", x[i], sum);

      if (inverse && (sum != 0 && sum != 1) && (x[i] == 0 || sum * 1e-4 < GSL_FN_EVAL(f,x[i]) * fabs(x[i])))
        printf (fmt, "_Pinv", sum, x[i]);
    }

  printf("\n");

  sum=0;
  sumerr=0;

  for (i = N-1; i >= 0; i--)
    {
      res = 0;
      err = 0;

      if (x[N-1] > 1000)
        {
          if (x[i] > 0)
            {
              res = integrate_upper (f, x[i]);
            }
          else
            {
              res = 1-integrate_lower (f, x[i]);
            }

          sum = res;
        }
      else
        {
          if (i == N-1)
            sum += 0;
          else
            sum += integrate(f, x[i], x[i+1]);
        }

      printf (fmt, "_Q", x[i], sum);

      if (inverse && (sum != 0 && sum != 1) && (x[i] == 0 || sum * 1e-4 < GSL_FN_EVAL(f,x[i]) * fabs(x[i])))
        printf (fmt, "_Qinv", sum, x[i]);
    }

  printf ("}\n\n");

  gsl_integration_workspace_free (w);

}

int
main (void)
{
  const int N = 43;

  double xall[] = {
    -1e10, -1e9, -1e8, -1e7, -1e6, -1e5, -1e4, -1e3, -1e2, -1e1,
    -1,
    -1e-1, -1e-2, -1e-3, -1e-4, -1e-5, -1e-6, -1e-7, -1e-8, -1e-9, -1e-10,
    0,
    1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
    1,
    1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10
  };

  const int Npos = 22;

  double xpos[] = {
    0,
    1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
    1,
    1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10
  };

  const int N01 = 23;

  double x01[] = {
    0,
    1e-10, 1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
    0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 0.9,
    0.99, 0.999, 0.9999, 0.99999,
    1.0
  };

#define TEST(name,params,range,n) { double p[] = { params } ; gsl_function f = {&name, p}; test(#name, &f, range, n, "    TEST(gsl_cdf_" #name "%s, (%.16e," #params "), %.12e, TEST_TOL6);\n"); }

#define TEST2(name,p1,p2,range,n) { double p[] = { p1,p2 } ; gsl_function f = {&name, p}; test(#name, &f, range, n, "    TEST(gsl_cdf_" #name "%s, (%.16e," #p1 "," #p2 "), %.12e, TEST_TOL6);\n"); }

#define TEST2A(desc,name,p1,p2,range,n) { double p[] = { p1,p2 } ; gsl_function f = {&name, p}; test(#desc, &f, range, n, "    TEST(gsl_cdf_" #name "%s, (%.16e," #p1 "," #p2 "), %.12e, TEST_TOL6);\n"); }

  inverse = 0;
  TEST2(beta, 1.3, 2.7, x01, N01);
  TEST2(fdist, 5.3, 2.7, xpos, Npos);

  inverse = 1;
  TEST(cauchy, 1.3, xall, N);
  TEST(gaussian, 1.3, xall, N);
  TEST(laplace, 1.3, xall, N);
  TEST(rayleigh, 1.3, xpos, Npos);

  TEST2(flat, 1.3, 750.0, xpos, Npos);
  TEST2(lognormal, 1.3,2.7, xpos, Npos);
  TEST2(gamma, 1.3,2.7, xpos, Npos);

  TEST(chisq, 1.3, xpos, Npos);
  TEST(tdist, 1.3, xall, N);

  TEST2(gumbel1, 1.3, 2.7, xall, N);
  TEST2(gumbel2, 1.3, 2.7, xpos, Npos);
  TEST2(weibull, 1.3, 2.7, xpos, Npos);

  TEST2(pareto, 1.3, 2.7, xpos, Npos);
  TEST(logistic, 1.3, xall, N);

  TEST2A(gammalarge, gamma, 1.3,123.0, xpos, Npos);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_permute_double.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

unsigned int p5[120][5] = {
  {0, 1, 2, 3, 4}, {0, 1, 2, 4, 3}, {0, 1, 3, 2, 4}, {0, 1, 3, 4, 2},
  {0, 1, 4, 2, 3}, {0, 1, 4, 3, 2}, {0, 2, 1, 3, 4}, {0, 2, 1, 4, 3},
  {0, 2, 3, 1, 4}, {0, 2, 3, 4, 1}, {0, 2, 4, 1, 3}, {0, 2, 4, 3, 1},
  {0, 3, 1, 2, 4}, {0, 3, 1, 4, 2}, {0, 3, 2, 1, 4}, {0, 3, 2, 4, 1},
  {0, 3, 4, 1, 2}, {0, 3, 4, 2, 1}, {0, 4, 1, 2, 3}, {0, 4, 1, 3, 2},
  {0, 4, 2, 1, 3}, {0, 4, 2, 3, 1}, {0, 4, 3, 1, 2}, {0, 4, 3, 2, 1},
  {1, 0, 2, 3, 4}, {1, 0, 2, 4, 3}, {1, 0, 3, 2, 4}, {1, 0, 3, 4, 2},
  {1, 0, 4, 2, 3}, {1, 0, 4, 3, 2}, {1, 2, 0, 3, 4}, {1, 2, 0, 4, 3},
  {1, 2, 3, 0, 4}, {1, 2, 3, 4, 0}, {1, 2, 4, 0, 3}, {1, 2, 4, 3, 0},
  {1, 3, 0, 2, 4}, {1, 3, 0, 4, 2}, {1, 3, 2, 0, 4}, {1, 3, 2, 4, 0},
  {1, 3, 4, 0, 2}, {1, 3, 4, 2, 0}, {1, 4, 0, 2, 3}, {1, 4, 0, 3, 2},
  {1, 4, 2, 0, 3}, {1, 4, 2, 3, 0}, {1, 4, 3, 0, 2}, {1, 4, 3, 2, 0},
  {2, 0, 1, 3, 4}, {2, 0, 1, 4, 3}, {2, 0, 3, 1, 4}, {2, 0, 3, 4, 1},
  {2, 0, 4, 1, 3}, {2, 0, 4, 3, 1}, {2, 1, 0, 3, 4}, {2, 1, 0, 4, 3},
  {2, 1, 3, 0, 4}, {2, 1, 3, 4, 0}, {2, 1, 4, 0, 3}, {2, 1, 4, 3, 0},
  {2, 3, 0, 1, 4}, {2, 3, 0, 4, 1}, {2, 3, 1, 0, 4}, {2, 3, 1, 4, 0},
  {2, 3, 4, 0, 1}, {2, 3, 4, 1, 0}, {2, 4, 0, 1, 3}, {2, 4, 0, 3, 1},
  {2, 4, 1, 0, 3}, {2, 4, 1, 3, 0}, {2, 4, 3, 0, 1}, {2, 4, 3, 1, 0},
  {3, 0, 1, 2, 4}, {3, 0, 1, 4, 2}, {3, 0, 2, 1, 4}, {3, 0, 2, 4, 1},
  {3, 0, 4, 1, 2}, {3, 0, 4, 2, 1}, {3, 1, 0, 2, 4}, {3, 1, 0, 4, 2},
  {3, 1, 2, 0, 4}, {3, 1, 2, 4, 0}, {3, 1, 4, 0, 2}, {3, 1, 4, 2, 0},
  {3, 2, 0, 1, 4}, {3, 2, 0, 4, 1}, {3, 2, 1, 0, 4}, {3, 2, 1, 4, 0},
  {3, 2, 4, 0, 1}, {3, 2, 4, 1, 0}, {3, 4, 0, 1, 2}, {3, 4, 0, 2, 1},
  {3, 4, 1, 0, 2}, {3, 4, 1, 2, 0}, {3, 4, 2, 0, 1}, {3, 4, 2, 1, 0},
  {4, 0, 1, 2, 3}, {4, 0, 1, 3, 2}, {4, 0, 2, 1, 3}, {4, 0, 2, 3, 1},
  {4, 0, 3, 1, 2}, {4, 0, 3, 2, 1}, {4, 1, 0, 2, 3}, {4, 1, 0, 3, 2},
  {4, 1, 2, 0, 3}, {4, 1, 2, 3, 0}, {4, 1, 3, 0, 2}, {4, 1, 3, 2, 0},
  {4, 2, 0, 1, 3}, {4, 2, 0, 3, 1}, {4, 2, 1, 0, 3}, {4, 2, 1, 3, 0},
  {4, 2, 3, 0, 1}, {4, 2, 3, 1, 0}, {4, 3, 0, 1, 2}, {4, 3, 0, 2, 1},
  {4, 3, 1, 0, 2}, {4, 3, 1, 2, 0}, {4, 3, 2, 0, 1}, {4, 3, 2, 1, 0}
} ;

unsigned int c5[120][5] = {
  {4, 3, 2, 1, 0}, {3, 4, 2, 1, 0}, {4, 2, 3, 1, 0}, {2, 3, 4, 1, 0},
  {2, 4, 3, 1, 0}, {3, 2, 4, 1, 0}, {4, 3, 1, 2, 0}, {3, 4, 1, 2, 0},
  {4, 1, 2, 3, 0}, {1, 2, 3, 4, 0}, {1, 2, 4, 3, 0}, {3, 1, 2, 4, 0},
  {4, 1, 3, 2, 0}, {1, 3, 4, 2, 0}, {4, 2, 1, 3, 0}, {2, 1, 3, 4, 0},
  {2, 4, 1, 3, 0}, {1, 3, 2, 4, 0}, {1, 4, 3, 2, 0}, {3, 1, 4, 2, 0},
  {2, 1, 4, 3, 0}, {3, 2, 1, 4, 0}, {1, 4, 2, 3, 0}, {2, 3, 1, 4, 0},
  {4, 3, 2, 0, 1}, {3, 4, 2, 0, 1}, {4, 2, 3, 0, 1}, {2, 3, 4, 0, 1},
  {2, 4, 3, 0, 1}, {3, 2, 4, 0, 1}, {4, 3, 0, 1, 2}, {3, 4, 0, 1, 2},
  {4, 0, 1, 2, 3}, {0, 1, 2, 3, 4}, {0, 1, 2, 4, 3}, {3, 0, 1, 2, 4},
  {4, 0, 1, 3, 2}, {0, 1, 3, 4, 2}, {4, 2, 0, 1, 3}, {2, 0, 1, 3, 4},
  {2, 4, 0, 1, 3}, {0, 1, 3, 2, 4}, {0, 1, 4, 3, 2}, {3, 0, 1, 4, 2},
  {2, 0, 1, 4, 3}, {3, 2, 0, 1, 4}, {0, 1, 4, 2, 3}, {2, 3, 0, 1, 4},
  {4, 3, 0, 2, 1}, {3, 4, 0, 2, 1}, {4, 0, 2, 3, 1}, {0, 2, 3, 4, 1},
  {0, 2, 4, 3, 1}, {3, 0, 2, 4, 1}, {4, 3, 1, 0, 2}, {3, 4, 1, 0, 2},
  {4, 1, 0, 2, 3}, {1, 0, 2, 3, 4}, {1, 0, 2, 4, 3}, {3, 1, 0, 2, 4},
  {4, 1, 3, 0, 2}, {1, 3, 4, 0, 2}, {4, 0, 2, 1, 3}, {0, 2, 1, 3, 4},
  {0, 2, 4, 1, 3}, {1, 3, 0, 2, 4}, {1, 4, 3, 0, 2}, {3, 1, 4, 0, 2},
  {0, 2, 1, 4, 3}, {3, 0, 2, 1, 4}, {1, 4, 0, 2, 3}, {0, 2, 3, 1, 4},
  {4, 0, 3, 2, 1}, {0, 3, 4, 2, 1}, {4, 2, 0, 3, 1}, {2, 0, 3, 4, 1},
  {2, 4, 0, 3, 1}, {0, 3, 2, 4, 1}, {4, 1, 0, 3, 2}, {1, 0, 3, 4, 2},
  {4, 2, 1, 0, 3}, {2, 1, 0, 3, 4}, {2, 4, 1, 0, 3}, {1, 0, 3, 2, 4},
  {4, 0, 3, 1, 2}, {0, 3, 4, 1, 2}, {4, 1, 2, 0, 3}, {1, 2, 0, 3, 4},
  {1, 2, 4, 0, 3}, {0, 3, 1, 2, 4}, {0, 3, 1, 4, 2}, {1, 4, 0, 3, 2},
  {1, 4, 2, 0, 3}, {0, 3, 2, 1, 4}, {2, 1, 4, 0, 3}, {2, 0, 3, 1, 4},
  {0, 4, 3, 2, 1}, {3, 0, 4, 2, 1}, {2, 0, 4, 3, 1}, {3, 2, 0, 4, 1},
  {0, 4, 2, 3, 1}, {2, 3, 0, 4, 1}, {1, 0, 4, 3, 2}, {3, 1, 0, 4, 2},
  {2, 1, 0, 4, 3}, {3, 2, 1, 0, 4}, {1, 0, 4, 2, 3}, {2, 3, 1, 0, 4},
  {0, 4, 3, 1, 2}, {3, 0, 4, 1, 2}, {1, 2, 0, 4, 3}, {3, 1, 2, 0, 4},
  {0, 4, 1, 2, 3}, {1, 2, 3, 0, 4}, {1, 3, 0, 4, 2}, {0, 4, 1, 3, 2},
  {0, 4, 2, 1, 3}, {1, 3, 2, 0, 4}, {2, 0, 4, 1, 3}, {2, 1, 3, 0, 4}
} ;

unsigned int cycles[120] = {
  5, 4, 4, 3, 3, 4, 4, 3, 3, 2,
  2, 3, 3, 2, 4, 3, 3, 2, 2, 3,
  3, 4, 2, 3, 4, 3, 3, 2, 2, 3,
  3, 2, 2, 1, 1, 2, 2, 1, 3, 2,
  2, 1, 1, 2, 2, 3, 1, 2, 3, 2,
  2, 1, 1, 2, 4, 3, 3, 2, 2, 3,
  3, 2, 2, 1, 1, 2, 2, 3, 1, 2,
  2, 1, 2, 1, 3, 2, 2, 1, 3, 2,
  4, 3, 3, 2, 2, 1, 3, 2, 2, 1,
  1, 2, 2, 1, 3, 2, 1, 2, 2, 3,
  1, 2, 2, 3, 3, 4, 2, 3, 1, 2,
  2, 3, 1, 2, 2, 1, 1, 2, 2, 3
} ;

unsigned int inversions[120] = {
  0, 1, 1, 2, 2, 3, 1, 2, 2, 3,
  3, 4, 2, 3, 3, 4, 4, 5, 3, 4,
  4, 5, 5, 6, 1, 2, 2, 3, 3, 4,
  2, 3, 3, 4, 4, 5, 3, 4, 4, 5,
  5, 6, 4, 5, 5, 6, 6, 7, 2, 3,
  3, 4, 4, 5, 3, 4, 4, 5, 5, 6,
  4, 5, 5, 6, 6, 7, 5, 6, 6, 7,
  7, 8, 3, 4, 4, 5, 5, 6, 4, 5,
  5, 6, 6, 7, 5, 6, 6, 7, 7, 8,
  6, 7, 7, 8, 8, 9, 4, 5, 5, 6,
  6, 7, 5, 6, 6, 7, 7, 8, 6, 7,
  7, 8, 8, 9, 7, 8, 8, 9, 9, 10
} ;

int
main (void)
{
  gsl_ieee_env_setup ();

  {
    int i = 0, j, status = 0;

    gsl_permutation * p ;

    p = gsl_permutation_alloc (5);

    gsl_permutation_init (p);

    do
      {
        for (j = 0; j < 5; j++)
          {
            status |= (p->data[j] != p5[i][j]);
          }

        i++;
      }
    while (gsl_permutation_next(p) == GSL_SUCCESS);

    gsl_test(status, "gsl_permutation_next, 5-th order permutation, 120 steps");

    do
      {
        i--;

        for (j = 0; j < 5; j++)
          {
            status |= (p->data[j] != p5[i][j]);
          }
      }
    while (gsl_permutation_prev(p) == GSL_SUCCESS);

    gsl_test(status, "gsl_permutation_prev, 5-th order permutation, 120 steps");

    gsl_permutation_free (p);
  }

#ifdef JUNK
  {
    int i;
    int status = 0 ;

    gsl_permutation * p1 = gsl_permutation_alloc (5);
    gsl_permutation * p2 = gsl_permutation_alloc (5);
    gsl_permutation * p = gsl_permutation_alloc (5);

    double v[5] = { 100.0, 101.0, 102.0, 103.0, 104.0 } ;

    gsl_permutation_init (p1);

    do
      {
        gsl_permutation_init (p2);

        do
          {
            double x[5], y[5];

            /* Compute x= p1 p2 v */
            memcpy (x, v, 5 * sizeof(double));
            gsl_permute (p2->data, x, 1, 5);
            gsl_permute (p1->data, x, 1, 5);

            /* Compute P= p1 p2, y = P v */
            gsl_permutation_mul (p, p1, p2);
            memcpy (y, v, 5 * sizeof(double));
            gsl_permute (p->data, y, 1, 5);

            for (i = 0; i < 5; i++)
              {
                if (x[i] != y[i])
                  status = 1;
              }

            if (status == 1)
              break;
          }
        while (gsl_permutation_next(p2) == GSL_SUCCESS);

        if (status == 1)
          break;
      }
    while (gsl_permutation_next(p1) == GSL_SUCCESS);

    gsl_permutation_free (p1);
    gsl_permutation_free (p2);
    gsl_permutation_free (p);

    gsl_test(status, "gsl_permutation_mul, all 5-th order combinations");
  }
#endif

  /* testing cycles representations */
  {
    int i = 0, j, status = 0;

    gsl_permutation * p = gsl_permutation_alloc (5);

    gsl_permutation * plin = gsl_permutation_alloc (5);
    gsl_permutation * pcan = gsl_permutation_alloc (5);

    gsl_permutation_init (p);

    do
      {
        gsl_permutation_memcpy (plin, p);

        for (j = 0; j < 5; j++)
          {
            pcan->data[j] = 0;
          }

        gsl_permutation_linear_to_canonical (pcan, plin);

        for (j = 0; j < 5; j++)
          {
            status |= (pcan->data[j] != c5[i][j]);
          }

        status |= (gsl_permutation_canonical_cycles (pcan) != cycles[i]);

        status |= (gsl_permutation_linear_cycles (plin) != cycles[i]);

        for (j = 0; j < 5; j++)
          {
            plin->data[j] = 0;
          }

        gsl_permutation_canonical_to_linear (plin, pcan);

        for (j = 0; j < 5; j++)
          {
            status |= (plin->data[j] != p5[i][j]);
          }

        i++;
      }
    while (gsl_permutation_next(p) == GSL_SUCCESS);

    gsl_permutation_free (p);
    gsl_permutation_free (plin);
    gsl_permutation_free (pcan);

    gsl_test (status, "gsl_permutation canonical conversion, 5-th order permutation, 120 steps");
  }

  /* testing number of inversions */
  {
    int i = 0, status = 0;

    gsl_permutation * p = gsl_permutation_alloc (5);

    gsl_permutation_init (p);

    do
      {
        status |= gsl_permutation_inversions (p) != inversions[i];
        i++;
      }
    while (gsl_permutation_next(p) == GSL_SUCCESS);

    gsl_permutation_free (p);

    gsl_test (status, "gsl_permutation_inversions, 5-th order permutation, 120 steps");
  }

  exit (gsl_test_summary());
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/**
 *
 * @file al4san/codelet_dlacpy.c
 *
 * @copyright 2009-2014 The University of Tennessee and The University of
 *                      Tennessee Research Foundation. All rights reserved.
 * @copyright 2012-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      Univ. Bordeaux. All rights reserved.
 * @copyright 2017-2018 King Abdullah University of Science and Technology (KAUST).
 *                      All rights reserved.
 * @brief Chameleon dlacpy AL4SAN codelet
 *
 * @version 1.1.0
 * @author Rabab Alomairy
 * @date 2018-09-19
 * @generated d Wed Jan 16 11:13:56 2019
 *
 */
/**
 *
 * @ingroup CORE_CHAMELEON_Complex64_t
 *
 */
#include <lapacke.h>
#include <al4san.h>
#include <runtime/al4san_quark.h>
#include <runtime/al4san_starpu.h>

AL4SAN_TASK_CPU(lacpy1, lacpy_cpu_func1)

void EIG_AL4SAN_CORE_dlacpy1(AL4SAN_option_t *options,
                         int uplo, int m, int n, int nb,
                         const AL4SAN_desc_t *A, int Am, int An, int lda,
                         const AL4SAN_desc_t *B, int Bm, int Bn, int ldb )
{
    (void)nb;

    AL4SAN_Insert_Task(AL4SAN_TASK(lacpy1), (AL4SAN_option_t*)options,
        AL4SAN_VALUE,                       &uplo,                                        sizeof(int),
        AL4SAN_VALUE,                       &m,                                           sizeof(int),
        AL4SAN_VALUE,                       &n,                                           sizeof(int),
        AL4SAN_INPUT,                       AL4SAN_ADDR(A, double, Am, An),  AL4SAN_DEP,
        AL4SAN_VALUE,                       &lda,                                         sizeof(int),
        AL4SAN_OUTPUT,                      AL4SAN_ADDR(B, double, Bm, Bn),  AL4SAN_DEP,
        AL4SAN_VALUE,                       &ldb,                                         sizeof(int),
        AL4SAN_PRIORITY,                    options->priority,                     sizeof(int),
        AL4SAN_LABEL,                       "dlacpy",                                     sizeof(char),
        AL4SAN_COLOR,                       "white",                                      sizeof(char),
        ARG_END);
}

#if !defined(CHAMELEON_SIMULATION)
 void lacpy_cpu_func1(AL4SAN_arg_list *al4san_arg)
{
    int uplo;
    int M;
    int N;
    int displA;
    int displB;
    const double *A;
    int LDA;
    double *B;
    int LDB;

    AL4SAN_Unpack_Arg(al4san_arg, &uplo, &M, &N, &A, &LDA, &B, &LDB);

    CORE_dlacpy(uplo, M, N, A, LDA, B , LDB);
}
#endif /* !defined(CHAMELEON_SIMULATION) */


text:
#pragma once

#include "Common.h"
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/depth_first_search.hpp>
#include <boost/graph/edge_list.hpp>
#include <boost/graph/graph_selectors.hpp>
#include <boost/graph/properties.hpp>
#include <boost/graph/visitors.hpp>
#include <any>
#include <forward_list>
#include <functional>
#include <gsl/span>
#include <set>
#include <string_view>
#include <variant>

namespace bcos::scheduler
{
class GraphKeyLocks
{
public:
    using Ptr = std::shared_ptr<GraphKeyLocks>;
    using KeyLock = std::tuple<std::string, std::string>;
    using KeyLockView = std::tuple<std::string_view, std::string_view>;

    using ContractView = std::string_view;
    using KeyView = std::string_view;

    GraphKeyLocks() = default;
    GraphKeyLocks(const GraphKeyLocks&) = delete;
    GraphKeyLocks(GraphKeyLocks&&) = delete;
    GraphKeyLocks& operator=(const GraphKeyLocks&) = delete;
    GraphKeyLocks& operator=(GraphKeyLocks&&) = delete;

    bool batchAcquireKeyLock(std::string_view contract, gsl::span<std::string const> keyLocks,
        ContextID contextID, Seq seq);

    bool acquireKeyLock(
        std::string_view contract, std::string_view key, ContextID contextID, Seq seq);

    std::vector<std::string> getKeyLocksNotHoldingByContext(
        std::string_view contract, ContextID excludeContextID) const;

    void releaseKeyLocks(ContextID contextID, Seq seq);

    bool detectDeadLock(ContextID contextID);

    struct Vertex : public std::variant<ContextID, KeyLock>
    {
        using std::variant<ContextID, KeyLock>::variant;

        bool operator==(const KeyLockView& rhs) const
        {
            if (index() != 1)
            {
                return false;
            }

            auto view = std::make_tuple(std::string_view(std::get<0>(std::get<1>(*this))),
                std::string_view(std::get<1>(std::get<1>(*this))));
            return view == rhs;
        }
    };

private:
    struct VertexIterator
    {
        using kind = boost::vertex_property_tag;
    };
    using VertexProperty = boost::property<VertexIterator, const Vertex*>;
    struct EdgeSeq
    {
        using kind = boost::edge_property_tag;
    };
    using EdgeProperty = boost::property<EdgeSeq, int64_t>;

    using Graph = boost::adjacency_list<boost::multisetS, boost::multisetS, boost::bidirectionalS,
        VertexProperty, EdgeProperty>;
    using VertexPropertyTag = boost::property_map<Graph, VertexIterator>::const_type;
    using EdgePropertyTag = boost::property_map<Graph, EdgeSeq>::const_type;
    using VertexID = Graph::vertex_descriptor;
    using EdgeID = Graph::edge_descriptor;

    Graph m_graph;
    std::map<Vertex, VertexID, std::less<>> m_vertexes;

    VertexID touchContext(ContextID contextID);
    VertexID touchKeyLock(KeyLockView keylockView);

    void addEdge(VertexID source, VertexID target, Seq seq);
    void removeEdge(VertexID source, VertexID target, Seq seq);
};

}  // namespace bcos::scheduler

namespace std
{
inline bool operator<(const bcos::scheduler::GraphKeyLocks::Vertex& lhs,
    const bcos::scheduler::GraphKeyLocks::KeyLockView& rhs)
{
    if (lhs.index() != 1)
    {
        return true;
    }

    auto view = std::make_tuple(std::string_view(std::get<0>(std::get<1>(lhs))),
        std::string_view(std::get<1>(std::get<1>(lhs))));
    return view < rhs;
}

inline bool operator<(const bcos::scheduler::GraphKeyLocks::KeyLockView& lhs,
    const bcos::scheduler::GraphKeyLocks::Vertex& rhs)
{
    if (rhs.index() != 1)
    {
        return false;
    }

    auto view = std::make_tuple(std::string_view(std::get<0>(std::get<1>(rhs))),
        std::string_view(std::get<1>(std::get<1>(rhs))));
    return lhs < view;
}
}  // namespace std


text:
/* filter/rmedian.c
 *
 * Copyright (C) 2018 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.
 */

/*
 * This module contains routines related to the recursive median filter. The
 * algorithm is based on this paper,
 *
 * [1] S-J Ko, Y. H. Lee, and A. T. Fam, Efficient Implementation of One-Dimensional
 * Recursive Median Filters, IEEE Transactions on Circuits and Systems, Vol 37,
 * No 11, 1990.
 */

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_movstat.h>

typedef struct
{
  const gsl_movstat_accum * minmax_acc; /* minimum/maximum accumulator */
  void *minmax_state;                   /* minimum/maximum accumulator workspace */
} rmedian_state_t;

static size_t rmedian_size(const size_t n);
static int rmedian_init(const size_t n, void * vstate);
static int rmedian_insert(const double x, void * vstate);
static int rmedian_delete(void * vstate);
static int rmedian_get(void * params, double * result, const void * vstate);

static const gsl_movstat_accum rmedian_accum_type;

gsl_filter_rmedian_workspace *
gsl_filter_rmedian_alloc(const size_t K)
{
  gsl_filter_rmedian_workspace *w;
  size_t state_size;

  w = calloc(1, sizeof(gsl_filter_rmedian_workspace));
  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->H = K / 2;
  w->K = 2*w->H + 1;
  w->minmaxacc = gsl_movstat_accum_minmax;

  w->window = malloc(w->K * sizeof(double));
  if (w->window == NULL)
    {
      gsl_filter_rmedian_free(w);
      GSL_ERROR_NULL ("failed to allocate space for window", GSL_ENOMEM);
    }

  state_size = rmedian_size(w->H + 1);

  w->state = malloc(state_size);
  if (w->state == NULL)
    {
      gsl_filter_rmedian_free(w);
      GSL_ERROR_NULL ("failed to allocate space for min/max state", GSL_ENOMEM);
    }

  w->movstat_workspace_p = gsl_movstat_alloc_with_size(state_size, 0, w->H);
  if (!w->movstat_workspace_p)
    {
      gsl_filter_rmedian_free(w);
      GSL_ERROR_NULL ("failed to allocate space for movstat workspace", GSL_ENOMEM);
    }

  return w;
}

void
gsl_filter_rmedian_free(gsl_filter_rmedian_workspace * w)
{
  if (w->state)
    free(w->state);

  if (w->window)
    free(w->window);

  if (w->movstat_workspace_p)
    gsl_movstat_free(w->movstat_workspace_p);

  free(w);
}

/*
gsl_filter_rmedian()
  Recursive median filter

Inputs: endtype - end point handling
        x       - input vector
        y       - output vector
        w       - workspace
*/

int
gsl_filter_rmedian(const gsl_filter_end_t endtype, const gsl_vector * x, gsl_vector * y, gsl_filter_rmedian_workspace * w)
{
  if (x->size != y->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else
    {
      int status = GSL_SUCCESS;
      const size_t n = x->size;
      const int H = (int) w->H;
      double yprev;
      int wsize;

      /* find median of first window to initialize filter */
      wsize = gsl_movstat_fill(endtype, x, 0, H, H, w->window);
      yprev = gsl_stats_median(w->window, 1, wsize);
      gsl_vector_set(y, 0, yprev);

      if (x->size > 1)
        {
          gsl_vector_const_view xv = gsl_vector_const_subvector(x, 1, n - 1);
          gsl_vector_view yv = gsl_vector_subvector(y, 1, n - 1);

          /* apply recursive median filter to x[2:end] */
          status = gsl_movstat_apply_accum(endtype, &xv.vector, &rmedian_accum_type, (void *) &yprev, &yv.vector,
                                           NULL, w->movstat_workspace_p);
        }

      return status;
    }
}

static size_t
rmedian_size(const size_t n)
{
  size_t size = 0;
  const gsl_movstat_accum * acc = gsl_movstat_accum_minmax;

  size += sizeof(rmedian_state_t);
  size += (acc->size)(n);

  return size;
}

static int
rmedian_init(const size_t n, void * vstate)
{
  rmedian_state_t * state = (rmedian_state_t *) vstate;

  state->minmax_acc = gsl_movstat_accum_minmax;
  state->minmax_state = (void *) ((unsigned char *) vstate + sizeof(rmedian_state_t));

  (state->minmax_acc->init)(n, state->minmax_state);

  return GSL_SUCCESS;
}

static int
rmedian_insert(const double x, void * vstate)
{
  rmedian_state_t * state = (rmedian_state_t *) vstate;
  return (state->minmax_acc->insert)(x, state->minmax_state);
}

static int
rmedian_delete(void * vstate)
{
  rmedian_state_t * state = (rmedian_state_t *) vstate;
  return (state->minmax_acc->delete_oldest)(state->minmax_state);
}

static int
rmedian_get(void * params, double * result, const void * vstate)
{
  const rmedian_state_t * state = (const rmedian_state_t *) vstate;
  double *yprev = (double *) params; /* previous filter output */
  double y;                          /* new filter output */
  double xminmax[2];

  /* get minimum/maximum values of {x_i,...,x_{i+H}} */
  (state->minmax_acc->get)(NULL, xminmax, state->minmax_state);

  /* y = median [ yprev, xmin, xmax ] */
  if (*yprev <= xminmax[0])
    y = xminmax[0];
  else if (*yprev <= xminmax[1])
    y = *yprev;
  else
    y = xminmax[1];

  *result = y;
  *yprev = y;

  return GSL_SUCCESS;
}

static const gsl_movstat_accum rmedian_accum_type =
{
  rmedian_size,
  rmedian_init,
  rmedian_insert,
  rmedian_delete,
  rmedian_get
};

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <gsl/gsl_rng.h>

#include "allvars.h"

// file names
char fname_param[100];
char fname_con[100], fname_line[100], fname_results[100];

FILE *fp_results;

// transfer function parameters
double tau_lim_low, tau_lim_up;
int nc, nc_lim_low, nc_lim_up, nc_best;
double *grid_tau;
const int ntau=200;
double *TF_tau, *TF;

// mcmc
const int ntheta_max = 200;
int ntheta;
int nmcmc, nbuilt;
double *workspace;
int flag_detrend, flag_sim, flag_mcmc;
double *Cmat, *ICmat, *Smat, *Nmat, *USmat, *ASmat;
double *ICvmat, *Tmat1, *Tmat2, *Tmat3, *Tmat4;

double *cov_matrix;
double **theta_range, *theta_best, *theta_best_var;
double *theta_input, *sigma_input;
int *theta_fixed;

double *theta_best_con, *theta_best_var_con;

const gsl_rng_type * gsl_T;
const gsl_rng * gsl_r;

// light curve data
const int ndata_max=500;
int ncon_data, nline_data, nall_data;
double *Tcon_data, *Tline_data, *Fcon_data, *Fcerrs_data, *Fline_data, *Flerrs_data;
double *Fall_data;
double scale_con, scale_line;
double len_con, len_line, cad_con, cad_line;

int ncon=500, nline=500, nall;
double *Tcon, *Fcon, *Fcerrs, *Tline, *Fline, *Flerrs;

// error exit
char str_error_exit[200];

// mathematic functions
int *workspace_ipiv;

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
//STARTWHOLE
static char help[] = "Solve a tridiagonal system of arbitrary size.\n"
"Option prefix = tri_.\n";

#include <petsc.h>

int main(int argc,char **args) {
    PetscErrorCode ierr;
    Vec    x, b, xexact;
    Mat    A;
    KSP    ksp;
    int    m = 4, i, Istart, Iend, j[3];
    double v[3], xval, errnorm;

    PetscInitialize(&argc,&args,NULL,help);

    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"tri_","options for tri",""); CHKERRQ(ierr);
    ierr = PetscOptionsInt("-m","dimension of linear system","tri.c",m,&m,NULL); CHKERRQ(ierr);
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);

    ierr = VecCreate(PETSC_COMM_WORLD,&x); CHKERRQ(ierr);
    ierr = VecSetSizes(x,PETSC_DECIDE,m); CHKERRQ(ierr);
    ierr = VecSetFromOptions(x); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&b); CHKERRQ(ierr);
    ierr = VecDuplicate(x,&xexact); CHKERRQ(ierr);

    ierr = MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
    ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,m); CHKERRQ(ierr);
    ierr = MatSetOptionsPrefix(A,"a_"); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A); CHKERRQ(ierr);
    ierr = MatSetUp(A); CHKERRQ(ierr);
    ierr = MatGetOwnershipRange(A,&Istart,&Iend); CHKERRQ(ierr);
    for (i=Istart; i<Iend; i++) {
        if (i == 0) {
            v[0] = 3.0;  v[1] = -1.0;
            j[0] = 0;    j[1] = 1;
            ierr = MatSetValues(A,1,&i,2,j,v,INSERT_VALUES); CHKERRQ(ierr);
        } else {
            v[0] = -1.0;  v[1] = 3.0;  v[2] = -1.0;
            j[0] = i-1;   j[1] = i;    j[2] = i+1;
            if (i == m-1) {
                ierr = MatSetValues(A,1,&i,2,j,v,INSERT_VALUES); CHKERRQ(ierr);
            } else {
                ierr = MatSetValues(A,1,&i,3,j,v,INSERT_VALUES); CHKERRQ(ierr);
            }
        }
        xval = exp(cos(i));
        ierr = VecSetValues(xexact,1,&i,&xval,INSERT_VALUES); CHKERRQ(ierr);
    }
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(xexact); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(xexact); CHKERRQ(ierr);
    ierr = MatMult(A,xexact,b); CHKERRQ(ierr);

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);

    ierr = VecAXPY(x,-1.0,xexact); CHKERRQ(ierr);
    ierr = VecNorm(x,NORM_2,&errnorm); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
    "error for m = %d system is |x-xexact|_2 = %.1e\n",m,errnorm); CHKERRQ(ierr);

    KSPDestroy(&ksp);  MatDestroy(&A);
    VecDestroy(&x);  VecDestroy(&b);  VecDestroy(&xexact);
    return PetscFinalize();
}
//ENDWHOLE


text:
/*
 * Copyright 2009-2019 The VOTCA-MPIP Development Team (http://www.votca.org)
 *
 * 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.
 *
 * author: Denis Andrienko
 */

#ifndef __VOTCA_CTP_GSL_H
#define __VOTCA_CTP_GSL_H

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_eigen.h>

#include <boost/numeric/conversion/cast.hpp>

#define DEBUG_LINALG
#ifdef DEBUG_LINALG
#include <chrono>
#endif

using namespace std;
namespace ub = boost::numeric::ublas;

namespace votca { namespace ctp { namespace linalg {

    /**
     * \brief inverts A
     * @param A symmetric positive definite matrix
     * @param V inverse matrix
     */
inline void invert_symm( const ub::matrix<double> &A, ub::matrix<double> &V){
        // matrix inversion using gsl

    #ifdef DEBUG_LINALG
        auto start = std::chrono::system_clock::now();
    #endif

        gsl_error_handler_t *handler = gsl_set_error_handler_off();
	const size_t N = A.size1();
	// signum s (for LU decomposition)
	int s;
        //make copy of A as A is destroyed by GSL
        ub::matrix<double> work=A;
        V.resize(N, N, false);

	// Define all the used matrices
        gsl_matrix_view A_view = gsl_matrix_view_array(&work(0,0), N, N);
        gsl_matrix_view V_view = gsl_matrix_view_array(&V(0,0), N, N);
	gsl_permutation * perm = gsl_permutation_alloc (N);

	// Make LU decomposition of matrix A_view
	gsl_linalg_LU_decomp (&A_view.matrix, perm, &s);

	// Invert the matrix A_view
	(void)gsl_linalg_LU_invert (&A_view.matrix, perm, &V_view.matrix);

        gsl_set_error_handler(handler);

    #ifdef DEBUG_LINALG
        auto end = std::chrono::system_clock::now();

        std::chrono::duration<double> elapsed_seconds = end - start;

        std::cout << "\n... ... ... \x1b[0;34mGSL Invert"
                  << "[" << A.size1() << "x" << A.size2() << "]"
                  << " time: " << elapsed_seconds.count()
                  << "s\x1b[0;39m";
    #endif

}

/**
 * \brief eigenvalues of a symmetric matrix A*x=E*x
 * @param E vector of eigenvalues
 * @param V input: matrix to diagonalize
 * @param V output: eigenvectors
 *
 * This function wraps gsl_eigen_symmv / DSYEV
 *
 */
inline void eigenvalues_symm(const ub::matrix<double> &A,
                                     ub::vector<double> &E,
                                     ub::matrix<double> &V)
{

   #ifdef DEBUG_LINALG
        auto start = std::chrono::system_clock::now();
    #endif

	gsl_error_handler_t *handler = gsl_set_error_handler_off();
	const size_t N = A.size1();

        // gsl does not handle conversion of a symmetric_matrix
        ub::matrix<double> _A( N,N );
        _A = A;

	E.resize(N, false);
	V.resize(N, N, false);
	gsl_matrix_view A_view = gsl_matrix_view_array(&_A(0,0), N, N);
	gsl_vector_view E_view = gsl_vector_view_array(&E(0), N);
	gsl_matrix_view V_view = gsl_matrix_view_array(&V(0,0), N, N);
	gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(N);

#ifdef NDEBUG
	(void)gsl_eigen_symmv(&A_view.matrix, &E_view.vector, &V_view.matrix, w);
#else
	int status = gsl_eigen_symmv(&A_view.matrix, &E_view.vector, &V_view.matrix, w);
        assert(status == 0);
#endif
	gsl_eigen_symmv_sort(&E_view.vector, &V_view.matrix, GSL_EIGEN_SORT_VAL_ASC);
	gsl_eigen_symmv_free(w);
	gsl_set_error_handler(handler);
	//return (status == 0);

    #ifdef DEBUG_LINALG
        auto end = std::chrono::system_clock::now();
        std::chrono::duration<double> elapsed_seconds = end - start;
        std::cout << "\n... ... ... \x1b[0;34mGSL Eigenvalues"
                  << "[" << A.size1() << "x" << A.size2() << "]"
                  << " time: " << elapsed_seconds.count()
                  << "s\x1b[0;39m";
    #endif

};

}}}

/*
 * Since ublas has very slow prod() function, we overwrite it here
 * by the gsl (blas) product function. GSL has separate implementations
 * for floats and doubles hence we have two specializations for these types
 * Separate templates are provided for diagonal and symmetric matreces
 * and their combinations. To check  whether a specific routine is called,
 * define the GSLDEBUG flag
 */
namespace boost { namespace numeric { namespace ublas {

// specialization for double precision (dgemm)
template<class F, class A>
inline matrix<double,F,A>    // prod( m1, m2 )
prod(const matrix<double,F,A> &m1, const matrix<double,F,A> &m2)
{

    #ifdef DEBUG_LINALG
        auto start = std::chrono::system_clock::now();
    #endif

       gsl_matrix_const_view mA = gsl_matrix_const_view_array (&m1(0,0), m1.size1(), m1.size2());
       gsl_matrix_const_view mB = gsl_matrix_const_view_array (&m2(0,0), m2.size1(), m2.size2());

       boost::numeric::ublas::matrix<double,F,A> AxB( m1.size1(), m2.size2() );
       gsl_matrix_view mC = gsl_matrix_view_array (&AxB(0,0), AxB.size1(), AxB.size2());

#ifdef NDEBUG
       (void)gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                  1.0, &mA.matrix, &mB.matrix,
                  0.0, &mC.matrix);
#else
       int status = gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                  1.0, &mA.matrix, &mB.matrix,
                  0.0, &mC.matrix);
       assert(status == 0);
#endif

    #ifdef DEBUG_LINALG
        auto end = std::chrono::system_clock::now();

        std::chrono::duration<double> elapsed_seconds = end - start;

        std::cout << "\n... ... ... \x1b[0;34mGSL "
                  << "[" << m1.size1() << "x" << m1.size2() << "]"
                  << "*[" << m2.size1() << "x" << m2.size2() << "]"
                  << " time: " << elapsed_seconds.count()
                  << "s\x1b[0;39m";
    #endif

    return AxB;

}

    // specialization for single precision (sgemm)
    template<class F, class A>
    inline matrix<float,F,A>    // prod( m1, m2 )
    prod(const matrix<float,F,A> &m1, const matrix<float,F,A> &m2)
    {
    #ifdef DEBUG_LINALG
        auto start = std::chrono::system_clock::now();
    #endif

       gsl_matrix_float_const_view mA = gsl_matrix_float_const_view_array (&m1(0,0), m1.size1(), m1.size2());
       gsl_matrix_float_const_view mB = gsl_matrix_float_const_view_array (&m2(0,0), m2.size1(), m2.size2());

       boost::numeric::ublas::matrix<float,F,A> AxB( m1.size1(), m2.size2() );
       gsl_matrix_float_view mC = gsl_matrix_float_view_array (&AxB(0,0), AxB.size1(), AxB.size2());

       gsl_blas_sgemm (CblasNoTrans, CblasNoTrans,
                  1.0, &mA.matrix, &mB.matrix,
                  0.0, &mC.matrix);

    #ifdef DEBUG_LINALG
        auto end = std::chrono::system_clock::now();
        std::chrono::duration<double> elapsed_seconds = end - start;
        std::cout << "\n... ... ... \x1b[0;34mGSL "
                  << "[" << m1.size1() << "x" << m1.size2() << "]"
                  << "*[" << m2.size1() << "x" << m2.size2() << "]"
                  << " time: " << elapsed_seconds.count()
                  << "s\x1b[0;39m";
    #endif

       return AxB;
    }

}}}

#endif  // __VOTCA_CTP_GSL_H

max_stars_repo_name: jimbach/ctp

max_issues_repo_name: jimbach/ctp

max_forks_repo_name: jimbach/ctp


text:
#ifndef TRANS_UTIL
#define TRANS_UTIL

#define prefetch(x) __builtin_prefetch(x)
#define prefetch_for_r(x) __builtin_prefetch(x, 0, 3)
#define prefetch_for_w(x) __builtin_prefetch(x, 1, 3)

#include <stdarg.h>
#include <gsl/gsl_vector.h>
#include "../include/type.h"

int setulb_(integer *, integer *, doublereal *,
            doublereal *, doublereal *, integer *, doublereal *, doublereal *,
            doublereal *, doublereal *, doublereal *, integer *, char *,
            integer *, char *, logical *, integer *, doublereal *, integer *);

double cuscumsum(gsl_vector *vec, double (*fun) (double, va_list), int argc, ...);

int parse_title(char* str, char** output);

int str_to_vector(char* str, gsl_vector** vec);

int parse_line(char* str, char** id, size_t* id_n, gsl_vector** inc1,
               gsl_vector** skp1, gsl_vector** inc2, gsl_vector** skp2,
               int* inclu_len, int* skip_len);

int parse_file(const char* filename, diff_list_node* list, char** title);

double logit(double i);

double cumprod(const gsl_vector* vec);

double sum_for_multivar(const double i, va_list argv);

double sum_for_multivar_der(const double i, va_list argv);

double myfunc_marginal_2der(const double x, const double I, const double S,
                            const double beta, const double var,
                            const int inclu_len, const int skip_len);

double sum_for_marginal(const double i, va_list argv);

double sum_for_marginal_der(const double i, va_list argv);

odiff* diff_alloc(gsl_vector* inc1, gsl_vector* inc2,
                  gsl_vector* skp1, gsl_vector* skp2,
                  int inclu_len, int skip_len, int flag, char* id);

int diff_append(diff_list_node* header, odiff* data);

void mp_threadpool(int nthread, int ntask, void* (*func)(void *), void** datum, void **ret);

double l_bfgs_b_wrapper(integer n, integer m, doublereal x[], doublereal l[],
                        doublereal u[], integer nbd[],
                        double (*fp) (const double x[], va_list argv),
                        void (*gp) (const double x[], double res[], va_list argv),
                        doublereal factr, doublereal pgtol, integer iprint,
                        int maxfun, int maxiter, int argc, ...);

#endif

max_stars_repo_name: yharlne/rmats-turbo

max_issues_repo_name: yharlne/rmats-turbo

max_forks_repo_name: yharlne/rmats-turbo


text:
// MIT License
//
// Copyright (c) 2020 SunnyCase
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
#pragma once
#include <chino/ddk/kernel.h>
#include <cstdint>
#include <gsl/gsl-lite.hpp>

namespace chino::arch
{
// No asynchronous interrupt is supported by this target
// Only save callee saved registers in context
struct armv7m_thread_context
{
    uintptr_t r4;
    uintptr_t r5;
    uintptr_t r6;
    uintptr_t r7;
    uintptr_t r8;
    uintptr_t r9;
    uintptr_t r10;
    uintptr_t r11;
    uintptr_t lr;
    uintptr_t sp;
    uintptr_t stack_bottom;
};

using thread_context_t = armv7m_thread_context;

struct armv7m_arch
{
    static constexpr size_t ALLOCATE_ALIGNMENT = 8;

    static uint32_t current_processor() noexcept { return 0; }
    static void yield_processor() noexcept;

    static uintptr_t disable_irq() noexcept;
    static void restore_irq(uintptr_t state) noexcept;

    static void init_thread_context(thread_context_t &context, gsl::span<uintptr_t> stack, kernel::thread_thunk_t start, void *arg0, void *arg1) noexcept;
    [[noreturn]] static void start_schedule(thread_context_t &context) noexcept;
    static void yield(thread_context_t &old_context, thread_context_t &new_context) noexcept;

    static void init_stack_check() noexcept;
};

using arch_t = armv7m_arch;
}

max_stars_repo_name: chino-os/chino-os

max_issues_repo_name: dotnetGame/chino-os

max_forks_repo_name: dotnetGame/chino-os


text:
#ifndef UTIL_H
#define UTIL_H

#include <iostream>
#include <stdio.h>
#include <vector>
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/Householder>
#include <assert.h>
#include <time.h>
#include <sys/time.h>
#include <random>
#include <atomic>
#ifdef USE_MKL
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#include <mkl_lapack.h>
#else
#include <cblas.h>
#include <lapacke.h>
#endif

typedef Eigen::SparseMatrix<double, 0, int> SpMat;
typedef Eigen::VectorBlock<Eigen::Matrix<double, -1,  1, 0, -1,  1>, -1> Segment;
typedef Eigen::Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, -1, -1> MatrixBlock;

bool are_connected(Eigen::VectorXi &a, Eigen::VectorXi &b, SpMat &A);
bool should_be_disconnected(int lvl1, int lvl2, int sep1, int sep2);
double elapsed(timeval& start, timeval& end);
void swap2perm(Eigen::VectorXi* swap, Eigen::VectorXi* perm);
bool isperm(Eigen::VectorXi* perm);
SpMat symmetric_graph(SpMat& A);

typedef timeval timer;
timer wctime();

/* Concatenate the matrices vertically in the vector H */
void concatenate(std::vector<Eigen::MatrixXd*> H, Eigen::MatrixXd* V);
Eigen::MatrixXd Vconcatenate(std::vector<Eigen::MatrixXd*> H);
void rVconcatenate(std::vector<Eigen::MatrixXd*> H, Eigen::MatrixXd& Hc);

Eigen::VectorXd Vconcatenate(std::vector<Eigen::VectorXd*> H);

/* Concatenate the matrices horizontally in the vector H */
Eigen::MatrixXd Hconcatenate(std::vector<Eigen::MatrixXd*> H);
/* Reverse the concatenation*/
void rHconcatenate(std::vector<Eigen::MatrixXd*> H, Eigen::MatrixXd& Hc);

/**
 * C <- alpha A   * B   + beta C
 * C <- alpha A^T * B   + beta C
 * C <- alpha A   * B^T + beta C
 * C <- alpha A^T * B^T + beta C
 * Gemm
 */
void gemm(Eigen::MatrixXd* A, Eigen::MatrixXd* B, Eigen::MatrixXd* C, CBLAS_TRANSPOSE tA, CBLAS_TRANSPOSE tB, double alpha, double beta);

/** Return a new
 * C <- alpha A^(/T) * B^(/T)
 **/
Eigen::MatrixXd* gemm_new(Eigen::MatrixXd* A, Eigen::MatrixXd* B, CBLAS_TRANSPOSE tA, CBLAS_TRANSPOSE tB, double alpha);

/**
 * C <- C - A * A^T
 */
void syrk(Eigen::MatrixXd* A, Eigen::MatrixXd* C);

/**
 * A <- L, L L^T = A
 * Return != 0 if potf failed (not spd)
 */
int potf(Eigen::MatrixXd* A);

/**
 * A <- [L\U] (lower and upper)
 * p <- swap (NOT a permutation)
 * A[p] = L*U
 * L is unit diagonal
 * U is not
 * Return != 0 if getf failed (singular)
 */
int getf(Eigen::MatrixXd* A, Eigen::VectorXi* swap);

/**
 * Compute an estimated 1-norm condition number of A using its LU or Cholesky factorization
 */
double rcond_1_getf(Eigen::MatrixXd* A_LU, double A_1_norm);
double rcond_1_potf(Eigen::MatrixXd* A_LLT, double A_1_norm);

/**
 * B <- B * L^(-1)
 * B <- B * L^(-T)
 * B <- B * U^(-1)
 * B <- B * U^(-T)
 */
void trsm_right(Eigen::MatrixXd* L, Eigen::MatrixXd* B, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans, CBLAS_DIAG diag);

/**
 * B <- L^(-1) * B
 * B <- L^(-T) * B
 * B <- U^(-1) * B
 * B <- U^(-T) * B
 */
void trsm_left(Eigen::MatrixXd* L, Eigen::MatrixXd* B, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans, CBLAS_DIAG diag);

/**
 * x <- L^(-1) * x
 * x <- L^(-T) * x
 * x <- U^(-1) * x
 * x <- U^(-T) * x
 */
void trsv(Eigen::MatrixXd* L, Segment* x, CBLAS_UPLO uplo, CBLAS_TRANSPOSE trans, CBLAS_DIAG diag);

/**
 * x <- L^T * x
 */
void trmv_trans(Eigen::MatrixXd* L, Segment* x);

/**
 * A <- L^T * A
 */
void trmm_trans(Eigen::MatrixXd* L, Eigen::MatrixXd* A);

/**
 * x2 <- x2 - A21 * x1
 */
void gemv_notrans(Eigen::MatrixXd* A21, Segment* x1, Segment* x2);

/**
 * x2 <- x2 - A12^T * x1
 */
void gemv_trans(Eigen::MatrixXd* A12, Segment* x1, Segment* x2);

/**
 * AP = QR
 */
void geqp3(Eigen::MatrixXd* A, Eigen::VectorXi* jpvt, Eigen::VectorXd* tau);

/**
 * Form householder vector from x
 */
double house(Eigen::VectorXd& x);

/**
 * RRQR with truncation when the max R_ii < tol
 */
void rrqr(Eigen::MatrixXd* A, Eigen::VectorXi* jpvt, Eigen::VectorXd* tau, double& tol, int& rank);

// template<typename T>
// Eigen::MatrixXd get_gaussian(const int rows, const int cols, T* gen);
template<typename T>
Eigen::MatrixXd get_gaussian(const int rows, const int cols, T* gen) {
    std::normal_distribution<double> norm_dist(0.0, 1.0);
    Eigen::MatrixXd W(rows, cols);
    for(int j = 0; j < cols; j++){
        for(int i = 0; i < rows; i++){
            W(i,j) = norm_dist(*gen);
        }
    }
    Eigen::VectorXd col_norms = W.colwise().norm();
    assert(col_norms.size() == cols);
    assert(col_norms.minCoeff() >= 0);
    if (col_norms.minCoeff() > 0) {
        return W.normalized();
    }
    return W;
}

Eigen::MatrixXd get_uniform(const int rows, const int cols);
void random_rrqr(const Eigen::MatrixXd* A, Eigen::MatrixXd* v, Eigen::VectorXd* h,  double tol, std::function<Eigen::MatrixXd(int,int)> gen_gaussian);

// void laqps(Eigen::MatrixXd* A, Eigen::VectorXi* jpvt, Eigen::VectorXd* tau, double& tol, int block_size, int& rank);
void laqps(Eigen::MatrixXd* A, Eigen::VectorXi* jpvt, Eigen::VectorXd* tau, double& tol, int block_size, int& rank);
/**
 * A = U S U^T
 * Compute the full EVD, where if A is mxm, U is mxm and S is mxm
 */

void geevd(Eigen::MatrixXd* A, Eigen::MatrixXd* U, Eigen::VectorXd* S);

/**
 * A = U S VT
 * Compute the full SVD, where if A is mxn, U is mxm, V is nxn, and S is min(M,N)
 * VT is V^T, *not* V.
 */
void gesvd(Eigen::MatrixXd* A, Eigen::MatrixXd* U, Eigen::VectorXd* S, Eigen::MatrixXd* VT);

/**
 * Compute the singular values of a matrix A
 * A = m x n matrix
 * S = min(M,N) vector
 */
void gesvd_values(Eigen::MatrixXd* A,  Eigen::VectorXd* S);

/**
 * x <- Q * x
 * A <- Q * A
 */
void ormqr_notrans(Eigen::MatrixXd* v, Eigen::VectorXd* h, Segment* x);
void ormqr_notrans_left(Eigen::MatrixXd* v, Eigen::VectorXd* h, Eigen::MatrixXd* A);

/**
 * A <- A * Q
 * */
void ormqr_notrans_right(Eigen::MatrixXd* v, Eigen::VectorXd* h, Eigen::MatrixXd* A);
/**
 * x <- Q^T * x
 * A <- Q^T * A
 */
void ormqr_trans(Eigen::MatrixXd* v, Eigen::VectorXd* h, Segment* x);
void ormqr_trans_left(Eigen::MatrixXd* v, Eigen::VectorXd* h, Eigen::MatrixXd* A);

/**
 * A <- A   * Q
 * A <- A   * Q^T
 * A <- Q   * A
 * A <- Q^T * A
 */
void ormqr(Eigen::MatrixXd* v, Eigen::VectorXd* h, Eigen::MatrixXd* A, char side, char trans);

/**
 * Create the thin Q
 */
void orgqr(Eigen::MatrixXd* v, Eigen::VectorXd* h);

/* Finds upper triangular matrix T such that,
    H  =  I - V * T * V^T
    using Householder vectors V and their norms tau
    H = Householder reflector matrix
*/
void larft(Eigen::MatrixXd* V, Eigen::VectorXd* tau, Eigen::MatrixXd* T);

/* Apply householder vectors on a rectangular matrix
V = [1    *     *   * *
     v(1) 1     *   * *
     v(1) v(2)  *   * *
     v(1) v(2) v(3) * *];
H = H(1) H(2) H(3) ... H(k)
C <- H^T C
*/
void larfb(Eigen::MatrixXd* V, Eigen::MatrixXd* T, Eigen::MatrixXd* C);
void larfb(Eigen::MatrixXd* V, Eigen::MatrixXd* T, Eigen::VectorXd* C, char side = 'L', char trans = 'T', char direct = 'F', char storev = 'C');

/**
 * A = QR
 */
void geqrf(Eigen::MatrixXd* A, Eigen::VectorXd* tau);

int choose_rank(Eigen::VectorXd& s, double tol, bool rel = true);

std::size_t hashv(std::vector<size_t> vals);

// Hash function for Eigen matrix and vector.
// The code is from `hash_combine` function of the Boost library. See
// http://www.boost.org/doc/libs/1_55_0/doc/html/hash/reference.html#boost.hash_combine .
template<typename T>
struct matrix_hash : std::unary_function<T, size_t> {
  std::size_t operator()(T const& matrix) const {
    // Note that it is oblivious to the storage order of Eigen matrix (column- or
    // row-major). It will give you the same hash value for two different matrices if they
    // are the transpose of each other in different storage order.
    size_t seed = 0;
    for (size_t i = 0; i < matrix.size(); ++i) {
      auto elem = *(matrix.data() + i);
      seed ^= std::hash<typename T::Scalar>()(elem) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
    }
    return seed;
  }
};

void block2dense(Eigen::VectorXi &rowval, Eigen::VectorXi &colptr, Eigen::VectorXd &nnzval, int i, int j, int li, int lj, Eigen::MatrixXd *dst, bool transpose);

Eigen::MatrixXd linspace_nd(int n, int dim);

// Returns A[p,p]
SpMat symm_perm(SpMat &A, Eigen::VectorXi &p);

// Permute the columns of non-square matrix A
SpMat col_perm(SpMat &A, Eigen::VectorXi &p);

// Random vector with seed
Eigen::VectorXd random(int size, int seed);

Eigen::MatrixXd random(int rows, int cols, int seed);

// Print vector
template<typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& v) {
    for(auto v_ : v) {
        os << v_ << " " ;
    }
    os << std::endl;
    return os;
}

#endif