text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

int
main (void)
{
  int i;
  double xi, yi, x[10], y[10];

  printf ("#m=0,S=17\n");

  for (i = 0; i < 10; i++)
    {
      x[i] = i + 0.5 * sin (i);
      y[i] = i + cos (i * i);
      printf ("%g %g\n", x[i], y[i]);
    }

  printf ("#m=1,S=0\n");

  {
    gsl_interp_accel *acc
      = gsl_interp_accel_alloc ();
    gsl_spline *spline
      = gsl_spline_alloc (gsl_interp_cspline, 10);

    gsl_spline_init (spline, x, y, 10);

    for (xi = x[0]; xi < x[9]; xi += 0.01)
      {
        yi = gsl_spline_eval (spline, xi, acc);
        printf ("%g %g\n", xi, yi);
      }
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);
  }
  return 0;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
 * Copyright (c) 2009-2020 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 *
 */
#ifndef _TESTSCOMMON_H
#define _TESTSCOMMON_H

/* parsec things */
#include "parsec.h"

/* system and io */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
/* Plasma and math libs */
#include <math.h>
#include <cblas.h>
#include <lapacke.h>

#include "parsec/profiling.h"
#include "parsec/parsec_internal.h"
#include "parsec/utils/debug.h"
#include "dplasma.h"
#include "dplasmatypes.h"

#include <core_blas.h>
/* timings */
#include "common_timing.h"

#include "flops.h"

/* these are globals in common.c */
extern char *PARSEC_SCHED_NAME[];
extern int unix_timestamp;
extern char cwd[];

/* Update PASTE_CODE_PROGRESS_KERNEL below if you change this list */
enum iparam_t {
  IPARAM_RANK,         /* Rank                              */
  IPARAM_NNODES,       /* Number of nodes                   */
  IPARAM_NCORES,       /* Number of cores                   */
  IPARAM_NGPUS,        /* Number of GPUs                    */
  IPARAM_P,            /* Rows in the process grid          */
  IPARAM_Q,            /* Columns in the process grid       */
  IPARAM_M,            /* Number of rows of the matrix      */
  IPARAM_N,            /* Number of columns of the matrix   */
  IPARAM_K,            /* RHS or K                          */
  IPARAM_LDA,          /* Leading dimension of A            */
  IPARAM_LDB,          /* Leading dimension of B            */
  IPARAM_LDC,          /* Leading dimension of C            */
  IPARAM_IB,           /* Inner-blocking size               */
  IPARAM_NB,           /* Number of columns in a tile       */
  IPARAM_MB,           /* Number of rows in a tile          */
  IPARAM_KQ,           /* Number of columns repetitions in a k-cyclic distribution */
  IPARAM_KP,           /* Number of rows repetititions in a k-cyclic distribution */
  IPARAM_HMB,          /* Small MB for recursive hdags */
  IPARAM_HNB,          /* Small NB for recursive hdags */
  IPARAM_CHECK,        /* Checking activated or not         */
  IPARAM_CHECKINV,     /* Inverse Checking activated or not */
  IPARAM_ASYNC,        /* Bench the asynchronous version    */
  IPARAM_VERBOSE,      /* How much noise do we want?        */
  IPARAM_LOWLVL_TREE,  /* Tree used for reduction inside nodes  (specific to xgeqrf_param) */
  IPARAM_HIGHLVL_TREE, /* Tree used for reduction between nodes (specific to xgeqrf_param) */
  IPARAM_QR_TS_SZE,    /* Size of TS domain                     (specific to xgeqrf_param) */
  IPARAM_QR_HLVL_SZE,  /* Size of the high level tree           (specific to xgeqrf_param) */
  IPARAM_RANDOM_SEED,  /* Seed for the pseudo-random generators */
  IPARAM_MATRIX_INIT,  /* Matrix generator type */
  IPARAM_QR_DOMINO,    /* Enable/disable the domino between the upper and the lower tree (specific to xgeqrf_param) */
  IPARAM_QR_TSRR,      /* Enable/disable the round-robin on TS domain */
  IPARAM_BUT_LEVEL,    /* Butterfly level */
  IPARAM_SCHEDULER,    /* User-selected scheduler */
  IPARAM_SIZEOF
};

#define PARSEC_SCHEDULER_DEFAULT 0
#define PARSEC_SCHEDULER_LFQ 1
#define PARSEC_SCHEDULER_LTQ 2
#define PARSEC_SCHEDULER_AP  3
#define PARSEC_SCHEDULER_LHQ 4
#define PARSEC_SCHEDULER_GD  5
#define PARSEC_SCHEDULER_PBQ 6
#define PARSEC_SCHEDULER_IP  7
#define PARSEC_SCHEDULER_RND 8

void iparam_default_facto(int* iparam);
void iparam_default_solve(int* iparam);
void iparam_default_gemm(int* iparam);
void iparam_default_ibnbmb(int* iparam, int ib, int nb, int mb);

#define PASTE_CODE_IPARAM_LOCALS(iparam)                                \
    int rank  = iparam[IPARAM_RANK];                                    \
    int nodes = iparam[IPARAM_NNODES];                                  \
    int cores = iparam[IPARAM_NCORES];                                  \
    int gpus  = iparam[IPARAM_NGPUS];                                   \
    int P     = iparam[IPARAM_P];                                       \
    int Q     = iparam[IPARAM_Q];                                       \
    int M     = iparam[IPARAM_M];                                       \
    int N     = iparam[IPARAM_N];                                       \
    int K     = iparam[IPARAM_K];                                       \
    int NRHS  = K;                                                      \
    int LDA   = max(M, iparam[IPARAM_LDA]);                             \
    int LDB   = max(N, iparam[IPARAM_LDB]);                             \
    int LDC   = max(K, iparam[IPARAM_LDC]);                             \
    int IB    = iparam[IPARAM_IB];                                      \
    int MB    = iparam[IPARAM_MB];                                      \
    int NB    = iparam[IPARAM_NB];                                      \
    int SMB   = iparam[IPARAM_KP];                                      \
    int SNB   = iparam[IPARAM_KQ];                                      \
    int HMB   = iparam[IPARAM_HMB];                                     \
    int HNB   = iparam[IPARAM_HNB];                                     \
    int MT    = (M%MB==0) ? (M/MB) : (M/MB+1);                          \
    int NT    = (N%NB==0) ? (N/NB) : (N/NB+1);                          \
    int KT    = (K%MB==0) ? (K/MB) : (K/MB+1);                          \
    int check = iparam[IPARAM_CHECK];                                   \
    int check_inv = iparam[IPARAM_CHECKINV];                            \
    int loud  = iparam[IPARAM_VERBOSE];                                 \
    int scheduler = iparam[IPARAM_SCHEDULER];                           \
    int random_seed = iparam[IPARAM_RANDOM_SEED];                       \
    int matrix_init = iparam[IPARAM_MATRIX_INIT];                       \
    int butterfly_level = iparam[IPARAM_BUT_LEVEL];                     \
    int async = iparam[IPARAM_ASYNC];                                   \
    (void)rank;(void)nodes;(void)cores;(void)gpus;(void)P;(void)Q;(void)M;(void)N;(void)K;(void)NRHS; \
    (void)LDA;(void)LDB;(void)LDC;(void)IB;(void)MB;(void)NB;(void)MT;(void)NT;(void)KT; \
    (void)SMB;(void)SNB;(void)HMB;(void)HNB;(void)check;(void)loud;(void)async; \
    (void)scheduler;(void)butterfly_level;(void)check_inv;(void)random_seed;(void)matrix_init;

/* Define a double type which not pass through the precision generation process */
typedef double DagDouble_t;
#define PASTE_CODE_FLOPS( FORMULA, PARAMS ) \
  double gflops = -1.0, flops = FORMULA PARAMS;

#if defined(PRECISION_z) || defined(PRECISION_c)
#define PASTE_CODE_FLOPS_COUNT(FADD,FMUL,PARAMS) \
  double gflops = -1.0, flops = (2. * FADD PARAMS + 6. * FMUL PARAMS);
#else
#define PASTE_CODE_FLOPS_COUNT(FADD,FMUL,PARAMS) \
  double gflops = -1.0, flops = (FADD PARAMS + FMUL PARAMS);
#endif

/*******************************
 * globals values
 *******************************/

#if defined(PARSEC_HAVE_MPI)
extern MPI_Datatype SYNCHRO;
#endif  /* PARSEC_HAVE_MPI */

extern const int side[2];
extern const int uplo[2];
extern const int diag[2];
extern const int trans[3];
extern const int norms[4];
extern const char *sidestr[2];
extern const char *uplostr[2];
extern const char *diagstr[2];
extern const char *transstr[3];
extern const char *normsstr[4];

void print_usage(void);

parsec_context_t *setup_parsec(int argc, char* argv[], int *iparam);
void cleanup_parsec(parsec_context_t* parsec, int *iparam);

/**
 * No macro with the name max or min is acceptable as there is
 * no way to correctly define them without borderline effects.
 */
#undef max
#undef min
static inline int max(int a, int b) { return a > b ? a : b; }
static inline int min(int a, int b) { return a < b ? a : b; }

/* Paste code to allocate a matrix in desc if cond_init is true */
#define PASTE_CODE_ALLOCATE_MATRIX(DC, COND, TYPE, INIT_PARAMS)      \
    TYPE##_t DC;                                                     \
    if(COND) {                                                          \
        TYPE##_init INIT_PARAMS;                                        \
        DC.mat = parsec_data_allocate((size_t)DC.super.nb_local_tiles * \
                                        (size_t)DC.super.bsiz *      \
                                        (size_t)parsec_datadist_getsizeoftype(DC.super.mtype)); \
        parsec_data_collection_set_key((parsec_data_collection_t*)&DC, #DC);          \
    }

#define PASTE_CODE_ENQUEUE_KERNEL(PARSEC, KERNEL, PARAMS)               \
    SYNC_TIME_START();                                                  \
    parsec_taskpool_t* PARSEC_##KERNEL = dplasma_##KERNEL##_New PARAMS;   \
    PARSEC_CHECK_ERROR(parsec_context_add_taskpool(PARSEC, PARSEC_##KERNEL), "parsec_context_add_taskpool");      \
    if( loud > 2 ) SYNC_TIME_PRINT(rank, ( #KERNEL "\tDAG created\n"));

#define PASTE_CODE_PROGRESS_KERNEL(PARSEC, KERNEL)                      \
    SYNC_TIME_START();                                                  \
    PARSEC_CHECK_ERROR(parsec_context_start(PARSEC), "parsec_context_start"); \
    TIME_START();                                                       \
    PARSEC_CHECK_ERROR(parsec_context_wait(PARSEC), "parsec_context_wait"); \
    SYNC_TIME_PRINT(rank, (#KERNEL "\tPxQ= %3d %-3d NB= %4d N= %7d : %14f gflops\n", \
                           P, Q, NB, N,                                 \
                           gflops=(flops/1e9)/sync_time_elapsed));      \
    PROFILING_SAVE_dINFO("TIME_ELAPSED", time_elapsed);                 \
    PROFILING_SAVE_dINFO("SYNC_TIME_ELAPSED", sync_time_elapsed);       \
    PROFILING_SAVE_dINFO("GFLOPS", gflops);                             \
    PROFILING_SAVE_iINFO("PARAM_RANK", iparam[IPARAM_RANK]);            \
    PROFILING_SAVE_iINFO("PARAM_NNODES", iparam[IPARAM_NNODES]);        \
    PROFILING_SAVE_iINFO("PARAM_NCORES", iparam[IPARAM_NCORES]);        \
    PROFILING_SAVE_iINFO("PARAM_NGPUS", iparam[IPARAM_NGPUS]);          \
    PROFILING_SAVE_iINFO("PARAM_P", iparam[IPARAM_P]);                  \
    PROFILING_SAVE_iINFO("PARAM_Q", iparam[IPARAM_Q]);                  \
    PROFILING_SAVE_iINFO("PARAM_M", iparam[IPARAM_M]);                  \
    PROFILING_SAVE_iINFO("PARAM_N", iparam[IPARAM_N]);                  \
    PROFILING_SAVE_iINFO("PARAM_K", iparam[IPARAM_K]);                  \
    PROFILING_SAVE_iINFO("PARAM_LDA", iparam[IPARAM_LDA]);              \
    PROFILING_SAVE_iINFO("PARAM_LDB", iparam[IPARAM_LDB]);              \
    PROFILING_SAVE_iINFO("PARAM_LDC", iparam[IPARAM_LDC]);              \
    PROFILING_SAVE_iINFO("PARAM_IB", iparam[IPARAM_IB]);                \
    PROFILING_SAVE_iINFO("PARAM_NB", iparam[IPARAM_NB]);                \
    PROFILING_SAVE_iINFO("PARAM_MB", iparam[IPARAM_MB]);                \
    PROFILING_SAVE_iINFO("PARAM_KP", iparam[IPARAM_KQ]);                \
    PROFILING_SAVE_iINFO("PARAM_KQ", iparam[IPARAM_KP]);                \
    PROFILING_SAVE_iINFO("PARAM_HNB", iparam[IPARAM_HNB]);              \
    PROFILING_SAVE_iINFO("PARAM_CHECK", iparam[IPARAM_CHECK]);          \
    PROFILING_SAVE_iINFO("PARAM_CHECKINV", iparam[IPARAM_CHECKINV]);    \
    PROFILING_SAVE_iINFO("PARAM_VERBOSE", iparam[IPARAM_VERBOSE]);      \
    PROFILING_SAVE_iINFO("PARAM_LOWLVL_TREE", iparam[IPARAM_LOWLVL_TREE]); \
    PROFILING_SAVE_iINFO("PARAM_HIGHLVL_TREE", iparam[IPARAM_HIGHLVL_TREE]); \
    PROFILING_SAVE_iINFO("PARAM_QR_TS_SZE", iparam[IPARAM_QR_TS_SZE]);  \
    PROFILING_SAVE_iINFO("PARAM_QR_HLVL_SZE", iparam[IPARAM_QR_HLVL_SZE]); \
    PROFILING_SAVE_iINFO("PARAM_QR_DOMINO", iparam[IPARAM_QR_DOMINO]);  \
    PROFILING_SAVE_iINFO("PARAM_QR_TSRR", iparam[IPARAM_QR_TSRR]);      \
    PROFILING_SAVE_iINFO("PARAM_BUT_LEVEL", iparam[IPARAM_BUT_LEVEL]);  \
    PROFILING_SAVE_iINFO("PARAM_SCHEDULER", iparam[IPARAM_SCHEDULER]);  \
    if(loud >= 5 && rank == 0) {                                        \
        printf("<DartMeasurement name=\"performance\" type=\"numeric/double\"\n" \
               "                 encoding=\"none\" compression=\"none\">\n" \
               "%g\n"                                                   \
               "</DartMeasurement>\n",                                  \
               gflops);                                                 \
    }                                                                   \
    (void)gflops;

#endif /* _TESTSCOMMON_H */

ext: h

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
#pragma once

#include "halley/text/halleystring.h"
#include "halley/utils/utils.h"
#include <gsl/gsl>
#include <array>

namespace Halley {
	class Deserializer;
	class Serializer;

	class UUID {
    public:
        UUID();
        UUID(std::array<Byte, 16> bytes);
        explicit UUID(gsl::span<const gsl::byte> bytes);
		explicit UUID(const Bytes& bytes);
        explicit UUID(const String& str);

        bool operator==(const UUID& other) const;
        bool operator!=(const UUID& other) const;
		bool operator<(const UUID& other) const;

		String toString() const;

        static UUID generate();
        static UUID generateFromUUIDs(const UUID& one, const UUID& two);
    	bool isValid() const;

        gsl::span<const gsl::byte> getBytes() const;
		gsl::span<gsl::byte> getWriteableBytes();
        gsl::span<const uint64_t> getUint64Bytes() const;

    	void serialize(Serializer& s) const;
		void deserialize(Deserializer& s);

    private:
        std::array<uint64_t, 2> bytes;
    };
}

template<>
struct std::hash<Halley::UUID>
{
    std::size_t operator() (const Halley::UUID& uuid) const noexcept
    {
        const auto& bytes = uuid.getUint64Bytes();
        return bytes[0] ^ bytes[1];
    }
};

namespace natvis {
    struct x4lo {
    	uint8_t v: 4;
    	uint8_t _: 4;
    };
    struct x4hi {
	    uint8_t _ : 4;
    	uint8_t v : 4;
    };
    struct x8 {
	    uint8_t _;
    };
	struct x16 {
	    uint16_t _;
    };
    struct x32 {
	    uint32_t _;
    };
	struct x48 {
	    uint32_t b0;
		uint16_t b1;
    };
}

max_stars_repo_name: code-disaster/halley

max_issues_repo_name: code-disaster/halley

max_forks_repo_name: code-disaster/halley


text:
/// @file
///
/// Class to handle Taylor-term fitting to a spectrum
///
/// @copyright (c) 2017 CSIRO
/// Australia Telescope National Facility (ATNF)
/// Commonwealth Scientific and Industrial Research Organisation (CSIRO)
/// PO Box 76, Epping NSW 1710, Australia
/// atnf-enquiries@csiro.au
///
/// This file is part of the ASKAP software distribution.
///
/// The ASKAP software distribution 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
///
/// @author Matthew Whiting <Matthew.Whiting@csiro.au>
///
#ifndef ASKAP_ANALYSISUTILS_FIT_TT_H_
#define ASKAP_ANALYSISUTILS_FIT_TT_H_

#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <casacore/casa/Arrays/Array.h>

namespace askap {

namespace analysisutilities {

/// Data structure used within the fitting as a way of accessing the
/// various arrays within the GSL functions.
struct data {
    size_t ndata;
    double * xdat;
    double * ydat;
    double * weights;
};

/// @brief Class to perform non-linear fitting to find alpha & beta.
/// @details A non-linear function describing a
/// power-law-with-curvature function is fitted to a spectrum using
/// non-linear Levenberg-Marquardt algorithm. The function is F(nu) =
/// F_0 (nu/nu_0) ^ (alpha + beta * log(nu/nu_0)), where F_0 is the
/// flux at the reference frequency nu_0. It is a weighted fit,
/// taking into account a provided noise array. The number of terms is
/// user-selectable, from 1 (constant flux), 2 (constant flux +
/// spectral index) or 3 (add spectral curvature).
class FitTT {
    public:
        /// Initialise, setting the number of terms that should be fit.
        FitTT(unsigned int nterms);
        virtual ~FitTT() {};

        /// @brief Fit to provided arrays.
        /// @details For spectral-index fitting, x should be the
        /// frequencies divided by the reference frequency, and y should
        /// be the flux values. No weights array is provided, so weights
        /// are implicitly set to 1.  Front end to fit(size_t ndata,
        /// double *xdata, double *ydata, double *weights).
        void fit(casa::Array<float> &x, casa::Array<float> &y);

        /// @brief Fit to provided arrays, with associated weights array.
        /// @details For spectral-index fitting, x should be the
        /// frequencies divided by the reference frequency, y should be
        /// the flux values, and w the weights on the flux values
        /// (ie. noise). Front end to fit(size_t ndata, double *xdata,
        /// double *ydata, double *weights).
        void fit(casa::Array<float> &x, casa::Array<float> &y, casa::Array<float> &w);

        /// @brief Fit to provided arrays - main function that does the fitting.
        /// @details For spectral-index fitting, xdata should be the
        /// frequencies divided by the reference frequency, ydata should be
        /// the flux values, and weights the weights on the flux values
        /// (ie. noise), with ndata the size of the arrays.
        void fit(size_t ndata, double *xdata, double *ydata, double *weights);

        /// @brief Model function, for nterms=1
        static int taylor_f1(const gsl_vector * p, void *data,  gsl_vector * f);
        /// @brief Jacobian function, for nterms=1
        static int taylor_df1(const gsl_vector * p, void *data,  gsl_matrix * J);
        /// @brief Combined model & Jacobian function, for nterms=1
        static int taylor_fdf1(const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J);
        /// @brief Model function, for nterms=2
        static int taylor_f2(const gsl_vector * p, void *data,  gsl_vector * f);
        /// @brief Jacobian function, for nterms=2
        static int taylor_df2(const gsl_vector * p, void *data,  gsl_matrix * J);
        /// @brief Combined model & Jacobian function, for nterms=2
        static int taylor_fdf2(const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J);
        /// @brief Model function, for nterms=3
        static int taylor_f3(const gsl_vector * p, void *data,  gsl_vector * f);
        /// @brief Jacobian function, for nterms=3
        static int taylor_df3(const gsl_vector * p, void *data,  gsl_matrix * J);
        /// @brief Combined model & Jacobian function, for nterms=3
        static int taylor_fdf3(const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J);

        /// @brief Return the fitted value of the flux at the reference pixel
        float fluxZero() {return itsFluxZero;};
        /// @brief Return the fitted value of the spectral index
        float alpha() {return itsAlpha;};
        /// @brief Return the fitted value of the spectral curvature
        float beta() {return itsBeta;};

        /// @brief Return the error on the fitted value of the flux at the reference frequency
        float fluxZeroErr() {return itsFluxZeroErr;};
        /// @brief Return the error on the fitted value of the spectral index
        float alphaErr() {return itsAlphaErr;};
        /// @brief Return the error on the fitted value of the spectral curvature
        float betaErr() {return itsBetaErr;};

    protected:
        /// @brief How many terms to fit?
        unsigned int itsNterms;

        /// @brief Flux at reference frequency
        float itsFluxZero;
        /// @brief Spectral index
        float itsAlpha;
        /// @brief Spectral curvature
        float itsBeta;
        /// @brief Error on flux at reference frequency
        float itsFluxZeroErr;
        /// @brief Error on spectral index
        float itsAlphaErr;
        /// @brief Error on spectral curvature
        float itsBetaErr;

};

}
}
#endif

max_stars_repo_name: steve-ord/askapsoft

max_issues_repo_name: ATNF/askapsoft

max_forks_repo_name: ATNF/askapsoft


text:
/********************************************************
  Stanford Driving Software
  Copyright (c) 2011 Stanford University
  All rights reserved.

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

* Redistributions of source code must retain the above
  copyright notice, this list of conditions and the
  following disclaimer.
* 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.
* The names of the contributors may not be used to endorse
  or promote products derived from this software
  without specific prior written permission.

  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
  CONTRIBUTORS "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 THE
  COPYRIGHT OWNER 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.
 ********************************************************/

#ifndef GSL_SUPPORT_H
#define GSL_SUPPORT_H

#include <gsl/gsl_blas.h>

namespace vlr {

#define matrix_get gsl_matrix_get
#define vector_get gsl_vector_get
#define matrix_set gsl_matrix_set
#define vector_set gsl_vector_set
#define vector_copy gsl_vector_memcpy
#define matrix_copy gsl_matrix_memcpy
#define matrix_set_identity gsl_matrix_set_identity
#define vector_set_zero gsl_vector_set_zero
#define matrix_set_zero gsl_matrix_set_zero
#define vector_scale gsl_vector_scale

gsl_vector *vector_init(int length);
void vector_fill(gsl_vector *v, ...);
void vector_free(gsl_vector *vector);
gsl_matrix *matrix_init(int rows, int cols);
void matrix_fill(gsl_matrix *m, ...);
void matrix_free(gsl_matrix *matrix);
void vector_add(gsl_vector *result, gsl_vector *a, gsl_vector *b);
void vector_add_ip(gsl_vector *a, gsl_vector *b);
void vector_elem_add(gsl_vector *v, int elem, double x);
void vector_subtract(gsl_vector *result, gsl_vector *a,
                                gsl_vector *b);
void vector_subtract_ip(gsl_vector *a, gsl_vector *b);
void matrix_add(gsl_matrix *result, gsl_matrix *a, gsl_matrix *b);
void matrix_add_ip(gsl_matrix *a, gsl_matrix *b);
void matrix_elem_add(gsl_matrix *m, int row, int col, double x);
void matrix_subtract(gsl_matrix *result, gsl_matrix *a,
                                gsl_matrix *b);
void matrix_subtract_ip(gsl_matrix *a, gsl_matrix *b);
void matrix_multiply(gsl_matrix *result, gsl_matrix *a, int a_t,
                     gsl_matrix *b, int b_t);
void matrix_invert(gsl_matrix *inverse, gsl_matrix *m);
void matrix_invert_ip(gsl_matrix *m);
void matrix_vector_multiply(gsl_vector *result, gsl_matrix *a,
                                       gsl_vector *x);
void matrix_linear_transform(gsl_matrix *result, gsl_matrix *m,
                                        gsl_matrix *a);
void matrix_inverse_linear_transform(gsl_matrix *result,
                                                gsl_matrix *m, gsl_matrix *a);
void matrix_scale(gsl_matrix *m, double scale);
void matrix_print(gsl_matrix *m, const char *str);
void vector_print(gsl_vector *v, const char *str);
void matrix_print_signs(gsl_matrix *m, const char *str);
void matrix_to_file(gsl_matrix *m, const char *filename);
void vector_to_file(gsl_vector *v, const char *filename);
double vector_inner_product(gsl_vector *a, gsl_vector *b);
void gaussian_sample(gsl_vector *sample, gsl_vector *mean,
                                gsl_matrix *cov);
double matrix_determinant(gsl_matrix *A);

} // namespace vlr

#endif


text:
/*
**  Matrix Projection,
**   Ref:  Lohmann, Loktyushin, Stelzer, Scheffler (2018) bioRXiv
**
** G.Lohmann, MPI-KYB, Nov 2018
*/

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

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_linalg.h>

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

#ifdef _OPENMP
#include <omp.h>
#endif /*_OPENMP*/

#define SQR(x) ((x) * (x))
#define IMIN(x,y) ((x) < (y) ? (x) : (y))
#define IMAX(x,y) ((x) > (y) ? (x) : (y))

extern double ECMcorrelation(const float *arr1,const float *arr2,size_t nt,int type);
extern void NormVec(float *x,size_t);

void VMatrixProjection(gsl_matrix_float *X,float *ev,int type,int seed)
{
  size_t i,j;
  size_t nvox = X->size1;
  size_t nt = X->size2;
  size_t p = 32;

  fprintf(stderr," VMatrixProjection...\n");

  /* generate random matrix */
  const gsl_rng_type *Tx = gsl_rng_default;
  gsl_rng *rx = gsl_rng_alloc (Tx);
  gsl_rng_set(rx,(unsigned long int)seed);
  gsl_matrix *R = gsl_matrix_calloc(nvox,p);
  for (i=0; i<nvox; i++) {
    for (j=0; j<p; j++) {
      gsl_matrix_set(R,i,j,gsl_ran_ugaussian(rx));
    }
  }

  /* matrix product 1 */
  gsl_matrix *Y = gsl_matrix_calloc(nvox,p);
  size_t progress=0;
  size_t step=4;  /* sparse */
  fprintf(stderr," pass 1:\n");

#pragma omp parallel for shared(progress) schedule(guided) firstprivate(X,Y,R)
  for (i=0; i<nvox; i++) {
    if (i%1000 == 0) fprintf(stderr," %d000  of  %lu\r",(int)(++progress),nvox);

    size_t j=0,k=0;
    double *tmp1 = (double *)VCalloc(nvox,sizeof(double));
    double *tmp2 = NULL;
    double sum=0;
    const float *arr1 = gsl_matrix_float_const_ptr(X,i,0);

    for (j=0; j<nvox; j++) {
      const float *arr2 = gsl_matrix_float_const_ptr(X,j,0);
      tmp1[j] = ECMcorrelation(arr1,arr2,nt,type);
    }

    for (j=0; j<p; j++) {
      tmp2 = gsl_matrix_ptr(R,j,0);
      sum=0;
      for (k=0; k<nvox; k+=step) {
	sum += tmp1[k]*tmp2[k];
      }
      gsl_matrix_set(Y,i,j,sum);
    }
    VFree(tmp1);
  }
  fprintf(stderr,"\n");
  gsl_matrix_free(R);

  /* get orthonormal basis for Y  (SVD) */
  gsl_vector *sv = gsl_vector_calloc(p);
  gsl_vector *work = gsl_vector_calloc(p);
  gsl_matrix *V = gsl_matrix_calloc(p,p);
  gsl_linalg_SV_decomp(Y,V,sv,work);

  /* matrix product 2 */
  fprintf(stderr," pass 2:\n");
  gsl_matrix *D = gsl_matrix_calloc(nvox,p);
  progress=0;

  #pragma omp parallel for shared(progress) schedule(guided) firstprivate(X,Y,D)
  for (i=0; i<nvox; i++) {
    if (i%1000 == 0) fprintf(stderr," %d000  of  %lu\r",(int)(++progress),nvox);

    size_t j=0,k=0;
    double *tmp1 = (double *)VCalloc(nvox,sizeof(double));
    double sum=0;
    const float *arr1 = gsl_matrix_float_const_ptr(X,i,0);

    for (j=0; j<nvox; j++) {
      const float *arr2 = gsl_matrix_float_const_ptr(X,j,0);
      tmp1[j] = ECMcorrelation(arr1,arr2,nt,type);
    }

    for (j=0; j<p; j++) {
      sum=0;
      for (k=0; k<nvox; k+=step) {
	sum += tmp1[k] * gsl_matrix_get(Y,k,j);
      }
      gsl_matrix_set(D,i,j,sum);
    }
    VFree(tmp1);
  }
  fprintf(stderr,"\n");
  gsl_matrix_free(Y);

  /* final SVD */
  gsl_linalg_SV_decomp(D,V,sv,work);

  /* output */
  size_t k=0;
  for (i=0; i<nvox; i++) {
    ev[i] = gsl_matrix_get(D,i,0);
    if (ev[i] < 0) k++;
  }

  /* sign switch if needed */
  if (k > nvox/2) {
    for (i=0; i<nvox; i++) {
      ev[i] = -ev[i];
    }
  }

  /* normalize */
  NormVec(ev,nvox);
  float kx = sqrt((float)(nvox));
  for (i=0; i<nvox; i++)  ev[i] *= kx;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:43:52 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -real2hc 1 */

/*
 * This function contains 0 FP additions, 0 FP multiplications,
 * (or, 0 additions, 0 multiplications, 0 fused multiply/add),
 * 1 stack variables, and 2 memory accesses
 */

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_real2hc_1(const fftw_real *input, fftw_real *real_output, fftw_real *imag_output, int istride, int real_ostride, int imag_ostride)
{
     fftw_real tmp1;
     ASSERT_ALIGNED_DOUBLE;
     tmp1 = input[0];
     real_output[0] = tmp1;
}

fftw_codelet_desc fftw_real2hc_1_desc =
{
     "fftw_real2hc_1",
     (void (*)()) fftw_real2hc_1,
     1,
     FFTW_FORWARD,
     FFTW_REAL2HC,
     24,
     0,
     (const int *) 0,
};

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


text:
#ifndef gigp_header

#define gigp_header

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define GIGP_GAMMA_ROW 0
#define GIGP_B_ROW 1
#define GIGP_C_ROW 2
#define GIGP_N_ROW 3
#define GIGP_s_ROW 4
#define GIGP_S_ROW 5
#define GIGP_MAXCOUNT_ROW 6
#define GIGP_FIT_ROW 7

struct gigp_param{
        double gamma;
        double b;
        double c;
        double N;
        double s;
        double S;
        double max_count;
        double fit;
};

extern struct gigp_param* fit_exact_gigp(struct gigp_param* gigp_param,double* fx_unique_transcript_count, int len);
extern double gigp_dist (double x, void * p);
extern double sichel_function(const gsl_vector *v, void *params);
extern struct gigp_param* init_gigp_param(void);
extern int copy_gigp_param(struct gigp_param* source, struct gigp_param* target);
extern double give_me_sichel_p0(struct gigp_param* gigp_param ,double* p);

/* housekeeping stuff */
extern struct gigp_param* init_gigp_param(void);
extern int clear_gigp_param(struct gigp_param* tmp);
extern int copy_gigp_param(struct gigp_param* source, struct gigp_param* target);

#endif

max_stars_repo_name: TimoLassmann/SingleCellSim

max_issues_repo_name: TimoLassmann/SingleCellSim

max_forks_repo_name: TimoLassmann/SingleCellSim


text:
#ifndef MATH_FUNCTIONS_H
#define MATH_FUNCTIONS_H

#include "string.h"
#include <math.h>

extern "C" {
#include <cblas.h>
}

void caffe_copy(const int N, const float *X, float *Y);

void caffe_set(const int N, const float alpha, float *X);

void caffe_cpu_gemm(const CBLAS_TRANSPOSE TransA,
	const CBLAS_TRANSPOSE TransB, const int M, const int N, const int K,
	const float alpha, const float* A, const float* B, const float beta,
	float* C);

void caffe_sqr(const int N, const float* a, float* y);

float caffe_cpu_asum(const int n, const float* x);

void caffe_cpu_scale(const int n, const float alpha, const float *x, float* y);

void caffe_cpu_gemv(const CBLAS_TRANSPOSE TransA, const int M,
	const int N, const float alpha, const float* A, const float* x,
	const float beta, float* y);

void caffe_powx(const int n, const float* a, const float b, float* y);

void caffe_div(const int n, const float* a, const float* b, float* y);

void caffe_scal(const int N, const float alpha, float *X);

void caffe_mul(const int n, const float* a, const float* b, float* y);

void caffe_exp(const int n, const float* a, float* y);

#endif


text:
/* min/fsolver.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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>

#include "min.h"

static int
compute_f_values (gsl_function * f, double x_minimum, double * f_minimum,
                  double x_lower, double * f_lower,
                  double x_upper, double * f_upper);

static int
compute_f_values (gsl_function * f, double x_minimum, double * f_minimum,
                  double x_lower, double * f_lower,
                  double x_upper, double * f_upper)
{
  SAFE_FUNC_CALL(f, x_lower, f_lower);
  SAFE_FUNC_CALL(f, x_upper, f_upper);
  SAFE_FUNC_CALL(f, x_minimum, f_minimum);

  return GSL_SUCCESS;
}

int
gsl_min_fminimizer_set (gsl_min_fminimizer * s,
                        gsl_function * f,
                        double x_minimum, double x_lower, double x_upper)
{
  int status ;

  double f_minimum, f_lower, f_upper;

  status = compute_f_values (f, x_minimum, &f_minimum,
                             x_lower, &f_lower,
                             x_upper, &f_upper);

  if (status != GSL_SUCCESS)
    {
      return status ;
    }

  status = gsl_min_fminimizer_set_with_values (s, f, x_minimum, f_minimum,
                                               x_lower, f_lower,
                                               x_upper, f_upper);
  return status;
}

gsl_min_fminimizer *
gsl_min_fminimizer_alloc (const gsl_min_fminimizer_type * T)
{
  gsl_min_fminimizer * s =
    (gsl_min_fminimizer *) malloc (sizeof (gsl_min_fminimizer));

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

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

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

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

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

  return s;
}

int
gsl_min_fminimizer_set_with_values (gsl_min_fminimizer * s, gsl_function * f,
                                    double x_minimum, double f_minimum,
                                    double x_lower, double f_lower,
                                    double x_upper, double f_upper)
{
  s->function = f;
  s->x_minimum = x_minimum;
  s->x_lower = x_lower;
  s->x_upper = x_upper;

  if (x_lower > x_upper)
    {
      GSL_ERROR ("invalid interval (lower > upper)", GSL_EINVAL);
    }

  if (x_minimum >= x_upper || x_minimum <= x_lower)
    {
      GSL_ERROR ("x_minimum must lie inside interval (lower < x < upper)",
                 GSL_EINVAL);
    }

  s->f_lower = f_lower;
  s->f_upper = f_upper;
  s->f_minimum = f_minimum;

  if (f_minimum >= f_lower || f_minimum >= f_upper)
    {
      GSL_ERROR ("endpoints do not enclose a minimum", GSL_EINVAL);
    }

  return (s->type->set) (s->state, s->function,
                         x_minimum, f_minimum,
                         x_lower, f_lower,
                         x_upper, f_upper);
}

int
gsl_min_fminimizer_iterate (gsl_min_fminimizer * s)
{
  return (s->type->iterate) (s->state, s->function,
                             &(s->x_minimum), &(s->f_minimum),
                             &(s->x_lower), &(s->f_lower),
                             &(s->x_upper), &(s->f_upper));
}

void
gsl_min_fminimizer_free (gsl_min_fminimizer * s)
{
  free (s->state);
  free (s);
}

const char *
gsl_min_fminimizer_name (const gsl_min_fminimizer * s)
{
  return s->type->name;
}

/* Deprecated, use x_minimum instead */
double
gsl_min_fminimizer_minimum (const gsl_min_fminimizer * s)
{
  return s->x_minimum;
}

double
gsl_min_fminimizer_x_minimum (const gsl_min_fminimizer * s)
{
  return s->x_minimum;
}

double
gsl_min_fminimizer_x_lower (const gsl_min_fminimizer * s)
{
  return s->x_lower;
}

double
gsl_min_fminimizer_x_upper (const gsl_min_fminimizer * s)
{
  return s->x_upper;
}

double
gsl_min_fminimizer_f_minimum (const gsl_min_fminimizer * s)
{
  return s->f_minimum;
}

double
gsl_min_fminimizer_f_lower (const gsl_min_fminimizer * s)
{
  return s->f_lower;
}

double
gsl_min_fminimizer_f_upper (const gsl_min_fminimizer * s)
{
  return s->f_upper;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <mpi.h>

/* Pull in the CBLAS header. */
#ifdef _MACOSX
#include <Accelerate/Accelerate.h>
#else
#ifdef _ATLAS
#include <cblas.h>
#else
#include <gsl_cblas.h>
#endif /* _ATLAS */
#endif /* _MACOSX */

#include "ScaleME.h"

#include "precision.h"

#include "io.h"
#include "mlfma.h"
#include "direct.h"
#include "util.h"
#include "fsgreen.h"
#include "integrate.h"

fmadesc fmaconf;

static int farmatrow (cplx *, real, real *, real, int);
static int farmatcol (cplx *, real, real *, int, int);
static int acabuild (cplx **, real, real, int, int, real, int);
static int fullbuild (cplx **, real, int, int, real, int);
static int recaca (cplx *, cplx *, int, int, int, real);

/* Computes the far-field pattern for the specified group with the specified
 * center, and stores the output in a provided vector. sgn is positive for
 * radiation pattern and negative for receiving pattern. The list of "basis
 * functions" and center are ignored, since the calling program ensures each
 * FMM basis function is actually a single group of gridded elements with the
 * same affine grid. The center of the finest-level FMM group coincides with
 * the single FMM basis function contained therein. */
void farpattern (int nbs, int *bsl, void *vcrt, void *vpat, real *cen, int sgn) {
	cplx fact, beta = 1.0, *crt = (cplx *)vcrt,
		*pat = *((cplx **)vpat);

	if (sgn >= 0) {
		/* Scalar factors for the matrix multiplication. */
		fact = fmaconf.k0;

		/* Don't add in the existing output vector since it should
		 * have been zeroed anyway. */
		beta = 0.0;

		/* Perform the matrix-vector product. */
		GEMV (CblasColMajor, CblasNoTrans, fmaconf.nsamp,
				fmaconf.bspboxvol, &fact, fmaconf.radpats,
				fmaconf.nsamp, crt, 1, &beta, pat, 1);
	} else {
		/* Distribute the far-field patterns to the basis functions. */
		/* Scalar factors for the matrix multiplication. */
		fact = I * fmaconf.k0 * fmaconf.k0 / (4 * M_PI);

		/* Perform the matrix-vector product. */
		GEMV (CblasColMajor, CblasConjTrans, fmaconf.nsamp,
				fmaconf.bspboxvol, &fact, fmaconf.radpats,
				fmaconf.nsamp, pat, 1, &beta, crt, 1);
	}
}

/* Computes the far-field pattern for the specified group with the specified
 * center, and stores the output in a provided vector. sgn is positive for
 * radiation pattern and negative for receiving pattern. The list of "basis
 * functions" and center are ignored, since the calling program ensures each
 * FMM basis function is actually a single group of gridded elements with the
 * same affine grid. The center of the finest-level FMM group coincides with
 * the single FMM basis function contained therein. ACA is used to approximate
 * the matrix for efficient computations. */
void acafarpattern (int nbs, int *bsl, void *vcrt, void *vpat, real *cen, int sgn) {
	cplx fact, beta = 1.0, *crt = (cplx *)vcrt,
		*pat = *((cplx **)vpat), *work, *u, *v;

	u = fmaconf.radpats;
	v = fmaconf.radpats + fmaconf.acarank * fmaconf.nsamp;

	work = malloc (fmaconf.acarank * sizeof(cplx));

	if (sgn >= 0) {
		beta = 0.0;
		fact = 1.0;

		GEMV (CblasColMajor, CblasConjTrans, fmaconf.bspboxvol,
				fmaconf.acarank, &fact, v, fmaconf.bspboxvol,
				crt, 1, &beta, work, 1);

		/* Scalar factors for the matrix multiplication. */
		fact = fmaconf.k0;

		/* Perform the matrix-vector product. */
		GEMV (CblasColMajor, CblasNoTrans, fmaconf.nsamp,
				fmaconf.acarank, &fact, u, fmaconf.nsamp,
				work, 1, &beta, pat, 1);
	} else {
		beta = 0.0;
		fact = 1.0;

		GEMV (CblasColMajor, CblasConjTrans, fmaconf.nsamp,
				fmaconf.acarank, &fact, u, fmaconf.nsamp,
				pat, 1, &beta, work, 1);

		/* Distribute the far-field patterns to the basis functions. */
		/* Scalar factors for the matrix multiplication. */
		fact = I * fmaconf.k0 * fmaconf.k0 / (4 * M_PI);
		beta = 1.0;

		/* Perform the matrix-vector product. */
		GEMV (CblasColMajor, CblasNoTrans, fmaconf.bspboxvol,
				fmaconf.acarank, &fact, v, fmaconf.bspboxvol,
				work, 1, &beta, crt, 1);
	}

	free (work);
}

/* Build a column of the far-field signature for a point a distance rmc
 * from the center of the parent box. */
static int farmatcol (cplx *col, real k, real *rmc, int ntheta, int nphi) {
	int nsamp = nphi * (ntheta - 2) + 2;
	real cell[3] = { fmaconf.cell, fmaconf.cell, fmaconf.cell };

#pragma omp parallel default(shared)
{
	real s[3];
	int i;
#pragma omp for
	for (i = 0; i < nsamp; ++i) {
		/* Compute the Cartesian coordinates of the far-field sample. */
		sampcoords (s, i, ntheta, nphi);
		/* Compute the far-field sample. */
		col[i] = srcint(k, rmc, s, cell, fsplane);
	}
}

	return 0;
}

/* Build a row of the far-field signature for a fixed field sample s. */
static int farmatrow (cplx *row, real k, real *s, real dx, int bpd) {
	int bptot = bpd * bpd * bpd;
	real cell[3] = { fmaconf.cell, fmaconf.cell, fmaconf.cell };

#pragma omp parallel default(shared)
{
	int l;
	real dist[3];

#pragma omp for
	for (l = 0; l < bptot; ++l) {
		/* The distance from the basis to the box center. */
		cellcoords (dist, l, bpd, dx);

		/* The value of the far-field sample of the cell. */
		row[l] = srcint(k, dist, s, cell, fsplane);
	}
}
	return 0;
}

/* Construct an ACA approximation to the far-field signature matrix. */
static int acabuild (cplx **mats, real k0, real tol,
		int ntheta, int nphi, real dx, int bpd) {
	int maxrank, *irow, *icol, i, j, k,
	    nsamp = (ntheta - 2) * nphi + 2, nelt = bpd * bpd * bpd;
	cplx *u, *v, *uptr, *vptr, *row, *col, dpr, dpc;
	real dist[3], s[3], err = 0, tolsq;

	/* Compute the square of the tolerance and decrease by two orders
	 * of magnitude for recompression. */
	tolsq = tol * tol * 1e-2;

	maxrank = MIN(nelt, nsamp);

	/* Allocate the row and column index arrays. The extra element should
	 * only be stored (but never recalled) in the limiting case when the
	 * rank cannot be reduced. */
	irow = malloc((2 * maxrank + 1) * sizeof(int));
	icol = irow + maxrank + 1;

	/* Allocate the workspace for the row and column matrices. */
	u = malloc(maxrank * (nelt + nsamp) * sizeof(cplx));
	v = u + nsamp * maxrank;

	/* Start with the first row of the matrix. */
	irow[0] = 0;

	for (i = 0, row = v, col = u; i < maxrank; ++i, row += nelt, col += nsamp) {
		/* Find the coordinate of the observer element. */
		sampcoords (s, irow[i], ntheta, nphi);

		/* Fill the row for the selected observer element. */
		farmatrow (row, k0, s, dx, bpd);

		/* Subtract existing contributions from earlier ranks. */
#pragma omp parallel for default(shared) private(vptr, uptr, j, k)
		for (j = 0; j < nelt; ++j) {
			for (k = 0; k < i; ++k, uptr += nsamp, vptr += nelt) {
				uptr = u + k * nsamp;
				vptr = v + k * nelt;
				row[j] -= vptr[j] * uptr[irow[i]];
			}
		}

		/* Find the index and value of the maximum value. */
		icol[i] = maxind (row, nelt, icol, i);
		dpr = row[icol[i]];

		/* Scale the row. */
#pragma omp parallel for default(shared) private(j)
		for (j = 0; j < nelt; ++j) row[j] /= dpr;

		/* Find the coordinates of the source cell. */
		cellcoords (dist, icol[i], bpd, dx);

		/* Construct the radiation pattern of the source cell. */
		farmatcol (col, k0, dist, ntheta, nphi);

		/* Subtract existing contributions from earlier ranks. */
#pragma omp parallel for default(shared) private(vptr, uptr, j, k)
		for (j = 0; j < nsamp; ++j) {
			for (k = 0; k < i; ++k) {
				vptr = v + k * nelt;
				uptr = u + k * nsamp;
				col[j] -= vptr[icol[i]] * uptr[j];
			}
		}

		/* Update the error approximation. */
		for (k = 0, uptr = u, vptr = v; k < i; ++k) {
			DOTU_SUB (nelt, vptr, 1, row, 1, &dpr);
			DOTU_SUB (nsamp, uptr, 1, col, 1, &dpc);
			uptr += nsamp;
			vptr += nelt;

			err += 2.0 * cabs(dpr) * cabs(dpc);
		}

		DOTC_SUB (nelt, row, 1, row, 1, &dpr);
		DOTC_SUB (nsamp, col, 1, col, 1, &dpc);

		err += creal(dpr) * creal(dpc);

		if (creal(dpr) * creal(dpc) <= tolsq * err) break;

		/* Find the next row index. */
		irow[i + 1] = maxind (col, nsamp, irow, i + 1);
	}

	maxrank = i;

	/* Conjugate the matrix v. */
	for (i = 0, k = maxrank * nelt; i < k; ++i) v[i] = conj(v[i]);

	/* Recompress (in place) the matrices with the actual tolerance. */
	maxrank = recaca (u, v, nsamp, nelt, maxrank, tol);

	/* Allocate the final matrix storage. */
	*mats = malloc (maxrank * (nelt + nsamp) * sizeof(cplx));
	/* Copy the colum matrix in first, then the row matrix. */
	memcpy (*mats, u, maxrank * nsamp * sizeof(cplx));
	memcpy (*mats + maxrank * nsamp, v, maxrank * nelt * sizeof(cplx));

	/* Free the work arrays. */
	free (irow);
	free (u);

	return maxrank;
}

/* Recompress an ACA approximation using truncated singular values. */
static int recaca (cplx *u, cplx *v, int m, int n, int k, real tol) {
	int rank, lwork, info, i, j, l;
	real *rwork, *ss;
	cplx *qu, *qv, *rp, *work, *tu, *tv, *us, *vs, alpha = 1.0, beta = 0.0;

	/* Allocate space for the orthogonal matrices and the outer product. */
	qu = malloc (k * (m + n + 3 * k) * sizeof(cplx));
	qv = qu + k * m;
	rp = qv + k * n;
	us = rp + k * k;
	vs = us + k * k;

	/* Allocate space for the reflector coefficients. */
	tu = malloc (2 * k * sizeof(cplx));
	tv = tu + k;

	/* Allocate space for the singular values and real work. */
	ss = malloc (6 * k * sizeof(real));
	rwork = ss + k;

	/* Figure out the maximum size of the work array. */
	lwork = -1;
	GEQRF (&m, &k, qu, &m, tu, qu, &lwork, &info);
	GEQRF (&n, &k, qv, &n, tv, qu + 1, &lwork, &info);
	UNGQR (&m, &k, &k, qu, &m, tu, qu + 2, &lwork, &info);
	UNGQR (&n, &k, &k, qv, &n, tv, qu + 3, &lwork, &info);
	GESVD ("S", "S", &k, &k, rp, &k, ss, us, &k, vs, &k, qu + 4, &lwork, rwork, &info);

	/* Find the maximum workspace size. */
	for (i = 0; i < 5; ++i) lwork = MAX(lwork, (int)creal(qu[i]));

	work = malloc(lwork * sizeof(cplx));

	/* Copy the input matrices into the workspaces. */
	memcpy (qu, u, (m * k) * sizeof(cplx));
	memcpy (qv, v, (n * k) * sizeof(cplx));

	/* Compute the QR factorizations of the row and column matrices. */
	GEQRF (&m, &k, qu, &m, tu, work, &lwork, &info);
	GEQRF (&n, &k, qv, &n, tv, work, &lwork, &info);

	/* Compute the outer product of the triangular matrices. */
	for (i = 0; i < k; ++i) {
		for (j = 0; j < k; ++j) {
			rp[i + j * k] = 0.0;
			for (l = MAX(i,j); l < k; ++l)
				rp[i + j * k] += qu[i + l * m] * conj(qv[j + l * n]);
		}
	}

	/* Compute the SVD of the triangular matrix product. */
	GESVD ("S", "S", &k, &k, rp, &k, ss, us, &k, vs, &k, work, &lwork, rwork, &info);

	/* Find the truncated rank. */
	for (rank = 0; rank < k; ++rank)
		if (fabs(ss[rank] / ss[0]) < tol) break;

	/* Multiply the singular values by the column matrix. */
	for (i = 0; i < rank; ++i)
		for (j = 0; j < k; ++j)
			vs[i + k * j] *= ss[i];

	/* Expand the orthogonal matrices. */
	UNGQR (&m, &k, &k, qu, &m, tu, work, &lwork, &info);
	UNGQR (&n, &k, &k, qv, &n, tv, work, &lwork, &info);

	/* Compute the new row and column matrices. */
	GEMM (CblasColMajor, CblasNoTrans, CblasNoTrans, m, rank, k,
			&alpha, qu, m, us, k, &beta, u, m);
	GEMM (CblasColMajor, CblasNoTrans, CblasConjTrans, n, rank, k,
			&alpha, qv, n, vs, k, &beta, v, n);

	free (qu);
	free (tu);
	free (ss);
	free (work);

	return rank;
}

static int svdbuild (cplx **mats, real k0, real tol,
		int ntheta, int nphi, real dx, int bpd) {
	int nsamp = (ntheta - 2) * nphi + 2, nelt = bpd * bpd * bpd,
	    lwork, info, mindim, rank, l, i;
	cplx *col, *a, *u, *vt, *work, *ptr;
	real dist[3], *s, *rwork;

	/* The minimum dimension. */
	mindim = MIN(nsamp, nelt);

	/* Allocate the full matrix. */
	a = malloc(nsamp * nelt * sizeof(cplx));

	/* Allocate the singular vector matrices. */
	u = malloc(mindim * (nsamp + nelt) * sizeof(cplx));
	vt = u + mindim * nsamp;

	/* Allocate the real workspace and the singular value array. */
	s = malloc(6 * mindim * sizeof(real));
	rwork = s + mindim;

	/* Perform a workspace query. */
	lwork = -1;
	GESVD ("S", "S", &nsamp, &nelt, a, &nsamp, s, u, &nsamp, vt, &mindim,
			a, &lwork, rwork, &info);

	/* Allocate the workspace. */
	lwork = creal(a[0]);
	work = malloc(lwork * sizeof(cplx));

	/* Loop through all columns (source grid elements) and build the matrix. */
	for (l = 0, col = a; l < nelt; ++l, col += nsamp) {
		/* The relative position of the source grid element. */
		cellcoords (dist, l, bpd, dx);

		/* Build the corresponding matrix column. */
		farmatcol (col, k0, dist, ntheta, nphi);
	}

	/* Perform an SVD on the matrix. */
	GESVD ("S", "S", &nsamp, &nelt, a, &nsamp, s, u, &nsamp, vt, &mindim,
			work, &lwork, rwork, &info);

	/* Find the first rank below the desired tolerance. */
	for (rank = 0; rank < mindim; ++rank)
		if (fabs(s[rank] / s[0]) < tol) break;

	*mats = malloc (rank * (nelt + nsamp) * sizeof(cplx));

	/* Multiply the left singular vector matrix by the singular values. */
	for (l = 0, col = u, ptr = *mats; l < rank; ++l)
		for (i = 0; i < nsamp; ++i, ++col, ++ptr)
			*ptr = s[l] * (*col);

	/* Tranpose and conjugate the right singular vector matrix. */
	for (l = 0, ptr = *mats + rank * nsamp; l < rank; ++l)
		for (i = 0; i < nelt; ++i)
			ptr[i + l * nelt] = conj(vt[l + i * mindim]);

	free (a);
	free (u);
	free (s);
	free (work);

	return rank;
}

static int fullbuild (cplx **mats, real k0, int ntheta, int nphi, real dx, int bpd) {
	int nsamp = (ntheta - 2) * nphi + 2, nelt = bpd * bpd * bpd, l;
	cplx *col;
	real dist[3];

	/* Allocate the full far-field matrix. */
	*mats = malloc(nsamp * nelt * sizeof(cplx));

	/* Loop through all columns (source grid elements). */
	for (l = 0, col = *mats; l < nelt; ++l, col += nsamp) {
		/* The relative position of the source grid element. */
		cellcoords (dist, l, bpd, dx);

		/* Build the corresponding matrix column. */
		farmatcol (col, k0, dist, ntheta, nphi);
	}

	return nelt * nsamp;
}

/* Precomputes the near interactions for redundant calculations and sets up
 * the wave vector directions to be used for fast calculation of far-field patterns. */
int fmmprecalc (real acatol, int useaca) {
	int ntheta, nphi, rank, i;

	MPI_Comm_rank (MPI_COMM_WORLD, &rank);

	/* Get the finest level parameters. */
	ScaleME_getFinestLevelParams (&(fmaconf.nsamp), &ntheta, &nphi, NULL);

	if (!useaca) {
		/* Build the direct far-field matrices. */
		fmaconf.acarank = 0;
		i = fullbuild (&(fmaconf.radpats), fmaconf.k0, ntheta,
				nphi, fmaconf.cell, fmaconf.bspbox);
		fprintf (stderr, "Rank %d: Far-field matrix element count: %d\n", rank, i);
	} else {
		/* Use ACA if the tolerance is positive, otherwise use SVD. */
		if (acatol > 0)
			fmaconf.acarank = acabuild (&(fmaconf.radpats),
					fmaconf.k0, acatol, ntheta, nphi,
					fmaconf.cell, fmaconf.bspbox);
		else
			fmaconf.acarank = svdbuild (&(fmaconf.radpats),
					fmaconf.k0, -acatol, ntheta,
					nphi, fmaconf.cell, fmaconf.bspbox);

		fprintf (stderr, "Rank %d: Far-field matrix rank: %d\n", rank, fmaconf.acarank);
	}

	return fmaconf.acarank;
}

/* initialisation and finalisation routines for ScaleME */
int ScaleME_preconf (int useaca) {
	int error;
	real len, cen[3];

	/* The problem and tree are both three-dimensional. */
	ScaleME_setDimen (3);
	ScaleME_setTreeType (3);

	/* The fields are scalar-valued. */
	ScaleME_setFields (1);

	/* The wave number is real-valued. */
	ScaleME_setWaveNumber (fmaconf.k0);

	/* Set some MLFMA parameters. */
	ScaleME_setNumBasis (fmaconf.gnumbases);
	ScaleME_setMaxLevel (fmaconf.maxlev);
	ScaleME_setPrecision (fmaconf.precision);
	ScaleME_setMAC (fmaconf.numbuffer);
	ScaleME_setInterpOrder (fmaconf.interpord);

	ScaleME_setTopComputeLevel (fmaconf.toplev);

	ScaleME_setRHSDataWidth (fmaconf.bspboxvol);

	/* Enable fast O2I if desired. */
	if (fmaconf.fo2itxlev > 0 && fmaconf.fo2ibclev > 0) {
		ScaleME_useFastO2IGhosts (1);
		ScaleME_selectFastO2I (fmaconf.fo2itxlev, fmaconf.fo2ibclev,
				fmaconf.fo2iord, fmaconf.fo2iosr);
	}

	/* Set the root box length. */
	len = (1 <<  fmaconf.maxlev) * fmaconf.bspbox * fmaconf.cell;
	/* Position the root box properly. */
	cen[0] = fmaconf.min[0] + 0.5 * len;
	cen[1] = fmaconf.min[1] + 0.5 * len;
	cen[2] = fmaconf.min[2] + 0.5 * len;
	ScaleME_setRootBox (len, cen);

	/*  let all processes start the initialisation together */
	MPI_Barrier(MPI_COMM_WORLD);

	/* open the std files */
	MPFMA_stdout = stdout;
	MPFMA_stderr = stderr;

	/* Use the external near-field interactions. */
	ScaleME_setBlockDirInterFunc (blockinteract);
	if (useaca) ScaleME_useExternFarField (acafarpattern);
	else ScaleME_useExternFarField (farpattern);

	/* Finish the setup with the external interactions. */
	error = ScaleME_initSetUp (MPI_COMM_WORLD, NULL, NULL, NULL, bscenter);

	if (error) {
		fprintf(stdout, "ERROR: ScaleME pre-init failed.\n");
		ScaleME_finalizeParHostFMA();
		return -1;
	}

	return 0;
}

int ScaleME_postconf (void) {
	if (ScaleME_completeSetUp()) {
		fprintf(stdout, "ERROR: ScaleME setup routine failed.\n");
		goto error_handle;
	}

	if (ScaleME_initParHostDataStructs()) {
		fprintf(stdout, "ERROR: ScaleME parallel host init failed.\n");
		goto error_handle;
	}

	return 0;

error_handle:
	ScaleME_finalizeParHostFMA();
	return -1;
}

max_stars_repo_name: ahesford/afma

max_issues_repo_name: ahesford/afma

max_forks_repo_name: ahesford/afma


text:
/* rng/vax.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, 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_rng.h>

/* This is the old vax generator MTH$RANDOM. The sequence is,

   x_{n+1} = (a x_n + c) mod m

   with a = 69069, c = 1 and m = 2^32. The seed specifies the initial
   value, x_1.

   The theoretical value of x_{10001} is 3051034865.

   The period of this generator is 2^32. */

static inline unsigned long int vax_get (void *vstate);
static double vax_get_double (void *vstate);
static void vax_set (void *state, unsigned long int s);

typedef struct
  {
    unsigned long int x;
  }
vax_state_t;

static inline unsigned long int
vax_get (void *vstate)
{
  vax_state_t *state = (vax_state_t *) vstate;

  state->x = (69069 * state->x + 1) & 0xffffffffUL;

  return state->x;
}

static double
vax_get_double (void *vstate)
{
  return vax_get (vstate) / 4294967296.0 ;
}

static void
vax_set (void *vstate, unsigned long int s)
{
  vax_state_t *state = (vax_state_t *) vstate;

  /* default seed is 0. The constant term c stops the series from
     collapsing to 0,0,0,0,0,... */

  state->x = s;

  return;
}

static const gsl_rng_type vax_type =
{"vax",				/* name */
 0xffffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (vax_state_t),
 &vax_set,
 &vax_get,
 &vax_get_double};

const gsl_rng_type *gsl_rng_vax = &vax_type;


text:
/*
 * Matrix.h
 *
 *  Created on: Dec 24, 2014
 *      Author: Phong
 */

#ifndef MATRIX_H_
#define MATRIX_H_

#include "Default.h"
#include <cblas.h>
#include <fstream>
#include <random>

using namespace std;

// row-major

class Matrix {
public:
	int rows, cols, length;
	real* data;

    ~Matrix();

	Matrix(int rows, int cols);

	// vector
	Matrix(int rows);

	Matrix(int rows, int cols, real* data);

	Matrix(int rows, real* data);

	real get(int r, int c);

	void put(int r, int c, real value);

	Matrix* dup();

	void copy(Matrix* other);

	void errorIfNotSameSize(Matrix* other);

	void fill(real value);

    void save(ofstream& f);

    static Matrix* load(ifstream& f);

    void print(bool trans = false);

	// y <- alpha * x + y
	Matrix* addi(real alpha, Matrix* x);

	// y <- x + y
	Matrix* addi(Matrix* x);

	// y[i] <- y[i] + a
	Matrix* addi(real a);

	// z <- alpha * x + y
	Matrix* add(real alpha, Matrix* x);

	// y <- x + y
	Matrix* add(Matrix* x);

	// z[i] <- y[i] + a
	Matrix* add(real a);

    // z[i] <- x[i] * y[i]
    Matrix* mul(Matrix* x);

	// y[i] <- x[i] * y[i]
	Matrix* muli(Matrix* x);

	// y[i] <- a * y[i]
	Matrix* muli(real a);

	// z[i] <- a * y[i]
	Matrix* mul(real a);

	Matrix* divi(real a);

    Matrix* divi(Matrix* x);

	Matrix* div(real a);

    Matrix* sqrt();

    Matrix* sqrti();

	Matrix* getColumn(int col);

	Matrix* getRows(int* ids, int n);

	Matrix* addColumn(int c, Matrix* x);

	real sum();

	static void free(Matrix** &As, int n);

	static void free(Matrix*** &As, int n1, int n2);

	static Matrix* zeros(int rows, int cols);

	static Matrix* zeros(int dim);

	static Matrix* eye(int dim);

	static Matrix* uniform(int rows, int cols, real a, real b);

	static Matrix* uniform(int rows, real a, real b);

	static Matrix* bernoulli(int rows, int cols, real p);

	static Matrix* bernoulli(int rows, real p);

	static Matrix* normal(int rows, int cols, real mean, real dev);

	static Matrix* normal(int rows, real mean, real dev);

	static void gemv(real alpha, enum CBLAS_TRANSPOSE transA, Matrix* A, Matrix* x, real beta, Matrix* y);

	static real dot(Matrix* X, Matrix* Y);

	static void ger(real alpha, Matrix* x, Matrix* y, Matrix* A);

	static void axpy(real alpha, Matrix* x, Matrix* y);

	static real nrm2(Matrix* x);

};

#endif /* MATRIX_H_ */

max_stars_repo_name: zhangmeishan/lstm-rnn

max_issues_repo_name: zhangmeishan/lstm-rnn

max_forks_repo_name: zhangmeishan/lstm-rnn