text:
#pragma once
#include "icommand.h"
#include "optioninfo.h"
#include "configaccess.h"
#include "config.h"
#include "format.h"
#include "streamreader.h"
#include <gsl/gsl>
#include <cmdlime/errors.h>
#include <cmdlime/customnames.h>
#include <cmdlime/detail/configaccess.h>
#include <cmdlime/detail/flag.h>
#include <sstream>
#include <functional>
#include <memory>

namespace cmdlime::detail{

template <typename TConfig>
class Command : public ICommand{
public:
    enum class Type{
        Normal,
        SubCommand
    };

    Command(const std::string& name,
            std::function<std::optional<TConfig>&()> commandGetter,
            Type type)
        : info_(name, {}, {})
        , commandGetter_(commandGetter)
        , type_(type)
    {
    }

    OptionInfo& info() override
    {
        return info_;
    }

    const OptionInfo& info() const override
    {
        return info_;
    }

private:
    void read(const std::vector<std::string>& commandLine) override
    {
        auto& commandCfg = commandGetter_();
        if (!commandCfg.has_value()){
            commandCfg.emplace();
            if (helpFlag_){
                detail::ConfigAccess<TConfig>(*commandCfg).addFlag(std::move(helpFlag_));
                detail::ConfigAccess<TConfig>(*commandCfg).addHelpFlagToCommands(programName_);
            }
        }
        commandCfg->read(commandLine);
    }

    void enableHelpFlag(const std::string& programName) override
    {
        programName_ = programName + " " + info_.name();
        using NameProvider = typename detail::Format<detail::ConfigAccess<TConfig>::format()>::nameProvider;
        helpFlag_ = std::make_unique<detail::Flag>(NameProvider::name("help"),
                                                   std::string{},
                                                   [this]()->bool&{return helpFlagValue_;},
                                                   detail::Flag::Type::Exit);
        helpFlag_->info().addDescription("show usage info and exit");
    }

    bool isHelpFlagSet() const override
    {
        return helpFlagValue_;
    }

    bool isSubCommand() const override
    {
        return type_ == Type::SubCommand;
    }

    std::string usageInfo() const override
    {
        auto& commandCfg = commandGetter_();
        if (!commandCfg.has_value())
            return {};
        return commandCfg->usageInfo(programName_);
    }

    std::string usageInfoDetailed() const override
    {
        auto& commandCfg = commandGetter_();
        if (!commandCfg.has_value())
            return {};
        return commandCfg->usageInfoDetailed(programName_);
    }

    std::vector<not_null<ICommand*>> commandList() override
    {
        auto& commandCfg = commandGetter_();
        if (!commandCfg.has_value())
            return {};
        return detail::ConfigAccess<TConfig>(*commandCfg).commandList();
    }

private:
    OptionInfo info_;
    std::function<std::optional<TConfig>&()> commandGetter_;
    Type type_;
    std::string programName_;
    std::unique_ptr<IFlag> helpFlag_;
    bool helpFlagValue_ = false;
};

template<typename T, typename TConfig>
class CommandCreator{
    using NameProvider = typename Format<ConfigAccess<TConfig>::format()>::nameProvider;
public:
    CommandCreator(TConfig& cfg,
                   const std::string& varName,
                   std::function<std::optional<T>&()> commandGetter,
                   typename Command<T>::Type type = Command<T>::Type::Normal)
        : cfg_(cfg)
    {
        static_assert (std::is_base_of_v<Config<ConfigAccess<TConfig>::format()>, T>,
                       "Command's type must be a subclass of Config<FormatType> and have the same format as its parent config.");
        Expects(!varName.empty());
        command_ = std::make_unique<Command<T>>(NameProvider::fullName(varName), commandGetter, type);
    }

    CommandCreator<T, TConfig>& operator<<(const std::string& info)
    {
        command_->info().addDescription(info);
        return *this;
    }

    CommandCreator<T, TConfig>& operator<<(const Name& customName)
    {
        command_->info().resetName(customName.value());
        return *this;
    }

    operator std::optional<T>()
    {
        ConfigAccess<TConfig>{cfg_}.addCommand(std::move(command_));
        return std::optional<T>{};
    }

private:
    std::unique_ptr<Command<T>> command_;
    TConfig& cfg_;
};

template <typename T, typename TConfig>
CommandCreator<T, TConfig> makeCommandCreator(TConfig& cfg,
                                              const std::string& varName,
                                              std::function<std::optional<T>&()> commandGetter,
                                              typename Command<T>::Type type = Command<T>::Type::Normal)
{
    return CommandCreator<T, TConfig>{cfg, varName, commandGetter, type};
}

}

max_stars_repo_name: GerHobbelt/hypertextcpp

max_issues_repo_name: GerHobbelt/hypertextcpp

max_forks_repo_name: GerHobbelt/hypertextcpp


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpMath.h>
#include <gbpCosmo_linear_theory.h>
#include <gsl/gsl_multimin.h>

typedef struct init_gbpCosmo2gbpCosmo_integrand_params_struct init_gbpCosmo2gbpCosmo_integrand_params;
struct init_gbpCosmo2gbpCosmo_integrand_params_struct {
    double                     inv_s;
    double                     z_source;
    double                     z_target;
    double                     R_1;
    double                     R_2;
    cosmo_info **              cosmo_source;
    cosmo_info **              cosmo_target;
    int                        n_int;
    gsl_integration_workspace *wspace;
};

double init_gbpCosmo2gbpCosmo_integrand(double R, void *params_in) {
    init_gbpCosmo2gbpCosmo_integrand_params *params       = (init_gbpCosmo2gbpCosmo_integrand_params *)params_in;
    double                                   inv_s        = params->inv_s;
    double                                   z_source     = params->z_source;
    double                                   z_target     = params->z_target;
    cosmo_info **                            cosmo_source = params->cosmo_source;
    cosmo_info **                            cosmo_target = params->cosmo_target;
    return (pow(1. - (sigma_R(cosmo_source, inv_s * R, z_source, PSPEC_LINEAR_TF, PSPEC_ALL_MATTER)) /
                         (sigma_R(cosmo_target, R, z_target, PSPEC_LINEAR_TF, PSPEC_ALL_MATTER)),
                2.) /
            R);
}

double init_gbpCosmo2gbpCosmo_minimize_function(const gsl_vector *v_i, void *params_in) {
    // Set variable integrand parameters
    init_gbpCosmo2gbpCosmo_integrand_params *params = (init_gbpCosmo2gbpCosmo_integrand_params *)params_in;
    params->inv_s                                   = gsl_vector_get(v_i, 0);
    params->z_target                                = gsl_vector_get(v_i, 1);

    // Perform integral to minimize
    gsl_function integrand;
    double       delta_i;
    double       abs_error;
    integrand.function = init_gbpCosmo2gbpCosmo_integrand;
    integrand.params   = params_in;
    gsl_integration_qag(&integrand, params->R_1, params->R_2, 0, 1e-3, params->n_int, GSL_INTEG_GAUSS61, params->wspace, &delta_i, &abs_error);
    return (delta_i / take_ln(params->R_2 / params->R_1));
}

void init_gbpCosmo2gbpCosmo(cosmo_info **           cosmo_source,
                            cosmo_info **           cosmo_target,
                            double                  z_min,
                            double                  M_min,
                            double                  M_max,
                            gbpCosmo2gbpCosmo_info *gbpCosmo2gbpCosmo) {
    SID_log("Initializing cosmology scaling...", SID_LOG_OPEN | SID_LOG_TIMER);
    SID_set_verbosity(SID_SET_VERBOSITY_RELATIVE, -1);

    // Store some infor in the gbpCosmo2gbpCosmo_info structure
    gbpCosmo2gbpCosmo->M_min        = M_min;
    gbpCosmo2gbpCosmo->M_max        = M_max;
    gbpCosmo2gbpCosmo->z_min        = z_min;
    gbpCosmo2gbpCosmo->cosmo_source = (*cosmo_source);
    gbpCosmo2gbpCosmo->cosmo_target = (*cosmo_target);

    // Perform minimization
    // const gsl_multimin_fminimizer_type *T=gsl_multimin_fminimizer_nmsimplex2;
    const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
    gsl_multimin_fminimizer *           s = NULL;
    gsl_vector *                        ss, *x;
    gsl_multimin_function               minex_func;

    // Starting point
    x = gsl_vector_alloc(2);
    gsl_vector_set(x, 0, 1.);    // inv_s
    gsl_vector_set(x, 1, z_min); // z_scaled

    // Set initial step sizes to 1
    ss = gsl_vector_alloc(2);
    gsl_vector_set_all(ss, 1.0);

    // Set parameters
    init_gbpCosmo2gbpCosmo_integrand_params params;
    params.cosmo_source = cosmo_source;
    params.cosmo_target = cosmo_target;
    params.z_source     = z_min;
    params.R_1          = R_of_M(M_min, *cosmo_source);
    params.R_2          = R_of_M(M_max, *cosmo_source);
    params.inv_s        = gsl_vector_get(x, 0);
    params.z_target     = gsl_vector_get(x, 1);
    params.n_int        = 100;
    params.wspace       = gsl_integration_workspace_alloc(params.n_int);

    // Initialize method
    minex_func.n      = 2;
    minex_func.f      = init_gbpCosmo2gbpCosmo_minimize_function;
    minex_func.params = (void *)(&params);
    s                 = gsl_multimin_fminimizer_alloc(T, 2);
    gsl_multimin_fminimizer_set(s, &minex_func, x, ss);

    // Perform minimization
    double size;
    int    status;
    size_t iter     = 0;
    size_t iter_max = 200;
    do {
        iter++;
        status = gsl_multimin_fminimizer_iterate(s);
        if(status)
            SID_exit_error("Error encountered during minimisation in init_gbpCosmo2gbpCosmo() (status=%d).",
                           SID_ERROR_LOGIC, status);
        size   = gsl_multimin_fminimizer_size(s);
        status = gsl_multimin_test_size(size, 1e-2);
    } while(status == GSL_CONTINUE && iter <= iter_max);
    if(status != GSL_SUCCESS)
        SID_exit_error("Failed to converge during minimisation in init_gbpCosmo2gbpCosmo() (status=%d,iter=%d).",
                       SID_ERROR_LOGIC, status, iter);

    // Finalize results
    double Omega_M_source           = ((double *)ADaPS_fetch(*cosmo_source, "Omega_M"))[0];
    double H_Hubble_source          = 1e2 * ((double *)ADaPS_fetch(*cosmo_source, "h_Hubble"))[0];
    double Omega_M_target           = ((double *)ADaPS_fetch(*cosmo_target, "Omega_M"))[0];
    double H_Hubble_target          = 1e2 * ((double *)ADaPS_fetch(*cosmo_target, "h_Hubble"))[0];
    gbpCosmo2gbpCosmo->s_L          = 1. / gsl_vector_get(s->x, 0);
    gbpCosmo2gbpCosmo->s_M          = (Omega_M_target * H_Hubble_target) / (Omega_M_source * H_Hubble_source) * pow((gbpCosmo2gbpCosmo->s_L), 3.);
    gbpCosmo2gbpCosmo->z_min_scaled = gsl_vector_get(s->x, 1);
    ;

    // Calculate growth factors needed for
    //    determining redshift mappings
    gbpCosmo2gbpCosmo->D_prime_z_min = linear_growth_factor(z_min, cosmo_target);
    gbpCosmo2gbpCosmo->D_z_scaled    = linear_growth_factor(gbpCosmo2gbpCosmo->z_min_scaled, cosmo_source);
    gbpCosmo2gbpCosmo->D_ratio       = gbpCosmo2gbpCosmo->D_prime_z_min / gbpCosmo2gbpCosmo->D_z_scaled;

    // Clean-up
    gsl_vector_free(x);
    gsl_vector_free(ss);
    gsl_multimin_fminimizer_free(s);
    gsl_integration_workspace_free(params.wspace);
    SID_set_verbosity(SID_SET_VERBOSITY_DEFAULT);
    SID_log("Done.", SID_LOG_CLOSE);
}

ext: c

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#pragma once

#include <winrt\Windows.Foundation.h>
#include <d3d11.h>
#include <DirectXMath.h>
#include <gsl\gsl>

namespace Library
{
	class SamplerStates final
	{
	public:
		inline static winrt::com_ptr<ID3D11SamplerState> TrilinearWrap;
		inline static winrt::com_ptr<ID3D11SamplerState> TrilinearMirror;
		inline static winrt::com_ptr<ID3D11SamplerState> TrilinearClamp;
		inline static winrt::com_ptr<ID3D11SamplerState> TrilinearBorder;
		inline static winrt::com_ptr<ID3D11SamplerState> PointClamp;
		inline static winrt::com_ptr<ID3D11SamplerState> DepthMap;
		inline static winrt::com_ptr<ID3D11SamplerState> ShadowMap;
		inline static winrt::com_ptr<ID3D11SamplerState> PcfShadowMap;

		static DirectX::XMVECTORF32 BorderColor;
		static DirectX::XMVECTORF32 ShadowMapBorderColor;

		static void Initialize(gsl::not_null<ID3D11Device*> direct3DDevice);
		static void Shutdown();

		SamplerStates() = delete;
		SamplerStates(const SamplerStates&) = delete;
		SamplerStates& operator=(const SamplerStates&) = delete;
		SamplerStates(SamplerStates&&) = delete;
		SamplerStates& operator=(SamplerStates&&) = delete;
		~SamplerStates() = default;
	};
}

ext: h


text:
/*
 * Copyright 2009-2018 The VOTCA Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#ifndef _GSL_BOOST_UBLAS_MATRIX_PROD_
#define _GSL_BOOST_UBLAS_MATRIX_PROD_

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

#include <gsl/gsl_blas.h>

namespace boost { namespace numeric { namespace ublas {

    /*
     * GSL has separate implementations for floats and doubles
     * hence we first have specializations for these two types
     * TSeparate templates for diagonal and symmetric matreces
     * and their combinations. To check  whether a specific
     * routine is called, define GSLDEBUG
     */

    // partial specialization for double precision (dgemm))
    template<class F, class A>
    inline matrix<double,F,A>    // prod( m1, m2 )
    prod(const matrix<double,F,A> &m1, const matrix<double,F,A> &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;34mGSL [CMD,CM]\x1b[0;39m\n" << std::flush;
       #endif

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

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

       gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                  1.0, &mA.matrix, &mB.matrix,
                  0.0, &mC.matrix);
       return AxB;
    }

    // partial specialization for single precision (sgemm))
    template<class F, class A>
    inline matrix<float,F,A>    // prod( m1, m2 )
    prod(const matrix<float,F,A> &m1, const matrix<float,F,A> &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;34mGSL [CMS,CM]\x1b[0;39m\n" << std::flush;
       #endif

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

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

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

    // full specialization
    template<class T, class F, class A>
    inline matrix<T,F,A>    // prod( m1, m2 )
    prod(const matrix<T,F,A> &m1, const matrix<T,F,A> &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;34mGSL [CM,CM]\x1b[0;39m -> " << std::flush;
       #endif

       return prod(m1,m2);
    }

    /// transpose products
    template<class T, class F, class A>
    inline matrix<T,F,A>    // prod( trans(m1), m2 )
    prod(const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m1, const matrix<T,F,A> &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;32mGSL [CTM,CM]\x1b[0;39m -> " << std::flush;
       #endif
       boost::numeric::ublas::matrix<T,F,A> _m1 = m1;
       return prod(_m1,m2);
    }

    template<class T, class F, class A>
    inline matrix<T,F,A>    // prod( m1, trans(m2) )
    prod(const matrix<T,F,A> &m1, const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;32mGSL [CM,CTM]\x1b[0;39m -> " << std::flush;
       #endif

       const boost::numeric::ublas::matrix<T,F,A> _m2 = m2;
       return prod(m1,_m2);
    }

    template<class T, class F, class A>
    inline matrix<T,F,A>    // prod( trans(m1), trans(m2) )
    prod(const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m1, const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;32mGSL [CTM,CTM]\x1b[0;39m -> " << std::flush;
       #endif

       boost::numeric::ublas::matrix<T,F,A> _m1 = m1;
       boost::numeric::ublas::matrix<T,F,A> _m2 = m2;

       return prod(_m1,_m2);
    }

    /// diagonal matrix
    template<class T, class F, class L, class A>
    inline matrix<T,F,A>    // prod( diagonal m1, m2 )
    prod(const diagonal_matrix<T,L,A> &m1, const matrix<T,F,A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CDM, CM]\x1b[0;39m -> " << std::flush;
       #endif
       const matrix<T,F,A> _m1 = m1;
       return prod(_m1,m2);
    }

    template<class T, class F, class L, class A>
    inline matrix<T,F,A>    // prod( m1, diagonal m2 )
    prod(const matrix<T,F,A> &m1, const diagonal_matrix<T,L,A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CM, CDM]\x1b[0;39m -> " << std::flush;
       #endif
       const matrix<T,F,A> _m2 = m2;
       return prod(m1,_m2);
    }

    template<class T, class F, class L, class A>
    inline matrix<T,F,A>    // prod( diagonal m1, diagonal m2 )
    prod(const diagonal_matrix<T,L,A> &m1, const diagonal_matrix<T,L,A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CDM, CM]\x1b[0;39m -> " << std::flush;
       #endif
       const matrix<T,F,A> _m1 = m1;
       const matrix<T,F,A> _m2 = m2;
       return prod(_m1,_m2);
    }

    template<class T, class F, class L, class A>
    inline matrix<T,F,A>    // prod( diagonal m1, transpose m2 )
    prod(const diagonal_matrix<T,L,A> &m1, const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CDM, CTM]\x1b[0;39m -> " << std::flush;
       #endif

       const matrix<T,F,A> _m1 = m1;
       const matrix<T,F,A> _m2 = m2;

       return prod(_m1,_m2);

    }

   template<class T, class F, class L, class A>
    inline matrix<T,F,A>    // prod( transpose m1, diagonal m2 )
    prod(const matrix_unary2<matrix<T,F,A>,scalar_identity<T> > &m1, const diagonal_matrix<T,L,A> &m2)
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CTM, CDM]\x1b[0;39m -> " << std::flush;
       #endif

       const matrix<T,F,A> _m1 = m1;
       const matrix<T,F,A> _m2 = m2;

       return prod(_m1,_m2);

    }

    // symmetric matrix
    template<class T, class F, class A, class TRI, class L>
    inline matrix<T,F,A>    // prod( symmetric m1, m2 )
    prod(symmetric_matrix<T, TRI, L, A> &m1, matrix<T,F,A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CSM, CM]\x1b[0;39m -> " << std::flush;
       #endif

       assert( m1.size1() == m2.size1() );
       const matrix<T,F,A> _m1 = m1;

       return prod(_m1, m2 );
    }

    template<class T, class F, class A, class TRI, class L>
    inline matrix<T,F,A>    // prod( m1, symmetric m2 )
    prod( matrix<T,F,A> &m1, symmetric_matrix<T, TRI, L, A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CSM, CM]\x1b[0;39m -> " << std::flush;
       #endif

       assert( m1.size1() == m2.size1() );
       const matrix<T,F,A> _m2 = m2;

       return prod(m1, _m2 );
    }

    template<class T, class F, class A, class TRI, class L>
    inline matrix<T,F,A>    // prod( symmetric m1, symmetric m2 )
    prod(symmetric_matrix<T, TRI, L, A> &m1, symmetric_matrix<T, TRI, L, A> &m2 )
    {
       #ifdef GSLDEBUG
          std::cout << "\x1b[0;33mGSL [CSM, CM]\x1b[0;39m -> " << std::flush;
       #endif

       assert( m1.size1() == m2.size1() );
       const matrix<T,F,A> _m1 = m1;
       const matrix<T,F,A> _m2 = m2;

       return prod(_m1, _m2 );
    }

}}}

#endif  // _GSL_BOOST_UBLAS_MATRIX_PROD_

max_stars_repo_name: mbarbry/ctp

max_issues_repo_name: mbarbry/ctp

max_forks_repo_name: mbarbry/ctp


text:
#ifndef GSLMATRIX_H
#define GSLMATRIX_H

#include <gsl/gsl_matrix_double.h> //Science Library header file
#include <gsl/gsl_linalg.h>
#include "Matrix.h"

template<typename T>
class GSLMatrix : public Matrix<T>
{
private:
	gsl_matrix *_matrix; //Pointer to the gsl_matrix struct
	void setMatrix(gsl_matrix* newMatrix);

public:
	GSLMatrix(int nRow, int nCol);
	~GSLMatrix();

	void resize(int nRow, int nCol); //If the matrix has elements in it should it keep those elements or is it ok if they are lost?

	T getElement(int i, int j) const;
	void setElement(int i, int j, T value);
	gsl_matrix* getMatrix() const;

	int getSize() const;
	void setZero();
	void invert();

	Matrix<T>* operator*(const Matrix<T>* RHS);
	Matrix<T>* operator+(const Matrix<T>* RHS);
	Matrix<T>* operator-(const Matrix<T>* RHS);
	void operator=(const GSLMatrix<T>* RHS);

};

template<typename T>
GSLMatrix<T>::GSLMatrix(int nRow, int nCol)
{
	_matrix = gsl_matrix_alloc(nRow, nCol);
}

template<typename T>
GSLMatrix<T>::~GSLMatrix()
{
	if(_matrix != NULL)
		gsl_matrix_free(_matrix);
}

template<typename T>
void GSLMatrix<T>::setMatrix(gsl_matrix* newMatrix)
{
	//Free the old matrix, if there is one
	if(_matrix != NULL)
		gsl_matrix_free(_matrix);

	_matrix = newMatrix;
}

template<typename T>
void GSLMatrix<T>::resize(int nRow, int nCol)
{
	if(_matrix != NULL)
		gsl_matrix_free(_matrix);

	_matrix = NULL;

	_matrix = gsl_matrix_alloc(nRow, nCol);
}

template<typename T>
T GSLMatrix<T>::getElement(int i, int j) const
{
	return gsl_matrix_get(_matrix, i, j);
}

template<typename T>
void GSLMatrix<T>::setElement(int i, int j, T value)
{
	gsl_matrix_set(_matrix, i, j, value);
}

template<typename T>
void GSLMatrix<T>::invert()
{
	int sigNum;

	//gsl_matrix for the output
	gsl_matrix* outputMatrix = gsl_matrix_alloc(_matrix->size1, _matrix->size1 );

	//Allocate permutation matrix
	gsl_permutation* permutationM = gsl_permutation_alloc(_matrix->size1);

	//get the LU decomposition of the matrix
	gsl_linalg_LU_decomp((gsl_matrix*) _matrix, permutationM, &sigNum);

	//invert the matrix
	gsl_linalg_LU_invert(_matrix, permutationM, outputMatrix);

	this->setMatrix(outputMatrix);

	gsl_permutation_free(permutationM);
}

template<typename T>
gsl_matrix* GSLMatrix<T>::getMatrix() const
{
	return _matrix;
}

template<typename T>
int GSLMatrix<T>::getSize() const
{
	return _matrix->size1;
}

template<typename T>
void GSLMatrix<T>::setZero()
{
	gsl_matrix_set_zero (_matrix);
}

template<typename T>
Matrix<T>* GSLMatrix<T>::operator*(const Matrix<T>* RHS)
{
  //Check if its the same size
  if(_matrix->size1 != RHS->getSize() || _matrix->size1 != _matrix->size1)
    return NULL;

  //Check if its the same type (GSL Matrix)
  if(dynamic_cast<const GSLMatrix<T>*> (RHS) == NULL)
    throw Error (InvalidParam, "GSLMatrix<T>::operator*",
		 "RHS is not a GSLMatrix");

  //Create the new result matrix
  GSLMatrix<T> *result = new GSLMatrix<T>(_matrix->size1, _matrix->size1);
  result->setZero();

  //For each row
  for(int row = 0; row < _matrix->size1; row++)
  {
    //for each element in that row
    for(int element = 0; element < _matrix->size2; element++)
    {
      //calculate the value for the current element
      for(int i = 0; i < _matrix->size1; i++)
      {
	result->setElement(row, element, (result->getElement(row, element) ) + this->getElement(row, i) * RHS->getElement(i,element));
      }
    }
  }

  return result;
}

template<typename T>
Matrix<T>* GSLMatrix<T>::operator+(const Matrix<T>* RHS)
{
  //Check if its the same size
  if(_matrix->size1 != RHS->getSize() || _matrix->size1 != _matrix->size1)
    return NULL;

  //Check if its the same type (GSL Matrix)
  if(dynamic_cast<const GSLMatrix<T>*> (RHS) == NULL)
    throw Error (InvalidParam, "GSLMatrix<T>::operator*",
		 "RHS is not a GSLMatrix");

	GSLMatrix<T> *result = new GSLMatrix<T>(_matrix->size1, _matrix->size1); //Create the new result matrix
	result->setZero();

	//For each row
	for(int row = 0; row < _matrix->size1; row++)
	{
		//for each element in that row
		for(int element = 0; element < _matrix->size1; element++)
		{
			result->setElement(row, element, this->getElement(row, element) + RHS->getElement(row, element) );
		}
	}

	return result;
}

template<typename T>
Matrix<T>* GSLMatrix<T>::operator-(const Matrix<T>* RHS)
{
	if(_matrix->size1 != RHS->getSize() || _matrix->size1 != _matrix->size1) //Check if its the same size
			return NULL;

  //Check if its the same type (GSL Matrix)
  if(dynamic_cast<const GSLMatrix<T>*> (RHS) == NULL)
    throw Error (InvalidParam, "GSLMatrix<T>::operator*",
		 "RHS is not a GSLMatrix");

	GSLMatrix<T> *result = new GSLMatrix<T>(RHS->getSize(), RHS->getSize()); //Create the new result matrix
	result->setZero();

	//For each row
	for(int row = 0; row < _matrix->size1; row++)
	{
		//for each element in that row
		for(int element = 0; element < _matrix->size1; element++)
		{
			result->setElement(row, element, this->getElement(row, element) - RHS->getElement(row, element) );
		}
	}

	return result;
}

template<typename T>
void GSLMatrix<T>::operator=(const GSLMatrix<T>* RHS)
{
	if(_matrix->size1 != RHS->getSize() )
	{
		gsl_matrix_free(_matrix);
		gsl_matrix_alloc(RHS->getSize(), RHS->getSize());
	}

	gsl_matrix_memcpy(_matrix, RHS->getMatrix());
}

#endif


text:
#ifndef __LLVUTILS_H__
#define __LLVUTILS_H__ 1

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_cdf.h>

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"
#include "wip.h"
#include "likelihood.h"
#include "splinecoeffs.h"
#include "LLVFDresponse.h"
#include "LLVnoise.h"

/***************** Structures for parameters *****************/

/* Parameters for the generation of a LLV waveform (in the form of a list of modes) */
typedef struct tagLLVParams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double m1;                 /* mass of companion 1 (solar masses) */
  double m2;                 /* mass of companion 2 (solar masses) */
  double distance;           /* distance of source (Mpc) */
  double ra;                 /* right ascension of the source (rad) */
  double dec;                /* declination of the source (rad) */
  double inclination;        /* inclination of L relative to line of sight (rad) */
  double polarization;       /* polarization angle (rad) */
  int nbmode;                /* number of modes to generate (starting with 22) - defaults to 5 (all modes) */
} LLVParams;

/* Global parameters for the waveform generation and overlap computation */
typedef struct tagLLVGlobalParams {
  double fRef;               /* reference frequency (Hz, default 0 which is interpreted as Mf=0.14) */
  double minf;               /* Minimal frequency (Hz) - when set to 0 (default), use the first frequency covered by the noise data of the detector */
  double maxf;               /* Maximal frequency (Hz) - when set to 0 (default), use the last frequency covered by the noise data of the detector */
  int setphiRefatfRef;       /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) (default=1) */
  int nbmodeinj;             /* number of modes to include in the injection (starting with 22) - defaults to 5 (all modes) */
  int nbmodetemp;            /* number of modes to include in the templates (starting with 22) - defaults to 5 (all modes) */
  int tagint;                /* Tag choosing the integrator: 0 for wip (default), 1 for linear integration */
  int tagnetwork;            /* Tag choosing the network of detectors to use */
  int nbptsoverlap;          /* Number of points to use in loglinear overlaps (default 32768) */
  int constL;                /* set all logLikelihood to 0 - allows to sample from the prior for testing */
} LLVGlobalParams;

typedef struct tagLLVPrior {
	double deltaT;             /* width of time prior centered on injected value (s) (default 0.1) */
	double comp_min;           /* minimum component mass (solar masses) (default 4) */
	double comp_max;           /* maximum component mass (solar masses) (default 50) */
	double mtot_min;           /* minimum total mass (solar masses) (default 8) */
	double mtot_max;           /* maximum total mass (solar masses) (default 100) */
	double qmax;               /* maximum asymmetric mass ratio (>=1) (default 12) */
	double dist_min;           /* minimum distance of source (Mpc) (default 1) */
	double dist_max;           /* maximum distance of source (Mpc) (default 1e4) */
  double ra_min;            /* minimum ra (rad, default 0) - for testing */
  double ra_max;            /* maximum ra (rad, default 2pi) - for testing */
  double dec_min;           /* minimum dec (rad, default 0) - for testing */
  double dec_max;           /* maximum dec (rad, default pi) - for testing */
  double phase_min;          /* minimum phase (rad, default 0) - for testing */
  double phase_max;          /* maximum phase (rad, default 2pi) - for testing */
  double pol_min;            /* minimum polarization (rad, default 0) - for testing */
  double pol_max;            /* maximum polarization (rad, default 2pi) - for testing */
  double inc_min;            /* minimum inclination (rad, default 0) - for testing */
  double inc_max;            /* maximum inclination (rad, default pi) - for testing */
  double fix_m1;
  double fix_m2;
  double fix_time;
  double fix_ra;
  double fix_dec;
  double fix_phase;
  double fix_pol;
  double fix_dist;
  double fix_inc;
  int pin_m1;
  int pin_m2;
  int pin_time;
  int pin_ra;
  int pin_dec;
  int pin_phase;
  int pin_pol;
  int pin_dist;
  int pin_inc;
  double snr_target;
  int rescale_distprior;
  int flat_distprior;
} LLVPrior;

typedef struct tagLLVRunParams {
	double eff;                /* target efficiency (default 0.1) */
	double tol;                /* logZ tolerance (default 0.5) */
	int    nlive;              /* number of live points (default 1000) */
	char   outroot[200];       /* output root (default "chains/LLVinference_") */
	int    bambi;              /* run BAMBI? (default 0) */
	int    resume;             /* resume form previous run? (default 0) */
	int    maxiter;            /* max number of iterations (default 0 - ignore) */
	char   netfile[200];       /* NN settings file (default "LLVinference.inp") */
  int    mmodal;             /* use multimodal decomposition ? */
  int    maxcls;             /* max number of modes in multimodal decomposition */
  int    nclspar;            /* number of parameters to use for multimodal decomposition - in the order of the cube */
  double ztol;               /* in multimodal decomposition, modes with lnZ lower than Ztol are ignored */
  int    seed;               /* seed the inference by setting one of the live points to the injection ? */
} LLVRunParams;

/************ Structures for signals and injections ************/

// typedef struct tagLLVSignal
// {
//   struct tagListmodesCAmpPhaseFrequencySeries* LHOSignal;   /* Signal in LHO, in the form of a list of the contribution of each mode */
//   struct tagListmodesCAmpPhaseFrequencySeries* LLOSignal;   /* Signal in LLO, in the form of a list of the contribution of each mode */
//   struct tagListmodesCAmpPhaseFrequencySeries* VIRGOSignal; /* Signal in VIRGO, in the form of a list of the contribution of each mode */
//   double LHOhh;                                             /* Inner product (h|h) for LHO */
//   double LLOhh;                                             /* Inner product (h|h) for LLO */
//   double VIRGOhh;                                           /* Inner product (h|h) for VIRGO */
// } LLVSignal;

typedef struct tagLLVSignalCAmpPhase
{
  struct tagListmodesCAmpPhaseFrequencySeries* LHOSignal;   /* Signal in LHO, in the form of a list of the contribution of each mode */
  struct tagListmodesCAmpPhaseFrequencySeries* LLOSignal;   /* Signal in LLO, in the form of a list of the contribution of each mode */
  struct tagListmodesCAmpPhaseFrequencySeries* VIRGOSignal; /* Signal in VIRGO, in the form of a list of the contribution of each mode */
  double LLVhh;                                             /* Combined Inner product (h|h) for dectectors LHV */
} LLVSignalCAmpPhase;

typedef struct tagLLVInjectionCAmpPhase
{
  struct tagListmodesCAmpPhaseSpline* LHOSplines;   /* Signal in LHO, in the form of a list of splines for the contribution of each mode */
  struct tagListmodesCAmpPhaseSpline* LLOSplines;   /* Signal in LLO, in the form of a list of splines for the contribution of each mode */
  struct tagListmodesCAmpPhaseSpline* VIRGOSplines;   /* Signal in VIRGO, in the form of a list of splines for the contribution of each mode */
  double LLVss;                                   /* Combined Inner product (s|s) for dectectors LHV */
} LLVInjectionCAmpPhase;

typedef struct tagLLVSignalReIm /* We don't store the SNRs here, as we will use -1/2(h-s|h-s) for the likelihood */
{
  struct tagReImFrequencySeries* LHOSignal;   /* Signal in LHO, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* LLOSignal;   /* Signal in LLO, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* VIRGOSignal;   /* Signal in VIRGO, in the form of a Re/Im frequency series where the modes have been summed */
} LLVSignalReIm;

typedef struct tagLLVInjectionReIm /* Storing the vectors of frequencies and noise values - We don't store the SNRs here, as we will use -1/2(h-s|h-s) for the likelihood */
{
  struct tagReImFrequencySeries* LHOSignal;   /* Signal in LHO, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* LLOSignal;   /* Signal in LLO, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* VIRGOSignal;   /* Signal in VIRGO, in the form of a Re/Im frequency series where the modes have been summed */
  gsl_vector* freq;                            /* Vector of frequencies of the injection (assumed to be the same for LHO, LLO, VIRGO) */
  gsl_vector* noisevaluesLHO;                    /* Vector of noise values on freq LHO */
  gsl_vector* noisevaluesLLO;                    /* Vector of noise values on freq LLO */
  gsl_vector* noisevaluesVIRGO;                  /* Vector of noise values on freq VIRGO */
} LLVInjectionReIm;

/************ Functions for LLV parameters, injection, likelihood, prior ************/

/* Parsing parameters for the generation of a LLV waveform, from the command line */
/* Masses are input in solar masses and distances in Mpc - converted in SI for the internals */
void parse_args_LLV(ssize_t argc, char **argv,
    LLVParams* params,
    LLVGlobalParams* globalparams,
    LLVPrior* prior,
    LLVRunParams* run,
    LLVParams* addparams);

/* Functions to print the parameters of the run in files for reference */
int print_parameters_to_file_LLV(
  LLVParams* params,
  LLVGlobalParams* globalparams,
  LLVPrior* prior,
  LLVRunParams* run);
int print_rescaleddist_to_file_LLV(
  LLVParams* params,
  LLVGlobalParams* globalparams,
  LLVPrior* prior,
  LLVRunParams* run);

/* Initialization and clean-up for LLVSignal structures */
void LLVSignalCAmpPhase_Cleanup(LLVSignalCAmpPhase* signal);
void LLVSignalCAmpPhase_Init(LLVSignalCAmpPhase** signal);
void LLVInjectionCAmpPhase_Cleanup(LLVInjectionCAmpPhase* signal);
void LLVInjectionCAmpPhase_Init(LLVInjectionCAmpPhase** signal);
void LLVSignalReIm_Cleanup(LLVSignalReIm* signal);
void LLVSignalReIm_Init(LLVSignalReIm** signal);
void LLVInjectionReIm_Cleanup(LLVInjectionReIm* signal);
void LLVInjectionReIm_Init(LLVInjectionReIm** signal);

/* Function generating a LLV signal as a list of modes in CAmp/Phase form, from LLV parameters */
int LLVGenerateSignalCAmpPhase(
  struct tagLLVParams* params,                 /* Input: set of LLV parameters of the signal */
  struct tagLLVSignalCAmpPhase* signal);  /* Output: structure for the generated signal */
/* Function generating a LLV injection as a list of modes, given as preinterpolated splines, from LLV parameters */
int LLVGenerateInjectionCAmpPhase(
  struct tagLLVParams* injectedparams,    /* Input: set of LLV parameters of the signal */
  struct tagLLVInjectionCAmpPhase* signal);  /* Output: structure for the generated signal */
/* Function generating a LLV signal as a frequency series in Re/Im form where the modes have been summed, from LLV parameters - takes as argument the frequencies on which to evaluate */
int LLVGenerateSignalReIm(
  struct tagLLVParams* params,       /* Input: set of LLV parameters of the template */
  gsl_vector* freq,                   /* Input: frequencies on which evaluating the waveform (from the injection) */
  struct tagLLVSignalReIm* signal);  /* Output: structure for the generated signal */
/* Function generating a LLV injection as a frequency series in Re/Im form where the modes have been summed, from LLV parameters - frequencies on which to evaluate are to be determined internally */
int LLVGenerateInjectionReIm(
  struct tagLLVParams* injectedparams,       /* Input: set of LLV parameters of the injection */
  double fLow,                                /* Input: starting frequency (from argument minf) */
  double fHigh,                               /* Input: upper frequency (from argument maxf) */
  int nbpts,                                  /* Input: number of frequency samples */
  int tagsampling,                            /* Input: tag for using linear (0) or logarithmic (1) sampling */
  struct tagLLVInjectionReIm* signal);       /* Output: structure for the generated signal */

// checks prior boundaires
int PriorBoundaryCheck(LLVPrior *prior, double *Cube);

// Prior functions from Cube to physical parameters
// x1 is min, x2 is max when specified
// r is Cube value
double CubeToFlatPrior(double r, double x1, double x2);
double CubeToLogFlatPrior(double r, double x1, double x2);
double CubeToPowerPrior(double p, double r, double x1, double x2);
double CubeToGaussianPrior(double r, double mean, double sigma);
double CubeToSinPrior(double r, double x1, double x2);
double CubeToCosPrior(double r, double x1, double x2);

/* Prior functions from physical parameters to Cube
   x1 is min, x2 is max when specified
   y is physical value */
double FlatPriorToCube(double y, double x1, double x2);
double LogFlatPriorToCube(double y, double x1, double x2);
double PowerPriorToCube(double p, double y, double x1, double x2);
double SinPriorToCube(double y, double x1, double x2);
double CosPriorToCube(double y, double x1, double x2);

/* log-Likelihood functions */
double CalculateLogLCAmpPhase(LLVParams *params, LLVInjectionCAmpPhase* injection);
double CalculateLogLReIm(LLVParams *params, LLVInjectionReIm* injection);

/************ Global Parameters ************/

extern LLVParams* injectedparams;
extern LLVGlobalParams* globalparams;
extern LLVPrior* priorParams;
double logZdata;

#endif

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
#include <sys/sysinfo.h>
#include <sys/time.h>

#define PRINT_MAX 100 // don't print more than 100 elements
#define MEAN 0
#define STANDARD_DEVIATION 1

struct arguments {
        double *array;
        int start;
        int end;
        unsigned long seed;
};

void * initialize_array (void *arguments) {

        /* Get arguments */
        struct arguments *args = (struct arguments *) arguments;
        double *array = args->array;
        int start = args->start;
        int end = args->end;
        unsigned long seed = args->seed;

        /* Create a new random number generator */
        const gsl_rng_type *T;
        gsl_rng *r;

        gsl_rng_env_setup(); // set the default values for the random number generator variables

        T = gsl_rng_default;
        r = gsl_rng_alloc (T);
        gsl_rng_set(r, seed); // setup the generator using the seed we got from the main function

        /* Using the GSL library, initialize the array with normal random values */
        for (int i = start; i < end; i++) {
                array[i] = MEAN + gsl_ran_gaussian (r, STANDARD_DEVIATION);;
        }

        /* Free memory */
        gsl_rng_free (r);

        /* Exit thread when finished */
        pthread_exit(NULL);
}

int main (int argc, char* argv[])
{
        int i;
        int number_of_elements;
        int number_of_threads;
        if (argc == 3) {
                number_of_elements = atoi(argv[1]);
                number_of_threads = atoi(argv[2]);
        } else {
                fprintf(stderr, "%s:\n", "Wrong number of arguments.");
                fprintf(stderr, "%s:\n", "Use: ./avoid_race_condition number_of_elements number_of_threads");
                exit(EXIT_FAILURE);
        }

        int elements_per_thread = number_of_elements / number_of_threads;

        double *array = (double *) malloc(number_of_elements * sizeof(double));
        pthread_t *threads = (pthread_t *) malloc(number_of_threads * sizeof(pthread_t));
        struct arguments *args = (struct arguments *) malloc(number_of_threads * sizeof(struct arguments));;

        srand(time(NULL)); // need to to use a new seed each time we run the program

        // create threads
	for (i = 0; i < number_of_threads; i++) {
                args[i].array = array;
                args[i].start = i * elements_per_thread;
                if (i == number_of_threads - 1) { // in the last thread, process the remaining elements in case
                                                  // that NUMBER_OF_ELEMENTS is not divisible by NUMBER_OF_THREADS
                        args[i].end = number_of_elements;
                } else {
                        args[i].end = args[i]. start + elements_per_thread;
                }
                args[i].seed = random();

                // create the threads
                if(pthread_create(&threads[i], NULL, initialize_array, (void *) &args[i]) == -1) {
                        fprintf(stderr, "%s %d\n", "Can't create thread number", i);
                }
	}

	// wait until every thread ends
        for (i = 0; i < number_of_threads; i++) {
		if (pthread_join(threads[i], NULL) == -1) {
			fprintf(stderr, "%s %d\n", "Can't join thread number", i);
		}
	}

        // print numbers
        if(number_of_elements <= PRINT_MAX) {
                for (i = 0; i < number_of_elements; i++) {
                        printf (" %f", array[i]);
                }
                printf("\n");
        } else {
                printf("(Array size = %d. Printing only its first and last elements)\n", number_of_elements);
                printf("%f %f ... %f %f\n", array[0], array[1], array[number_of_elements -2], array[number_of_elements -1]);
        }

        // free memory
        free(array);
        free(threads);
	free(args);

        return 0;
}

max_stars_repo_name: sernamar/random-numbers

max_issues_repo_name: sernamar/random-numbers

max_forks_repo_name: sernamar/random-numbers


text:
/**
 *
 * @file qwrapper_zlaset.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_zlaset(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       PLASMA_Complex64_t alpha, PLASMA_Complex64_t beta,
                       PLASMA_Complex64_t *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_zlaset_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(PLASMA_Complex64_t),         &alpha, VALUE,
        sizeof(PLASMA_Complex64_t),         &beta,  VALUE,
        sizeof(PLASMA_Complex64_t)*LDA*N,    A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlaset_quark = PCORE_zlaset_quark
#define CORE_zlaset_quark PCORE_zlaset_quark
#endif
void CORE_zlaset_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    PLASMA_Complex64_t alpha;
    PLASMA_Complex64_t beta;
    PLASMA_Complex64_t *A;
    int LDA;

    quark_unpack_args_7(quark, uplo, M, N, alpha, beta, A, LDA);
    LAPACKE_zlaset_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, alpha, beta, A, LDA);
}


text:
/**
 * Copyright (c) 2019-2020 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Imported from:
 *
 * @file core_blas.h
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Jakub Kurzak
 * @author Hatem Ltaief
 * @date 2010-11-15
 *
 **/
#ifndef _PLASMA_CORE_BLAS_H_
#define _PLASMA_CORE_BLAS_H_

#include "cores/dplasma_plasmatypes.h"
#include "cores/descriptor.h"
#include <cblas.h>
#include <lapacke.h>

#include "cores/core_zblas.h"
#include "cores/core_dblas.h"
#include "cores/core_cblas.h"
#include "cores/core_sblas.h"
#include "cores/core_zcblas.h"
#include "cores/core_dsblas.h"

#ifdef __cplusplus
extern "C" {
#endif

  /*
   * Coreblas Error
   */
#define coreblas_error(k, str) fprintf(stderr, "%s: Parameter %d / %s\n", __func__, k, str);

/* CBLAS requires for scalar arguments to be passed by address rather than by value */
#ifndef CBLAS_SADDR
#define CBLAS_SADDR( _val_ ) &(_val_)
#endif

 /** ****************************************************************************
  *  External interface of the GKK algorithm for InPlace Layout Translation
  **/
int  GKK_minloc(int n, int *T);
void GKK_BalanceLoad(int thrdnbr, int *Tp, int *leaders, int nleaders, int L);
int  GKK_getLeaderNbr(int me, int ne, int *nleaders, int **leaders);

void CORE_pivot_update(int m, int n, int *ipiv, int *indices,
                       int offset, int init);

#ifdef __cplusplus
}
#endif

#endif /* _PLASMA_CORE_BLAS_H_ */

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
#pragma once
#include <gsl/gsl>
#include <map>
#include "halley/text/halleystring.h"
#include "halley/maths/vector4.h"

namespace Halley
{
	class Path;
	class Image;
	struct ImageData;

	class AsepriteExternalReader
	{
	public:
		static std::vector<ImageData> importAseprite(String baseName, gsl::span<const gsl::byte> fileData, bool trim);

	private:
		static std::vector<ImageData> loadImagesFromPath(Path tmp, bool crop);
		static std::map<int, int> getSpriteDurations(Path jsonPath);
		static void processFrameData(String baseName, std::vector<ImageData>& frameData, std::map<int, int> durations);
	};

	class AsepriteReader
	{
	public:
		static std::vector<ImageData> importAseprite(String baseName, gsl::span<const gsl::byte> fileData, bool trim);
	};
}

max_stars_repo_name: lye/halley

max_issues_repo_name: lye/halley

max_forks_repo_name: lye/halley


text:
/* cdf/betainv.c
 *
 * Copyright (C) 2004 Free Software Foundation, Inc.
 * Copyright (C) 2006, 2007 Brian Gough
 * Written by Jason H. Stover.
 * Modified for GSL by 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307, USA.
 */

/*
 * Invert the Beta distribution.
 *
 * References:
 *
 * Roger W. Abernathy and Robert P. Smith. "Applying Series Expansion
 * to the Inverse Beta Distribution to Find Percentiles of the
 * F-Distribution," ACM Transactions on Mathematical Software, volume
 * 19, number 4, December 1993, pages 474-480.
 *
 * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions of a
 * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39,
 * number 8, August 1968, pages 1264-1273.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_randist.h>

#include "error.h"

static double
bisect (double x, double P, double a, double b, double xtol, double Ptol)
{
  double x0 = 0, x1 = 1, Px;

  while (fabs(x1 - x0) > xtol) {
    Px = gsl_cdf_beta_P (x, a, b);
    if (fabs(Px - P) < Ptol) {
      /* return as soon as approximation is good enough, including on
         the first iteration */
      return x;
    } else if (Px < P) {
      x0 = x;
    } else if (Px > P) {
      x1 = x;
    }
    x = 0.5 * (x0 + x1);
  }
  return x;
}

double
gsl_cdf_beta_Pinv (const double P, const double a, const double b)
{
  double x, mean;

  if (P < 0.0 || P > 1.0)
    {
      CDF_ERROR ("P must be in range 0 < P < 1", GSL_EDOM);
    }

  if (a < 0.0)
    {
      CDF_ERROR ("a < 0", GSL_EDOM);
    }

  if (b < 0.0)
    {
      CDF_ERROR ("b < 0", GSL_EDOM);
    }

  if (P == 0.0)
    {
      return 0.0;
    }

  if (P == 1.0)
    {
      return 1.0;
    }

  if (P > 0.5)
    {
      return gsl_cdf_beta_Qinv (1 - P, a, b);
    }

  mean = a / (a + b);

  if (P < 0.1)
    {
      /* small x */

      double lg_ab = gsl_sf_lngamma (a + b);
      double lg_a = gsl_sf_lngamma (a);
      double lg_b = gsl_sf_lngamma (b);

      double lx = (log (a) + lg_a + lg_b - lg_ab + log (P)) / a;
      if (lx <= 0) {
        x = exp (lx);             /* first approximation */
        x *= pow (1 - x, -(b - 1) / a);   /* second approximation */
      } else {
        x = mean;
      }

      if (x > mean)
        x = mean;
    }
  else
    {
      /* Use expected value as first guess */
      x = mean;
    }

  /* Do bisection to get closer */
  x = bisect (x, P, a, b, 0.01, 0.01);

  {
    double lambda, dP, phi;
    unsigned int n = 0;

  start:
    dP = P - gsl_cdf_beta_P (x, a, b);
    phi = gsl_ran_beta_pdf (x, a, b);

    if (dP == 0.0 || n++ > 64)
      goto end;

    lambda = dP / GSL_MAX (2 * fabs (dP / x), phi);

    {
      double step0 = lambda;
      double step1 = -((a - 1) / x - (b - 1) / (1 - x)) * lambda * lambda / 2;

      double step = step0;

      if (fabs (step1) < fabs (step0))
        {
          step += step1;
        }
      else
        {
          /* scale back step to a reasonable size when too large */
          step *= 2 * fabs (step0 / step1);
        };

      if (x + step > 0 && x + step < 1)
        {
          x += step;
        }
      else
        {
          x = sqrt (x) * sqrt (mean);   /* try a new starting point */
        }

      if (fabs (step0) > 1e-10 * x)
        goto start;
    }

  end:

    if (fabs(dP) > GSL_SQRT_DBL_EPSILON * P)
      {
        GSL_ERROR_VAL("inverse failed to converge", GSL_EFAILED, GSL_NAN);
      }

    return x;
  }
}

double
gsl_cdf_beta_Qinv (const double Q, const double a, const double b)
{

  if (Q < 0.0 || Q > 1.0)
    {
      CDF_ERROR ("Q must be inside range 0 < Q < 1", GSL_EDOM);
    }

  if (a < 0.0)
    {
      CDF_ERROR ("a < 0", GSL_EDOM);
    }

  if (b < 0.0)
    {
      CDF_ERROR ("b < 0", GSL_EDOM);
    }

  if (Q == 0.0)
    {
      return 1.0;
    }

  if (Q == 1.0)
    {
      return 0.0;
    }

  if (Q > 0.5)
    {
      return gsl_cdf_beta_Pinv (1 - Q, a, b);
    }
  else
    {
      return 1 - gsl_cdf_beta_Pinv (Q, b, a);
    };
}


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 rho_NFW(double r, double M_vir, double z, int mode, cosmo_info **cosmo) {
    double c_vir;
    double R_vir;
    double g_c;
    double rho_o;
    double x;
    set_NFW_params(M_vir, z, mode, cosmo, &c_vir, &R_vir);
    g_c   = 1. / (log(1. + c_vir) - c_vir / (1. + c_vir));
    rho_o = M_vir * g_c / (FOUR_PI * pow(R_vir / c_vir, 3.));
    x     = r * c_vir / R_vir;
    return (rho_o / (x * pow(1. + x, 2.)));
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2020-Present Couchbase, Inc.
 *
 *   Use of this software is governed by the Business Source License included
 *   in the file licenses/BSL-Couchbase.txt.  As of the Change Date specified
 *   in that file, in accordance with the Business Source License, use of this
 *   software will be governed by the Apache License, Version 2.0, included in
 *   the file licenses/APL2.txt.
 */

#pragma once

#include <folly/lang/Assume.h>
#include <gsl/gsl-lite.hpp>

#include <array>
#include <ratio>
#include <string>
#include <string_view>

namespace cb::stats {

// Enum of the relevant base units used by convention in Prometheus
enum class BaseUnit {
    None, // used for text stats where units aren't relevant
    Count, // Generic count which does not meet a better unit
    Seconds,
    Bytes,
    Ratio,
};

/**
 * Type representing a unit (kilobytes, microseconds etc.).
 * This is encoded as a BaseUnit (see above) and a scaling factor stored as
 * two integers (see std::ratio).
 *
 * Thus,
 *     kilobytes -> bytes, 1000/1
 *     microseconds -> seconds, 1/1000000
 *
 * Given a numerical value in this unit, the scaling factor indicates
 * how to convert the value back to the base unit.
 *
 *     10 kilobytes
 *     10 * 1000 / 1
 *     -> 10000 bytes
 *
 * This can be performed using the toBaseUnit method -
 *
 * Unit(std::kilo{}, BaseUnit::Bytes).toBaseUnit(10) == 10000
 *
 * Units will usually be used through the constexpr values
 * in namespace units
 *
 * units::kilobytes.toBaseUnit(10) == 10000
 */
class Unit {
public:
    explicit constexpr Unit(BaseUnit baseUnit)
        : Unit(std::ratio<1>{}, baseUnit) {
    }

    template <class RatioType>
    constexpr Unit(RatioType ratio, BaseUnit baseUnit)
        : numerator(gsl::narrow_cast<int64_t>(RatioType::num)),
          denominator(gsl::narrow_cast<int64_t>(RatioType::den)),
          baseUnit(baseUnit) {
    }

    /**
     * Scale a value of the current unit (e.g., milliseconds) to the base
     * unit (e.g., seconds).
     */
    [[nodiscard]] constexpr double toBaseUnit(double value) const {
        return (value * numerator) / denominator;
    }

    [[nodiscard]] std::string_view getSuffix() const {
        using namespace std::string_view_literals;
        switch (baseUnit) {
        case BaseUnit::None:
        case BaseUnit::Count:
            return ""sv;
        case BaseUnit::Seconds:
            return "_seconds"sv;
        case BaseUnit::Bytes:
            return "_bytes"sv;
        case BaseUnit::Ratio:
            return "_ratio"sv;
        }
        folly::assume_unreachable();
    }

private:
    int64_t numerator;
    int64_t denominator;
    BaseUnit baseUnit;
};

namespace units {
constexpr Unit none{std::ratio<1>{}, BaseUnit::None};

constexpr Unit count{std::ratio<1>{}, BaseUnit::Count};

// floating point between 0 and 1, this is already in the correct base unit
constexpr Unit ratio{std::ratio<1>{}, BaseUnit::Ratio};
// percents should be scaled down to ratios for Prometheus
constexpr Unit percent{std::ratio<1, 100>{}, BaseUnit::Ratio};

// time units
constexpr Unit nanoseconds{std::nano{}, BaseUnit::Seconds};
constexpr Unit microseconds{std::micro{}, BaseUnit::Seconds};
constexpr Unit milliseconds{std::milli{}, BaseUnit::Seconds};
constexpr Unit seconds{std::ratio<1>{}, BaseUnit::Seconds};
constexpr Unit minutes{std::ratio<60>{}, BaseUnit::Seconds};
constexpr Unit hours{std::ratio<60 * 60>{}, BaseUnit::Seconds};
constexpr Unit days{std::ratio<60 * 60 * 24>{}, BaseUnit::Seconds};

// Note: std ratio definitions are powers of 10, not 2
// e.g., kilo = std::ratio<1000, 1>
// so 1 kilobyte = 1000 bytes
// if powers of 2 are explicitly needed, IEC prefixes
// (kibi, gibi, tebi) could easily be defined.

// byte units
constexpr Unit bits{std::ratio<1, 8>{}, BaseUnit::Bytes};
constexpr Unit bytes{std::ratio<1>{}, BaseUnit::Bytes};
constexpr Unit kilobytes{std::kilo{}, BaseUnit::Bytes};
constexpr Unit megabytes{std::mega{}, BaseUnit::Bytes};
constexpr Unit gigabytes{std::giga{}, BaseUnit::Bytes};
} // namespace units
} // namespace cb::stats


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 painterport_eca0917d_5109_4de2_ae6f_18d5c163d02c_h
#define painterport_eca0917d_5109_4de2_ae6f_18d5c163d02c_h

#include <gslib/rtree.h>
#include <ariel/painter.h>

__ariel_begin__

class painter_obj;

typedef rtree_entity<painter_obj*> painter_obj_entity;
typedef rtree_node<painter_obj_entity> painter_obj_node;
typedef _tree_allocator<painter_obj_node> painter_obj_alloc;
typedef tree<painter_obj_entity, painter_obj_node, painter_obj_alloc> painter_obj_tree;
typedef rtree<painter_obj_entity, quadratic_split_alg<25, 10, painter_obj_tree>, painter_obj_node, painter_obj_alloc> painter_obj_rtree;
typedef vector<painter_obj*> painter_objs;

class __gs_novtable painter_obj abstract
{
public:
    enum type
    {
        po_line,
        po_rect,
        po_path,
        po_text,
    };

public:
    painter_obj(const painter_context& ctx): _ctx(ctx) {}
    virtual ~painter_obj() {}
    virtual type get_type() const = 0;
    virtual rectf& get_rect(rectf& rc) const = 0;

protected:
    painter_context         _ctx;

public:
    const painter_context& get_context() const { return _ctx; }
};

class painter_line_obj:
    public painter_obj
{
public:
    painter_line_obj(const painter_context& ctx, const pointf& p1, const pointf& p2): painter_obj(ctx), _p1(p1), _p2(p2) {}
    virtual type get_type() const override { return po_line; }
    virtual rectf& get_rect(rectf& rc) const override;

private:
    pointf                  _p1;
    pointf                  _p2;
};

class painter_rect_obj:
    public painter_obj
{
public:
    painter_rect_obj(const painter_context& ctx,  const rectf& rc): painter_obj(ctx), _rc(rc) {}
    virtual type get_type() const override { return po_rect; }
    virtual rectf& get_rect(rectf& rc) const override { return rc = _rc; }

private:
    rectf                   _rc;
};

class painter_path_obj:
    public painter_obj
{
public:
    painter_path_obj(const painter_context& ctx, painter_path& take_me): painter_obj(ctx) { _path.swap(take_me); }
    painter_path_obj(const painter_context& ctx, const painter_path& path): painter_obj(ctx), _path(path) {}
    virtual type get_type() const override { return po_path; }
    virtual rectf& get_rect(rectf& rc) const override;
    rectf& rebuild_rect(rectf& rc) const;

private:
    painter_path            _path;
    mutable bool            _rc_valid = false;
    mutable rectf           _rc;
};

class painter_text_obj:
    public painter_obj
{
public:
    painter_text_obj(const painter_context& ctx, const string& txt, const pointf& p, const color& cr): painter_obj(ctx), _text(txt), _pos(p), _cr(cr) {}
    virtual type get_type() const override { return po_text; }
    virtual rectf& get_rect(rectf& rc) const override;

private:
    string                  _text;
    pointf                  _pos;
    color                   _cr;
};

class painterport:
    public painter
{
public:
    virtual ~painterport();
    virtual void resize(int w, int h) override {}
    virtual void draw_path(const painter_path& path) override;
    virtual void draw_line(const vec2& p1, const vec2& p2) override;
    virtual void draw_rect(const rectf& rc) override;
    virtual void draw_text(const gchar* str, float x, float y, const color& cr, int length) override;
    virtual void on_draw_begin() override;

protected:
    painter_obj_rtree       _rtree;

public:
    rectf& get_area_rect(rectf& rc) const;
    void query_objs(painter_objs& objs, const rectf& rc);

protected:
    void add_obj(painter_obj* obj);
    void clear_objs();
};

__ariel_end__

#endif

max_stars_repo_name: lymastee/gslib

max_issues_repo_name: lymastee/gslib

max_forks_repo_name: lymastee/gslib


text:
#include <ceed.h>
#include <petsc.h>
#include "../problems/mooney-rivlin.h"

// Build libCEED context object
PetscErrorCode PhysicsContext_MR(MPI_Comm comm, Ceed ceed, Units *units,
                                 CeedQFunctionContext *ctx) {
  PetscErrorCode ierr;
  Physics_MR phys;

  PetscFunctionBegin;

  ierr = PetscMalloc1(1, units); CHKERRQ(ierr);
  ierr = PetscMalloc1(1, &phys); CHKERRQ(ierr);
  ierr = ProcessPhysics_MR(comm, phys, *units); CHKERRQ(ierr);
  CeedQFunctionContextCreate(ceed, ctx);
  CeedQFunctionContextSetData(*ctx, CEED_MEM_HOST, CEED_COPY_VALUES,
                              sizeof(*phys), phys);
  ierr = PetscFree(phys); CHKERRQ(ierr);

  PetscFunctionReturn(0);
}

// Build libCEED smoother context object
PetscErrorCode PhysicsSmootherContext_MR(MPI_Comm comm, Ceed ceed,
    CeedQFunctionContext ctx, CeedQFunctionContext *ctx_smoother) {
  PetscErrorCode ierr;
  PetscScalar nu_smoother = 0;
  PetscBool nu_flag = PETSC_FALSE;
  Physics_MR phys, phys_smoother;

  PetscFunctionBegin;

  ierr = PetscOptionsBegin(comm, NULL,
                           "Mooney Rivlin physical parameters for smoother", NULL);
  CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-nu_smoother", "Poisson's ratio for smoother",
                            NULL, nu_smoother, &nu_smoother, &nu_flag);
  CHKERRQ(ierr);
  if (nu_smoother < 0 ||
      nu_smoother >= 0.5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                                    "Mooney-Rivlin model requires Poisson ratio -nu option in [0, .5)");

  ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics

  if (nu_flag) {
    // Copy context
    CeedQFunctionContextGetData(ctx, CEED_MEM_HOST, &phys);
    ierr = PetscMalloc1(1, &phys_smoother); CHKERRQ(ierr);
    ierr = PetscMemcpy(phys_smoother, phys, sizeof(*phys)); CHKERRQ(ierr);
    CeedQFunctionContextRestoreData(ctx, &phys);
    // Create smoother context
    CeedQFunctionContextCreate(ceed, ctx_smoother);
    phys_smoother->lambda = 2 * (phys_smoother->mu_1 + phys_smoother->mu_2) *
                            nu_smoother / (1 - 2*nu_smoother);
    CeedQFunctionContextSetData(*ctx_smoother, CEED_MEM_HOST, CEED_COPY_VALUES,
                                sizeof(*phys_smoother), phys_smoother);
    ierr = PetscFree(phys_smoother); CHKERRQ(ierr);
  } else {
    *ctx_smoother = NULL;
  }

  PetscFunctionReturn(0);
}

// Process physics options - Mooney-Rivlin
PetscErrorCode ProcessPhysics_MR(MPI_Comm comm, Physics_MR phys, Units units) {
  PetscErrorCode ierr;
  PetscReal nu = -1;
  phys->mu_1 = -1;
  phys->mu_2 = -1;
  phys->lambda = -1;
  units->meter     = 1;        // 1 meter in scaled length units
  units->second    = 1;        // 1 second in scaled time units
  units->kilogram  = 1;        // 1 kilogram in scaled mass units

  PetscFunctionBeginUser;

  ierr = PetscOptionsBegin(comm, NULL, "Mooney Rivlin physical parameters", NULL);
  CHKERRQ(ierr);

  ierr = PetscOptionsScalar("-mu_1", "Material Property mu_1", NULL,
                            phys->mu_1, &phys->mu_1, NULL); CHKERRQ(ierr);
  if (phys->mu_1 < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                                "Mooney-Rivlin model requires non-negative -mu_1 option (Pa)");

  ierr = PetscOptionsScalar("-mu_2", "Material Property mu_2", NULL,
                            phys->mu_2, &phys->mu_2, NULL); CHKERRQ(ierr);
  if (phys->mu_2 < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                                "Mooney-Rivlin model requires non-negative -mu_2 option (Pa)");

  ierr = PetscOptionsScalar("-nu", "Poisson ratio", NULL,
                            nu, &nu, NULL); CHKERRQ(ierr);
  if (nu < 0 || nu >= 0.5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                                     "Mooney-Rivlin model requires Poisson ratio -nu option in [0, .5)");
  phys->lambda = 2 * (phys->mu_1 + phys->mu_2) * nu / (1 - 2*nu);

  ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units",
                            NULL, units->meter, &units->meter, NULL);
  CHKERRQ(ierr);
  units->meter = fabs(units->meter);

  ierr = PetscOptionsScalar("-units_second", "1 second in scaled time units",
                            NULL, units->second, &units->second, NULL);
  CHKERRQ(ierr);
  units->second = fabs(units->second);

  ierr = PetscOptionsScalar("-units_kilogram", "1 kilogram in scaled mass units",
                            NULL, units->kilogram, &units->kilogram, NULL);
  CHKERRQ(ierr);
  units->kilogram = fabs(units->kilogram);

  ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics

  // Define derived units
  units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second));

  // Scale material parameters based on units of Pa
  phys->mu_1 *= units->Pascal;
  phys->mu_2 *= units->Pascal;
  phys->lambda *= units->Pascal;

  PetscFunctionReturn(0);
};

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED