text:
/**
 * \file RIAAFilter.h
 * Filters for RIAA correction
 */

#ifndef ATK_EQ_RIAAFILTER_H
#define ATK_EQ_RIAAFILTER_H

#include <gsl/gsl>

#include <ATK/EQ/SecondOrderFilter.h>

namespace ATK
{
  /// RIAA coefficients used for RIAA correction
  template<typename DataType_>
  class RIAACoefficients : public SecondOrderCoreCoefficients<DataType_>
  {
  public:
    /// Simplify parent calls
    using Parent = SecondOrderCoreCoefficients<DataType_>;
    using typename Parent::AlignedScalarVector;
    using typename Parent::DataType;
    using CoeffDataType = typename TypeTraits<DataType>::Scalar;
    using Parent::input_sampling_rate;
    using Parent::output_sampling_rate;
    using Parent::setup;
  protected:
    using Parent::in_order;
    using Parent::out_order;
    using Parent::coefficients_in;
    using Parent::coefficients_out;

    void setup() override;

  public:
    /*!
     * @brief Constructor
     * @param nb_channels is the number of input and output channels
     */
    explicit RIAACoefficients(gsl::index nb_channels = 1);
  };

  /// RIAA coefficients used to create the inverse compensation (for mastering engineers)
  template<typename DataType_>
  class InverseRIAACoefficients : public SecondOrderCoreCoefficients<DataType_>
  {
  public:
    /// Simplify parent calls
    using Parent = SecondOrderCoreCoefficients<DataType_>;
    using typename Parent::AlignedScalarVector;
    using typename Parent::DataType;
    using CoeffDataType = typename TypeTraits<DataType>::Scalar;
    using Parent::input_sampling_rate;
    using Parent::output_sampling_rate;
    using Parent::setup;
  protected:
    using Parent::in_order;
    using Parent::out_order;
    using Parent::coefficients_in;
    using Parent::coefficients_out;

    void setup() override;

  public:
    /*!
     * @brief Constructor
     * @param nb_channels is the number of input and output channels
     */
    explicit InverseRIAACoefficients(gsl::index nb_channels = 1);
  };
}

#endif

ext: h

max_stars_repo_name: D-J-Roberts/AudioTK

max_issues_repo_name: D-J-Roberts/AudioTK

max_forks_repo_name: D-J-Roberts/AudioTK


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/

#ifndef __StgFEM_Discretisation_Discretisation_h__
#define __StgFEM_Discretisation_Discretisation_h__

	#include "types.h"

	#include <petsc.h>
	#include <petscvec.h>
	#include <petscmat.h>
	#include <petscksp.h>
#if( (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR==0) || (PETSC_VERSION_MAJOR<3) )
	#include <petscmg.h>
#endif
	#include <petscsnes.h>
	#include "PETScErrorChecking.h"

	#include "FeMesh_Algorithms.h"
	#include "FeMesh_ElementType.h"
	#include "ElementType.h"
	#include "ElementType_Register.h"
	#include "ConstantElementType.h"
	#include "LinearElementType.h"
	#include "BilinearElementType.h"
	#include "TrilinearElementType.h"
	#include "Biquadratic.h"
	#include "Triquadratic.h"
	#include "LinearTriangleElementType.h"
	#include "BilinearInnerElType.h"
	#include "TrilinearInnerElType.h"
    #include "dQ13DElementType.h"
    #include "dQ12DElementType.h"
	#include "dQ1Generator.h"

    #include "Element.h"
	#include "FeMesh.h"
	#include "C0Generator.h"
	#include "C2Generator.h"
	#include "Inner2DGenerator.h"
	#include "LinkedDofInfo.h"
	#include "FeEquationNumber.h"
	#include "FeVariable.h"

    #include "IrregularMeshGaussLayout.h"

	#include "Init.h"
	#include "Finalise.h"

#endif /* __StgFEM_Discretisation_Discretisation_h__ */


text:
#include <gsl/gsl_vector.h>

#include "icp.h"
#include "../csm_all.h"

int compatible(struct sm_params*params, int i, int j);

int compatible(struct sm_params*params, int i, int j) {
	if(!params->do_alpha_test) return 1;

	double theta0 = 0; /* FIXME */
	if((params->laser_sens->alpha_valid[i]==0) ||
	 (params->laser_ref->alpha_valid[j]==0))
	 return 1;

	double alpha_i = params->laser_sens->alpha[i];
	double alpha_j = params->laser_ref->alpha[j];
	double tolerance = deg2rad(params->do_alpha_test_thresholdDeg);

	/** FIXME remove alpha test */
	double theta = angleDiff(alpha_j, alpha_i);
	if(fabs(angleDiff(theta,theta0))>
		tolerance+deg2rad(params->max_angular_correction_deg)) {
	 return 0;
	} else {
	 return 1;
	}
}

void find_correspondences(struct sm_params*params) {
	const LDP laser_ref  = params->laser_ref;
	const LDP laser_sens = params->laser_sens;

	int i;
	for(i=0;i<laser_sens->nrays;i++) {
		if(!ld_valid_ray(laser_sens,i)) {
/*			sm_debug("dumb: i %d is invalid \n", i);*/
			ld_set_null_correspondence(laser_sens, i);
			continue;
		}

		double *p_i_w = laser_sens->points_w[i].p;

		int j1 = -1;
		double best_dist = 10000;

		int from; int to; int start_cell;
		possible_interval(p_i_w, laser_ref, params->max_angular_correction_deg,
			params->max_linear_correction, &from, &to, &start_cell);

/*		sm_debug("dumb: i %d from  %d to %d \n", i, from, to); */
		int j;
		for(j=from;j<=to;j++) {
			if(!ld_valid_ray(laser_ref,j)) {
/*				sm_debug("dumb: i %d      j %d invalid\n", i, j);*/
				continue;
			}
			double dist = distance_squared_d(p_i_w, laser_ref->points[j].p);
/*			sm_debug("dumb: i %d j1 %d j %d d %f\n", i,j1,j,dist);*/
			if(dist>square(params->max_correspondence_dist)) continue;

			if( (-1 == j1) || (dist < best_dist) ) {
				if(compatible(params, i, j)) {
					j1 = j;
					best_dist = dist;
				}
			}
		}

		if(j1==-1) {/* no match */
			ld_set_null_correspondence(laser_sens, i);
			continue;
		}
		/* Do not match with extrema*/
		if(j1==0 || (j1 == (laser_ref->nrays-1))) {/* no match */
			ld_set_null_correspondence(laser_sens, i);
			continue;
		}

		int j2;
		int j2up   = ld_next_valid_up   (laser_ref, j1);
		int j2down = ld_next_valid_down (laser_ref, j1);
		if((j2up==-1)&&(j2down==-1)) {
			ld_set_null_correspondence(laser_sens, i);
			continue;
		}
		if(j2up  ==-1) { j2 = j2down; } else
		if(j2down==-1) { j2 = j2up; } else {
			double dist_up   = distance_squared_d(p_i_w, laser_ref->points[j2up  ].p);
			double dist_down = distance_squared_d(p_i_w, laser_ref->points[j2down].p);
			j2 = dist_up < dist_down ? j2up : j2down;
		}

		ld_set_correspondence(laser_sens, i, j1, j2);
		laser_sens->corr[i].dist2_j1 = best_dist;
		laser_sens->corr[i].type =
			params->use_point_to_line_distance ? corr_pl : corr_pp;

	}

}

ext: c


text:
/*
#
# * The source code in this file is developed independently by NEC Corporation.
#
# # NLCPy License #
#
#     Copyright (c) 2020-2021 NEC Corporation
#     All rights reserved.
#
#     Redistribution and use in source and binary forms, with or without
#     modification, are permitted provided that the following conditions are met:
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright notice,
#       this list of conditions and the following disclaimer in the documentation
#       and/or other materials provided with the distribution.
#     * Neither NEC Corporation 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 <cblas.h>
#include "nlcpy.h"

uint64_t wrapper_cblas_sdot(ve_arguments *args, int32_t *psw)
{
#ifdef _OPENMP
#pragma omp single
#endif /* _OPENMP */
{
    ve_array *x = &(args->binary.x);
    ve_array *y = &(args->binary.y);
    ve_array *z = &(args->binary.z);

    float *px = (float *)x->ve_adr;
    if (px == NULL) {
        px = (float *)nlcpy__get_scalar(x);
        if (px == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    float *py = (float *)y->ve_adr;
    if (py == NULL) {
        py = (float *)nlcpy__get_scalar(y);
        if (py == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    float *pz = (float *)z->ve_adr;
    if (pz == NULL) {
       return (uint64_t)NLCPY_ERROR_MEMORY;
    }
    assert(x->ndim <= 1);
    assert(y->ndim <= 1);
    assert(z->ndim <= 1);
    assert(x->size == y->size);

    *pz = cblas_sdot(x->size, px, x->strides[0] / x->itemsize,
                            py, y->strides[0] / y->itemsize);
} /* omp single */
    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_ddot(ve_arguments *args, int32_t *psw)
{
#ifdef _OPENMP
#pragma omp single
#endif /* _OPENMP */
{
    ve_array *x = &(args->binary.x);
    ve_array *y = &(args->binary.y);
    ve_array *z = &(args->binary.z);

    double *px = (double *)x->ve_adr;
    if (px == NULL) {
        px = (double *)nlcpy__get_scalar(x);
        if (px == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    double *py = (double *)y->ve_adr;
    if (py == NULL) {
        py = (double *)nlcpy__get_scalar(y);
        if (py == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    double *pz = (double *)z->ve_adr;
    if (pz == NULL) {
       return (uint64_t)NLCPY_ERROR_MEMORY;
    }
    assert(x->ndim <= 1);
    assert(y->ndim <= 1);
    assert(z->ndim <= 1);
    assert(x->size == y->size);

    *pz = cblas_ddot(x->size, px, x->strides[0] / x->itemsize,
                            py, y->strides[0] / y->itemsize);
} /* omp single */
    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_cdotu_sub(ve_arguments *args, int32_t *psw)
{
#ifdef _OPENMP
#pragma omp single
#endif /* _OPENMP */
{
    ve_array *x = &(args->binary.x);
    ve_array *y = &(args->binary.y);
    ve_array *z = &(args->binary.z);

    float _Complex *px = (float _Complex *)x->ve_adr;
    if (px == NULL) {
        px = (float _Complex *)nlcpy__get_scalar(x);
        if (px == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    float _Complex *py = (float _Complex *)y->ve_adr;
    if (py == NULL) {
        py = (float _Complex *)nlcpy__get_scalar(y);
        if (py == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    float _Complex *pz = (float _Complex *)z->ve_adr;
    if (pz == NULL) {
       return (uint64_t)NLCPY_ERROR_MEMORY;
    }
    assert(x->ndim <= 1);
    assert(y->ndim <= 1);
    assert(z->ndim <= 1);
    assert(x->size == y->size);

    cblas_cdotu_sub(x->size, px, x->strides[0] / x->itemsize,
                        py, y->strides[0] / y->itemsize, pz);
} /* omp single */
    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_zdotu_sub(ve_arguments *args, int32_t *psw)
{
#ifdef _OPENMP
#pragma omp single
#endif /* _OPENMP */
{
    ve_array *x = &(args->binary.x);
    ve_array *y = &(args->binary.y);
    ve_array *z = &(args->binary.z);

    double _Complex *px = (double _Complex *)x->ve_adr;
    if (px == NULL) {
        px = (double _Complex *)nlcpy__get_scalar(x);
        if (px == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    double _Complex *py = (double _Complex *)y->ve_adr;
    if (py == NULL) {
        py = (double _Complex *)nlcpy__get_scalar(y);
        if (py == NULL) {
            return NLCPY_ERROR_MEMORY;
        }
    }
    double _Complex *pz = (double _Complex *)z->ve_adr;
    if (pz == NULL) {
       return (uint64_t)NLCPY_ERROR_MEMORY;
    }
    assert(x->ndim <= 1);
    assert(y->ndim <= 1);
    assert(z->ndim <= 1);
    assert(x->size == y->size);

    cblas_zdotu_sub(x->size, px, x->strides[0] / x->itemsize,
                        py, y->strides[0] / y->itemsize, pz);
} /* omp single */
    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_sgemm(ve_arguments *args, int32_t *psw)
{
    const int32_t order = args->gemm.order;
    const int32_t transA = args->gemm.transA;
    const int32_t transB = args->gemm.transB;
    const int32_t m = args->gemm.m;
    const int32_t n = args->gemm.n;
    const int32_t k = args->gemm.k;
    const float alpha = *((float *)nlcpy__get_scalar(&(args->gemm.alpha)));
    float* const a = (float *)args->gemm.a.ve_adr;
    const int32_t lda = args->gemm.lda;
    float* const b = (float *)args->gemm.b.ve_adr;
    const int32_t ldb = args->gemm.ldb;
    const float beta = *((float *)nlcpy__get_scalar(&(args->gemm.beta)));
    float* const c = (float *)args->gemm.c.ve_adr;
    const int32_t ldc = args->gemm.ldc;

    if (a == NULL || b == NULL || c == NULL) {
        return NLCPY_ERROR_MEMORY;
    }

#ifdef _OPENMP
    const int32_t nt = omp_get_num_threads();
    const int32_t it = omp_get_thread_num();
#else
    const int32_t nt = 1;
    const int32_t it = 0;
#endif /* _OPENMP */

    const int32_t m_s = m * it / nt;
    const int32_t m_e = m * (it + 1) / nt;
    const int32_t m_d = m_e - m_s;
    const int32_t n_s = n * it / nt;
    const int32_t n_e = n * (it + 1) / nt;
    const int32_t n_d = n_e - n_s;

    int32_t mode = 1;
    if ( n > nt ) {
        mode = 2;
    }
    int32_t iar, iac, ibr, ibc, icr, icc;
    if (transA == CblasNoTrans ) {
        iar = 1;
        iac = lda;
    } else {
        iar = lda;
        iac = 1;
    }
    if (transB == CblasNoTrans ) {
        ibr = 1;
        ibc = ldb;
    } else {
        ibr = ldb;
        ibc = 1;
    }
    if (order == CblasColMajor ) {
        icr = 1;
        icc = ldc;
    } else {
        icr = ldc;
        icc = 1;
    }

    if (order == CblasColMajor) {
        if ( mode == 1 ) {
            // split 'm'
            cblas_sgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iar, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_sgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibc, ldb, beta, c + n_s * icc, ldc);
        }
    } else {
        if ( mode == 1 ) {
            // split 'm'
            cblas_sgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iac, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_sgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibr, ldb, beta, c + n_s * icc, ldc);
        }
    }

    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_dgemm(ve_arguments *args, int32_t *psw)
{
    const int32_t order = args->gemm.order;
    const int32_t transA = args->gemm.transA;
    const int32_t transB = args->gemm.transB;
    const int32_t m = args->gemm.m;
    const int32_t n = args->gemm.n;
    const int32_t k = args->gemm.k;
    const double alpha = *((double *)nlcpy__get_scalar(&(args->gemm.alpha)));
    double* const a = (double *)args->gemm.a.ve_adr;
    const int32_t lda = args->gemm.lda;
    double* const b = (double *)args->gemm.b.ve_adr;
    const int32_t ldb = args->gemm.ldb;
    const double beta = *((double *)nlcpy__get_scalar(&(args->gemm.beta)));
    double* const c = (double *)args->gemm.c.ve_adr;
    const int32_t ldc = args->gemm.ldc;

    if (a == NULL || b == NULL || c == NULL) {
        return NLCPY_ERROR_MEMORY;
    }

#ifdef _OPENMP
    const int32_t nt = omp_get_num_threads();
    const int32_t it = omp_get_thread_num();
#else
    const int32_t nt = 1;
    const int32_t it = 0;
#endif /* _OPENMP */

    const int32_t m_s = m * it / nt;
    const int32_t m_e = m * (it + 1) / nt;
    const int32_t m_d = m_e - m_s;
    const int32_t n_s = n * it / nt;
    const int32_t n_e = n * (it + 1) / nt;
    const int32_t n_d = n_e - n_s;

    int32_t mode = 1;
    if ( n > nt ) {
        mode = 2;
    }
    int32_t iar, iac, ibr, ibc, icr, icc;
    if (transA == CblasNoTrans ) {
        iar = 1;
        iac = lda;
    } else {
        iar = lda;
        iac = 1;
    }
    if (transB == CblasNoTrans ) {
        ibr = 1;
        ibc = ldb;
    } else {
        ibr = ldb;
        ibc = 1;
    }
    if (order == CblasColMajor ) {
        icr = 1;
        icc = ldc;
    } else {
        icr = ldc;
        icc = 1;
    }

    if (order == CblasColMajor) {
        if ( mode == 1 ) {
            // split 'm'
            cblas_dgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iar, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_dgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibc, ldb, beta, c + n_s * icc, ldc);
        }
    } else {
        if ( mode == 1 ) {
            // split 'm'
            cblas_dgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iac, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_dgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibr, ldb, beta, c + n_s * icc, ldc);
        }
    }

    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_cgemm(ve_arguments *args, int32_t *psw)
{
    const int32_t order = args->gemm.order;
    const int32_t transA = args->gemm.transA;
    const int32_t transB = args->gemm.transB;
    const int32_t m = args->gemm.m;
    const int32_t n = args->gemm.n;
    const int32_t k = args->gemm.k;
    const void *alpha = (void *)nlcpy__get_scalar(&(args->gemm.alpha));
    if (alpha == NULL) return (uint64_t)NLCPY_ERROR_MEMORY;
    float  _Complex* const a = (float  _Complex *)args->gemm.a.ve_adr;
    const int32_t lda = args->gemm.lda;
    float  _Complex* const b = (float  _Complex *)args->gemm.b.ve_adr;
    const int32_t ldb = args->gemm.ldb;
    const void *beta = (void *)nlcpy__get_scalar(&(args->gemm.beta));
    if (beta == NULL) return (uint64_t)NLCPY_ERROR_MEMORY;
    float  _Complex* const c = (float  _Complex *)args->gemm.c.ve_adr;
    const int32_t ldc = args->gemm.ldc;

    if (a == NULL || b == NULL || c == NULL) {
        return NLCPY_ERROR_MEMORY;
    }

#ifdef _OPENMP
    const int32_t nt = omp_get_num_threads();
    const int32_t it = omp_get_thread_num();
#else
    const int32_t nt = 1;
    const int32_t it = 0;
#endif /* _OPENMP */

    const int32_t m_s = m * it / nt;
    const int32_t m_e = m * (it + 1) / nt;
    const int32_t m_d = m_e - m_s;
    const int32_t n_s = n * it / nt;
    const int32_t n_e = n * (it + 1) / nt;
    const int32_t n_d = n_e - n_s;

    int32_t mode = 1;
    if ( n > nt ) {
        mode = 2;
    }
    int32_t iar, iac, ibr, ibc, icr, icc;
    if (transA == CblasNoTrans ) {
        iar = 1;
        iac = lda;
    } else {
        iar = lda;
        iac = 1;
    }
    if (transB == CblasNoTrans ) {
        ibr = 1;
        ibc = ldb;
    } else {
        ibr = ldb;
        ibc = 1;
    }
    if (order == CblasColMajor ) {
        icr = 1;
        icc = ldc;
    } else {
        icr = ldc;
        icc = 1;
    }

    if (order == CblasColMajor) {
        if ( mode == 1 ) {
            // split 'm'
            cblas_cgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iar, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_cgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibc, ldb, beta, c + n_s * icc, ldc);
        }
    } else {
        if ( mode == 1 ) {
            // split 'm'
            cblas_cgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iac, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_cgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibr, ldb, beta, c + n_s * icc, ldc);
        }
    }

    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

uint64_t wrapper_cblas_zgemm(ve_arguments *args, int32_t *psw)
{
    const int32_t order = args->gemm.order;
    const int32_t transA = args->gemm.transA;
    const int32_t transB = args->gemm.transB;
    const int32_t m = args->gemm.m;
    const int32_t n = args->gemm.n;
    const int32_t k = args->gemm.k;
    const void *alpha = (void *)nlcpy__get_scalar(&(args->gemm.alpha));
    if (alpha == NULL) return (uint64_t)NLCPY_ERROR_MEMORY;
    double _Complex* const a = (double _Complex *)args->gemm.a.ve_adr;
    const int32_t lda = args->gemm.lda;
    double _Complex* const b = (double _Complex *)args->gemm.b.ve_adr;
    const int32_t ldb = args->gemm.ldb;
    const void *beta = (void *)nlcpy__get_scalar(&(args->gemm.beta));
    if (beta == NULL) return (uint64_t)NLCPY_ERROR_MEMORY;
    double _Complex* const c = (double _Complex *)args->gemm.c.ve_adr;
    const int32_t ldc = args->gemm.ldc;

    if (a == NULL || b == NULL || c == NULL) {
        return NLCPY_ERROR_MEMORY;
    }

#ifdef _OPENMP
    const int32_t nt = omp_get_num_threads();
    const int32_t it = omp_get_thread_num();
#else
    const int32_t nt = 1;
    const int32_t it = 0;
#endif /* _OPENMP */

    const int32_t m_s = m * it / nt;
    const int32_t m_e = m * (it + 1) / nt;
    const int32_t m_d = m_e - m_s;
    const int32_t n_s = n * it / nt;
    const int32_t n_e = n * (it + 1) / nt;
    const int32_t n_d = n_e - n_s;

    int32_t mode = 1;
    if ( n > nt ) {
        mode = 2;
    }
    int32_t iar, iac, ibr, ibc, icr, icc;
    if (transA == CblasNoTrans ) {
        iar = 1;
        iac = lda;
    } else {
        iar = lda;
        iac = 1;
    }
    if (transB == CblasNoTrans ) {
        ibr = 1;
        ibc = ldb;
    } else {
        ibr = ldb;
        ibc = 1;
    }
    if (order == CblasColMajor ) {
        icr = 1;
        icc = ldc;
    } else {
        icr = ldc;
        icc = 1;
    }

    if (order == CblasColMajor) {
        if ( mode == 1 ) {
            // split 'm'
            cblas_zgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iar, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_zgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibc, ldb, beta, c + n_s * icc, ldc);
        }
    } else {
        if ( mode == 1 ) {
            // split 'm'
            cblas_zgemm(order, transA, transB, m_d, n, k, alpha, a + m_s * iac, lda, b, ldb, beta, c + m_s * icr, ldc);
        } else {
            // split 'n'
            cblas_zgemm(order, transA, transB, m, n_d, k, alpha, a, lda, b + n_s * ibr, ldb, beta, c + n_s * icc, ldc);
        }
    }

    retrieve_fpe_flags(psw);
    return (uint64_t)NLCPY_ERROR_OK;
}

max_stars_repo_name: SX-Aurora/nlcpy

max_issues_repo_name: SX-Aurora/nlcpy

max_forks_repo_name: SX-Aurora/nlcpy


text:
/* Modified Cholesky Decomposition
 *
 * Copyright (C) 2016 Patrick Alken
 *
 * This 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, or (at your option) any
 * later version.
 *
 * This source 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.
 */

#include <config.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_permute_vector.h>

#include "cholesky_common.c"

/*
 * This module contains routines related to the Modified Cholesky
 * Decomposition, which factors a symmetric indefinite matrix A as
 *
 * P (A + E) P^T = L D L^T
 *
 * where:
 * P: permutation matrix
 * E: small, non-negative diagonal matrix
 * L: unit lower triangular matrix
 * D: strictly positive diagonal matrix
 *
 * These routines follow these works closely:
 *
 * [1] P. E. Gill, W. Murray, M. H. Wright, Practical Optimization,
 *     Academic Press, 1981.
 *
 * [2] Dennis and Schnabel, Numerical Methods for Unconstrained Optimization
 *     and Nonlinear Equations, SIAM, 1996
 */

static size_t mcholesky_maxabs(const gsl_vector * v, double *maxabs);

/*
gsl_linalg_mcholesky_decomp()
  Perform Pivoted Modified Cholesky LDLT decomposition of a symmetric positive
indefinite matrix:

P (A + E) P^T = L D L^T

Inputs: A - (input) symmetric, positive indefinite matrix,
                    stored in lower triangle
            (output) lower triangle contains L; diagonal contains D
        p - (output) permutation matrix P
        E - (output) perturbation matrix E

Return: success/error

Notes:
1) Based on algorithm 4.2.2 (Outer Product LDLT with Pivoting) of
Golub and Van Loan, Matrix Computations (4th ed), with modifications
described in [1] and [2]

2) E can be set to NULL if not required
*/

int
gsl_linalg_mcholesky_decomp (gsl_matrix * A, gsl_permutation * p,
                             gsl_vector * E)
{
  const size_t N = A->size1;

  if (N != A->size2)
    {
      GSL_ERROR("LDLT decomposition requires square matrix", GSL_ENOTSQR);
    }
  else if (p->size != N)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else
    {
      const double delta = GSL_DBL_EPSILON;
      double beta;
      double gamma = 0.0;
      double xi = 0.0;
      gsl_vector_view diag = gsl_matrix_diagonal(A);
      size_t i, j;

      /* save a copy of A in upper triangle (for later rcond calculation) */
      gsl_matrix_transpose_tricpy('L', 0, A, A);

      gsl_permutation_init(p);

      /* compute:
       * gamma = max | A_{ii} |
       * xi = max_{i \ne j} | A_{ij} |
       */
      for (i = 0; i < N; ++i)
        {
          double aii = gsl_matrix_get(A, i, i);

          gamma = GSL_MAX(gamma, fabs(aii));

          for (j = 0; j < i; ++j)
            {
              double aij = gsl_matrix_get(A, i, j);
              xi = GSL_MAX(xi, fabs(aij));
            }
        }

      /* compute:
       * beta = sqrt[ max { gamma, xi/nu, eps } ]
       * with: nu = max{ sqrt(N^2 - 1), 1 }
       */
      if (N == 1)
        {
          beta = GSL_MAX(GSL_MAX(gamma, xi), GSL_DBL_EPSILON);
        }
      else
        {
          double nu = sqrt(N*N - 1.0);
          beta = GSL_MAX(GSL_MAX(gamma, xi / nu), GSL_DBL_EPSILON);
        }

      beta = sqrt(beta);

      for (j = 0; j < N; ++j)
        {
          double ajj, thetaj, u, alpha, alphainv;
          gsl_vector_view w;
          size_t q;

          /* compute q = max_idx { A_jj, ..., A_nn } */
          w = gsl_vector_subvector(&diag.vector, j, N - j);
          q = mcholesky_maxabs(&w.vector, NULL) + j;

          gsl_permutation_swap(p, q, j);
          cholesky_swap_rowcol(A, q, j);

          /* theta_j = max_{j+1 <= i <= n} |A_{ij}| */
          if (j < N - 1)
            {
              w = gsl_matrix_subcolumn(A, j, j + 1, N - j - 1);
              mcholesky_maxabs(&w.vector, &thetaj);
            }
          else
            {
              thetaj = 0.0;
            }

          u = thetaj / beta;

          /* compute alpha = d_j */
          ajj = gsl_matrix_get(A, j, j);
          alpha = GSL_MAX(GSL_MAX(delta, fabs(ajj)), u * u);
          alphainv = 1.0 / alpha;

          if (j < N - 1)
            {
              /* v = A(j+1:n, j) */
              gsl_vector_view v = gsl_matrix_subcolumn(A, j, j + 1, N - j - 1);

              /* m = A(j+1:n, j+1:n) */
              gsl_matrix_view m = gsl_matrix_submatrix(A, j + 1, j + 1, N - j - 1, N - j - 1);

              /* m = m - v v^T / alpha */
              gsl_blas_dsyr(CblasLower, -alphainv, &v.vector, &m.matrix);

              /* v = v / alpha */
              gsl_vector_scale(&v.vector, alphainv);

            }

          if (E)
            gsl_vector_set(E, j, alpha - ajj);

          gsl_matrix_set(A, j, j, alpha);
        }

      if (E)
        {
          /* we currently have: P A P^T + E = L D L^T, permute E
           * so that we have: P (A + E) P^T = L D L^T */
          gsl_permute_vector_inverse(p, E);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_mcholesky_solve(const gsl_matrix * LDLT,
                           const gsl_permutation * p,
                           const gsl_vector * b,
                           gsl_vector * x)
{
  int status = gsl_linalg_pcholesky_solve(LDLT, p, b, x);
  return status;
}

int
gsl_linalg_mcholesky_svx(const gsl_matrix * LDLT,
                         const gsl_permutation * p,
                         gsl_vector * x)
{
  int status = gsl_linalg_pcholesky_svx(LDLT, p, x);
  return status;
}

int
gsl_linalg_mcholesky_rcond (const gsl_matrix * LDLT, const gsl_permutation * p,
                            double * rcond, gsl_vector * work)
{
  int status = gsl_linalg_pcholesky_rcond(LDLT, p, rcond, work);
  return status;
}

int
gsl_linalg_mcholesky_invert(const gsl_matrix * LDLT, const gsl_permutation * p,
                            gsl_matrix * Ainv)
{
  int status = gsl_linalg_pcholesky_invert(LDLT, p, Ainv);
  return status;
}

/*
mcholesky_maxabs()
  Compute:

val = max_i |v_i|

Inputs: v      - vector
        maxabs - (output) max abs value

Return: index corresponding to max_i |v_i|
*/

static size_t
mcholesky_maxabs(const gsl_vector * v, double *maxabs)
{
  const size_t n = v->size;
  size_t i;
  size_t idx = 0;
  double max = gsl_vector_get(v, idx);

  for (i = 1; i < n; ++i)
    {
      double vi = gsl_vector_get(v, i);
      double absvi = fabs(vi);

      if (absvi > max)
        {
          max = absvi;
          idx = i;
        }
    }

  if (maxabs)
    *maxabs = max;

  return idx;
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#ifndef __GSL_VECTOR_DOUBLE_H__
#define __GSL_VECTOR_DOUBLE_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/block/gsl_check_range.h>
#include <gsl/block/gsl_block_double.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;		//the number of elements in this vector
  size_t stride;	//in the block pointer, this indicates the spacing between elements
  double *data;		//a pointer to the data in this vector
  gsl_block *block;	//a pointer to the block that the data belongs to
  int owner;		//if this is 1 it means that this vector is the owner of this block, otherwise it means it inherited it from somewhere else and it shouldn't be deallocated with the vector
}
gsl_vector;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_view;

typedef _gsl_vector_view gsl_vector_view;

typedef struct
{
  gsl_vector vector;
} _gsl_vector_const_view;

typedef const _gsl_vector_const_view gsl_vector_const_view;

/* Allocation */

GSL_FUNC gsl_vector *gsl_vector_alloc (const size_t n);
GSL_FUNC gsl_vector *gsl_vector_calloc (const size_t n);

GSL_FUNC gsl_vector *gsl_vector_alloc_from_block (gsl_block * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUNC gsl_vector *gsl_vector_alloc_from_vector (gsl_vector * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUNC void gsl_vector_free (gsl_vector * v);

/* Views */

GSL_FUNC _gsl_vector_view
gsl_vector_view_array (double *v, size_t n);

GSL_FUNC _gsl_vector_view
gsl_vector_view_array_with_stride (double *base,
                                         size_t stride,
                                         size_t n);

GSL_FUNC _gsl_vector_const_view
gsl_vector_const_view_array (const double *v, size_t n);

GSL_FUNC _gsl_vector_const_view
gsl_vector_const_view_array_with_stride (const double *base,
                                               size_t stride,
                                               size_t n);

GSL_FUNC _gsl_vector_view
gsl_vector_subvector (gsl_vector *v,
                            size_t i,
                            size_t n);

GSL_FUNC _gsl_vector_view
gsl_vector_subvector_with_stride (gsl_vector *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUNC _gsl_vector_const_view
gsl_vector_const_subvector (const gsl_vector *v,
                                  size_t i,
                                  size_t n);

GSL_FUNC _gsl_vector_const_view
gsl_vector_const_subvector_with_stride (const gsl_vector *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUNC double gsl_vector_get (const gsl_vector * v, const size_t i);
GSL_FUNC void gsl_vector_set (gsl_vector * v, const size_t i, double x);

GSL_FUNC double *gsl_vector_ptr (gsl_vector * v, const size_t i);
GSL_FUNC const double *gsl_vector_const_ptr (const gsl_vector * v, const size_t i);

GSL_FUNC void gsl_vector_set_zero (gsl_vector * v);
GSL_FUNC void gsl_vector_set_all (gsl_vector * v, double x);
GSL_FUNC int gsl_vector_set_basis (gsl_vector * v, size_t i);

GSL_FUNC int gsl_vector_fread (FILE * stream, gsl_vector * v);
GSL_FUNC int gsl_vector_fwrite (FILE * stream, const gsl_vector * v);
GSL_FUNC int gsl_vector_fscanf (FILE * stream, gsl_vector * v);
GSL_FUNC int gsl_vector_fprintf (FILE * stream, const gsl_vector * v,
                              const char *format);

GSL_FUNC int gsl_vector_memcpy (gsl_vector * dest, const gsl_vector * src);

GSL_FUNC int gsl_vector_reverse (gsl_vector * v);

GSL_FUNC int gsl_vector_swap (gsl_vector * v, gsl_vector * w);
GSL_FUNC int gsl_vector_swap_elements (gsl_vector * v, const size_t i, const size_t j);

GSL_FUNC double gsl_vector_max (const gsl_vector * v);
GSL_FUNC double gsl_vector_min (const gsl_vector * v);
GSL_FUNC void gsl_vector_minmax (const gsl_vector * v, double * min_out, double * max_out);

GSL_FUNC size_t gsl_vector_max_index (const gsl_vector * v);
GSL_FUNC size_t gsl_vector_min_index (const gsl_vector * v);
GSL_FUNC void gsl_vector_minmax_index (const gsl_vector * v, size_t * imin, size_t * imax);

GSL_FUNC int gsl_vector_add (gsl_vector * a, const gsl_vector * b);
GSL_FUNC int gsl_vector_sub (gsl_vector * a, const gsl_vector * b);
GSL_FUNC int gsl_vector_mul (gsl_vector * a, const gsl_vector * b);
GSL_FUNC int gsl_vector_div (gsl_vector * a, const gsl_vector * b);
GSL_FUNC int gsl_vector_scale (gsl_vector * a, const double x);
GSL_FUNC int gsl_vector_add_constant (gsl_vector * a, const double x);

GSL_FUNC int gsl_vector_isnull (const gsl_vector * v);
GSL_FUNC int gsl_vector_ispos (const gsl_vector * v);
GSL_FUNC int gsl_vector_isneg (const gsl_vector * v);

#ifdef HAVE_INLINE

extern inline
double
gsl_vector_get (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_set (gsl_vector * v, const size_t i, double x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
double *
gsl_vector_ptr (gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (double *) (v->data + i * v->stride);
}

extern inline
const double *
gsl_vector_const_ptr (const gsl_vector * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const double *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_DOUBLE_H__ */

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#define _GNU_SOURCE
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>

#include "../core_allvars.h"

#define LUV_LOOKUPTABLE_MASS 1e6

int32_t init_UVlookup(void)
{

  // For calcUVmag == 1 we must explicitly track the stellar ages of a galaxy.
  // Then we determine the UV magnitude using the age of the stellar population.

  // Using STARBURST99 it was determined that the 1600A Luminosty depends on the mass
  // of stars formed and the time since the starburst. Furthermore, the UV luminosity linearly with the mass of stars formed.
  // That is, a starburst that forms 6.0e11Msun worth of stars will exhibit 10 times the UV luminosity as a starburst that forms 6.0e10Msun worth of stars.

  // So if we read in a table that contains the number of ionizing photons emitted from a starburst for a 6.0e10Msun episode, then we can scale our values to this lookup table
  // using log10 LUV(Msun, t) = (log10 M* - 6.0) + log10 LUV(6.0, t).

#define MAXBINS 10000

  char fname[MAX_STRING_LEN];
  FILE *LUVtable;
  int32_t num_lines = 0;
  float t, lambda, LUV, norm_spec;

  snprintf(fname, MAX_STRING_LEN - 1, ROOT_DIR "/extra/LUV_table.txt");
  LUVtable = fopen(fname, "r");
  if (LUVtable == NULL)
  {
    fprintf(stderr, "Could not open file %s\n", fname);
    return EXIT_FAILURE;
  }

  stars_LUV = calloc(MAXBINS, sizeof(*(stars_LUV)));
  if (stars_LUV == NULL)
  {
    fprintf(stderr, "Could not allocate memory for the UV luminosity bins for the tracking of stellar populations.\n");
    return EXIT_FAILURE;
  }

  while (fscanf(LUVtable, "%f %f %f %f", &t, &lambda, &LUV, &norm_spec) == 4)
  {

    stars_LUV[num_lines] = LUV;

    ++num_lines;
    if (num_lines == MAXBINS - 1)
    {
      fprintf(stderr, "Exceeding the maximum bins for the tracking of stellar populations.\n");
      return EXIT_FAILURE;
    }
  }
  fclose(LUVtable);

  // Check that the Nion lookup table had enough datapoints to cover the time we're tracking the ages for.
  if (t / 1.0e6 < STELLAR_TRACKING_TIME)
  {
    fprintf(stderr, "The final time specified in the UV luminosity lookup table is %.4f Myr. However we specified to track stellar ages over %d Myr.\n", t / 1.0e6, STELLAR_TRACKING_TIME);
    fprintf(stderr, "Either update the UV luminosity lookup table or reduce the value of STELLAR_TRACKING_TIME in `core_allvars.h`.\n");
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;

#undef MAXBINS
}

/*
Calculates the UV Luminosity (i.e., at 1600A) for a given galaxy at a specified snapshot.

**Important** See Units.

Parameters
----------

g: struct GALAXY pointer.  See `core_allvars.h` for the struct architecture.
  Pointer to the galaxy that we are calculating UV luminosity.

*LUV: Float Pointer.
  Pointer that will store the UV luminosity.

Returns
----------

EXIT_SUCCESS or EXIT_FAILURE.
  If the UV luminosity is negative, EXIT_FAILURE is returned.  Otherwise EXIT_SUCCESS is returned.

Pointer Updates
----------

*LUV.

Units
----------

The UV luminosty is returned in units of 1.0e50 erg s^-1 A^-1
*/

int32_t calc_LUV(struct GALAXY *g, float *LUV)
{

  double t;
  int32_t i, lookup_idx;

  // We ran STARBURST99 for a single mass and scale our results.
  // First convert the mass to internal code units.
  const double code_LUV_lookuptable_mass = LUV_LOOKUPTABLE_MASS*1.0e-10*Hubble_h;

  *LUV = 0.0;

  // Then go through the previous 100Myr worth of SF and sum up the UV luminosity.
  for (i = 0; i < StellarTracking_Len - 1; ++i)
  {
    if (g->Stellar_Stars[i] < 1e-10)
      continue;

    t = (i + 1) * TimeResolutionStellar; // (i + 1) because 0th entry will be at TimeResolutionSN.
    lookup_idx = t; // Find the index in the lookup table.

    *LUV += exp10(log10(g->Stellar_Stars[i] / code_LUV_lookuptable_mass) + stars_LUV[lookup_idx] - 50.0);
    //printf("t %.4f\tlookup_idx %d\tg->Stellar_Stars[i] %.4e\tstars_LUV[lookup_idx] %.4e\tRunningTotal %.4e\n", t, lookup_idx, g->Stellar_Stars[i], stars_LUV[lookup_idx], *LUV);
  }

  // The units of LUV are 1.0e50 erg s^-1 A^-1.  Hence a negative value is not allowed.
  if (*LUV < 0.0)
  {
    fprintf(stderr, "Got an LUV value of %.4e.  This MUST be a positive value.\nPrinting out information for every element used to calculate LUV.\n", *LUV);

    // Print out information for every element of the array so we can try identify the problem.
    for (i = 0; i < StellarTracking_Len; ++i)
    {
      t = (i + 1) * TimeResolutionStellar; // (i + 1) because 0th entry will be at TimeResolutionSN.
      lookup_idx = (t / 0.1); // Find the index in the lookup table.

      double total = 0.0;
      total += exp10(log10(g->Stellar_Stars[i] / code_LUV_lookuptable_mass) + stars_LUV[lookup_idx] - 50.0);
      printf("t %.4f\tlookup_idx %d\tg->Stellar_Stars[i] %.4e\tstars_LUV[lookup_idx] %.4e\tRunningTotal %.4e\n", t, lookup_idx, g->Stellar_Stars[i], stars_LUV[lookup_idx], total);
    }
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

max_stars_repo_name: jacobseiler/rsage

max_issues_repo_name: jacobseiler/rsage

max_forks_repo_name: jacobseiler/rsage


text:
/*
 *  Copyright (C) 2021 FISCO BCOS.
 *  SPDX-License-Identifier: Apache-2.0
 *  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.
 *
 * @brief interface of Table
 * @file Table.h
 * @author: xingqiangbai
 * @date: 2021-04-07
 */
#pragma once

#include "Common.h"
#include "StorageInterface.h"
#include <future>
#include <gsl/span>

namespace bcos
{
namespace storage
{
class Table
{
public:
    Table(StorageInterface* _db, TableInfo::ConstPtr _tableInfo)
      : m_storage(_db), m_tableInfo(std::move(_tableInfo))
    {}

    Table(const Table&) = default;
    Table(Table&&) = default;
    Table& operator=(const Table&) = default;
    Table& operator=(Table&&) = default;
    ~Table() {}

    std::optional<Entry> getRow(std::string_view _key);
    std::vector<std::optional<Entry>> getRows(
        const std::variant<const gsl::span<std::string_view const>,
            const gsl::span<std::string const>>& _keys);
    std::vector<std::string> getPrimaryKeys(const std::optional<const Condition>& _condition);

    void setRow(std::string_view _key, Entry _entry);

    void asyncGetPrimaryKeys(std::optional<const Condition> const& _condition,
        std::function<void(Error::UniquePtr, std::vector<std::string>)> _callback) noexcept;

    void asyncGetRow(std::string_view _key,
        std::function<void(Error::UniquePtr, std::optional<Entry>)> _callback) noexcept;

    void asyncGetRows(const std::variant<const gsl::span<std::string_view const>,
                          const gsl::span<std::string const>>& _keys,
        std::function<void(Error::UniquePtr, std::vector<std::optional<Entry>>)>
            _callback) noexcept;

    void asyncSetRow(
        std::string_view key, Entry entry, std::function<void(Error::UniquePtr)> callback) noexcept;

    TableInfo::ConstPtr tableInfo() const { return m_tableInfo; }
    Entry newEntry() { return Entry(m_tableInfo); }
    Entry newDeletedEntry()
    {
        auto deletedEntry = newEntry();
        deletedEntry.setStatus(Entry::DELETED);
        return deletedEntry;
    }

protected:
    StorageInterface* m_storage;
    TableInfo::ConstPtr m_tableInfo;
};

}  // namespace storage
}  // namespace bcos


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

/* Implicit Euler a.k.a backward Euler method. */

/* Reference: Ascher, U.M., Petzold, L.R., Computer methods for
   ordinary differential and differential-algebraic equations, SIAM,
   Philadelphia, 1998.
*/

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

#include "odeiv_util.h"
#include "rksubs.c"
#include "modnewton1.c"

/* Stage of method */
#define RK1IMP_STAGE 1

typedef struct
{
  gsl_matrix *A;                /* Runge-Kutta coefficients */
  double *y_onestep;            /* Result with one step */
  double *y_twostep;            /* Result with two half steps */
  double *ytmp;                 /* Temporary work space */
  double *y_save;               /* Backup space */
  double *YZ;                   /* Runge-Kutta points */
  double *fYZ;                  /* Derivatives at YZ */
  gsl_matrix *dfdy;             /* Jacobian matrix */
  double *dfdt;                 /* time derivative of f */
  modnewton1_state_t *esol;     /* nonlinear equation solver */
  double *errlev;               /* desired error level of y */
  const gsl_odeiv2_driver *driver;      /* pointer to driver object */
}
rk1imp_state_t;

static void *
rk1imp_alloc (size_t dim)
{
  rk1imp_state_t *state = (rk1imp_state_t *) malloc (sizeof (rk1imp_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for rk1imp_state",
                      GSL_ENOMEM);
    }

  state->A = gsl_matrix_alloc (RK1IMP_STAGE, RK1IMP_STAGE);

  if (state->A == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for A", GSL_ENOMEM);
    }

  state->y_onestep = (double *) malloc (dim * sizeof (double));

  if (state->y_onestep == 0)
    {
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->y_twostep = (double *) malloc (dim * sizeof (double));

  if (state->y_twostep == 0)
    {
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_onestep", GSL_ENOMEM);
    }

  state->ytmp = (double *) malloc (dim * sizeof (double));

  if (state->ytmp == 0)
    {
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for ytmp", GSL_ENOMEM);
    }

  state->y_save = (double *) malloc (dim * sizeof (double));

  if (state->y_save == 0)
    {
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for y_save", GSL_ENOMEM);
    }

  state->YZ = (double *) malloc (dim * RK1IMP_STAGE * sizeof (double));

  if (state->YZ == 0)
    {
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for YZ", GSL_ENOMEM);
    }

  state->fYZ = (double *) malloc (dim * RK1IMP_STAGE * sizeof (double));

  if (state->fYZ == 0)
    {
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for fYZ", GSL_ENOMEM);
    }

  state->dfdt = (double *) malloc (dim * sizeof (double));

  if (state->dfdt == 0)
    {
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdt", GSL_ENOMEM);
    }

  state->dfdy = gsl_matrix_alloc (dim, dim);

  if (state->dfdy == 0)
    {
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for dfdy", GSL_ENOMEM);
    }

  state->esol = modnewton1_alloc (dim, RK1IMP_STAGE);

  if (state->esol == 0)
    {
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for esol", GSL_ENOMEM);
    }

  state->errlev = (double *) malloc (dim * sizeof (double));

  if (state->errlev == 0)
    {
      modnewton1_free (state->esol);
      gsl_matrix_free (state->dfdy);
      free (state->dfdt);
      free (state->fYZ);
      free (state->YZ);
      free (state->y_save);
      free (state->ytmp);
      free (state->y_twostep);
      free (state->y_onestep);
      gsl_matrix_free (state->A);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for errlev", GSL_ENOMEM);
    }

  state->driver = NULL;

  return state;
}

static int
rk1imp_apply (void *vstate, size_t dim, double t, double h,
              double y[], double yerr[],
              const double dydt_in[], double dydt_out[],
              const gsl_odeiv2_system * sys)
{
  /* Makes an implicit Euler step with size h and estimates the local
     error of the step by step doubling.
   */

  rk1imp_state_t *state = (rk1imp_state_t *) vstate;

  double *const y_onestep = state->y_onestep;
  double *const y_twostep = state->y_twostep;
  double *const ytmp = state->ytmp;
  double *const y_save = state->y_save;
  double *const YZ = state->YZ;
  double *const fYZ = state->fYZ;
  gsl_matrix *const dfdy = state->dfdy;
  double *const dfdt = state->dfdt;
  double *const errlev = state->errlev;

  const modnewton1_state_t *esol = state->esol;

  /* Runge-Kutta coefficients */

  gsl_matrix *A = state->A;
  const double b[] = { 1.0 };
  const double c[] = { 1.0 };
  gsl_matrix_set (A, 0, 0, 1.0);

  if (esol == NULL)
    {
      GSL_ERROR ("no non-linear equation solver speficied", GSL_EINVAL);
    }

  /* Get desired error levels via gsl_odeiv2_control object through driver
     object, which is a requirement for this stepper.
   */

  if (state->driver == NULL)
    {
      return GSL_EFAULT;
    }
  else
    {
      size_t i;

      for (i = 0; i < dim; i++)
        {
          if (dydt_in != NULL)
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           dydt_in[i], h, i, &errlev[i]);
            }
          else
            {
              gsl_odeiv2_control_errlevel (state->driver->c, y[i],
                                           0.0, h, i, &errlev[i]);
            }
        }
    }

  /* Evaluate Jacobian for modnewton1 */

  {
    int s = GSL_ODEIV_JA_EVAL (sys, t, y, dfdy->data, dfdt);

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

  /* Calculate a single step with size h */

  {
    int s = modnewton1_init ((void *) esol, A, h, dfdy, sys);

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

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h, y,
                              sys, YZ, errlev);

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

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + c[0] * h, YZ, fYZ);

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

  {
    int s = rksubs (y_onestep, h, y, fYZ, b, RK1IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      return s;
  }

  /* Error estimation by step doubling */

  {
    int s = modnewton1_init ((void *) esol, A, h / 2.0, dfdy, sys);

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

  /* 1st half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t, h / 2.0, y,
                              sys, YZ, errlev);

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

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + c[0] * h / 2.0, YZ, fYZ);
    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    int s = rksubs (ytmp, h / 2.0, y, fYZ, b, RK1IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      return s;
  }

  /* Save original y values in case of error */

  DBL_MEMCPY (y_save, y, dim);

  /* 2nd half step */

  {
    int s = modnewton1_solve ((void *) esol, A, c, t + h / 2.0, h / 2.0,
                              ytmp, sys, YZ, errlev);

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

  {
    int s = GSL_ODEIV_FN_EVAL (sys, t + h / 2.0 + c[0] * h / 2.0, YZ, fYZ);
    if (s != GSL_SUCCESS)
      {
        return s;
      }
  }

  {
    /* Note: rk1imp returns y using the results from two half steps
       instead of the single step since the results are freely
       available and more precise.
     */

    int s = rksubs (y_twostep, h / 2.0, ytmp, fYZ, b, RK1IMP_STAGE, dim);

    if (s != GSL_SUCCESS)
      {
        DBL_MEMCPY (y, y_save, dim);
        return s;
      }
  }

  DBL_MEMCPY (y, y_twostep, dim);

  /* Error estimation */

  {
    size_t i;
    for (i = 0; i < dim; i++)
      {
        yerr[i] = ODEIV_ERR_SAFETY * 0.5 * fabs (y_twostep[i] - y_onestep[i]);
      }
  }

  /* Derivatives at output */

  if (dydt_out != NULL)
    {
      int s = GSL_ODEIV_FN_EVAL (sys, t + h, y, dydt_out);

      if (s != GSL_SUCCESS)
        {
          /* Restore original values */
          DBL_MEMCPY (y, y_save, dim);

          return s;
        }
    }

  return GSL_SUCCESS;
}

static int
rk1imp_set_driver (void *vstate, const gsl_odeiv2_driver * d)
{
  rk1imp_state_t *state = (rk1imp_state_t *) vstate;

  state->driver = d;

  return GSL_SUCCESS;
}

static int
rk1imp_reset (void *vstate, size_t dim)
{
  rk1imp_state_t *state = (rk1imp_state_t *) vstate;

  DBL_ZERO_MEMSET (state->y_onestep, dim);
  DBL_ZERO_MEMSET (state->y_twostep, dim);
  DBL_ZERO_MEMSET (state->ytmp, dim);
  DBL_ZERO_MEMSET (state->y_save, dim);
  DBL_ZERO_MEMSET (state->YZ, dim);
  DBL_ZERO_MEMSET (state->fYZ, dim);

  return GSL_SUCCESS;
}

static unsigned int
rk1imp_order (void *vstate)
{
  rk1imp_state_t *state = (rk1imp_state_t *) vstate;
  state = 0;                    /* prevent warnings about unused parameters */
  return 1;
}

static void
rk1imp_free (void *vstate)
{
  rk1imp_state_t *state = (rk1imp_state_t *) vstate;

  free (state->errlev);
  modnewton1_free (state->esol);
  gsl_matrix_free (state->dfdy);
  free (state->dfdt);
  free (state->fYZ);
  free (state->YZ);
  free (state->y_save);
  free (state->ytmp);
  free (state->y_twostep);
  free (state->y_onestep);
  gsl_matrix_free (state->A);
  free (state);
}

static const gsl_odeiv2_step_type rk1imp_type = {
  "rk1imp",                     /* name */
  1,                            /* can use dydt_in? */
  1,                            /* gives exact dydt_out? */
  &rk1imp_alloc,
  &rk1imp_apply,
  &rk1imp_set_driver,
  &rk1imp_reset,
  &rk1imp_order,
  &rk1imp_free
};

const gsl_odeiv2_step_type *gsl_odeiv2_step_rk1imp = &rk1imp_type;

max_stars_repo_name: ielomariala/Hex-Game


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

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

#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_matrix.h>

/*
 * This module computes the eigenvalues and eigenvectors of a real
 * nonsymmetric matrix.
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license. The LAPACK
 * routines used are DTREVC and DLALN2.
 */

#define GSL_NONSYMMV_SMLNUM (2.0 * GSL_DBL_MIN)
#define GSL_NONSYMMV_BIGNUM ((1.0 - GSL_DBL_EPSILON) / GSL_NONSYMMV_SMLNUM)

static void nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                            gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec,
                                            gsl_eigen_nonsymmv_workspace *w);
static void nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
                                            gsl_matrix_complex *evec);

/*
gsl_eigen_nonsymmv_alloc()

Allocate a workspace for solving the nonsymmetric eigenvalue problem.
The size of this workspace is O(5n).

Inputs: n - size of matrices

Return: pointer to workspace
*/

gsl_eigen_nonsymmv_workspace *
gsl_eigen_nonsymmv_alloc(const size_t n)
{
  gsl_eigen_nonsymmv_workspace *w;

  if (n == 0)
    {
      GSL_ERROR_NULL ("matrix dimension must be positive integer",
                      GSL_EINVAL);
    }

  w = (gsl_eigen_nonsymmv_workspace *)
      calloc (1, sizeof (gsl_eigen_nonsymmv_workspace));

  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->size = n;
  w->Z = NULL;
  w->nonsymm_workspace_p = gsl_eigen_nonsymm_alloc(n);

  if (w->nonsymm_workspace_p == 0)
    {
      gsl_eigen_nonsymmv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for nonsymm workspace", GSL_ENOMEM);
    }

  /*
   * set parameters to compute the full Schur form T and balance
   * the matrices
   */
  gsl_eigen_nonsymm_params(1, 1, w->nonsymm_workspace_p);

  w->work = gsl_vector_alloc(n);
  w->work2 = gsl_vector_alloc(n);
  w->work3 = gsl_vector_alloc(n);
  if (w->work == 0 || w->work2 == 0 || w->work3 == 0)
    {
      gsl_eigen_nonsymmv_free(w);
      GSL_ERROR_NULL ("failed to allocate space for nonsymmv additional workspace", GSL_ENOMEM);
    }

  return (w);
} /* gsl_eigen_nonsymmv_alloc() */

/*
gsl_eigen_nonsymmv_free()
  Free workspace w
*/

void
gsl_eigen_nonsymmv_free (gsl_eigen_nonsymmv_workspace * w)
{
  RETURN_IF_NULL (w);

  if (w->nonsymm_workspace_p)
    gsl_eigen_nonsymm_free(w->nonsymm_workspace_p);

  if (w->work)
    gsl_vector_free(w->work);

  if (w->work2)
    gsl_vector_free(w->work2);

  if (w->work3)
    gsl_vector_free(w->work3);

  free(w);
} /* gsl_eigen_nonsymmv_free() */

/*
gsl_eigen_nonsymmv()

Solve the nonsymmetric eigensystem problem

A x = \lambda x

for the eigenvalues \lambda and right eigenvectors x

Inputs: A    - general real matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        w    - workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv (gsl_matrix * A, gsl_vector_complex * eval,
                    gsl_matrix_complex * evec,
                    gsl_eigen_nonsymmv_workspace * w)
{
  const size_t N = A->size1;

  /* check matrix and vector sizes */

  if (N != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues", GSL_ENOTSQR);
    }
  else if (eval->size != N)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != N)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else
    {
      int s;
      gsl_matrix Z;

      /*
       * We need a place to store the Schur vectors, so we will
       * treat evec as a real matrix and store them in the left
       * half - the factor of 2 in the tda corresponds to the
       * complex multiplicity
       */
      Z.size1 = N;
      Z.size2 = N;
      Z.tda = 2 * N;
      Z.data = evec->data;
      Z.block = 0;
      Z.owner = 0;

      /* compute eigenvalues, Schur form, and Schur vectors */
      s = gsl_eigen_nonsymm_Z(A, eval, &Z, w->nonsymm_workspace_p);

      if (w->Z)
        {
          /*
           * save the Schur vectors in user supplied matrix, since
           * they will be destroyed when computing eigenvectors
           */
          gsl_matrix_memcpy(w->Z, &Z);
        }

      /* only compute eigenvectors if we found all eigenvalues */
      if (s == GSL_SUCCESS)
        {
          /* compute eigenvectors */
          nonsymmv_get_right_eigenvectors(A, &Z, eval, evec, w);

          /* normalize so that Euclidean norm is 1 */
          nonsymmv_normalize_eigenvectors(eval, evec);
        }

      return s;
    }
} /* gsl_eigen_nonsymmv() */

/*
gsl_eigen_nonsymmv_Z()
  Compute eigenvalues and eigenvectors of a real nonsymmetric matrix
and also save the Schur vectors. See comments in gsl_eigen_nonsymm_Z
for more information.

Inputs: A    - real nonsymmetric matrix
        eval - where to store eigenvalues
        evec - where to store eigenvectors
        Z    - where to store Schur vectors
        w    - nonsymmv workspace

Return: success or error
*/

int
gsl_eigen_nonsymmv_Z (gsl_matrix * A, gsl_vector_complex * eval,
                      gsl_matrix_complex * evec, gsl_matrix * Z,
                      gsl_eigen_nonsymmv_workspace * w)
{
  /* check matrix and vector sizes */

  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix must be square to compute eigenvalues/eigenvectors", GSL_ENOTSQR);
    }
  else if (eval->size != A->size1)
    {
      GSL_ERROR ("eigenvalue vector must match matrix size", GSL_EBADLEN);
    }
  else if (evec->size1 != evec->size2)
    {
      GSL_ERROR ("eigenvector matrix must be square", GSL_ENOTSQR);
    }
  else if (evec->size1 != A->size1)
    {
      GSL_ERROR ("eigenvector matrix has wrong size", GSL_EBADLEN);
    }
  else if ((Z->size1 != Z->size2) || (Z->size1 != A->size1))
    {
      GSL_ERROR ("Z matrix has wrong dimensions", GSL_EBADLEN);
    }
  else
    {
      int s;

      w->Z = Z;

      s = gsl_eigen_nonsymmv(A, eval, evec, w);

      w->Z = NULL;

      return s;
    }
} /* gsl_eigen_nonsymmv_Z() */

/********************************************
 *           INTERNAL ROUTINES              *
 ********************************************/

/*
nonsymmv_get_right_eigenvectors()
  Compute the right eigenvectors of the Schur form T and then
backtransform them using the Schur vectors to get right eigenvectors of
the original matrix.

Inputs: T    - Schur form
        Z    - Schur vectors
        eval - where to store eigenvalues (to ensure that the
               correct eigenvalue is stored in the same position
               as the eigenvectors)
        evec - where to store eigenvectors
        w    - nonsymmv workspace

Return: none

Notes: 1) based on LAPACK routine DTREVC - the algorithm used is
          backsubstitution on the upper quasi triangular system T
          followed by backtransformation by Z to get vectors of the
          original matrix.

       2) The Schur vectors in Z are destroyed and replaced with
          eigenvectors stored with the same storage scheme as DTREVC.
          The eigenvectors are also stored in 'evec'

       3) The matrix T is unchanged on output

       4) Each eigenvector is normalized so that the element of
          largest magnitude has magnitude 1; here the magnitude of
          a complex number (x,y) is taken to be |x| + |y|
*/

static void
nonsymmv_get_right_eigenvectors(gsl_matrix *T, gsl_matrix *Z,
                                gsl_vector_complex *eval,
                                gsl_matrix_complex *evec,
                                gsl_eigen_nonsymmv_workspace *w)
{
  const size_t N = T->size1;
  const double smlnum = GSL_DBL_MIN * N / GSL_DBL_EPSILON;
  const double bignum = (1.0 - GSL_DBL_EPSILON) / smlnum;
  int i;              /* looping */
  size_t iu,          /* looping */
         ju,
         ii;
  gsl_complex lambda; /* current eigenvalue */
  double lambda_re,   /* Re(lambda) */
         lambda_im;   /* Im(lambda) */
  gsl_matrix_view Tv, /* temporary views */
                  Zv;
  gsl_vector_view y,  /* temporary views */
                  y2,
                  ev,
                  ev2;
  double dat[4],      /* scratch arrays */
         dat_X[4];
  double scale;       /* scale factor */
  double xnorm;       /* |X| */
  gsl_vector_complex_view ecol, /* column of evec */
                          ecol2;
  int complex_pair;   /* complex eigenvalue pair? */
  double smin;

  /*
   * Compute 1-norm of each column of upper triangular part of T
   * to control overflow in triangular solver
   */

  gsl_vector_set(w->work3, 0, 0.0);
  for (ju = 1; ju < N; ++ju)
    {
      gsl_vector_set(w->work3, ju, 0.0);
      for (iu = 0; iu < ju; ++iu)
        {
          gsl_vector_set(w->work3, ju,
                         gsl_vector_get(w->work3, ju) +
                         fabs(gsl_matrix_get(T, iu, ju)));
        }
    }

  for (i = (int) N - 1; i >= 0; --i)
    {
      iu = (size_t) i;

      /* get current eigenvalue and store it in lambda */
      lambda_re = gsl_matrix_get(T, iu, iu);

      if (iu != 0 && gsl_matrix_get(T, iu, iu - 1) != 0.0)
        {
          lambda_im = sqrt(fabs(gsl_matrix_get(T, iu, iu - 1))) *
                      sqrt(fabs(gsl_matrix_get(T, iu - 1, iu)));
        }
      else
        {
          lambda_im = 0.0;
        }

      GSL_SET_COMPLEX(&lambda, lambda_re, lambda_im);

      smin = GSL_MAX(GSL_DBL_EPSILON * (fabs(lambda_re) + fabs(lambda_im)),
                     smlnum);
      smin = GSL_MAX(smin, GSL_NONSYMMV_SMLNUM);

      if (lambda_im == 0.0)
        {
          int k, l;
          gsl_vector_view bv, xv;

          /* real eigenvector */

          /*
           * The ordering of eigenvalues in 'eval' is arbitrary and
           * does not necessarily follow the Schur form T, so store
           * lambda in the right slot in eval to ensure it corresponds
           * to the eigenvector we are about to compute
           */
          gsl_vector_complex_set(eval, iu, lambda);

          /*
           * We need to solve the system:
           *
           * (T(1:iu-1, 1:iu-1) - lambda*I)*X = -T(1:iu-1,iu)
           */

          /* construct right hand side */
          for (k = 0; k < i; ++k)
            {
              gsl_vector_set(w->work,
                             (size_t) k,
                             -gsl_matrix_get(T, (size_t) k, iu));
            }

          gsl_vector_set(w->work, iu, 1.0);

          for (l = i - 1; l >= 0; --l)
            {
              size_t lu = (size_t) l;

              if (lu == 0)
                complex_pair = 0;
              else
                complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;

              if (!complex_pair)
                {
                  double x;

                  /*
                   * 1-by-1 diagonal block - solve the system:
                   *
                   * (T_{ll} - lambda)*x = -T_{l(iu)}
                   */

                  Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
                  bv = gsl_vector_view_array(dat, 1);
                  gsl_vector_set(&bv.vector, 0,
                                 gsl_vector_get(w->work, lu));
                  xv = gsl_vector_view_array(dat_X, 1);

                  gsl_schur_solve_equation(1.0,
                                           &Tv.matrix,
                                           lambda_re,
                                           1.0,
                                           1.0,
                                           &bv.vector,
                                           &xv.vector,
                                           &scale,
                                           &xnorm,
                                           smin);

                  /* scale x to avoid overflow */
                  x = gsl_vector_get(&xv.vector, 0);
                  if (xnorm > 1.0)
                    {
                      if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
                        {
                          x /= xnorm;
                          scale /= xnorm;
                        }
                    }

                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  gsl_vector_set(w->work, lu, x);

                  if (lu > 0)
                    {
                      gsl_vector_view v1, v2;

                      /* update right hand side */

                      v1 = gsl_matrix_subcolumn(T, lu, 0, lu);
                      v2 = gsl_vector_subvector(w->work, 0, lu);
                      gsl_blas_daxpy(-x, &v1.vector, &v2.vector);
                    } /* if (l > 0) */
                } /* if (!complex_pair) */
              else
                {
                  double x11, x21;

                  /*
                   * 2-by-2 diagonal block
                   */

                  Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
                  bv = gsl_vector_view_array(dat, 2);
                  gsl_vector_set(&bv.vector, 0,
                                 gsl_vector_get(w->work, lu - 1));
                  gsl_vector_set(&bv.vector, 1,
                                 gsl_vector_get(w->work, lu));
                  xv = gsl_vector_view_array(dat_X, 2);

                  gsl_schur_solve_equation(1.0,
                                           &Tv.matrix,
                                           lambda_re,
                                           1.0,
                                           1.0,
                                           &bv.vector,
                                           &xv.vector,
                                           &scale,
                                           &xnorm,
                                           smin);

                  /* scale X(1,1) and X(2,1) to avoid overflow */
                  x11 = gsl_vector_get(&xv.vector, 0);
                  x21 = gsl_vector_get(&xv.vector, 1);

                  if (xnorm > 1.0)
                    {
                      double beta;

                      beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
                                     gsl_vector_get(w->work3, lu));
                      if (beta > bignum / xnorm)
                        {
                          x11 /= xnorm;
                          x21 /= xnorm;
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  gsl_vector_set(w->work, lu - 1, x11);
                  gsl_vector_set(w->work, lu, x21);

                  /* update right hand side */
                  if (lu > 1)
                    {
                      gsl_vector_view v1, v2;

                      v1 = gsl_matrix_subcolumn(T, lu - 1, 0, lu - 1);
                      v2 = gsl_vector_subvector(w->work, 0, lu - 1);
                      gsl_blas_daxpy(-x11, &v1.vector, &v2.vector);

                      v1 = gsl_matrix_subcolumn(T, lu, 0, lu - 1);
                      gsl_blas_daxpy(-x21, &v1.vector, &v2.vector);
                    }

                  --l;
                } /* if (complex_pair) */
            } /* for (l = i - 1; l >= 0; --l) */

          /*
           * At this point, w->work is an eigenvector of the
           * Schur form T. To get an eigenvector of the original
           * matrix, we multiply on the left by Z, the matrix of
           * Schur vectors
           */

          ecol = gsl_matrix_complex_column(evec, iu);
          y = gsl_matrix_column(Z, iu);

          if (iu > 0)
            {
              gsl_vector_view x;

              Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu);

              x = gsl_vector_subvector(w->work, 0, iu);

              /* compute Z * w->work and store it in Z(:,iu) */
              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work, iu),
                             &y.vector);
            } /* if (iu > 0) */

          /* store eigenvector into evec */

          ev = gsl_vector_complex_real(&ecol.vector);
          ev2 = gsl_vector_complex_imag(&ecol.vector);

          scale = 0.0;
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y.vector, ii);

              /* store real part of eigenvector */
              gsl_vector_set(&ev.vector, ii, a);

              /* set imaginary part to 0 */
              gsl_vector_set(&ev2.vector, ii, 0.0);

              if (fabs(a) > scale)
                scale = fabs(a);
            }

          if (scale != 0.0)
            scale = 1.0 / scale;

          /* scale by magnitude of largest element */
          gsl_blas_dscal(scale, &ev.vector);
        } /* if (GSL_IMAG(lambda) == 0.0) */
      else
        {
          gsl_vector_complex_view bv, xv;
          size_t k;
          int l;
          gsl_complex lambda2;

          /* complex eigenvector */

          /*
           * Store the complex conjugate eigenvalues in the right
           * slots in eval
           */
          GSL_SET_REAL(&lambda2, GSL_REAL(lambda));
          GSL_SET_IMAG(&lambda2, -GSL_IMAG(lambda));
          gsl_vector_complex_set(eval, iu - 1, lambda);
          gsl_vector_complex_set(eval, iu, lambda2);

          /*
           * First solve:
           *
           * [ T(i:i+1,i:i+1) - lambda*I ] * X = 0
           */

          if (fabs(gsl_matrix_get(T, iu - 1, iu)) >=
              fabs(gsl_matrix_get(T, iu, iu - 1)))
            {
              gsl_vector_set(w->work, iu - 1, 1.0);
              gsl_vector_set(w->work2, iu,
                             lambda_im / gsl_matrix_get(T, iu - 1, iu));
            }
          else
            {
              gsl_vector_set(w->work, iu - 1,
                             -lambda_im / gsl_matrix_get(T, iu, iu - 1));
              gsl_vector_set(w->work2, iu, 1.0);
            }
          gsl_vector_set(w->work, iu, 0.0);
          gsl_vector_set(w->work2, iu - 1, 0.0);

          /* construct right hand side */
          for (k = 0; k < iu - 1; ++k)
            {
              gsl_vector_set(w->work, k,
                             -gsl_vector_get(w->work, iu - 1) *
                             gsl_matrix_get(T, k, iu - 1));
              gsl_vector_set(w->work2, k,
                             -gsl_vector_get(w->work2, iu) *
                             gsl_matrix_get(T, k, iu));
            }

          /*
           * We must solve the upper quasi-triangular system:
           *
           * [ T(1:i-2,1:i-2) - lambda*I ] * X = s*(work + i*work2)
           */

          for (l = i - 2; l >= 0; --l)
            {
              size_t lu = (size_t) l;

              if (lu == 0)
                complex_pair = 0;
              else
                complex_pair = gsl_matrix_get(T, lu, lu - 1) != 0.0;

              if (!complex_pair)
                {
                  gsl_complex bval;
                  gsl_complex x;

                  /*
                   * 1-by-1 diagonal block - solve the system:
                   *
                   * (T_{ll} - lambda)*x = work + i*work2
                   */

                  Tv = gsl_matrix_submatrix(T, lu, lu, 1, 1);
                  bv = gsl_vector_complex_view_array(dat, 1);
                  xv = gsl_vector_complex_view_array(dat_X, 1);

                  GSL_SET_COMPLEX(&bval,
                                  gsl_vector_get(w->work, lu),
                                  gsl_vector_get(w->work2, lu));
                  gsl_vector_complex_set(&bv.vector, 0, bval);

                  gsl_schur_solve_equation_z(1.0,
                                             &Tv.matrix,
                                             &lambda,
                                             1.0,
                                             1.0,
                                             &bv.vector,
                                             &xv.vector,
                                             &scale,
                                             &xnorm,
                                             smin);

                  if (xnorm > 1.0)
                    {
                      if (gsl_vector_get(w->work3, lu) > bignum / xnorm)
                        {
                          gsl_blas_zdscal(1.0/xnorm, &xv.vector);
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                      wv = gsl_vector_subvector(w->work2, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }

                  x = gsl_vector_complex_get(&xv.vector, 0);
                  gsl_vector_set(w->work, lu, GSL_REAL(x));
                  gsl_vector_set(w->work2, lu, GSL_IMAG(x));

                  /* update the right hand side */
                  if (lu > 0)
                    {
                      gsl_vector_view v1, v2;

                      v1 = gsl_matrix_subcolumn(T, lu, 0, lu);
                      v2 = gsl_vector_subvector(w->work, 0, lu);
                      gsl_blas_daxpy(-GSL_REAL(x), &v1.vector, &v2.vector);

                      v2 = gsl_vector_subvector(w->work2, 0, lu);
                      gsl_blas_daxpy(-GSL_IMAG(x), &v1.vector, &v2.vector);
                    } /* if (lu > 0) */
                } /* if (!complex_pair) */
              else
                {
                  gsl_complex b1, b2, x1, x2;

                  /*
                   * 2-by-2 diagonal block - solve the system
                   */

                  Tv = gsl_matrix_submatrix(T, lu - 1, lu - 1, 2, 2);
                  bv = gsl_vector_complex_view_array(dat, 2);
                  xv = gsl_vector_complex_view_array(dat_X, 2);

                  GSL_SET_COMPLEX(&b1,
                                  gsl_vector_get(w->work, lu - 1),
                                  gsl_vector_get(w->work2, lu - 1));
                  GSL_SET_COMPLEX(&b2,
                                  gsl_vector_get(w->work, lu),
                                  gsl_vector_get(w->work2, lu));
                  gsl_vector_complex_set(&bv.vector, 0, b1);
                  gsl_vector_complex_set(&bv.vector, 1, b2);

                  gsl_schur_solve_equation_z(1.0,
                                             &Tv.matrix,
                                             &lambda,
                                             1.0,
                                             1.0,
                                             &bv.vector,
                                             &xv.vector,
                                             &scale,
                                             &xnorm,
                                             smin);

                  x1 = gsl_vector_complex_get(&xv.vector, 0);
                  x2 = gsl_vector_complex_get(&xv.vector, 1);

                  if (xnorm > 1.0)
                    {
                      double beta;

                      beta = GSL_MAX(gsl_vector_get(w->work3, lu - 1),
                                     gsl_vector_get(w->work3, lu));
                      if (beta > bignum / xnorm)
                        {
                          gsl_blas_zdscal(1.0/xnorm, &xv.vector);
                          scale /= xnorm;
                        }
                    }

                  /* scale if necessary */
                  if (scale != 1.0)
                    {
                      gsl_vector_view wv;

                      wv = gsl_vector_subvector(w->work, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                      wv = gsl_vector_subvector(w->work2, 0, iu + 1);
                      gsl_blas_dscal(scale, &wv.vector);
                    }
                  gsl_vector_set(w->work, lu - 1, GSL_REAL(x1));
                  gsl_vector_set(w->work, lu, GSL_REAL(x2));
                  gsl_vector_set(w->work2, lu - 1, GSL_IMAG(x1));
                  gsl_vector_set(w->work2, lu, GSL_IMAG(x2));

                  /* update right hand side */
                  if (lu > 1)
                    {
                      gsl_vector_view v1, v2, v3, v4;

                      v1 = gsl_matrix_subcolumn(T, lu - 1, 0, lu - 1);
                      v4 = gsl_matrix_subcolumn(T, lu, 0, lu - 1);
                      v2 = gsl_vector_subvector(w->work, 0, lu - 1);
                      v3 = gsl_vector_subvector(w->work2, 0, lu - 1);

                      gsl_blas_daxpy(-GSL_REAL(x1), &v1.vector, &v2.vector);
                      gsl_blas_daxpy(-GSL_REAL(x2), &v4.vector, &v2.vector);
                      gsl_blas_daxpy(-GSL_IMAG(x1), &v1.vector, &v3.vector);
                      gsl_blas_daxpy(-GSL_IMAG(x2), &v4.vector, &v3.vector);
                    } /* if (lu > 1) */

                  --l;
                } /* if (complex_pair) */
            } /* for (l = i - 2; l >= 0; --l) */

          /*
           * At this point, work + i*work2 is an eigenvector
           * of T - backtransform to get an eigenvector of the
           * original matrix
           */

          y = gsl_matrix_column(Z, iu - 1);
          y2 = gsl_matrix_column(Z, iu);

          if (iu > 1)
            {
              gsl_vector_view x;

              /* compute real part of eigenvectors */

              Zv = gsl_matrix_submatrix(Z, 0, 0, N, iu - 1);
              x = gsl_vector_subvector(w->work, 0, iu - 1);

              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work, iu - 1),
                             &y.vector);

              /* now compute the imaginary part */
              x = gsl_vector_subvector(w->work2, 0, iu - 1);

              gsl_blas_dgemv(CblasNoTrans,
                             1.0,
                             &Zv.matrix,
                             &x.vector,
                             gsl_vector_get(w->work2, iu),
                             &y2.vector);
            }
          else
            {
              gsl_blas_dscal(gsl_vector_get(w->work, iu - 1), &y.vector);
              gsl_blas_dscal(gsl_vector_get(w->work2, iu), &y2.vector);
            }

          /*
           * Now store the eigenvectors into evec - the real parts
           * are Z(:,iu - 1) and the imaginary parts are
           * +/- Z(:,iu)
           */

          /* get views of the two eigenvector slots */
          ecol = gsl_matrix_complex_column(evec, iu - 1);
          ecol2 = gsl_matrix_complex_column(evec, iu);

          /*
           * save imaginary part first as it may get overwritten
           * when copying the real part due to our storage scheme
           * in Z/evec
           */
          ev = gsl_vector_complex_imag(&ecol.vector);
          ev2 = gsl_vector_complex_imag(&ecol2.vector);
          scale = 0.0;
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y2.vector, ii);

              scale = GSL_MAX(scale,
                              fabs(a) + fabs(gsl_vector_get(&y.vector, ii)));

              gsl_vector_set(&ev.vector, ii, a);
              gsl_vector_set(&ev2.vector, ii, -a);
            }

          /* now save the real part */
          ev = gsl_vector_complex_real(&ecol.vector);
          ev2 = gsl_vector_complex_real(&ecol2.vector);
          for (ii = 0; ii < N; ++ii)
            {
              double a = gsl_vector_get(&y.vector, ii);

              gsl_vector_set(&ev.vector, ii, a);
              gsl_vector_set(&ev2.vector, ii, a);
            }

          if (scale != 0.0)
            scale = 1.0 / scale;

          /* scale by largest element magnitude */

          gsl_blas_zdscal(scale, &ecol.vector);
          gsl_blas_zdscal(scale, &ecol2.vector);

          /*
           * decrement i since we took care of two eigenvalues at
           * the same time
           */
          --i;
        } /* if (GSL_IMAG(lambda) != 0.0) */
    } /* for (i = (int) N - 1; i >= 0; --i) */
} /* nonsymmv_get_right_eigenvectors() */

/*
nonsymmv_normalize_eigenvectors()
  Normalize eigenvectors so that their Euclidean norm is 1

Inputs: eval - eigenvalues
        evec - eigenvectors
*/

static void
nonsymmv_normalize_eigenvectors(gsl_vector_complex *eval,
                                gsl_matrix_complex *evec)
{
  const size_t N = evec->size1;
  size_t i;     /* looping */
  gsl_complex ei;
  gsl_vector_complex_view vi;
  gsl_vector_view re, im;
  double scale; /* scaling factor */

  for (i = 0; i < N; ++i)
    {
      ei = gsl_vector_complex_get(eval, i);
      vi = gsl_matrix_complex_column(evec, i);

      re = gsl_vector_complex_real(&vi.vector);

      if (GSL_IMAG(ei) == 0.0)
        {
          scale = 1.0 / gsl_blas_dnrm2(&re.vector);
          gsl_blas_dscal(scale, &re.vector);
        }
      else if (GSL_IMAG(ei) > 0.0)
        {
          im = gsl_vector_complex_imag(&vi.vector);

          scale = 1.0 / gsl_hypot(gsl_blas_dnrm2(&re.vector),
                                  gsl_blas_dnrm2(&im.vector));
          gsl_blas_zdscal(scale, &vi.vector);

          vi = gsl_matrix_complex_column(evec, i + 1);
          gsl_blas_zdscal(scale, &vi.vector);
        }
    }
} /* nonsymmv_normalize_eigenvectors() */


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_odeiv2.h>

BC_INLINE4(GSL_ODEIV_FN_EVAL,gsl_odeiv2_system*,double,double*,double*,int)
BC_INLINE5(GSL_ODEIV_JA_EVAL,gsl_odeiv2_system*,double,double*,double*,double*,int)


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

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_spmatrix.h>

/*
create_random_sparse()
  Create a random sparse matrix with approximately
M*N*density non-zero entries

Inputs: M       - number of rows
        N       - number of columns
        density - sparse density \in [0,1]
                  0 = no non-zero entries
                  1 = all m*n entries are filled
        r       - random number generator

Return: pointer to sparse matrix in triplet format (must be freed by caller)

Notes:
1) non-zero matrix entries are uniformly distributed in [0,1]
*/

static gsl_spmatrix *
create_random_sparse(const size_t M, const size_t N, const double density,
                     const gsl_rng *r)
{
  size_t nnzwanted = (size_t) floor(M * N * GSL_MIN(density, 1.0));
  gsl_spmatrix *m = gsl_spmatrix_alloc_nzmax(M, N,
                                             nnzwanted,
                                             GSL_SPMATRIX_TRIPLET);

  while (gsl_spmatrix_nnz(m) < nnzwanted)
    {
      /* generate a random row and column */
      size_t i = gsl_rng_uniform(r) * M;
      size_t j = gsl_rng_uniform(r) * N;

      /* generate random m_{ij} and add it */
      double x = gsl_rng_uniform(r);
      gsl_spmatrix_set(m, i, j, x);
    }

  return m;
} /* create_random_sparse() */

static void
test_getset(const size_t M, const size_t N, const gsl_rng *r)
{
  int status;
  size_t i, j;

  /* test triplet versions of _get and _set */
  {
    size_t k = 0;
    gsl_spmatrix *m = gsl_spmatrix_alloc(M, N);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double x = (double) ++k;
            double y;

            gsl_spmatrix_set(m, i, j, x);
            y = gsl_spmatrix_get(m, i, j);
            if (x != y)
              status = 1;
          }
      }

    gsl_test(status, "test_getset: M=%zu N=%zu _get != _set", M, N);

    /* test setting an element to 0 */
    gsl_spmatrix_set(m, 0, 0, 1.0);
    gsl_spmatrix_set(m, 0, 0, 0.0);

    status = gsl_spmatrix_get(m, 0, 0) != 0.0;
    gsl_test(status, "test_getset: M=%zu N=%zu m(0,0) = %f",
             M, N, gsl_spmatrix_get(m, 0, 0));

    /* test gsl_spmatrix_set_zero() */
    gsl_spmatrix_set(m, 0, 0, 1.0);
    gsl_spmatrix_set_zero(m);
    status = gsl_spmatrix_get(m, 0, 0) != 0.0;
    gsl_test(status, "test_getset: M=%zu N=%zu set_zero m(0,0) = %f",
             M, N, gsl_spmatrix_get(m, 0, 0));

    /* resassemble matrix to ensure nz is calculated correctly */
    k = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double x = (double) ++k;
            gsl_spmatrix_set(m, i, j, x);
          }
      }

    status = gsl_spmatrix_nnz(m) != M * N;
    gsl_test(status, "test_getset: M=%zu N=%zu set_zero nz = %zu",
             M, N, gsl_spmatrix_nnz(m));

    gsl_spmatrix_free(m);
  }

  /* test duplicate values are handled correctly */
  {
    size_t min = GSL_MIN(M, N);
    size_t expected_nnz = min;
    size_t nnz;
    size_t k = 0;
    gsl_spmatrix *m = gsl_spmatrix_alloc(M, N);

    status = 0;
    for (i = 0; i < min; ++i)
      {
        for (j = 0; j < 5; ++j)
          {
            double x = (double) ++k;
            double y;

            gsl_spmatrix_set(m, i, i, x);
            y = gsl_spmatrix_get(m, i, i);
            if (x != y)
              status = 1;
          }
      }

    gsl_test(status, "test_getset: duplicate test M=%zu N=%zu _get != _set", M, N);

    nnz = gsl_spmatrix_nnz(m);
    status = nnz != expected_nnz;
    gsl_test(status, "test_getset: duplicate test M=%zu N=%zu nnz=%zu, expected=%zu",
             M, N, nnz, expected_nnz);

    gsl_spmatrix_free(m);
  }

  /* test compressed version of gsl_spmatrix_get() */
  {
    gsl_spmatrix *T = create_random_sparse(M, N, 0.3, r);
    gsl_spmatrix *C = gsl_spmatrix_compcol(T);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Tij = gsl_spmatrix_get(T, i, j);
            double Cij = gsl_spmatrix_get(C, i, j);

            if (Tij != Cij)
              status = 1;
          }
      }

    gsl_test(status, "test_getset: M=%zu N=%zu compressed _get", M, N);

    gsl_spmatrix_free(T);
    gsl_spmatrix_free(C);
  }
} /* test_getset() */

static void
test_memcpy(const size_t M, const size_t N, const gsl_rng *r)
{
  int status;

  {
    gsl_spmatrix *at = create_random_sparse(M, N, 0.2, r);
    gsl_spmatrix *ac = gsl_spmatrix_compcol(at);
    gsl_spmatrix *bt, *bc;

    bt = gsl_spmatrix_alloc(M, N);
    gsl_spmatrix_memcpy(bt, at);

    status = gsl_spmatrix_equal(at, bt) != 1;
    gsl_test(status, "test_memcpy: _memcpy M=%zu N=%zu triplet format", M, N);

    bc = gsl_spmatrix_alloc_nzmax(M, N, ac->nzmax, GSL_SPMATRIX_CCS);
    gsl_spmatrix_memcpy(bc, ac);

    status = gsl_spmatrix_equal(ac, bc) != 1;
    gsl_test(status, "test_memcpy: _memcpy M=%zu N=%zu compressed column format", M, N);

    gsl_spmatrix_free(at);
    gsl_spmatrix_free(ac);
    gsl_spmatrix_free(bt);
    gsl_spmatrix_free(bc);
  }

  /* test transpose_memcpy */
  {
    gsl_spmatrix *A = create_random_sparse(M, N, 0.3, r);
    gsl_spmatrix *B = gsl_spmatrix_compcol(A);
    gsl_spmatrix *AT = gsl_spmatrix_alloc(N, M);
    gsl_spmatrix *BT = gsl_spmatrix_alloc_nzmax(N, M, 1, GSL_SPMATRIX_CCS);
    size_t i, j;

    gsl_spmatrix_transpose_memcpy(AT, A);
    gsl_spmatrix_transpose_memcpy(BT, B);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double Aij = gsl_spmatrix_get(A, i, j);
            double ATji = gsl_spmatrix_get(AT, j, i);
            double Bij = gsl_spmatrix_get(B, i, j);
            double BTji = gsl_spmatrix_get(BT, j, i);

            if ((Aij != ATji) || (Bij != BTji) || (Aij != Bij))
              status = 1;
          }
      }

    gsl_test(status, "test_memcpy: _transpose_memcpy M=%zu N=%zu triplet format", M, N);

    gsl_spmatrix_free(A);
    gsl_spmatrix_free(AT);
    gsl_spmatrix_free(B);
    gsl_spmatrix_free(BT);
  }
} /* test_memcpy() */

static void
test_ops(const size_t M, const size_t N, const gsl_rng *r)
{
  size_t i, j;
  int status;

  /* test gsl_spmatrix_add */
  {
    gsl_spmatrix *Ta = create_random_sparse(M, N, 0.2, r);
    gsl_spmatrix *Tb = create_random_sparse(M, N, 0.2, r);
    gsl_spmatrix *a = gsl_spmatrix_compcol(Ta);
    gsl_spmatrix *b = gsl_spmatrix_compcol(Tb);
    gsl_spmatrix *c = gsl_spmatrix_alloc_nzmax(M, N, 1, GSL_SPMATRIX_CCS);

    gsl_spmatrix_add(c, a, b);

    status = 0;
    for (i = 0; i < M; ++i)
      {
        for (j = 0; j < N; ++j)
          {
            double aij = gsl_spmatrix_get(a, i, j);
            double bij = gsl_spmatrix_get(b, i, j);
            double cij = gsl_spmatrix_get(c, i, j);

            if (aij + bij != cij)
              status = 1;
          }
      }

    gsl_test(status, "test_ops: _add M=%zu N=%zu compressed format", M, N);

    gsl_spmatrix_free(Ta);
    gsl_spmatrix_free(Tb);
    gsl_spmatrix_free(a);
    gsl_spmatrix_free(b);
    gsl_spmatrix_free(c);
  }
} /* test_ops() */

int
main()
{
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  test_memcpy(10, 10, r);
  test_memcpy(10, 15, r);
  test_memcpy(53, 213, r);
  test_memcpy(920, 2, r);
  test_memcpy(2, 920, r);

  test_getset(20, 20, r);
  test_getset(30, 20, r);
  test_getset(15, 210, r);

  test_ops(20, 20, r);
  test_ops(50, 20, r);
  test_ops(20, 50, r);
  test_ops(76, 43, r);

  gsl_rng_free(r);

  exit (gsl_test_summary());
} /* main() */

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc