text:
/**
 * @copyright (c) 2017 King Abdullah University of Science and Technology (KAUST).
 *                     All rights reserved.
 **/
/**
 * @file codelet_zhagcm.c
 *
 *  HiCMA codelets kernel
 *  HiCMA is a software package provided by King Abdullah University of Science and Technology (KAUST)
 *
 * @version 0.1.1
 * @author Kadir Akbudak
 * @date 2018-11-08
 * @precisions normal z -> c d s
 **/
#include "morse.h"
#include "runtime/starpu/chameleon_starpu.h"
//#include "runtime/starpu/include/runtime_codelet_z.h"

#include "runtime/starpu/runtime_codelets.h"
ZCODELETS_HEADER(hagcm)

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>  //FIXME remove mallocs from this code
#include <sys/time.h>//FIXME for gettimeofday

#include "hicma.h"
#include "starsh.h"
#include "starsh-spatial.h"
#include "starsh-randtlr.h"
#ifdef MKL
  #include <mkl.h>
  #include <mkl_lapack.h>
  //#pragma message("MKL is used")
#else
  #ifdef ARMPL
    #include <armpl.h>
  #else
    #include <cblas.h>
  #endif
  #ifdef LAPACKE_UTILS
    #include <lapacke_utils.h>
  #endif
  #include <lapacke.h>
  //#pragma message("MKL is NOT used")
#endif

//#warning "An experimental feature is enabled!!!"
int steal_lrtile = 0; //non-zero values are for experimental reasons. Otherwise set to 0!!!!!

extern void _printmat(double * A, int m, int n, int ld);

void zhagcm( int m, int n, /*dimension of squareAD*/
        double *AU,
        double *AV,
        double *Ark,
        int ldu,
        int ldv,
        int tile_row_index, int tile_column_index,
        int maxrank, double tol,
        int A_mt)
{
    int ii = tile_row_index;
    int jj = tile_column_index;
    if(steal_lrtile == 1 && ii == jj){
        if(ii == A_mt-1) { // steal tile above
            ii = ii - 1;
        } else { // still tile below
            ii = ii + 1;
        }
    }
    int64_t i, j;
    int lda = ldu; //FIXME ASSUMPTION
    struct timeval tvalBefore, tvalAfter;  // removed comma
    gettimeofday (&tvalBefore, NULL);

    int shape[2];
    int rank = 0;
    int oversample = 10;
    double *work;
    int *iwork;
    STARSH_cluster *RC = HICMA_get_starsh_format()->row_cluster, *CC = RC;
    void *RD = RC->data, *CD = RD;
    double *AD;
    AD = malloc(sizeof(double) * lda * n);

    HICMA_get_starsh_format()->problem->kernel(m, n, RC->pivot+RC->start[ii], CC->pivot+CC->start[jj],
            RD, CD, AD, lda);
    int mn = m;
    int mn2 = maxrank+oversample;
    if(mn2 > mn)
        mn2 = mn;
    // Get size of temporary arrays
    size_t lwork = n, lwork_sdd = (4*mn2+7)*mn2;
    if(lwork_sdd > lwork)
        lwork = lwork_sdd;
    lwork += (size_t)mn2*(2*n+m+mn2+1);
    size_t liwork = 8*mn2;
    // Allocate temporary arrays
    //STARSH_MALLOC(iwork, liwork);
    iwork = malloc(sizeof(*iwork) * liwork);
    if(iwork == NULL) {
        fprintf(stderr, "%s %s %d:\t Allocation failed. No memory! liwork:%d", __FILE__, __func__, __LINE__, liwork);
        exit(-1);
    }
    //STARSH_MALLOC(work, lwork);
    work = malloc(sizeof(*work) * lwork);
    if(work == NULL) {
        fprintf(stderr, "%s %s %d:\t Allocation failed. No memory! lwork:%d", __FILE__, __func__, __LINE__, lwork);
        exit(-1);
    }
    starsh_dense_dlrrsdd(m, n, AD, lda, AU, ldu,  AV, ldv, &rank, maxrank, oversample, tol, work, lwork, iwork);

    Ark[0] = rank;

    free(work);
    free(iwork);
}
/**
 * HICMA_TASK_zhagcm - Generate compressed matrix from a problem determined according to current global setting of HiCMA library
 */
void HICMA_TASK_zhagcm( const MORSE_option_t *options,
                        int m, int n,
                        const MORSE_desc_t *AUV,
                        const MORSE_desc_t *Ark,
                        int Am, int An,
                        int ldu,
                        int ldv,
                        int maxrank, double tol,
                        int A_mt
                        )
{
    struct starpu_codelet *codelet = &cl_zhagcm;
    //void (*callback)(void*) = options->profiling ? cl_zhagcm_callback : NULL;
    void (*callback)(void*) = NULL;
    int nAUV = AUV->nb;

    MORSE_BEGIN_ACCESS_DECLARATION;
    MORSE_ACCESS_W(AUV, Am, An);
    MORSE_ACCESS_W(Ark, Am, An);
    MORSE_END_ACCESS_DECLARATION;

    //printf("%s:%d: Am:%d An:%d lda:%d bigM:%d m0:%d n0:%d\n ", __FILE__, __LINE__, Am, An, lda, bigM, m0, n0);
        //printf("%s %d: Am:%d An:%d ADm:%d ADn:%d ptr:%p\n", __func__, __LINE__, Am, An, ADm, ADn, ptr);
    starpu_insert_task(
            starpu_mpi_codelet(codelet),
            STARPU_VALUE,    &m,                      sizeof(int),
            STARPU_VALUE,    &n,                      sizeof(int),
            STARPU_VALUE,    &nAUV,                      sizeof(int),
            STARPU_W,         RTBLKADDR(AUV, double, Am, An),
            STARPU_W,         RTBLKADDR(Ark, double, Am, An),
            STARPU_VALUE,  &ldu,                      sizeof(int),
            STARPU_VALUE,  &ldv,                      sizeof(int),
            STARPU_VALUE,   &Am,                      sizeof(int),
            STARPU_VALUE,   &An,                      sizeof(int),
            STARPU_VALUE,   &maxrank,                      sizeof(int),
            STARPU_VALUE,   &tol,                      sizeof(double),
            STARPU_VALUE,   &A_mt,                      sizeof(int),
            STARPU_PRIORITY,    options->priority,
            STARPU_CALLBACK,    callback,
#if defined(CHAMELEON_CODELETS_HAVE_NAME)
            STARPU_NAME, "zhagcm",
#endif
            0);
}

/*   cl_zhagcm_cpu_func - Generate a tile for random matrix. */

#if !defined(CHAMELEON_SIMULATION)
static void cl_zhagcm_cpu_func(void *descr[], void *cl_arg)
{
    int m;
    int n;
    int nAUV;
    double *AUV;
    double *Ark;
    int ldu;
    int ldv;
    int tile_row_index;
    int tile_column_index;
    int maxrank;
    double tol;
    int A_mt;

    AUV = (double *)STARPU_MATRIX_GET_PTR(descr[0]);
    Ark = (double *)STARPU_MATRIX_GET_PTR(descr[1]);

    starpu_codelet_unpack_args(cl_arg, &m, &n, &nAUV, &ldu, &ldv, &tile_row_index, &tile_column_index, &maxrank, &tol, &A_mt);

    double *AU = AUV;
    int nAU = nAUV/2;
    assert(ldu == ldv);
    size_t nelm_AU = (size_t)ldu * (size_t)nAU;
    double *AV = &(AUV[nelm_AU]);

    //printf("(%d,%d)%d %s %d %d\n", m0/m,n0/n,MORSE_My_Mpi_Rank(), __func__, __LINE__, AD == Dense);
    zhagcm( m, n,
            AU,
            AV,
            Ark,
            ldu,
            ldv,
            tile_row_index, tile_column_index,
            maxrank, tol, A_mt
            );
}
#endif /* !defined(CHAMELEON_SIMULATION) */

/*
 * Codelet definition
 */
CODELETS_CPU(zhagcm, 2, cl_zhagcm_cpu_func)


text:
/**
 * @file batchv_zgemv.c
 *
 * Part of API test for Batched BLAS routines.
 *
 * @author Samuel D. Relton
 * @author Pedro   V. Lara
 * @author Mawussi Zounon
 * @date 2016-06-01
 *
 * @precisions normal z -> c d s
 *
 **/

#include <cblas.h>
#include "bblas.h"

#define COMPLEX
void batchv_zgemv(
	const enum BBLAS_TRANS *trans,
	const int *m, const int *n,
	const BBLAS_Complex64_t *alpha,
	const BBLAS_Complex64_t **arrayA, const int *lda,
	const BBLAS_Complex64_t **arrayx, const int *incx,
	const BBLAS_Complex64_t *beta,
	BBLAS_Complex64_t **arrayy, const int *incy,
	const int batch_count, int* info)
{
	/* Local variables */
//	int first_index = 0;
	int batch_iter = 0;
	char func_name[15] = "batchv_zgemv";

	if (batch_count < 0)
	{
		xerbla_batch(func_name, BBLAS_ERR_BATCH_COUNT, -1);
	}

	for (batch_iter = 0; batch_iter < batch_count; batch_iter++)
	{
		/* Check input arguments */
		if ((trans[batch_iter] != BblasTrans) &&
			(trans[batch_iter] != BblasNoTrans) &&
			(trans[batch_iter] != BblasConjTrans))
		{
			xerbla_batch(func_name, BBLAS_ERR_TRANS, batch_iter);
			info[batch_iter] = BBLAS_ERR_TRANS;
		}

		if (m[batch_iter] < 0)
		{
			xerbla_batch(func_name, BBLAS_ERR_M, batch_iter);
			info[batch_iter] = BBLAS_ERR_M;
		}

		if (n[batch_iter] < 0)
		{
			xerbla_batch(func_name, BBLAS_ERR_N, batch_iter);
			info[batch_iter] = BBLAS_ERR_N;
		}

		/* Column major */
		if ((lda[batch_iter] < 1) && (lda[batch_iter] < m[batch_iter]))
		{
			xerbla_batch(func_name, BBLAS_ERR_LDA, batch_iter);
			info[batch_iter] = BBLAS_ERR_LDA;
		}

		if (incx[batch_iter] < 1)
		{
			xerbla_batch(func_name, BBLAS_ERR_INCX, batch_iter);
			info[batch_iter] = BBLAS_ERR_INCX;
		}

		if (incy[batch_iter] < 1)
		{
			xerbla_batch(func_name, BBLAS_ERR_INCY, batch_iter);
			info[batch_iter] = BBLAS_ERR_INCY;
		}

		/* Call CBLAS */
		cblas_zgemv(
			BblasColMajor,
			trans[batch_iter],
			m[batch_iter], n[batch_iter],
			CBLAS_SADDR( alpha[batch_iter] ),
			arrayA[batch_iter], lda[batch_iter],
			arrayx[batch_iter], incx[batch_iter],
			CBLAS_SADDR( beta[batch_iter] ),
			arrayy[batch_iter], incy[batch_iter]);
		/* Successful */
		info[batch_iter] = BBLAS_SUCCESS;
	} /* End variable size for loop */
}
#undef COMPLEX

ext: c

max_stars_repo_name: mawussi/BBLAS-group


text:
#ifndef GrowthFunctionBase_H
#define GrowthFunctionBase_H

#include <stdio.h>
#include <iostream>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <array>
#include <vector>

class GrowthFunctionBase{
private:

public:
	GrowthFunctionBase(int id, int type, float initTime, float endTime, bool applyToColumnarLayer, bool applyToPeripodialMembrane, bool applyToBasalECM, bool applyToLateralECM){
		/**
		 *  integer id will set GrowthFunctionBase#Id. \n
		 *  integer type will set GrowthFunctionBase#Type. \n
		 *  floats initTime and endTime will set GrowthFunctionBase#initTime and GrowthFunctionBase#endTime respectively. \n
		 *	booleans applyToColumnarLayer and applyToPeripodialMembrane will set GrowthFunctionBase#applyToColumnarLayer and GrowthFunctionBase#applyToPeripodialMembrane, respectively. \n
		 *   */
		this->Id = id;
		this->Type = type;
		this->initTime = initTime;
		this->endTime = endTime;
		this->applyToColumnarLayer = applyToColumnarLayer;
		this->applyToPeripodialMembrane = applyToPeripodialMembrane;
		this->applyToBasalECM = applyToBasalECM;
		this->applyToLateralECM = applyToLateralECM;
		zMin = 0;
		zMax = 1.0;
	} ///< The constructor of GrowthFunctionBase. Different growth functions will be derived from this class
    ~GrowthFunctionBase(){}

	int Type; 						///< The type of the growth function, 1: uniform growth, 2: Ring shaped growth, 3: Grid based growth, where growth rates read from a separate input file
	int Id;							///< The unique identification number of the growth function
	float initTime;					///< The initiation time of the growth, in seconds
	float endTime;					///< The end time of the growth, in seconds.
	bool applyToColumnarLayer;		///< Boolean stating if the growth should be applied to columnar layer
	bool applyToPeripodialMembrane; ///< Boolean stating if the growth should be applied to peripodial membrane
	bool applyToBasalECM;			///< Boolean stating if the growth should be applied to basal explicit ECM layer
	bool applyToLateralECM ;		///< Boolean stating if the growth should be applied to lateral explicit ECM layer
    bool applyTissueApical;         ///< Boolean stating if the growth should be applied to apical side of the tissue
    bool applyTissueBasal;          ///< Boolean stating if the growth should be applied to basal side of the tissue
    bool applyTissueMidLine;        ///< Boolean stating if the growth should be applied to mid-line side of the tissue
    double zMin;                    ///< Double giving the minimum z fraction of the tissue height that the growth is applied to
    double zMax;                    ///< Double giving the maximum z fraction of the tissue height that the growth is applied to
    std::vector <int> appliedEllipseBandIds;    ///< The vector of integers holding the iIDs of the marker ellipse bands that the growth isapplied to

    void 	ParentErrorMessage(std::string functionName){
        std::cerr<<"You are calling the function: "<<functionName<<" from a parent here, check declaration is via pointers"<<std::endl;
	}
    double 	ParentErrorMessage(std::string functionName, double returnValue){
        std::cerr<<"You are calling the function: "<<functionName<<" from a parent here, check declaration is via pointers"<<std::endl;
		return returnValue;
	}
    int 	ParentErrorMessage(std::string functionName, int returnValue){
        std::cerr<<"You are calling the function: "<<functionName<<" from a parent here, check declaration is via pointers"<<std::endl;
		return returnValue;
	}
    virtual void		writeSummary(std::ofstream&/*saveFileSimulationSummary*/){ParentErrorMessage("writeSummary_no_dt");} /// Virtual function of the parent to write the summary of growthFunction.
    virtual void		writeSummary(std::ofstream&/*saveFileSimulationSummary*/, double /*dt*/){ParentErrorMessage("writeSummary_with_dt");} /// Virtual function of the parent to write the summary of growthFunction.
    virtual void 		getCentre(float &/*centreX*/, float &/*centreY*/){ParentErrorMessage("getCentre");} /// Virtual function of the parent to get centre of the growth function definition.
    virtual float 		getInnerRadius(){return ParentErrorMessage("getInnerRadius",0.0);} /// Virtual function of the parent to get inner radius of the growth function definition.
    virtual float 		getOuterRadius(){return ParentErrorMessage("getOuterRadius",0.0);} /// Virtual function of the parent to get outer radius of the growth function definition.
    virtual std::array<double,3> 		getGrowthRate(){ParentErrorMessage("getGrowthRate");return std::array<double,3>{0.0};} /// Virtual function of the parent to get growth rate of the growth function definition.
    virtual std::array<double,3> 		getShapeChangeRateRate(){ParentErrorMessage("getShapeChangeRateRate");return std::array<double,3>{0.0};} /// Virtual function of the parent to get shape change rate of the growth function definition.
    virtual gsl_matrix* getShearAngleRotationMatrix(){ParentErrorMessage("getShearAngleRotationMatrix"); gsl_matrix* dummy = gsl_matrix_calloc(0,0); return dummy;}  /// Virtual function of the parent to get rotation angle matrix of the growth function definition.
    virtual double 		getShearAngle(){ParentErrorMessage("getShearAngle");return 0.0;} /// Virtual function of the parent to get rotation angle of the growth function definition.
    virtual size_t      getGridX(){return ParentErrorMessage("getGridX",0);} /// Virtual function of the parent to get gridX of the growth function definition.
    virtual size_t      getGridY(){return ParentErrorMessage("getGridY",0);} /// Virtual function of the parent to get gridY of the growth function definition.
    virtual	double 		getGrowthMatrixElement(int /*i*/, int /*j*/, int /*k*/){return ParentErrorMessage("getGrowthMatrixElement",0.0);} /// Virtual function of the parent to get grid matrix element of the growth function definition.
    virtual	double 		getXyShearAngleMatrixElement(int /*i*/, int /*j*/){return ParentErrorMessage("getXyShearhMatrixElement",0.0);} /// Virtual function of the parent to get grid matrix element of the growth function definition.
    virtual bool 		isAspectRatioOverOne(int /*i*/, int /*j*/){return ParentErrorMessage("isAspectRatioOverOne",0);}/// Virtual function of the parent to check aspect ratioof the growth function definition.
    virtual gsl_matrix* getXyShearRotationsMatrixElement(int /*i*/, int /*j*/){ParentErrorMessage("getShearAngleRotationMatrixElement");gsl_matrix* dummy = gsl_matrix_calloc(0,0); return dummy;} //this is used by grid based growth
    virtual void 		getGrowthProfileAt4Corners(int /*IndexX*/, int /*IndexY*/, std::array<double,3>& /*growth0*/, std::array<double,3>&/*growth1*/, std::array<double,3>&/*growth2*/, std::array<double,3>&/*growth3*/, std::array<double,4>&/*angles*/, std::array<bool,4>&/*anglesEliminated*/){ParentErrorMessage("getGrowthProfileAt4Corners");}

    virtual void		setGrowtRate(double /*ex*/, double /*ey*/, double /*ez*/){ParentErrorMessage("setGrowtRate");}/// Virtual function of the parent to set growth rateof the growth function definition.
    virtual void		setGrowthMatrixElement(double /*ex*/, double /*ey*/, double /*ez*/, int /*i*/, int /*j*/){ParentErrorMessage("setGrowtRate");} /// Virtual function of the parent to set growth rateof the growth function definition.
};
#endif

ext: h


text:
#ifndef ITIMING_H
#define ITIMING_H

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <fcntl.h>
#include <chameleon.h>
#include <al4san.h>
#include <runtime/al4san_quark.h>
#define _TYPE  double
#define _PREC  double
#define _LAMCH LAPACKE_dlamch_work

#define _NAME  "sygv"
/* See Lawn 41 page 120 */
/* cholesky + 2 trsm's + trd */
#define _FMULS (((2. / 3.) * ((double)N * (double)N * (double)N)) + (N * (1.0 / 6.0 * N + 0.5) * N) + 2 * (N * N * (double)((N + 1) / 2.0) ))
#define _FADDS (((2. / 3.) * ((double)N * (double)N * (double)N)) + (N * (1.0 / 6.0 * N )      * N) + 2 * (N * N * (double)((N + 1) / 2.0) ))

/*
#define _FMULS ((double)N * (((1. / 6.) * (double)N + 0.5) * (double)N + (1. / 3.)))
#define _FADDS ((double)N * (((1. / 6.) * (double)N      ) * (double)N - (1. / 6.)))
 */
#include "../include/flops.h"

#include<mkl_lapacke.h>
//#include <lapacke.h>
//#include <core_blas.h>
//#include <cblas.h>
#include<mkl.h>
#if defined(EIG_SCHED_STARPU)
#include <starpu.h>
#endif /* defined(EIG_SCHED_STARPU) */
#if defined(EIG_SCHED_QUARK)
#include <quark.h>
#endif /* defined(EIG_SCHED_QUARK) */
#if defined(EIG_SCHED_PARSEC)
#include <parsec.h>
#endif /* defined(EIG_SCHED_PARSEC) */

//#include <plasma.h>
#include <quark.h>
//#include "plasmatypes.h"
//#include "descriptor.h"

typedef double real_Double_t;

#define EIG_INPLACE    1
#define EIG_OUTOFPLACE 2

#if defined( _WIN32 ) || defined( _WIN64 )
#include <windows.h>
#include <time.h>
#include <sys/timeb.h>
#if defined(_MSC_VER) || defined(_MSC_EXTENSIONS)
#define DELTA_EPOCH_IN_MICROSECS  11644473600000000Ui64
#else
#define DELTA_EPOCH_IN_MICROSECS  11644473600000000ULL
#endif

double cWtime(void);

struct timezone
{
	int  tz_minuteswest; /* minutes W of Greenwich */
	int  tz_dsttime;     /* type of dst correction */
};

int gettimeofday(struct timeval* tv, struct timezone* tz)
{
	FILETIME         ft;
	unsigned __int64 tmpres = 0;
	static int       tzflag;

	if (NULL != tv)
	{
		GetSystemTimeAsFileTime(&ft);
		tmpres |=  ft.dwHighDateTime;
		tmpres <<= 32;
		tmpres |=  ft.dwLowDateTime;

		/*converting file time to unix epoch*/
		tmpres /= 10;  /*convert into microseconds*/
		tmpres -= DELTA_EPOCH_IN_MICROSECS;

		tv->tv_sec  = (long)(tmpres / 1000000UL);
		tv->tv_usec = (long)(tmpres % 1000000UL);
	}
	if (NULL != tz)
	{
		if (!tzflag)
		{
			_tzset();
			tzflag++;
		}
		tz->tz_minuteswest = _timezone / 60;
		tz->tz_dsttime     = _daylight;
	}
	return 0;
}

#else  /* Non-Windows */
#include <unistd.h>
#include <sys/time.h>
#include <sys/resource.h>
#endif

enum iparam_t {
	IPARAM_THRDNBR,        /* Number of cores                            */
	IPARAM_THRDNBR_SUBGRP, /* Number of cores in a subgroup (NUMA node)  */
	IPARAM_NCUDAS,
	IPARAM_NMPI,
	IPARAM_P,
	IPARAM_Q,
	IPARAM_SCHEDULER,      /* What scheduler do we choose (dyn, stat)    */
	IPARAM_M,              /* Number of rows of the matrix               */
	IPARAM_N,              /* Number of columns of the matrix            */
	IPARAM_K,              /* RHS or K                                   */
	IPARAM_LDA,            /* Leading dimension of A                     */
	IPARAM_LDB,            /* Leading dimension of B                     */
	IPARAM_LDC,            /* Leading dimension of C                     */
	IPARAM_IB,             /* Inner-blocking size                        */
	IPARAM_NB,             /* Number of columns in a tile                */
	IPARAM_MB,             /* Number of rows in a tile                   */
	IPARAM_NITER,          /* Number of iteration of each test           */
	IPARAM_WARMUP,         /* Run one test to load dynamic libraries     */
	IPARAM_CHECK,          /* Checking activated or not                  */
	IPARAM_VERBOSE,        /* How much noise do we want?                 */
	IPARAM_AUTOTUNING,     /* Disable/enable autotuning                  */
	IPARAM_INPUTFMT,       /* Input format (Use only for getmi/gecfi)    */
	IPARAM_OUTPUTFMT,      /* Output format (Use only for getmi/gecfi)   */
	IPARAM_TRACE,          /* Generate trace on the first non warmup run */
	IPARAM_DAG,            /* Do we require to output the DOT file?      */
	IPARAM_ASYNC,          /* Asynchronous calls                         */
	IPARAM_MX,             /* */
	IPARAM_NX,             /* */
	IPARAM_RHBLK,          /* Householder reduction parameter for QR/LQ  */
	IPARAM_TNTPIV_MODE,    /* Tournament pivoting LU mode (LU or QR)     */
	IPARAM_TNTPIV_SIZE,    /* Tournament pivoting LU group size          */
	IPARAM_INPLACE,        /* InPlace/OutOfPlace translation mode        */
	IPARAM_MODE,           /* Eigenvalue generation mode                 */
	IPARAM_LAPACK,         /* Eigenvalue dstedc using or not LAPACK      */
	IPARAM_MATRIX,         /* FD matrix entry used for tridiagonal solvers */
	IPARAM_MATRIX_CHECK,   /* FD for checking results                    */
	IPARAM_SIZEOF,
	IPARAM_GEMM3M
};

enum dparam_timing {
	IPARAM_TIME,
	IPARAM_ANORM,
	IPARAM_BNORM,
	IPARAM_XNORM,
	IPARAM_RES,
	IPARAM_DNBPARAM
};
/*
   int min( int a, int b ){
   if (a<b) return(a); else return(b);
   }

   int max( int a, int b ){
   if (a>b) return(a); else return(b);
   }*/

#define PASTE_CODE_IPARAM_LOCALS(iparam)           \
	double  t;                                     \
double t1, t2;                                 \
int64_t M     = iparam[IPARAM_M];              \
int64_t N     = iparam[IPARAM_N];              \
int64_t K     = iparam[IPARAM_K];              \
int64_t NRHS  = K;                             \
int64_t LDA   = N;    \
int64_t LDB   = max(N, iparam[IPARAM_LDB]);    \
int64_t LDC   = max(K, iparam[IPARAM_LDC]);    \
int64_t IB    = iparam[IPARAM_IB];             \
int64_t MB    = iparam[IPARAM_MB];             \
int64_t NB    = iparam[IPARAM_NB];             \
int64_t MT    = (M%MB==0) ? (M/MB) : (M/MB+1); \
int64_t NT    = (N%NB==0) ? (N/NB) : (N/NB+1); \
int64_t P     = iparam[IPARAM_P];              \
int64_t Q     = iparam[IPARAM_Q];              \
int check     = iparam[IPARAM_CHECK];          \
int loud      = iparam[IPARAM_VERBOSE];        \
(void)M;(void)N;(void)K;(void)NRHS;            \
(void)LDA;(void)LDB;(void)LDC;                 \
(void)IB;(void)MB;(void)NB;(void)MT;(void)NT;  \
(void)check;(void)loud;

#define PASTE_CODE_IPARAM_LOCALS2(iparam)           \
M     = iparam[IPARAM_M];              \
N     = iparam[IPARAM_N];              \
K     = iparam[IPARAM_K];              \
NRHS  = K;                             \
LDA   = N;    \
LDB   = max(N, iparam[IPARAM_LDB]);    \
LDC   = max(K, iparam[IPARAM_LDC]);    \
IB    = iparam[IPARAM_IB];             \
MB    = iparam[IPARAM_MB];             \
NB    = iparam[IPARAM_NB];             \
MT    = (M%MB==0) ? (M/MB) : (M/MB+1); \
NT    = (N%NB==0) ? (N/NB) : (N/NB+1); \
P     = iparam[IPARAM_P];              \
Q     = iparam[IPARAM_Q];              \
check     = iparam[IPARAM_CHECK];          \
loud      = iparam[IPARAM_VERBOSE];        \
(void)M;(void)N;(void)K;(void)NRHS;            \
(void)LDA;(void)LDB;(void)LDC;                 \
(void)IB;(void)MB;(void)NB;(void)MT;(void)NT;  \
(void)check;(void)loud;

/* Paste code to allocate a matrix in desc if cond_init is true */
/*#define PASTE_CODE_ALLOCATE_MATRIX_TILE(_desc_, _cond_, _type_, _type2_, _lda_, _m_, _n_) \
	PLASMA_desc *_desc_ = NULL;                                         \
if( _cond_ ) {                                                      \
	_type_ *ptr = NULL;                                             \
	ptr = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_) );      \
	if ( ! ptr ) {                                                  \
		fprintf(stderr, "Our of Memory for %s\n", #_desc_);         \
		return -1;                                                  \
	}                                                               \
	PLASMA_Desc_Create(&(_desc_), ptr, _type2_, MB, NB, MB*NB, _lda_, _n_, 0, 0, _m_, _n_); \
}

#define PASTE_CODE_FREE_MATRIX(_desc_)                                  \
	if ( _desc_ != NULL ) {                                             \
		free(_desc_->mat);                                              \
	}                                                                   \
PLASMA_Desc_Destroy( &_desc_ );

#define PASTE_TILE_TO_LAPACK(_desc_, _name_, _cond_, _type_, _lda_, _n_) \
	_type_ *_name_ = NULL;                                              \
if ( _cond_ ) {                                                     \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_));    \
	if ( ! _name_ ) {                                               \
		fprintf(stderr, "Our of Memory for %s\n", #_name_);         \
		return -1;                                                  \
	}                                                               \
	PLASMA_Tile_to_Lapack(_desc_, (void*)_name_, _lda_);            \
}

#define PASTE_CODE_ALLOCATE_MATRIX(_name_, _cond_, _type_, _lda_, _n_)  \
	_type_ *_name_ = NULL;                                              \
if( _cond_ ) {                                                      \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_) );   \
	if ( ! _name_ ) {                                               \
		fprintf(stderr, "Our of Memory for %s\n", #_name_);         \
		return -1;                                                  \
	}                                                               \
}

#define PASTE_CODE_ALLOCATE_COPY(_name_, _cond_, _type_, _orig_, _lda_, _n_) \
	_type_ *_name_ = NULL;                                              \
if( _cond_ ) {                                                      \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_) );   \
	if ( ! _name_ ) {                                               \
		fprintf(stderr, "Our of Memory for %s\n", #_name_);         \
		return -1;                                                  \
	}                                                               \
	memcpy(_name_, _orig_, (_lda_) * (_n_) * sizeof(_type_) );      \
}
*/
/* Paste code to allocate a matrix in desc if cond_init is true */
#define CHAM_PASTE_CODE_ALLOCATE_MATRIX_TILE(_desc_, _cond_, _type_, _type2_, _lda_, _m_, _n_) \
	CHAM_desc_t *_desc_ = NULL;                                        \
int status ## _desc_ ;                                              \
if( _cond_ ) {                                                      \
	if (ooc)                                                        \
	status ## _desc_ = CHAMELEON_Desc_Create_OOC(&(_desc_), _type2_, MB, NB, MB*NB, _lda_, _n_, 0, 0, _m_, _n_, \
			P, Q);             \
	else if (!bigmat)                                               \
	status ## _desc_ = CHAMELEON_Desc_Create_User(&(_desc_), NULL, _type2_, MB, NB, MB*NB, _lda_, _n_, 0, 0, _m_, _n_, \
			P, Q, chameleon_getaddr_null, NULL, NULL); \
	else                                                            \
	status ## _desc_ = CHAMELEON_Desc_Create(&(_desc_), NULL, _type2_, MB, NB, MB*NB, _lda_, _n_, 0, 0, _m_, _n_, \
			P, Q);                 \
	if (status ## _desc_ != CHAMELEON_SUCCESS) return (status ## _desc_); \
}

#define CHAM_PASTE_CODE_FREE_MATRIX(_desc_)                                  \
	CHAMELEON_Desc_Destroy( &_desc_ );

#define CHAM_PASTE_TILE_TO_LAPACK(_desc_, _name_, _cond_, _type_, _lda_, _n_) \
	_type_ *_name_ = NULL;                                               \
if ( _cond_ ) {                                                      \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_));     \
	if ( ! _name_ ) {                                                \
		fprintf(stderr, "Out of Memory for %s\n", #_name_);          \
		return -1;                                                   \
	}                                                                \
	CHAMELEON_Tile_to_Lapack(_desc_, (void*)_name_, _lda_);              \
}

#define CHAM_PASTE_CODE_ALLOCATE_MATRIX(_name_, _cond_, _type_, _lda_, _n_)  \
	_type_ *_name_ = NULL;                                              \
if( _cond_ ) {                                                      \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_) );   \
	if ( ! _name_ ) {                                               \
		fprintf(stderr, "Out of Memory for %s\n", #_name_);         \
		return -1;                                                  \
	}                                                               \
}

#define CHAM_PASTE_CODE_ALLOCATE_COPY(_name_, _cond_, _type_, _orig_, _lda_, _n_) \
	_type_ *_name_ = NULL;                                                   \
if( _cond_ ) {                                                           \
	_name_ = (_type_*)malloc( (_lda_) * (_n_) * sizeof(_type_) );        \
	if ( ! _name_ ) {                                                    \
		fprintf(stderr, "Out of Memory for %s\n", #_name_);              \
		return -1;                                                       \
	}                                                                    \
	memcpy(_name_, _orig_, (_lda_) * (_n_) * sizeof(_type_) );           \
}

/*********************
 *
 * General Macros for timing
 *
 */
#define START_TIMING1()                          \
        t1=0;                                    \
	t1 = -cWtime();

#define STOP_TIMING1()                          \
        t1 += cWtime();

#define START_TIMING2()                          \
        t2=0;                                    \
        t2 = -cWtime();

#define STOP_TIMING2()                           \
	t2 += cWtime();                                \
*t_ = t1+t2;
#endif /* ITIMING_H */


text:
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sort_double.h>

int
main ()
{
  gsl_rng * r;

  int i, k = 5, N = 100000;

  double * x = malloc (N * sizeof(double));
  double * small = malloc (k * sizeof(double));

  gsl_rng_env_setup();

  r = gsl_rng_alloc (gsl_rng_default);

  for (i = 0; i < N; i++)
    {
      x[i] = gsl_rng_uniform(r);
    }

  gsl_sort_smallest (small, k, x, 1, N);

  printf("%d smallest values from %d\n", k, N);

  for (i = 0; i < k; i++)
    {
      printf ("%d: %.18f\n", i, small[i]);
    }

}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include "ccv.h"
#include "ccv_internal.h"
#ifdef HAVE_CBLAS
#include <cblas.h>
#elif HAVE_ACCELERATE_FRAMEWORK
#include <Accelerate/Accelerate.h>
#endif

double ccv_trace(ccv_matrix_t* mat)
{
	return 0;
}

double ccv_norm(ccv_matrix_t* mat, int type)
{
	return 0;
}

double ccv_normalize(ccv_matrix_t* a, ccv_matrix_t** b, int btype, int flag)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	assert(CCV_GET_CHANNEL(da->type) == CCV_C1);
	ccv_declare_derived_signature(sig, da->sig != 0, ccv_sign_with_format(20, "ccv_normalize(%d)", flag), da->sig, CCV_EOF_SIGN);
	btype = (btype == 0) ? CCV_GET_DATA_TYPE(da->type) | CCV_C1 : CCV_GET_DATA_TYPE(btype) | CCV_C1;
	ccv_dense_matrix_t* db = *b = ccv_dense_matrix_renew(*b, da->rows, da->cols, CCV_ALL_DATA_TYPE | CCV_C1, btype, sig);
	assert(db);
	ccv_object_return_if_cached(db->tag.f64, db);
	double sum = 0, inv;
	int i, j;
	unsigned char* a_ptr = da->data.u8;
	unsigned char* b_ptr = db->data.u8;
	switch (flag)
	{
		case CCV_L1_NORM:
#define for_block(_for_set, _for_get) \
			for (i = 0; i < da->rows; i++) \
			{ \
				for (j = 0; j < da->cols; j++) \
					sum += _for_get(a_ptr, j, 0); \
				a_ptr += da->step; \
			} \
			inv = 1.0 / sum; \
			a_ptr = da->data.u8; \
			for (i = 0; i < da->rows; i++) \
			{ \
				for (j = 0; j < da->cols; j++) \
					_for_set(b_ptr, j, _for_get(a_ptr, j, 0) * inv, 0); \
				a_ptr += da->step; \
				b_ptr += db->step; \
			}
			ccv_matrix_setter(db->type, ccv_matrix_getter, da->type, for_block);
#undef for_block
			break;
		case CCV_L2_NORM:
#define for_block(_for_set, _for_get) \
			for (i = 0; i < da->rows; i++) \
			{ \
				for (j = 0; j < da->cols; j++) \
					sum += _for_get(a_ptr, j, 0) * _for_get(a_ptr, j, 0); \
				a_ptr += da->step; \
			} \
			sum = sqrt(sum); \
			inv = 1.0 / sum; \
			a_ptr = da->data.u8; \
			for (i = 0; i < da->rows; i++) \
			{ \
				for (j = 0; j < da->cols; j++) \
					_for_set(b_ptr, j, _for_get(a_ptr, j, 0) * inv, 0); \
				a_ptr += da->step; \
				b_ptr += db->step; \
			}
			ccv_matrix_setter(db->type, ccv_matrix_getter, da->type, for_block);
#undef for_block
			break;
	}
	return db->tag.f64 = sum;
}

void ccv_sat(ccv_dense_matrix_t* a, ccv_dense_matrix_t** b, int type, int padding_pattern)
{
	ccv_declare_derived_signature(sig, a->sig != 0, ccv_sign_with_format(20, "ccv_sat(%d)", padding_pattern), a->sig, CCV_EOF_SIGN);
	int safe_type = (a->type & CCV_8U) ? ((a->rows * a->cols >= 0x808080) ? CCV_64S : CCV_32S) : ((a->type & CCV_32S) ? CCV_64S : a->type);
	type = (type == 0) ? CCV_GET_DATA_TYPE(safe_type) | CCV_GET_CHANNEL(a->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(a->type);
	int ch = CCV_GET_CHANNEL(a->type);
	int i, j;
	unsigned char* a_ptr = a->data.u8;
	ccv_dense_matrix_t* db;
	unsigned char* b_ptr;
	switch (padding_pattern)
	{
		case CCV_NO_PADDING:
			db = *b = ccv_dense_matrix_renew(*b, a->rows, a->cols, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(a->type), type, sig);
			ccv_object_return_if_cached(, db);
			b_ptr = db->data.u8;
#define for_block(_for_set_b, _for_get_b, _for_get) \
			for (j = 0; j < ch; j++) \
				_for_set_b(b_ptr, j, _for_get(a_ptr, j, 0), 0); \
			for (j = ch; j < a->cols * ch; j++) \
				_for_set_b(b_ptr, j, _for_get_b(b_ptr, j - ch, 0) + _for_get(a_ptr, j, 0), 0); \
			a_ptr += a->step; \
			b_ptr += db->step; \
			for (i = 1; i < a->rows; i++) \
			{ \
				for (j = 0; j < ch; j++) \
					_for_set_b(b_ptr, j, _for_get_b(b_ptr - db->step, j, 0) + _for_get(a_ptr, j, 0), 0); \
				for (j = ch; j < a->cols * ch; j++) \
					_for_set_b(b_ptr, j, _for_get_b(b_ptr, j - ch, 0) - _for_get_b(b_ptr - db->step, j - ch, 0) + _for_get_b(b_ptr - db->step, j, 0) + _for_get(a_ptr, j, 0), 0); \
				a_ptr += a->step; \
				b_ptr += db->step; \
			}
			ccv_matrix_setter_getter(db->type, ccv_matrix_getter, a->type, for_block);
#undef for_block
			break;
		case CCV_PADDING_ZERO:
			db = *b = ccv_dense_matrix_renew(*b, a->rows + 1, a->cols + 1, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(a->type), type, sig);
			ccv_object_return_if_cached(, db);
			b_ptr = db->data.u8;
#define for_block(_for_set_b, _for_get_b, _for_get) \
			for (j = 0; j < db->cols * ch; j++) \
				_for_set_b(b_ptr, j, 0, 0); \
			b_ptr += db->step; \
			for (i = 0; i < a->rows; i++) \
			{ \
				for (j = 0; j < ch; j++) \
					_for_set_b(b_ptr, j, 0, 0); \
				for (j = ch; j < db->cols * ch; j++) \
					_for_set_b(b_ptr, j, _for_get_b(b_ptr, j - ch, 0) - _for_get_b(b_ptr - db->step, j - ch, 0) + _for_get_b(b_ptr - db->step, j, 0) + _for_get(a_ptr, j - ch, 0), 0); \
				a_ptr += a->step; \
				b_ptr += db->step; \
			}
			ccv_matrix_setter_getter(db->type, ccv_matrix_getter, a->type, for_block);
#undef for_block
			break;
	}
}

double ccv_sum(ccv_matrix_t* mat, int flag)
{
	ccv_dense_matrix_t* dmt = ccv_get_dense_matrix(mat);
	double sum = 0;
	unsigned char* m_ptr = dmt->data.u8;
	int i, j, ch = CCV_GET_CHANNEL(dmt->type);
#define for_block(_, _for_get) \
	switch (flag) \
	{ \
		case CCV_UNSIGNED: \
			for (i = 0; i < dmt->rows; i++) \
			{ \
				for (j = 0; j < dmt->cols * ch; j++) \
					sum += fabs(_for_get(m_ptr, j, 0)); \
				m_ptr += dmt->step; \
			} \
			break; \
		case CCV_SIGNED: \
		default: \
			for (i = 0; i < dmt->rows; i++) \
			{ \
				for (j = 0; j < dmt->cols * ch; j++) \
					sum += _for_get(m_ptr, j, 0); \
				m_ptr += dmt->step; \
			} \
	}
	ccv_matrix_getter(dmt->type, for_block);
#undef for_block
	return sum;
}

double ccv_variance(ccv_matrix_t* mat)
{
	ccv_dense_matrix_t* dmt = ccv_get_dense_matrix(mat);
	double mean = 0, variance = 0;
	unsigned char* m_ptr = dmt->data.u8;
	int i, j, ch = CCV_GET_CHANNEL(dmt->type);
#define for_block(_, _for_get) \
	for (i = 0; i < dmt->rows; i++) \
	{ \
		for (j = 0; j < dmt->cols * ch; j++) \
		{ \
			mean += _for_get(m_ptr, j, 0); \
			variance += _for_get(m_ptr, j, 0) * _for_get(m_ptr, j, 0); \
		} \
		m_ptr += dmt->step; \
	}
	ccv_matrix_getter(dmt->type, for_block);
#undef for_block
	mean = mean / (dmt->rows * dmt->cols * ch);
	variance = variance / (dmt->rows * dmt->cols * ch);
	return variance - mean * mean;
}

void ccv_multiply(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** c, int type)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
	assert(da->rows == db->rows && da->cols == db->cols && CCV_GET_DATA_TYPE(da->type) == CCV_GET_DATA_TYPE(db->type) && CCV_GET_CHANNEL(da->type) == CCV_GET_CHANNEL(db->type));
	ccv_declare_derived_signature(sig, da->sig != 0 && db->sig != 0, ccv_sign_with_literal("ccv_multiply"), da->sig, db->sig, CCV_EOF_SIGN);
	int no_8u_type = (da->type & CCV_8U) ? CCV_32S : da->type;
	type = (type == 0) ? CCV_GET_DATA_TYPE(no_8u_type) | CCV_GET_CHANNEL(da->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(da->type);
	ccv_dense_matrix_t* dc = *c = ccv_dense_matrix_renew(*c, da->rows, da->cols, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(da->type), type, sig);
	ccv_object_return_if_cached(, dc);
	int i, j, ch = CCV_GET_CHANNEL(da->type);
	unsigned char* aptr = da->data.u8;
	unsigned char* bptr = db->data.u8;
	unsigned char* cptr = dc->data.u8;
#define for_block(_for_get, _for_set) \
	for (i = 0; i < da->rows; i++) \
	{ \
		for (j = 0; j < da->cols * ch; j++) \
			_for_set(cptr, j, _for_get(aptr, j, 0) * _for_get(bptr, j, 0), 0); \
		aptr += da->step; \
		bptr += db->step; \
		cptr += dc->step; \
	}
	ccv_matrix_getter(da->type, ccv_matrix_setter, dc->type, for_block);
#undef for_block
}

void ccv_add(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** c, int type)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
	assert(da->rows == db->rows && da->cols == db->cols && CCV_GET_CHANNEL(da->type) == CCV_GET_CHANNEL(db->type));
	ccv_declare_derived_signature(sig, da->sig != 0 && db->sig != 0, ccv_sign_with_literal("ccv_add"), da->sig, db->sig, CCV_EOF_SIGN);
	int no_8u_type = (da->type & CCV_8U) ? CCV_32S : da->type;
	type = (type == 0) ? CCV_GET_DATA_TYPE(no_8u_type) | CCV_GET_CHANNEL(da->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(da->type);
	ccv_dense_matrix_t* dc = *c = ccv_dense_matrix_renew(*c, da->rows, da->cols, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(da->type), type, sig);
	ccv_object_return_if_cached(, dc);
	int i, j, ch = CCV_GET_CHANNEL(da->type);
	unsigned char* aptr = da->data.u8;
	unsigned char* bptr = db->data.u8;
	unsigned char* cptr = dc->data.u8;
#define for_block(_for_get_a, _for_get_b, _for_set) \
	for (i = 0; i < da->rows; i++) \
	{ \
		for (j = 0; j < da->cols * ch; j++) \
			_for_set(cptr, j, _for_get_a(aptr, j, 0) + _for_get_b(bptr, j, 0), 0); \
		aptr += da->step; \
		bptr += db->step; \
		cptr += dc->step; \
	}
	ccv_matrix_getter_a(da->type, ccv_matrix_getter_b, db->type, ccv_matrix_setter, dc->type, for_block);
#undef for_block
}

void ccv_subtract(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** c, int type)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
	assert(da->rows == db->rows && da->cols == db->cols && CCV_GET_DATA_TYPE(da->type) == CCV_GET_DATA_TYPE(db->type) && CCV_GET_CHANNEL(da->type) == CCV_GET_CHANNEL(db->type));
	ccv_declare_derived_signature(sig, da->sig != 0 && db->sig != 0, ccv_sign_with_literal("ccv_subtract"), da->sig, db->sig, CCV_EOF_SIGN);
	int no_8u_type = (da->type & CCV_8U) ? CCV_32S : da->type;
	type = (type == 0) ? CCV_GET_DATA_TYPE(no_8u_type) | CCV_GET_CHANNEL(da->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(da->type);
	ccv_dense_matrix_t* dc = *c = ccv_dense_matrix_renew(*c, da->rows, da->cols, CCV_ALL_DATA_TYPE | CCV_GET_CHANNEL(da->type), type, sig);
	ccv_object_return_if_cached(, dc);
	int i, j, ch = CCV_GET_CHANNEL(da->type);
	unsigned char* aptr = da->data.u8;
	unsigned char* bptr = db->data.u8;
	unsigned char* cptr = dc->data.u8;
#define for_block(_for_get, _for_set) \
	for (i = 0; i < da->rows; i++) \
	{ \
		for (j = 0; j < da->cols * ch; j++) \
			_for_set(cptr, j, _for_get(aptr, j, 0) - _for_get(bptr, j, 0), 0); \
		aptr += da->step; \
		bptr += db->step; \
		cptr += dc->step; \
	}
	ccv_matrix_getter(da->type, ccv_matrix_setter, dc->type, for_block);
#undef for_block
}

void ccv_gemm(ccv_matrix_t* a, ccv_matrix_t* b, double alpha, ccv_matrix_t* c, double beta, int transpose, ccv_matrix_t** d, int type)
{
	ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
	ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
	ccv_dense_matrix_t* dc = (c == 0) ? 0 : ccv_get_dense_matrix(c);

	assert(CCV_GET_DATA_TYPE(da->type) == CCV_GET_DATA_TYPE(db->type) && CCV_GET_CHANNEL(da->type) == 1 && CCV_GET_CHANNEL(db->type) == 1 && ((transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols) == ((transpose & CCV_B_TRANSPOSE) ? db->cols : db->rows));

	if (dc != 0)
		assert(CCV_GET_DATA_TYPE(dc->type) == CCV_GET_DATA_TYPE(da->type) && CCV_GET_CHANNEL(dc->type) == 1 && ((transpose & CCV_A_TRANSPOSE) ? da->cols : da->rows) == dc->rows && ((transpose & CCV_B_TRANSPOSE) ? db->rows : db->cols) == dc->cols);

	ccv_declare_derived_signature_case(sig, ccv_sign_with_format(20, "ccv_gemm(%d)", transpose), ccv_sign_if(dc == 0 && da->sig != 0 && db->sig != 0, da->sig, db->sig, CCV_EOF_SIGN), ccv_sign_if(dc != 0 && da->sig != 0 && db->sig != 0 && dc->sig != 0, da->sig, db->sig, dc->sig, CCV_EOF_SIGN));
	type = CCV_GET_DATA_TYPE(da->type) | CCV_GET_CHANNEL(da->type);
	ccv_dense_matrix_t* dd = *d = ccv_dense_matrix_renew(*d, (transpose & CCV_A_TRANSPOSE) ? da->cols : da->rows, (transpose & CCV_B_TRANSPOSE) ? db->rows : db->cols, type, type, sig);
	ccv_object_return_if_cached(, dd);

	if (dd != dc && dc != 0)
		memcpy(dd->data.u8, dc->data.u8, dc->step * dc->rows);
	else if (dc == 0) // clean up dd if dc is not provided
		memset(dd->data.u8, 0, dd->step * dd->rows);

#if (defined HAVE_CBLAS || defined HAVE_ACCELERATE_FRAMEWORK)
	switch (CCV_GET_DATA_TYPE(dd->type))
	{
		case CCV_32F:
			cblas_sgemm(CblasRowMajor, (transpose & CCV_A_TRANSPOSE) ? CblasTrans : CblasNoTrans, (transpose & CCV_B_TRANSPOSE) ? CblasTrans : CblasNoTrans, dd->rows, dd->cols, (transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols, alpha, da->data.f32, da->cols, db->data.f32, db->cols, beta, dd->data.f32, dd->cols);
			break;
		case CCV_64F:
			cblas_dgemm(CblasRowMajor, (transpose & CCV_A_TRANSPOSE) ? CblasTrans : CblasNoTrans, (transpose & CCV_B_TRANSPOSE) ? CblasTrans : CblasNoTrans, dd->rows, dd->cols, (transpose & CCV_A_TRANSPOSE) ? da->rows : da->cols, alpha, da->data.f64, da->cols, db->data.f64, db->cols, beta, dd->data.f64, dd->cols);
			break;
	}
#else
	assert(0 && "You need a BLAS compatible library for this function, e.g. libatlas.");
#endif
}

max_stars_repo_name: mrgloom/ccv

max_issues_repo_name: mrgloom/ccv

max_forks_repo_name: mrgloom/ccv


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

#define NUMBER_OF_ELEMENTS 1
#define MEAN 0
#define STANDARD_DEVIATION 1

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

        if (argc == 1) {
                number_of_elements = NUMBER_OF_ELEMENTS;
        } else if (argc == 2) {
                number_of_elements = atoi(argv[1]);
        } else {
                fprintf(stderr, "%s\n", "Wrong number of arguments.");
                fprintf(stderr, "%s\n", "Use: ./normal_numbers number_of_elements");
                exit(EXIT_FAILURE);
        }

        double *numbers = malloc (number_of_elements * sizeof(double));

        double mean = MEAN;
        double standard_deviation = STANDARD_DEVIATION;

        // declare the necessary random number generator variables
        const gsl_rng_type *T;
        gsl_rng *r;

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

        // create a seed based on time
        struct timeval tv;
        gettimeofday(&tv,0);
        unsigned long seed = tv.tv_sec + tv.tv_usec;

        // setup the generator using the seed just created
        T = gsl_rng_default;
        r = gsl_rng_alloc (T);
        gsl_rng_set(r, seed);

        // generate number_of_elements normal random numbers
        for (int i = 0; i < number_of_elements; i++) {
                numbers[i] = mean + gsl_ran_gaussian (r, standard_deviation);
        }

        // print numbers
        if(number_of_elements < 10) {
                for (int i = 0; i < number_of_elements; i++) {
                        printf (" %f", numbers[i]);
                }
                printf("\n");
        } else {
                printf("Done.\n");
        }

        // free memory
        gsl_rng_free (r);
        free(numbers);

        return 0;
}

max_stars_repo_name: sernamar/random-numbers

max_issues_repo_name: sernamar/random-numbers

max_forks_repo_name: sernamar/random-numbers


text:
/* Foundation and miscellenea for the statistics modules.
 *
 * Contents:
 *   1. Summary statistics (means, variances)
 *   2. Special functions
 *   3. Standard statistical tests
 *   4. Data fitting.
 *   5. Unit tests.
 *   6. Test driver.
 *   7. Examples.
 *      - driver for linear regression
 *      - driver for G-test
 *   8. License and copyright information.
 *
 */
#include "esl_config.h"

#include <math.h>

#include "easel.h"
#include "esl_stats.h"

/*****************************************************************
 * 1. Summary statistics calculations (means, variances)
 *****************************************************************/

/* Function:  esl_stats_DMean()
 * Synopsis:  Calculates mean and $\sigma^2$ for samples $x_i$.
 *
 * Purpose:   Calculates the sample mean and $s^2$, the unbiased
 *            estimator of the population variance, for a
 *            sample of <n> numbers <x[0]..x[n-1]>, and optionally
 *            returns either or both through <ret_mean> and
 *            <ret_var>.
 *
 *            <esl_stats_FMean()> and <esl_stats_IMean()> do the same,
 *            for float and integer vectors.
 *
 * Args:      x        - samples x[0]..x[n-1]
 *            n        - number of samples
 *            opt_mean - optRETURN: mean
 *            opt_var  - optRETURN: estimate of population variance
 *
 * Returns:   <eslOK> on success.
 */
int
esl_stats_DMean(const double *x, int n, double *opt_mean, double *opt_var)
{
  double sum   = 0.;
  double sqsum = 0.;
  int i;

  for (i = 0; i < n; i++)
    {
      sum   += x[i];
      sqsum += x[i]*x[i];
    }
  if (opt_mean != NULL)  *opt_mean = sum / (double) n;
  if (opt_var  != NULL)  *opt_var  = (sqsum - sum*sum/(double)n) / ((double)n-1);
  return eslOK;
}
int
esl_stats_FMean(const float *x, int n, double *opt_mean, double *opt_var)
{
  double sum   = 0.;
  double sqsum = 0.;
  int i;

  for (i = 0; i < n; i++)
    {
      sum   += x[i];
      sqsum += x[i]*x[i];
    }
  if (opt_mean != NULL)  *opt_mean = sum / (double) n;
  if (opt_var  != NULL)  *opt_var  = (sqsum - sum*sum/(double)n) / ((double)n-1);
  return eslOK;
}
int
esl_stats_IMean(const int *x, int n, double *opt_mean, double *opt_var)
{
  double sum   = 0.;
  double sqsum = 0.;
  int i;

  for (i = 0; i < n; i++)
    {
      sum   += x[i];
      sqsum += x[i]*x[i];
    }
  if (opt_mean != NULL)  *opt_mean = sum / (double) n;
  if (opt_var  != NULL)  *opt_var  = (sqsum - sum*sum/(double)n) / ((double)n-1);
  return eslOK;
}
/*--------------- end, summary statistics -----------------------*/

/*****************************************************************
 * 2. Special functions.
 *****************************************************************/

/* Function:  esl_stats_LogGamma()
 * Synopsis:  Calculates $\log \Gamma(x)$.
 *
 * Purpose:   Returns natural log of $\Gamma(x)$, for $x > 0$.
 *
 * Credit:    Adapted from a public domain implementation in the
 *            NCBI core math library. Thanks to John Spouge and
 *            the NCBI. (According to NCBI, that's Dr. John
 *            "Gammas Galore" Spouge to you, pal.)
 *
 * Args:      x          : argument, x > 0.0
 *            ret_answer : RETURN: the answer
 *
 * Returns:   Put the answer in <ret_answer>; returns <eslOK>.
 *
 * Throws:    <eslERANGE> if $x <= 0$.
 */
int
esl_stats_LogGamma(double x, double *ret_answer)
{
  int i;
  double xx, tx;
  double tmp, value;
  static double cof[11] = {
     4.694580336184385e+04,
    -1.560605207784446e+05,
     2.065049568014106e+05,
    -1.388934775095388e+05,
     5.031796415085709e+04,
    -9.601592329182778e+03,
     8.785855930895250e+02,
    -3.155153906098611e+01,
     2.908143421162229e-01,
    -2.319827630494973e-04,
     1.251639670050933e-10
  };

  /* Protect against invalid x<=0  */
  if (x <= 0.0)  ESL_EXCEPTION(eslERANGE, "invalid x <= 0 in esl_stats_LogGamma()");

  xx       = x - 1.0;
  tx = tmp = xx + 11.0;
  value    = 1.0;
  for (i = 10; i >= 0; i--)	/* sum least significant terms first */
    {
      value += cof[i] / tmp;
      tmp   -= 1.0;
    }
  value  = log(value);
  tx    += 0.5;
  value += 0.918938533 + (xx+0.5)*log(tx) - tx;
  *ret_answer = value;
  return eslOK;
}

/* Function:  esl_stats_Psi()
 * Synopsis:  Calculates $\Psi(x)$ (the digamma function).
 *
 * Purpose:   Computes $\Psi(x)$ (the "digamma" function), which is
 *            the derivative of log of the Gamma function:
 *            $d/dx \log \Gamma(x) = \frac{\Gamma'(x)}{\Gamma(x)} = \Psi(x)$.
 *            Argument $x$ is $> 0$.
 *
 *            This is J.M. Bernardo's "Algorithm AS103",
 *            Appl. Stat. 25:315-317 (1976).
 */
int
esl_stats_Psi(double x, double *ret_answer)
{
  double answer = 0.;
  double x2;

  if (x <= 0.0) ESL_EXCEPTION(eslERANGE, "invalid x <= 0 in esl_stats_Psi()");

  /* For small x, Psi(x) ~= -0.5772 - 1/x + O(x), we're done.
   */
  if (x <= 1e-5) {
    *ret_answer = -eslCONST_EULER - 1./x;
    return eslOK;
  }

  /* For medium x, use Psi(1+x) = \Psi(x) + 1/x to c.o.v. x,
   * big enough for Stirling approximation to work...
   */
  while (x < 8.5) {
    answer = answer - 1./x;
    x += 1.;
  }

  /* For large X, use Stirling approximation
   */
  x2 = 1./x;
  answer += log(x) - 0.5 * x2;
  x2 = x2*x2;
  answer -= (1./12.)*x2;
  answer += (1./120.)*x2*x2;
  answer -= (1./252.)*x2*x2*x2;

  *ret_answer = answer;
  return eslOK;
}

/* Function: esl_stats_IncompleteGamma()
 * Synopsis: Calculates the incomplete Gamma function.
 *
 * Purpose:  Returns $P(a,x)$ and $Q(a,x)$ where:
 *
 *           \begin{eqnarray*}
 *             P(a,x) & = & \frac{1}{\Gamma(a)} \int_{0}^{x} t^{a-1} e^{-t} dt \\
 *                    & = & \frac{\gamma(a,x)}{\Gamma(a)} \\
 *             Q(a,x) & = & \frac{1}{\Gamma(a)} \int_{x}^{\infty} t^{a-1} e^{-t} dt\\
 *                    & = & 1 - P(a,x) \\
 *           \end{eqnarray*}
 *
 *           $P(a,x)$ is the CDF of a gamma density with $\lambda = 1$,
 *           and $Q(a,x)$ is the survival function.
 *
 *           For $x \simeq 0$, $P(a,x) \simeq 0$ and $Q(a,x) \simeq 1$; and
 *           $P(a,x)$ is less prone to roundoff error.
 *
 *           The opposite is the case for large $x >> a$, where
 *           $P(a,x) \simeq 1$ and $Q(a,x) \simeq 0$; there, $Q(a,x)$ is
 *           less prone to roundoff error.
 *
 * Method:   Based on ideas from Numerical Recipes in C, Press et al.,
 *           Cambridge University Press, 1988.
 *
 * Args:     a          - for instance, degrees of freedom / 2     [a > 0]
 *           x          - for instance, chi-squared statistic / 2  [x >= 0]
 *           ret_pax    - RETURN: P(a,x)
 *           ret_qax    - RETURN: Q(a,x)
 *
 * Return:   <eslOK> on success.
 *
 * Throws:   <eslERANGE> if <a> or <x> is out of accepted range.
 *           <eslENOHALT> if approximation fails to converge.
 */
int
esl_stats_IncompleteGamma(double a, double x, double *ret_pax, double *ret_qax)
{
  int    iter;			/* iteration counter */
  double pax;			/* P(a,x) */
  double qax;			/* Q(a,x) */

  if (a <= 0.) ESL_EXCEPTION(eslERANGE, "esl_stats_IncompleteGamma(): a must be > 0");
  if (x <  0.) ESL_EXCEPTION(eslERANGE, "esl_stats_IncompleteGamma(): x must be >= 0");

  /* For x > a + 1 the following gives rapid convergence;
   * calculate Q(a,x) = \frac{\Gamma(a,x)}{\Gamma(a)},
   * using a continued fraction development for \Gamma(a,x).
   */
  if (x > a+1)
    {
      double oldp;		/* previous value of p    */
      double nu0, nu1;		/* numerators for continued fraction calc   */
      double de0, de1;		/* denominators for continued fraction calc */

      nu0 = 0.;			/* A_0 = 0       */
      de0 = 1.;			/* B_0 = 1       */
      nu1 = 1.;			/* A_1 = 1       */
      de1 = x;			/* B_1 = x       */

      oldp = nu1;
      for (iter = 1; iter < 100; iter++)
	{
	  /* Continued fraction development:
	   * set A_j = b_j A_j-1 + a_j A_j-2
	   *     B_j = b_j B_j-1 + a_j B_j-2
           * We start with A_2, B_2.
	   */
				/* j = even: a_j = iter-a, b_j = 1 */
				/* A,B_j-2 are in nu0, de0; A,B_j-1 are in nu1,de1 */
	  nu0 = nu1 + ((double)iter - a) * nu0;
	  de0 = de1 + ((double)iter - a) * de0;
				/* j = odd: a_j = iter, b_j = x */
				/* A,B_j-2 are in nu1, de1; A,B_j-1 in nu0,de0 */
	  nu1 = x * nu0 + (double) iter * nu1;
	  de1 = x * de0 + (double) iter * de1;
				/* rescale */
	  if (de1 != 0.)
	    {
	      nu0 /= de1;
	      de0 /= de1;
	      nu1 /= de1;
	      de1 =  1.;
	    }
				/* check for convergence */
	  if (fabs((nu1-oldp)/nu1) < 1.e-7)
	    {
	      esl_stats_LogGamma(a, &qax);
	      qax = nu1 * exp(a * log(x) - x - qax);

	      if (ret_pax != NULL) *ret_pax = 1 - qax;
	      if (ret_qax != NULL) *ret_qax = qax;
	      return eslOK;
	    }

	  oldp = nu1;
	}
      ESL_EXCEPTION(eslENOHALT,
		"esl_stats_IncompleteGamma(): fraction failed to converge");
    }
  else /* x <= a+1 */
    {
      double p;			/* current sum               */
      double val;		/* current value used in sum */

      /* For x <= a+1 we use a convergent series instead:
       *   P(a,x) = \frac{\gamma(a,x)}{\Gamma(a)},
       * where
       *   \gamma(a,x) = e^{-x}x^a \sum_{n=0}{\infty} \frac{\Gamma{a}}{\Gamma{a+1+n}} x^n
       * which looks appalling but the sum is in fact rearrangeable to
       * a simple series without the \Gamma functions:
       *   = \frac{1}{a} + \frac{x}{a(a+1)} + \frac{x^2}{a(a+1)(a+2)} ...
       * and it's obvious that this should converge nicely for x <= a+1.
       */
      p = val = 1. / a;
      for (iter = 1; iter < 10000; iter++)
	{
	  val *= x / (a+(double)iter);
	  p   += val;

	  if (fabs(val/p) < 1.e-7)
	    {
	      esl_stats_LogGamma(a, &pax);
	      pax = p * exp(a * log(x) - x - pax);

	      if (ret_pax != NULL) *ret_pax = pax;
	      if (ret_qax != NULL) *ret_qax = 1. - pax;
	      return eslOK;
	    }
	}
      ESL_EXCEPTION(eslENOHALT,
		"esl_stats_IncompleteGamma(): series failed to converge");
    }
  /*NOTREACHED*/
  return eslOK;
}
/*----------------- end, special functions ----------------------*/

/*****************************************************************
 * 3. Standard statistical tests.
 *****************************************************************/

/* Function:  esl_stats_GTest()
 * Synopsis:  Calculates a G-test on 2 vs. 1 binomials.
 *
 * Purpose:   In experiment a, we've drawn <ca> successes in <na> total
 *            trials; in experiment b, we've drawn <cb> successes in
 *            <nb> total trials. Are the counts different enough to
 *            conclude that the two experiments are different? The
 *            null hypothesis is that the successes in both experiments
 *            were drawn from the same binomial distribution with
 *            per-trial probability $p$. The tested hypothesis is that
 *            experiments a,b have different binomial probabilities
 *            $p_a,p_b$. The G-test is a log-likelihood-ratio statistic,
 *            assuming maximum likelihood values for $p,p_a,p_b$.
 *            $2G$ is distributed approximately as $X^2(1)$,
 *            %"X" is "Chi"
 *            which we use to calculate a P-value for the G statistic.
 *
 * Args:      ca    - number of positives in experiment a
 *            na    - total number in experiment a
 *            cb    - number of positives in experiment b
 *            nb    - total number in experiment b
 *            ret_G - RETURN: G statistic, a log likelihood ratio, in nats
 *            ret_P - RETURN: P-value for the G-statistic
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    (no abnormal error conditions)
 *
 * Xref:      Archive1999/0906-sagescore/sagescore.c
 */
int
esl_stats_GTest(int ca, int na, int cb, int nb, double *ret_G, double *ret_P)
{
  double a,b,c,d,n;
  double G = 0.;

  a = (double) ca;
  b = (double) (na - ca);
  c = (double) cb;
  d = (double) (nb - cb);
  n = (double) na+nb;

  /* Yes, the calculation here is correct; algebraic
   * rearrangement of the log-likelihood-ratio with
   * p_a = ca/na, p_b = cb/nb, and p = (ca+cb)/(na+nb).
   * Guard against 0 probabilities; assume 0 log 0 => 0.
   */
  if (a   > 0.) G  = a * log(a);
  if (b   > 0.) G += b * log(b);
  if (c   > 0.) G += c * log(c);
  if (d   > 0.) G += d * log(d);
  if (n   > 0.) G += n * log(n);
  if (a+b > 0.) G -= (a+b) * log(a+b);
  if (c+d > 0.) G -= (c+d) * log(c+d);
  if (a+c > 0.) G -= (a+c) * log(a+c);
  if (b+d > 0.) G -= (b+d) * log(b+d);

  *ret_G = G;
  return esl_stats_IncompleteGamma( 0.5, G, NULL, ret_P);
}

/* Function:  esl_stats_ChiSquaredTest()
 * Synopsis:  Calculates a $\chi^2$ P-value.
 * Incept:    SRE, Tue Jul 19 11:39:32 2005 [St. Louis]
 *
 * Purpose:   Calculate the probability that a chi-squared statistic
 *            with <v> degrees of freedom would exceed the observed
 *            chi-squared value <x>; return it in <ret_answer>. If
 *            this probability is less than some small threshold (say,
 *            0.05 or 0.01), then we may reject the hypothesis we're
 *            testing.
 *
 * Args:      v          - degrees of freedom
 *            x          - observed chi-squared value
 *            ret_answer - RETURN: P(\chi^2 > x)
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslERANGE> if <v> or <x> are out of valid range.
 *            <eslENOHALT> if iterative calculation fails.
 */
int
esl_stats_ChiSquaredTest(int v, double x, double *ret_answer)
{
  return esl_stats_IncompleteGamma((double)v/2, x/2, NULL, ret_answer);
}
/*----------------- end, statistical tests  ---------------------*/

/*****************************************************************
 * 4. Data fitting.
 *****************************************************************/

/* Function:  esl_stats_LinearRegression()
 * Synopsis:  Fit data to a straight line.
 * Incept:    SRE, Sat May 26 11:33:46 2007 [Janelia]
 *
 * Purpose:   Fit <n> points <x[i]>, <y[i]> to a straight line
 *            $y = a + bx$ by linear regression.
 *
 *            The $x_i$ are taken to be known, and the $y_i$ are taken
 *            to be observed quantities associated with a sampling
 *            error $\sigma_i$. If known, the standard deviations
 *            $\sigma_i$ for $y_i$ are provided in the <sigma> array.
 *            If they are unknown, pass <sigma = NULL>, and the
 *            routine will proceed with the assumption that $\sigma_i
 *            = 1$ for all $i$.
 *
 *            The maximum likelihood estimates for $a$ and $b$ are
 *            optionally returned in <opt_a> and <opt_b>.
 *
 *            The estimated standard deviations of $a$ and $b$ and
 *            their estimated covariance are optionally returned in
 *            <opt_sigma_a>, <opt_sigma_b>, and <opt_cov_ab>.
 *
 *            The Pearson correlation coefficient is optionally
 *            returned in <opt_cc>.
 *
 *            The $\chi^2$ P-value for the regression fit is
 *            optionally returned in <opt_Q>. This P-value may only be
 *            obtained when the $\sigma_i$ are known. If <sigma> is
 *            passed as <NULL> and <opt_Q> is requested, <*opt_Q> is
 *            set to 1.0.
 *
 *            This routine follows the description and algorithm in
 *            \citep[pp.661-666]{Press93}.
 *
 *            <n> must be greater than 2; at least two x[i] must
 *            differ; and if <sigma> is provided, all <sigma[i]> must
 *            be $>0$. If any of these conditions isn't met, the
 *            routine throws <eslEINVAL>.
 *
 * Args:      x            - x[0..n-1]
 *            y            - y[0..n-1]
 *            sigma        - sample error in observed y_i
 *            n            - number of data points
 *            opt_a        - optRETURN: intercept estimate
 *            opt_b        - optRETURN: slope estimate
 *            opt_sigma_a  - optRETURN: error in estimate of a
 *            opt_sigma_b  - optRETURN: error in estimate of b
 *            opt_cov_ab   - optRETURN: covariance of a,b estimates
 *            opt_cc       - optRETURN: Pearson correlation coefficient for x,y
 *            opt_Q        - optRETURN: X^2 P-value for linear fit
 *
 * Returns:   <eslOK> on success.
 *
 * Throws:    <eslEMEM> on allocation error;
 *            <eslEINVAL> if a contract condition isn't met;
 *            <eslENORESULT> if the chi-squared test fails.
 *            In these cases, all optional return values are set to 0.
 */
int
esl_stats_LinearRegression(const double *x, const double *y, const double *sigma, int n,
			   double *opt_a,       double *opt_b,
			   double *opt_sigma_a, double *opt_sigma_b, double *opt_cov_ab,
			   double *opt_cc,      double *opt_Q)
{
  int     status;
  double *t      = NULL;
  double  S, Sx, Sy, Stt;
  double  Sxy, Sxx, Syy;
  double  a, b, sigma_a, sigma_b, cov_ab, cc, X2, Q;
  double  xdev, ydev;
  double  tmp;
  int     i;

  /* Contract checks. */
  if (n <= 2) ESL_XEXCEPTION(eslEINVAL, "n must be > 2 for linear regression fitting");
  if (sigma != NULL)
    for (i = 0; i < n; i++) if (sigma[i] <= 0.) ESL_XEXCEPTION(eslEINVAL, "sigma[%d] <= 0", i);
  status = eslEINVAL;
  for (i = 0; i < n; i++) if (x[i] != 0.) { status = eslOK; break; }
  if (status != eslOK) ESL_XEXCEPTION(eslEINVAL, "all x[i] are 0.");

  /* Allocations */
  ESL_ALLOC(t, sizeof(double) * n);

  /* S = \sum_{i=1}{n} \frac{1}{\sigma_i^2}.  (S > 0.) */
  if (sigma != NULL) { for (S = 0., i = 0; i < n; i++) S += 1./ (sigma[i] * sigma[i]);  }
  else S = (double) n;

  /* S_x = \sum_{i=1}{n} \frac{x[i]}{ \sigma_i^2}  (Sx real.) */
  for (Sx = 0., i = 0; i < n; i++) {
    if (sigma == NULL) Sx += x[i];
    else               Sx += x[i] / (sigma[i] * sigma[i]);
  }

  /* S_y = \sum_{i=1}{n} \frac{y[i]}{\sigma_i^2}  (Sy real.) */
  for (Sy = 0., i = 0; i < n; i++) {
    if (sigma == NULL) Sy += y[i];
    else               Sy += y[i] / (sigma[i] * sigma[i]);
  }

  /* t_i = \frac{1}{\sigma_i} \left( x_i - \frac{S_x}{S} \right)   (t_i real) */
  for (i = 0; i < n; i++) {
    t[i] = x[i] - Sx/S;
    if (sigma != NULL) t[i] /= sigma[i];
  }

  /* S_{tt} = \sum_{i=1}^n t_i^2  (if at least one x is != 0, Stt > 0) */
  for (Stt = 0., i = 0; i < n; i++) { Stt += t[i] * t[i]; }

  /* b = \frac{1}{S_{tt}} \sum_{i=1}^{N} \frac{t_i y_i}{\sigma_i}  */
  for (b = 0., i = 0; i < n; i++) {
    if (sigma != NULL) { b += t[i]*y[i] / sigma[i]; }
    else               { b += t[i]*y[i]; }
  }
  b /= Stt;

  /* a = \frac{ S_y - S_x b } {S}   */
  a = (Sy - Sx * b) / S;

  /* \sigma_a^2 = \frac{1}{S} \left( 1 + \frac{ S_x^2 }{S S_{tt}} \right) */
  sigma_a = sqrt ((1. + (Sx*Sx) / (S*Stt)) / S);

  /* \sigma_b = \frac{1}{S_{tt}} */
  sigma_b = sqrt (1. / Stt);

  /* Cov(a,b) = - \frac{S_x}{S S_{tt}}    */
  cov_ab = -Sx / (S * Stt);

  /* Pearson correlation coefficient */
  Sxy = Sxx = Syy = 0.;
  for (i = 0; i < n; i++) {
    if (sigma != NULL) {
      xdev = (x[i] / (sigma[i] * sigma[i])) - (Sx / n);
      ydev = (y[i] / (sigma[i] * sigma[i])) - (Sy / n);
    } else {
      xdev = x[i] - (Sx / n);
      ydev = y[i] - (Sy / n);
    }
    Sxy += xdev * ydev;
    Sxx += xdev * xdev;
    Syy += ydev * ydev;
  }
  cc = Sxy / (sqrt(Sxx) * sqrt(Syy));

  /* \chi^2 */
  for (X2 = 0., i = 0; i < n; i++) {
    tmp =  y[i] - a - b*x[i];
    if (sigma != NULL) tmp /= sigma[i];
    X2 += tmp*tmp;
  }

  /* We can calculate a goodness of fit if we know the \sigma_i */
  if (sigma != NULL) {
    if (esl_stats_ChiSquaredTest(n-2, X2, &Q) != eslOK) { status = eslENORESULT; goto ERROR; }
  } else Q = 1.0;

  /* If we didn't use \sigma_i, adjust the sigmas for a,b */
  if (sigma == NULL) {
    tmp = sqrt(X2 / (double)(n-2));
    sigma_a *= tmp;
    sigma_b *= tmp;
  }

  /* Done. Set up for normal return.
   */
  free(t);
  if (opt_a       != NULL) *opt_a       = a;
  if (opt_b       != NULL) *opt_b       = b;
  if (opt_sigma_a != NULL) *opt_sigma_a = sigma_a;
  if (opt_sigma_b != NULL) *opt_sigma_b = sigma_b;
  if (opt_cov_ab  != NULL) *opt_cov_ab  = cov_ab;
  if (opt_cc      != NULL) *opt_cc      = cc;
  if (opt_Q       != NULL) *opt_Q       = Q;
  return eslOK;

 ERROR:
  if (t != NULL) free(t);
  if (opt_a       != NULL) *opt_a       = 0.;
  if (opt_b       != NULL) *opt_b       = 0.;
  if (opt_sigma_a != NULL) *opt_sigma_a = 0.;
  if (opt_sigma_b != NULL) *opt_sigma_b = 0.;
  if (opt_cov_ab  != NULL) *opt_cov_ab  = 0.;
  if (opt_cc      != NULL) *opt_cc      = 0.;
  if (opt_Q       != NULL) *opt_Q       = 0.;
  return status;
}
/*------------------- end, data fitting -------------------------*/

/*****************************************************************
 * 5. Unit tests.
 *****************************************************************/
#ifdef eslSTATS_TESTDRIVE
#include "esl_random.h"
#include "esl_stopwatch.h"
#ifdef HAVE_LIBGSL
#include <gsl/gsl_sf_gamma.h>
#endif

/* The LogGamma() function is rate-limiting in hmmbuild, because it is
 * used so heavily in mixture Dirichlet calculations.
 *    ./configure --with-gsl; [compile test driver]
 *    ./stats_utest -v
 * runs a comparison of time/precision against GSL.
 * SRE, Sat May 23 10:04:41 2009, on home Mac:
 *     LogGamma       = 1.29u  / N=1e8  =  13 nsec/call
 *     gsl_sf_lngamma = 1.43u  / N=1e8  =  14 nsec/call
 */
static void
utest_LogGamma(ESL_RANDOMNESS *r, int N, int be_verbose)
{
  char          *msg = "esl_stats_LogGamma() unit test failed";
  ESL_STOPWATCH *w   = esl_stopwatch_Create();
  double        *x   = malloc(sizeof(double) * N);
  double        *lg  = malloc(sizeof(double) * N);
  double        *lg2 = malloc(sizeof(double) * N);
  int            i;

  for (i = 0; i < N; i++)
    x[i] = esl_random(r) * 100.;

  esl_stopwatch_Start(w);
  for (i = 0; i < N; i++)
    if (esl_stats_LogGamma(x[i], &(lg[i])) != eslOK) esl_fatal(msg);
  esl_stopwatch_Stop(w);

  if (be_verbose) esl_stopwatch_Display(stdout, w, "esl_stats_LogGamma() timing: ");

#ifdef HAVE_LIBGSL
  esl_stopwatch_Start(w);
  for (i = 0; i < N; i++) lg2[i] = gsl_sf_lngamma(x[i]);
  esl_stopwatch_Stop(w);

  if (be_verbose) esl_stopwatch_Display(stdout, w, "gsl_sf_lngamma() timing:     ");

  for (i = 0; i < N; i++)
    if (esl_DCompare(lg[i], lg2[i], 1e-2) != eslOK) esl_fatal(msg);
#endif

  free(lg2);
  free(lg);
  free(x);
  esl_stopwatch_Destroy(w);
}

/* The test of esl_stats_LinearRegression() is a statistical test,
 * so we can't be too aggressive about testing results.
 *
 * Args:
 *    r          - a source of randomness
 *    use_sigma  - TRUE to pass sigma to the regression fit.
 *    be_verbose - TRUE to print results (manual, not automated test mode)
 */
static void
utest_LinearRegression(ESL_RANDOMNESS *r, int use_sigma, int be_verbose)
{
  char msg[] = "linear regression unit test failed";
  double a     = -3.;
  double b     = 1.;
  int    n     = 100;
  double xori  = -20.;
  double xstep = 1.0;
  double setsigma = 1.0;		/* sigma on all points */
  int    i;
  double *x     = NULL;
  double *y     = NULL;
  double *sigma = NULL;
  double  ae, be, siga, sigb, cov_ab, cc, Q;

  if ((x     = malloc(sizeof(double) * n)) == NULL) esl_fatal(msg);
  if ((y     = malloc(sizeof(double) * n)) == NULL) esl_fatal(msg);
  if ((sigma = malloc(sizeof(double) * n)) == NULL) esl_fatal(msg);

  /* Simulate some linear data */
  for (i = 0; i < n; i++)
    {
      sigma[i] = setsigma;
      x[i]     = xori + i*xstep;
      y[i]     = esl_rnd_Gaussian(r, a + b*x[i], sigma[i]);
    }

  if (use_sigma) {
    if (esl_stats_LinearRegression(x, y, sigma, n, &ae, &be, &siga, &sigb, &cov_ab, &cc, &Q) != eslOK) esl_fatal(msg);
  } else {
    if (esl_stats_LinearRegression(x, y,  NULL, n, &ae, &be, &siga, &sigb, &cov_ab, &cc, &Q) != eslOK) esl_fatal(msg);
  }

  if (be_verbose) {
    printf("Linear regression test:\n");
    printf("estimated intercept a = %8.4f   [true = %8.4f]\n", ae, a);
    printf("estimated slope b     = %8.4f   [true = %8.4f]\n", be, b);
    printf("estimated sigma on a  = %8.4f\n",                  siga);
    printf("estimated sigma on b  = %8.4f\n",                  sigb);
    printf("estimated cov(a,b)    = %8.4f\n",                  cov_ab);
    printf("correlation coeff     = %8.4f\n",                  cc);
    printf("P-value               = %8.4f\n",                  Q);
  }

  /* The following tests are statistical.
   */
  if ( fabs(ae-a) > 2*siga ) esl_fatal(msg);
  if ( fabs(be-b) > 2*sigb ) esl_fatal(msg);
  if ( cc < 0.95)            esl_fatal(msg);
  if (use_sigma) {
    if (Q < 0.001)           esl_fatal(msg);
  } else {
    if (Q != 1.0)            esl_fatal(msg);
  }

  free(x);
  free(y);
  free(sigma);
}
#endif /*eslSTATS_TESTDRIVE*/
/*-------------------- end of unit tests ------------------------*/

/*****************************************************************
 * 6. Test driver.
 *****************************************************************/
#ifdef eslSTATS_TESTDRIVE
/* gcc -g -Wall -o stats_utest  -L. -I. -DeslSTATS_TESTDRIVE esl_stats.c -leasel -lm
 * gcc -DHAVE_LIBGSL -O2 -o stats_utest -L. -I. -DeslSTATS_TESTDRIVE esl_stats.c -leasel -lgsl -lm
 */
#include <stdio.h>
#include "easel.h"
#include "esl_getopts.h"
#include "esl_random.h"
#include "esl_stats.h"

static ESL_OPTIONS options[] = {
  /* name  type         default  env   range togs  reqs  incomp  help                docgrp */
  {"-h",  eslARG_NONE,    FALSE, NULL, NULL, NULL, NULL, NULL, "show help and usage",                   0},
  {"-s",  eslARG_INT,      "42", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>",         0},
  {"-v",  eslARG_NONE,    FALSE, NULL, NULL, NULL, NULL, NULL, "verbose: show verbose output",          0},
  {"-N",  eslARG_INT,"10000000", NULL, NULL, NULL, NULL, NULL, "number of trials in LogGamma test",     0},
  { 0,0,0,0,0,0,0,0,0,0},
};
static char usage[]  = "[-options]";
static char banner[] = "test driver for stats special functions";

int
main(int argc, char **argv)
{
  ESL_GETOPTS    *go         = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
  ESL_RANDOMNESS *r          = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
  int             be_verbose = esl_opt_GetBoolean(go, "-v");
  int             N          = esl_opt_GetInteger(go, "-N");

  if (be_verbose) printf("seed = %" PRIu32 "\n", esl_randomness_GetSeed(r));

  utest_LogGamma(r, N, be_verbose);
  utest_LinearRegression(r, TRUE,  be_verbose);
  utest_LinearRegression(r, FALSE, be_verbose);

  esl_getopts_Destroy(go);
  esl_randomness_Destroy(r);
  exit(0);
}
#endif /*eslSTATS_TESTDRIVE*/
/*------------------- end of test driver ------------------------*/

/*****************************************************************
 * 7. Examples.
 *****************************************************************/

/* Compile:  gcc -g -Wall -o example -I. -DeslSTATS_EXAMPLE esl_stats.c esl_random.c easel.c -lm
 * or        gcc -g -Wall -o example -I. -L. -DeslSTATS_EXAMPLE esl_stats.c -leasel -lm
 */
#ifdef eslSTATS_EXAMPLE
/*::cexcerpt::stats_example::begin::*/
/* gcc -g -Wall -o example -I. -DeslSTATS_EXAMPLE esl_stats.c esl_random.c easel.c -lm  */
#include <stdio.h>
#include "easel.h"
#include "esl_random.h"
#include "esl_stats.h"

int main(void)
{
  ESL_RANDOMNESS *r   = esl_randomness_Create(0);
  double a            = -3.;
  double b            = 1.;
  double xori         = -20.;
  double xstep        = 1.0;
  double setsigma     = 1.0;		/* sigma on all points */
  int    n            = 100;
  double *x           = malloc(sizeof(double) * n);
  double *y           = malloc(sizeof(double) * n);
  double *sigma       = malloc(sizeof(double) * n);
  int    i;
  double  ae, be, siga, sigb, cov_ab, cc, Q;

  /* Simulate some linear data, with Gaussian noise added to y_i */
  for (i = 0; i < n; i++) {
    sigma[i] = setsigma;
    x[i]     = xori + i*xstep;
    y[i]     = esl_rnd_Gaussian(r, a + b*x[i], sigma[i]);
  }

  if (esl_stats_LinearRegression(x, y, sigma, n, &ae, &be, &siga, &sigb, &cov_ab, &cc, &Q) != eslOK)
    esl_fatal("linear regression failed");

  printf("estimated intercept a = %8.4f   [true = %8.4f]\n", ae, a);
  printf("estimated slope b     = %8.4f   [true = %8.4f]\n", be, b);
  printf("estimated sigma on a  = %8.4f\n",                  siga);
  printf("estimated sigma on b  = %8.4f\n",                  sigb);
  printf("estimated cov(a,b)    = %8.4f\n",                  cov_ab);
  printf("correlation coeff     = %8.4f\n",                  cc);
  printf("P-value               = %8.4f\n",                  Q);

  free(x);  free(y);  free(sigma);
  esl_randomness_Destroy(r);
  exit(0);
}
/*::cexcerpt::stats_example::end::*/
#endif /* eslSTATS_EXAMPLE */

#ifdef eslSTATS_EXAMPLE2

#include <stdlib.h>

#include "easel.h"
#include "esl_getopts.h"
#include "esl_stats.h"

static ESL_OPTIONS options[] = {
  /* name           type      default  env  range toggles reqs incomp  help                                       docgroup*/
  { "-h",        eslARG_NONE,   FALSE,  NULL, NULL,  NULL,  NULL, NULL, "show brief help on version and usage",    0 },
  {  0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
};
static char usage[]  = "[-options] <ca> <na> <cb> <nb>";
static char banner[] = "example from the stats module: using a G-test";

int
main(int argc, char **argv)
{
  ESL_GETOPTS  *go  = esl_getopts_CreateDefaultApp(options, 4, argc, argv, banner, usage);
  int           ca  = strtol(esl_opt_GetArg(go, 1), NULL, 10);
  int           na  = strtol(esl_opt_GetArg(go, 2), NULL, 10);
  int           cb  = strtol(esl_opt_GetArg(go, 3), NULL, 10);
  int           nb  = strtol(esl_opt_GetArg(go, 4), NULL, 10);
  double        G, P;
  int           status;

  if (ca > na || cb > nb) esl_fatal("argument order wrong? expect ca, na, cb, nb for ca/na, cb/nb");

  if ( (status = esl_stats_GTest(ca, na, cb, nb, &G, &P)) != eslOK) esl_fatal("G-test failed?");
  printf("%-10.3g %12.2f\n", P, G);
  exit(0);
}
#endif /* eslSTATS_EXAMPLE2 */
/*--------------------- end of examples -------------------------*/

/*****************************************************************
 * Easel - a library of C functions for biological sequence analysis
 * Version h3.1b2; February 2015
 * Copyright (C) 2015 Howard Hughes Medical Institute.
 * Other copyrights also apply. See the COPYRIGHT file for a full list.
 *
 * Easel is distributed under the Janelia Farm Software License, a BSD
 * license. See the LICENSE file for more details.
 *
 * SVN $Id: esl_stats.c 854 2013-02-25 22:00:19Z wheelert $
 * SVN $URL: https://svn.janelia.org/eddylab/eddys/easel/branches/hmmer/3.1/esl_stats.c $
 *****************************************************************/