text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_core.h>
#include <gbpCosmo_NFW_etc.h>
#include <gsl/gsl_sf_expint.h>

double R_vir_NFW(double M_vir, double z, int mode, cosmo_info **cosmo) {
    double c_vir;
    double R_vir;
    set_NFW_params(M_vir, z, mode, cosmo, &c_vir, &R_vir);
    return (R_vir);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
/*
 * Gsl_blas_dgemv: test gsl_blas_dgemv (matrix . vector)
 *
 * Copyright (c) 2012 Jérémie Decock
 *
 * Required: GSL library (libgsl0-dev)
 * Usage: gcc gsl_blas_dgemv.c -lgsl -lgslcblas -lm
 *    or: gcc gsl_blas_dgemv.c $(pkg-config --libs gsl)
 */

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

void main() {

    gsl_vector * v = gsl_vector_alloc(3);
    gsl_vector * r = gsl_vector_alloc(3);
    gsl_matrix * m = gsl_matrix_alloc(3, 3);

    gsl_vector_set_all(v, 2.0);
    gsl_matrix_set_all(m, 3.0);

    gsl_blas_dgemv(CblasNoTrans,
                   1.0, m, v,
                   0.0, r);

    gsl_vector_fprintf(stdout, r, "%f");

    gsl_vector_free(v);
    gsl_vector_free(r);
    gsl_matrix_free(m);

}

max_stars_repo_name: jeremiedecock/snippets

max_issues_repo_name: jeremiedecock/snippets

max_forks_repo_name: jeremiedecock/snippets


text:
/**
 *
 * @file testing_zgemm.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_zmain.h"

#undef REAL
#define COMPLEX

static int check_solution(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                          PLASMA_Complex64_t alpha, PLASMA_Complex64_t *A, int LDA,
                          PLASMA_Complex64_t *B, int LDB,
                          PLASMA_Complex64_t beta, PLASMA_Complex64_t *Cref, PLASMA_Complex64_t *Cplasma, int LDC);

int testing_zgemm(int argc, char **argv)
{
    /* Check for number of arguments*/
    if ( argc != 8) {
        USAGE("GEMM", "alpha beta M N K LDA LDB LDC",
              "   - alpha  : alpha coefficient\n"
              "   - beta   : beta coefficient\n"
              "   - M      : number of rows of matrices A and C\n"
              "   - N      : number of columns of matrices B and C\n"
              "   - K      : number of columns of matrix A / number of rows of matrix B\n"
              "   - LDA    : leading dimension of matrix A\n"
              "   - LDB    : leading dimension of matrix B\n"
              "   - LDC    : leading dimension of matrix C\n");
        return -1;
    }

    PLASMA_Complex64_t alpha = (PLASMA_Complex64_t) atol(argv[0]);
    PLASMA_Complex64_t beta = (PLASMA_Complex64_t) atol(argv[1]);
    int M     = atoi(argv[2]);
    int N     = atoi(argv[3]);
    int K     = atoi(argv[4]);
    int LDA   = atoi(argv[5]);
    int LDB   = atoi(argv[6]);
    int LDC   = atoi(argv[7]);

    double eps;
    int info_solution;
    int i, j, ta, tb;
    int LDAxK = LDA*max(M,K);
    int LDBxN = LDB*max(K,N);
    int LDCxN = LDC*N;

    PLASMA_Complex64_t *A      = (PLASMA_Complex64_t *)malloc(LDAxK*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *B      = (PLASMA_Complex64_t *)malloc(LDBxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *C      = (PLASMA_Complex64_t *)malloc(LDCxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Cinit  = (PLASMA_Complex64_t *)malloc(LDCxN*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *Cfinal = (PLASMA_Complex64_t *)malloc(LDCxN*sizeof(PLASMA_Complex64_t));

    /* Check if unable to allocate memory */
    if ((!A)||(!B)||(!Cinit)||(!Cfinal)){
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_dlamch_work('e');

    printf("\n");
    printf("------ TESTS FOR PLASMA ZGEMM ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, N);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
     *  TESTING ZGEMM
     */

    /* Initialize A, B, C */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxK, A);
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxN, B);
    LAPACKE_zlarnv_work(IONE, ISEED, LDCxN, C);

#ifdef COMPLEX
    for (ta=0; ta<3; ta++) {
        for (tb=0; tb<3; tb++) {
#else
    for (ta=0; ta<2; ta++) {
        for (tb=0; tb<2; tb++) {
#endif
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cinit[LDC*j+i] = C[LDC*j+i];
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cfinal[LDC*j+i] = C[LDC*j+i];

            /* PLASMA ZGEMM */
            PLASMA_zgemm(trans[ta], trans[tb], M, N, K, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);

            /* Check the solution */
            info_solution = check_solution(trans[ta], trans[tb], M, N, K,
                                           alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

            if (info_solution == 0) {
                printf("***************************************************\n");
                printf(" ---- TESTING ZGEMM (%s, %s) ............... PASSED !\n", transstr[ta], transstr[tb]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING ZGEMM (%s, %s) ... FAILED !\n", transstr[ta], transstr[tb]);
                printf("************************************************\n");
            }
        }
    }
#ifdef _UNUSED_
    }}
#endif
    free(A); free(B); free(C);
    free(Cinit); free(Cfinal);

    return 0;
}

/*--------------------------------------------------------------
 * Check the solution
 */

static int check_solution(PLASMA_enum transA, PLASMA_enum transB, int M, int N, int K,
                          PLASMA_Complex64_t alpha, PLASMA_Complex64_t *A, int LDA,
                          PLASMA_Complex64_t *B, int LDB,
                          PLASMA_Complex64_t beta, PLASMA_Complex64_t *Cref, PLASMA_Complex64_t *Cplasma, int LDC)
{
    int info_solution;
    double Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm, Rnorm, result;
    double eps;
    PLASMA_Complex64_t beta_const;

    double *work = (double *)malloc(max(K,max(M, N))* sizeof(double));
    int Am, An, Bm, Bn;

    beta_const  = -1.0;

    if (transA == PlasmaNoTrans) {
        Am = M; An = K;
    } else {
        Am = K; An = M;
    }
    if (transB == PlasmaNoTrans) {
        Bm = K; Bn = N;
    } else {
        Bm = N; Bn = K;
    }

    Anorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), Am, An, A,       LDA, work);
    Bnorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), Bm, Bn, B,       LDB, work);
    Cinitnorm   = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,  N,  Cref,    LDC, work);
    Cplasmanorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M,  N,  Cplasma, LDC, work);

    cblas_zgemm(CblasColMajor, (CBLAS_TRANSPOSE)transA, (CBLAS_TRANSPOSE)transB, M, N, K,
                CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Cref, LDC);

    Clapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    cblas_zaxpy(LDC * N, CBLAS_SADDR(beta_const), Cplasma, 1, Cref, 1);

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Cref, LDC, work);

    eps = LAPACKE_dlamch_work('e');

    printf("Rnorm %e, Anorm %e, Bnorm %e, Cinitnorm %e, Cplasmanorm %e, Clapacknorm %e\n",
           Rnorm, Anorm, Bnorm, Cinitnorm, Cplasmanorm, Clapacknorm);

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);
    printf("============\n");
    printf("Checking the norm of the difference against reference ZGEMM \n");
    printf("-- ||Cplasma - Clapack||_oo/((||A||_oo+||B||_oo+||C||_oo).N.eps) = %e \n",
           result);

    if (  isnan(Rnorm) || isinf(Rnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
         printf("-- The solution is suspicious ! \n");
         info_solution = 1;
    }
    else {
         printf("-- The solution is CORRECT ! \n");
         info_solution= 0 ;
    }

    free(work);

    return info_solution;
}


text:
#ifndef AMICI_VECTOR_H
#define AMICI_VECTOR_H

#include <vector>
#include <type_traits>

#include <amici/exception.h>

#include <nvector/nvector_serial.h>

#include <gsl/gsl-lite.hpp>

namespace amici {

/** Since const N_Vector is not what we want */
using const_N_Vector =
    std::add_const<typename std::remove_pointer<N_Vector>::type *>::type;

/** AmiVector class provides a generic interface to the NVector_Serial struct */
class AmiVector {
  public:
    /**
     * @brief Default constructor
     */
    AmiVector() = default;

    /** Creates an std::vector<realtype> and attaches the
     * data pointer to a newly created N_Vector_Serial.
     * Using N_VMake_Serial ensures that the N_Vector
     * module does not try to deallocate the data vector
     * when calling N_VDestroy_Serial
     * @brief empty constructor
     * @param length number of elements in vector
     */
    explicit AmiVector(const long int length)
        : vec_(static_cast<decltype(vec_)::size_type>(length), 0.0),
          nvec_(N_VMake_Serial(length, vec_.data())) {}

    /** Moves data from std::vector and constructs an nvec that points to the
     * data
     * @brief constructor from std::vector,
     * @param rvec vector from which the data will be moved
     */
    explicit AmiVector(std::vector<realtype> rvec)
        : vec_(std::move(rvec)),
          nvec_(N_VMake_Serial(static_cast<long int>(vec_.size()), vec_.data())) {}

    /** Copy data from gsl::span and constructs a vector
     * @brief constructor from gsl::span,
     * @param rvec vector from which the data will be copied
     */
    explicit AmiVector(gsl::span<realtype> rvec)
        : AmiVector(std::vector<realtype>(rvec.begin(), rvec.end())) {}

    /**
     * @brief copy constructor
     * @param vold vector from which the data will be copied
     */
    AmiVector(const AmiVector &vold) : vec_(vold.vec_) {
        nvec_ =
            N_VMake_Serial(static_cast<long int>(vold.vec_.size()), vec_.data());
    }

    /**
     * @brief move constructor
     * @param other vector from which the data will be moved
     */
    AmiVector(AmiVector&& other) noexcept : nvec_(nullptr) {
        vec_ = std::move(other.vec_);
        synchroniseNVector();
    }

    /**
     * @brief destructor
     */
    ~AmiVector();

    /**
     * @brief copy assignment operator
     * @param other right hand side
     * @return left hand side
     */
    AmiVector &operator=(AmiVector const &other);

    /**
     * @brief data accessor
     * @return pointer to data array
     */
    realtype *data();

    /**
     * @brief const data accessor
     * @return const pointer to data array
     */
    const realtype *data() const;

    /**
     * @brief N_Vector accessor
     * @return N_Vector
     */
    N_Vector getNVector();

    /**
     * @brief N_Vector accessor
     * @return N_Vector
     */
    const_N_Vector getNVector() const;

    /**
     * @brief Vector accessor
     * @return Vector
     */
    std::vector<realtype> const &getVector() const;

    /**
     * @brief returns the length of the vector
     * @return length
     */
    int getLength() const;

    /**
     * @brief resets the Vector by filling with zero values
     */
    void reset();

    /**
     * @brief changes the sign of data elements
     */
    void minus();

    /**
     * @brief sets all data elements to a specific value
     * @param val value for data elements
     */
    void set(realtype val);

    /**
     * @brief accessor to data elements of the vector
     * @param pos index of element
     * @return element
     */
    realtype &operator[](int pos);
    /**
     * @brief accessor to data elements of the vector
     * @param pos index of element
     * @return element
     */
    realtype &at(int pos);

    /**
     * @brief accessor to data elements of the vector
     * @param pos index of element
     * @return element
     */
    const realtype &at(int pos) const;

    /**
     * @brief copies data from another AmiVector
     * @param other data source
     */
    void copy(const AmiVector &other);

  private:
    /** main data storage */
    std::vector<realtype> vec_;

    /** N_Vector, will be synchronised such that it points to data in vec */
    N_Vector nvec_ {nullptr};

    /**
     * @brief reconstructs nvec such that data pointer points to vec data array
     */
    void synchroniseNVector();
};

/**
 * @brief AmiVectorArray class.
 *
 * Provides a generic interface to arrays of NVector_Serial structs
 */
class AmiVectorArray {
  public:
    /**
     * @brief Default constructor
     */
    AmiVectorArray() = default;

    /**
     * Creates an std::vector<realype> and attaches the
     * data pointer to a newly created N_VectorArray
     * using CloneVectorArrayEmpty ensures that the N_Vector
     * module does not try to deallocate the data vector
     * when calling N_VDestroyVectorArray_Serial
     * @brief empty constructor
     * @param length_inner length of vectors
     * @param length_outer number of vectors
     */
    AmiVectorArray(long int length_inner, long int length_outer);

    /**
     * @brief copy constructor
     * @param vaold object to copy from
     */
    AmiVectorArray(const AmiVectorArray &vaold);

    ~AmiVectorArray() = default;

    /**
     * @brief copy assignment operator
     * @param other right hand side
     * @return left hand side
     */
    AmiVectorArray &operator=(AmiVectorArray const &other);

    /**
     * @brief accessor to data of AmiVector elements
     * @param pos index of AmiVector
     * @return pointer to data array
     */
    realtype *data(int pos);

    /**
     * @brief const accessor to data of AmiVector elements
     * @param pos index of AmiVector
     * @return const pointer to data array
     */
    const realtype *data(int pos) const;

    /**
     * @brief accessor to elements of AmiVector elements
     * @param ipos inner index in AmiVector
     * @param jpos outer index in AmiVectorArray
     * @return element
     */
    realtype &at(int ipos, int jpos);

    /**
     * @brief const accessor to elements of AmiVector elements
     * @param ipos inner index in AmiVector
     * @param jpos outer index in AmiVectorArray
     * @return element
     */
    const realtype &at(int ipos, int jpos) const;

    /**
     * @brief accessor to NVectorArray
     * @return N_VectorArray
     */
    N_Vector *getNVectorArray();

    /**
     * @brief accessor to NVector element
     * @param pos index of corresponding AmiVector
     * @return N_Vector
     */
    N_Vector getNVector(int pos);

    /**
     * @brief const accessor to NVector element
     * @param pos index of corresponding AmiVector
     * @return N_Vector
     */
    const_N_Vector getNVector(int pos) const;

    /**
     * @brief accessor to AmiVector elements
     * @param pos index of AmiVector
     * @return AmiVector
     */
    AmiVector &operator[](int pos);

    /**
     * @brief const accessor to AmiVector elements
     * @param pos index of AmiVector
     * @return const AmiVector
     */
    const AmiVector &operator[](int pos) const;

    /**
     * @brief length of AmiVectorArray
     * @return length
     */
    int getLength() const;

    /**
     * @brief resets every AmiVector in AmiVectorArray
     */
    void reset();

    /**
     * @brief flattens the AmiVectorArray to a vector in row-major format
     * @param vec vector into which the AmiVectorArray will be flattened. Must
     * have length equal to number of elements.
     */
    void flatten_to_vector(std::vector<realtype> &vec) const;

    /**
     * @brief copies data from another AmiVectorArray
     * @param other data source
     */
    void copy(const AmiVectorArray &other);

  private:
    /** main data storage */
    std::vector<AmiVector> vec_array_;

    /**
     * N_Vector array, will be synchronised such that it points to
     * respective elements in the vec_array
     */
    std::vector<N_Vector> nvec_array_;
};

} // namespace amici

namespace gsl {
/**
 * @brief Create span from N_Vector
 * @param nv
 * @return
 */
inline span<realtype> make_span(N_Vector nv)
{
    return span<realtype>(N_VGetArrayPointer(nv), N_VGetLength_Serial(nv));
}
} // namespace gsl

#endif /* AMICI_VECTOR_H */

ext: h

max_stars_repo_name: lcontento/AMICI

max_issues_repo_name: lcontento/AMICI

max_forks_repo_name: lcontento/AMICI


text:
#ifndef BAYESIAN_MERGE_H
#define BAYESIAN_MERGE_H

#include <vector>
#include <multi/GaussianDistribution.h>
#include <gsl/gsl_linalg.h>

using namespace std;

namespace openworld {
  class BayesianMerge {
  public:
    // modifies means in-place to assure maximum likelihood
    static void mergeGaussianHierarchy(vector<GaussianDistribution*> lows, GaussianDistribution* high) {
      vector<double> A, b;
      for (unsigned ii = 0; ii < lows.size(); ii++) {
        for (unsigned jj = 0; jj < lows.size(); jj++) {
          if (ii == jj)
            A.push_back(1 / lows[ii]->getVariance() + 1 / high->getVariance());
          else
            A.push_back(1 / high->getVariance());
        }

        b.push_back(lows[ii]->getMean() / lows[ii]->getVariance() + high->getMean() / high->getVariance());
      }

      gsl_matrix_view Amat = gsl_matrix_view_array(A.data(), lows.size(), lows.size());
      gsl_vector_view bmat = gsl_vector_view_array(b.data(), lows.size());
      gsl_vector* x = gsl_vector_alloc(lows.size());

      int s;

      gsl_permutation* p = gsl_permutation_alloc(lows.size());

      gsl_linalg_LU_decomp(&Amat.matrix, p, &s);
      gsl_linalg_LU_solve(&Amat.matrix, p, &bmat.vector, x);

      double meansum = 0;
      for (unsigned ii = 0; ii < lows.size(); ii++) {
        meansum += gsl_vector_get(x, ii);
        lows[ii]->setMean(gsl_vector_get(x, ii));
      }

      high->setMean(meansum);

      gsl_permutation_free(p);
      gsl_vector_free(x);
    }
  };
}

#endif

max_stars_repo_name: jrising/openworld

max_issues_repo_name: jrising/openworld

max_forks_repo_name: jrising/openworld


text:
/*
 *  Copyright 2008-2014 NVIDIA Corporation
 *
 *  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.
 */

#include <cusp/complex.h>
#include <cblas.h>

#define CUSP_CBLAS_COMPLEX_DEFS_1(FUNC_MACRO)                                               \
  FUNC_MACRO(cusp::complex<float> , float , c)                                              \
  FUNC_MACRO(cusp::complex<double>, double, z)

#define CUSP_CBLAS_COMPLEX_DEFS_2(FUNC_MACRO)                                               \
  FUNC_MACRO(cusp::complex<float> , float , sc)                                             \
  FUNC_MACRO(cusp::complex<double>, double, dz)

#define CUSP_CBLAS_COMPLEX_DEFS_3(FUNC_MACRO)                                               \
  FUNC_MACRO(cusp::complex<float> , float , cs)                                             \
  FUNC_MACRO(cusp::complex<double>, double, zd)

#define CUSP_CBLAS_COMPLEX_AMAX(T,V,name)                                                   \
  template<int dummy>                                                                       \
  int amax( const int n, const T* X, const int incX )                                       \
  {                                                                                \
    return cblas_i##name##amax(n, (const V*) X, incX);                                      \
  }

#define CUSP_CBLAS_COMPLEX_ASUM(T,V,name)                                                   \
  template<int dummy>                                                                       \
  V asum( const int n, const T* X, const int incX )                                         \
  {                                                                                \
    return cblas_##name##asum(n, (const V*) X, incX);                                       \
  }

#define CUSP_CBLAS_COMPLEX_AXPY(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void axpy( const int n, const T alpha, const T* X, const int incX, T* Y, const int incY ) \
  {                                                                                \
    cblas_##name##axpy(n, (const V*) &alpha, (const V*) X, incX, (V*) Y, incY);             \
  }

#define CUSP_CBLAS_COMPLEX_COPY(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void copy( const int n, const T* X, const int incX, T* Y, const int incY )                \
  {                                                                                \
    cblas_##name##copy(n, (const V*) X, incX, (V*) Y, incY);                                \
  }

// #define CUSP_CBLAS_COMPLEX_DOTC(T,V,name)
//   T dotc( const int n, const T* X, const int incX, const T* Y, const int incY )
//   {
//     typedef typename thrust::detail::eval_if<
//       thrust::detail::is_same<V,float>::value,
//           thrust::detail::identity_<openblas_complex_float>,
//           thrust::detail::identity_<openblas_complex_double>
//     >::type Cmplx;
//     Cmplx z;
//     cblas_##name##dotc_sub(n, (const V*) X, incX, (const V*) Y, incY, &z);
//     return T(z.real, z.imag);
//   }
//
// #define CUSP_CBLAS_COMPLEX_DOTU(T,V,name)
//   T dotu( const int n, const T* X, const int incX, const T* Y, const int incY )
//   {
//     typedef typename thrust::detail::eval_if<
//       thrust::detail::is_same<V,float>::value,
//           thrust::detail::identity_<openblas_complex_float>,
//           thrust::detail::identity_<openblas_complex_double>
//     >::type Cmplx;
//     Cmplx z;
//     cblas_##name##dotu_sub(n, (const V*) X, incX, (const V*) Y, incY, &z);
//     return T(z.real, z.imag);
//   }

#define CUSP_CBLAS_COMPLEX_NRM2(T,V,name)                                                   \
  template<int dummy>                                                                       \
  V nrm2( const int n, const T* X, const int incX )                                         \
  {                                                                                \
    return cblas_##name##nrm2(n, (const V*) X, incX);                                       \
  }

#define CUSP_CBLAS_COMPLEX_SCAL(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void scal( const int n, const V alpha, T* X, const int incX )                             \
  {                                                                                \
    cblas_##name##scal(n, alpha, (V*) X, incX);                                             \
  }

#define CUSP_CBLAS_COMPLEX_SWAP(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void swap( const int n, T* X, const int incX, T* Y, const int incY )                      \
  {                                                                                \
    cblas_##name##swap(n, (V*) X, incX, (V*) Y, incY);                                      \
  }

#define CUSP_CBLAS_COMPLEX_GEMV(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void gemv(  CBLAS_ORDER order,  CBLAS_TRANSPOSE trans,                                    \
           int m, int n, T alpha, const T* A, int lda,                                      \
           const T* x, int incx, T beta, T* y, int incy)                                    \
{                                                                                \
    cblas_##name##gemv(order, trans, m, n, (const V*) &alpha, (const V*) A, lda,            \
                       (const V*) x, incx, (const V*) &beta, (V*) y, incy);                 \
}

#define CUSP_CBLAS_COMPLEX_GERC(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void ger(  CBLAS_ORDER order, int m, int n, T alpha, const T* x, int incx,                \
            const T* y, int incy, T* A, int lda)                                            \
{                                                                                \
    cblas_##name##gerc(order, m, n, (const V*) &alpha,                                      \
                      (const V*) x, incx, (const V*) y, incy,                               \
                      (V*) A, lda);                                                         \
}

#define CUSP_CBLAS_COMPLEX_HEMV(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void hemv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
             int n, T alpha, const T* A, int lda,                                           \
             const T* x, int incx, T beta, T* y, int incy)                                  \
{                                                                                \
    cblas_##name##hemv(order, uplo, n, (const V*) &alpha, (const V*) A, lda,                \
                       (const V*) x, incx, (const V*) &beta, (V*) y, incy);                 \
}

#define CUSP_CBLAS_COMPLEX_HER(T,V,name)                                                    \
  template<int dummy>                                                                       \
  void her(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                           \
            int n, const V alpha, const T* x, int incx, T* A, int lda)                      \
{                                                                                \
    cblas_##name##her(order, uplo, n, alpha,                                                \
                      (const V*) x, incx, (V*) A, lda);                                     \
}

#define CUSP_CBLAS_COMPLEX_TRMV(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void trmv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int n, const T* A, int lda, T* x, int incx)                                    \
{                                                                                \
    cblas_##name##trmv(order, uplo, trans, diag, n,                                         \
                       (const V*) A, lda, (V*) x, incx);                                    \
}

#define CUSP_CBLAS_COMPLEX_TRSV(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void trsv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int n, const T* A, int lda, T* x, int incx)                                    \
{                                                                                \
    cblas_##name##trsv(order, uplo, trans, diag, n,                                         \
                       (const V*) A, lda, (V*) x, incx);                                    \
}

#define CUSP_CBLAS_COMPLEX_GEMM(T,V,name)                                                   \
  template<int dummy>                                                                       \
void gemm(  CBLAS_ORDER order,                                                              \
            CBLAS_TRANSPOSE transa,  CBLAS_TRANSPOSE transb,                                \
           int m, int n, int k, T alpha, const T* A, int lda,                               \
           const T* B, int ldb, T beta, T* C, int ldc)                                      \
{                                                                                \
    cblas_##name##gemm(order, transa, transb,                                               \
                       m, n, k, (const V*) &alpha, (const V*) A, lda,                       \
                       (const V*) B, ldb, (const V*) &beta, (V*) C, ldc);                   \
}

#define CUSP_CBLAS_COMPLEX_SYMM(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void symm(  CBLAS_ORDER order,                                                            \
              CBLAS_SIDE side,  CBLAS_UPLO uplo,                                            \
             int m, int n, T alpha, const T* A, int lda,                                    \
             const T* B, int ldb, T beta, T* C, int ldc)                                    \
{                                                                                \
    cblas_##name##symm(order, side, uplo, m, n,                                             \
                       (const V*) &alpha, (const V*) A, lda,                                \
                       (const V*) B, ldb, (const V*) &beta, (V*) C, ldc);                   \
}

#define CUSP_CBLAS_COMPLEX_SYRK(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void syrk(  CBLAS_ORDER order,                                                            \
              CBLAS_UPLO uplo,  CBLAS_TRANSPOSE trans,                                      \
             int n, int k, T alpha, const T* A, int lda,                                    \
             T beta, T* C, int ldc)                                                         \
{                                                                                \
    cblas_##name##syrk(order, uplo, trans, n, k,                                            \
                       (const V*) &alpha, (const V*) A, lda,                                \
                       (const V*) &beta, (V*) C, ldc);                                      \
}

#define CUSP_CBLAS_COMPLEX_SYR2K(T,V,name)                                                  \
  template<int dummy>                                                                       \
  void syr2k(  CBLAS_ORDER order,                                                           \
               CBLAS_UPLO uplo,  CBLAS_TRANSPOSE trans,                                     \
              int n, int k, T& alpha, const T* A, int lda,                                  \
              const T* B, int ldb, T& beta, T* C, int ldc)                                  \
{                                                                                \
    cblas_##name##syr2k(order, uplo, trans, n, k,                                           \
                        (const V*) &alpha, (const V*) A, lda,                               \
                        (const V*) B, ldb, (const V*) &beta, (V*) C, ldc);                  \
}

#define CUSP_CBLAS_COMPLEX_TRMM(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void trmm(  CBLAS_ORDER order,                                                            \
              CBLAS_SIDE side,  CBLAS_UPLO uplo,                                            \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int m, int n, T alpha, const T* A, int lda,                                    \
             T* B, int ldb)                                                                 \
{                                                                                \
    cblas_##name##trmm(order, side, uplo, trans, diag, m, n,                                \
                      (const V*) &alpha, (const V*) A, lda, (V*) B, ldb);                   \
}

#define CUSP_CBLAS_COMPLEX_TRSM(T,V,name)                                                   \
  template<int dummy>                                                                       \
  void trsm(  CBLAS_ORDER order,                                                            \
              CBLAS_SIDE side,  CBLAS_UPLO uplo,                                            \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int m, int n, T alpha, const T* A, int lda,                                    \
             T* B, int ldb)                                                                 \
{                                                                                \
    cblas_##name##trsm(order, side, uplo, trans, diag, m, n,                                \
                       (const V*) &alpha, (const V*) A, lda, (V*) B, ldb);                  \
}

namespace cusp
{
namespace system
{
namespace cpp
{
namespace detail
{
namespace cblas
{

// LEVEL 1
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_AMAX);
CUSP_CBLAS_COMPLEX_DEFS_2(CUSP_CBLAS_COMPLEX_ASUM);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_AXPY);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_COPY);
// CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_DOTC);
CUSP_CBLAS_COMPLEX_DEFS_2(CUSP_CBLAS_COMPLEX_NRM2);
CUSP_CBLAS_COMPLEX_DEFS_3(CUSP_CBLAS_COMPLEX_SCAL);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_SWAP);

// LEVEL 2
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_GEMV);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_GERC);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_HEMV);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_HER);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_TRMV);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_TRSV);

// LEVEL 3
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_GEMM);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_SYMM);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_SYRK);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_SYR2K);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_TRMM);
CUSP_CBLAS_COMPLEX_DEFS_1(CUSP_CBLAS_COMPLEX_TRSM);

} // end namespace cblas
} // end namespace detail
} // end namespace cpp
} // end namespace system
} // end namespace cusp

max_stars_repo_name: Raman-sh/cusplibrary


text:
/////////////////
//example25.5.c
/////////////////
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <math.h>

typedef struct{
  double g;//gravity
  double M;//mass of the car
  double m;//mass of the load
  double r;//length of the rope
  double C;//damping coefficient

  double F;//input force proportional to torque
  double T;//tension to the rope

  double x,x0;//x position of the load
  double y,y0;//y position of the load
  double X;//X position of the car
  double a;//angle (theta)

  double dx,dx0;
  double dy,dy0;
  double dX,dX0;
  double dr;
  double da,da0;

  double ddx;
  double ddy;
  double ddX;
  double ddr;
  double dda;
  double dXmax;
  double ddXmax;
  double Fmax;//Force

  double h;
  int nh;
  double _h;
  int _dim;
  double *_y;
  double *_y_err;
  double *_dydt_in;
  double *_dydt_out;
  const gsl_odeiv_step_type *_T;
  gsl_odeiv_step *_s;
  gsl_odeiv_system _sys;
  double _t;
} CRANE;

#define dim_crane 4
int cranefunc (double t, const double y[], double f[], void *params)
{
  CRANE *c = (CRANE *)params;
  //y[0]=crane.p, y[1]=crane.da, y[2]=crane.X, y[3]=crane.dX
  f[0] = y[1];
  f[1] = (-(2*c->dr+c->C)*y[1] -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  //  f[1] = (-2*c->dr -c->ddX*cos(y[0]) -c->g*sin(y[0]))/c->r;
  f[2] = y[3];
  //  f[3] = c->ddX;
  f[3] = (c->F+c->T*sin(y[0]))/c->M;
  return GSL_SUCCESS;
}
////////////////
#define square(x) ((x)*(x))
CRANE crane;
int initialize()
{

  crane.Fmax=20.0;//check//  crane.Fmax=30;//check
  //  crane.dXmax=1.0;
  crane.dXmax=1.0;
  crane.ddXmax=0.2;
  crane.g=9.8;
  crane.M=100;
  crane.m=20;//10;//10;//20
  crane.r=5;//5;
  crane.C=1;//5;
  crane.dr=crane.ddr=0;
  crane.x0=0;
  crane.y0=crane.r;
  crane.dX0=crane.dx0=crane.dy0=crane.da0=0;
#ifdef CRANESUB
  //  crane.Fmax=_crane_Fmax;//check//  crane.Fmax=30;//check
  crane.m=_crane_m;
  crane.C=_crane_C;
  crane.r=_crane_r;
  crane.h=AP_tS;
  crane.dXmax=_crane_dXmax;
  if(_AP_umax>0) AP_u_max=crane.Fmax=_AP_umax;
  else AP_u_max=crane.Fmax;
  AP_u_min=-AP_u_max;
  //  starttime=0;
  //  totaltime=100;//  totaltime=40;
  //  totaltime=50;//  totaltime=40;
  //  int kmax=totaltime/AP_tS;
  //  _rr=(double*)malloc(kmax*sizeof(double));
  rr=AP_r=_AP_r;//10
  //  int k;for(k=0;k<kmax;k++) _rr[k]=AP_r;//10;
  rr_kyoyou=_rr_kyoyou;
  //  p=(double *)malloc(buffsize*sizeof(double));
  C_MODE=11;
  //  iteration=_iteration;
  //  fprintf(stdout,"?????????iteration=%d\n",_iteration);
#else
  crane.h=0.01;
#endif
  crane.nh=10;
  crane.F=crane.a=crane.da=crane.dda=0;
  crane.T=crane.m*crane.g;
  crane._h=crane.h/crane.nh;//0.001
  crane._dim=4;
  if(crane._y==NULL){
    crane._y=(double*)malloc(crane._dim*sizeof(double));
    crane._y_err=(double*)malloc(crane._dim*sizeof(double));
    crane._dydt_in=(double*)malloc(crane._dim*sizeof(double));
    crane._dydt_out=(double*)malloc(crane._dim*sizeof(double));
  }
  crane._T= gsl_odeiv_step_rk4;
  crane._s= gsl_odeiv_step_alloc (crane._T, crane._dim);
  crane._sys= (gsl_odeiv_system){cranefunc, NULL, crane._dim, &crane};
  crane._t=crane._y[0]=crane._y[1]=crane._y[2]=crane._y[3]=0;
  int i;for(i=0;i<crane._dim;i++) crane._dydt_in[i]=0;
  GSL_ODEIV_FN_EVAL(&(crane._sys), crane._t, crane._y, crane._dydt_in);
  return(0);
}
#ifndef CRANESUB
char *fn="crane2io.dat";
FILE *fp;
#endif//#ifndef CRANESUB

double plant(double uu)
{
  //  crane.F=uu;
  if(uu>crane.Fmax) uu=crane.Fmax;
  else if(uu<-crane.Fmax) uu=-crane.Fmax;
  int n;
  for(n=0;n<crane.nh;n++){
    crane.F=uu;//
    if(crane.dX>=crane.dXmax && crane.F>0) crane.F=0;
    else if(crane.dX<=-crane.dXmax && crane.F<0) crane.F=0;
    int status = gsl_odeiv_step_apply (crane._s, crane._t, crane._h,
				       crane._y, crane._y_err,
				       crane._dydt_in,
				       crane._dydt_out,
				       &crane._sys);
    if (status != GSL_SUCCESS)   break;

    int i;for(i=0;i<dim_crane;i++) crane._dydt_in[i]=crane._dydt_out[i];
    crane._t+=crane._h;

    crane.a  =crane._y[0];
    crane.da =crane._y[1];
    crane.X  =crane._y[2];
    crane.dX =crane._y[3];

    crane.x  =crane.X+crane.r * sin(crane.a);
    crane.y  =crane.r*cos(crane.a);
    crane.dx=(crane.x-crane.x0)/crane._h;
    crane.dy=(crane.y-crane.y0)/crane._h;
    crane.ddx=(crane.dx-crane.dx0)/crane._h;
    crane.ddy=(crane.dy-crane.dy0)/crane._h;
    crane.T=crane.m*sqrt(square(crane.ddx)+square(crane.ddy-crane.g));
    crane.ddX=(crane.dX-crane.dX0)/crane._h;
    //    crane.F=crane.M*crane.ddX-crane.T*sin(crane.a);
    crane.dda=(crane.da-crane.da0)/crane._h;
    /*    crane.dda=(-2*crane.dr*crane.da-crane.g*sin(crane.a)-ddX[n]*cos(crane.a))/crane.r; */
    crane.x0=crane.x;
    crane.y0=crane.y;
    crane.dx0=crane.dx;
    crane.dy0=crane.dy;
    crane.da0=crane.da;
    crane.dX0=crane.dX;
  }
#ifndef CRANESUB
  fprintf(fp,"%.7e %.7e %.7e", crane._t,crane._y[0],crane._y[1]);//crane.a,crane.da
  fprintf(fp," %.7e %.7e %.7e" ,crane.X,crane.dX,crane.ddX);
  fprintf(fp," %.7e %.7e %.7e" ,crane.x,crane.dx,crane.ddx);
  fprintf(fp," %.7e %.7e %.7e" ,crane.y,crane.dy,crane.ddy);
  fprintf(fp," %.7e %.7e" ,crane.T,crane.F);
  fprintf(fp,"\n");
#endif
  return(crane.x);
  //  return(crane.X);
}

#ifndef CRANESUB
int main (void)
{
  //  CRANE crane;
  crane.g=9.8;
  crane.M=100;
  crane.m=20;
  crane.r=5;
  crane.dr=crane.ddr=0;

  ////////////////////////////////////////////////////
  /// method 2                                     ///
  /// input F                                     ///
  /// output p,x,y,X,...                          ///
  ////////////////////////////////////////////////////
  double h;//,hh=0.1;
  //result for different h
  /*gnuplot
   set style data lines;n=1;
   n=n+1;plot "crane1io.dat" using 1:n, "crane2io.dat" using 1:n;print "n=",n
   #strange at n=9(ddx),12(ddy),13(T)
  */
  ////////////////////////////////////////////////////
  ///   read F                                    ///
  ////////////////////////////////////////////////////
  char *fn1="crane1io.dat";
  FILE *fp1=fopen(fn1,"r");
  int n4;
#define buffsize 512
  char buff[buffsize];
  int n0;
  fgets(buff,buffsize,fp1);
  sscanf(buff,"#%lf%d%d%lf%lf%lf%lf",&h,&n0,&n4,&crane.M,&crane.m,&crane.r,&crane.C);
  crane.h=h;
  double t4=n4*h;
  double *F  =(double*)malloc(sizeof(double)*(n4+1));
  //  int M=n4+1;double *F  =(double*)malloc(sizeof(double)*M);
  //  double *ddX=(double*)malloc(sizeof(double)*M);

  int n;
  double t;
  //  double dx0=0,dy0=0,dX0=0,da0=0;
  //  double x0=0,y0=crane.r;

  fgets(buff,buffsize,fp1);
  for(n=0;n<n4;n++){//
    double _x1,_x2,_x3,_x4,_x5,_x6,_x7,_x8,_x9,_x10,_x11,_x12,_x13,_x14;
    fgets(buff,buffsize,fp1);
    sscanf(buff,"%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf",
	   &_x1,&_x2,&_x3,&_x4,&_x5,&_x6,&_x7,&_x8,&_x9,&_x10,&_x11,&_x12,&_x13,&_x14);
    F[n]=_x14;
  }
  ////////////////////////////////////////////////////
  ///   Solve by the Runge Kutta Method of GSL     ///
  ////////////////////////////////////////////////////
  initialize();
  fp=fopen(fn,"w");
  for(t=0;t<t4;t+=crane.h){
    n=t/crane.h;
    plant(F[n]);
  }
  fclose(fp);
  fprintf(stdout,"Results are stored in '%s'.\n",fn);

  gsl_odeiv_step_free (crane._s);
  return 0;
}
#endif


text:
/** @file */
#ifndef __CCL_UTILS_H_INCLUDED__
#define __CCL_UTILS_H_INCLUDED__

#include <gsl/gsl_spline.h>

#define CCL_MIN(a, b)  (((a) < (b)) ? (a) : (b))
#define CCL_MAX(a, b)  (((a) > (b)) ? (a) : (b))

CCL_BEGIN_DECLS

/**
 * Compute bin edges of N-1 linearly spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_linear_spacing(double xmin, double xmax, int N);

/**
 * Compute bin edges of N-1 logarithmically and then linearly spaced bins on the interval [xmin,xmax]
 * @param xminlog minimum value of spacing
 * @param xmin value when logarithmical ends and linear spacing begins
 * @param xmax maximum value of spacing
 * @param Nlin number of linear bins plus one (number of bin edges)
 * @param Nlog number of logarithmic bins plus one (number of bin edges)
 * @return x, bin edges in range [xminlog, xmax]
 */
double * ccl_linlog_spacing(double xminlog, double xmin, double xmax, int Nlin, int Nlog);

/**
 * Compute bin edges of N-1 logarithmically spaced bins on the interval [xmin,xmax]
 * @param xmin minimum value of spacing
 * @param xmax maximum value of spacing
 * @param N number of bins plus one (number of bin edges)
 * @return x, bin edges in range [xmin, xmax]
 */
double * ccl_log_spacing(double xmin, double xmax, int N);
//Returns array of N logarithmically-spaced values between xmin and xmax

double ccl_j_bessel(int l,double x);
//Spherical Bessel function of order l (adapted from CAMB)

CCL_END_DECLS

#endif

max_stars_repo_name: benediktdiemer/CCL

max_issues_repo_name: benediktdiemer/CCL

max_forks_repo_name: benediktdiemer/CCL


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

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram2d.h>

#include "find.c"

double
gsl_histogram2d_get (const gsl_histogram2d * h, const size_t i, const size_t j)
{
  const size_t nx = h->nx;
  const size_t ny = h->ny;

  if (i >= nx)
    {
      GSL_ERROR_VAL ("index i lies outside valid range of 0 .. nx - 1",
			GSL_EDOM, 0);
    }

  if (j >= ny)
    {
      GSL_ERROR_VAL ("index j lies outside valid range of 0 .. ny - 1",
			GSL_EDOM, 0);
    }

  return h->bin[i * ny + j];
}

int
gsl_histogram2d_get_xrange (const gsl_histogram2d * h, const size_t i,
			    double *xlower, double *xupper)
{
  const size_t nx = h->nx;

  if (i >= nx)
    {
      GSL_ERROR ("index i lies outside valid range of 0 .. nx - 1", GSL_EDOM);
    }

  *xlower = h->xrange[i];
  *xupper = h->xrange[i + 1];

  return GSL_SUCCESS;
}

int
gsl_histogram2d_get_yrange (const gsl_histogram2d * h, const size_t j,
			    double *ylower, double *yupper)
{
  const size_t ny = h->ny;

  if (j >= ny)
    {
      GSL_ERROR ("index j lies outside valid range of 0 .. ny - 1", GSL_EDOM);
    }

  *ylower = h->yrange[j];
  *yupper = h->yrange[j + 1];

  return GSL_SUCCESS;
}

int
gsl_histogram2d_find (const gsl_histogram2d * h,
		      const double x, const double y,
		      size_t * i, size_t * j)
{
  int status = find (h->nx, h->xrange, x, i);

  if (status)
    {
      GSL_ERROR ("x not found in range of h", GSL_EDOM);
    }

  status = find (h->ny, h->yrange, y, j);

  if (status)
    {
      GSL_ERROR ("y not found in range of h", GSL_EDOM);
    }

  return GSL_SUCCESS;
}


text:
/*
 * Copyright (c) 2016-2021 lymastee, All rights reserved.
 * Contact: lymastee@hotmail.com
 *
 * This file is part of the gslib project.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#pragma once

#ifndef scenegraph_979fc7df_3f68_4a38_9dea_5eed6d400c98_h
#define scenegraph_979fc7df_3f68_4a38_9dea_5eed6d400c98_h

#include <gslib/type.h>
#include <gslib/tree.h>
#include <ariel/config.h>

__ariel_begin__

class scene_node
{
public:
    //typedef gs::string string;
    typedef const string* sntype;
    typedef void* snptr;
    typedef snptr (*fnaccess)(int);
    typedef void (*fndestroy)(snptr);

protected:
    sntype          _tag;
    snptr           _ptr;
    string          _key;
    fndestroy       _del;

public:
    /* sntype declaration */
#define declare_sntype(snt) \
    static sntype snt;
#include <ariel\snt.h>
#undef declare_sntype

    static void register_sntypes();

public:
    scene_node() { reset(); }
    scene_node(const gchar* name)
    {
        reset();
        set_key(name);
    }
    scene_node(const gchar* name, sntype tag, snptr ptr, fndestroy del = 0)
    {
        set_key(name);
        bind(tag, ptr, del);
    }
    ~scene_node()
    {
        destroy();
        reset();
    }
    void reset()
    {
        _tag = sn_void;
        _ptr = 0;
        _del = 0;
    }
    void destroy() { if(_ptr && _del) _del(_ptr); }
    sntype get_type() const { return _tag; }
    snptr get_ptr() const { return _ptr; }
    const gchar* get_name() const { return _key.c_str(); }
    const string& get_key() const { return _key; }
    void set_destroy(fndestroy del) { _del = del; }
    void set_nodestroy() { set_destroy(0); }
    void set_key(const gchar* key) { _key.assign(key); }
    void bind(sntype tag, snptr ptr, fndestroy del = 0)
    {
        _tag = tag;
        _ptr = ptr;
        set_destroy(del);
    }
    void rebind(sntype tag, snptr ptr, fndestroy del = 0)
    {
        destroy();
        bind(tag, ptr, del);
    }
    scene_node& operator=(const scene_node& that)
    {
        _tag = that._tag;
        _ptr = that._ptr;
        _key = that._key;
        _del = that._del;
        const_cast<scene_node&>(that)._del = 0;
        return *this;
    }
};

struct scene_key
{
    const gchar*    _key;
    scene_key() { _key = 0; }
    scene_key(const gchar* k) { _key = k; }
    scene_key(const scene_node* n) { _key = n->get_name(); }
};

class scene_graph
{
public:
    struct snkey_hash
    {
        size_t operator()(const scene_key& key) const
        { return string_hash(key._key); }
    };
    struct snkey_equal
    {
        bool operator()(const scene_key& k1, const scene_key& k2) const
        { return string_hash(k1._key) == string_hash(k2._key); }
    };

    typedef tree<scene_node> sntree;
    typedef sntree::wrapper wrapper;
    typedef sntree::iterator iterator;
    typedef sntree::const_iterator const_iterator;
    typedef unordered_map<scene_key, scene_node*, snkey_hash, snkey_equal> sntable;
    typedef scene_node::sntype sntype;
    typedef scene_node::snptr snptr;
    typedef scene_node::fndestroy fndestroy;

public:
    template<class castop, int bias, class castp>
    static castop fbop_cast(castp ptr) { return ptr ? reinterpret_cast<castop>(((int)ptr) + bias) : 0; }
    template<class castop, class castp>
    static castop fbop_cast(castp ptr, int bias) { return ptr ? reinterpret_cast<castop>(((int)ptr) + bias) : 0; }
    static int wrapper_node_bias()
    {
        static wrapper inst;
        static const int bias = (int)&inst - ((int)inst.get_ptr());
        return bias;
    }
    static wrapper* node_to_wrapper(scene_node* ptr) { return fbop_cast<wrapper*>(ptr, -wrapper_node_bias()); }
    static const wrapper* node_to_wrapper(const scene_node* ptr) { return fbop_cast<const wrapper*>(ptr, -wrapper_node_bias()); }
    static scene_node* iter_to_node(iterator i) { return i.get_ptr(); }
    static const scene_node* iter_to_node(const_iterator i) { return i.get_ptr(); }
    static iterator node_to_iter(scene_node* ptr) { return iterator(node_to_wrapper(ptr)); }
    static const_iterator node_to_iter(const scene_node* ptr) { return const_iterator(node_to_wrapper(ptr)); }

public:
    scene_graph();
    virtual ~scene_graph() {}

public:
    scene_node* get_root_node() { return iter_to_node(_sntree.get_root()); }
    const scene_node* get_root_node() const { return iter_to_node(_sntree.get_root()); }
    bool is_node_attached(const scene_node* node) const { return node && _sntree.is_mine(node_to_iter(node)); }
    sntree* get_node_tree(const scene_node* node);
    const sntree* get_node_tree(const scene_node* node) const;
    void fix_node(scene_node* old, scene_node* new1);
    void map_node(scene_node* node, const gchar* name);
    void unmap_node(scene_node* node);

public:
    virtual scene_node* create_node(sntype tag = scene_node::sn_void, snptr ptr = 0, fndestroy del = 0);
    virtual scene_node* create_node(const gchar* name, sntype tag = scene_node::sn_void, snptr ptr = 0, fndestroy del = 0);
    virtual scene_node* attach_node(scene_node* parent, scene_node* node);
    virtual scene_node* detach_node(scene_node* node);
    virtual void destroy_node(scene_node* node);
    virtual void clear_miscs();

protected:
    sntree          _sntree;
    sntree          _miscs;
    sntable         _sntable;
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
/* rng/r250.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 a shift-register random number generator. The sequence is

   x_n = x_{n-103} ^ x_{n-250}        ("^" means XOR)

   defined on 32-bit words.

   The first 250 elements x_1 .. x_250 are first initialized as x_n =
   s_n, where s_n = (69069*s_{n-1}) mod 2^32 and s_0=s is the
   user-supplied seed. To ensure that the sequence does not lie on a
   subspace we force 32 of the entries to be linearly independent.  We
   take the 32 elements x[3], x[10], x[17], x[24], ..., 213 and apply
   the following operations,

   x[3]   &= 11111111111111111111111111111111
   x[3]   |= 10000000000000000000000000000000
   x[10]  &= 01111111111111111111111111111111
   x[10]  |= 01000000000000000000000000000000
   x[17]  &= 00111111111111111111111111111111
   x[17]  |= 00100000000000000000000000000000
   ....      ...
   x[206] &= 00000000000000000000000000000111
   x[206] |= 00000000000000000000000000000100
   x[213] &= 00000000000000000000000000000011
   x[213] |= 00000000000000000000000000000010
   x[220] &= 00000000000000000000000000000001
   x[220] |= 00000000000000000000000000000001

   i.e. if we consider the bits of the 32 elements as forming a 32x32
   array then we are setting the diagonal bits of the array to one and
   masking the lower triangle below the diagonal to zero.

   With this initialization procedure the theoretical value of
   x_{10001} is 1100653588 for s = 1 (Actually I got this by running
   the original code). The subscript 10001 means (1) seed the
   generator with s = 1 and then do 10000 actual iterations.

   The period of this generator is about 2^250.

   The algorithm works for any number of bits. It is implemented here
   for 32 bits.

   From: S. Kirkpatrick and E. Stoll, "A very fast shift-register
   sequence random number generator", Journal of Computational Physics,
   40, 517-526 (1981). */

static inline unsigned long int r250_get (void *vstate);
static double r250_get_double (void *vstate);
static void r250_set (void *state, unsigned long int s);

typedef struct
  {
    int i;
    unsigned long x[250];
  }
r250_state_t;

static inline unsigned long int
r250_get (void *vstate)
{
  r250_state_t *state = (r250_state_t *) vstate;
  unsigned long int k;
  int j;

  int i = state->i;

  if (i >= 147)
    {
      j = i - 147;
    }
  else
    {
      j = i + 103;
    }

  k = state->x[i] ^ state->x[j];
  state->x[i] = k;

  if (i >= 249)
    {
      state->i = 0;
    }
  else
    {
      state->i = i + 1;
    }

  return k;
}

static double
r250_get_double (void *vstate)
{
  return r250_get (vstate) /  4294967296.0 ;
}

static void
r250_set (void *vstate, unsigned long int s)
{
  r250_state_t *state = (r250_state_t *) vstate;

  int i;

  if (s == 0)
    s = 1;	/* default seed is 1 */

  state->i = 0;

#define LCG(n) ((69069 * n) & 0xffffffffUL)

  for (i = 0; i < 250; i++)	/* Fill the buffer  */
    {
      s = LCG (s);
      state->x[i] = s;
    }

  {
    /* Masks for turning on the diagonal bit and turning off the
       leftmost bits */

    unsigned long int msb = 0x80000000UL;
    unsigned long int mask = 0xffffffffUL;

    for (i = 0; i < 32; i++)
      {
	int k = 7 * i + 3;	/* Select a word to operate on        */
	state->x[k] &= mask;	/* Turn off bits left of the diagonal */
	state->x[k] |= msb;	/* Turn on the diagonal bit           */
	mask >>= 1;
	msb >>= 1;
      }
  }

  return;
}

static const gsl_rng_type r250_type =
{"r250",			/* name */
 0xffffffffUL,			/* RAND_MAX */
 0,			        /* RAND_MIN */
 sizeof (r250_state_t),
 &r250_set,
 &r250_get,
 &r250_get_double};

const gsl_rng_type *gsl_rng_r250 = &r250_type;


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

#include <config.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

/* Factorise a general M x N matrix A into
 *
 *   A = Q L
 *
 * where Q is orthogonal (M x M) and L is lower triangular (M x N).
 *
 * Q is stored as a packed set of Householder transformations in the
 * strict upper triangular part of the input matrix.
 *
 * L is stored in the diagonal and lower triangle of the input matrix.
 *
 * The full matrix for Q can be obtained as the product
 *
 *       Q = Q_k .. Q_2 Q_1
 *
 * where k = MIN(M,N) and
 *
 *       Q_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [A(1,N-k+i), A(2,N-k+i), ... , A(M-k+i,N-k+i), 1, 0, ..., 0]
 *
 * This storage scheme is the same as in LAPACK.  */

/*
gsl_linalg_QL_decomp()
  Perform QL decomposition of a matrix A

Inputs: A   - M-by-N matrix
        tau - (output) Householder coefficients, length N

Notes:
1) The K = MIN(M, N) Householder scalars are stored in tau(N-K+1:N)
on output; the rest of tau is used as temporary workspace
*/

int
gsl_linalg_QL_decomp (gsl_matrix * A, gsl_vector * tau)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (tau->size != N)
    {
      GSL_ERROR ("size of tau must be N", GSL_EBADLEN);
    }
  else
    {
      const size_t K = GSL_MIN(M, N);
      size_t i;

      for (i = 0; i < K; i++)
        {
          /* compute the Householder transformation to annihilate the (N-K+i)-th
             column of the matrix */

          gsl_vector_view c = gsl_matrix_subcolumn (A, N - i - 1, 0, M - i);
          double * alpha = gsl_matrix_ptr(A, M - i - 1, N - i - 1);
          double tau_j = gsl_linalg_householder_transform2 (alpha, &(c.vector));

          /* apply the transformation to A(1:M-i,1:N-i-2) from the left */
          if (i + 1 < N)
            {
              gsl_vector_view work = gsl_vector_subvector(tau, 0, N - i - 1);
              gsl_matrix_view m = gsl_matrix_submatrix (A, 0, 0, M - i, N - i - 1);
              double tmp = *alpha;

              *alpha = 1.0;
              gsl_linalg_householder_left (tau_j, &(c.vector), &(m.matrix), &work.vector);
              *alpha = tmp;
            }

          gsl_vector_set (tau, N - i - 1, tau_j);
        }

      return GSL_SUCCESS;
    }
}

/* form the orthogonal matrix Q from the packed QL matrix */
int
gsl_linalg_QL_unpack (const gsl_matrix * QL, const gsl_vector * tau, gsl_matrix * Q, gsl_matrix * L)
{
  const size_t M = QL->size1;
  const size_t N = QL->size2;

  if (Q->size1 != M || Q->size2 != M)
    {
      GSL_ERROR ("Q matrix must be M x M", GSL_ENOTSQR);
    }
  else if (L->size1 != M || L->size2 != N)
    {
      GSL_ERROR ("L matrix must be M x N", GSL_ENOTSQR);
    }
  else if (tau->size != N)
    {
      GSL_ERROR ("size of tau must be N", GSL_EBADLEN);
    }
  else
    {
      const size_t K = GSL_MIN(M, N);
      size_t i;

      /* initialize Q to the identity */
      gsl_matrix_set_identity (Q);

      for (i = 0; i < K; ++i)
        {
          gsl_vector_const_view h = gsl_matrix_const_subcolumn (QL, N - K + i, 0, M - K + i + 1);
          gsl_matrix_view m = gsl_matrix_submatrix (Q, 0, 0, M - K + i + 1, M - K + i + 1);
          gsl_vector_view work = gsl_matrix_subcolumn(L, 0, 0, M - K + i + 1);
          double ti = gsl_vector_get (tau, N - K + i);
          double * ptr = gsl_matrix_ptr((gsl_matrix *) QL, M - K + i, N - K + i);
          double tmp = *ptr;

          *ptr = 1.0;
          gsl_linalg_householder_left (ti, &h.vector, &m.matrix, &work.vector);
          *ptr = tmp;
        }

      /* form the left triangular matrix L from a packed QL matrix */
      gsl_matrix_set_zero(L);

      if (M >= N)
        {
          gsl_matrix_const_view src = gsl_matrix_const_submatrix(QL, M - N, 0, N, N);
          gsl_matrix_view dest = gsl_matrix_submatrix(L, M - N, 0, N, N);
          gsl_matrix_tricpy(CblasLower, CblasNonUnit, &dest.matrix, &src.matrix);
        }
      else
        {
          gsl_matrix_const_view src1 = gsl_matrix_const_submatrix(QL, 0, 0, M, N - M);
          gsl_matrix_view dest1 = gsl_matrix_submatrix(L, 0, 0, M, N - M);

          gsl_matrix_const_view src2 = gsl_matrix_const_submatrix(QL, 0, N - M, M, M);
          gsl_matrix_view dest2 = gsl_matrix_submatrix(L, 0, N - M, M, M);

          gsl_matrix_memcpy(&dest1.matrix, &src1.matrix);
          gsl_matrix_tricpy(CblasLower, CblasNonUnit, &dest2.matrix, &src2.matrix);
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/*==================BEGIN ASF AUTO-GENERATED DOCUMENTATION==================*/
/*
ABOUT EDITING THIS DOCUMENTATION:
If you wish to edit the documentation for this program, you need to
change the following defines. For the short ones (like
ASF_NAME_STRING) this is no big deal. However, for some of the longer
ones, such as ASF_COPYRIGHT_STRING, it can be a daunting task to get
all the newlines in correctly, etc. In order to help you with this
task, there is a tool, edit_man_header. The tool *only* works with
this portion of the code, so fear not. It will scan in defines of the
format #define ASF_<something>_STRING between the two auto-generated
documentation markers, format them for a text editor, run that editor,
allow you to edit the text in a clean manner, and then automatically
generate these defines, formatted appropriately. The only warning is
that any text between those two markers and not part of one of those
defines will not be preserved, and that all of this auto-generated
code will be at the top of the source file. Save yourself the time and
trouble, and use edit_man_header. :)
*/

#define ASF_NAME_STRING \
"asf_geocode"

#define ASF_USAGE_STRING \
"   "ASF_NAME_STRING" -p <projection name> <<projection parameters>>\n"\
"             [-force] [-resample-method <method>] [-height <height>]\n"\
"             [-datum <datum>] [-pixel-size <pixel size>] [-band <band_id | all>]\n"\
"             [-log <file>] [-write-proj-file <file>] [-read-proj-file <file>]\n"\
"             [-save-mapping] [-background <value>] [-quiet] [-license]\n"\
"             [-version] [-help]\n"\
"             <in_base_name> <out_base_name>\n"\
"\n"\
"   Use the -help option for more projection parameter controls.\n"

#define ASF_DESCRIPTION_STRING \
"     This program takes a map projected or an unprojected (ground\n"\
"     range) image in the ASF internal format and geocodes it,\n"\
"     i.e. swizzles it around into one of the standard projections used\n"\
"     for maps (universal transverse mercator, polar stereo, etc).  The\n"\
"     output is a new image in ASF internal format.\n"

#define ASF_INPUT_STRING \
"     Most of the \"options\" are actually required.  The specification\n"\
"     of a certain projection type implies that all the parameters\n"\
"     required to fully specify a projection of that type be included.\n"\
"\n"\
"     This must be an ASF internal format image base name.\n"

#define ASF_OUTPUT_STRING \
"     The base name of the geocoded image to produce.\n"

#define ASF_OPTIONS_STRING \
"%s"\
"\n"\
"     -save-mapping\n"\
"          Creates two additional files during geocoding.  One will contain\n"\
"          the line number from which the pixel was obtained from the\n"\
"          original file, the other the sample numbers.  Together, these\n"\
"          define the mapping of pixels performed by the geocoding.\n"\
"\n"\
"     -log <log file>\n"\
"          Output will be written to a specified log file.\n"\
"\n"\
"     -quiet\n"\
"          Supresses all non-essential output.\n"\
"\n"\
"     -license\n"\
"          Print copyright and license for this software then exit.\n"\
"\n"\
"     -version\n"\
"          Print version and copyright then exit.\n"\
"\n"\
"     -help\n"\
"          Print a help page and exit.\n"

#define ASF_EXAMPLES_STRING \
"     To map project an image with centerpoint at -147 degrees\n"\
"     longitude and average height 466 meters into universal transverse\n"\
"     mercator projection, with one pixel 50 meters on a side:\n"\
"\n"\
"     "ASF_NAME_STRING" -p utm --central-meridian -147.0 --height 466\n"\
"                 input_image output_image\n"\
"\n"\
"     To geocode one band within an image file, you specify the selected band\n"\
"     with the -band option, and the selected band MUST be one of which appears\n"\
"     in the list of available bands as noted in the 'bands' item found in the\n"\
"     'general' (first) block in the metadata file.  For example, if 'bands'\n"\
"     contains \"01,02,03,04\", then you could specify a band_id, e.g. \"-band 02\"\n"\
"     etc on the command line.  The same applies to band lists, such\n"\
"     as \"HH,HV,VH,VV\" or just \"03\" etc.\n"\
"\n"\
"     "ASF_NAME_STRING" -p utm -band HV file outfile_HV\n"

#define ASF_LIMITATIONS_STRING \
"     May fail badly if bad projection parameters are supplied for the\n"\
"     area in the image.\n"

#define ASF_SEE_ALSO_STRING \
"     asf_import, asf_export\n"

/*===================END ASF AUTO-GENERATED DOCUMENTATION===================*/
#include <asf_contact.h>
#include <asf_license.h>

// Standard libraries.
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// Libraries from packages outside ASF.
#include <glib.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_statistics_double.h>

// Libraries developed at ASF.
#include <asf.h>
#include <asf_nan.h>
#include <asf_meta.h>
#include <asf_raster.h>
#include "float_image.h"
#include <libasf_proj.h>
#include <spheroids.h>
#include <asf_contact.h>

// Headers used by this program.
#include "asf_geocode.h"

// Print minimalistic usage info & exit
static void print_usage(void)
{
  asfPrintStatus("\n"
      "Usage:\n"
      ASF_USAGE_STRING
      "\n");
  exit(EXIT_FAILURE);
}

// Print the help info & exit
static void print_help(void)
{
  asfPrintStatus(
      "\n"
      "Tool name:\n   " ASF_NAME_STRING "\n\n"
      "Usage:\n" ASF_USAGE_STRING "\n"
      "Description:\n" ASF_DESCRIPTION_STRING "\n"
      "Input:\n" ASF_INPUT_STRING "\n"
      "Output:\n"ASF_OUTPUT_STRING "\n"
      "Options:\n" ASF_OPTIONS_STRING "\n"
      "Examples:\n" ASF_EXAMPLES_STRING "\n"
      "Limitations:\n" ASF_LIMITATIONS_STRING "\n"
      "See also:\n" ASF_SEE_ALSO_STRING "\n"
      "Contact:\n" ASF_CONTACT_STRING "\n"
      "Version:\n   " SVN_REV " (part of " TOOL_SUITE_NAME " " MAPREADY_VERSION_STRING ")\n\n",
        geocode_projection_options_help());
  exit(EXIT_FAILURE);
}

// Main routine.
int
main (int argc, char **argv)
{
  int force_flag = FALSE;
  int debug_dump = FALSE;
  char band_id[256]="";
  char *in_base_name, *out_base_name;

  in_base_name = (char *) MALLOC(sizeof(char)*255);
  out_base_name = (char *) MALLOC(sizeof(char)*255);

  // Get the projection parameters from the command line.
  projection_type_t projection_type;
  // Terrain height to assume.  Defaults to 0 (no height correction)
  double average_height = 0.0;
  // Pixel size to use for output image, in projection coordinate
  // units.  This variable corresponds to a "private"
  // (i.e. undocumented, so users don't fiddle with it) option.
  // NOTE: Setting the pixel size to a negative number results in the
  // code to pick a pixel size from the metadata
  double pixel_size = -1;
  // Datum to use in the target projection
  datum_type_t datum;
  spheroid_type_t spheroid;
  // Method to use to resample images.
  resample_method_t resample_method;
  // Value to put in the region outside the image
  double background_val = 0.0;
  // Should we save the mapping files?
  int save_map_flag;

  if (detect_flag_options(argc, argv, "-help", "--help", "-h", NULL)) {
    print_help();
  }

  // Detect & Process logging arguments
  if ((logflag = detect_string_options(argc, argv, logFile,
                      "-log", "--log", NULL))) {
      fLog = fopen (logFile, "a");
      if ( fLog == NULL ) {
    // Couldn't open the log file, so just don't do logging.
    logflag = FALSE;
      }
  }
  quietflag = detect_flag_options(argc, argv, "-quiet", "--quiet", NULL);
  save_map_flag = extract_flag_options(&argc, &argv, "-save-mapping", "--save_mapping", NULL);

  handle_license_and_version_args(argc, argv, ASF_NAME_STRING);

  asfSplashScreen(argc, argv);

  project_parameters_t *pp
    = get_geocode_options (&argc, &argv, &projection_type, &average_height,
			   &pixel_size, &datum, &spheroid, &resample_method,
			   &force_flag, band_id);

  if (!pp && projection_type != LAT_LONG_PSEUDO_PROJECTION) {
      print_usage();
  }

  // The argument at which the filenames start
  int arg_num = 1;

  if (detect_flag_options(argc, argv, "-debug", NULL)) {
    debug_dump=TRUE;
    ++arg_num;
  }

  if (!extract_double_options(&argc, &argv, &background_val, "--background",
			      "-background", NULL)) {
    strcpy(in_base_name, argv[arg_num]);
    meta_parameters *meta = meta_read(in_base_name);
    background_val = meta->general->no_data;
    asfPrintStatus("Extracted no data value (%f) out of the metadata.\n",
		   meta->general->no_data);
    meta_free(meta);
  }
  if (ISNAN(background_val)) background_val = DEFAULT_NO_DATA_VALUE;

  // Get non-option command line arguments.
  if ( argc != 3 && !debug_dump ) {
    int ii;
    int bad_arg = FALSE;

    for (ii = 0; ii < argc; ++ii) {
      if (argv[ii][0] == '-') {
    bad_arg = TRUE;
    asfPrintStatus("Unrecognized argument: %s\n", argv[ii]);
      }
    }

    if (!bad_arg)
      fprintf (stderr, "Wrong number of arguments\n");

    print_usage ();
  }

  strcpy (in_base_name, argv[arg_num]);
  strcpy (out_base_name, argv[arg_num + 1]);

  // Strip .img and also check to make sure input and output
  // base names are not the same
  char *ext = findExt(in_base_name);
  if (ext && strncmp(ext, ".img", 4) == 0) *ext = '\0';
  ext = findExt(out_base_name);
  if (ext && strncmp(ext, ".img", 4) == 0) *ext = '\0';
  char msg[512];
  sprintf(msg,"Input and output basenames cannot be the same:\n"
      "    Input base name : %s\n"
          "    Output base name: %s\n",
      in_base_name, out_base_name);
  if (strcmp(in_base_name, out_base_name) == 0)
    asfPrintError(msg);

  // Call library function that does the actual work
  asf_geocode(pp, projection_type, force_flag, resample_method, average_height,
          datum, pixel_size, band_id, in_base_name, out_base_name,
              (float)background_val, save_map_flag);

  // Close Log, if needed
  if (logflag)
    FCLOSE (fLog);

  exit (EXIT_SUCCESS);
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


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

#ifndef __GSL_MULTIFIT_NLIN_H__
#define __GSL_MULTIFIT_NLIN_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_permutation.h>

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

__BEGIN_DECLS

int gsl_multifit_gradient (const gsl_matrix * J, const gsl_vector * f,
                           gsl_vector * g);

int gsl_multifit_covar (const gsl_matrix * J, const double epsrel, gsl_matrix * covar);
int gsl_multifit_covar_QRPT (gsl_matrix * r, gsl_permutation * perm,
                             const double epsrel, gsl_matrix * covar);

/* Definition of vector-valued functions with parameters based on gsl_vector */

struct gsl_multifit_function_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  size_t n;   /* number of functions */
  size_t p;   /* number of independent variables */
  void * params;
};

typedef struct gsl_multifit_function_struct gsl_multifit_function ;

#define GSL_MULTIFIT_FN_EVAL(F,x,y) (*((F)->f))(x,(F)->params,(y))

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n, size_t p);
    int (*set) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    int (*iterate) (void *state, gsl_multifit_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
    void (*free) (void *state);
  }
gsl_multifit_fsolver_type;

typedef struct
  {
    const gsl_multifit_fsolver_type * type;
    gsl_multifit_function * function ;
    gsl_vector * x ;
    gsl_vector * f ;
    gsl_vector * dx ;
    void *state;
  }
gsl_multifit_fsolver;

gsl_multifit_fsolver *
gsl_multifit_fsolver_alloc (const gsl_multifit_fsolver_type * T,
                            size_t n, size_t p);

void gsl_multifit_fsolver_free (gsl_multifit_fsolver * s);

int gsl_multifit_fsolver_set (gsl_multifit_fsolver * s,
                              gsl_multifit_function * f,
                              const gsl_vector * x);

int gsl_multifit_fsolver_iterate (gsl_multifit_fsolver * s);

int gsl_multifit_fsolver_driver (gsl_multifit_fsolver * s,
                                 const size_t maxiter,
                                 const double epsabs, const double epsrel);

const char * gsl_multifit_fsolver_name (const gsl_multifit_fsolver * s);
gsl_vector * gsl_multifit_fsolver_position (const gsl_multifit_fsolver * s);

/* Definition of vector-valued functions and gradient with parameters
   based on gsl_vector */

struct gsl_multifit_function_fdf_struct
{
  int (* f) (const gsl_vector * x, void * params, gsl_vector * f);
  int (* df) (const gsl_vector * x, void * params, gsl_matrix * df);
  int (* fdf) (const gsl_vector * x, void * params, gsl_vector * f, gsl_matrix *df);
  size_t n;       /* number of functions */
  size_t p;       /* number of independent variables */
  void * params;  /* user parameters */
  size_t nevalf;  /* number of function evaluations */
  size_t nevaldf; /* number of Jacobian evaluations */
};

typedef struct gsl_multifit_function_fdf_struct gsl_multifit_function_fdf ;

typedef struct
  {
    const char *name;
    size_t size;
    int (*alloc) (void *state, size_t n, size_t p);
    int (*set) (void *state, const gsl_vector * wts,
                gsl_multifit_function_fdf * fdf, gsl_vector * x,
                gsl_vector * f, gsl_vector * dx);
    int (*iterate) (void *state, const gsl_vector * wts,
                    gsl_multifit_function_fdf * fdf, gsl_vector * x,
                    gsl_vector * f, gsl_vector * dx);
    int (*gradient) (void *state, gsl_vector * g);
    int (*jac) (void *state, gsl_matrix * J);
    void (*free) (void *state);
  }
gsl_multifit_fdfsolver_type;

typedef struct
  {
    const gsl_multifit_fdfsolver_type * type;
    gsl_multifit_function_fdf * fdf ;
    gsl_vector * x;        /* parameter values x */
    gsl_vector * f;        /* residual vector f(x) */
    gsl_vector * dx;       /* step dx */
    gsl_vector * g;        /* gradient J^T f */
    gsl_vector * sqrt_wts; /* sqrt(wts) */
    size_t niter;          /* number of iterations performed */
    void *state;
  }
gsl_multifit_fdfsolver;

gsl_multifit_fdfsolver *
gsl_multifit_fdfsolver_alloc (const gsl_multifit_fdfsolver_type * T,
                              size_t n, size_t p);

int
gsl_multifit_fdfsolver_set (gsl_multifit_fdfsolver * s,
                            gsl_multifit_function_fdf * fdf,
                            const gsl_vector * x);
int gsl_multifit_fdfsolver_wset (gsl_multifit_fdfsolver * s,
                                 gsl_multifit_function_fdf * f,
                                 const gsl_vector * x,
                                 const gsl_vector * wts);

int
gsl_multifit_fdfsolver_iterate (gsl_multifit_fdfsolver * s);

int gsl_multifit_fdfsolver_driver (gsl_multifit_fdfsolver * s,
                                   const size_t maxiter,
                                   const double xtol,
                                   const double gtol,
                                   const double ftol,
                                   int *info);

int gsl_multifit_fdfsolver_jac (gsl_multifit_fdfsolver * s,
                                gsl_matrix * J);

void
gsl_multifit_fdfsolver_free (gsl_multifit_fdfsolver * s);

const char * gsl_multifit_fdfsolver_name (const gsl_multifit_fdfsolver * s);
gsl_vector * gsl_multifit_fdfsolver_position (const gsl_multifit_fdfsolver * s);
gsl_vector * gsl_multifit_fdfsolver_residual (const gsl_multifit_fdfsolver * s);
size_t gsl_multifit_fdfsolver_niter (const gsl_multifit_fdfsolver * s);
int gsl_multifit_eval_wf(gsl_multifit_function_fdf *fdf,
                         const gsl_vector *x, const gsl_vector *wts,
                         gsl_vector *y);
int gsl_multifit_eval_wdf(gsl_multifit_function_fdf *fdf,
                          const gsl_vector *x, const gsl_vector *wts,
                          gsl_matrix *dy);

int gsl_multifit_fdfsolver_test (const gsl_multifit_fdfsolver * s,
                                 const double xtol,
                                 const double gtol,
                                 const double ftol, int *info);
int gsl_multifit_test_delta (const gsl_vector * dx, const gsl_vector * x,
                             double epsabs, double epsrel);

int gsl_multifit_test_gradient (const gsl_vector * g, double epsabs);

int gsl_multifit_fdfsolver_dif_df(const gsl_vector *x,
                                  const gsl_vector *wts,
                                  gsl_multifit_function_fdf *fdf,
                                  const gsl_vector *f, gsl_matrix *J);
int gsl_multifit_fdfsolver_dif_fdf(const gsl_vector *x, gsl_multifit_function_fdf *fdf,
                                   gsl_vector *f, gsl_matrix *J);

typedef struct
{
  size_t n;                         /* number of (original) residuals */
  size_t p;                         /* number of model parameters */
  double lambda;                    /* damping parameter */
  const gsl_vector *L_diag;         /* diagonal damping matrix or NULL */
  const gsl_matrix *L;              /* general damping matrix or NULL */
  gsl_vector *f;                    /* function values for finite diff J */
  gsl_vector *wts;                  /* weight vector for augmented system */
  gsl_multifit_fdfsolver * s;
  gsl_multifit_function_fdf *fdf;   /* user defined fdf */
  gsl_multifit_function_fdf fdftik; /* Tikhonov modified fdf */
} gsl_multifit_fdfridge;

gsl_multifit_fdfridge *
gsl_multifit_fdfridge_alloc (const gsl_multifit_fdfsolver_type * T,
                             const size_t n, const size_t p);
void gsl_multifit_fdfridge_free(gsl_multifit_fdfridge *work);
const char *gsl_multifit_fdfridge_name(const gsl_multifit_fdfridge * w);
gsl_vector *gsl_multifit_fdfridge_position (const gsl_multifit_fdfridge * w);
gsl_vector *gsl_multifit_fdfridge_residual (const gsl_multifit_fdfridge * w);
size_t gsl_multifit_fdfridge_niter (const gsl_multifit_fdfridge * w);
int gsl_multifit_fdfridge_set (gsl_multifit_fdfridge * w,
                               gsl_multifit_function_fdf * f,
                               const gsl_vector * x,
                               const double lambda);
int gsl_multifit_fdfridge_wset (gsl_multifit_fdfridge * w,
                                gsl_multifit_function_fdf * f,
                                const gsl_vector * x,
                                const double lambda,
                                const gsl_vector * wts);
int gsl_multifit_fdfridge_set2 (gsl_multifit_fdfridge * w,
                                gsl_multifit_function_fdf * f,
                                const gsl_vector * x,
                                const gsl_vector * lambda);
int gsl_multifit_fdfridge_wset2 (gsl_multifit_fdfridge * w,
                                 gsl_multifit_function_fdf * f,
                                 const gsl_vector * x,
                                 const gsl_vector * lambda,
                                 const gsl_vector * wts);
int gsl_multifit_fdfridge_set3 (gsl_multifit_fdfridge * w,
                                gsl_multifit_function_fdf * f,
                                const gsl_vector * x,
                                const gsl_matrix * L);
int gsl_multifit_fdfridge_wset3 (gsl_multifit_fdfridge * w,
                                 gsl_multifit_function_fdf * f,
                                 const gsl_vector * x,
                                 const gsl_matrix * L,
                                 const gsl_vector * wts);
int gsl_multifit_fdfridge_iterate (gsl_multifit_fdfridge * w);
int gsl_multifit_fdfridge_driver (gsl_multifit_fdfridge * w,
                                  const size_t maxiter,
                                  const double xtol,
                                  const double gtol,
                                  const double ftol,
                                  int *info);

/* extern const gsl_multifit_fsolver_type * gsl_multifit_fsolver_gradient; */

GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmsder;
GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmder;
GSL_VAR const gsl_multifit_fdfsolver_type * gsl_multifit_fdfsolver_lmniel;

__END_DECLS

#endif /* __GSL_MULTIFIT_NLIN_H__ */

max_stars_repo_name: shi-bash-cmd/qtTest

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_linalg.h>

#include "../include/run_svd.h"

int main(void) {
  FILE *fp = fopen("./data/my_data.csv", "r");
  if (fp == NULL) {
    perror("Unable to open file!");

    exit(1);
  }
  int m = 4;
  int n = 5;

  char chunk[128];

  const char s[2] = ",";
  int i = 0;
  int j = 0;
  gsl_matrix *mat = gsl_matrix_alloc(m, n);
  while (fgets(chunk, sizeof(chunk), fp) != NULL) {
    char *token;

    token = strtok(chunk, s);
    while (token != NULL) {
       double x = atof(token);
       gsl_matrix_set(mat, i, j, x);

       j = (j + 1) % n;
       if (j == 0) {
         i = (i + 1) % m;
       }

       token = strtok(NULL, s);
    }
  }
  printf("a_matrix\n");
  pretty_print(mat);

  int result = run_svd(mat);

  fclose(fp);

  return result;
}

max_stars_repo_name: paul-reiners/getting-smaller

max_issues_repo_name: paul-reiners/getting-smaller

max_forks_repo_name: paul-reiners/getting-smaller


text:
#ifndef S3D_CV_DISPARITY_DISPARITY_ANALYZER_STAN_H
#define S3D_CV_DISPARITY_DISPARITY_ANALYZER_STAN_H

#include <s3d/disparity/disparity_analyzer.h>

#include "s3d/cv/features/match_finder_surf.h"
#include "s3d/multiview/sampson_distance_function.h"
#include "s3d/multiview/stan_fundamental_matrix_solver.h"
#include "s3d/multiview/stan_results.h"
#include "s3d/robust/ransac.h"
#include "s3d/robust/lmeds.h"
#include "s3d/utilities/math.h"

#include "s3d/utilities/eigen.h"

#include <opencv2/core/types.hpp>

#include <gsl/gsl>

namespace s3d {
struct StanAlignment;
class MatchFinderCV;
class MatchFinderSurf;

class DisparityAnalyzerSTAN : public DisparityAnalyzer {
 public:
  // can alternatively use Lmeds
  using RansacAlgorithmSTAN =
      s3d::robust::Lmeds<s3d::StanFundamentalMatrixSolver, s3d::SampsonDistanceFunction>;

  struct Results {
    explicit Results();

    void updateParameters(double minDisparity,
                          double maxDisparity,
                          float widthRatio,
                          RansacAlgorithmSTAN::ModelType model);

    void updatePoints(const std::vector<Eigen::Vector2d>& bestPtsLeft,
                      const std::vector<Eigen::Vector2d>& bestPtsRight,
                      std::vector<double> disparities,
                      float widthRatio,
                      float resizeRatio);

    StanResults stan;
    double minDisparityPercent{};
    double maxDisparityPercent{};

    std::vector<float> disparitiesPercent{};
  };

  DisparityAnalyzerSTAN();

  gsl::owner<DisparityAnalyzerSTAN*> clone() const override;

  void setMatchFinder(std::unique_ptr<MatchFinderCV> matchFinder);
  void setMaxNumberOfFeatures(int maxNumberOfFeatures);

  bool analyze(const cv::Mat& left, const cv::Mat& right);
  bool analyze(const Image<uint8_t>& left, const Image<uint8_t>& right) override;
  const std::vector<float>& getDisparitiesPercent() const override;
  const std::vector<Eigen::Vector2f>& getFeaturePointsLeft() const override;
  const std::vector<Eigen::Vector2f>& getFeaturePointsRight() const override;
  void setMinimumNumberOfInliers(int minNbInliers);

  // outputs
  Results results;

 private:
  MatchFinder::Matches findMatches(const cv::Mat& left, const cv::Mat& right);
  bool enoughMatches(int nbOfMatches);
  RansacAlgorithmSTAN createRansac(Size imageSize);

  std::unique_ptr<MatchFinderCV> matchFinder_{std::make_unique<MatchFinderCV>()};
  int minNbInliers_{4};
};
}  // namespace s3d

#endif  // S3D_CV_DISPARITY_DISPARITY_ANALYZER_STAN_H