text:
/**
 *
 * @file example_dpotrf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

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

static int check_factorization(int, int, double*, double*, int, int*);

int testing_dgetrf(int argc, char **argv)
{
	int M	  = 1000;
    int N     = 500;
    int LDA   = 1000;
    int info_factorization;

    double *A1   = (double *)malloc(LDA*N*sizeof(double));
    double *A2   = (double *)malloc(LDA*N*sizeof(double));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return 0;
    }

    /* Initialize A1 and A2 */
	LAPACKE_dlarnv(3,ISEED, LDA*N, A1);

	int i;
	for(i = 0; i < N*LDA; ++i){
		A2[i] = A1[i];
	}
    /* Plasma routines */
	int *IPIV = calloc(LDA, sizeof(int));
    PLASMA_dgetrf(M, N, A2, LDA, IPIV);

    /* Check the factorization */
    info_factorization = check_factorization( M, N, A1, A2, LDA, IPIV);

    if ( info_factorization != 0 )
       printf("-- Error in DGETRF example ! \n");
    else
       printf("-- Run of DGETRF example successful ! \n");

    free(A1); free(A2);

    return 0;
}

static int check_factorization(int M, int N, double *A1, double *A2, int LDA, int *IPIV)
{
    double Anorm, Rnorm;
	double Anorm1, Rnorm1;
    double alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

//    double *Residual = (double *)malloc(M*N*sizeof(double));
//    double *L1       = (double *)malloc(M*N*sizeof(double));
//    double *L2       = (double *)malloc(M*N*sizeof(double));
    double *work              = (double *)malloc(N*sizeof(double));
	int *ipiv0 = calloc(LDA, sizeof(int));
	int dff = 0;

	LAPACKE_dgetrf(LAPACK_COL_MAJOR, M, N, A1, M, ipiv0);

	for(j=0; j<N; j++)
		dff += IPIV[j] - ipiv0[j];
	printf("dff: %d\n", dff);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A1, M, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A2, LDA, work);

	Rnorm1 = PLASMA_dlange(PlasmaInfNorm, M, N, A1, LDA);
	Anorm1 = PLASMA_dlange(PlasmaInfNorm, M, N, A2, LDA);

	printf("|Rnorm-Rnorm1|: %e, |Anorm-Anorm1|: %e\n", fabs(Rnorm-Rnorm1), fabs(Anorm-Anorm1));
    printf("============\n");
    printf("Checking the LU Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",fabs(Rnorm-Anorm)/(Anorm));

    if ( isnan(fabs(Rnorm-Anorm)/(Anorm)) || (fabs(Rnorm-Anorm)/(Anorm) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

//    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}


text:
#pragma once
#include <gsl/gsl_fit.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_vector.h>
#include <cstdlib>
#include "DataFrame.h"
#include "Ops.h"

namespace adsl {

	// Evaluate the result of a fit at a given point based on the DataFrame description
	// double <- DataFrame
	/*
	auto evalFit = [](double t) {
		auto retFunc = [t](DataFrame& df) {
			if (df.getDesc() == "gsl_fit_linear") {
				double m = df + select({ "slope" }) + getFirst + single();
				double b = df + select({ "intercept" }) + getFirst + single();
				return (m * t) + b;
			}
			else {
				return (double)NAN;
			}
		};
		return retFunc;
	};
	*/

	// Perform a linear fit using gsl_fit_linear
	// df must have only 2 columns
	// DataFrame <- DataFrame
	auto fitLinear = [](DataFrame& df) {
		if (df.getCols() == 2) {
			if (df.getData()[0].type != DBL || df.getData()[1].type != DBL) {
				std::cout << "[fitLinear] <<ERROR>> At least one of the columns is not of type DBL " << std::endl;
				exit(1);
				DataFrame empty_df;
				return empty_df;
			}
			// Setup
			DataFrame ret;
			auto xData = df.getData()[0].toVec_dbl();
			auto yData = df.getData()[1].toVec_dbl();
			size_t length = xData.size();
			double* x;
			double* y;
			x = (double*)malloc(length * sizeof(double));
			y = (double*)malloc(length * sizeof(double));
			std::copy(xData.begin(), xData.end(), x);
			std::copy(yData.begin(), yData.end(), y);

			// Fitting & r-squared
			double c0, c1, cov00, cov01, cov11, sumsq;
			gsl_fit_linear(x, 1, y, 1, length, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
			vd tmpSlope = {c1};
			DataList slope(&tmpSlope, DataType::DBL, "slope");
			vd tmpIntercept = {c0};
			DataList intercept(&tmpIntercept, DataType::DBL, "intercept");
			vd tmpRsqrd = { gsl_stats_correlation(x, 1, y, 1, length) };
			DataList rsqrd(&tmpRsqrd, DataType::DBL, "r-squared");

			// Clean up
			free(x);
			free(y);

			// Return
			ret.addCol(slope);
			ret.addCol(intercept);
			ret.addCol(rsqrd);
			ret.setDesc("gsl_fit_linear");
			return ret;
		}
		else {
			std::cout << "[fitLinear] <<ERROR>> DataFrame doesn't have two columns; has " << df.getCols() << std::endl;
			exit(1);
			DataFrame empty_df;
			return empty_df;
		}
	};

	// TODO: add more fit types

	// Calculate the moving average of a DataList
	// DataList <- DataList <- double
	/*
	auto SMA = [](double period) {
		auto retFunc = [period](DataList& dl) {
			DataFrame ret;
			double* input;
			double* output;
			// TODO
		};
		return retFunc;
	};
	*/

 }

ext: h

max_stars_repo_name: ianfr/adsl-cpp

max_issues_repo_name: ianfr/adsl-cpp

max_forks_repo_name: ianfr/adsl-cpp


text:
/**
 * @file bblas_ztesting.c
 *
 * @brief BBLAS testing  for double _Complex routines.
 *
 * BBLAS is a software package provided by Univ. of Manchester,
 * Univ. of Tennessee.
 *
 * @version 1.0.0
 * @author  Samuel  D. Relton
 * @author  Pedro   V. Lara
 * @author  Mawussi Zounon
 * @date    2016-02-20
 *
 */

#ifndef DOXYGEN_SHOULD_SKIP_THIS
/**
 * Code generation
 * @precisions normal z -> c d s
 **/
#endif

#include "bblas_common.h"
#if defined(BBLAS_WITH_MKL)
#include <mkl_lapacke.h>
#else
#include <lapacke.h>
#endif
#include <cblas.h>

#define COMPLEX

/**
 * Initialize test parameters to their default values.
 **/

void bblas_zinit_config (bblas_ztest_t *test)
{
    test->gen_uplo       =1;
    test->gen_transA     =1;
    test->gen_transB     =1;
    test->gen_trans      =1;
    test->gen_side       =1;
    test->gen_diag       =1;
    test->minM           = 0;
    test->minN           = 0;
    test->minK           = 0;
    test->maxM           = 0;
    test->maxN           = 0;
    test->maxK           = 0;
    test->minbatch_count = 1;
    test->maxbatch_count = 1;
    test->batch_opts     = 0;
    test->routine        = 1;
    test->nb_test        = 1;
    test->set_error      = 0;
    test->global_error   = 0;
    test->faulty_iter    = 0;
    test->mkl_sequential = 0;
    test->new_accuracy   = 1;
}

/**
 * Set the values of all the BBLAS parameters inside the test structure.
 **/

void bblas_zsettest(bblas_ztest_t *test)
{

    enum BBLAS_ROUTINE routine = test->routine;

    /*
     * Set the value of batch count
     */
    bblas_zset_batch_count(test);

    /*
     * Allocate memory and set values for uplo
     */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K)||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	bblas_zsetuplo(test);
    }

    /*
     * Allocate memory and set values for transA
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_TRMM) ||
	(routine == BBLAS_TRSM))
    {
	bblas_zsettransA(test);
    }

    /*
     * Allocate memory and set values for transB
     */
    if (routine == BBLAS_GEMM)
    {
	bblas_zsettransB(test);
    }

    /*
     * Allocate memory and set values for trans
     */
    if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
    {
	bblas_zsettrans(test);
    }

    /*
     * Allocate memory and set values for side
     */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	bblas_zsetside(test);
    }

    /*
     * Allocate memory and set values for diag
     */
    if  ((routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	bblas_zsetdiag(test);
    }

    /*
     * Allocate memory and set values for  M
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	(routine == BBLAS_TRSM))
    {
	bblas_zsetM(test);
    }

    /*
     * Allocate memory and set values for N, all routines
     */
    bblas_zsetN(test);

    /*
     * Allocate memory and set values for K
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K))
    {
	bblas_zsetK(test);
    }

    /*
     * Allocate memory and set values for lda, all routines
     */
    bblas_zsetlda(test);

    /*
     * Allocate memory and set values for ldb, all routines
     */
    if ((routine == BBLAS_GEMM)  || (routine == BBLAS_SYMM)  ||
	(routine == BBLAS_HEMM)  || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K) || (routine == BBLAS_TRMM)  ||
	(routine == BBLAS_TRSM))
    {
	bblas_zsetldb(test);
    }

    /*
     * Allocate memory and set values for ldc, all routines
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	(routine == BBLAS_HER2K))
    {
	bblas_zsetldc(test);
    }

    /*
     * Allocate memory and set values for alpha, all routines
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK)  ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K)||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	bblas_zsetalpha(test);
    }

    /*
     * Allocate memory and set values for alpha
     */
    if (routine == BBLAS_HERK)
    {
	bblas_zsetalpha_herk(test);
    }

    /*
     * Allocate memory and set values for beta
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
        (routine == BBLAS_SYR2K))
    {
	bblas_zsetbeta(test);
    }

    /*
     * Allocate memory and set values for beta_herk
     */
    if ((routine == BBLAS_HERK) || (routine == BBLAS_HER2K))
    {
	bblas_zsetbeta_herk(test);
    }

    /*
     * Allocate memory and set values for arrayA
     */
    bblas_zsetarrayA(test);

    /*
     * Allocate memory and set values for arrayB
     */
    if ((routine == BBLAS_GEMM)  || (routine == BBLAS_SYMM)  ||
	(routine == BBLAS_HEMM)  || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K) || (routine == BBLAS_TRMM)  ||
	(routine == BBLAS_TRSM))
    {
	bblas_zsetarrayB(test);
    }

    /*
     * Allocate memory and set values for arrayC
     */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	(routine == BBLAS_HER2K))
    {
	bblas_zsetarrayC(test);
    }

    /* Memory allocation for result and error variables */
    bblas_zmalloc_result_error(test);
}

/**
 * Allocate memory and set the values of uplo
 **/

void bblas_zsetuplo(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "uplo";
    int random_number;
    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation */
    test->uplo =  (enum BBLAS_UPLO*) malloc(nb_data*sizeof(enum BBLAS_UPLO));

    /*Malloc checking */
    bblas_malloc_check(test->uplo, ptr_name);

    /*set UPLO values */
    switch (test->gen_uplo)
    {
	case UPLO_LOWER:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->uplo[batch_iter] = BblasLower;
	    }
	    break;

	case UPLO_UPPER:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->uplo[batch_iter] = BblasUpper;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50 )
		{
		    test->uplo[batch_iter] = BblasLower;

		}else
		{
		    test->uplo[batch_iter] = BblasUpper;
		}
		break;
	    }
    }
}

/**
 * Allocate memory and set the values of transA
 **/

void bblas_zsettransA(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "transA";
    int random_number;

    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation */
    test->transA =  (enum BBLAS_TRANS*) malloc(nb_data*sizeof(enum BBLAS_TRANS));

    /*Malloc checking */
    bblas_malloc_check(test->transA, ptr_name);

    /*set transA */
    switch (test->gen_transA)
    {
	case NO_TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transA[batch_iter] = BblasNoTrans;
	    }
	    break;

	case TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transA[batch_iter] = BblasTrans;
	    }
	    break;

	case CONJ:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transA[batch_iter] = BblasConjTrans;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50)
		{
		    test->transA[batch_iter] = BblasNoTrans;

		}else if (random_number < 80)
		{
		    test->transA[batch_iter] = BblasTrans;
		}else
		{
		    test->transA[batch_iter] = BblasConjTrans;
		}
	    }
	    break;
    }
}

/**
 * Allocate memory and set the values of transB
 **/

void bblas_zsettransB(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "transB";
    int random_number;

    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation for transB */
    test->transB =  (enum BBLAS_TRANS*) malloc(nb_data*sizeof(enum BBLAS_TRANS));

    /*checking memory allocation */
    bblas_malloc_check(test->transB, ptr_name);

    /*set transB */
    switch (test->gen_transB)
    {
	case NO_TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transB[batch_iter] = BblasNoTrans;
	    }
	    break;

	case TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transB[batch_iter] = BblasTrans;
	    }
	    break;

	case CONJ:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->transB[batch_iter] = BblasConjTrans;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50)
		{
		    test->transB[batch_iter] = BblasNoTrans;

		}else if (random_number < 80)
		{
		    test->transB[batch_iter] = BblasTrans;
		}else
		{
		    test->transB[batch_iter] = BblasConjTrans;
		}
	    }
	    break;
    }

}

/**
 * Allocate memory and set the values of trans
 **/

void bblas_zsettrans(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "trans";
    int random_number;

    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation for trans */
    test->trans =  (enum BBLAS_TRANS*) malloc(nb_data*sizeof(enum BBLAS_TRANS));

    /*checking memory allocation */
    bblas_malloc_check(test->trans, ptr_name);

    /*Set the values of trans */
    switch (test->gen_trans)
    {
	case NO_TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->trans[batch_iter] = BblasNoTrans;
	    }
	    break;

	case TRANS:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->trans[batch_iter] = BblasTrans;
	    }
	    break;

	case CONJ:

	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		test->trans[batch_iter] = BblasConjTrans;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50)
		{
		    test->trans[batch_iter] = BblasNoTrans;

		}else if (random_number < 80)
		{
		    test->trans[batch_iter] = BblasTrans;
		}else
		{
		    test->trans[batch_iter] = BblasConjTrans;
		}
	    }
	    break;
    }
}

/**
 * Allocate memory and set the values of side
 **/

void bblas_zsetside(bblas_ztest_t *test)
{

    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "side";
    int random_number;

    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation for side */
    test->side =  (enum BBLAS_SIDE*) malloc(nb_data*sizeof(enum BBLAS_SIDE));

    /*checking memory allocation */
    bblas_malloc_check(test->side, ptr_name);

    /*Set the values of side */
    switch (test->gen_side)
    {
	case SIDE_LEFT:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->side[batch_iter] = BblasLeft;
	    }
	    break;

	case SIDE_RIGHT:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->side[batch_iter] = BblasRight;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50 )
		{
		    test->side[batch_iter] = BblasLeft;

		}else
		{
		    test->side[batch_iter] = BblasRight;
		}
		break;
	    }
    }
}

/**
 * Allocate memory and set the values of diag
 **/

void bblas_zsetdiag(bblas_ztest_t *test)
{

    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "diag";
    int random_number;

    /*Initialize random number generation */
    srand ( time(NULL) ) ;

    /*Memory allocation for diag */
    test->diag =  (enum BBLAS_DIAG*) malloc(nb_data*sizeof(enum BBLAS_DIAG));

    /*checking memory allocation */
    bblas_malloc_check(test->diag, ptr_name);

    /*Set the values of diag */
    switch (test->gen_diag)
    {
	case DIAG_NO_U:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->diag[batch_iter] = BblasNonUnit;
	    }
	    break;

	case DIAG_U:

	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		test->diag[batch_iter] = BblasUnit;
	    }
	    break;

	default:
	    for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	    {
		/*Generate a random number */
		random_number = rand() % 100;

		if (random_number < 50 )
		{
		    test->diag[batch_iter] = BblasNonUnit;

		}else
		{
		    test->diag[batch_iter] = BblasUnit;
		}
		break;
	    }
    }
}

/**
 * Allocate memory and set the values of M
 **/

void bblas_zsetM(bblas_ztest_t *test)
{

    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "M";

    /*Memory allocation for M */
    test->M = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->M, ptr_name);

    /*Set the values of M */
    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
    {
	test->M[batch_iter] = irandRange(test->minM, test->maxM);
    }
}

/**
 * Allocate memory and set the values of N
 **/

void bblas_zsetN(bblas_ztest_t *test)
{

    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "N";

    /*Memory allocation for N */
    test->N = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->N, ptr_name);

    /*Set the values of N */
    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
    {
	test->N[batch_iter] = irandRange(test->minN, test->maxN);
    }
}

/**
 * Allocate memory and set the values of K
 **/

void bblas_zsetK(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "K";

    /*Memory allocation for K */
    test->K = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->K, ptr_name);

    /*Set the values of K */
    for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
    {
	test->K[batch_iter] = irandRange(test->minK, test->maxK);
    }
}

/**
 * Allocate memory and set the values of lda
 **/

void bblas_zsetlda(bblas_ztest_t *test)
{
    int batch_iter;
    int nb_data                 = bblas_znbdata(test);
    int routine                 = test->routine;
    char ptr_name[NAME_LENGTH]  = "lda";

    /*Memory allocation for lda */
    test->lda = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->lda, ptr_name);

    /*LDA for GEMM */
    if (routine == BBLAS_GEMM)
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    if (test->transA[batch_iter] == BblasNoTrans)
	    {
		test->lda[batch_iter] = test->M[batch_iter];
	    }else
	    {
		test->lda[batch_iter] =test->K[batch_iter] ;
	    }
	}
    }

    /*LDA for SYMM, HEMM, TRMM AND TRSM */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    if (test->side[batch_iter] == BblasLeft)
	    {
		test->lda[batch_iter] = test->M[batch_iter];
	    }else
	    {
		test->lda[batch_iter] =test->N[batch_iter] ;
	    }
	}
    }

    /*LDA for SYRK, HERK, SYR2K, HER2K */
    if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    if (test->trans[batch_iter] == BblasNoTrans)
	    {
		test->lda[batch_iter] = test->N[batch_iter];
	    }else
	    {
		test->lda[batch_iter] =test->K[batch_iter] ;
	    }
	}
    }
}

/**
 * Allocate memory and set the values of ldb
 **/

void bblas_zsetldb(bblas_ztest_t *test)
{
    int batch_iter;
    int routine                 = test->routine;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "ldb";

    /*Memory allocation for ldb */
    test->ldb = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->ldb, ptr_name);

    /*LDB for GEMM */
    if ((routine == BBLAS_GEMM) )
    {
	for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	{
	    if (test->transB[batch_iter] == BblasNoTrans)
	    {
		test->ldb[batch_iter] = test->K[batch_iter];
	    }else
	    {
		test->ldb[batch_iter] =test->N[batch_iter] ;
	    }
	}
    }

    /*LDB SYMM, HEMM, TRMM AND TRSM */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	{
	    test->ldb[batch_iter] = test->M[batch_iter];
	}
    }

    /*LDB for SYR2K, HER2K */
    if ((routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
    {
	for( batch_iter =0; batch_iter < nb_data; batch_iter++)
	{
	    if (test->trans[batch_iter] == BblasNoTrans)
	    {
		test->ldb[batch_iter] = test->N[batch_iter];
	    }else
	    {
		test->ldb[batch_iter] =test->K[batch_iter] ;
	    }
	}
    }
}

/**
 * Allocate memory and set the values of ldc
 **/

void bblas_zsetldc(bblas_ztest_t *test)
{
    int batch_iter;
    int routine                 = test->routine;
    int nb_data                 = bblas_znbdata(test);
    char ptr_name[NAME_LENGTH]  = "ldc";

    /*Memory allocation for ldc */
    test->ldc = (int*) malloc(nb_data*sizeof(int));

    /*checking memory allocation */
    bblas_malloc_check(test->ldc, ptr_name);

    /*LDC for GEMM */
    if (routine == BBLAS_GEMM)
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    test->ldc[batch_iter] = test->M[batch_iter];
	}
    }

    /*LDC for SYMM, HEMM */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM))
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    test->ldc[batch_iter] = test->M[batch_iter];
	}
    }

    /*LDC for SYRK, HERK, SYR2K, HER2K */
    if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
    {
	for( batch_iter =0; batch_iter < nb_data ; batch_iter++)
	{
	    test->ldc[batch_iter] = test->N[batch_iter];
	}
    }
}

/**
 * Allocate memory and set the values of alpha
 **/

void bblas_zsetalpha(bblas_ztest_t *test)
{
    int batch_iter;
    int batch_count             = test->batch_count;
    char ptr_name[NAME_LENGTH]  = "alpha";

    /*Memory allocation for alpha */
    test->alpha = (BBLAS_Complex64_t*) malloc(batch_count*sizeof(BBLAS_Complex64_t));

    /*checking memory allocation */
    bblas_malloc_check(test->alpha, ptr_name);

    /* Set value of alpha */
    for( batch_iter =0; batch_iter < batch_count; batch_iter++)
    {
	test->alpha[batch_iter] = ((BBLAS_Complex64_t)rand()/(BBLAS_Complex64_t)RAND_MAX);
    }
}

/**
 * Allocate memory and set the values of alpha specifically for herk.
 **/

void bblas_zsetalpha_herk(bblas_ztest_t *test)
{

    int batch_iter;
    int batch_count             = test->batch_count;
    char ptr_name[NAME_LENGTH]  = "alpha_herk";

    /*Memory allocation for alpha_herk */
    test->alpha_herk = (double*) malloc(batch_count*sizeof(double));

    /*checking memory allocation */
    bblas_malloc_check(test->alpha_herk, ptr_name);

    /* Set value of alpha_herk */
    for( batch_iter =0; batch_iter < batch_count; batch_iter++)
    {
	test->alpha_herk[batch_iter] = ((double)rand()/(double)RAND_MAX);
    }
}

/**
 * Allocate memory and set the values of beta.
 **/

void bblas_zsetbeta(bblas_ztest_t *test)
{
    int batch_iter;
    int batch_count             = test->batch_count;
    char ptr_name[NAME_LENGTH]  = "beta";

    /*Memory allocation for beta */
    test->beta = (BBLAS_Complex64_t*) malloc(batch_count*sizeof(BBLAS_Complex64_t));

    /*checking memory allocation */
    bblas_malloc_check(test->beta, ptr_name);

    /* Set value of beta */
    for( batch_iter =0; batch_iter < batch_count; batch_iter++)
    {
	test->beta[batch_iter] = ((BBLAS_Complex64_t)rand()/(BBLAS_Complex64_t)RAND_MAX);
    }
}

/**
 * Allocate memory and set the values of beta specifically for herk.
 **/

void bblas_zsetbeta_herk(bblas_ztest_t *test)
{
     int batch_iter;
     int batch_count             = test->batch_count;
     char ptr_name[NAME_LENGTH]  = "beta_herk";

     /*Memory allocation for beta_herk */
     test->beta_herk = (double*) malloc(batch_count*sizeof(double));

     /*checking memory allocation */
     bblas_malloc_check(test->beta_herk, ptr_name);

     /* Set value of beta_herk */
     for( batch_iter =0; batch_iter < batch_count; batch_iter++)
     {
	 test->beta_herk[batch_iter] = ((double)rand()/(double)RAND_MAX);
     }
}

/**
 * Allocate memory and set the values of arrayA.
 **/

void bblas_zsetarrayA(bblas_ztest_t *test)
{
    int batch_iter, nb_row, nb_col;
     int first_index             = 0;
     int routine                 = test->routine;
     int batch_count             = test->batch_count;
     char ptr_name[NAME_LENGTH]  = "arrayA";

     int IONE     = 1;
     int ISEED[4] ={0,0,0,1};

     /*Memory allocation for **arrayA  */
     test->arrayA = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
     bblas_malloc_check(test->arrayA, ptr_name);

     if( test->batch_opts == BBLAS_VARIABLE )
     {
	 for( batch_iter =0; batch_iter < batch_count; batch_iter++)
	 {
	     nb_row = test->lda[batch_iter];

	     /* nb_col  for GEMM */
	     if (routine == BBLAS_GEMM)
	     {
		 if (test->transA[batch_iter] == BblasNoTrans)
		 {
		     nb_col =test->K[batch_iter];
		 }else
		 {
		     nb_col =test->M[batch_iter];
		 }
	     }

	     /*  nb_col  SYMM, HEMM, TRMM AND TRSM */
	     if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
		 (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	     {
		 if(test->side[batch_iter] == BblasLeft )
		 {
		     nb_col = test->M[batch_iter];
		 }else
		 {
		     nb_col = test->N[batch_iter];
		 }
	     }

	     /* nb_col for SYRK, HERK, SYR2K, HER2K */
	     if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
		 (routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	     {
		 if (test->trans[batch_iter] == BblasNoTrans)
		 {
		     nb_col = test->K[batch_iter];
		 }else
		 {
		     nb_col = test->N[batch_iter];
		 }
	     }

	     /*Matrix filling */
	     test->arrayA[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
	     bblas_malloc_check(test->arrayA[batch_iter], ptr_name);

#if defined(BBLAS_WITH_MKL)
	     LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayA[batch_iter]);
#else
		 LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayA[batch_iter]);
#endif

	     if( (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	       {
		 for(int i=0; i<max(nb_row,nb_col); i++)
		   {
		     test->arrayA[batch_iter][nb_col*i+i] = test->arrayA[batch_iter][nb_col*i+i] + 1.0;
		   }
	       }

	     if(routine == BBLAS_HEMM )
	       {
		 for(int i=0; i< nb_row; i++)
		   {
		     test->arrayA[batch_iter][nb_col*i+i] = creal(test->arrayA[batch_iter][nb_col*i+i]);
		   }
	       }
	 }

     }else if( test->batch_opts == BBLAS_FIXED )
     {
	 nb_row = test->lda[first_index];

	 /* nb_col  for GEMM */
	 if (routine == BBLAS_GEMM)
	 {
	     if (test->transA[first_index] == BblasNoTrans)
	     {
		 nb_col =test->K[first_index];
	     }else
	     {
		 nb_col =test->M[first_index];
	     }
	 }

	 /*  nb_col  SYMM, HEMM, TRMM AND TRSM */
	 if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	     (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	 {
	     if(test->side[first_index] == BblasLeft )
	     {
		 nb_col = test->M[first_index];
	     }else
	     {
		 nb_col = test->N[first_index];
	     }
	 }

	 /* nb_col for SYRK, HERK, SYR2K, HER2K */
	 if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	     (routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	 {
	     if (test->trans[first_index] == BblasNoTrans)
	     {
		 nb_col = test->K[first_index];
	     }else
	     {
		 nb_col = test->N[first_index];
	     }
	 }

	 /*Matrix filling */
	 for( batch_iter =0; batch_iter < batch_count; batch_iter++)
	 {
	     test->arrayA[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
	     bblas_malloc_check(test->arrayA[batch_iter], ptr_name);

#if defined(BBLAS_WITH_MKL)
	     LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayA[batch_iter]);
#else
		 LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayA[batch_iter]);
#endif

	     if( (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	       {
		 for(int i=0; i<max(nb_row,nb_col); i++)
		   {
		     test->arrayA[batch_iter][nb_col*i+i] = test->arrayA[batch_iter][nb_col*i+i] + 1.0;
		   }
	       }

	     if(routine == BBLAS_HEMM )
	       {
		 for(int i=0; i<max(nb_row,nb_col); i++)
		   {
		     test->arrayA[batch_iter][nb_col*i+i] = creal(test->arrayA[batch_iter][nb_col*i+i]);
		   }
	       }
	 }
     }else
     {
	 bblas_error("bblas_ztesting.c", "wrong batch_opts value");
     }
}

/**
 * Allocate memory and set the values of arrayB.
 **/

void bblas_zsetarrayB(bblas_ztest_t *test)
{
    int batch_iter, nb_row, nb_col;
    int first_index             = 0;
    int routine                 = test->routine;
    int batch_count             = test->batch_count;
    char ptr_name[NAME_LENGTH]  = "arrayB";
    int max_work_size = max(test->maxK, max(test->maxM, test->maxN));
    int IONE     = 1;
    int ISEED[4] ={0,0,0,1};

    /*Memory allocation for arrayB */
    test->arrayB = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
    bblas_malloc_check(test->arrayB, ptr_name);

    test->Binitnorm = (double *)malloc(batch_count*sizeof(double));
    bblas_malloc_check(test->Binitnorm, "Binitnorm");

    double *work = (double *)malloc(max_work_size*sizeof(double));
    bblas_malloc_check(work, "work");

    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    nb_row =test->ldb[batch_iter];

	    /* nb_col  for GEMM */
	    if (routine == BBLAS_GEMM)
	    {
		if (test->transB[batch_iter] == BblasNoTrans)
		{
		    nb_col = test->N[batch_iter];
		}else
		{
		    nb_col = test->K[batch_iter];
		}
	    }

	    /* nb_col for SYMM, HEMM, TRMM AND TRSM */
	    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
		(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	    {
		nb_col =test->N[batch_iter];
	    }

	    /* nb_col for SYR2K, HER2K */
	    if ((routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	    {
		if (test->trans[batch_iter] == BblasNoTrans)
		{
		    nb_col = test->K[batch_iter];
		}else
		{
		    nb_col = test->N[batch_iter];
		}
	    }

	    test->arrayB[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
	    bblas_malloc_check(test->arrayB[batch_iter], ptr_name);

#if defined(BBLAS_WITH_MKL)
	    LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayB[batch_iter]);
#else
		LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayB[batch_iter]);
#endif

	    /*Compute the infinity norm of B  */
	    if( (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	    {
#if defined(BBLAS_WITH_MKL)
		test->Binitnorm[batch_iter] = (BBLAS_Complex64_t) LAPACKE_zlange_work(LAPACK_COL_MAJOR,
								  'I', nb_row, nb_col,
								  (MKL_Complex16*) test->arrayB[batch_iter],
								  test->ldb[batch_iter], work);
#else
		test->Binitnorm[batch_iter] = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
								  'I', nb_row, nb_col,
								  test->arrayB[batch_iter],
								  test->ldb[batch_iter], work);
#endif
	    }
	}

    }else if( test->batch_opts == BBLAS_FIXED )
    {
	nb_row =test->ldb[first_index];

	/* nb_col  for GEMM */
	if (routine == BBLAS_GEMM)
	{
	    if (test->transB[first_index] == BblasNoTrans)
	    {
		nb_col = test->N[first_index];
	    }else
	    {
		nb_col = test->K[first_index];
	    }
	}

	/* nb_col for SYMM, HEMM, TRMM AND TRSM */
	if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	    (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	{
	    nb_col =test->N[first_index];
	}

	/* nb_col for SYR2K, HER2K */
	if ((routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
	{
	    if (test->trans[first_index] == BblasNoTrans)
	    {
		nb_col = test->K[first_index];
	    }else
	    {
		nb_col = test->N[first_index];
	    }
	}
	/*Matrix filling */
	for( batch_iter =0; batch_iter < batch_count; batch_iter++)
	{
	     test->arrayB[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
#if defined(BBLAS_WITH_MKL)
	     LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayB[batch_iter]);
#else
		 LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayB[batch_iter]);
#endif

	     /*Compute the infinity norm of B  */
	     if( (routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
	     {
#if defined(BBLAS_WITH_MKL)
		 test->Binitnorm[batch_iter] = (BBLAS_Complex64_t) LAPACKE_zlange_work(LAPACK_COL_MAJOR,
								   'I', nb_row, nb_col,
								   (MKL_Complex16*) test->arrayB[batch_iter],
								   test->ldb[first_index], work);
#else
		 test->Binitnorm[batch_iter] = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
								   'I', nb_row, nb_col,
								   test->arrayB[batch_iter],
								   test->ldb[first_index], work);
#endif
	     }
	 }
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }
    /*Free work */
    free(work);
}

/**
 * Allocate memory and set the values of arrayC.
 **/

void bblas_zsetarrayC(bblas_ztest_t *test)
{
    int batch_iter, nb_row, nb_col;
    int first_index             = 0;
    int batch_count             = test->batch_count;
    int max_work_size = max(test->maxK, max(test->maxM, test->maxN));
    char ptr_name[NAME_LENGTH]  = "arrayC";

    int IONE     = 1;
    int ISEED[4] ={0,0,0,1};

    test->arrayC = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
    bblas_malloc_check(test->arrayC, ptr_name);

    test->Cinitnorm = (double *)malloc(batch_count*sizeof(double));
    bblas_malloc_check(test->Cinitnorm, "Cinitnorm");

    double *work = (double *)malloc(max_work_size*sizeof(double));
    bblas_malloc_check(work, "work");

    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    nb_row =test->ldc[batch_iter];
	    nb_col =test->N[batch_iter];

	    test->arrayC[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
	    bblas_malloc_check(test->arrayC[batch_iter], ptr_name);
#if defined(BBLAS_WITH_MKL)
	    LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayC[batch_iter]);
#else
		LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayC[batch_iter]);
#endif

	    if( (test->routine == BBLAS_HERK) || (test->routine == BBLAS_HER2K))
	      {
		for(int i=0; i<max(nb_row,nb_col); i++)
		  {
		    test->arrayC[batch_iter][nb_col*i+i] = creal(test->arrayC[batch_iter][nb_col*i+i]);
		  }
	      }

	    /*Compuptation of the norm of C */
#if defined(BBLAS_WITH_MKL)
	    test->Cinitnorm[batch_iter] = (BBLAS_Complex64_t) LAPACKE_zlange_work(LAPACK_COL_MAJOR,
								  'I', nb_row, nb_col, (MKL_Complex16*) test->arrayC[batch_iter],
								  test->ldc[batch_iter], work);
#else
		test->Cinitnorm[batch_iter] = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
							      'I', nb_row, nb_col, test->arrayC[batch_iter],
								  test->ldc[batch_iter], work);
#endif
	}

    }else if( test->batch_opts == BBLAS_FIXED )
    {
	nb_row =test->ldc[first_index];
	nb_col =test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    test->arrayC[batch_iter] = (BBLAS_Complex64_t *) malloc(nb_row*nb_col* sizeof(BBLAS_Complex64_t ));
	    bblas_malloc_check(test->arrayC[batch_iter], ptr_name);

#if defined(BBLAS_WITH_MKL)
	    LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, (MKL_Complex16*) test->arrayC[batch_iter]);
#else
		LAPACKE_zlarnv_work(IONE, ISEED, nb_row*nb_col, test->arrayC[batch_iter]);
#endif

	    if( (test->routine == BBLAS_HERK) || (test->routine == BBLAS_HER2K))
	      {
		for(int i=0; i<max(nb_row,nb_col); i++)
		  {
		    test->arrayC[batch_iter][nb_col*i+i] = creal(test->arrayC[batch_iter][nb_col*i+i]);
		  }
	      }

	    /*Compuptation of the norm of C */
#if defined(BBLAS_WITH_MKL)
	    test->Cinitnorm[batch_iter] = (BBLAS_Complex64_t) LAPACKE_zlange_work(LAPACK_COL_MAJOR,
									  'I', nb_row, nb_col,
									  (MKL_Complex16*) test->arrayC[batch_iter],
									  test->ldc[first_index], work);
		#else
		test->Cinitnorm[batch_iter] = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
							      'I',nb_row, nb_col, test->arrayC[batch_iter],
								  test->ldc[first_index], work);
		#endif
	}
    } else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

    /*Free work */
    free(work);
}

/**
 * Allocate memory for error checking.
 **/

void bblas_zmalloc_result_error(bblas_ztest_t *test)
{

    int batch_count             = test->batch_count;

    /*Memory for error computation  */
    switch(test->target)
    {
	case BBLAS_MKL:
	    test->mkl_result  = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
	    bblas_malloc_check(test->mkl_result, "mkl_result");

	    test->mkl_error = (double*) malloc(batch_count*sizeof(double));
	    bblas_malloc_check(test->mkl_error, "mkl_error");
	    break;

	case BBLAS_CUBLAS:
	case BBLAS_MAGMA:
	    test->device_result = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
	    bblas_malloc_check(test->device_result, "device_result");

	    test->device_error = (double*) malloc(batch_count*sizeof(double));
	    bblas_malloc_check(test->device_error, "device_error");
	    break;

	case BBLAS_OTHER:
	    test->other_result  = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
	    bblas_malloc_check(test->other_result, "other_result");

	    test->other_error = (double*) malloc(batch_count*sizeof(double));
	    bblas_malloc_check(test->other_error, "other_error");
	    break;

	case BBLAS_CUMKL:
	    test->mkl_result    = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
	    bblas_malloc_check(test->mkl_result, "mkl_result");

	    test->device_result = (BBLAS_Complex64_t**) malloc(batch_count*sizeof(BBLAS_Complex64_t*));
	    bblas_malloc_check(test->device_result, "device_result");

	    test->mkl_error = (double*) malloc(batch_count*sizeof(double));
	    bblas_malloc_check(test->mkl_error, "mkl_error");

	    test->device_error = (double*) malloc(batch_count*sizeof(double));
	    bblas_malloc_check(test->device_error, "device_error");
	    break;

	default:
	    printf("Memory alloation for error: Target no defined\n");
	    exit(EXIT_FAILURE);
    }

    /*Memory for info  */
    test->info = (int*) malloc(batch_count*sizeof(int));
}

/**
 * Allocate memory and copy the initial arrayC values.
 * This is required to take norms after the computation is complete.
 **/

void bblas_zcopy_Cinit(bblas_ztest_t *test, BBLAS_Complex64_t **C_copy)
{
    /*Local variables */
    enum BBLAS_ROUTINE routine  = test->routine;
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int ldc, N;

    if (!((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	  (routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	  (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	  (routine == BBLAS_HER2K)))
    {
	printf("BBLAS FATAL ERROR: bblas_ztesting.c():\n");
	printf("\t bblas_zcopy_Cinit not defined for %s \n", bblas_getroutine(test->routine));
	exit(EXIT_FAILURE);
    }

    if( test->batch_opts == BBLAS_VARIABLE ) // Varible size
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    ldc = test->ldc[batch_iter];
	    N   = test->N[batch_iter];
	    C_copy[batch_iter] = (BBLAS_Complex64_t *) malloc(ldc*N* sizeof(BBLAS_Complex64_t ));

	    /*Copy the matrix {C}_i */
	    cblas_zcopy (ldc*N, test->arrayC[batch_iter], 1, C_copy[batch_iter], 1);
	}

    }else if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	ldc = test->ldc[first_index];
	N = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    C_copy[batch_iter] = (BBLAS_Complex64_t *) malloc(ldc*N* sizeof(BBLAS_Complex64_t ));

	    /*Copy the matrix {C}_i */
	    cblas_zcopy (ldc*N, test->arrayC[batch_iter], 1,  C_copy[batch_iter], 1);
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value\n");
    }
}

/**
 * Allocate memory and copy the initial arrayB values.
 * This is required to take norms after the computation is complete.
 **/

void bblas_zcopy_Binit(bblas_ztest_t *test, BBLAS_Complex64_t **B_copy)
{
    /*Local variables */
    enum BBLAS_ROUTINE routine  = test->routine;
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int ldb, N;

    if (!((routine == BBLAS_TRMM) || (routine == BBLAS_TRSM)))
    {
	printf("BBLAS FATAL ERROR: bblas_ztesting.c():\n");
	printf("\t bblas_zcopy_Binit not defined for %s \n", bblas_getroutine(test->routine));
	exit(EXIT_FAILURE);
    }

    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    ldb = test->ldb[batch_iter];
	    N   = test->N[batch_iter];

	    B_copy[batch_iter] = (BBLAS_Complex64_t *) malloc(ldb*N* sizeof(BBLAS_Complex64_t ));

	    /*Copy the matrix {B}_i */
	    cblas_zcopy (ldb*N, test->arrayB[batch_iter], 1, B_copy[batch_iter], 1);
	}

    }else if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	ldb = test->ldb[first_index];
	N = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    B_copy[batch_iter] = (BBLAS_Complex64_t *) malloc(ldb*N* sizeof(BBLAS_Complex64_t ));

	    /*Copy the matrix {B}_i */
	    cblas_zcopy (ldb*N, test->arrayB[batch_iter], 1, B_copy[batch_iter], 1);
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value\n");
    }
}

/**
 * Compute the relative error of each batch operation.
 **/

void bblas_zcheck_Cfinal(bblas_ztest_t *test, BBLAS_Complex64_t **C_final)
{

    /*Local variables */
    enum BBLAS_ROUTINE routine  = test->routine;
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int ldc, N;

    /* Error calculation variables */
    double Cnorm, Error_norm;
    BBLAS_Complex64_t alpha =-1;

    /*Check if the call has been made by the correct routine */
    if (!((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	  (routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	  (routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	  (routine == BBLAS_HER2K)))
    {
	printf("BBLAS FATAL ERROR: bblas_ztesting.c():\n");
	printf("\t bblas_zcheck_Cfinal not defined for %s \n", bblas_getroutine(test->routine));
	exit(EXIT_FAILURE);
    }

    /*Temporary buffer to save (test-arrayC -C_final) */
    BBLAS_Complex64_t **C_diff;
    C_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayC in C_diff */
    bblas_zcopy_Cinit(test, C_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    ldc = test->ldc[batch_iter];
	    N   = test->N[batch_iter];

	    /*Computation of the Frobenus norm of  {C}_batch_iter */
	    Cnorm = cblas_dznrm2(ldc*N, test->arrayC[batch_iter], 1);

	    /*Compute the error */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the norm assoicated  with the error */
	    Error_norm = cblas_dznrm2(ldc*N, C_diff[batch_iter], 1);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = Error_norm/Cnorm;
		    break;

		case BBLAS_CUBLAS:
                case BBLAS_MAGMA:
  		    test->device_error[batch_iter] = Error_norm/Cnorm;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = Error_norm/Cnorm;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal, Variable: Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED ) //fixed  size
    {
	ldc = test->ldc[first_index];
	N   = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Computation of the Frobenus norm of  {C}_batch_iter */
	    Cnorm = cblas_dznrm2(ldc*N, test->arrayC[batch_iter], 1);

	    /*Compute the error */
	    cblas_zaxpy (ldc*N, CBLAS_SADDR(alpha), C_final[batch_iter], 1, C_diff[batch_iter], 1);

	    /*Compute  the norm assoicated  with the error */
	    Error_norm = cblas_dznrm2(ldc*N, C_diff[batch_iter], 1);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = Error_norm/Cnorm;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = Error_norm/Cnorm;
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = Error_norm/Cnorm;
		    break;
		default:
		    printf("In bblas_zcheck_Cfinal: Fixed, Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

    /*Free  C_diff */
    for(batch_iter=0; batch_iter < batch_count ; batch_iter++)
    {
	free(C_diff[batch_iter]);
    }
    free(C_diff);
}

/**
 * Compute the relative error of each batch operation.
 **/

void bblas_zcheck_Bfinal(bblas_ztest_t *test, BBLAS_Complex64_t **B_final)
{

    /*Local variables */
    enum BBLAS_ROUTINE routine  = test->routine;
    int batch_count = test->batch_count;
    int batch_iter, first_index = 0;
    int ldb, N;

    /* Error calculation variables */
    double Bnorm, Error_norm;
    BBLAS_Complex64_t alpha =-1;

    /*Check if the call has been made by the correct routine */
    if (!((routine == BBLAS_TRSM) || (routine == BBLAS_TRMM)))
    {
	printf("BBLAS FATAL ERROR: bblas_ztesting.c():\n");
	printf("\t bblas_zcheck_Bfinal not defined for %s \n", bblas_getroutine(test->routine));
	exit(EXIT_FAILURE);
    }

    /*Temporary buffer to save (test-arrayB -B_final) */
    BBLAS_Complex64_t **B_diff;
    B_diff = (BBLAS_Complex64_t **) malloc(batch_count*sizeof(BBLAS_Complex64_t *));

    /*Make a copy of test->arrayB in B_diff */
    bblas_zcopy_Binit(test, B_diff);

    /*Variable size */
    if( test->batch_opts == BBLAS_VARIABLE )
    {
	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    ldb = test->ldb[batch_iter];
	    N   = test->N[batch_iter];

	    /*Computation of the Frobenus norm of  {B}_batch_iter */
	    Bnorm = cblas_dznrm2(ldb*N, test->arrayB[batch_iter], 1);

	    /*Compute the error */
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

	    /*Compute  the norm assoicated  with the error */
	    Error_norm = cblas_dznrm2(ldb*N, B_diff[batch_iter], 1);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = Error_norm/Bnorm;
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = Error_norm/(Bnorm);
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = Error_norm/(Bnorm);
		    break;
		default:
		    printf("In bblas_zcheck_Bfinal(): Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}

    }else  if( test->batch_opts == BBLAS_FIXED )
    {
	ldb = test->ldb[first_index];
	N = test->N[first_index];

	for( batch_iter =0; batch_iter < batch_count ; batch_iter++)
	{
	    /*Computation of the Frobenus norm of  {B}_batch_iter */
	    Bnorm = cblas_dznrm2(ldb*N, test->arrayB[batch_iter], 1);

	    /*Compute the error */
	    cblas_zaxpy (ldb*N, CBLAS_SADDR(alpha), B_final[batch_iter], 1, B_diff[batch_iter], 1);

	    /*Compute  the norm assoicated  with the error */
	    Error_norm = cblas_dznrm2(ldb*N, B_diff[batch_iter], 1);

	    /*Compute the relative error */
	    switch(test->target)
	    {
		case BBLAS_MKL:
		    test->mkl_error[batch_iter] = Error_norm/(Bnorm);
		    break;

		case BBLAS_CUBLAS:
		case BBLAS_MAGMA:
		    test->device_error[batch_iter] = Error_norm/(Bnorm);
		    break;

		case BBLAS_OTHER:
		    test->other_error[batch_iter] = Error_norm/(Bnorm);
		    break;
		default:
		    printf("In bblas_zcheck_Bfinal(): Target no defined\n");
		    exit(EXIT_FAILURE);
	    }
	}
    }else
    {
	bblas_error("bblas_ztesting.c", "wrong batch_opts value");
    }

    /*Free  B_diff */
    for(batch_iter=0; batch_iter < batch_count ; batch_iter++)
    {
	free(B_diff[batch_iter]);
    }
    free(B_diff);
}

/**
 * Free the memory associated with the test structure.
 **/

void bblas_zfreetest( bblas_ztest_t *test )
{
    enum BBLAS_ROUTINE routine  = test->routine;
    /*Free uplo */
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K)||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	free(test->uplo);
    }

    /*Free TRANSA */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_TRMM) ||
	(routine == BBLAS_TRSM))
    {
	free( test->transA );
    }

    /*Free TRANSB */
    if (routine == BBLAS_GEMM)
    {
	free( test->transB );
    }

    /*Free TRANS */
    if ((routine == BBLAS_SYRK) || (routine == BBLAS_HERK) ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K))
    {
	free( test->trans );
    }

    /*Free SIDE*/
    if ((routine == BBLAS_SYMM) || (routine == BBLAS_HEMM) ||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	free( test->side );
    }

    /*Free DIAG*/
    if  ((routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	free(test->diag);
    }

    /* Free M  memory*/
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_TRMM) ||
	(routine == BBLAS_TRSM))
    {
	free( test->M );
    }

    /* Free N memory*/
    free( test->N );

    /* Free K memory*/
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYRK)  ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K))
    {
	free( test->K );
    }

    /*Free LDA memory */
    free( test->lda );

    /*Free LDB memory */
    if ((routine == BBLAS_GEMM)  || (routine == BBLAS_SYMM)  ||
	(routine == BBLAS_HEMM)  || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K) || (routine == BBLAS_TRMM)  ||
	(routine == BBLAS_TRSM))
    {
	free( test->ldb );
    }

    /*Free LDC memory */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	(routine == BBLAS_HER2K))
    {
	free( test->ldc );
    }

    /*Free ALPHA memory */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) ||(routine == BBLAS_SYRK)  ||
	(routine == BBLAS_SYR2K)|| (routine == BBLAS_HER2K)||
	(routine == BBLAS_TRMM) || (routine == BBLAS_TRSM))
    {
	free(test->alpha);

    }else if (routine == BBLAS_HERK)
    {
	free(test->alpha_herk);
    }

    /*Free ALPHA BETA memory */
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_SYR2K))
    {
	free( test->beta  );

    }else if ((routine == BBLAS_HERK) || (routine == BBLAS_HER2K))
    {
	free(test->beta_herk);
    }

    /*Free matrices {A}_i  */
    for ( int batch_iter = 0; batch_iter < test->batch_count; batch_iter++ )
    {
	free( test->arrayA[batch_iter] );
    }
    free( test->arrayA );

    /*Free matrices  {B}_i*/
    if ((routine == BBLAS_GEMM)  || (routine == BBLAS_SYMM)  ||
	(routine == BBLAS_HEMM)  || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K) || (routine == BBLAS_TRMM)  ||
	(routine == BBLAS_TRSM))
    {
	for ( int batch_iter = 0; batch_iter < test->batch_count; batch_iter++ )
	{
	    free( test->arrayB[batch_iter] );
	}
	free( test->arrayB );
    }

    /*Free matrices  {C}_i*/
    if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	(routine == BBLAS_HER2K))
    {
	for ( int batch_iter = 0; batch_iter < test->batch_count; batch_iter++ )
	{
	    free( test->arrayC[batch_iter] );
	}
	free( test->arrayC );
    }

    /*Free memory allocated for error computation */
    switch(test->target)
    {
	case BBLAS_MKL:
	    free(test->mkl_error);
	    //free(test->group_size);
	    for ( int batch_iter = 0; batch_iter < test->batch_count; batch_iter++ )
	    {
		free( test->mkl_result[batch_iter] );
	    }
	    free( test->mkl_result );

	    break;

	case BBLAS_CUBLAS:
	case BBLAS_MAGMA:
	    free(test->device_error);

	    for (int batch_iter =0; batch_iter < test->batch_count; batch_iter++)
	    {
		free(test->device_result[batch_iter]);
	    }

	    free(test->device_result);

	    break;

	case BBLAS_OTHER:

	    free(test->other_error);

	    for (int batch_iter =0; batch_iter < test->batch_count; batch_iter++)
	    {
		free(test->other_result[batch_iter]);
	    }
	    free(test->other_result);

	    break;

	case BBLAS_CUMKL:
	    free(test->mkl_error);
	    free(test->device_error);
	    free(test->group_size);

	    for (int batch_iter =0; batch_iter < test->batch_count; batch_iter++)
	    {
		free(test->device_result[batch_iter]);
		free(test->mkl_result[batch_iter]);
	    }
	    free(test->device_result);
	    free(test->mkl_result);

	    break;

	default:
	    printf("In bblas_zfreetest(): Target no defined\n");
	    exit(EXIT_FAILURE);
    }
    /*Free INFO */
    free( test->info );

    /*Free cuda memory */
    if( (test->target == BBLAS_CUBLAS) || (test->target == BBLAS_MAGMA))
    {
	bblas_zcudaFree(test);
    }

	if ((routine == BBLAS_GEMM)  || (routine == BBLAS_SYMM)  ||
	(routine == BBLAS_HEMM)  || (routine == BBLAS_SYR2K) ||
	(routine == BBLAS_HER2K) || (routine == BBLAS_TRMM)  ||
	(routine == BBLAS_TRSM))
	{
		free(test->Binitnorm);
    }

    /*Free memory allocated for norm computing */
	if ((routine == BBLAS_GEMM) || (routine == BBLAS_SYMM) ||
	(routine == BBLAS_HEMM) || (routine == BBLAS_SYRK) ||
	(routine == BBLAS_HERK) || (routine == BBLAS_SYR2K)||
	(routine == BBLAS_HER2K))
    {
		free(test->Cinitnorm);
    }
}

#undef COMPLEX

ext: c

max_stars_repo_name: NLAFET/BBLAS-ref

max_issues_repo_name: NLAFET/BBLAS-ref

max_forks_repo_name: NLAFET/BBLAS-ref


text:
//===--- Sudoku/Solver_set_option.h                                     ---===//
//
// Helper functions
//===----------------------------------------------------------------------===//
#pragma once

#include "Board.h"
#include "Location.h"
#include "Location_Utilities.h" // is_same_*
#include "Options.h"
#include "Size.h"
#include "Solvers_find.h"
#include "Solvers_remove_option.h"
#include "Value.h"
#include "exceptions.h"
#include "traits.h"

#include <gsl/gsl>

#include <vector>

#include <algorithm>   // find_if
#include <iterator>    // next
#include <stdexcept>   // logic_error
#include <type_traits> // is_base_of

#include "Solver.fwd.h" // Forward declaration - single_option

#include <cassert>

namespace Sudoku
{
//===----------------------------------------------------------------------===//
template<int N, typename Options = Options<elem_size<N>>>
int set_Value(Board<Options, N>&, Location<N>, Value);
template<int N, typename Options = Options<elem_size<N>>, typename ItrT>
int set_Value(Board<Options, N>&, ItrT begin, ItrT end);

template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int set_unique(Board<Options, N>&, SectionT, Value);
template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int set_uniques(Board<Options, N>&, SectionT, Options worker);

template<int N, typename Options = Options<elem_size<N>>, typename SectionT>
int set_section_locals(
	Board<Options, N>&, SectionT, int rep_count, Options values);
template<int N, typename Options = Options<elem_size<N>>>
int set_section_locals(
	Board<Options, N>&,
	Board_Section::Block<Options, N>,
	int rep_count,
	Options values);

//===----------------------------------------------------------------------===//

// IF a possible option, Make [value] the answer for [loc]
template<int N, typename Options>
inline int set_Value(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Location<N> loc,
	const Value value)
{
	assert(is_valid(loc));
	assert(is_valid<N>(value));

	int changes{0};
	auto& elem = board.at(loc);
	if (not elem.test(value)) // value is option nor answer
		throw error::invalid_Board();

	if (not is_answer(elem))
	{
		changes = gsl::narrow_cast<int>(elem.count_all());
		elem.set_nocheck(value);
	}
	return changes;
}

// set board_ using a transferable container of values
template<int N, typename Options, typename ItrT>
int set_Value(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const ItrT begin,
	const ItrT end)
{
	{
		static_assert(traits::is_forward<ItrT>);
		assert(end - begin == full_size<N>);
	}
	int changes{0};
	int n{0};
	for (auto itr = begin; itr != end; ++itr)
	{
		const Location<N> loc(n++); // start at 0!

		// handle different input types
		Value value{};
		if constexpr (traits::iterator_to<ItrT, Value>)
		{
			value = *itr;
		}
		else
		{
			value = to_Value<N>(*itr);
		}

		if (value != Value{0})
		{
			if (not is_valid<N>(value))
			{
				throw std::domain_error{"Invalid Value"};
			}

			if (is_option(board.at(loc), value))
			{ // update options on board
				changes += single_option(board, loc, value);
			}
			assert(is_answer(board.at(loc), value));
		}
	}
	assert(n == full_size<N>);
	return changes;
}

// Set unique values in section as answer
template<int N, typename Options, typename SectionT>
inline int set_uniques(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const Options worker)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
	}
	int changes{0};

	if (worker.count_all() > 0)
	{
		for (Value val{1}; val < Value{worker.size()}; ++val)
		{
			if (worker[val])
			{
				changes += set_unique(board, section, val);
			}
		}
	}
	return changes;
}

// Set unique value in section as answer
template<int N, typename Options, typename SectionT>
inline int set_unique(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const Value value)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
		static_assert(traits::is_input<typename SectionT::iterator>);
		assert(is_valid<N>(value));
	}
	const auto end       = section.cend();
	const auto condition = [value](const Options& O) {
		return O.test(value);
	};

	const auto itr = std::find_if(section.cbegin(), end, condition);

	if (itr == end)
	{ // option not available in section
		throw error::invalid_Board();
	}
	assert(&(*itr) == &board[itr.location()]); // section must be part of board
	assert(std::find_if(std::next(itr), end, condition) == end); // unique

	return single_option(board, itr.location(), value);
}

// for [row/col] per value: if all in same block, remove from rest block
template<int N, typename Options, typename SectionT>
inline int set_section_locals(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const SectionT section,
	const int rep_count,
	const Options values)
{
	{
		static_assert(Board_Section::traits::is_Section_v<SectionT>);
		static_assert(std::is_same_v<typename SectionT::value_type, Options>);
		assert(rep_count > 1);  // should have been caught by caller
								// use the set_uniques specialization
		assert(rep_count <= N); // won't fit in single block-row/col
		assert(values.count_all() > 0);
	}
	int changes{0};
	std::vector<Location<N>> locations{};

	// start at 1, to skip the answer-bit
	for (Value value{1}; value < Value{values.size()}; ++value)
	{
		if (values[value])
		{
			locations = list_where_option<N>(section, value, rep_count);
			if (locations.size() != gsl::narrow_cast<size_t>(rep_count))
			{
				assert(changes > 0); // changed by earlier value in values
			}
			else if (is_same_block<N>(locations.cbegin(), locations.cend()))
			{ // remove from rest of block
				changes += remove_option_section(
					board, board.block(locations[0]), locations, value);
			}
		}
	}
	return changes;
}

// per value in [block]: if all in same row/col, remove from rest row/col
template<int N, typename Options>
inline int set_section_locals(
	Board<Options, N>& board, // NOLINT(runtime/references)
	const Board_Section::Block<Options, N> block,
	const int rep_count,
	const Options values)
{
	{
		assert(rep_count > 1);  // should have been caught by caller
								// use the set_uniques specialization
		assert(rep_count <= N); // won't fit in single block-row/col
		assert(values.count_all() > 0);
	}
	int changes{0};
	std::vector<Location<N>> locations{};

	// start at 1, to skip the answer-bit
	for (Value value{1}; value < Value{values.size()}; ++value)
	{
		if (values[value])
		{
			locations = list_where_option<N>(block, value, rep_count);
			if (locations.size() != gsl::narrow_cast<size_t>(rep_count))
			{
				assert(changes > 0); // changed by earlier value in values
			}
			else if (is_same_row<N>(locations.cbegin(), locations.cend()))
			{ // remove from rest of row
				changes += remove_option_outside_block(
					board, board.row(locations[0]), locations[0], value);
			}
			else if (is_same_col<N>(locations.cbegin(), locations.cend()))
			{ // remove from rest of col
				changes += remove_option_outside_block(
					board, board.col(locations[0]), locations[0], value);
			}
		}
	}
	return changes;
}

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#include "urand.c"

int
main (void)
{
  size_t i, j;
  int status = 0;
  double total = 0;
  size_t N = 200000;

  gsl_histogram2d *h;

  gsl_ieee_env_setup ();

  h = gsl_histogram2d_calloc_uniform (10, 10, 0.0, 1.0, 0.0, 1.0);

  for (i = 0; i < 10; i++)
    {
      for (j = 0; j < 10; j++)
	{
	  double w = 10.0 * i + j;
	  total += w;
	  gsl_histogram2d_accumulate (h, 0.1 * i, 0.1 * i, w);
	}
    }

  {
    gsl_histogram2d_pdf *p = gsl_histogram2d_pdf_alloc (10,10);

    gsl_histogram2d *hh = gsl_histogram2d_calloc_uniform (20, 20,
							  0.0, 1.0,
							  0.0, 1.0);

    gsl_histogram2d_pdf_init (p, h);

    for (i = 0; i < N; i++)
      {
	double u = urand();
	double v = urand();
	double x, y;
	status = gsl_histogram2d_pdf_sample (p, u, v, &x, &y);
	status = gsl_histogram2d_increment (hh, x, y);
      }

    status = 0;
    for (i = 0; i < 20; i++)
      {
	for (j = 0; j < 20; j++)
	  {
	    double z = 4 * total * gsl_histogram2d_get (hh, i, j) / (double) N;
	    size_t k1, k2;
	    double ya;
	    double x, xmax, y, ymax;

	    gsl_histogram2d_get_xrange (hh, i, &x, &xmax);
	    gsl_histogram2d_get_yrange (hh, j, &y, &ymax);

	    gsl_histogram2d_find (h, x, y, &k1, &k2);
	    ya = gsl_histogram2d_get (h, k1, k2);

	    if (ya == 0)
	      {
		if (z != 0)
		  {
		    status = 1;
		    printf ("(%d,%d): %g vs %g\n", (int)i, (int)j, z, ya);
		  }
	      }
	    else
	      {
		double err = 1 / sqrt (gsl_histogram2d_get (hh, i, j));
		double sigma = fabs ((z - ya) / (ya * err));
		if (sigma > 3)
		  {
		    status = 1;
		    printf ("%g vs %g err=%g sigma=%g\n", z, ya, err, sigma);
		  }
	      }
	  }
      }

    gsl_histogram2d_pdf_free (p) ;
    gsl_histogram2d_free (hh) ;

    gsl_test (status, "gsl_histogram2d_pdf_sample within statistical errors");
  }

  gsl_histogram2d_free (h) ;

  exit (gsl_test_summary ());
}


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

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

#include "recurse.h"

static int LU_decomp_L2 (gsl_matrix_complex * A, gsl_vector_uint * ipiv);
static int LU_decomp_L3 (gsl_matrix_complex * A, gsl_vector_uint * ipiv);
static int singular (const gsl_matrix_complex * LU);
static int apply_pivots(gsl_matrix_complex * A, const gsl_vector_uint * ipiv);

/* Factorise a general N x N complex matrix A into,
 *
 *   P A = L U
 *
 * where P is a permutation matrix, L is unit lower triangular and U
 * is upper triangular.
 *
 * L is stored in the strict lower triangular part of the input
 * matrix. The diagonal elements of L are unity and are not stored.
 *
 * U is stored in the diagonal and upper triangular part of the
 * input matrix.
 *
 * P is stored in the permutation p. Column j of P is column k of the
 * identity matrix, where k = permutation->data[j]
 *
 * signum gives the sign of the permutation, (-1)^n, where n is the
 * number of interchanges in the permutation.
 *
 * See Golub & Van Loan, Matrix Computations, Algorithm 3.4.1 (Gauss
 * Elimination with Partial Pivoting).
 */

int
gsl_linalg_complex_LU_decomp (gsl_matrix_complex * A, gsl_permutation * p, int *signum)
{
  const size_t M = A->size1;

  if (p->size != M)
    {
      GSL_ERROR ("permutation length must match matrix size1", GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t N = A->size2;
      const size_t minMN = GSL_MIN(M, N);
      gsl_vector_uint * ipiv = gsl_vector_uint_alloc(minMN);
      gsl_matrix_complex_view AL = gsl_matrix_complex_submatrix(A, 0, 0, M, minMN);
      size_t i;

      status = LU_decomp_L3 (&AL.matrix, ipiv);

      /* process remaining right matrix */
      if (M < N)
        {
          gsl_matrix_complex_view AR = gsl_matrix_complex_submatrix(A, 0, M, M, N - M);

          /* apply pivots to AR */
          apply_pivots(&AR.matrix, ipiv);

          /* AR = AL^{-1} AR */
          gsl_blas_ztrsm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, GSL_COMPLEX_ONE, &AL.matrix, &AR.matrix);
        }

      /* convert ipiv array to permutation */

      gsl_permutation_init(p);
      *signum = 1;

      for (i = 0; i < minMN; ++i)
        {
          unsigned int pivi = gsl_vector_uint_get(ipiv, i);

          if (p->data[pivi] != p->data[i])
            {
              size_t tmp = p->data[pivi];
              p->data[pivi] = p->data[i];
              p->data[i] = tmp;
              *signum = -(*signum);
            }
        }

      gsl_vector_uint_free(ipiv);

      return status;
    }
}

/*
LU_decomp_L2
  LU decomposition with partial pivoting using Level 2 BLAS

Inputs: A    - on input, matrix to be factored; on output, L and U factors
        ipiv - (output) array containing row swaps

Notes:
1) Based on LAPACK ZGETF2
*/

static int
LU_decomp_L2 (gsl_matrix_complex * A, gsl_vector_uint * ipiv)
{
  const size_t M = A->size1;
  const size_t N = A->size2;
  const size_t minMN = GSL_MIN(M, N);

  if (ipiv->size != minMN)
    {
      GSL_ERROR ("ipiv length must equal MIN(M,N)", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      for (j = 0; j < minMN; ++j)
        {
          /* find maximum in the j-th column */
          gsl_vector_complex_view v = gsl_matrix_complex_subcolumn(A, j, j, M - j);
          size_t j_pivot = j + gsl_blas_izamax(&v.vector);
          gsl_vector_complex_view v1, v2;

          gsl_vector_uint_set(ipiv, j, j_pivot);

          if (j_pivot != j)
            {
              /* swap rows j and j_pivot */
              v1 = gsl_matrix_complex_row(A, j);
              v2 = gsl_matrix_complex_row(A, j_pivot);
              gsl_blas_zswap(&v1.vector, &v2.vector);
            }

          if (j < M - 1)
            {
              gsl_complex Ajj = gsl_matrix_complex_get(A, j, j);
              gsl_complex Ajjinv = gsl_complex_inverse(Ajj);

              if (gsl_complex_abs(Ajj) >= GSL_DBL_MIN)
                {
                  v1 = gsl_matrix_complex_subcolumn(A, j, j + 1, M - j - 1);
                  gsl_blas_zscal(Ajjinv, &v1.vector);
                }
              else
                {
                  for (i = 1; i < M - j; ++i)
                    {
                      gsl_complex * ptr = gsl_matrix_complex_ptr(A, j + i, j);
                      *ptr = gsl_complex_mul(*ptr, Ajjinv);
                    }
                }
            }

          if (j < minMN - 1)
            {
              gsl_matrix_complex_view A22 = gsl_matrix_complex_submatrix(A, j + 1, j + 1, M - j - 1, N - j - 1);
              v1 = gsl_matrix_complex_subcolumn(A, j, j + 1, M - j - 1);
              v2 = gsl_matrix_complex_subrow(A, j, j + 1, N - j - 1);

              gsl_blas_zgeru(GSL_COMPLEX_NEGONE, &v1.vector, &v2.vector, &A22.matrix);
            }
        }

      return GSL_SUCCESS;
    }
}

/*
LU_decomp_L3
  LU decomposition with partial pivoting using Level 3 BLAS

Inputs: A    - on input, matrix to be factored; on output, L and U factors
        ipiv - (output) array containing row swaps

Notes:
1) Based on ReLAPACK DGETRF
*/

static int
LU_decomp_L3 (gsl_matrix_complex * A, gsl_vector_uint * ipiv)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M < N)
    {
      GSL_ERROR ("matrix must have M >= N", GSL_EBADLEN);
    }
  else if (ipiv->size != GSL_MIN(M, N))
    {
      GSL_ERROR ("ipiv length must equal MIN(M,N)", GSL_EBADLEN);
    }
  else if (N <= CROSSOVER_LU)
    {
      /* use Level 2 algorithm */
      return LU_decomp_L2(A, ipiv);
    }
  else
    {
      /*
       * partition matrix:
       *
       *       N1  N2
       * N1  [ A11 A12 ]
       * M2  [ A21 A22 ]
       *
       * and
       *      N1  N2
       * M  [ AL  AR  ]
       */
      int status;
      const size_t N1 = GSL_LINALG_SPLIT_COMPLEX(N);
      const size_t N2 = N - N1;
      const size_t M2 = M - N1;
      gsl_matrix_complex_view A11 = gsl_matrix_complex_submatrix(A, 0, 0, N1, N1);
      gsl_matrix_complex_view A12 = gsl_matrix_complex_submatrix(A, 0, N1, N1, N2);
      gsl_matrix_complex_view A21 = gsl_matrix_complex_submatrix(A, N1, 0, M2, N1);
      gsl_matrix_complex_view A22 = gsl_matrix_complex_submatrix(A, N1, N1, M2, N2);

      gsl_matrix_complex_view AL = gsl_matrix_complex_submatrix(A, 0, 0, M, N1);
      gsl_matrix_complex_view AR = gsl_matrix_complex_submatrix(A, 0, N1, M, N2);

      /*
       * partition ipiv = [ ipiv1 ] N1
       *                  [ ipiv2 ] N2
       */
      gsl_vector_uint_view ipiv1 = gsl_vector_uint_subvector(ipiv, 0, N1);
      gsl_vector_uint_view ipiv2 = gsl_vector_uint_subvector(ipiv, N1, N2);

      size_t i;

      /* recursion on (AL, ipiv1) */
      status = LU_decomp_L3(&AL.matrix, &ipiv1.vector);
      if (status)
        return status;

      /* apply ipiv1 to AR */
      apply_pivots(&AR.matrix, &ipiv1.vector);

      /* A12 = A11^{-1} A12 */
      gsl_blas_ztrsm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, GSL_COMPLEX_ONE, &A11.matrix, &A12.matrix);

      /* A22 = A22 - A21 * A12 */
      gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, GSL_COMPLEX_NEGONE, &A21.matrix, &A12.matrix, GSL_COMPLEX_ONE, &A22.matrix);

      /* recursion on (A22, ipiv2) */
      status = LU_decomp_L3(&A22.matrix, &ipiv2.vector);
      if (status)
        return status;

      /* apply pivots to A21 */
      apply_pivots(&A21.matrix, &ipiv2.vector);

      /* shift pivots */
      for (i = 0; i < N2; ++i)
        {
          unsigned int * ptr = gsl_vector_uint_ptr(&ipiv2.vector, i);
          *ptr += N1;
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_solve (const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size2 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;

      /* copy x <- b */
      gsl_vector_complex_memcpy (x, b);

      /* solve for x */
      status = gsl_linalg_complex_LU_svx (LU, p, x);

      return status;
    }
}

int
gsl_linalg_complex_LU_svx (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_vector_complex * x)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution/rhs size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      /* apply permutation to RHS */
      gsl_permute_vector_complex (p, x);

      /* solve for c using forward-substitution, L c = P b */
      gsl_blas_ztrsv (CblasLower, CblasNoTrans, CblasUnit, LU, x);

      /* perform back-substitution, U x = c */
      gsl_blas_ztrsv (CblasUpper, CblasNoTrans, CblasNonUnit, LU, x);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_complex_LU_refine (const gsl_matrix_complex * A, const gsl_matrix_complex * LU, const gsl_permutation * p, const gsl_vector_complex * b, gsl_vector_complex * x, gsl_vector_complex * work)
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("matrix a must be square", GSL_ENOTSQR);
    }
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (A->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be decomposition of a", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (LU->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LU->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else if (LU->size1 != work->size)
    {
      GSL_ERROR ("matrix size must match workspace size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;

      /* Compute residual = (A * x  - b) */

      gsl_vector_complex_memcpy (work, b);

      {
        gsl_complex one = GSL_COMPLEX_ONE;
        gsl_complex negone = GSL_COMPLEX_NEGONE;
        gsl_blas_zgemv (CblasNoTrans, one, A, x, negone, work);
      }

      /* Find correction, delta = - (A^-1) * residual, and apply it */

      status = gsl_linalg_complex_LU_svx (LU, p, work);

      {
        gsl_complex negone= GSL_COMPLEX_NEGONE;
        gsl_blas_zaxpy (negone, work, x);
      }

      return status;
    }
}

int
gsl_linalg_complex_LU_invert (const gsl_matrix_complex * LU, const gsl_permutation * p, gsl_matrix_complex * inverse)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (inverse->size1 != LU->size1 || inverse->size2 != LU->size2)
    {
      GSL_ERROR ("inverse matrix must match LU matrix dimensions", GSL_EBADLEN);
    }
  else
    {
      gsl_matrix_complex_memcpy(inverse, LU);
      return gsl_linalg_complex_LU_invx(inverse, p);
    }
}

int
gsl_linalg_complex_LU_invx (gsl_matrix_complex * LU, const gsl_permutation * p)
{
  if (LU->size1 != LU->size2)
    {
      GSL_ERROR ("LU matrix must be square", GSL_ENOTSQR);
    }
  else if (LU->size1 != p->size)
    {
      GSL_ERROR ("permutation length must match matrix size", GSL_EBADLEN);
    }
  else if (singular (LU))
    {
      GSL_ERROR ("matrix is singular", GSL_EDOM);
    }
  else
    {
      int status;
      const size_t N = LU->size1;
      size_t i;

      /* compute U^{-1} */
      status = gsl_linalg_complex_tri_invert(CblasUpper, CblasNonUnit, LU);
      if (status)
        return status;

      /* compute L^{-1} */
      status = gsl_linalg_complex_tri_invert(CblasLower, CblasUnit, LU);
      if (status)
        return status;

      /* compute U^{-1} L^{-1} */
      status = gsl_linalg_complex_tri_UL(LU);
      if (status)
        return status;

      /* apply permutation to columns of A^{-1} */
      for (i = 0; i < N; ++i)
        {
          gsl_vector_complex_view v = gsl_matrix_complex_row(LU, i);
          gsl_permute_vector_complex_inverse(p, &v.vector);
        }

      return GSL_SUCCESS;
    }
}

gsl_complex
gsl_linalg_complex_LU_det (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex det = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex zi = gsl_matrix_complex_get (LU, i, i);
      det = gsl_complex_mul (det, zi);
    }

  return det;
}

double
gsl_linalg_complex_LU_lndet (gsl_matrix_complex * LU)
{
  size_t i, n = LU->size1;

  double lndet = 0.0;

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);
      lndet += log (gsl_complex_abs (z));
    }

  return lndet;
}

gsl_complex
gsl_linalg_complex_LU_sgndet (gsl_matrix_complex * LU, int signum)
{
  size_t i, n = LU->size1;

  gsl_complex phase = gsl_complex_rect((double) signum, 0.0);

  for (i = 0; i < n; i++)
    {
      gsl_complex z = gsl_matrix_complex_get (LU, i, i);

      double r = gsl_complex_abs(z);

      if (r == 0)
        {
          phase = gsl_complex_rect(0.0, 0.0);
          break;
        }
      else
        {
          z = gsl_complex_div_real(z, r);
          phase = gsl_complex_mul(phase, z);
        }
    }

  return phase;
}

static int
singular (const gsl_matrix_complex * LU)
{
  size_t i, n = LU->size1;

  for (i = 0; i < n; i++)
    {
      gsl_complex u = gsl_matrix_complex_get (LU, i, i);
      if (GSL_REAL(u) == 0 && GSL_IMAG(u) == 0) return 1;
    }

 return 0;
}

static int
apply_pivots(gsl_matrix_complex * A, const gsl_vector_uint * ipiv)
{
  if (0)
    {
    }
  else
    {
      size_t i;

      for (i = 0; i < ipiv->size; ++i)
        {
          size_t pi = gsl_vector_uint_get(ipiv, i);

          if (i != pi)
            {
              /* swap rows i and pi */
              gsl_vector_complex_view v1 = gsl_matrix_complex_row(A, i);
              gsl_vector_complex_view v2 = gsl_matrix_complex_row(A, pi);
              gsl_blas_zswap(&v1.vector, &v2.vector);
            }
        }

      return GSL_SUCCESS;
    }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_LUc_decomp_eps(const gsl_matrix_complex * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;

  gsl_matrix_complex * A = gsl_matrix_complex_alloc(M, N);
  gsl_matrix_complex * PLU = gsl_matrix_complex_calloc(M, N);
  gsl_permutation * p = gsl_permutation_alloc(M);
  size_t i, j;
  int signum;

  gsl_matrix_complex_memcpy(A, m);
  gsl_linalg_complex_LU_decomp(A, p, &signum);

  if (M >= N)
    {
      gsl_matrix_complex_view U = gsl_matrix_complex_submatrix(A, 0, 0, N, N);

      /* copy L factor into PLU */
      for (j = 0; j < GSL_MIN(M, N); ++j)
        {
          gsl_vector_complex_view v1 = gsl_matrix_complex_subcolumn(A, j, j, M - j);
          gsl_vector_complex_view v2 = gsl_matrix_complex_subcolumn(PLU, j, j, M - j);

          gsl_vector_complex_memcpy(&v2.vector, &v1.vector);
          gsl_matrix_complex_set(PLU, j, j, GSL_COMPLEX_ONE);
        }

      /* PLU := L * U */
      gsl_blas_ztrmm(CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, GSL_COMPLEX_ONE, &U.matrix, PLU);
    }
  else
    {
      gsl_matrix_complex_view L = gsl_matrix_complex_submatrix(A, 0, 0, M, M);

      /* copy U factor into PLU */
      for (i = 0; i < GSL_MIN(M, N); ++i)
        {
          gsl_vector_complex_view v1 = gsl_matrix_complex_subrow(A, i, i, N - i);
          gsl_vector_complex_view v2 = gsl_matrix_complex_subrow(PLU, i, i, N - i);
          gsl_vector_complex_memcpy(&v2.vector, &v1.vector);
        }

      /* PLU := L * U */
      gsl_blas_ztrmm(CblasLeft, CblasLower, CblasNoTrans, CblasUnit, GSL_COMPLEX_ONE, &L.matrix, PLU);
    }

  /* PLU := P * L * U */
  for (j = 0; j < N; ++j)
    {
      gsl_vector_complex_view v = gsl_matrix_complex_column(PLU, j);
      gsl_permute_vector_complex_inverse(p, &v.vector);
    }

  /* now test m = PLU */
  for (i = 0; i < M; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex aij = gsl_matrix_complex_get(PLU, i, j);
          gsl_complex mij = gsl_matrix_complex_get(m, i, j);

          gsl_test_rel(GSL_REAL(aij), GSL_REAL(mij), eps,
                       "%s real: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, j, GSL_REAL(aij), GSL_REAL(mij));

          gsl_test_rel(GSL_IMAG(aij), GSL_IMAG(mij), eps,
                       "%s imag: (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i, j, GSL_IMAG(aij), GSL_IMAG(mij));
        }
    }

  gsl_matrix_complex_free(A);
  gsl_matrix_complex_free(PLU);
  gsl_permutation_free(p);

  return s;
}

static int
test_LUc_decomp(gsl_rng * r)
{
  int s = 0;
  size_t n;

  /* test square matrices */
  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix_complex * m = gsl_matrix_complex_alloc(n, n);

      create_random_complex_matrix(m, r);
      test_LUc_decomp_eps(m, 1.0e4 * n * GSL_DBL_EPSILON, "complex_LU_decomp random");

      gsl_matrix_complex_free(m);
    }

  {
    gsl_matrix_complex * m = gsl_matrix_complex_alloc(100, 50);
    create_random_complex_matrix(m, r);
    test_LUc_decomp_eps(m, 256.0 * n * GSL_DBL_EPSILON, "complex_LU_decomp rect1");
    gsl_matrix_complex_free(m);
  }

  {
    gsl_matrix_complex * m = gsl_matrix_complex_alloc(50, 100);
    create_random_complex_matrix(m, r);
    test_LUc_decomp_eps(m, 1.0e3 * n * GSL_DBL_EPSILON, "complex_LU_decomp rect2");
    gsl_matrix_complex_free(m);
  }

  {
    gsl_matrix_complex * m = gsl_matrix_complex_alloc(80, 100);
    create_random_complex_matrix(m, r);
    test_LUc_decomp_eps(m, 1.0e3 * n * GSL_DBL_EPSILON, "complex_LU_decomp rect3");
    gsl_matrix_complex_free(m);
  }

  return s;
}

static int
test_LUc_solve_eps(const gsl_matrix_complex * m, const gsl_vector_complex * rhs, const gsl_vector_complex * sol, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = m->size1;
  int signum;
  size_t i;

  gsl_permutation * perm = gsl_permutation_alloc(N);
  gsl_matrix_complex * lu  = gsl_matrix_complex_alloc(N, N);
  gsl_vector_complex * x = gsl_vector_complex_alloc(N);
  gsl_vector_complex * residual = gsl_vector_complex_alloc(N);

  gsl_matrix_complex_memcpy(lu, m);

  s += gsl_linalg_complex_LU_decomp(lu, perm, &signum);
  s += gsl_linalg_complex_LU_solve(lu, perm, rhs, x);

  for (i = 0; i < N; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

      gsl_test_rel(GSL_REAL(xi), GSL_REAL(yi), eps, "%s: real %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, GSL_REAL(xi), GSL_REAL(yi));

      gsl_test_rel(GSL_IMAG(xi), GSL_IMAG(yi), eps, "%s: imag %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, GSL_IMAG(xi), GSL_IMAG(yi));
    }

  s += gsl_linalg_complex_LU_refine(m, lu, perm, rhs, x, residual);

  for (i = 0; i < N; i++)
    {
      gsl_complex xi = gsl_vector_complex_get(x, i);
      gsl_complex yi = gsl_vector_complex_get(sol, i);

      gsl_test_rel(GSL_REAL(xi), GSL_REAL(yi), eps, "%s: improved real %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, GSL_REAL(xi), GSL_REAL(yi));

      gsl_test_rel(GSL_IMAG(xi), GSL_IMAG(yi), eps, "%s: improved imag %3lu[%lu]: %22.18g    %22.18g\n",
                   desc, N, i, GSL_IMAG(xi), GSL_IMAG(yi));
    }

  gsl_vector_complex_free(residual);
  gsl_vector_complex_free(x);
  gsl_matrix_complex_free(lu);
  gsl_permutation_free(perm);

  return s;
}

static int
test_LUc_solve(gsl_rng * r)
{
  int s = 0;
  size_t n;

  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix_complex * m = gsl_matrix_complex_alloc(n, n);
      gsl_vector_complex * rhs = gsl_vector_complex_alloc(n);
      gsl_vector_complex * sol = gsl_vector_complex_alloc(n);

      create_random_complex_matrix(m, r);
      create_random_complex_vector(sol, r);
      gsl_blas_zgemv(CblasNoTrans, GSL_COMPLEX_ONE, m, sol, GSL_COMPLEX_ZERO, rhs);
      test_LUc_solve_eps(m, rhs, sol, 1.0e5 * n * GSL_DBL_EPSILON, "complex_LU_solve random");

      gsl_matrix_complex_free(m);
      gsl_vector_complex_free(rhs);
      gsl_vector_complex_free(sol);
    }

  return s;
}

int
test_LUc_invert_eps(const gsl_matrix_complex * m, const double eps, const char *desc)
{
  int s = 0;
  size_t i, j, N = m->size1;

  gsl_matrix_complex * lu  = gsl_matrix_complex_alloc(N, N);
  gsl_matrix_complex * c  = gsl_matrix_complex_alloc(N, N);
  gsl_permutation * perm = gsl_permutation_alloc(N);
  int signum;

  gsl_matrix_complex_memcpy(lu, m);

  s += gsl_linalg_complex_LU_decomp(lu, perm, &signum);
  s += gsl_linalg_complex_LU_invx(lu, perm);

  /* c = m m^{-1} */
  gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, GSL_COMPLEX_ONE, m, lu, GSL_COMPLEX_ZERO, c);

  /* c should be the identity matrix */

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < N; ++j)
        {
          gsl_complex cij = gsl_matrix_complex_get(c, i, j);
          double expected = (i == j) ? 1.0 : 0.0;

          gsl_test_rel(GSL_REAL(cij), expected, eps, "%s real (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_REAL(cij), expected);

          gsl_test_rel(GSL_IMAG(cij), 0.0, eps, "%s imag (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, N, i, j, GSL_IMAG(cij), 0.0);
        }
    }

  gsl_matrix_complex_free(lu);
  gsl_matrix_complex_free(c);
  gsl_permutation_free(perm);

  return s;
}

static int
test_LUc_invert(gsl_rng * r)
{
  int s = 0;
  size_t n;

  for (n = 1; n <= 50; ++n)
    {
      gsl_matrix_complex * m = gsl_matrix_complex_alloc(n, n);

      create_random_complex_matrix(m, r);
      test_LUc_invert_eps(m, 1.0e5 * n * GSL_DBL_EPSILON, "complex_LU_invert random");

      gsl_matrix_complex_free(m);
    }

  return s;
}

max_stars_repo_name: igormcoelho/optstats

max_issues_repo_name: igormcoelho/optstats

max_forks_repo_name: igormcoelho/optstats


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <time.h>

#define N_OF_RUNS 10

#define RUSAGE_DIFF(r1, r2) \
    ((double) ((r1).tv_sec - (r2).tv_sec) + (double) ((r1).tv_nsec - (r2).tv_nsec)/(double)1E9)

gsl_vector* random_vector(size_t);
gsl_matrix* random_matrix(size_t);
double* ddot_test(size_t);
double* dgemv_test(size_t);
void test_suite(size_t, size_t, size_t, const char*);

int main() {
    srand48((unsigned int)time(NULL));
    const char filename[] = "results.csv";
    test_suite(1, 15000, 1000, filename);
    return 0;
}

void test_suite(const size_t from, const size_t to, const size_t step, const char* filename){
    FILE* file = fopen(filename, "w");
    if(file == NULL){
        printf("Error opening file, exiting \r\n");
        exit(EXIT_FAILURE);
    }
    fprintf(file, "LEVEL,SIZE,TIME\r\n");

    double* ddot_result, * dgemv_result;

    for(size_t i = from; i < to; i += step){
        for(int j = 0; j < N_OF_RUNS; j++) {
            ddot_result = ddot_test(i);
            dgemv_result = dgemv_test(i);
            fprintf(file, "%d,%d,%lf\r\n", 1, (int) i, *ddot_result);
            fprintf(file, "%d,%d,%lf\r\n", 2, (int) i, *dgemv_result);
            free(ddot_result);
            free(dgemv_result);
        }
    }

    fclose(file);
}

gsl_vector* random_vector(const size_t size) {
    gsl_vector* a = gsl_vector_calloc(size);
    for(int i = 0; i < size; i++)
        gsl_vector_set(a, (size_t)i, drand48()*10);
    return a;
}

gsl_matrix* random_matrix(const size_t size) {
    gsl_matrix* a = gsl_matrix_calloc(size, size);
    for(int i = 0; i < size; i++)
        for(int j = 0; j < size; j++)
            gsl_matrix_set(a, (size_t)i, (size_t)j, drand48()*10);
    return a;
}

double* ddot_test(const size_t size) {
    gsl_vector* v_1 = random_vector(size);
    gsl_vector* v_2 = random_vector(size);
    printf("Multiplying 2 vectors of size %d \r\n", (int)size);
    double * result = calloc(1, sizeof(double));
    struct timespec real_start, real_end;
    clock_gettime(CLOCK_REALTIME, &real_start);
    for(int i = 0; i < N_OF_RUNS; i++) {
        gsl_blas_ddot(v_1, v_2, result);
    }
    clock_gettime(CLOCK_REALTIME, &real_end);
    double* time = malloc(sizeof(double));
    *time = RUSAGE_DIFF(real_end, real_start)/10;
    printf("Real time = %lf \r\n", *time);
    free(result);
    gsl_vector_free(v_1); gsl_vector_free(v_2);
    return time;
}

double* dgemv_test(const size_t size) {
    gsl_vector* v_1 = random_vector(size);
    gsl_vector* v_2 = random_vector(size);
    gsl_matrix* m_1 = random_matrix(size);
    printf("Multiplying vector times matrix of size %d \r\n", (int)size);

    struct timespec real_start, real_end;
    clock_gettime(CLOCK_REALTIME, &real_start);

    gsl_blas_dgemv(CblasNoTrans, 1, m_1, v_1, 0, v_2);

    clock_gettime(CLOCK_REALTIME, &real_end);

    double* time = malloc(sizeof(double));
    *time = RUSAGE_DIFF(real_end, real_start);
    printf("Real time = %lf \r\n", *time);
    gsl_vector_free(v_1); gsl_matrix_free(m_1);
    return time;
}


text:
#pragma once

#include "Graphics/BufferObject.h"
#include "Graphics/GraphicsResource.h"
#include "Math/Bounds.h"

#include <glad/gl.h>
#include <glm/glm.hpp>
#include <gsl/span>

#include <vector>

struct DrawingContext;

template<typename T>
struct MeshAttributeData
{
   gsl::span<T> values;
   GLint valueSize = 0;
};

struct MeshData
{
   gsl::span<GLuint> indices;
   MeshAttributeData<GLfloat> positions;
   MeshAttributeData<GLfloat> normals;
   MeshAttributeData<GLfloat> texCoords;
   MeshAttributeData<GLfloat> tangents;
   MeshAttributeData<GLfloat> bitangents;
   MeshAttributeData<GLfloat> colors;
};

class MeshSection : public GraphicsResource
{
public:
   MeshSection();
   MeshSection(const MeshSection& other) = delete;
   MeshSection(MeshSection&& other);
   ~MeshSection();
   MeshSection& operator=(const MeshSection& other) = delete;
   MeshSection& operator=(MeshSection&& other);

private:
   void move(MeshSection&& other);
   void release();

public:
   void setData(const MeshData& data);
   void draw(const DrawingContext& context) const;

   const Bounds& getBounds() const
   {
      return bounds;
   }

   void setLabel(std::string newLabel);

private:
   void bind() const;

   BufferObject elementBufferObject;
   VertexBufferObject positionBufferObject;
   VertexBufferObject normalBufferObject;
   VertexBufferObject texCoordBufferObject;
   VertexBufferObject tangentBufferObject;
   VertexBufferObject bitangentBufferObject;
   VertexBufferObject colorBufferObject;

   GLsizei numIndices;

   Bounds bounds;
};

class Mesh
{
public:
   Mesh(std::vector<MeshSection>&& meshSections);

   void draw(const DrawingContext& context) const;

   std::vector<MeshSection>& getSections()
   {
      return sections;
   }

   const std::vector<MeshSection>& getSections() const
   {
      return sections;
   }

private:
   std::vector<MeshSection> sections;
};

max_stars_repo_name: aaronmjacobs/Swap

max_issues_repo_name: aaronmjacobs/Swap

max_forks_repo_name: aaronmjacobs/Swap