text:
/**
 *
 * @file core_dparfb.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Dulceneia Becker
 * @date 2011-06-14
 * @generated d Tue Jan  7 11:44:48 2014
 *
 **/
#ifdef USE_MKL
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#else
#include <cblas.h>
#include <lapacke.h>
#endif
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_double
 *
 *  CORE_dparfb applies a complex upper triangular block reflector H
 *  or its transpose H' to a complex rectangular matrix formed by
 *  coupling two tiles A1 and A2. Matrix V is:
 *
 *          COLUMNWISE                    ROWWISE
 *
 *         |     K     |                 |      N2-L     |   L  |
 *      __ _____________ __           __ _________________        __
 *         |    |      |                 |               | \
 *         |    |      |                 |               |   \    L
 *    M2-L |    |      |              K  |_______________|_____\  __
 *         |    |      | M2              |                      |
 *      __ |____|      |                 |                      | K-L
 *         \    |      |              __ |______________________| __
 *       L   \  |      |
 *      __     \|______| __              |          N2          |
 *
 *         | L |  K-L  |
 *
 *******************************************************************************
 *
 * @param[in] side
 *         @arg PlasmaLeft  : apply Q or Q**T from the Left;
 *         @arg PlasmaRight : apply Q or Q**T from the Right.
 *
 * @param[in] trans
 *         @arg PlasmaNoTrans   : No transpose, apply Q;
 *         @arg PlasmaTrans : ConjTranspose, apply Q**T.
 *
 * @param[in] direct
 *         Indicates how H is formed from a product of elementary
 *         reflectors
 *         @arg PlasmaForward  : H = H(1) H(2) . . . H(k) (Forward)
 *         @arg PlasmaBackward : H = H(k) . . . H(2) H(1) (Backward)
 *
 * @param[in] storev
 *         Indicates how the vectors which define the elementary
 *         reflectors are stored:
 *         @arg PlasmaColumnwise
 *         @arg PlasmaRowwise
 *
 * @param[in] M1
 *         The number of columns of the tile A1. M1 >= 0.
 *
 * @param[in] N1
 *         The number of rows of the tile A1. N1 >= 0.
 *
 * @param[in] M2
 *         The number of columns of the tile A2. M2 >= 0.
 *
 * @param[in] N2
 *         The number of rows of the tile A2. N2 >= 0.
 *
 * @param[in] K
 *         The order of the matrix T (= the number of elementary
 *         reflectors whose product defines the block reflector).
 *
 * @param[in] L
 *         The size of the triangular part of V
 *
 * @param[in,out] A1
 *         On entry, the M1-by-N1 tile A1.
 *         On exit, A1 is overwritten by the application of Q.
 *
 * @param[in] LDA1
 *         The leading dimension of the array A1. LDA1 >= max(1,N1).
 *
 * @param[in,out] A2
 *         On entry, the M2-by-N2 tile A2.
 *         On exit, A2 is overwritten by the application of Q.
 *
 * @param[in] LDA2
 *         The leading dimension of the tile A2. LDA2 >= max(1,N2).
 *
 * @param[in] V
 *         (LDV,K) if STOREV = 'C'
 *         (LDV,M2) if STOREV = 'R' and SIDE = 'L'
 *         (LDV,N2) if STOREV = 'R' and SIDE = 'R'
 *         Matrix V.
 *
 * @param[in] LDV
 *         The leading dimension of the array V.
 *         If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M2);
 *         if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N2);
 *         if STOREV = 'R', LDV >= K.
 *
 * @param[out] T
 *         The triangular K-by-K matrix T in the representation of the
 *         block reflector.
 *         T is upper triangular by block (economic storage);
 *         The rest of the array is not referenced.
 *
 * @param[in] LDT
 *         The leading dimension of the array T. LDT >= K.
 *
 * @param[in,out] WORK
 *
 * @param[in] LDWORK
 *         The dimension of the array WORK.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
/* This kernel is never traced so return type on previous line for convert2eztrace.pl script */
int
CORE_dparfb(PLASMA_enum side, PLASMA_enum trans, PLASMA_enum direct, PLASMA_enum storev,
            int M1, int N1, int M2, int N2, int K, int L,
                  double *A1, int LDA1,
                  double *A2, int LDA2,
            const double *V, int LDV,
            const double *T, int LDT,
                  double *WORK, int LDWORK)
{
    static double zone  =  1.0;
    static double mzone = -1.0;

    int j;

    /* Check input arguments */
    if ((side != PlasmaLeft) && (side != PlasmaRight)) {
        coreblas_error(1, "Illegal value of side");
        return -1;
    }
    if ((trans != PlasmaNoTrans) && (trans != PlasmaTrans)) {
        coreblas_error(2, "Illegal value of trans");
        return -2;
    }
    if ((direct != PlasmaForward) && (direct != PlasmaBackward)) {
        coreblas_error(3, "Illegal value of direct");
        return -3;
    }
    if ((storev != PlasmaColumnwise) && (storev != PlasmaRowwise)) {
        coreblas_error(4, "Illegal value of storev");
        return -4;
    }
    if (M1 < 0) {
        coreblas_error(5, "Illegal value of M1");
        return -5;
    }
    if (N1 < 0) {
        coreblas_error(6, "Illegal value of N1");
        return -6;
    }
    if ((M2 < 0) ||
        ( (side == PlasmaRight) && (M1 != M2) ) ) {
        coreblas_error(7, "Illegal value of M2");
        return -7;
    }
    if ((N2 < 0) ||
        ( (side == PlasmaLeft) && (N1 != N2) ) ) {
        coreblas_error(8, "Illegal value of N2");
        return -8;
    }
    if (K < 0) {
        coreblas_error(9, "Illegal value of K");
        return -9;
    }

    /* Quick return */
    if ((M1 == 0) || (N1 == 0) || (M2 == 0) || (N2 == 0) || (K == 0))
        return PLASMA_SUCCESS;

    if (direct == PlasmaForward) {

        if (side == PlasmaLeft) {

            /*
             * Column or Rowwise / Forward / Left
             * ----------------------------------
             *
             * Form  H * A  or  H' * A  where  A = ( A1 )
             *                                     ( A2 )
             */

            /* W = A1 + op(V) * A2 */
            CORE_dpamm(
                    PlasmaW, PlasmaLeft, storev,
                    K, N1, M2, L,
                    A1, LDA1,
                    A2, LDA2,
                    V, LDV,
                    WORK, LDWORK);

            /* W = op(T) * W */
            cblas_dtrmm(
                CblasColMajor, CblasLeft, CblasUpper,
                (CBLAS_TRANSPOSE)trans, CblasNonUnit, K, N2,
                (zone), T, LDT, WORK, LDWORK);

            /* A1 = A1 - W */
            for(j = 0; j < N1; j++) {
                cblas_daxpy(
                        K, (mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
            }

            /* A2 = A2 - op(V) * W  */
            /* W also changes: W = V * W, A2 = A2 - W */
            CORE_dpamm(
                    PlasmaA2, PlasmaLeft, storev,
                    M2, N2, K, L,
                    A1, LDA1,
                    A2, LDA2,
                    V, LDV,
                    WORK, LDWORK);
        }
        else {
            /*
             * Column or Rowwise / Forward / Right
             * -----------------------------------
             *
             * Form  H * A  or  H' * A  where A  = ( A1 A2 )
             *
             */

            /* W = A1 + A2 * op(V) */
            CORE_dpamm(
                    PlasmaW, PlasmaRight, storev,
                    M1, K, N2, L,
                    A1, LDA1,
                    A2, LDA2,
                    V, LDV,
                    WORK, LDWORK);

            /* W = W * op(T) */
            cblas_dtrmm(
                CblasColMajor, CblasRight, CblasUpper,
                (CBLAS_TRANSPOSE)trans, CblasNonUnit, M2, K,
                (zone), T, LDT, WORK, LDWORK);

            /* A1 = A1 - W */
            for(j = 0; j < K; j++) {
                cblas_daxpy(
                        M1, (mzone),
                        &WORK[LDWORK*j], 1,
                        &A1[LDA1*j], 1);
            }

            /* A2 = A2 - W * op(V) */
            /* W also changes: W = W * V', A2 = A2 - W */
            CORE_dpamm(
                    PlasmaA2, PlasmaRight, storev,
                    M2, N2, K, L,
                    A1, LDA1,
                    A2, LDA2,
                    V, LDV,
                    WORK, LDWORK);
        }
    }
    else {
        coreblas_error(3, "Not implemented (Backward / Left or Right)");
        return PLASMA_ERR_NOT_SUPPORTED;
    }

    return PLASMA_SUCCESS;
}


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

#include <gsl.h>
#include <nlohmann/json_fwd.hpp>
#include <cstdint>
#include <string>

namespace cb {
namespace crypto {
enum class Algorithm { MD5, SHA1, SHA256, SHA512 };

bool isSupported(Algorithm algorithm);

const int MD5_DIGEST_SIZE = 16;
const int SHA1_DIGEST_SIZE = 20;
const int SHA256_DIGEST_SIZE = 32;
const int SHA512_DIGEST_SIZE = 64;

/**
 * Generate a HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string HMAC(Algorithm algorithm,
                 std::string_view key,
                 std::string_view data);

/**
 * Generate a PBKDF2_HMAC digest of the key and data by using the given
 * algorithm
 *
 * @throws std::invalid_argument - unsupported algorithm
 *         std::runtime_error - Failures generating the HMAC
 */
std::string PBKDF2_HMAC(Algorithm algorithm,
                        const std::string& pass,
                        std::string_view salt,
                        unsigned int iterationCount);

/**
 * Generate a digest by using the requested algorithm
 */
std::string digest(Algorithm algorithm, std::string_view data);

enum class Cipher { AES_256_cbc };

Cipher to_cipher(const std::string& str);

/**
 * Encrypt the specified data by using a given cipher
 *
 * @param cipher The cipher to use
 * @param key The key used for encryption
 * @param ivec The IV to use for encryption
 * @param data The Pointer to the data to encrypt
 * @return The encrypted data
 */
std::string encrypt(Cipher cipher,
                    std::string_view key,
                    std::string_view iv,
                    std::string_view data);

/**
 * Encrypt the specified data
 *
 * {
 *    "cipher" : "name of the cipher",
 *    "key" : "base64 of the raw key",
 *    "iv" : "base64 of the raw iv"
 * }
 *
 * @param meta the json description of the encryption to use
 * @param data Pointer to the data to encrypt
 * @param length the length of the data to encrypt
 * @return The encrypted data
 */
std::string encrypt(const nlohmann::json& json, std::string_view data);

/**
 * Decrypt the specified data by using a given cipher
 *
 * @param cipher The cipher to use
 * @param key The key used for encryption
 * @param ivec The IV to use for encryption
 * @param data The data to decrypt
 * @return The decrypted data
 */
std::string decrypt(Cipher cipher,
                    std::string_view key,
                    std::string_view iv,
                    std::string_view data);

/**
 * Decrypt the specified data
 *
 * {
 *    "cipher" : "name of the cipher",
 *    "key" : "base64 of the raw key",
 *    "iv" : "base64 of the raw iv"
 * }
 *
 * @param meta the json description of the encryption to use
 * @param data Pointer to the data to decrypt
 * @param length the length of the data to decrypt
 * @return The decrypted data
 */
std::string decrypt(const nlohmann::json& json, std::string_view data);

} // namespace crypto
} // namespace cb

ext: h


text:
#ifndef PCH_H
#define PCH_H

#define ISOLATION_AWARE_ENABLED 1

/*
    Only put files that never change or are rarely changed in this file.
    An include graph will be provided for some headers that were not explicit
    added. VC++ is a requirement.

    We don't care about warnings for library files.
*/
#pragma warning(push)
#pragma warning(disable : 4091 4191 4365 4464 4571 4619 4623 4625 4626 4768 4774 5026 5027 5039 5045)
/* Windows Stuff */
#include "windows_nodefines.h"
#include <ShlObj.h> // CommCtrl
#include <WindowsX.h>
#include <atlbase.h> // atldef.h (Windows.h), atlcore.h (tchar.h), Shlwapi.h
#include <direct.h>
#include <io.h>
#include <wincodec.h>

#if WIN32_LEAN_AND_MEAN
#include <CommDlg.h>
#include <MMSystem.h>
#include <ShellAPI.h>
#endif // WIN32_LEAN_AND_MEAN

/*
    Nuke WinAPI #defines - anything including this file should be using the
    ra versions.
*/
#undef CreateDirectory
#undef GetMessage

/* C Stuff */
#include <cctype>

/* STL Stuff */
#include <array> // algorithm, iterator, tuple
#include <atomic>
#include <fstream>
#include <iomanip>
#include <map>
#include <memory>
#include <mutex> // chrono (time.h), functional, thread, utility
#include <queue> // deque, vector, algorithm
#include <set>
#include <sstream> // string
#include <stack>
#include <unordered_map>
#include <variant>

#pragma warning(push)
// Look at this file if you want to know what they are and what they mean
#include "pch_cppcorecheck_suppressions.h"
#include "pch_microsoft_suppressions.h"

/* RapidJSON Stuff */
#define RAPIDJSON_HAS_STDSTRING 1
#define RAPIDJSON_NOMEMBERITERATORCLASS 1
#include <rapidjson\document.h> // has reader.h
#include <rapidjson\error\en.h>
#include <rapidjson\istreamwrapper.h>
#include <rapidjson\ostreamwrapper.h>
#include <rapidjson\writer.h> // has stringbuffer.h

/* gsl stuff */
#define GSL_THROW_ON_CONTRACT_VIOLATION
#include <gsl\gsl>

// the default GSL Expects() macro throws an exception, which we don't try to catch anywhere.
// replace it with our custom handler, which will log and report the error before throwing it.
#ifndef RA_UTEST
 #ifdef NDEBUG
  extern void __gsl_contract_handler(const char* const file, unsigned int line);
  #undef Expects
  #define Expects(cond) (GSL_LIKELY(cond) ? static_cast<void>(0) : __gsl_contract_handler(__FILE__, __LINE__))
 #else
  extern void __gsl_contract_handler(const char* const file, unsigned int line, const char* const error);
  #undef Expects
  #define Expects(cond) (GSL_LIKELY(cond) ? static_cast<void>(0) : __gsl_contract_handler(__FILE__, __LINE__, GSL_STRINGIFY(cond)))
 #endif
#endif

#pragma warning(pop)

#if RA_UTEST
#include <CodeAnalysis\Warnings.h>
#pragma warning(push)
#pragma warning(disable : ALL_CPPCORECHECK_WARNINGS)
#include "CppUnitTest.h"
#pragma warning(pop)
#endif /* RA_UTEST */

/* rcheevos stuff */
#pragma warning(push)
#pragma warning(disable : 4201) // nameless struct
#include <rcheevos.h>
#pragma warning(pop)
#pragma warning(pop)

#endif /* !PCH_H */

max_stars_repo_name: Jamiras/RAIntegration

max_issues_repo_name: Jamiras/RAIntegration

max_forks_repo_name: Jamiras/RAIntegration


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

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

#include "ccl.h"

/*----- ROUTINE: dc_NakamuraSuto -----
INPUT: cosmology, scale factor
TASK: Computes the peak threshold: delta_c(z) assuming LCDM.
Cosmology dependence of the critical linear density according to the spherical-collapse model.
Fitting function from Nakamura & Suto (1997; arXiv:astro-ph/9710107).
*/
double dc_NakamuraSuto(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double dc0 = (3./20.)*pow(12.*M_PI,2./3.);
  double dc = dc0*(1.+0.012299*log10(Om_mz));

  return dc;

}

/*----- ROUTINE: Dv_BryanNorman -----
INPUT: cosmology, scale factor
TASK: Computes the virial collapse density contrast with respect to the matter density assuming LCDM.
Cosmology dependence of the virial collapse density according to the spherical-collapse model
Fitting function from Bryan & Norman (1998; arXiv:astro-ph/9710107)
*/
double Dv_BryanNorman(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double x = Om_mz-1.;
  double Dv0 = 18.*pow(M_PI,2);
  double Dv = (Dv0+82.*x-39.*pow(x,2))/Om_mz;

  return Dv;
}

static double sigmaM_m2r(ccl_cosmology *cosmo, double halomass, int *status)
{
  double rho_m, smooth_radius;

  // Comoving matter density
  rho_m = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

  smooth_radius = pow((3.0*halomass) / (4*M_PI*rho_m), (1.0/3.0));

  return smooth_radius;
}

void ccl_cosmology_compute_sigma(ccl_cosmology *cosmo, ccl_f2d_t *psp, int *status)
{
  if(cosmo->computed_sigma)
    return;

  int na = cosmo->spline_params.A_SPLINE_NA_SM + cosmo->spline_params.A_SPLINE_NLOG_SM - 1;
  int nm = cosmo->spline_params.LOGM_SPLINE_NM;
  double *m = NULL;
  double *y = NULL;
  double *aa = NULL;

  // create linearly-spaced values of log-mass.
  m = ccl_linear_spacing(cosmo->spline_params.LOGM_SPLINE_MIN,
                         cosmo->spline_params.LOGM_SPLINE_MAX, nm);
  if (m == NULL ||
      (fabs(m[0]-cosmo->spline_params.LOGM_SPLINE_MIN)>1e-5) ||
      (fabs(m[nm-1]-cosmo->spline_params.LOGM_SPLINE_MAX)>1e-5) ||
      (m[nm-1]>10E17)) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,
                                     "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
                                     "Error creating linear spacing in m\n");
  }

  // create scale factor array
  if (*status == 0) {
    aa = ccl_linlog_spacing(cosmo->spline_params.A_SPLINE_MINLOG_SM,
                            cosmo->spline_params.A_SPLINE_MIN_SM,
                            cosmo->spline_params.A_SPLINE_MAX,
                            cosmo->spline_params.A_SPLINE_NLOG_SM,
                            cosmo->spline_params.A_SPLINE_NA_SM);
    if (aa == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
                                       "Error creating scale factor array\n");
    }
  }

  // create space for y, to be filled with sigma
  if (*status == 0) {
    y = malloc(sizeof(double)*nm*na);
    if (y == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
                                       "memory allocation\n");
    }
  }

  // fill in sigma, if no errors have been triggered at this time.
  if (*status == 0) {
    #pragma omp parallel shared(na, aa, nm, m, y, status, cosmo, psp) \
                         default(none)
    {
      int i, j;
      double a_sf, smooth_radius;
      int local_status = *status;

      #pragma omp for
      for (j=0; j<na; j++) {
        a_sf = aa[j];
        for (i=0; i<nm; i++) {
          smooth_radius = sigmaM_m2r(cosmo, pow(10,m[i]), &local_status);
          y[j*nm + i] = log(ccl_sigmaR(cosmo, smooth_radius, a_sf,
                                       psp, &local_status));
        }
      } //end omp for
      if (local_status) {
        #pragma omp atomic write
        *status = local_status;
      }
    } //end omp parallel
  }

  gsl_spline2d *lsM = NULL;
  if (*status == 0) {
    lsM = gsl_spline2d_alloc(gsl_interp2d_bicubic, nm, na);
    if (lsM == NULL) {
      *status = CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
                                       "error allocating 2D spline\n");
    }
  }

  if(*status == 0) {
    int s2dstatus=gsl_spline2d_init(lsM, m, aa, y, nm, na);
    if (s2dstatus) {
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,
                                       "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
                                       "error initializing spline\n");
    }
  }

  if (*status == 0) {
    cosmo->computed_sigma = true;
    cosmo->data.logsigma = lsM;
  }
  else
    gsl_spline2d_free(lsM);

  free(aa);
  free(m);
  free(y);
}

/*----- ROUTINE: ccl_sigma_M -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns sigma from the sigmaM interpolation. Also computes the sigma interpolation if
necessary.
*/
double ccl_sigmaM(ccl_cosmology *cosmo, double log_halomass, double a, int *status)
{
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(cosmo,
                                     "ccl_massfunc.c: ccl_sigmaM(): "
                                     "sigma(M) spline has not been computed!");
    return NAN;
  }

  double lgsigmaM;
  int gslstatus = gsl_spline2d_eval_e(cosmo->data.logsigma, log_halomass,
                                      a, NULL, NULL, &lgsigmaM);

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_sigmaM():");
    *status |= gslstatus;
  }

  return exp(lgsigmaM);
}

/*----- ROUTINE: ccl_dlnsigM_dlogM -----
INPUT: ccl_cosmology *cosmo, double halo mass in units of Msun
TASK: returns the value of the derivative of ln(sigma^-1) with respect to log10 in halo mass.
*/
double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, double a, int *status)
{
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(cosmo,
                                     "ccl_massfunc.c: ccl_dlnsigM_dlogM(): "
                                     "sigma(M) spline has not been computed!");
    return NAN;
  }

  double dlsdlgm;
  int gslstatus = gsl_spline2d_eval_deriv_x_e(cosmo->data.logsigma,
                                              log_halomass, a,
                                              NULL, NULL, &dlsdlgm);
  if(gslstatus) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_dlnsigM_dlogM():");
    *status |= gslstatus;
  }
  return -dlsdlgm;
}

ext: c

max_stars_repo_name: Jappenn/CCL

max_issues_repo_name: Jappenn/CCL

max_forks_repo_name: Jappenn/CCL


text:
#ifndef S3D_CV_UTILITIES_CV_H
#define S3D_CV_UTILITIES_CV_H

#include "s3d/image/image.h"
#include "s3d/utilities/random_color_generator.h"
#include "s3d/video/video_frame.h"
#include "s3d/video/video_types.h"
#include "s3d/utilities/eigen.h"

#include <opencv/cxeigen.hpp>
#include <opencv2/opencv.hpp>

#include <cstdint>

#include <gsl/gsl>

#include <cassert>
#include <cmath>
#include <tuple>

namespace s3d {

inline Image<uint8_t> cv2image(cv::Mat cvImg) {
  Image<uint8_t> img(Size(cvImg.cols, cvImg.rows));
  cv::cv2eigen(cvImg, img);
  return img;
}

inline cv::Mat image2cv(const Image<uint8_t>& img) {
  cv::Mat cvImg;
  cv::eigen2cv(img, cvImg);
  return cvImg;
}

inline int nbChannelsToMatType(int nbChannels) {
  int imgType = CV_8U;
  switch (nbChannels) {
    case 1:
      imgType = CV_8U;
      break;
    case 2:
      imgType = CV_8UC2;
      break;
    case 3:
      imgType = CV_8UC3;
      break;
    case 4:
      imgType = CV_8UC4;
      break;
    default:
      break;
  }
  return imgType;
}

inline cv::Mat dataToMat(Size size, VideoPixelFormat pixelFormat, const gsl::span<const uint8_t>& data) {
  auto nbChannels = static_cast<int>(VideoFrame::NumBytesPerPixel(pixelFormat));
  return cv::Mat(size.getHeight(),
                 size.getWidth(),
                 nbChannelsToMatType(nbChannels),
                 const_cast<uchar*>(data.data()));
}
// does not copy data, to copy use resizeMat(...).clone();
inline void resizeMat(gsl::not_null<cv::Mat*> mat, float colsRatio, float rowsRatio) {
  float cols = mat->cols;
  float rows = mat->rows;
  // todo: use linear if upsampling, inter_area is good for down sampling
  cv::resize(*mat, *mat, cv::Size2f(cols * colsRatio, rows * rowsRatio), cv::INTER_AREA);
}

inline void resizeMat(gsl::not_null<cv::Mat*> mat, float ratio) {
  resizeMat(mat, ratio, ratio);
}

inline void displayInNewWindow(const std::string& name, cv::InputArray src) {
  cv::namedWindow(name, cv::WINDOW_AUTOSIZE);
  cv::imshow(name, src);
}

template <class T>
void displayMatches(const std::string& title,
                    const cv::Mat& left,
                    const cv::Mat& right,
                    const std::vector<Eigen::Matrix<T, 2, 1>>& featuresLeft,
                    const std::vector<Eigen::Matrix<T, 2, 1>>& featuresRight) {
  cv::Mat leftRight(left.rows, left.cols * 2, left.type());
  cv::hconcat(left, right, leftRight);

  // convert to color
  if (left.depth() == 1) {
    cv::cvtColor(leftRight, leftRight, cv::COLOR_GRAY2BGR);
  }

  // for random color
  RandomColorGenerator<cv::Scalar> colorGenerator;

  for (size_t i = 0; i < featuresLeft.size(); ++i) {
    // random color
    cv::Scalar color = colorGenerator.randomColor();

    // discretize features
    const Eigen::Vector2i featureLeft = featuresLeft[i].template cast<int>();
    const Eigen::Vector2i featureRight =
        featuresRight[i].template cast<int>() + Eigen::Vector2i(left.cols, 0.0f);

    const cv::Point ptLeft(featureLeft.x(), featureLeft.y());
    const cv::Point ptRight(featureRight.x(), featureRight.y());

    cv::line(leftRight, ptLeft, ptRight, color, 1, CV_AA);
    cv::circle(leftRight, ptLeft, 3, color, 1, CV_AA);
    cv::circle(leftRight, ptRight, 3, color, 1, CV_AA);
  }

  displayInNewWindow(title, leftRight);
}

template <class T>
std::vector<cv::Point_<T>> eigenPointsToCV(std::vector<Eigen::Matrix<T, 2, 1>> points) {
  std::vector<cv::Point_<T>> pointsCV;
  for (auto point : points) {
    pointsCV.emplace_back(point.x(), point.y());
  }
  return pointsCV;
}

template <class T, int N, int M>
cv::Matx<T, N, M> eigenMatToCV(Eigen::Matrix<T, N, M> mat) {
  cv::Matx<T, 3, 3> matCV;
  cv::eigen2cv(mat, matCV);
  return matCV;
};

template <typename T>
static float distancePointLine(const cv::Point_<T> point, const cv::Vec<T, 3>& line) {
  // Line is given as a*x + b*y + c = 0
  return std::abs(line(0) * point.x + line(1) * point.y + line(2)) /
         std::sqrt(line(0) * line(0) + line(1) * line(1));
}

/**
 * \brief Compute and draw the epipolar lines in two images
 *      associated to each other by a fundamental matrix
 *
 * \param title     Title of the window to display
 * \param F         Fundamental matrix
 * \param img1      First image
 * \param img2      Second image
 * \param points1   Set of points in the first image
 * \param points2   Set of points in the second image matching to the first set
 * \param inlierDistance      Points with a high distance to the epipolar lines are
 *                not displayed. If it is negative, all points are displayed
 **/
template <typename T1, typename T2>
static std::tuple<cv::Mat, cv::Mat> drawEpipolarLines(const cv::Matx<T1, 3, 3> F,
                                                      const cv::Mat& img1,
                                                      const cv::Mat& img2,
                                                      const std::vector<cv::Point_<T2>>& points1,
                                                      const std::vector<cv::Point_<T2>>& points2,
                                                      const float inlierDistance = -1) {
  assert(img1.size() == img2.size() && img1.type() == img2.type());
  assert(points1.size() == points2.size());
  assert(points1.size() > 0);

  cv::Mat outImgLeft(img1.rows, img1.cols, CV_8UC3);
  cv::Mat outImgRight(img1.rows, img1.cols, CV_8UC3);

  // for color drawing
  if (img1.type() == CV_8U) {
    cv::cvtColor(img1, outImgLeft, CV_GRAY2BGR);
    cv::cvtColor(img2, outImgRight, CV_GRAY2BGR);
  } else {
    img1.copyTo(outImgLeft);
    img2.copyTo(outImgRight);
  }

  // compute epilines for image 1 and 2
  std::vector<cv::Vec<T2, 3>> epilines1, epilines2;
  cv::computeCorrespondEpilines(points1, 1, F, epilines1);
  cv::computeCorrespondEpilines(points2, 2, F, epilines2);

  assert(points1.size() == points2.size() && points2.size() == epilines1.size() &&
         epilines1.size() == epilines2.size());

  s3d::RandomColorGenerator<cv::Scalar> colorGenerator;
  for (size_t i = 0; i < points1.size(); i++) {
    if (inlierDistance > 0) {
      if (distancePointLine(points1[i], epilines2[i]) > inlierDistance ||
          distancePointLine(points2[i], epilines1[i]) > inlierDistance) {
        // The point match is no inlier
        continue;
      }
    }

    // epipolar lines of the 1st point set are drawn in the 2nd image and vice-versa
    cv::Scalar color = colorGenerator.randomColor();

    cv::line(
        outImgRight,
        cv::Point(0, -epilines1[i][2] / epilines1[i][1]),
        cv::Point(img1.cols, -(epilines1[i][2] + epilines1[i][0] * img1.cols) / epilines1[i][1]),
        color,
        1,
        CV_AA);
    cv::circle(outImgLeft, points1[i], 3, color, -1, CV_AA);

    cv::line(
        outImgLeft,
        cv::Point(0, -epilines2[i][2] / epilines2[i][1]),
        cv::Point(img2.cols, -(epilines2[i][2] + epilines2[i][0] * img2.cols) / epilines2[i][1]),
        color,
        1,
        CV_AA);
    cv::circle(outImgRight, points2[i], 3, color, -1, CV_AA);
  }

  return std::tuple<cv::Mat, cv::Mat>(outImgLeft, outImgRight);
}

/**
 * \brief Compute and draw the epipolar lines in two images
 *      associated to each other by a fundamental matrix
 *
 * \param title     Title of the window to display
 * \param F         Fundamental matrix
 * \param img1      First image
 * \param img2      Second image
 * \param points1   Set of points in the first image
 * \param points2   Set of points in the second image matching to the first set
 **/
template <typename T1, typename T2>
static void displayEpipolarLines(const std::string& title,
                                 const cv::Matx<T1, 3, 3> F,
                                 const cv::Mat& img1,
                                 const cv::Mat& img2,
                                 const std::vector<cv::Point_<T2>>& points1,
                                 const std::vector<cv::Point_<T2>>& points2) {
  assert(img1.size() == img2.size() && img1.type() == img2.type());

  cv::Mat imgLeftEpilines, imgRightEpilines;
  std::tie(imgLeftEpilines, imgRightEpilines) =
      s3d::drawEpipolarLines(F, img1, img2, points1, points2);

  // horizontal concat
  cv::Mat outImg(img1.rows, img1.cols * 2, CV_8UC3);
  cv::hconcat(imgLeftEpilines, imgRightEpilines, outImg);

  cv::imshow(title, outImg);
  cv::waitKey(1);
}

}  // namespace s3d

#endif  // S3D_CV_UTILITIES_CV_H


text:
/**
 * VGMTrans (c) - 2002-2021
 * Licensed under the zlib license
 * See the included LICENSE for more information
 */

#pragma once

#include <vector>
#include <fluidsynth.h>
#include <gsl-lite.hpp>
#include <QObject>
#include <QString>

class VGMColl;

class MusicPlayer : public QObject {
  Q_OBJECT
public:
  static auto &the() {
    static MusicPlayer instance;
    return instance;
  }

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

  ~MusicPlayer();

  /**
   * Toggles the status of the player
   * @returns true if the player is playing
   */
  bool toggle();
  /**
   * Stops the player (unloads resources)
   */
  void stop();
  /**
   * Moves the player position
   * @param position relative to song start
   */
  void seek(int position);

  /**
   * Returns the title of the song currently playing
   * @return the song title
   */
  [[nodiscard]] QString songTitle() const;

  /**
   * Loads a VGMColl for playback
   * @param collection
   * @return true if data was loaded correctly
   */
  bool playCollection(VGMColl *collection);

  /**
   * Checks whether the player is playing
   * @return true if playing
   */
  [[nodiscard]] bool playing() const;
  /**
   * The number of MIDI ticks elapsed since the player was started
   * @return number of ticks relative to song start
   */
  [[nodiscard]] int elapsedTicks() const;
  /**
   * The total number of MIDI ticks in the song
   * @return total ticks in the song
   */
  [[nodiscard]] int totalTicks() const;

  /**
   * Gets all the audio driver the player supports
   * @return the driver list
   */
  [[nodiscard]] std::vector<const char *> getAvailableDrivers() const;
  /**
   * Changes the current audio driver and restarts playing
   * @param driver_name
   * @return true if the driver changed
   */
  bool setAudioDriver(const char *driver_name);

  /**
   * Checks the value of a setting of the player corresponds to the passed value
   * @param setting
   * @param value
   * @return true if the value matches
   */
  [[nodiscard]] bool checkSetting(const char *setting, const char *value) const;
  /**
   * Sets a parameter of the player to the corrisponding passed value
   * @param setting
   * @param value
   */
  void updateSetting(const char *setting, const char *value);

signals:
  void statusChange(bool playing);
  void playbackPositionChanged(int current, int max);

private:
  MusicPlayer();

  fluid_settings_t *m_settings = nullptr;
  fluid_synth_t *m_synth = nullptr;
  fluid_audio_driver_t *m_active_driver = nullptr;
  fluid_player_t *m_active_player = nullptr;
  VGMColl *m_active_coll = nullptr;

  void makeSettings();
  void makeSynth();
  void makePlayer();
};

max_issues_repo_name: DubiousDoggo/vgmtrans

max_forks_repo_name: DubiousDoggo/vgmtrans


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_chebyshev.h>

int gsl_cheb_calc_integ(gsl_cheb_series * integ, const gsl_cheb_series * f)
{
  const size_t n = f->order + 1;
  const double con = 0.25 * (f->b - f->a);

  if(integ->order != f->order)
    {
      GSL_ERROR ("order of chebyshev series must be equal", GSL_ENOMEM);
    }

  /* set the other parameters in the chebyshev struct */

  integ->a = f->a;
  integ->b = f->b;

  /* FIXME:  should probably set integ->f[] as well */

  if(n == 1) {
    integ->c[0] = 0.;
  }
  else if(n == 2) {
    integ->c[1] = con * f->c[0];
    integ->c[0] = 2.0 * integ->c[1];
  }
  else {
    double sum = 0.0;
    double fac = 1.0;
    size_t i;
    for(i=1; i<=n-2; i++) {
      integ->c[i] = con * (f->c[i-1] - f->c[i+1])/((double)i);
      sum += fac * integ->c[i];
      fac = -fac;
    }
    integ->c[n-1] = con * f->c[n-2]/(n-1.0);
    sum += fac * integ->c[n-1];
    integ->c[0] = 2.0 * sum;
  }

  return GSL_SUCCESS;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdlib.h>
#include <stdio.h>
#include "linterp.h"
#include "qspline.h"
#include "cspline.h"
#include <gsl/gsl_spline.h>

#define NUM_DATAPOINTS 9

/** Read a txt file with two columns seperated by tab and read all the items
 * into the arrays x and y.
 **/
void read_data(double *x, double *y){
    FILE *infile = fopen("xy_data.txt", "r");
    for (int i = 0; i < NUM_DATAPOINTS; ++i){
        int items = fscanf(infile, "%lf\t%lf\n", &x[i], &y[i]);
        if (items == EOF){
            break;
        }
    }
    fclose(infile);
}

void exerciseA(){
    FILE *outfile = fopen("out1.txt", "w");
    double *x, *y;
    x = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    y = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    read_data(x, y);

    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    gsl_spline *spline = gsl_spline_alloc (gsl_interp_linear, NUM_DATAPOINTS);

    gsl_spline_init (spline, x, y, NUM_DATAPOINTS);

    double val = 0.1;
    do {
        double interp_val = linterp(NUM_DATAPOINTS, x, y, val);
        double interp_val_integ = linterp_integ(NUM_DATAPOINTS, x, y, val);
        double y_gsl = gsl_spline_eval(spline, val, acc);
        fprintf(outfile, "%f\t%f\t%f\t%f\n", val, interp_val, interp_val_integ, y_gsl);
        val += 0.1;
    } while (val < 8);
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);
    free(x), free(y);
    fclose(outfile);
}

void exerciseB(){
    FILE *outfile = fopen("out2.txt", "w");
    double *x, *y;
    x = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    y = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    read_data(x, y);

    double val = 0.1;
    do {
        double interp_val = qspline_eval(NUM_DATAPOINTS, x, y, val);
        double interp_val_integ = qspline_integ(NUM_DATAPOINTS, x, y, val);
        double interp_val_deriv = qspline_deriv(NUM_DATAPOINTS, x, y, val);
        fprintf(outfile, "%f\t%f\t%f\t%f\n", val, interp_val, interp_val_integ, interp_val_deriv);
        val += 0.1;
    } while (val < 8);
    free(x), free(y);
    fclose(outfile);
}

void exerciseC(){
    FILE *outfile = fopen("out3.txt", "w");
    double *x, *y;
    x = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    y = (double *) malloc(sizeof(double *) * NUM_DATAPOINTS);
    read_data(x, y);

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

    gsl_spline_init (spline, x, y, NUM_DATAPOINTS);

    double val = 0.1;
    do {
        double interp_val = cspline_eval(NUM_DATAPOINTS, x, y, val);
        double interp_val_integ = cspline_integ(NUM_DATAPOINTS, x, y, val);
        double y_gsl = gsl_spline_eval(spline, val, acc);
        fprintf(outfile, "%f\t%f\t%f\t%f\n", val, interp_val, interp_val_integ, y_gsl);
        val += 0.1;
    } while (val < 8);
    gsl_spline_free (spline);
    gsl_interp_accel_free (acc);
    free(x), free(y);
    fclose(outfile);
}

int main(){
    exerciseA();
    exerciseB();
    exerciseC();
    return 0;
}


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:59 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 32 */

/*
 * This function contains 156 FP additions, 42 FP multiplications,
 * (or, 140 additions, 26 multiplications, 16 fused multiply/add),
 * 46 stack variables, and 64 memory accesses
 */
static const fftw_real K195090322 = FFTW_KONST(+0.195090322016128267848284868477022240927691618);
static const fftw_real K980785280 = FFTW_KONST(+0.980785280403230449126182236134239036973933731);
static const fftw_real K555570233 = FFTW_KONST(+0.555570233019602224742830813948532874374937191);
static const fftw_real K831469612 = FFTW_KONST(+0.831469612302545237078788377617905756738560812);
static const fftw_real K382683432 = FFTW_KONST(+0.382683432365089771728459984030398866761344562);
static const fftw_real K923879532 = FFTW_KONST(+0.923879532511286756128183189396788286822416626);
static const fftw_real K707106781 = FFTW_KONST(+0.707106781186547524400844362104849039284835938);

