text:
#ifndef __MAIN_H__
#define __MAIN_H__ 1

#ifdef PARALLEL
#include "mpi.h"
#endif
#include <stdio.h>
#include <math.h>
#include <iostream>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <gsl/gsl_randist.h>
//#include "NNopt.h"
#include "bambi.h"

#ifdef __INTEL_COMPILER 			// if the MultiNest library was compiled with ifort
       #define NESTRUN nested_mp_nestrun_
#elif defined __GNUC__ 				// if the MultiNest library was compiled with gfortran
       #define NESTRUN __nested_MOD_nestrun
#else
       #error Do not know how to link to Fortran libraries, check symbol table for your platform (nm libnest3.a | grep nestrun) & edit example_eggbox_C++/eggbox.cc
#endif

extern float thL[3],logLRange,tol;
extern float *omicron;
extern double *maxsigma;
extern double logZero;
extern int nNN,nNNerr,totpar,loglcalls,ncheck,myid,nproc;
extern char root[100],networkinputs[100];
extern bool likenetinit,converged,lastconverged,netres,firstrun,discardpts;
extern int ignoredbambicalls,counter;
extern size_t nlayers,nnodes[10];
extern int doBAMBI,useNN,whitenin,whitenout,resume;

/***************************************** C++ Interface to MultiNest **************************************************/

namespace nested
{

	// map the Fortran 90 entry points of libnest3.a to C++ functions

	// module nested, function nestRun maps to nested::run

	// the pass-by-reference nature of most of the Fortran is translated away
	// *apart* from the callbacks. The provided call back functions must still accept
	// references rather than values. There is also some confusion as to the type
	// of the first argument of LogLike.
	// Should it be a double * or an farray<double, 1> *? The former seems to
	// work and is simpler.

	// This structure is reverse engineered from looking
	// at gfortran stack traces. It is probably wrong

	template<typename type, int ndims> class farray_traits;

	template<> class farray_traits<double, 1> { public: static const int id = 537; };
	template<> class farray_traits<double, 2> { public: static const int id = 538; };
	template<> class farray_traits<int, 1> { public: static const int id = 265; };
	template<> class farray_traits<int, 2> { public: static const int id = 266; };

	// the extra data for f90 that defines how arrays are arranged.
	template<typename T, int ndim> class farray
	{
		public:
			farray(T *_data, int w, int h = 0) : data(_data), offset(0), type(farray_traits<T, ndim>::id),
			x_stride(1), x_lbound(1), x_ubound(w), y_stride(w), y_lbound(1), y_ubound(h) {};

			T *data;
			int offset;
			int type;
			int x_stride, x_lbound, x_ubound;
			int y_stride, y_lbound, y_ubound;
	};

	extern "C" {
		void NESTRUN(int &mmodal, int &ceff, int &nlive, double &tol, double &efr, int &ndims,
			int &nPar, int &nClsPar, int &maxModes, int &updInt, double &Ztol, char *root, int &seed,
			int *pWrap, int &fb, int &resume, int &outfile, int &initMPI, double &logZero, int &maxiter,
			void (*Loglike)(double *Cube, int &n_dim, int &n_par, double &lnew, void *),
			void (*dumper)(int &, int &, int &, double **, double **, double **, double &, double &, double &, void *),
			void (*bambi)(int &, int &, double **, double &), void *context, int &root_len);
	}

	static void run(bool mmodal, bool ceff, int nlive, double tol, double efr, int ndims, int nPar, int nClsPar, int maxModes,
		int updInt, double Ztol, const std::string & root, int seed, int *pWrap, bool fb, bool resume, bool outfile,
		bool initMPI, double logZero, int maxiter, void (*LogLike)(double *Cube, int &n_dim, int &n_par, double &lnew, void *),
		void (*dumper)(int &, int &, int &, double **, double **, double **, double &, double &, double &, void *),
		void (*bambi)(int &, int &, double **, double &), void *context)
	{
		char t_root[100];
		std::fill(t_root, t_root + 100, ' ');
		snprintf(t_root, 99, "%s", root.c_str());
		int root_len = strlen(t_root);
		t_root[strlen(t_root)] = ' ';

		int t_fb = fb;
		int t_resume = resume;
		int t_outfile = outfile;
		int t_initMPI = initMPI;
		int t_mmodal = mmodal;
		int t_ceff = ceff;

		NESTRUN(t_mmodal, t_ceff, nlive, tol, efr, ndims, nPar, nClsPar, maxModes, updInt, Ztol, t_root, seed, pWrap, t_fb,
		t_resume, t_outfile, t_initMPI, logZero, maxiter, LogLike, dumper, bambi, context, root_len);
	}
}

/***********************************************************************************************************************/

#endif


text:
#ifndef DELTAFREADER_H
#define DELTAFREADER_H

#include "ParameterReader.h"
#include "readindata.h"
#include "GaussThermal.h"
#include <fstream>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp.h>

using namespace std;

class Deltaf_Data
{
    private:
        ParameterReader * paraRdr;

        int hrg_eos; // type of pdg file for hadron resonance gas EoS
        int mode; //type of freezeout surface, VH or VAH
        int df_mode; // type of delta-f correction (e.g. 14-moment, CE, or modified distribution)
        int include_baryon;

        string hrg_eos_path;
        string urqmd = "deltaf_coefficients/vh/urqmd/"; // directories of df coefficient tables
        string smash = "deltaf_coefficients/vh/smash/";
        string smash_box = "deltaf_coefficients/vh/smash_box/";
    public:
        int points_T;
        int points_muB;

        double T_min;
        double muB_min;

        double dT;
        double dmuB;

        double * T_array;
        double * muB_array;
        //  Coefficients of 14 moment approximation (vhydro)
        // df ~ ((c0-c2)m^2 + b.c1(u.p) + (4c2-c0)(u.p)^2).Pi + (b.c3 + c4(u.p))p_u.V^u + c5.p_u.p_v.pi^uv
        double ** c0_data;
        double ** c1_data;
        double ** c2_data;
        double ** c3_data;
        double ** c4_data;

        //  Coefficients of Chapman Enskog expansion (vhydro)
        // df ~ ((c0-c2)m^2 + b.c1(u.p) + (4c2-c0)(u.p)^2).Pi + (b.c3 + c4(u.p))p_u.V^u + c5.p_u.p_v.pi^uv
        double ** F_data;
        double ** G_data;
        double ** betabulk_data;
        double ** betaV_data;
        double ** betapi_data;

        // cubic splines of the coefficients as function of temperature only (neglect muB, nB, Vmu)
        // (c1, G = 0) for muB = 0 and (c3, c4, betaV) aren't needed since they couple to baryon diffusion
        // so in the cubic spline evaluation: just set (G, c1, c3, c4) = 0 and betaV = 1 (betaV is in denominator)

        gsl_spline * c0_spline;
        gsl_spline * c2_spline;
        gsl_spline * c3_spline;

        gsl_spline * F_spline;
        gsl_spline * betabulk_spline;
        gsl_spline * betaV_spline;
        gsl_spline * betapi_spline;

        // Jonah coefficients
        const int jonah_points = 301;       // # lambda interpolation points
        const double lambda_min = -1.0;     // lambda min / max values
        const double lambda_max = 2.0;
        const double delta_lambda = (lambda_max - lambda_min) / ((double)jonah_points - 1.0);

        double * lambda_squared_array;      // squared isotropic momentum scale
        double * z_array;                   // renormalization factor (apart from detLambda)
        double * bulkPi_over_Peq_array;     // bulk pressure output
        double bulkPi_over_Peq_max;         // the maximum bulk pressure in the array

        gsl_spline * lambda_squared_spline; // cubic splines for lambda^2(bulkPi/Peq) and z(bulkPi/Peq)
        gsl_spline * z_spline;

        Deltaf_Data(ParameterReader * paraRdr_in);
        ~Deltaf_Data();

        void load_df_coefficient_data();    // read the data files in /deltaf_coefficients/vh

        void construct_cubic_splines();

        // I skip the photon because I think it breaks down for lambda = -1
        void compute_jonah_coefficients(particle_info * particle_data, int Nparticle);

        deltaf_coefficients evaluate_df_coefficients(double T, double muB, double E, double P, double bulkPi);

        deltaf_coefficients cubic_spline(double T, double E, double P, double bulkPi);

        double calculate_bilinear(double ** f_data, double T, double muB, double TL, double TR, double muBL, double muBR, int iTL, int iTR, int imuBL, int imuBR);

        deltaf_coefficients bilinear_interpolation(double T, double muB, double E, double P, double bulkPi);

        void test_df_coefficients(double bulkPi_over_P);

        void compute_particle_densities(particle_info * particle_data, int Nparticle);
};

#endif

ext: h


text:
/* Copyright 2019-2020 Canaan Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
#include "allocator.h"
#include "model.h"
#include "result.h"
#include "runtime_module.h"
#include <gsl/gsl-lite.hpp>
#include <memory>
#include <unordered_map>

BEGIN_NS_NNCASE_RUNTIME

class NNCASE_API options_dict
{
public:
    template <class T>
    result<T> get(const char *name)
    {
        auto it = values_.find(name);
        if (it != values_.end())
            return ok(it->second.as<T>());
        else
            return err(std::errc::result_out_of_range);
    }

    template <class T>
    result<void> set(const char *name, T value)
    {
        values_[name] = scalar(value);
        return ok();
    }

private:
    std::unordered_map<const char *, scalar> values_;
};

class NNCASE_API interpreter
{
public:
    interpreter() noexcept;
    interpreter(interpreter &) = delete;
    interpreter(interpreter &&) = default;

    NNCASE_NODISCARD result<void> load_model(gsl::span<const gsl::byte> buffer) noexcept;

    size_t inputs_size() const noexcept;
    size_t outputs_size() const noexcept;
    const memory_range &input_desc(size_t index) const noexcept;
    const memory_range &output_desc(size_t index) const noexcept;
    const runtime_shape_t &input_shape(size_t index) const noexcept;
    const runtime_shape_t &output_shape(size_t index) const noexcept;
    result<runtime_tensor> input_tensor(size_t index) noexcept;
    result<void> input_tensor(size_t index, runtime_tensor tensor) noexcept;
    result<runtime_tensor> output_tensor(size_t index) noexcept;
    result<void> output_tensor(size_t index, runtime_tensor tensor) noexcept;

    result<void> run() noexcept;

    result<runtime_module *> find_module_by_id(size_t index) noexcept;
    options_dict &options() noexcept;

private:
    std::vector<std::unique_ptr<runtime_module>> modules_;
    runtime_module *main_module_;
    options_dict options_;
};

END_NS_NNCASE_RUNTIME


text:
/*
 * The MIT License (MIT)
 *
 * Copyright (c) 2016 Josef Gajdusek
 *
 * 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.
 * */

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

#include "log.h"

#include "coords.h"

static projPJ proj_wgs;

void init_projs()
{
	proj_wgs = pj_init_plus("+init=epsg:4326");
}

cl_float2 wgs84_to_meters(cl_float2 wgs, projPJ proj_meters)
{
	cl_float2 ret;
	// The X and Ys are switched intentionally
	double ry = wgs.x * DEG_TO_RAD;
	double rx = wgs.y * DEG_TO_RAD;
	int err = pj_transform(proj_wgs, proj_meters, 1, 1, &ry, &rx, NULL);
	if (err) {
		log_error("Coordinate conversion failed: %s", pj_strerrno(err));
	}
	ret.x = rx;
	ret.y = ry;
	return ret;
}

void generate_translation_tile(int xtile, int ytile, int zoom, cl_float4 *out, projPJ proj_meters)
{
	int side = 20;
	int npoints = side * side;
	gsl_multifit_linear_workspace *gwsp = gsl_multifit_linear_alloc(npoints, 3);
	gsl_vector *yx = gsl_vector_alloc(npoints);
	gsl_vector *yy = gsl_vector_alloc(npoints);
	gsl_matrix *X = gsl_matrix_alloc(npoints, 3);
	gsl_vector *cx = gsl_vector_alloc(3);
	gsl_vector *cy = gsl_vector_alloc(3);
	gsl_matrix *cov = gsl_matrix_alloc(npoints, npoints);
	double chisq;

	// Generate side x side training points
	for (int x = 0; x < side; x++) {
		for (int y = 0; y < side; y++) {
			cl_float2 tile = {.x = xtile + (double)x / side,
							  .y = ytile + (double)y / side};
			cl_float2 wgs = tile_to_wgs84(tile, zoom);
			cl_float2 mets = wgs84_to_meters(wgs, proj_meters);

			int pt = x * side + y;
			// Inputs
			gsl_matrix_set(X, pt, 0, tile.x - xtile);
			gsl_matrix_set(X, pt, 1, tile.y - ytile);
			gsl_matrix_set(X, pt, 2,	 1);

			// Outputs
			gsl_vector_set(yx, pt, mets.x);
			gsl_vector_set(yy, pt, mets.y);

		}
	}

	// Fit
	gsl_multifit_linear(X, yx, cx, cov, &chisq, gwsp);
	gsl_multifit_linear(X, yy, cy, cov, &chisq, gwsp);

	gsl_matrix *mtrans = gsl_matrix_alloc(3, 2);
	gsl_matrix_set_col(mtrans, 0, cx);
	gsl_matrix_set_col(mtrans, 1, cy);

	// Output

	out[0].x = gsl_matrix_get(mtrans, 0, 0);
	out[0].y = gsl_matrix_get(mtrans, 1, 0);
	out[0].z = gsl_matrix_get(mtrans, 2, 0);
	out[1].x = gsl_matrix_get(mtrans, 0, 1);
	out[1].y = gsl_matrix_get(mtrans, 1, 1);
	out[1].z = gsl_matrix_get(mtrans, 2, 1);

	gsl_matrix_free(mtrans);
	gsl_matrix_free(cov);
	gsl_vector_free(cy);
	gsl_vector_free(cx);
	gsl_matrix_free(X);
	gsl_vector_free(yy);
	gsl_vector_free(yx);
	gsl_multifit_linear_free(gwsp);
}

ext: c

max_stars_repo_name: atalax/cl-heatmap

max_issues_repo_name: atalax/cl-heatmap

max_forks_repo_name: atalax/cl-heatmap


text:
#ifndef _ASA_H
#define _ASA_H

#include <iostream>
#include <memory>
#include <type_traits>
#include <vector>
#include <array>

#include <gsl/gsl_siman.h>

namespace gsl{

class SimulatedAnnealing{
public:
    SimulatedAnnealing();
    virtual ~SimulatedAnnealing();

    typedef gsl_siman_params_t Parameters;
    /* Struct defined in gsl/gsl_siman.h (gsl2.5) access on Jan 9, 2019 as:
     * 51 typedef struct {
     * 52   int n_tries;          // how many points to try for each step
     * 53   int iters_fixed_T;    // how many iterations at each temperature?
     * 54   double step_size;     // max step size in the random walk
     * 55   // the following parameters are for the Boltzmann distribution
     * 56   double k, t_initial, mu_t, t_min;
     * 57 } gsl_siman_params_t;
     */

protected:
    const gsl_rng_type * workspace;
    gsl_rng * randomNumberGenerator;

    double (*energy) (void*);
    double (*measure)(void*,
                      void*);
    void   (*step)   (const gsl_rng*,
                            void*,
                            double);
    void   (*print)  (void*);

    Parameters asaParameters;

public:
    Parameters& set_parameters(const Parameters&
                                     _asaParameters){
        asaParameters = _asaParameters;
        return asaParameters;
    }

    void run(void* init,
             size_t byteSize,
             size_t _n_tries = 1000);
    void run(std::vector<double>& init,
             size_t _n_tries = 1000);

    template<size_t N>
    void run(std::array<double,N>& init,
             size_t _n_tries = 1000){
        asaParameters.n_tries = _n_tries;
        void* data = static_cast<void*>(init.data());
        gsl_siman_solve(randomNumberGenerator, data,
                        energy, step,
                        measure, print,
                        NULL, NULL, NULL,
                        sizeof(double)*N, asaParameters);
    }

    template<class State>
    void run(State& init, unsigned _size){
        void* data = static_cast<void*>(&init);
        gsl_siman_solve(randomNumberGenerator, data,
                        energy, step,
                        measure, print,
                        NULL, NULL, NULL,
                        _size, asaParameters);
    }

    void set_energy (double (*)(void*));
    void set_measure(double (*)(void*,
                                void*));
    void set_step   (void   (*)(const gsl_rng*,
                                      void*,
                                      double));
    void set_print  (void   (*)(void*));
}; // end of class SimulatedAnnealing

template<typename T>
T myself(T t){
    return t;
}

template<typename R,typename ...T>
R pass(T... args __attribute__((unused))){
	if constexpr(std::is_void<R>::value) return;
	if constexpr(std::is_arithmetic<R>::value){
		if constexpr(sizeof...(args)==2) return R(1);
		else return R(0);
	}
	return R();
}

} // end of namespace gsl

#endif

max_stars_repo_name: adujovic/JorG

max_issues_repo_name: adujovic/JorG

max_forks_repo_name: adujovic/JorG


text:
/*
 * Copyright (c) 2011-2020 The University of Tennessee and The University
 *                         of Tennessee Research Foundation.  All rights
 *                         reserved.
 * Copyright (c) 2013      Inria. All rights reserved.
 *
 * @precisions normal z -> c d s
 *
 */

#include <lapacke.h>
#include "dplasma.h"
#include "dplasma/types.h"
#include "dplasmaaux.h"

static int
dplasma_zlaset_operator( parsec_execution_stream_t *es,
                         const parsec_tiled_matrix_t *descA,
                         void *_A,
                         dplasma_enum_t uplo, int m, int n,
                         void *args )
{
    int tempmm, tempnn, ldam;
    dplasma_complex64_t *alpha = (dplasma_complex64_t*)args;
    dplasma_complex64_t *A = (dplasma_complex64_t*)_A;
    (void)es;

    tempmm = ((m)==((descA->mt)-1)) ? ((descA->m)-(m*(descA->mb))) : (descA->mb);
    tempnn = ((n)==((descA->nt)-1)) ? ((descA->n)-(n*(descA->nb))) : (descA->nb);
    ldam = BLKLDD( descA, m );

    if (m == n) {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR, dplasma_lapack_const( uplo ), tempmm, tempnn,
            alpha[0], alpha[1], A, ldam);
    } else {
        LAPACKE_zlaset_work(
            LAPACK_COL_MAJOR, 'A', tempmm, tempnn,
            alpha[0], alpha[0], A, ldam);
    }
    return 0;
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlaset_New - Generates the taskpool that set the elements of the matrix
 * A on the diagonal to beta and the off-diagonals eklements to alpha.
 *
 * See parsec_apply_New() for further information.
 *
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is set:
 *          = dplasmaUpperLower: All matrix is referenced.
 *          = dplasmaUpper:      Only upper part is referenced.
 *          = dplasmaLower:      Only lower part is referenced.
 *
 * @param[in] alpha
 *         The constant to which the off-diagonal elements are to be set.
 *
 * @param[in] beta
 *         The constant to which the diagonal elements are to be set.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A. Any tiled matrix
 *          descriptor can be used.
 *          On exit, A has been set accordingly.
 *
 *******************************************************************************
 *
 * @return
 *          \retval NULL if incorrect parameters are given.
 *          \retval The parsec taskpool describing the operation that can be
 *          enqueued in the runtime with parsec_context_add_taskpool(). It, then, needs to be
 *          destroy with dplasma_zlaset_Destruct();
 *
 *******************************************************************************
 *
 * @sa dplasma_zlaset
 * @sa dplasma_zlaset_Destruct
 * @sa dplasma_claset_New
 * @sa dplasma_dlaset_New
 * @sa dplasma_slaset_New
 *
 ******************************************************************************/
parsec_taskpool_t*
dplasma_zlaset_New( dplasma_enum_t uplo,
                    dplasma_complex64_t alpha,
                    dplasma_complex64_t beta,
                    parsec_tiled_matrix_t *A )
{
    dplasma_complex64_t *params = (dplasma_complex64_t*)malloc(2 * sizeof(dplasma_complex64_t));

    params[0] = alpha;
    params[1] = beta;

    return parsec_apply_New( uplo, A, dplasma_zlaset_operator, params );
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 *  dplasma_zlaset_Destruct - Free the data structure associated to an taskpool
 *  created with dplasma_zlaset_New().
 *
 *******************************************************************************
 *
 * @param[in,out] taskpool
 *          On entry, the taskpool to destroy.
 *          On exit, the taskpool cannot be used anymore.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlaset_New
 * @sa dplasma_zlaset
 *
 ******************************************************************************/
void
dplasma_zlaset_Destruct( parsec_taskpool_t *tp )
{
    parsec_apply_Destruct(tp);
}

/**
 *******************************************************************************
 *
 * @ingroup dplasma_complex64
 *
 * dplasma_zlaset - Set the elements of the matrix
 * A on the diagonal to beta and the off-diagonals eklements to alpha.
 *
 *
 *******************************************************************************
 *
 * @param[in,out] parsec
 *          The parsec context of the application that will run the operation.
 *
 * @param[in] uplo
 *          Specifies which part of matrix A is set:
 *          = dplasmaUpperLower: All matrix is referenced.
 *          = dplasmaUpper:      Only upper part is refrenced.
 *          = dplasmaLower:      Only lower part is referenced.
 *
 * @param[in] alpha
 *         The constant to which the off-diagonal elements are to be set.
 *
 * @param[in] beta
 *         The constant to which the diagonal elements are to be set.
 *
 * @param[in,out] A
 *          Descriptor of the distributed matrix A. Any tiled matrix
 *          descriptor can be used.
 *          On exit, A has been set accordingly.
 *
 *******************************************************************************
 *
 * @return
 *          \retval -i if the ith parameters is incorrect.
 *          \retval 0 on success.
 *
 *******************************************************************************
 *
 * @sa dplasma_zlaset_New
 * @sa dplasma_zlaset_Destruct
 * @sa dplasma_claset
 * @sa dplasma_dlaset
 * @sa dplasma_slaset
 *
 ******************************************************************************/
int
dplasma_zlaset( parsec_context_t *parsec,
                dplasma_enum_t uplo,
                dplasma_complex64_t alpha,
                dplasma_complex64_t beta,
                parsec_tiled_matrix_t *A )
{
    parsec_taskpool_t *parsec_zlaset = NULL;

    /* Check input arguments */
    if ((uplo != dplasmaLower) &&
        (uplo != dplasmaUpper) &&
        (uplo != dplasmaUpperLower))
    {
        dplasma_error("dplasma_zlaset", "illegal value of type");
        return -1;
    }

    parsec_zlaset = dplasma_zlaset_New(uplo, alpha, beta, A);

    if ( parsec_zlaset != NULL ) {
        parsec_context_add_taskpool(parsec, (parsec_taskpool_t*)parsec_zlaset);
        dplasma_wait_until_completion(parsec);
        dplasma_zlaset_Destruct( parsec_zlaset );
    }
    return 0;
}

max_stars_repo_name: therault/dplasma

max_issues_repo_name: therault/dplasma

max_forks_repo_name: therault/dplasma


text:
#ifndef PETSC4PY_COMPAT_H
#define PETSC4PY_COMPAT_H

#include <petsc.h>
#include "compat/mpi.h"
#include "compat/hdf5.h"
#include "compat/tao.h"

#endif/*PETSC4PY_COMPAT_H*/


text:
/*
ODE: a program to get optime Runge-Kutta and multi-steps methods.

Copyright 2011-2019, Javier Burguete Tolosa.

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

	1. Redistributions of source code must retain the above copyright notice,
		this list of conditions and the following disclaimer.

	2. Redistributions in binary form must reproduce the above copyright notice,
		this list of conditions and the following disclaimer in the
		documentation and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT
SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file rk_5_4.c
 * \brief Source file to optimize Runge-Kutta 5 steps 4th order methods.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2011-2019.
 */
#define _GNU_SOURCE
#include <string.h>
#include <math.h>
#include <libxml/parser.h>
#include <glib.h>
#include <libintl.h>
#include <gsl/gsl_rng.h>
#include "config.h"
#include "utils.h"
#include "optimize.h"
#include "rk.h"
#include "rk_5_4.h"

#define DEBUG_RK_5_4 0          ///< macro to debug.

/**
 * Function to obtain the coefficients of a 5 steps 4th order Runge-Kutta
 * method.
 */
int
rk_tb_5_4 (Optimize * optimize) ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t5 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  b21 (tb) = r[2];
  t3 (tb) = r[3];
  b31 (tb) = r[4];
  b32 (tb) = r[5];
  t4 (tb) = r[6];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = 1.L / 3.L;
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = 0.25L;
  A[3] = D[3] = 0.L;
  B[3] = b21 (tb) * t1 (tb) * (t2 (tb) - t4 (tb));
  C[3] = (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb)) * (t3 (tb) - t4 (tb));
  E[3] = 0.125L - 1.L / 6.L * t4 (tb);
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b54 (tb) = E[3];
  b53 (tb) = E[2];
  b52 (tb) = E[1];
  b51 (tb) = E[0];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 1.L / 6.L - b52 (tb) * b21 (tb) * t1 (tb)
    - b53 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb));
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 12.L - b52 (tb) * b21 (tb) * sqr (t1 (tb))
    - b53 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = 1.L / 24.L - b53 (tb) * b32 (tb) * b21 (tb) * t1 (tb);
  solve_3 (A, B, C, D);
  b43 (tb) = D[2] / b54 (tb);
  if (isnan (b43 (tb)))
    return 0;
  b42 (tb) = D[1] / b54 (tb);
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = D[0] / b54 (tb);
  if (isnan (b41 (tb)))
    return 0;
  rk_b_5 (tb);
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 4th order, 5th order in
 * equations depending only on time, Runge-Kutta method.
 */
int
rk_tb_5_4t (Optimize * optimize)        ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4t: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t5 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  t3 (tb) = r[2];
  t4 (tb) = r[3];
  b31 (tb) = r[4];
  b21 (tb) = r[5];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = 1.L / 3.L;
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = 0.25L;
  A[3] = A[2] * t1 (tb);
  B[3] = B[2] * t2 (tb);
  C[3] = C[2] * t3 (tb);
  D[3] = D[2] * t4 (tb);
  E[3] = 0.2L;
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b54 (tb) = E[3];
  b53 (tb) = E[2];
  b52 (tb) = E[1];
  b51 (tb) = E[0];
  b32 (tb) = (1.L / 6.L * t4 (tb) - 0.125L
              - t1 (tb) * (b52 (tb) * b21 (tb) * (t4 (tb) - t2 (tb))
                           + b53 (tb) * b31 (tb) * (t4 (tb) - t3 (tb))))
    / (b53 (tb) * t2 (tb) * (t4 (tb) - t3 (tb)));
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 1.L / 6.L - b52 (tb) * b21 (tb) * t1 (tb)
    - b53 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb));
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 12.L - b52 (tb) * b21 (tb) * sqr (t1 (tb))
    - b53 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = 1.L / 24.L - b53 (tb) * b32 (tb) * b21 (tb) * t1 (tb);
  solve_3 (A, B, C, D);
  b43 (tb) = D[2] / b54 (tb);
  if (isnan (b43 (tb)))
    return 0;
  b42 (tb) = D[1] / b54 (tb);
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = D[0] / b54 (tb);
  if (isnan (b41 (tb)))
    return 0;
  rk_b_5 (tb);
#if DEBUG_RK_5_4
  rk_print_tb (optimize, "rk_tb_5_4t", stderr);
  fprintf (stderr, "rk_tb_5_4t: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 3rd-4th order Runge-Kutta
 * pair.
 */
int
rk_tb_5_4p (Optimize * optimize)        ///< Optimize struct.
{
  long double A[3], B[3], C[3], D[3];
  long double *tb;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4p: start\n");
#endif
  if (!rk_tb_5_4 (optimize))
    return 0;
  tb = optimize->coefficient;
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 3.L;
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = 1.L / 6.L;
  solve_3 (A, B, C, D);
  if (isnan (D[0]) || isnan (D[1]) || isnan (D[2]))
    return 0;
  e53 (tb) = D[2];
  e52 (tb) = D[1];
  e51 (tb) = D[0];
  rk_e_5 (tb);
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4p: end\n");
#endif
  return 1;
}

/**
 * Function to obtain the coefficients of a 5 steps 3th-4th order, 4th-5th order
 * in equations depending only on time, Runge-Kutta pair.
 */
int
rk_tb_5_4tp (Optimize * optimize)       ///< Optimize struct.
{
  long double A[4], B[4], C[4], D[4], E[4];
  long double *tb, *r;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_tb_5_4tp: start\n");
#endif
  tb = optimize->coefficient;
  r = optimize->random_data;
  t5 (tb) = 1.L;
  t1 (tb) = r[0];
  t2 (tb) = r[1];
  t3 (tb) = r[2];
  t4 (tb) = r[3];
  b31 (tb) = r[4];
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = t4 (tb);
  E[0] = 0.5L;
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = D[0] * t4 (tb);
  E[1] = 1.L / 3.L;
  A[2] = A[1] * t1 (tb);
  B[2] = B[1] * t2 (tb);
  C[2] = C[1] * t3 (tb);
  D[2] = D[1] * t4 (tb);
  E[2] = 0.25L;
  A[3] = A[2] * t1 (tb);
  B[3] = B[2] * t2 (tb);
  C[3] = C[2] * t3 (tb);
  D[3] = D[2] * t4 (tb);
  E[3] = 0.2L;
  solve_4 (A, B, C, D, E);
  if (isnan (E[0]) || isnan (E[1]) || isnan (E[2]) || isnan (E[3]))
    return 0;
  b54 (tb) = E[3];
  b53 (tb) = E[2];
  b52 (tb) = E[1];
  b51 (tb) = E[0];
  e53 (tb) = (0.25L - 1.L / 3.L * t1 (tb)
              - (1.L / 3.L - 0.5L * t1 (tb)) * t2 (tb))
    / (t3 (tb) * (t3 (tb) - t2 (tb)) * (t3 (tb) - t1 (tb)));
  if (isnan (e53 (tb)))
    return 0;
  e52 (tb) = (1.L / 3.L - 0.5L * t1 (tb)
              - t3 (tb) * (t3 (tb) - t1 (tb)) * e53 (tb))
    / (t2 (tb) * (t2 (tb) - t1 (tb)));
  if (isnan (e52 (tb)))
    return 0;
  e51 (tb) = (0.5L - t2 (tb) * e52 (tb) - t3 (tb) * e53 (tb)) / t1 (tb);
  if (isnan (e51 (tb)))
    return 0;
  b21 (tb) = (1.L / 6.L * b53 (tb) * (t4 (tb) - t3 (tb))
              + e53 (tb) * (0.125L - 1.L / 6.L * t4 (tb)))
    / (t1 (tb) * (e52 (tb) * b53 (tb) * (t4 (tb) - t3 (tb))
                  - e53 (tb) * b52 (tb) * (t4 (tb) - t2 (tb))));
  if (isnan (b21 (tb)))
    return 0;
  b32 (tb) = (1.L / 6.L * t4 (tb) - 0.125L
              - t1 (tb) * (b52 (tb) * b21 (tb) * (t4 (tb) - t2 (tb))
                           + b53 (tb) * b31 (tb) * (t4 (tb) - t3 (tb))))
    / (b53 (tb) * t2 (tb) * (t4 (tb) - t3 (tb)));
  if (isnan (b32 (tb)))
    return 0;
  A[0] = t1 (tb);
  B[0] = t2 (tb);
  C[0] = t3 (tb);
  D[0] = 1.L / 6.L - b52 (tb) * b21 (tb) * t1 (tb)
    - b53 (tb) * (b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb));
  A[1] = A[0] * t1 (tb);
  B[1] = B[0] * t2 (tb);
  C[1] = C[0] * t3 (tb);
  D[1] = 1.L / 12.L - b52 (tb) * b21 (tb) * sqr (t1 (tb))
    - b53 (tb) * (b31 (tb) * sqr (t1 (tb)) + b32 (tb) * sqr (t2 (tb)));
  A[2] = 0.L;
  B[2] = b21 (tb) * t1 (tb);
  C[2] = b31 (tb) * t1 (tb) + b32 (tb) * t2 (tb);
  D[2] = 1.L / 24.L - b53 (tb) * b32 (tb) * b21 (tb) * t1 (tb);
  solve_3 (A, B, C, D);
  b43 (tb) = D[2] / b54 (tb);
  if (isnan (b43 (tb)))
    return 0;
  b42 (tb) = D[1] / b54 (tb);
  if (isnan (b42 (tb)))
    return 0;
  b41 (tb) = D[0] / b54 (tb);
  if (isnan (b41 (tb)))
    return 0;
  rk_b_5 (tb);
  rk_e_5 (tb);
#if DEBUG_RK_5_4
  rk_print_tb (optimize, "rk_tb_5_4tp", stderr);
  fprintf (stderr, "rk_tb_5_4tp: end\n");
#endif
  return 1;
}

/**
 * Function to calculate the objective function of a 5 steps 4th order
 * Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_4 (RK * rk)   ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b54 (tb) < 0.L)
    o += b54 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_4: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 5 steps 4th order, 5th
 * order in equations depending only on time, Runge-Kutta method.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_4t (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4t: start\n");
#endif
  tb = rk->tb->coefficient;
#if DEBUG_RK_5_4
  rk_print_tb (optimize, "rk_objective_tb_5_4t", stderr);
#endif
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b54 (tb) < 0.L)
    o += b54 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4t: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_4t: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 5 steps 3rd-4th order
 * Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_4p (RK * rk)  ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4p: start\n");
#endif
  tb = rk->tb->coefficient;
  o = fminl (0.L, b20 (tb));
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b54 (tb) < 0.L)
    o += b54 (tb);
  if (e50 (tb) < 0.L)
    o += e50 (tb);
  if (e51 (tb) < 0.L)
    o += e51 (tb);
  if (e52 (tb) < 0.L)
    o += e52 (tb);
  if (e53 (tb) < 0.L)
    o += e53 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4p: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_4p: end\n");
#endif
  return o;
}

/**
 * Function to calculate the objective function of a 5 steps 3th-4th order,
 * 4th-5th order in equations depending only on time, Runge-Kutta pair.
 *
 * \return objective function value.
 */
long double
rk_objective_tb_5_4tp (RK * rk) ///< RK struct.
{
  long double *tb;
  long double o;
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4t: start\n");
#endif
  tb = rk->tb->coefficient;
#if DEBUG_RK_5_4
  rk_print_tb (optimize, "rk_objective_tb_5_4tp", stderr);
#endif
  o = fminl (0.L, b20 (tb));
  if (b21 (tb) < 0.L)
    o += b21 (tb);
  if (b30 (tb) < 0.L)
    o += b30 (tb);
  if (b32 (tb) < 0.L)
    o += b32 (tb);
  if (b40 (tb) < 0.L)
    o += b40 (tb);
  if (b41 (tb) < 0.L)
    o += b41 (tb);
  if (b42 (tb) < 0.L)
    o += b42 (tb);
  if (b43 (tb) < 0.L)
    o += b43 (tb);
  if (b50 (tb) < 0.L)
    o += b50 (tb);
  if (b51 (tb) < 0.L)
    o += b51 (tb);
  if (b52 (tb) < 0.L)
    o += b52 (tb);
  if (b53 (tb) < 0.L)
    o += b53 (tb);
  if (b54 (tb) < 0.L)
    o += b54 (tb);
  if (e50 (tb) < 0.L)
    o += e50 (tb);
  if (e51 (tb) < 0.L)
    o += e51 (tb);
  if (e52 (tb) < 0.L)
    o += e52 (tb);
  if (e53 (tb) < 0.L)
    o += e53 (tb);
  if (o < 0.L)
    {
      o = 40.L - o;
      goto end;
    }
  o = 30.L
    + fmaxl (1.L, fmaxl (t1 (tb), fmaxl (t2 (tb), fmaxl (t3 (tb), t4 (tb)))));
  if (rk->strong)
    {
      rk_bucle_ac (rk);
      o = fminl (o, *rk->ac0->optimal);
    }
end:
#if DEBUG_RK_5_4
  fprintf (stderr, "rk_objective_tb_5_4tp: optimal=%Lg\n", o);
  fprintf (stderr, "rk_objective_tb_5_4tp: end\n");
#endif
  return o;
}

max_stars_repo_name: jburguete/ode

max_issues_repo_name: jburguete/ode

max_forks_repo_name: jburguete/ode


text:
///
/// @file
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <starneig_test_config.h>
#include <starneig/configuration.h>
#include "solvers.h"
#include "../common/common.h"
#include "../common/parse.h"
#include "../common/local_pencil.h"
#ifdef STARNEIG_ENABLE_MPI
#include "../common/starneig_pencil.h"
#endif
#include "../common/threads.h"
#include <starneig/starneig.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <cblas.h>
#include <omp.h>

#ifdef MAGMA_FOUND
#include <cuda_runtime_api.h>
#include <cuda.h>
#include <magma_auxiliary.h>
#include <magma_d.h>
#endif

static hook_solver_state_t lapack_prepare(
    int argc, char * const *argv, struct hook_data_env *env)
{
    return (hook_solver_state_t) env->data;
}

static int lapack_finalize(hook_solver_state_t state, struct hook_data_env *env)
{
    return 0;
}

static int lapack_run(hook_solver_state_t state)
{
    pencil_t data = (pencil_t) state;

    extern void dgehrd_(
        int const *,    // the order of the matrix A
        int const *,    // left bound
        int const *,    // right bound
        double *,       // input/output matrix
        int const *,    // input/output matrix leading dimension
        double *,       // tau (scalar factors)
        double *,       // work space
        int const *,    // work space size
        int *);         // info

    extern void dormhr_(
        char const *,   // side
        char const *,   // transpose
        int const *,    // row count
        int const *,    // column count
        int const *,    // left bound
        int const *,    // right bound
        double const *, // elementary reflectors
        int const *,    // elementary reflector leading dimension
        double const *, // scalar factors (tau)
        double *,       // input/output matrix
        int const *,    // input/output matrix leading dimension
        double *,       // work space
        int const *,    // work space size
        int *);         // info

    extern void dgeqrf_(int const *, int const *, double *, int const *,
        double *, double *, int const *, int *);

    extern void dormqr_(char const *, char const *, int const *, int const *,
        int const *, double const *, int const *, double const *, double *,
        int const *, double*, const int *, int *);

    extern void dgghd3_(char const *, char const *, int const *, int const *,
        int const *, double *, int const *, double *, int const *, double *,
        int const *, double *, int const *, double *, int const *, int *);

    int n = LOCAL_MATRIX_N(data->mat_a);
    double *A = LOCAL_MATRIX_PTR(data->mat_a);
    int ldA = LOCAL_MATRIX_LD(data->mat_a);
    double *Q = LOCAL_MATRIX_PTR(data->mat_q);
    int ldQ = LOCAL_MATRIX_LD(data->mat_q);

    double *B = NULL;
    int ldB = 0;
    double *Z = NULL;
    int ldZ = 0;

    double *tau = NULL;
    double *work = NULL;
    int info, ilo = 1, ihi = n;

    threads_set_mode(THREADS_MODE_LAPACK);

    if (data->mat_b != NULL) {
        B = LOCAL_MATRIX_PTR(data->mat_b);
        ldB = LOCAL_MATRIX_LD(data->mat_b);
        Z = LOCAL_MATRIX_PTR(data->mat_z);
        ldZ = LOCAL_MATRIX_LD(data->mat_z);

        //
        // allocate workspace
        //

        int lwork = 0;

        {
            int _lwork = -1;
            double dlwork;

            dgeqrf_(&n, &n, B, &ldB, tau, &dlwork, &_lwork, &info);
            if (info != 0)
                goto cleanup;

            lwork = MAX(lwork, dlwork);
        }

        {
            int _lwork = -1;
            double dlwork;

            dormqr_("L", "T", &n, &n, &n,
                B, &ldB, tau, A, &ldA, &dlwork, &_lwork, &info);
            if (info != 0)
                goto cleanup;

            lwork = MAX(lwork, dlwork);
        }

        {
            int _lwork = -1;
            double dlwork;

            dormqr_("R", "N", &n, &n, &n,
                B, &ldB, tau, Q, &ldQ, &dlwork, &_lwork, &info);
            if (info != 0)
                goto cleanup;

            lwork = MAX(lwork, dlwork);
        }

        {
            int _lwork = -1;
            double dlwork;

            dgghd3_("V", "V", &n, &ilo, &ihi,
                A, &ldA, B, &ldB, Q, &ldQ, Z, &ldZ, &dlwork, &_lwork, &info);
            if (info != 0)
                goto cleanup;

            lwork = MAX(lwork, dlwork);
        }

        tau = malloc(n*sizeof(double));
        work = malloc(lwork*sizeof(double));

        //
        // reduce
        //

        // form B = ~Q * R
        dgeqrf_(&n, &n, B, &ldB, tau, work, &lwork, &info);
        if (info != 0)
            goto cleanup;

        // A <- ~Q^T * A
        dormqr_("L", "T", &n, &n, &n, B, &ldB, tau, A, &ldA, work, &lwork,
            &info);
        if (info != 0)
            goto cleanup;

        // Q <- Q * ~Q
        dormqr_("R", "N", &n, &n, &n, B, &ldB, tau, Q, &ldQ, work, &lwork,
            &info);
        if (info != 0)
            goto cleanup;

        // clean B (B <- R)
        for (int i = 0; i < n; i++)
            for (int j = i+1; j < n; j++)
                B[(size_t)i*ldB+j] = 0.0;

        // reduce (A,B) to Hessenberg-triangular form
        dgghd3_("V", "V", &n, &ilo, &ihi,
            A, &ldA, B, &ldB, Q, &ldQ, Z, &ldZ, work, &lwork, &info);
        if (info != 0)
            goto cleanup;
    }
    else {

        int lwork = -1;
        double dlwork;

        // request optimal work space size
        dgehrd_(&n, &ilo, &ihi, A, &ldA,
            tau, &dlwork, &lwork, &info);

        if (info != 0)
            goto cleanup;

        lwork = dlwork;
        work = malloc(lwork*sizeof(double));
        tau = malloc(n*sizeof(double));

        // reduce
        dgehrd_(&n, &ilo, &ihi, A, &ldA,
            tau, work, &lwork, &info);

        if (info != 0)
            goto cleanup;

        free(work);
        work = NULL;

        // request optimal work space size
        lwork = -1;
        dormhr_("Right", "No transpose", &n, &n, &ilo, &ihi, A, &ldA, tau,
            Q, &ldQ, &dlwork, &lwork, &info);

        if (info != 0)
            goto cleanup;

        lwork = dlwork;
        work = malloc(lwork*sizeof(double));

        // form Q
        dormhr_("Right", "No transpose", &n, &n, &ilo, &ihi, A, &ldA, tau,
            Q, &ldQ, work, &lwork, &info);

        if (info != 0)
            goto cleanup;

        for (int i = 0; i < n; i++)
            for (int j = i+2; j < n; j++)
                A[i*ldA+j] = 0.0;
    }

cleanup:

    threads_set_mode(THREADS_MODE_DEFAULT);

    free(work);
    free(tau);
    return info;
}

const struct hook_solver hessenberg_lapack_solver = {
    .name = "lapack",
    .desc = "LAPACK's dgehrd/dgghrd subroutine",
    .formats = (hook_data_format_t[]) { HOOK_DATA_FORMAT_PENCIL_LOCAL, 0 },
    .prepare = &lapack_prepare,
    .finalize = &lapack_finalize,
    .run = &lapack_run
};

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

#if defined(PDGEHRD_FOUND) && defined(PDORMHR_FOUND) && defined(PDLASET_FOUND)

static hook_solver_state_t scalapack_prepare(
    int argc, char * const *argv, struct hook_data_env *env)
{
    return env;
}

static int scalapack_finalize(
    hook_solver_state_t state, struct hook_data_env *env)
{
    return 0;
}

static int has_valid_descr(
    int matrix_size, int section_size, const starneig_blacs_descr_t *descr)
{
    if (descr->m != matrix_size || descr->n != matrix_size)
        return 0;
    if (descr->sm != section_size || descr->sn != section_size)
        return 0;
    return 1;
}

static int scalapack_run(hook_solver_state_t state)
{
    extern void pdgehrd_(int const *, int const *, int const *, double *,
        int const *, int const *, const starneig_blacs_descr_t *, double *,
        double *, int const *, int *);

    extern void pdormhr_(char const *, char const *, int const *, int const *,
        int const *, int const *, double *, int const *, int const *,
        const starneig_blacs_descr_t *, double *, double *, int const *,
        int const *, const starneig_blacs_descr_t *, double *, int const *,
        int *);

    extern void pdlaset_(char const *, int const *, int const *,
        double const *, double const *, double *, int const *, int const *,
        const starneig_blacs_descr_t *);

    threads_set_mode(THREADS_MODE_SCALAPACK);

    struct hook_data_env *env = state;
    pencil_t pencil = (pencil_t) env->data;

    if (pencil->mat_a == NULL) {
        fprintf(stderr, "Missing matrix A.\n");
        return -1;
    }

    if (pencil->mat_b != NULL) {
        fprintf(stderr, "Solver does not support generalized cases.\n");
        return -1;
    }

    int n = STARNEIG_MATRIX_N(pencil->mat_a);
    int sn = STARNEIG_MATRIX_BN(pencil->mat_a);

    starneig_distr_t distr = STARNEIG_MATRIX_DISTR(pencil->mat_a);
    starneig_blacs_context_t context = starneig_distr_to_blacs_context(distr);

    starneig_blacs_descr_t desc_a, desc_q;
    double *local_a, *local_q;
    STARNEIG_BLACS_MATRIX_DESCR_LOCAL(
        pencil->mat_a, context, &desc_a, (void **)&local_a);
    STARNEIG_BLACS_MATRIX_DESCR_LOCAL(
        pencil->mat_q, context, &desc_q, (void **)&local_q);

    if (!has_valid_descr(n, sn, &desc_a)) {
        fprintf(stderr, "Matrix A has invalid dimensions.\n");
        return -1;
    }

    if (pencil->mat_q != NULL && !has_valid_descr(n, sn, &desc_q)) {
        fprintf(stderr, "Matrix Q has invalid dimension.\n");
        return -1;
    }

    int ilo = 1, ihi = n, ia = 1, ja = 1, ic = 1, jc = 1, lwork, info;
    double *work = NULL, *tau = NULL, _work;

    lwork = -1;
    pdgehrd_(&n, &ilo, &ihi, NULL, &ia, &ja, &desc_a, NULL,
        &_work, &lwork, &info);

    if (info)
        goto cleanup;

    lwork = _work;
    work = malloc(lwork*sizeof(double));
    tau = malloc(n*sizeof(double));

    pdgehrd_(&n, &ilo, &ihi, local_a, &ia, &ja, &desc_a, tau,
        work, &lwork, &info);

    if (info)
        goto cleanup;

    lwork = -1;
    pdormhr_("Right", "No transpose", &n, &n, &ilo, &ihi, NULL,
        &ia, &ja, &desc_a, NULL, NULL, &ic, &jc,
        &desc_q, &_work, &lwork, &info);

    if (info)
        goto cleanup;

    free(work);
    lwork = _work;
    work = malloc(lwork*sizeof(double));

    pdormhr_("Right", "No transpose", &n, &n, &ilo, &ihi, local_a,
        &ia, &ja, &desc_a, tau, local_q, &ic, &jc,
        &desc_q, work, &lwork, &info);

    {
        int nm1 = n-2, one = 1, three = 3;
        double dzero = 0.0;
        pdlaset_("Lower", &nm1, &nm1, &dzero, &dzero, local_a, &three,
            &one, &desc_a);
    }

cleanup:

    threads_set_mode(THREADS_MODE_DEFAULT);
    starneig_blacs_gridexit(context);

    free(work);
    free(tau);

    return info;
}

const struct hook_solver hessenberg_scalapack_solver = {
    .name = "scalapack",
    .desc = "pdgehrd subroutine from scaLAPACK",
    .formats = (hook_data_format_t[]) {
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .prepare = &scalapack_prepare,
    .finalize = &scalapack_finalize,
    .run = &scalapack_run
};

#endif // PDGEHRD_FOUND

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

struct starpu_state {
    int argc;
    char * const *argv;
    struct hook_data_env *env;
};

static void starpu_print_usage(int argc, char * const *argv)
{
    printf(
        "  --cores [default,(num)} -- Number of CPU cores\n"
        "  --gpus [default,(num)} -- Number of GPUS\n"
        "  --tile-size [default,(num)] -- tile size\n"
        "  --panel-width [default,(num)] -- Panel width\n"
    );
}

static int starpu_check_args(int argc, char * const *argv, int *argr)
{
    struct multiarg_t arg_cores = read_multiarg(
        "--cores", argc, argv, argr, "default", NULL);
    struct multiarg_t arg_gpus = read_multiarg(
        "--gpus", argc, argv, argr, "default", NULL);
    struct multiarg_t tile_size = read_multiarg(
        "--tile-size", argc, argv, argr, "default", NULL);
    struct multiarg_t panel_width = read_multiarg(
        "--panel-width", argc, argv, argr, "default", NULL);

    if (arg_cores.type == MULTIARG_INVALID)
        return -1;

    if (arg_gpus.type == MULTIARG_INVALID)
        return -1;

    if (tile_size.type == MULTIARG_INVALID ||
    (tile_size.type == MULTIARG_INT && tile_size.int_value < 1)) {
        fprintf(stderr, "Invalid tile size.\n");
        return -1;
    }

    if (panel_width.type == MULTIARG_INVALID ||
    (panel_width.type == MULTIARG_INT && panel_width.int_value < 1)) {
        fprintf(stderr, "Invalid panel width.\n");
        return -1;
    }

    return 0;
}

static void starpu_print_args(int argc, char * const *argv)
{
    print_multiarg("--cores", argc, argv, "default", NULL);
    print_multiarg("--gpus", argc, argv, "default", NULL);
    print_multiarg("--tile-size", argc, argv, "default", NULL);
    print_multiarg("--panel-width", argc, argv, "default", NULL);
}

static hook_solver_state_t starpu_prepare(
    int argc, char * const *argv, struct hook_data_env *env)
{
    struct starpu_state *state = malloc(sizeof(struct starpu_state));

    state->argc = argc;
    state->argv = argv;
    state->env = env;

    struct multiarg_t arg_cores = read_multiarg(
        "--cores", argc, argv, NULL, "default", NULL);
    struct multiarg_t arg_gpus = read_multiarg(
        "--gpus", argc, argv, NULL, "default", NULL);

    int cores = STARNEIG_USE_ALL;
    if (arg_cores.type == MULTIARG_INT)
        cores = arg_cores.int_value;

    int gpus = STARNEIG_USE_ALL;
    if (arg_gpus.type == MULTIARG_INT)
        gpus = arg_gpus.int_value;

#ifdef STARNEIG_ENABLE_MPI
    if (env->format == HOOK_DATA_FORMAT_PENCIL_STARNEIG ||
    env->format == HOOK_DATA_FORMAT_PENCIL_BLACS)
        starneig_node_init(cores, gpus, STARNEIG_FAST_DM);
    else
#endif
        starneig_node_init(
            cores, gpus, STARNEIG_HINT_SM | STARNEIG_AWAKE_WORKERS);

    return state;
}

static int starpu_finalize(hook_solver_state_t state, struct hook_data_env *env)
{
    if (state == NULL)
        return 0;

    starneig_node_finalize();

    free(state);
    return 0;
}

static int starpu_run(hook_solver_state_t state)
{
    int argc = ((struct starpu_state *) state)->argc;
    char * const *argv = ((struct starpu_state *) state)->argv;
    struct hook_data_env *env = ((struct starpu_state *) state)->env;

    struct starneig_hessenberg_conf conf;
    starneig_hessenberg_init_conf(&conf);

    struct multiarg_t tile_size = read_multiarg(
        "--tile-size", argc, argv, NULL, "default", NULL);
    struct multiarg_t panel_width = read_multiarg(
        "--panel-width", argc, argv, NULL, "default", NULL);

    if (tile_size.type == MULTIARG_INT)
        conf.tile_size = tile_size.int_value;
    if (panel_width.type == MULTIARG_INT)
        conf.panel_width = panel_width.int_value;

    int ret = 0;

    if (env->format == HOOK_DATA_FORMAT_PENCIL_LOCAL) {
        pencil_t pencil = (pencil_t) env->data;
        if (pencil->mat_b != NULL) {
            ret = starneig_GEP_SM_HessenbergTriangular(
                LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_b), LOCAL_MATRIX_LD(pencil->mat_b),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q),
                LOCAL_MATRIX_PTR(pencil->mat_z), LOCAL_MATRIX_LD(pencil->mat_z)
            );
        }
        else {
            ret = starneig_SEP_SM_Hessenberg_expert(&conf,
                LOCAL_MATRIX_N(pencil->mat_a), 0, LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q)
            );
        }
    }
#ifdef STARNEIG_ENABLE_MPI
    if (env->format == HOOK_DATA_FORMAT_PENCIL_BLACS) {
        pencil_t pencil = (pencil_t) env->data;
        if (pencil->mat_b != NULL) {
#ifdef STARNEIG_GEP_DM_HESSENBERGTRIANGULAR
            ret = starneig_GEP_DM_HessenbergTriangular(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_b),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q),
                STARNEIG_MATRIX_HANDLE(pencil->mat_z));
#else
            fprintf(stderr,
                "Solver does not support generalized cases in distributed "
                "memory.\n");
            return -1;
#endif
        }
        else {
            ret = starneig_SEP_DM_Hessenberg(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q));
        }
    }
#endif

    return ret;
}

const struct hook_solver hessenberg_starpu_solver = {
    .name = "starneig",
    .desc = "StarPU based subroutine",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &starpu_print_usage,
    .print_args = &starpu_print_args,
    .check_args = &starpu_check_args,
    .prepare = &starpu_prepare,
    .finalize = &starpu_finalize,
    .run = &starpu_run
};

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

static void starpu_simple_print_usage(int argc, char * const *argv)
{
    printf(
        "  --cores [default,(num)} -- Number of CPU cores\n"
        "  --gpus [default,(num)} -- Number of GPUS\n"
    );
}

static void starpu_simple_print_args(int argc, char * const *argv)
{
    print_multiarg("--cores", argc, argv, "default", NULL);
    print_multiarg("--gpus", argc, argv, "default", NULL);
}

static int starpu_simple_check_args(int argc, char * const *argv, int *argr)
{
    struct multiarg_t arg_cores = read_multiarg(
        "--cores", argc, argv, argr, "default", NULL);
    struct multiarg_t arg_gpus = read_multiarg(
        "--gpus", argc, argv, argr, "default", NULL);

    if (arg_cores.type == MULTIARG_INVALID)
        return -1;

    if (arg_gpus.type == MULTIARG_INVALID)
        return -1;

    return 0;
}

static hook_solver_state_t starpu_simple_prepare(
    int argc, char * const *argv, struct hook_data_env *env)
{
    struct starpu_state *state = malloc(sizeof(struct starpu_state));

    state->argc = argc;
    state->argv = argv;
    state->env = env;

    struct multiarg_t arg_cores = read_multiarg(
        "--cores", argc, argv, NULL, "default", NULL);
    struct multiarg_t arg_gpus = read_multiarg(
        "--gpus", argc, argv, NULL, "default", NULL);

    int cores = STARNEIG_USE_ALL;
    if (arg_cores.type == MULTIARG_INT)
        cores = arg_cores.int_value;

    int gpus = STARNEIG_USE_ALL;
    if (arg_gpus.type == MULTIARG_INT)
        gpus = arg_gpus.int_value;

#ifdef STARNEIG_ENABLE_MPI
    if (env->format == HOOK_DATA_FORMAT_PENCIL_STARNEIG ||
    env->format == HOOK_DATA_FORMAT_PENCIL_BLACS)
        starneig_node_init(cores, gpus, STARNEIG_FAST_DM);
    else
#endif
        starneig_node_init(
            cores, gpus, STARNEIG_HINT_SM | STARNEIG_AWAKE_WORKERS);

    return state;
}

static int starpu_simple_finalize(
    hook_solver_state_t state, struct hook_data_env *env)
{
    if (state == NULL)
        return 0;

    starneig_node_finalize();

    free(state);
    return 0;
}

static int starpu_simple_run(hook_solver_state_t state)
{
    struct hook_data_env *env = ((struct starpu_state *) state)->env;

    int ret = 0;

    if (env->format == HOOK_DATA_FORMAT_PENCIL_LOCAL) {
        pencil_t pencil = (pencil_t) env->data;

        if (pencil->mat_b != NULL) {
            ret = starneig_GEP_SM_HessenbergTriangular(
                LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_b), LOCAL_MATRIX_LD(pencil->mat_b),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q),
                LOCAL_MATRIX_PTR(pencil->mat_z), LOCAL_MATRIX_LD(pencil->mat_z)
            );
        }
        else {
            ret = starneig_SEP_SM_Hessenberg(LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q)
            );
        }
    }
#ifdef STARNEIG_ENABLE_MPI
    if (env->format == HOOK_DATA_FORMAT_PENCIL_BLACS) {
        pencil_t pencil = (pencil_t) env->data;

        if (pencil->mat_b != NULL) {
#ifdef STARNEIG_GEP_DM_HESSENBERGTRIANGULAR
            ret = starneig_GEP_DM_HessenbergTriangular(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_b),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q),
                STARNEIG_MATRIX_HANDLE(pencil->mat_z));
#else
            fprintf(stderr,
                "Solver does not support distributed memory in generalized "
                "cases.\n");
            return -1;
#endif
        }
        else {
            ret = starneig_SEP_DM_Hessenberg(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q));
        }
    }
#endif

    return ret;
}

const struct hook_solver hessenberg_starpu_simple_solver = {
    .name = "starneig-simple",
    .desc = "StarPU based subroutine (simplified interface)",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &starpu_simple_print_usage,
    .print_args = &starpu_simple_print_args,
    .check_args = &starpu_simple_check_args,
    .prepare = &starpu_simple_prepare,
    .finalize = &starpu_simple_finalize,
    .run = &starpu_simple_run
};

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

#ifdef MAGMA_FOUND

static hook_solver_state_t magma_dgehrd_prepare(
    int argc, char * const *argv, struct hook_data_env *env)
{
    if (magma_init() != MAGMA_SUCCESS)
        return NULL;
    return (hook_solver_state_t) env->data;
}

static int magma_dgehrd_finalize(
    hook_solver_state_t state, struct hook_data_env *env)
{
    magma_finalize();
    return 0;
}

static int magma_dgehrd_run(hook_solver_state_t state)
{
    pencil_t data = (pencil_t) state;

    int n = LOCAL_MATRIX_N(data->mat_a);
    double *A = LOCAL_MATRIX_PTR(data->mat_a);
    double *Q = LOCAL_MATRIX_PTR(data->mat_q);
    int ldA = LOCAL_MATRIX_LD(data->mat_a);
    int ldQ = LOCAL_MATRIX_LD(data->mat_q);

    int info;

    double *tau = NULL;
    double *work = NULL;
    double *dT = NULL;

    threads_set_mode(THREADS_MODE_LAPACK);

    // request optimal work space size
    double _work;
    magma_dgehrd(n, 1, n, A, ldA, tau, &_work, -1, dT, &info);

    if (info != 0)
        goto finalize;

    tau = malloc(LOCAL_MATRIX_N(data->mat_a)*sizeof(double));

    int lwork = _work;
    work = malloc(lwork*sizeof(double));

    int nb = magma_get_dgehrd_nb(LOCAL_MATRIX_N(data->mat_a));
    cudaMalloc((void**)&dT, nb*n*sizeof(double));

    // reduce
    magma_dgehrd(n, 1, n, A, ldA, tau, work, lwork, dT, &info);

    if (info != 0)
        goto finalize;

    free(work);
    work = NULL;

    // copy A -> Q
    for (int i = 0; i < n; i++)
        memcpy(Q+i*ldQ, A+i*ldA, n*sizeof(double));

    // form Q
    magma_dorghr(n, 1, n, Q, ldQ, tau, dT, nb, &info);

    if (info != 0)
        goto finalize;

    // zero entries below the first sub-diagonal
    for (int i = 0; i < n-1; i++)
        memset(A+i*ldA+i+2, 0, (n-i-2)*sizeof(double));

finalize:

    threads_set_mode(THREADS_MODE_DEFAULT);

    cudaFree(dT);
    free(work);
    free(tau);
    return info;
}

const struct hook_solver hessenberg_magma_solver = {
    .name = "magma",
    .desc = "MAGMA's dgehrd subroutine",
    .formats = (hook_data_format_t[]) { HOOK_DATA_FORMAT_PENCIL_LOCAL, 0 },
    .prepare = &magma_dgehrd_prepare,
    .finalize = &magma_dgehrd_finalize,
    .run = &magma_dgehrd_run
};

#endif // MAGMA_FOUND

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
/**
  *
  * For reference: Pearson's rho...
  *
  *               sum( (x_i-E[x])(y_i-E[y]) )
  * rho = -------------------------------------------
  *       sqrt( sum((x_i-E[x])^2) sum((y_i-E[y])^2) )
  */

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <limits.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <alloca.h>

#include <gsl/gsl_cdf.h>

#include "stattest.h"
#include "mix.h"
#include "bvr.h"
#include "limits.h"

struct Pair {
	float cv; // Continuous Variable
	unsigned int dv; // Discrete Variable
};

static int _cmp_pair( const void *pvl, const void *pvr ) {
	const struct Pair *l = (const struct Pair *)pvl;
	const struct Pair *r = (const struct Pair *)pvr;
	const float delta = l->cv - r->cv;
	if( delta == 0.0 )
		return 0.0;
	else
		return  delta < 0.0 ? -1 : +1;
}

struct MixCovars {

	/**
	  * Maximum number of samples
	  */
	unsigned int SAMPLE_CAPACITY;

	/**
	  * All category labels pushed into this accumulator must
	  * be in [0,category_capacity).
	 */
	unsigned int CATEGORY_CAPACITY;

	/**
	  * Count of samples pushed since last cat_clear.
	  */
	unsigned int sample_count;

	/**
	  * This essentially bounds the allowed category LABELS.
	  * Any category label push'ed in must be in [0,expected_categories).
	  * Only this member, NOT category_capacity, should be used to size
	  * counts arrays!
	  */
	unsigned int expected_categories;

	/**
	  * This is an count of DISTINCT LABELS observed by push.
	  */
	unsigned int observed_categories;

	// Calculation state ///////////////////////////////////////////////////

	double mean_rank;  // necessarily of both covariates
	double sum_sq_dev; // used for Kruskal-Wallis and Spearman

	/**
	  * Following struct array is ONLY for Spearman rho calculation.
	  * Note that the number of non-empty categories after data entry
	  * says NOTHING about WHICH categories (indices) are non-zero.
	  * An earlier implementation was assuming category_count[0,1] > 0,
	  * BUT THIS NEED NOT BE TRUE! Thus, all this hack...
	  */
	struct {
		unsigned int index;
		unsigned int count;
		unsigned int meanRank;
	} edge[2];

	double sum_dev_prod; // used only for Spearman-rho

	// Buffer //////////////////////////////////////////////////////////////

	/**
	  * The TOTAL amount of malloc'ed space (for the purposes
	  * of fast clearing with memset).
	  */
	size_t SIZEOF_BUFFERS;

	/**
	  * Two buffers
	  */
	struct Pair *samples;

	/**
	  * A buffer allocated once in construction and reused for
	  * all subsequent analyses.
	  * Reallocation is not currently supported.
	  */
	unsigned int *category_count;
};

#ifdef _UNITTEST_MIX_
unsigned arg_min_mixb_count = 1;
#else
extern unsigned arg_min_mixb_count;
#endif

/**
  * Both minCat and maxCat assume that samples is non-empty,
  * but it's entirely possible that there were NO samples.
  */
static unsigned _minCat( struct MixCovars *co ) {
	unsigned int i = 0;
	while( i < co->expected_categories ) {
		if( co->category_count[i] > 0 )
			return i;
		i++;
	}
	return INT_MAX; // just to silence compiler.
}

static unsigned _maxCat( struct MixCovars *co ) {

	int i = co->expected_categories;
	while( i-- > 0 ) {
		if( co->category_count[i] > 0 )
			return (unsigned int)i;
	}
	return INT_MAX; // just to silence compiler
}

/**
  * Notice that all calcuations are performed without actually requiring
  * an array of ranks.
  */
static unsigned int _rank_sums( struct MixCovars *co,
	double *sums ) {

	const unsigned int N = co->sample_count;
	const bool PRECALC_SPEARMAN_PRODUCT
		= co->observed_categories == 2;
	const unsigned UPPER
		= co->edge[1].index;
	unsigned int ties = 0;
	unsigned int until = 0;
	float rank = 0;
	double diff;

	// The default pair::operator< automatically uses the first element.

	qsort( co->samples, co->sample_count, sizeof(struct Pair), _cmp_pair );

	for(unsigned int i = 0; i < N; i++ ) {

		const unsigned int cat
			= co->samples[i].dv;

		// The rank of the current sample is i+1 UNLESS we're in
		// the midst of a run of samples with the same value in
		// which case the last computed rank (a mean) is reused.

		if( ! ( i < until ) ) {

			until = i + 1;
			while( until < N
					&& ( co->samples[i].cv == co->samples[until].cv ) ) {
				ties++;
				until++;
			}

			if( i == 0 && until == N ) {
				// degenerate case: all ties.
			}

			// Mean of n integers starting on r:
			// [ (r+0) + (r+1) + (r+2) + ... + (r+(n-1)) ] / n
			// ... = \frac{nr + \sum_{i=0}^{n-1} i}{n}
			// ... = r + (1/n)[ (n-1)n / 2 ]
			// ... = r + (n-1)/2

			rank = (1.0+i) + ((until-i)-1)/2.0;

			// ...notice if the very next sample -is- different
			// then rank == i+1, the trivial case.
		}

		/**
		  * Warning: If the numeric feature is constant then rank==meanRank
		  * and diff == 0
		  */
		diff = rank - co->mean_rank;

		co->sum_sq_dev += ( diff * diff );

		if( PRECALC_SPEARMAN_PRODUCT ) {
			co->sum_dev_prod += ( diff * ( co->edge[UPPER==cat?1:0].meanRank - co->mean_rank ) );
		}

		sums[ cat ] += rank;
	}

	return ties;
}

#if defined(_UNITTEST_MIX_)
static void dbg_dump( struct MixCovars *co, FILE *fp ) {
	for(unsigned int i = 0; i < co->sample_count; i++ )
		fprintf( fp, "%.3f\t%d\n", co->samples[i].cv, co->samples[i].dv );
}
#endif

/***************************************************************************
  * Publics
  */

void mix_destroy( void *pv ) {

	if( pv ) {
		struct MixCovars *co = (struct MixCovars *)pv;
		if( co->samples )
			free( co->samples );
		free( pv );
	}
}

/**
  * Pre-allocate a set of working buffers large enough for all anticipated
  * calculations (max feature length) and a struct to wrap them.
  */
void *mix_create( unsigned int sample_capacity, unsigned int category_capacity ) {

	struct MixCovars *co
		= calloc( 1, sizeof(struct MixCovars) );
	if( co ) {
		co->SAMPLE_CAPACITY   = sample_capacity;
		co->CATEGORY_CAPACITY = category_capacity;
		co->SIZEOF_BUFFERS
			= ( sample_capacity * sizeof(struct Pair) )
			+ ( category_capacity * sizeof(unsigned int) );

		// Allocate one large buffer and partition it up.
		co->samples
			= calloc( co->SIZEOF_BUFFERS, sizeof(char) );
		co->category_count
			= (unsigned int*)(co->samples
			+ sample_capacity);
		// If -anything- failed clean up any successes.
		if( NULL == co->samples ) {
			mix_destroy( co );
			return NULL;
		}
		return co;
	}
	return NULL;
}

void mix_clear( void *pv, unsigned expcat ) {

	struct MixCovars *co = (struct MixCovars *)pv;

	assert( expcat <= co->CATEGORY_CAPACITY );

	co->sample_count = 0;
	co->expected_categories = expcat;
	co->observed_categories = 0;
	co->mean_rank    = 0.0;
	co->sum_sq_dev   = 0.0;
	co->sum_dev_prod = 0.0;
	memset( co->edge, 0, sizeof(co->edge) );
	co->sum_dev_prod = 0.0;

	memset( co->samples, 0, co->SIZEOF_BUFFERS );
}

void mix_push( void *pv, float num, unsigned int cat ) {

	struct MixCovars *co = (struct MixCovars *)pv;

	assert( cat              < co->expected_categories );
	assert( co->sample_count < co->SAMPLE_CAPACITY     );

	if( 0 == co->category_count[ cat ] )
		co->observed_categories++;

	co->category_count[ cat ] += 1;

	// Not updating edges in here because the number of conditionals
	// executed for sample counts > 32 exceeds the work to find the
	// edges post-sample accumulation.
	co->samples[ co->sample_count ].cv = num;
	co->samples[ co->sample_count ].dv = cat;
	co->sample_count++;
}

size_t mix_size( void *pv ) {
	struct MixCovars *co = (struct MixCovars *)pv;
	return co->sample_count;
}

/**
  * - Establish whether or not Spearman is even sensible (according to
  *   whether or not the categorical variable is binary).
  * - Compute the mean rank used throughout subsequent calcs.
  * - Compute the mean ranks of the categorical values assuming
  *   they -are- binary.
  */
bool mix_complete( void *pv ) {

	struct MixCovars *co = (struct MixCovars *)pv;
	if( mix_degenerate( co ) ) return false;

	co->mean_rank = (1.0 + co->sample_count)/2.0;

	co->edge[1].index = _maxCat( co ); // ...this is used by kruskal_wallis, but

	assert( co->edge[1].index < co->expected_categories );

	// ...the remainder of the edge[] struct is only used if...

	if( co->observed_categories == 2 ) {

		// All the following is relevant only to Spearman rho calculation.

		co->edge[0].index = _minCat( co );
		assert( co->edge[0].index < co->expected_categories );

		co->edge[0].count = co->category_count[ co->edge[0].index ];
		co->edge[0].meanRank = MEAN_RANK_OF_TIES( 0, co->edge[0].count );

		co->edge[1].count = co->category_count[ co->edge[1].index ];

		if( co->edge[0].count < arg_min_mixb_count ||
			co->edge[1].count < arg_min_mixb_count )
			return false; // another degeneracy class

		co->edge[1].meanRank = MEAN_RANK_OF_TIES( co->edge[0].count, co->sample_count );

		assert( (co->edge[0].count + co->edge[1].count) == co->sample_count );
	}
	return true;
}

bool mix_degenerate( void *pv ) {
	struct MixCovars *co = (struct MixCovars *)pv;
	return co->observed_categories < 2 || co->sample_count < 2;
}

bool mix_categoricalIsBinary( void *pv ) {
	struct MixCovars *co = (struct MixCovars *)pv;
	return co->observed_categories == 2;
}

/**
 * Calculates the Kruskal-Wallis statistic and optionally a p-value.
 * Importantly, it does it in a one pass iteration over the data.
 */
int mix_kruskal_wallis( void *pv, struct Statistic *result ) {

	struct MixCovars *co = (struct MixCovars *)pv;
	const unsigned int N
		= co->sample_count ;
	const size_t SIZEOF_SUMS
		= co->expected_categories * sizeof(double);
	double *rank_sum
		= (double*)alloca( SIZEOF_SUMS );
	double numerator = 0.0;

	memset( rank_sum, 0, SIZEOF_SUMS   );

	result->extra_value[0] = _rank_sums( co, rank_sum );

	for(unsigned int i = 0; i <= co->edge[1].index; i++ ) {
		if( co->category_count[i] > 0 ) {
			rank_sum[i] /= co->category_count[i];
			double delta
				= rank_sum[i] - co->mean_rank;
			numerator
				+= ( co->category_count[i] * delta * delta );
		}
	}

	result->name
		= "Kruskal-Wallis_K";
	result->sample_count
		= N;
	result->value
		= (N-1) * ( numerator / co->sum_sq_dev );
	result->probability
		= gsl_cdf_chisq_Q( result->value, co->observed_categories-1 );

	return 0;
}

/**
 * P-value returned is ONE-SIDED. Double it for two-sided.
 */
#ifdef HAVE_MANN_WHITNEY
int mix_mann_whitney( void *pv, struct Statistic *result ) {

	struct MixCovars *co = (struct MixCovars *)pv;
	const size_t SIZEOF_SUMS
		= co->expected_categories * sizeof(double);
	double *rank_sum
		= (double*)alloca( SIZEOF_SUMS );
	memset( rank_sum, 0, SIZEOF_SUMS   );

	ms->ties
		= _rank_sums( co, rank_sum );

	// TODO: If too many ties, the following method of calculating U
	//       becomes invalid. What to do?

	const unsigned int SQ
		= category_count[0]*category_count[1];

	// If either group count is 0 we're done...

	if( SQ ) {

		const double U1
			= rank_sum[0] - (category_count[0]*(category_count[0]+1.0))/2.0;
		const double U2
			= rank_sum[1] - (category_count[1]*(category_count[1]+1.0))/2.0;

		assert( U1+U2 == SQ );

		ms->U = U1 < U2 ? U1 : U2;

		const double mean_U
			= SQ / 2.0;
		const double sigma_U
			= sqrt( SQ*( category_count[0] + category_count[1] + 1.0 ) / 12.0 );

		cs->P = gsl_cdf_ugaussian_P( ( ms->U - mean_U ) / sigma_U );
	}
	return 0;
}
#endif

/**
  * This calculation depends on a quantities calculated within
  * kruskal_wallis (or mann_whitney), so must be called after
  * either of them.
  * (This is also the only reason this method can be const.)
  */
double mix_spearman_rho( void *pv ) {

	struct MixCovars *co = (struct MixCovars *)pv;
	/**
	  * denominator requires sum of square deviation for
	  * the binary variables which collapses to...
	  */
	double SQD0 = (co->edge[0].meanRank - co->mean_rank ); SQD0 *= SQD0;
	double SQD1 = (co->edge[1].meanRank - co->mean_rank ); SQD1 *= SQD1;

	const double SUMSQD
		= co->edge[0].count * SQD0 + co->edge[1].count * SQD1;
	/**
	  * And the sum of squared deviation of the numeric covariate's
	  * ranks (sum_sq_dev) was (necessarily) calculated earlier!
	  */
	return co->sum_dev_prod / sqrt( co->sum_sq_dev * SUMSQD );
}

#ifdef _UNITTEST_MIX_

/**
 * This is intended to be exercised with the following R script
 * that generates a small table of grouped floating-point values
 * executes a Kruskal-Wallis test on it, and dumps the table
 * to a tab-delimited file with the test results as a comment on
 * the first line.
 *
 * x <- data.frame(
 * 		cat=as.integer(gl(3,4))-1,
 * 		num=c( runif(4)*10, runif(4)*20, 10+runif(4)*10 ) );
 * k <- with( x, kruskal.test( num, cat ) );
 * cat( sprintf( "# K=%f p-value=%f\n", k$statistic, k$p.value ), file="foo.tab" );
 * write.table( x, 'foo.tab', quote=F, sep='\t', row.names=F, col.names=F, append=TRUE );
 */
#include <err.h>

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

	if( argc >= 3 ) {

		struct Statistic result;
		const int EXPCAT
			= atoi( argv[1] );
		struct MixCovars *accum
		   = mix_create( atoi( argv[2] ), MAX_CATEGORY_COUNT );
		FILE *fp
			= argc > 3
			? fopen( argv[3], "r" )
			: stdin;
		char *line = NULL;
		size_t n = 0;

		mix_clear( accum, EXPCAT );

		while( getline( &line, &n, fp ) > 0 ) {
			unsigned int cat;
			float        num;
			if( line[0] == '#' ) {
				fputs( line, stdout );
				continue;
			}
			if( 2 == sscanf( line, "%d\t%f\n", &cat, &num ) ) {
				mix_push( accum, num, cat );
			} else {
				fprintf( stderr, "Failure parsing line: %s", line );
			}
		}
		free( line );
		fclose( fp );

		// Always calculate Kruskal-Wallis since it's valid for any groups
		// count >= 2...

		if( mix_complete( accum ) ) {
			memset( &result, 0, sizeof(struct Statistic) );
			if( mix_kruskal_wallis( accum, &result ) == 0 ) {
				printf( "K=%f, p-value=%f", result.value, result.probability );
				if( mix_categoricalIsBinary( accum ) )
					printf( ", spearman=%f\n", mix_spearman_rho(accum) );
				else
					fputc( '\n', stdout );
			} else {
				printf( "error\n" );
			}
		}

#ifdef HAVE_MANN_WHITNEY
		// Do a Mann-Whitney, too, if there are only 2 groups.
		if( mix_categoricalIsBinary(accum) == 2 ) {
		}
#endif
		mix_destroy( accum );
	} else
		err( -1, "%s <categories> <sample count> [ <input file> ]", argv[0] );
	return 0;
}
#endif

max_stars_repo_name: IlyaLab/kramtools

max_issues_repo_name: IlyaLab/kramtools

max_forks_repo_name: IlyaLab/kramtools