/*
 * 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_32(const fftw_real *input, fftw_real *real_output, fftw_real *imag_output, int istride, int real_ostride, int imag_ostride)
{
     fftw_real tmp7;
     fftw_real tmp135;
     fftw_real tmp31;
     fftw_real tmp83;
     fftw_real tmp14;
     fftw_real tmp148;
     fftw_real tmp34;
     fftw_real tmp82;
     fftw_real tmp123;
     fftw_real tmp143;
     fftw_real tmp71;
     fftw_real tmp99;
     fftw_real tmp126;
     fftw_real tmp144;
     fftw_real tmp76;
     fftw_real tmp100;
     fftw_real tmp22;
     fftw_real tmp136;
     fftw_real tmp38;
     fftw_real tmp80;
     fftw_real tmp29;
     fftw_real tmp137;
     fftw_real tmp41;
     fftw_real tmp79;
     fftw_real tmp116;
     fftw_real tmp140;
     fftw_real tmp54;
     fftw_real tmp96;
     fftw_real tmp119;
     fftw_real tmp141;
     fftw_real tmp59;
     fftw_real tmp97;
     ASSERT_ALIGNED_DOUBLE;
     {
	  fftw_real tmp1;
	  fftw_real tmp2;
	  fftw_real tmp3;
	  fftw_real tmp4;
	  fftw_real tmp5;
	  fftw_real tmp6;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp1 = input[0];
	  tmp2 = input[16 * istride];
	  tmp3 = tmp1 + tmp2;
	  tmp4 = input[8 * istride];
	  tmp5 = input[24 * istride];
	  tmp6 = tmp4 + tmp5;
	  tmp7 = tmp3 + tmp6;
	  tmp135 = tmp3 - tmp6;
	  tmp31 = tmp1 - tmp2;
	  tmp83 = tmp4 - tmp5;
     }
     {
	  fftw_real tmp10;
	  fftw_real tmp32;
	  fftw_real tmp13;
	  fftw_real tmp33;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp8;
	       fftw_real tmp9;
	       fftw_real tmp11;
	       fftw_real tmp12;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp8 = input[4 * istride];
	       tmp9 = input[20 * istride];
	       tmp10 = tmp8 + tmp9;
	       tmp32 = tmp8 - tmp9;
	       tmp11 = input[28 * istride];
	       tmp12 = input[12 * istride];
	       tmp13 = tmp11 + tmp12;
	       tmp33 = tmp11 - tmp12;
	  }
	  tmp14 = tmp10 + tmp13;
	  tmp148 = tmp13 - tmp10;
	  tmp34 = K707106781 * (tmp32 + tmp33);
	  tmp82 = K707106781 * (tmp33 - tmp32);
     }
     {
	  fftw_real tmp63;
	  fftw_real tmp121;
	  fftw_real tmp75;
	  fftw_real tmp122;
	  fftw_real tmp66;
	  fftw_real tmp124;
	  fftw_real tmp69;
	  fftw_real tmp125;
	  fftw_real tmp70;
	  fftw_real tmp72;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp61;
	       fftw_real tmp62;
	       fftw_real tmp73;
	       fftw_real tmp74;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp61 = input[31 * istride];
	       tmp62 = input[15 * istride];
	       tmp63 = tmp61 - tmp62;
	       tmp121 = tmp61 + tmp62;
	       tmp73 = input[7 * istride];
	       tmp74 = input[23 * istride];
	       tmp75 = tmp73 - tmp74;
	       tmp122 = tmp73 + tmp74;
	  }
	  {
	       fftw_real tmp64;
	       fftw_real tmp65;
	       fftw_real tmp67;
	       fftw_real tmp68;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp64 = input[3 * istride];
	       tmp65 = input[19 * istride];
	       tmp66 = tmp64 - tmp65;
	       tmp124 = tmp64 + tmp65;
	       tmp67 = input[27 * istride];
	       tmp68 = input[11 * istride];
	       tmp69 = tmp67 - tmp68;
	       tmp125 = tmp67 + tmp68;
	  }
	  tmp123 = tmp121 + tmp122;
	  tmp143 = tmp121 - tmp122;
	  tmp70 = K707106781 * (tmp66 + tmp69);
	  tmp71 = tmp63 + tmp70;
	  tmp99 = tmp63 - tmp70;
	  tmp126 = tmp124 + tmp125;
	  tmp144 = tmp125 - tmp124;
	  tmp72 = K707106781 * (tmp69 - tmp66);
	  tmp76 = tmp72 - tmp75;
	  tmp100 = tmp75 + tmp72;
     }
     {
	  fftw_real tmp18;
	  fftw_real tmp36;
	  fftw_real tmp21;
	  fftw_real tmp37;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp16;
	       fftw_real tmp17;
	       fftw_real tmp19;
	       fftw_real tmp20;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp16 = input[2 * istride];
	       tmp17 = input[18 * istride];
	       tmp18 = tmp16 + tmp17;
	       tmp36 = tmp16 - tmp17;
	       tmp19 = input[10 * istride];
	       tmp20 = input[26 * istride];
	       tmp21 = tmp19 + tmp20;
	       tmp37 = tmp19 - tmp20;
	  }
	  tmp22 = tmp18 + tmp21;
	  tmp136 = tmp18 - tmp21;
	  tmp38 = (K923879532 * tmp36) - (K382683432 * tmp37);
	  tmp80 = (K382683432 * tmp36) + (K923879532 * tmp37);
     }
     {
	  fftw_real tmp25;
	  fftw_real tmp39;
	  fftw_real tmp28;
	  fftw_real tmp40;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp23;
	       fftw_real tmp24;
	       fftw_real tmp26;
	       fftw_real tmp27;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp23 = input[30 * istride];
	       tmp24 = input[14 * istride];
	       tmp25 = tmp23 + tmp24;
	       tmp39 = tmp23 - tmp24;
	       tmp26 = input[6 * istride];
	       tmp27 = input[22 * istride];
	       tmp28 = tmp26 + tmp27;
	       tmp40 = tmp26 - tmp27;
	  }
	  tmp29 = tmp25 + tmp28;
	  tmp137 = tmp25 - tmp28;
	  tmp41 = (K923879532 * tmp39) + (K382683432 * tmp40);
	  tmp79 = (K382683432 * tmp39) - (K923879532 * tmp40);
     }
     {
	  fftw_real tmp46;
	  fftw_real tmp114;
	  fftw_real tmp58;
	  fftw_real tmp115;
	  fftw_real tmp49;
	  fftw_real tmp117;
	  fftw_real tmp52;
	  fftw_real tmp118;
	  fftw_real tmp53;
	  fftw_real tmp55;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp44;
	       fftw_real tmp45;
	       fftw_real tmp56;
	       fftw_real tmp57;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp44 = input[istride];
	       tmp45 = input[17 * istride];
	       tmp46 = tmp44 - tmp45;
	       tmp114 = tmp44 + tmp45;
	       tmp56 = input[9 * istride];
	       tmp57 = input[25 * istride];
	       tmp58 = tmp56 - tmp57;
	       tmp115 = tmp56 + tmp57;
	  }
	  {
	       fftw_real tmp47;
	       fftw_real tmp48;
	       fftw_real tmp50;
	       fftw_real tmp51;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp47 = input[5 * istride];
	       tmp48 = input[21 * istride];
	       tmp49 = tmp47 - tmp48;
	       tmp117 = tmp47 + tmp48;
	       tmp50 = input[29 * istride];
	       tmp51 = input[13 * istride];
	       tmp52 = tmp50 - tmp51;
	       tmp118 = tmp50 + tmp51;
	  }
	  tmp116 = tmp114 + tmp115;
	  tmp140 = tmp114 - tmp115;
	  tmp53 = K707106781 * (tmp49 + tmp52);
	  tmp54 = tmp46 + tmp53;
	  tmp96 = tmp46 - tmp53;
	  tmp119 = tmp117 + tmp118;
	  tmp141 = tmp118 - tmp117;
	  tmp55 = K707106781 * (tmp52 - tmp49);
	  tmp59 = tmp55 - tmp58;
	  tmp97 = tmp58 + tmp55;
     }
     {
	  fftw_real tmp15;
	  fftw_real tmp30;
	  fftw_real tmp131;
	  fftw_real tmp132;
	  fftw_real tmp133;
	  fftw_real tmp134;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp15 = tmp7 + tmp14;
	  tmp30 = tmp22 + tmp29;
	  tmp131 = tmp15 + tmp30;
	  tmp132 = tmp116 + tmp119;
	  tmp133 = tmp123 + tmp126;
	  tmp134 = tmp132 + tmp133;
	  real_output[8 * real_ostride] = tmp15 - tmp30;
	  imag_output[8 * imag_ostride] = tmp133 - tmp132;
	  real_output[16 * real_ostride] = tmp131 - tmp134;
	  real_output[0] = tmp131 + tmp134;
     }
     {
	  fftw_real tmp113;
	  fftw_real tmp129;
	  fftw_real tmp128;
	  fftw_real tmp130;
	  fftw_real tmp120;
	  fftw_real tmp127;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp113 = tmp7 - tmp14;
	  tmp129 = tmp29 - tmp22;
	  tmp120 = tmp116 - tmp119;
	  tmp127 = tmp123 - tmp126;
	  tmp128 = K707106781 * (tmp120 + tmp127);
	  tmp130 = K707106781 * (tmp127 - tmp120);
	  real_output[12 * real_ostride] = tmp113 - tmp128;
	  real_output[4 * real_ostride] = tmp113 + tmp128;
	  imag_output[4 * imag_ostride] = tmp129 + tmp130;
	  imag_output[12 * imag_ostride] = tmp130 - tmp129;
     }
     {
	  fftw_real tmp139;
	  fftw_real tmp155;
	  fftw_real tmp149;
	  fftw_real tmp151;
	  fftw_real tmp146;
	  fftw_real tmp150;
	  fftw_real tmp154;
	  fftw_real tmp156;
	  fftw_real tmp138;
	  fftw_real tmp147;
	  ASSERT_ALIGNED_DOUBLE;
	  tmp138 = K707106781 * (tmp136 + tmp137);
	  tmp139 = tmp135 + tmp138;
	  tmp155 = tmp135 - tmp138;
	  tmp147 = K707106781 * (tmp137 - tmp136);
	  tmp149 = tmp147 - tmp148;
	  tmp151 = tmp148 + tmp147;
	  {
	       fftw_real tmp142;
	       fftw_real tmp145;
	       fftw_real tmp152;
	       fftw_real tmp153;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp142 = (K923879532 * tmp140) + (K382683432 * tmp141);
	       tmp145 = (K923879532 * tmp143) - (K382683432 * tmp144);
	       tmp146 = tmp142 + tmp145;
	       tmp150 = tmp145 - tmp142;
	       tmp152 = (K923879532 * tmp141) - (K382683432 * tmp140);
	       tmp153 = (K382683432 * tmp143) + (K923879532 * tmp144);
	       tmp154 = tmp152 + tmp153;
	       tmp156 = tmp153 - tmp152;
	  }
	  real_output[14 * real_ostride] = tmp139 - tmp146;
	  real_output[2 * real_ostride] = tmp139 + tmp146;
	  imag_output[6 * imag_ostride] = tmp149 + tmp150;
	  imag_output[10 * imag_ostride] = tmp150 - tmp149;
	  imag_output[2 * imag_ostride] = tmp151 + tmp154;
	  imag_output[14 * imag_ostride] = tmp154 - tmp151;
	  real_output[10 * real_ostride] = tmp155 - tmp156;
	  real_output[6 * real_ostride] = tmp155 + tmp156;
     }
     {
	  fftw_real tmp95;
	  fftw_real tmp111;
	  fftw_real tmp110;
	  fftw_real tmp112;
	  fftw_real tmp102;
	  fftw_real tmp106;
	  fftw_real tmp105;
	  fftw_real tmp107;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp93;
	       fftw_real tmp94;
	       fftw_real tmp108;
	       fftw_real tmp109;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp93 = tmp31 - tmp34;
	       tmp94 = tmp80 + tmp79;
	       tmp95 = tmp93 + tmp94;
	       tmp111 = tmp93 - tmp94;
	       tmp108 = (K831469612 * tmp97) - (K555570233 * tmp96);
	       tmp109 = (K555570233 * tmp99) + (K831469612 * tmp100);
	       tmp110 = tmp108 + tmp109;
	       tmp112 = tmp109 - tmp108;
	  }
	  {
	       fftw_real tmp98;
	       fftw_real tmp101;
	       fftw_real tmp103;
	       fftw_real tmp104;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp98 = (K831469612 * tmp96) + (K555570233 * tmp97);
	       tmp101 = (K831469612 * tmp99) - (K555570233 * tmp100);
	       tmp102 = tmp98 + tmp101;
	       tmp106 = tmp101 - tmp98;
	       tmp103 = tmp41 - tmp38;
	       tmp104 = tmp83 + tmp82;
	       tmp105 = tmp103 - tmp104;
	       tmp107 = tmp104 + tmp103;
	  }
	  real_output[13 * real_ostride] = tmp95 - tmp102;
	  real_output[3 * real_ostride] = tmp95 + tmp102;
	  imag_output[5 * imag_ostride] = tmp105 + tmp106;
	  imag_output[11 * imag_ostride] = tmp106 - tmp105;
	  imag_output[3 * imag_ostride] = tmp107 + tmp110;
	  imag_output[13 * imag_ostride] = tmp110 - tmp107;
	  real_output[11 * real_ostride] = tmp111 - tmp112;
	  real_output[5 * real_ostride] = tmp111 + tmp112;
     }
     {
	  fftw_real tmp43;
	  fftw_real tmp91;
	  fftw_real tmp90;
	  fftw_real tmp92;
	  fftw_real tmp78;
	  fftw_real tmp86;
	  fftw_real tmp85;
	  fftw_real tmp87;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp35;
	       fftw_real tmp42;
	       fftw_real tmp88;
	       fftw_real tmp89;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp35 = tmp31 + tmp34;
	       tmp42 = tmp38 + tmp41;
	       tmp43 = tmp35 + tmp42;
	       tmp91 = tmp35 - tmp42;
	       tmp88 = (K980785280 * tmp59) - (K195090322 * tmp54);
	       tmp89 = (K195090322 * tmp71) + (K980785280 * tmp76);
	       tmp90 = tmp88 + tmp89;
	       tmp92 = tmp89 - tmp88;
	  }
	  {
	       fftw_real tmp60;
	       fftw_real tmp77;
	       fftw_real tmp81;
	       fftw_real tmp84;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp60 = (K980785280 * tmp54) + (K195090322 * tmp59);
	       tmp77 = (K980785280 * tmp71) - (K195090322 * tmp76);
	       tmp78 = tmp60 + tmp77;
	       tmp86 = tmp77 - tmp60;
	       tmp81 = tmp79 - tmp80;
	       tmp84 = tmp82 - tmp83;
	       tmp85 = tmp81 - tmp84;
	       tmp87 = tmp84 + tmp81;
	  }
	  real_output[15 * real_ostride] = tmp43 - tmp78;
	  real_output[real_ostride] = tmp43 + tmp78;
	  imag_output[7 * imag_ostride] = tmp85 + tmp86;
	  imag_output[9 * imag_ostride] = tmp86 - tmp85;
	  imag_output[imag_ostride] = tmp87 + tmp90;
	  imag_output[15 * imag_ostride] = tmp90 - tmp87;
	  real_output[9 * real_ostride] = tmp91 - tmp92;
	  real_output[7 * real_ostride] = tmp91 + tmp92;
     }
}

fftw_codelet_desc fftw_real2hc_32_desc =
{
     "fftw_real2hc_32",
     (void (*)()) fftw_real2hc_32,
     32,
     FFTW_FORWARD,
     FFTW_REAL2HC,
     706,
     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:
// A rectangular block of data of type uint8_t.  A transparent block
// cache is used for each instance of this class if necessary, so the
// image can be much bigger than will fit in memory, and pixels can
// still be accessed and written efficiently, provided subsequent
// accesses are spatially correlated.  A variety of useful methods are
// implemented (filtering, subsetting, interpolating, etc.)
//
// Don't try to access the same instance concurrently.  Split your
// images up into separate instances if you must parallelize things.
//
// For many methods, arguments of type ssize_t are used, but are not
// allowed to be negative.  This is to help prevent people from
// shooting themselves in the foor by accidently passing negative
// values which don't yield warnings and can't be caught with
// assertions.

#ifndef UINT8_IMAGE_H
#define UINT8_IMAGE_H

#ifndef solaris
#  include <stdint.h>
#endif
#include "asf_meta.h"
#include <stdio.h>
#include <sys/types.h>

#include <glib.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_histogram.h>

// sometimes we don't have these - choose conservative values
#ifndef SSIZE_MAX
#define SSIZE_MAX 32767
#endif

#ifndef UINT8_MAX
#define UINT8_MAX 255
#endif

// Instance structure.  Everything here is private and need not be
// used or understood by client code, except for the size_x and size_y
// fields.
typedef struct {
  size_t size_x, size_y;	// Image dimensions.
  size_t cache_space;		// Memory cache space in bytes.
  size_t cache_area;		// Memory cache area in pixels.
  size_t tile_size;		// Tile size in pixels on a side.
  size_t cache_size_in_tiles;	// Number of tiles in cache.
  size_t tile_count_x;		// Number of tiles in image in x direction.
  size_t tile_count_y;		// Number of tiles in image in y direction.
  size_t tile_count;		// Total number of tiles in image.
  size_t tile_area;             // Area of a tile, in pixels.
  uint8_t *cache;		// Memory cache.
  uint8_t **tile_addresses;	// Addresss of individual tiles in the cache.
  GQueue *tile_queue;		// Queue of tile offsets kept in load order.
  FILE *tile_file;              // File with tiles stored contiguously.
  GString *tile_file_name;  // Filename of the tile file
} UInt8Image;

///////////////////////////////////////////////////////////////////////////////
//
// Creating New Instances
//
///////////////////////////////////////////////////////////////////////////////

// Thaw out a previously frozen instance (produced with
// uint8_image_freeze) using data pointed to by file_pointer.  Frozen
// instances aren't portable between platforms.  After thawing
// file_pointer points to the data immediately following the data from
// the thawed instance (or to the end of the file).
UInt8Image *
uint8_image_thaw (FILE *file_pointer);

// Create a new image filled with zero pixel values.
UInt8Image *
uint8_image_new (ssize_t size_x, ssize_t size_y);

// Create a new image with pixels initialized to value.
UInt8Image *
uint8_image_new_with_value (ssize_t size_x, ssize_t size_y, uint8_t value);

// Create a new image from memory.  This pixels are assumed to be
// layed out in memory in the usual way, i.e. contiguous rows of
// pixels in the x direction are contiguous in memory.
UInt8Image *
uint8_image_new_from_memory (ssize_t size_x, ssize_t size_y, uint8_t *buffer);

// Create a new independent copy of model.
UInt8Image *
uint8_image_copy (UInt8Image *model);

// Form reduced resolution version of the model.  The scale_factor
// must be positive and odd.  The new image will be round ((double)
// model->size_x / scale_factor) pixels by round ((double)
// model->size_y / scale_factor) pixels.  Scaling is performed by
// averaging blocks of pixels together, using odd pixel reflection
// around the image edges (see the description of the apply_kernel
// method).  The average for a block of pixels is first computed as a
// double precision floating point number, then rounded with the C99
// round() function.  The upper and leftmost blocks of pixels averaged
// together are always centered at the 0 index in the direction in
// question, so reflection is always used for these edges.  Whether
// reflection is used for the right and lower edges depends on the
// relationship between the model dimensions and the scale factor.
UInt8Image *
uint8_image_new_from_model_scaled (UInt8Image *model, ssize_t scale_factor);

// Create a new image by copying the portion of model with upper left
// corner at model coordinates (x, y), width size_x, and height
// size_y.
UInt8Image *
uint8_image_new_subimage (UInt8Image *model, ssize_t x, ssize_t y,
			  ssize_t size_x, ssize_t size_y);

// Create a new image from data at byte offset in file.  The pixel
// layout in the file is assumed to be the same as for the
// uint8_image_new_from_memory method.  The byte order of individual
// pixels in the file should be byte_order.
UInt8Image *
uint8_image_new_from_file (ssize_t size_x, ssize_t size_y, const char *file,
			   off_t offset);

// The method is like new_from_file, but takes a file pointer instead
// of a file name, and the offset argument is with respect to the
// current position in the file_pointer stream.
UInt8Image *
uint8_image_new_from_file_pointer (ssize_t size_x, ssize_t size_y,
				   FILE *file_pointer, off_t offset);

// Form a low quality reduced resolution version of the
// original_size_x by original_size_y image in file.  The new image
// will be size_x by size_y pixels.  This method is like new_from_file
// method, but gets its data by sampling in each dimension using
// bilinear interpolation, and rounding the interpolated values using
// the C99 round() routine.  This is a decent way of forming quick
// thumbnails of images, or of scaling images down just slightly (by a
// factor of say 1.5 or less), but not much else.
UInt8Image *
uint8_image_new_from_file_scaled (ssize_t size_x, ssize_t size_y,
				  ssize_t original_size_x,
				  ssize_t original_size_y,
				  const char *file, off_t offset);

// The function that does it all, generating an instance of UInt8Image
// from a file and the metadata
UInt8Image *
uint8_image_new_from_metadata(meta_parameters *meta, const char *file);

// For multi-band imagery the previous function needs to be more specific.
UInt8Image *
uint8_image_band_new_from_metadata(meta_parameters *meta,
				   int band, const char *file);

///////////////////////////////////////////////////////////////////////////////
//
// Getting and Setting Image Pixels and Regions
//
///////////////////////////////////////////////////////////////////////////////

// Get pixel at 0-indexed position x, y.  A cache is used so pixel
// access does not usually involve the hard disk, provided subsequent
// pixel lookups are spatially close together.  You can probably get
// away with trating this function just as if everything was in
// memory.  For details, see the cache control methods below.
//
// The x and y arguments should always be positive, ssize_t is used
// only so people can't as easily shoot themselves in the foot by
// accidently supplying negative arguments which get promoted in some
// strange way.
uint8_t
uint8_image_get_pixel (UInt8Image *self, ssize_t x, ssize_t y);

// Set pixel at 0-indexed position x, y to value.  A cache is used to
// make this fast, as for the uint8_image_get_pixel method.
void
uint8_image_set_pixel (UInt8Image *self, ssize_t x, ssize_t y, uint8_t value);

// Get rectangular image region of size_x, size_y having upper left
// corner at x, y and copy it into already allocated buffer.  There is
// not necessarily any caching help for this method, i.e. it may
// always involve disk access and always be slow.
void
uint8_image_get_region (UInt8Image *self, ssize_t x, ssize_t y,
			ssize_t size_x, ssize_t size_y, uint8_t *buffer);

// This method is analogous to uint8_image_get_region.
void
uint8_image_set_region (UInt8Image *self, size_t x, size_t y, size_t size_x,
			size_t size_y, uint8_t *buffer);

// Get a full row of pixels, copying the data into already allocated
// buffer.  This method will be fast on the average for calls with
// sequential row numbers.
void
uint8_image_get_row (UInt8Image *self, size_t row, uint8_t *buffer);

// Get a pixel, performing odd reflection at image edges if the pixel
// indicies fall outside the image.  See the description of the
// apply_kernel method for an explanation of reflection.
uint8_t
uint8_image_get_pixel_with_reflection (UInt8Image *self, ssize_t x, ssize_t y);

///////////////////////////////////////////////////////////////////////////////
//
// Image Analysis and Statistics
//
///////////////////////////////////////////////////////////////////////////////

// Default mask value when figuring image stats
#define UINT8_IMAGE_DEFAULT_MASK (0.0)

// Finds the minimum and maximum pixel values in the image, and the
// mean and standard deviation of all pixels.  This function considers
// every pixel in the image when mask is NAN, otherwise it discounts
// all values within .00000000001 of the mask
void
uint8_image_statistics (UInt8Image *self, uint8_t *min, uint8_t *max,
			double *mean, double *standard_deviation,
			gboolean use_mask_value, uint8_t mask_value);

// Does the same thing as uint8_image_statistics() except that the
// statistics are calculated for a particular band selection within
// a UInt8Image which contains multiple bands or channels in sequential
// order.  Returns 1 on error, otherwise 0
int
uint8_image_band_statistics (UInt8Image *self, meta_stats *stats,
                             int line_count, int band_no,
                             gboolean use_mask_value, uint8_t mask_value);

// This method works like the statistics method, except values in the
// interval [interval_start, interval_end] are not considered at all
// for the purposes of determining any of the outputs.
void
uint8_image_statistics_with_mask_interval (UInt8Image *self, uint8_t *min,
					   uint8_t *max, double *mean,
					   double *standard_deviation,
					   uint8_t interval_start,
					   uint8_t interval_end);

// Compute an efficient estimate of the mean and standard deviation of
// the pixels in the image, by sampling every stride th pixel in each
// dimension, beginning with pixel (0, 0). If the mask is a non-NAN
// value, this function will discount all values within .00000000001
// of the mask
void
uint8_image_approximate_statistics (UInt8Image *self, size_t stride,
                                    double *mean, double *standard_deviation,
				    gboolean use_mask_value,
				    uint8_t mask_value);

// This method is a logical combination of the
// statistics_with_mask_interval and approximate_statistics methods.
void
uint8_image_approximate_statistics_with_mask_interval
  (UInt8Image *self, size_t stride, double *mean, double *standard_deviation,
   uint8_t interval_start, uint8_t interval_end);

// Creates a gsl_histogram with bin_count bins evenly spaced over the
// interval [min, max].  This function considers every pixel in the
// image.
gsl_histogram *
uint8_image_gsl_histogram (UInt8Image *self, double min, double max,
                           size_t bin_count);

///////////////////////////////////////////////////////////////////////////////
//
// Kernels, Interpolation, and Sampling
//
///////////////////////////////////////////////////////////////////////////////

// Apply kernel centerd at pixel x, y and return the value.  The
// kernel matrix must be have equal odd dimensions.  The values in the
// kernel are multiplied by the pixels, and the sum of the products
// returned.  When part of the kernel would fall outside the image
// extents, the values used for the out-of-image pixels are the mirror
// images of the corresponding in-image pixels, with the edge pixels
// not duplicated, i.e. reflection about the middle of the edge pixels
// is used.
double
uint8_image_apply_kernel (UInt8Image *self, ssize_t x, ssize_t y,
			  gsl_matrix *kern);

// Type used to specify whether disk files should be in big or little
// endian byte order.

// Sample method types.  These dictate how nearby pixels are
// considered when we want to find the approximate value for a point
// which falls between pixel indicies.
typedef enum {
  // Nearest pixel.
  UINT8_IMAGE_SAMPLE_METHOD_NEAREST_NEIGHBOR,
  // Linearly weited average of four nearest pixels
  UINT8_IMAGE_SAMPLE_METHOD_BILINEAR,
  // Bicubic spline interpolation (which consideres the nearest 16 pixels).
  UINT8_IMAGE_SAMPLE_METHOD_BICUBIC
} uint8_image_sample_method_t;

double
uint8_image_sample (UInt8Image *self, double x, double y,
		    uint8_image_sample_method_t sample_method);

///////////////////////////////////////////////////////////////////////////////
//
// Comparing Images
//
///////////////////////////////////////////////////////////////////////////////

// Return true iff self and other have identical sizes and pixels.
gboolean
uint8_image_equals (UInt8Image *self, UInt8Image *other);

///////////////////////////////////////////////////////////////////////////////
//
// Manipulating Images
//
///////////////////////////////////////////////////////////////////////////////

// Flip an image about a horizontal line through the center of the image
void uint8_image_flip_y(UInt8Image *self);

// Flip an image about a vertical line through the center of the image
void uint8_image_flip_x(UInt8Image *self);

///////////////////////////////////////////////////////////////////////////////
//
// Storing Images in Files
//
///////////////////////////////////////////////////////////////////////////////

// Store instance self at position pointed to by file_pointer, for
// later retrieval using uint8_image_thaw.  The serialized version of
// self is not portable between platforms.
void
uint8_image_freeze (UInt8Image *self, FILE *file_pointer);

// Store image pixels in file.  The image is stored in the usual
// order, i.e. contiguous rows of pixels in the x direction are stored
// contiguously in memory.  Individual pixels are stored in byte order
// byte_order.  Returns 0 on success, nonzero on error.
int
uint8_image_store (UInt8Image *self, const char *file);

int
uint8_image_store_ext(UInt8Image *self, const char *file, int append_flag);

int
uint8_image_band_store(UInt8Image *self, const char *file,
		       meta_parameters *meta, int append_flag);

///////////////////////////////////////////////////////////////////////////////
//
// Exporting Images in Various Image File Formats
//
///////////////////////////////////////////////////////////////////////////////

// Export image to fila as a gray scaled jpeg image, with largest
// dimension no larger than max_dimension.  The max_dimension argument
// must be less than or equal to the largest dimension of the image.
// The image may be scaled st its largest dimension is considerably
// less than max_dimension.  Scaling is performed by averaging blocks
// of pixels together, using odd pixel reflection around the image
// edges (see the description of the apply_kernel method), and
// rounding pixel averages with the C99 round() routine.  If
// use_mask_value is true, then in determining the image statistics
// (mean and standard deviation), values equal to mask are not
// considered.  This routine slurps the whole image into memory, so
// beware.  Returns 0 on success, nonzero on error.
int
uint8_image_export_as_jpeg (UInt8Image *self, const char *file,
			    size_t max_dimension, gboolean use_mask_value,
			    uint8_t mask_value);

///////////////////////////////////////////////////////////////////////////////
//
// Controlling Image Data Caching
//
// It probably isn't necessary to use these methods.  They are
// provided largely to make it clear how the cache works.  The major
// tunable parameter is the size of the in-memory cache to use.
//
// When a new image is created, the following steps are performed:
//
//      1. The image is divided up into square tiles st two full rows
//         or columns of tiles will fit in the memory cache.
//
//      2. A copy of the image is created on disk with the memory
//         layout rearranged st individual tiles are contiguous in
//         memory.  This allows tiles to be quickly retrieved later.
//
// When a pixel is accessed (read or set), the following happens:
//
//      1. If the pixel is in a tile already loaded into the cache,
//         it is simply fetched or set.
//
//      2. Otherwise, the tile containing the pixel is loaded,
//         possibly displacing an already loaded tile, and then the
//         pixel is fetched or set.  The tile displaced is the one
//         loaded longest ago (there is no most-recently-accessed
//         heuristic, as this would make pixel access too slow).
//
// Thus, using a larger memory cache will result in larger tiles being
// used, and fewer tile loads being needed.  In general, the default
// behavior is pretty good, but if you know will be performing lots of
// widely (but not too widely) scattered accesses, you might want to
// make it bigger.
//
///////////////////////////////////////////////////////////////////////////////

// Get the image memory cache size setting, in bytes.  Note that this
// is the memory cache used per image, not the class-wide cache usage.
// If you will have a lot of objects instantiated simultaneously, you
// may find it necessary to use a smaller cache for each image.
size_t
uint8_image_get_cache_size (UInt8Image *self);

// Set the image memory cache to size bytes.  Changing the cache size
// requires the tiling to be recomputed, the on-disk tile cache to be
// regenerated, and the in memory cache to be flushed, so its slow.
void
uint8_image_set_cache_size (UInt8Image *self, size_t size);

///////////////////////////////////////////////////////////////////////////////
//
// Freeing Instances
//
///////////////////////////////////////////////////////////////////////////////

// Destroy self.
void
uint8_image_free (UInt8Image *self);

#endif // #ifndef UINT8_IMAGE_H

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


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

double f(double d, void* fdata){
    double * double_fdata = (double*) fdata;

    double par_obs = double_fdata[0];
    double var_par = double_fdata[1];
    double DM = double_fdata[2];
    double var_DM = double_fdata[3];
    double par_offset = double_fdata[4];
    double L = double_fdata[5];

    //printf("Data = %f, %f, %f, %f, %f\n", d, par_corrected, var_par, DM, var_DM);
    return pow(d, 2) * exp(-0.5 * (2*d/L + pow(par_obs - (1./d + par_offset), 2.)/var_par + pow(DM - (5*log10(d) - 5), 2.)/var_DM));
}

double integrate_likelihood(const double xmin, const double xmax, double *fdata, double rel_tolerance) {
    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);

    gsl_function F;
    F.function = &f;
    F.params = fdata;

    double result, err;

    gsl_integration_qag (&F, xmin, xmax, 0, rel_tolerance, 1000, 1, w, &result, &err);

    gsl_integration_workspace_free (w);

    return result;
}


text:
/* Siconos-Numerics, Copyright INRIA 2005-2018.
 * Siconos is a program dedicated to modeling, simulation and control
 * of non smooth dynamical systems.
 * Siconos is a 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 2 of the License, or
 * (at your option) any later version.
 * Siconos 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 Siconos; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact: Vincent ACARY, siconos-team@lists.gforge.inria.fr
 */

#ifndef SiconosLAPACKE_H
#define SiconosLAPACKE_H

#include "SiconosBlas.h"
#include <stdlib.h>
#include <assert.h>
#include <stdio.h>

// -------- Headers and routines naming conventions for the different Lapack implementations --------

// --- Intel MKL Header ---
#if defined(HAS_MKL_LAPACKE)
#include <mkl_lapacke.h>
#else
// Standard lapacke header
#include <lapacke.h>
#endif

// Name of the routines
#define LAPACK_NAME(N) LAPACKE_##N

#define LA_TRANS 'T'
#define LA_NOTRANS 'N'
#define LA_UP 'U'
#define LA_LO 'L'
#define LA_NONUNIT 'N'
#define LA_UNIT 'U'

#define INTEGER(X) X
#define INTEGERP(X) X
#define CHAR(X) X

#ifndef lapack_int
#define lapack_int int
#endif

// --- DGESVD ---
#if defined(HAS_LAPACK_DGESVD)
#define WRAP_DGESVD(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,INFO)      \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12)
#else
#define WRAP_DGESVD(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,A10,A11,A12,INFO)      \
  fprintf(stderr, "Your lapack version misses dgesvd function.\n");
#endif

// --- DGETRS ---
#define WRAP_DGETRS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)   \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8)

// --- DGESV ---
#define WRAP_DGESV(F,A1,A2,A3,A4,A5,A6,A7,INFO)   \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7)

// --- DGELS ---
#if defined(HAS_LAPACK_DGELS)
#define WRAP_DGELS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)    \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8)
#else
#define WRAP_DGELS(F,A1,A2,A3,A4,A5,A6,A7,A8,INFO)                      \
  fprintf(stderr, "Your lapack version misses dgels function.\n");
#endif

// --- DGETRI ---
#define WRAP_DGETRI(F,A1,A2,A3,A4,INFO)         \
  INFO = F(CblasColMajor,A1,A2,A3,A4)

// --- DPOTRF ---
#define WRAP_DPOTRF(F,A1,A2,A3,A4,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4)

// --- DGETRF ---
#define WRAP_DGETRF(F,A1,A2,A3,A4,A5,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5)

// --- DTRTRS ---
#if defined(HAS_LAPACK_DTRTRS)
#define WRAP_DTRTRS(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,INFO)  \
  INFO = F(CblasColMajor,A1,A2,A3,A4,A5,A6,A7,A8,A9)
#else
#define WRAP_DTRTRS(F,A1,A2,A3,A4,A5,A6,A7,A8,A9,INFO)                \
  fprintf(stderr, "Your lapack version misses dtrtrs function.\n");
#endif

#endif // SICONOSLAPACKE_H

max_stars_repo_name: ljktest/siconos

max_issues_repo_name: ljktest/siconos

max_forks_repo_name: ljktest/siconos


text:
/* rng/ranmar.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 RANMAR lagged fibonacci generator of Marsaglia, Zaman
   and Tsang.  The sequence is a series of 24-bit integers, x_n,

   x_n = (y_n - c_n + 2^24) mod 2^24

   where,

   y_n = (y_{n-97) - y_{n-33} + 2^24) mod 2^24
   c_n = (c_{n-1} - 7654321 + 2^24 - 3) mod (2^24 - 3)

   The period of this generator is 2^144.

   The generator provides about 900 million different subsequences
   each of length O(10^30). Thus each seed up to 900,000,000 gives an
   independent sequence.

   Although it was good in its day this generator now has known
   statistical defects and has been superseded by RANLUX.

   From: F. James, "A Review of Pseudorandom number generators",
   Computer Physics Communications 60, 329 (1990).

   G. Marsaglia, A. Zaman and W.W. Tsang, Stat. Prob. Lett. 9, 35 (1990)  */

static inline unsigned long int ranmar_get (void *vstate);
static double ranmar_get_double (void *vstate);
static void ranmar_set (void *state, unsigned long int s);

static const unsigned long int two24 = 16777216;	/* 2^24 */

typedef struct
  {
    unsigned int i;
    unsigned int j;
    long int carry;
    unsigned long int u[97];
  }
ranmar_state_t;

static inline unsigned long int
ranmar_get (void *vstate)
{
  ranmar_state_t *state = (ranmar_state_t *) vstate;

  unsigned int i = state->i;
  unsigned int j = state->j;
  long int carry = state->carry;

  long int delta = state->u[i] - state->u[j];

  if (delta < 0)
    delta += two24 ;

  state->u[i] = delta;

  if (i == 0)
    {
      i = 96;
    }
  else
    {
      i--;
    }

  state->i = i;

  if (j == 0)
    {
      j = 96;
    }
  else
    {
      j--;
    }

  state->j = j;

  carry += - 7654321 ;

  if (carry < 0)
    carry += two24 - 3;

  state->carry = carry ;

  delta += - carry ;

  if (delta < 0)
    delta += two24 ;

  return delta;
}

static double
ranmar_get_double (void *vstate)
{
  return ranmar_get (vstate) / 16777216.0 ;
}

static void
ranmar_set (void *vstate, unsigned long int s)
{
  ranmar_state_t *state = (ranmar_state_t *) vstate;

  unsigned long int ij = s / 30082 ;
  unsigned long int kl = s % 30082 ;

  int i = (ij / 177) % 177 + 2 ;
  int j = (ij % 177) + 2 ;
  int k = (kl / 169) % 178 + 1 ;
  int l = (kl % 169) ;

  int a, b;

  for (a = 0; a < 97; a++)
    {
      unsigned long int sum = 0 ;
      unsigned long int t = two24 ;

      for (b = 0; b < 24; b++)
	{
	  unsigned long int m = (((i * j) % 179) * k) % 179 ;
	  i = j ;
	  j = k ;
	  k = m ;
	  l = (53 * l + 1) % 169 ;
	  t >>= 1 ;

	  if ((l * m) % 64 >= 32)
	    sum += t ;
	}

      state->u[a] = sum ;
    }

  state->i = 96;
  state->j = 32;
  state->carry = 362436 ;

}

static const gsl_rng_type ranmar_type =
{"ranmar",			/* name */
 0x00ffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranmar_state_t),
 &ranmar_set,
 &ranmar_get,
 &ranmar_get_double};

const gsl_rng_type *gsl_rng_ranmar = &ranmar_type;


text:
/* linalg/tridiag.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2004 Gerard Jungman,
 * Brian Gough, David Necas
 *
 * 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.
 */

/* Author: G. Jungman */

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include "tridiag.h"
#include <gsl/gsl_linalg.h>

/* for description of method see [Engeln-Mullges + Uhlig, p. 92]
 *
 *     diag[0]  offdiag[0]             0   .....
 *  offdiag[0]     diag[1]    offdiag[1]   .....
 *           0  offdiag[1]       diag[2]
 *           0           0    offdiag[2]   .....
 */
static
int
solve_tridiag(
  const double diag[], size_t d_stride,
  const double offdiag[], size_t o_stride,
  const double b[], size_t b_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double *gamma = (double *) malloc (N * sizeof (double));
  double *alpha = (double *) malloc (N * sizeof (double));
  double *c = (double *) malloc (N * sizeof (double));
  double *z = (double *) malloc (N * sizeof (double));

  if (gamma == 0 || alpha == 0 || c == 0 || z == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      size_t i, j;

      /* Cholesky decomposition
         A = L.D.L^t
         lower_diag(L) = gamma
         diag(D) = alpha
       */
      alpha[0] = diag[0];
      gamma[0] = offdiag[0] / alpha[0];

      for (i = 1; i < N - 1; i++)
        {
          alpha[i] = diag[d_stride * i] - offdiag[o_stride*(i - 1)] * gamma[i - 1];
          gamma[i] = offdiag[o_stride * i] / alpha[i];
        }

      if (N > 1)
        {
          alpha[N - 1] = diag[d_stride * (N - 1)] - offdiag[o_stride*(N - 2)] * gamma[N - 2];
        }

      /* update RHS */
      z[0] = b[0];
      for (i = 1; i < N; i++)
        {
          z[i] = b[b_stride * i] - gamma[i - 1] * z[i - 1];
        }
      for (i = 0; i < N; i++)
        {
          c[i] = z[i] / alpha[i];
        }

      /* backsubstitution */
      x[x_stride * (N - 1)] = c[N - 1];
      if (N >= 2)
        {
          for (i = N - 2, j = 0; j <= N - 2; j++, i--)
            {
              x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)];
            }
        }

      status = GSL_SUCCESS;
    }

  if (z != 0)
    free (z);
  if (c != 0)
    free (c);
  if (alpha != 0)
    free (alpha);
  if (gamma != 0)
    free (gamma);

  return status;
}

/* plain gauss elimination, only not bothering with the zeroes
 *
 *       diag[0]  abovediag[0]             0   .....
 *  belowdiag[0]       diag[1]  abovediag[1]   .....
 *             0  belowdiag[1]       diag[2]
 *             0             0  belowdiag[2]   .....
 */
static
int
solve_tridiag_nonsym(
  const double diag[], size_t d_stride,
  const double abovediag[], size_t a_stride,
  const double belowdiag[], size_t b_stride,
  const double rhs[], size_t r_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double *alpha = (double *) malloc (N * sizeof (double));
  double *z = (double *) malloc (N * sizeof (double));

  if (alpha == 0 || z == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      size_t i, j;

      /* Bidiagonalization (eliminating belowdiag)
         & rhs update
         diag' = alpha
         rhs' = z
       */
      alpha[0] = diag[0];
      z[0] = rhs[0];

      for (i = 1; i < N; i++)
        {
          const double t = belowdiag[b_stride*(i - 1)]/alpha[i-1];
          alpha[i] = diag[d_stride*i] - t*abovediag[a_stride*(i - 1)];
          z[i] = rhs[r_stride*i] - t*z[i-1];
          /* FIXME!!! */
          if (alpha[i] == 0) {
            status = GSL_EZERODIV;
            goto solve_tridiag_nonsym_END;
          }
        }

      /* backsubstitution */
      x[x_stride * (N - 1)] = z[N - 1]/alpha[N - 1];
      if (N >= 2)
        {
          for (i = N - 2, j = 0; j <= N - 2; j++, i--)
            {
              x[x_stride * i] = (z[i] - abovediag[a_stride*i] * x[x_stride * (i + 1)])/alpha[i];
            }
        }

      status = GSL_SUCCESS;
    }

solve_tridiag_nonsym_END:
  if (z != 0)
    free (z);
  if (alpha != 0)
    free (alpha);

  return status;
}

/* for description of method see [Engeln-Mullges + Uhlig, p. 96]
 *
 *      diag[0]  offdiag[0]             0   .....  offdiag[N-1]
 *   offdiag[0]     diag[1]    offdiag[1]   .....
 *            0  offdiag[1]       diag[2]
 *            0           0    offdiag[2]   .....
 *          ...         ...
 * offdiag[N-1]         ...
 *
 */
static
int
solve_cyc_tridiag(
  const double diag[], size_t d_stride,
  const double offdiag[], size_t o_stride,
  const double b[], size_t b_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double * delta = (double *) malloc (N * sizeof (double));
  double * gamma = (double *) malloc (N * sizeof (double));
  double * alpha = (double *) malloc (N * sizeof (double));
  double * c = (double *) malloc (N * sizeof (double));
  double * z = (double *) malloc (N * sizeof (double));

  if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      size_t i, j;
      double sum = 0.0;

      /* factor */

      if (N == 1)
        {
          x[0] = b[0] / diag[0];
          return GSL_SUCCESS;
        }

      alpha[0] = diag[0];
      gamma[0] = offdiag[0] / alpha[0];
      delta[0] = offdiag[o_stride * (N-1)] / alpha[0];

      for (i = 1; i < N - 2; i++)
        {
          alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
          gamma[i] = offdiag[o_stride * i] / alpha[i];
          delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
        }

      for (i = 0; i < N - 2; i++)
        {
          sum += alpha[i] * delta[i] * delta[i];
        }

      alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];

      gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];

      alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2];

      /* update */
      z[0] = b[0];
      for (i = 1; i < N - 1; i++)
        {
          z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
        }
      sum = 0.0;
      for (i = 0; i < N - 2; i++)
        {
          sum += delta[i] * z[i];
        }
      z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
      for (i = 0; i < N; i++)
        {
          c[i] = z[i] / alpha[i];
        }

      /* backsubstitution */
      x[x_stride * (N - 1)] = c[N - 1];
      x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
      if (N >= 3)
        {
          for (i = N - 3, j = 0; j <= N - 3; j++, i--)
            {
              x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
            }
        }

      status = GSL_SUCCESS;
    }

  if (z != 0)
    free (z);
  if (c != 0)
    free (c);
  if (alpha != 0)
    free (alpha);
  if (gamma != 0)
    free (gamma);
  if (delta != 0)
    free (delta);

  return status;
}

/* solve following system w/o the corner elements and then use
 * Sherman-Morrison formula to compensate for them
 *
 *        diag[0]  abovediag[0]             0   .....  belowdiag[N-1]
 *   belowdiag[0]       diag[1]  abovediag[1]   .....
 *              0  belowdiag[1]       diag[2]
 *              0             0  belowdiag[2]   .....
 *            ...           ...
 * abovediag[N-1]           ...
 */
static
int solve_cyc_tridiag_nonsym(
  const double diag[], size_t d_stride,
  const double abovediag[], size_t a_stride,
  const double belowdiag[], size_t b_stride,
  const double rhs[], size_t r_stride,
  double x[], size_t x_stride,
  size_t N)
{
  int status;
  double *alpha = (double *) malloc (N * sizeof (double));
  double *zb = (double *) malloc (N * sizeof (double));
  double *zu = (double *) malloc (N * sizeof (double));
  double *w = (double *) malloc (N * sizeof (double));
  double beta;

  if (alpha == 0 || zb == 0 || zu == 0 || w == 0)
    {
      status = GSL_ENOMEM;
    }
  else
    {
      /* Bidiagonalization (eliminating belowdiag)
         & rhs update
         diag' = alpha
         rhs' = zb
         rhs' for Aq=u is zu
       */
      zb[0] = rhs[0];
      if (diag[0] != 0) beta = -diag[0]; else beta = 1;
      {
        const double q = 1 - abovediag[0]*belowdiag[0]/(diag[0]*diag[d_stride]);
        if (fabs(q/beta) > 0.5 && fabs(q/beta) < 2) {
          beta *= (fabs(q/beta) < 1) ? 0.5 : 2;
        }
      }
      zu[0] = beta;
      alpha[0] = diag[0] - beta;

      {
        size_t i;
        for (i = 1; i+1 < N; i++)
        {
          const double t = belowdiag[b_stride*(i - 1)]/alpha[i-1];
          alpha[i] = diag[d_stride*i] - t*abovediag[a_stride*(i - 1)];
          zb[i] = rhs[r_stride*i] - t*zb[i-1];
          zu[i] = -t*zu[i-1];
          /* FIXME!!! */
          if (alpha[i] == 0) {
            status = GSL_EZERODIV;
            goto solve_cyc_tridiag_nonsym_END;
          }
        }
      }

      {
        const size_t i = N-1;
        const double t = belowdiag[b_stride*(i - 1)]/alpha[i-1];
        alpha[i] = diag[d_stride*i]
                   - abovediag[a_stride*i]*belowdiag[b_stride*i]/beta
                   - t*abovediag[a_stride*(i - 1)];
        zb[i] = rhs[r_stride*i] - t*zb[i-1];
        zu[i] = abovediag[a_stride*i] - t*zu[i-1];
        /* FIXME!!! */
        if (alpha[i] == 0) {
          status = GSL_EZERODIV;
          goto solve_cyc_tridiag_nonsym_END;
        }
      }

      /* backsubstitution */
      {
        size_t i, j;
        w[N-1] = zu[N-1]/alpha[N-1];
        x[N-1] = zb[N-1]/alpha[N-1];
        for (i = N - 2, j = 0; j <= N - 2; j++, i--)
          {
            w[i] = (zu[i] - abovediag[a_stride*i] * w[i+1])/alpha[i];
            x[i*x_stride] = (zb[i] - abovediag[a_stride*i] * x[x_stride*(i + 1)])/alpha[i];
          }
      }

      /* Sherman-Morrison */
      {
        const double vw = w[0] + belowdiag[b_stride*(N - 1)]/beta * w[N-1];
        const double vx = x[0] + belowdiag[b_stride*(N - 1)]/beta * x[x_stride*(N - 1)];
        /* FIXME!!! */
        if (vw + 1 == 0) {
          status = GSL_EZERODIV;
          goto solve_cyc_tridiag_nonsym_END;
        }

        {
          size_t i;
          for (i = 0; i < N; i++)
            x[i] -= vx/(1 + vw)*w[i];
        }
      }

      status = GSL_SUCCESS;
    }

solve_cyc_tridiag_nonsym_END:
  if (zb != 0)
    free (zb);
  if (zu != 0)
    free (zu);
  if (w != 0)
    free (w);
  if (alpha != 0)
    free (alpha);

  return status;
}

int
gsl_linalg_solve_symm_tridiag(
  const gsl_vector * diag,
  const gsl_vector * offdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size)
    {
      GSL_ERROR ("size of diag must match rhs", GSL_EBADLEN);
    }
  else if (offdiag->size != rhs->size-1)
    {
      GSL_ERROR ("size of offdiag must match rhs-1", GSL_EBADLEN);
    }
  else if (solution->size != rhs->size)
    {
      GSL_ERROR ("size of solution must match rhs", GSL_EBADLEN);
    }
  else
    {
      return solve_tridiag(diag->data, diag->stride,
                           offdiag->data, offdiag->stride,
                           rhs->data, rhs->stride,
                           solution->data, solution->stride,
                           diag->size);
    }
}

int
gsl_linalg_solve_tridiag(
  const gsl_vector * diag,
  const gsl_vector * abovediag,
  const gsl_vector * belowdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size)
    {
      GSL_ERROR ("size of diag must match rhs", GSL_EBADLEN);
    }
  else if (abovediag->size != rhs->size-1)
    {
      GSL_ERROR ("size of abovediag must match rhs-1", GSL_EBADLEN);
    }
  else if (belowdiag->size != rhs->size-1)
    {
      GSL_ERROR ("size of belowdiag must match rhs-1", GSL_EBADLEN);
    }
  else if (solution->size != rhs->size)
    {
      GSL_ERROR ("size of solution must match rhs", GSL_EBADLEN);
    }
  else
    {
      return solve_tridiag_nonsym(diag->data, diag->stride,
                                  abovediag->data, abovediag->stride,
                                  belowdiag->data, belowdiag->stride,
                                  rhs->data, rhs->stride,
                                  solution->data, solution->stride,
                                  diag->size);
    }
}

int
gsl_linalg_solve_symm_cyc_tridiag(
  const gsl_vector * diag,
  const gsl_vector * offdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size)
    {
      GSL_ERROR ("size of diag must match rhs", GSL_EBADLEN);
    }
  else if (offdiag->size != rhs->size)
    {
      GSL_ERROR ("size of offdiag must match rhs", GSL_EBADLEN);
    }
  else if (solution->size != rhs->size)
    {
      GSL_ERROR ("size of solution must match rhs", GSL_EBADLEN);
    }
  else if (diag->size < 3)
    {
      GSL_ERROR ("size of cyclic system must be 3 or more", GSL_EBADLEN);
    }
  else
    {
      return solve_cyc_tridiag(diag->data, diag->stride,
                               offdiag->data, offdiag->stride,
                               rhs->data, rhs->stride,
                               solution->data, solution->stride,
                               diag->size);
    }
}

int
gsl_linalg_solve_cyc_tridiag(
  const gsl_vector * diag,
  const gsl_vector * abovediag,
  const gsl_vector * belowdiag,
  const gsl_vector * rhs,
  gsl_vector * solution)
{
  if(diag->size != rhs->size)
    {
      GSL_ERROR ("size of diag must match rhs", GSL_EBADLEN);
    }
  else if (abovediag->size != rhs->size)
    {
      GSL_ERROR ("size of abovediag must match rhs", GSL_EBADLEN);
    }
  else if (belowdiag->size != rhs->size)
    {
      GSL_ERROR ("size of belowdiag must match rhs", GSL_EBADLEN);
    }
  else if (solution->size != rhs->size)
    {
      GSL_ERROR ("size of solution must match rhs", GSL_EBADLEN);
    }
  else if (diag->size < 3)
    {
      GSL_ERROR ("size of cyclic system must be 3 or more", GSL_EBADLEN);
    }
  else
    {
      return solve_cyc_tridiag_nonsym(diag->data, diag->stride,
                                      abovediag->data, abovediag->stride,
                                      belowdiag->data, belowdiag->stride,
                                      rhs->data, rhs->stride,
                                      solution->data, solution->stride,
                                      diag->size);
    }
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified