text:
/**
 *
 * @file timing.c
 *
 *  PLASMA auxiliary routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author ???
 * @author Mathieu Faverge
 * @author Dulceneia Becker
 * @date 2010-11-15
 *
 **/

#if defined( _WIN32 ) || defined( _WIN64 )
#define int64_t __int64
#endif

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifdef PLASMA_EZTRACE
#include <eztrace.h>
#endif

#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

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

#include <cblas.h>
#include <lapacke.h>
#include <plasma.h>
#include <core_blas.h>
#include "flops.h"
#include "timing.h"
#include "auxiliary.h"
#include <omp.h>

static int RunTest(int *iparam, _PREC *dparam, double *t_);

double cWtime(void);

int ISEED[4] = {0,0,0,1};   /* initial seed for zlarnv() */

/*
 * struct timeval {time_t tv_sec; suseconds_t tv_usec;};
 */
double cWtime(void)
{
    struct timeval tp;
    gettimeofday( &tp, NULL );
    return tp.tv_sec + 1e-6 * tp.tv_usec;
}

static int
Test(int64_t n, int *iparam) {
    int      i, j, iter;
    int      thrdnbr, niter;
    int64_t  M, N, K, NRHS;
    double  *t;
    _PREC    eps = _LAMCH( 'e' );
    _PREC    dparam[IPARAM_DNBPARAM];
    double   fmuls, fadds, fp_per_mul, fp_per_add;
    double   sumgf, sumgf2, sumt, sd, flops, gflops;
    char    *s;
    char    *env[] = {
        "OMP_NUM_THREADS",
        "MKL_NUM_THREADS",
        "GOTO_NUM_THREADS",
        "ACML_NUM_THREADS",
        "ATLAS_NUM_THREADS",
        "BLAS_NUM_THREADS", ""
    };
    int gnuplot = 0;

    memset( &dparam, 0, IPARAM_DNBPARAM * sizeof(_PREC) );

    thrdnbr = iparam[IPARAM_THRDNBR];
    niter   = iparam[IPARAM_NITER];

    M    = iparam[IPARAM_M];
    N    = iparam[IPARAM_N];
    K    = iparam[IPARAM_K];
    NRHS = K;
    (void)M;(void)N;(void)K;(void)NRHS;

    if ( (n < 0) || (thrdnbr < 0 ) ) {
        if (gnuplot) {
            printf( "set title '%d_NUM_THREADS: ", thrdnbr );
            for (i = 0; env[i][0]; ++i) {
                s = getenv( env[i] );

                if (i) printf( " " ); /* separating space */

                for (j = 0; j < 5 && env[i][j] && env[i][j] != '_'; ++j)
                    printf( "%c", env[i][j] );

                if (s)
                    printf( "=%s", s );
                else
                    printf( "->%s", "?" );
            }
            printf( "'\n" );
            printf( "%s\n%s\n%s\n%s\n%s%s%s\n",
                    "set xlabel 'Matrix size'",
                    "set ylabel 'Gflop/s'",
                    "set key bottom",
                    gnuplot > 1 ? "set terminal png giant\nset output 'timeplot.png'" : "",
                    "plot '-' using 1:5 title '", _NAME, "' with linespoints" );
        }

        return 0;
    }

    printf( "%7d %7d %7d ", iparam[IPARAM_M], iparam[IPARAM_N], iparam[IPARAM_K] );
    fflush( stdout );

    t = (double*)malloc(niter*sizeof(double));
    memset(t, 0, niter*sizeof(double));

    if (sizeof(_TYPE) == sizeof(_PREC)) {
        fp_per_mul = 1;
        fp_per_add = 1;
    } else {
        fp_per_mul = 6;
        fp_per_add = 2;
    }

    fadds = (double)(_FADDS);
    fmuls = (double)(_FMULS);
    flops = 1e-9 * (fmuls * fp_per_mul + fadds * fp_per_add);
    gflops = 0.0;

    if ( iparam[IPARAM_WARMUP] ) {
        RunTest( iparam, dparam, &(t[0]));
    }

    sumgf  = 0.0;
    sumgf2 = 0.0;
    sumt   = 0.0;

    for (iter = 0; iter < niter; iter++)
    {
        if( iter == 0 ) {
          if ( iparam[IPARAM_TRACE] )
            iparam[IPARAM_TRACE] = 2;
          if ( iparam[IPARAM_DAG] )
            iparam[IPARAM_DAG] = 2;

          RunTest( iparam, dparam, &(t[iter]));

          iparam[IPARAM_TRACE] = 0;
          iparam[IPARAM_DAG] = 0;
        }
        else
            RunTest( iparam, dparam, &(t[iter]));

        gflops = flops / t[iter];

        sumt   += t[iter];
        sumgf  += gflops;
        sumgf2 += gflops*gflops;
    }

    gflops = sumgf / niter;
    sd = sqrt((sumgf2 - (sumgf*sumgf)/niter)/niter);

    if ( iparam[IPARAM_CHECK] )
        printf( "%9.3f %9.2f %9.2f %8.5e %8.5e %8.5e %8.5e %8.5e\n",
                sumt/niter, gflops, sd, dparam[IPARAM_RES], dparam[IPARAM_ANORM], dparam[IPARAM_XNORM], dparam[IPARAM_BNORM],
                dparam[IPARAM_RES] / n / eps / (dparam[IPARAM_ANORM] * dparam[IPARAM_XNORM] + dparam[IPARAM_BNORM] ));
    else{
        printf( "%9.3f %9.2f %9.2f\n", sumt/niter, gflops, sd );
    }
    fflush( stdout );
    free(t);

    return 0;
}

static int
startswith(const char *s, const char *prefix) {
    size_t n = strlen( prefix );
    if (strncmp( s, prefix, n ))
        return 0;
    return 1;
}

static int
get_range(char *range, int *start_p, int *stop_p, int *step_p) {
    char *s, *s1, buf[21];
    int colon_count, copy_len, nbuf=20, n;
    int start=1000, stop=10000, step=1000;

    colon_count = 0;
    for (s = strchr( range, ':'); s; s = strchr( s+1, ':'))
        colon_count++;

    if (colon_count == 0) { /* No colon in range. */
        if (sscanf( range, "%d", &start ) < 1 || start < 1)
            return -1;
        step = start / 10;
        if (step < 1) step = 1;
        stop = start + 10 * step;

    } else if (colon_count == 1) { /* One colon in range.*/
        /* First, get the second number (after colon): the stop value. */
        s = strchr( range, ':' );
        if (sscanf( s+1, "%d", &stop ) < 1 || stop < 1)
            return -1;

        /* Next, get the first number (before colon): the start value. */
        n = s - range;
        copy_len = n > nbuf ? nbuf : n;
        strncpy( buf, range, copy_len );
        buf[copy_len] = 0;
        if (sscanf( buf, "%d", &start ) < 1 || start > stop || start < 1)
            return -1;

        /* Let's have 10 steps or less. */
        step = (stop - start) / 10;
        if (step < 1)
            step = 1;
    } else if (colon_count == 2) { /* Two colons in range. */
        /* First, get the first number (before the first colon): the start value. */
        s = strchr( range, ':' );
        n = s - range;
        copy_len = n > nbuf ? nbuf : n;
        strncpy( buf, range, copy_len );
        buf[copy_len] = 0;
        if (sscanf( buf, "%d", &start ) < 1 || start < 1)
            return -1;

        /* Next, get the second number (after the first colon): the stop value. */
        s1 = strchr( s+1, ':' );
        n = s1 - (s + 1);
        copy_len = n > nbuf ? nbuf : n;
        strncpy( buf, s+1, copy_len );
        buf[copy_len] = 0;
        if (sscanf( buf, "%d", &stop ) < 1 || stop < start)
            return -1;

        /* Finally, get the third number (after the second colon): the step value. */
        if (sscanf( s1+1, "%d", &step ) < 1 || step < 1)
            return -1;
    } else

        return -1;

    *start_p = start;
    *stop_p = stop;
    *step_p = step;

    return 0;
}

static void
show_help(char *prog_name) {
    printf( "Usage:\n%s [options]\n\n", prog_name );
    printf( "Options are:\n"
            "  --help           Show this help\n"
            "\n"
            "  --threads=X      Number of CPU workers (default: _SC_NPROCESSORS_ONLN)\n"
            "\n"
            "  --[a]sync        Enable/Disable synchronous calls in wrapper function such as POTRI. (default: async)\n"
            "  --[no]check      Check result (default: nocheck)\n"
            "  --[no]inv        Check on inverse (default: noinv)\n"
            "  --[no]warmup     Perform a warmup run to pre-load libraries (default: warmup)\n"
#if defined(PLASMA_EZTRACE)
            "  --[no]trace      Enable/Disable trace generation (default: notrace)\n"
#endif
            "  --[no]dag        Enable/Disable DAG generation (default: nodag)\n"
            "                   Generates a dot_dag_file.dot.\n"
            "  --[no]atun       Activate autotuning (default: noatun) [Deprecated]\n"
            "  --inplace        Enable layout conversion inplace for lapack interface timers (default: enable)\n"
            "  --outplace       Enable layout conversion out of place for lapack interface timers (default: disable)\n"
            "\n"
            "  --n_range=R      Range of N values\n"
            "                   with R=Start:Stop:Step (default: 500:5000:500)\n"
            "  --m=X            dimension (M) of the matrices (default: N)\n"
            "  --k=X            dimension (K) of the matrices (default: 1)\n"
            "  --nrhs=X         Number of right-hand size (default: 1)\n"
            "  --nb=N           Nb size. (default: 128)\n"
            "  --ib=N           IB size. (default: 32)\n"
            "\n"
            "  --niter=N        Number of iterations performed for each test (default: 1)\n"
            "\n"
            "  --rhblk=N        If N > 0, enable Householder mode for QR and LQ factorization\n"
            "                   N is the size of each subdomain (default: 0)\n"
            "  --tntmode=x      Factorization used to select the pivots in the tournament pivoting algorithm (default: 1)\n"
            "                     1 - LU factorization with partial pivoting\n"
            "                     2 - QR rank revealing factorization\n"
            "  --tntsize=N      If N > 0, number of tiles grouped together during the tournament (default=4)\n"
            "\n"
            " Options specific to the conversion format timings xgetri and xgecfi:\n"
            "  --ifmt           Input format. (default: 0)\n"
            "  --ofmt           Output format. (default: 1)\n"
            "                   The possible values are:\n"
            "                     0 - PlasmaCM, Column major\n"
            "                     1 - PlasmaCCRB, Column-Colum rectangular block\n"
            "                     2 - PlasmaCRRB, Column-Row rectangular block\n"
            "                     3 - PlasmaRCRB, Row-Colum rectangular block\n"
            "                     4 - PlasmaRRRB, Row-Row rectangular block\n"
            "                     5 - PlasmaRM, Row Major\n"
            "  --thrdbypb       Number of threads per subproblem for inplace transformation (default: 1)\n"
            "\n");
}

static void
print_header(char *prog_name, int * iparam) {
    _PREC    eps = _LAMCH( 'e' );
    printf( "#\n"
            "# PLASMA %d.%d.%d, %s\n"
            "# Nb threads: %d\n"
            "# NB:         %d\n"
            "# IB:         %d\n"
            "# eps:        %e\n"
            "#\n",
            PLASMA_VERSION_MAJOR,
            PLASMA_VERSION_MINOR,
            PLASMA_VERSION_MICRO,
            prog_name,
            iparam[IPARAM_THRDNBR],
            iparam[IPARAM_NB],
            iparam[IPARAM_IB],
            eps );
    if  (iparam[IPARAM_CHECK] )
        printf( "#     M       N  K/NRHS   seconds   Gflop/s Deviation"
                "    ||Ax-b||       ||A||       ||x||       ||b|| ||Ax-b||/N/eps/(||A||||x||+||b||)\n" );
    else
        printf( "#     M       N  K/NRHS   seconds   Gflop/s Deviation\n" );

    return;
}

static void
get_thread_count(int *thrdnbr) {
#if defined WIN32 || defined WIN64
    sscanf( getenv( "NUMBER_OF_PROCESSORS" ), "%d", thrdnbr );
#else
    *thrdnbr = sysconf(_SC_NPROCESSORS_ONLN);
#endif
}

int
main(int argc, char *argv[]) {
    int i, m, mx, nx;
    int start =  500;
    int stop  = 5000;
    int step  =  500;
    int iparam[IPARAM_SIZEOF];

    memset(iparam, 0, IPARAM_SIZEOF*sizeof(int));

    iparam[IPARAM_THRDNBR       ] = 1;
    iparam[IPARAM_THRDNBR_SUBGRP] = 1;
    iparam[IPARAM_SCHEDULER     ] = 1;
    iparam[IPARAM_RUNTIME       ] = 2;
    iparam[IPARAM_M             ] = -1;
    iparam[IPARAM_N             ] = 500;
    iparam[IPARAM_K             ] = 500;
    iparam[IPARAM_LDA           ] = 500;
    iparam[IPARAM_LDB           ] = 500;
    iparam[IPARAM_LDC           ] = 500;
    iparam[IPARAM_MB            ] = 128;
    iparam[IPARAM_NB            ] = 128;
    iparam[IPARAM_IB            ] = 32;
    iparam[IPARAM_NITER         ] = 1;
    iparam[IPARAM_WARMUP        ] = 1;
    iparam[IPARAM_CHECK         ] = 1;
    iparam[IPARAM_VERBOSE       ] = 1;
    iparam[IPARAM_AUTOTUNING    ] = 1;
    iparam[IPARAM_INPUTFMT      ] = 0;
    iparam[IPARAM_OUTPUTFMT     ] = 0;
    iparam[IPARAM_TRACE         ] = 0;
    iparam[IPARAM_DAG           ] = 0;
    iparam[IPARAM_ASYNC         ] = 1;
    iparam[IPARAM_MX            ] = -1;
    iparam[IPARAM_NX            ] = -1;
    iparam[IPARAM_RHBLK         ] = 0;
    iparam[IPARAM_TNTPIV_MODE   ] = PLASMA_TOURNAMENT_LU;
    iparam[IPARAM_TNTPIV_SIZE   ] = 4;
    iparam[IPARAM_INPLACE       ] = PLASMA_OUTOFPLACE;

    get_thread_count( &(iparam[IPARAM_THRDNBR]) );

    for (i = 1; i < argc && argv[i]; ++i) {
        if (startswith( argv[i], "--help" )) {
            show_help( argv[0] );
            return EXIT_SUCCESS;
        } else if (startswith( argv[i], "--threads=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_THRDNBR]) );
        } else if (startswith( argv[i], "--check" )) {
            iparam[IPARAM_CHECK] = 1;
        } else if (startswith( argv[i], "--nocheck" )) {
            iparam[IPARAM_CHECK] = 0;
        } else if (startswith( argv[i], "--warmup" )) {
            iparam[IPARAM_WARMUP] = 1;
        } else if (startswith( argv[i], "--nowarmup" )) {
            iparam[IPARAM_WARMUP] = 0;
        } else if (startswith( argv[i], "--dyn" )) {
            iparam[IPARAM_SCHEDULER] = 1;
        } else if (startswith( argv[i], "--nodyn" )) {
            iparam[IPARAM_SCHEDULER] = 0;
        } else if (startswith( argv[i], "--quark" )) {
            iparam[IPARAM_RUNTIME] = 1;
        } else if (startswith( argv[i], "--ompss" )) {
            iparam[IPARAM_RUNTIME] = 2;
        /* } else if (startswith( argv[i], "--atun" )) { */
        /*     iparam[IPARAM_AUTOTUNING] = 1; */
        /* } else if (startswith( argv[i], "--noatun" )) { */
        /*     iparam[IPARAM_AUTOTUNING] = 0; */
        } else if (startswith( argv[i], "--trace" )) {
            iparam[IPARAM_TRACE] = 1;
        } else if (startswith( argv[i], "--notrace" )) {
            iparam[IPARAM_TRACE] = 0;
        } else if (startswith( argv[i], "--dag" )) {
            iparam[IPARAM_DAG] = 1;
        } else if (startswith( argv[i], "--nodag" )) {
            iparam[IPARAM_DAG] = 0;
        } else if (startswith( argv[i], "--sync" )) {
            iparam[IPARAM_ASYNC] = 0;
        } else if (startswith( argv[i], "--async" )) {
            iparam[IPARAM_ASYNC] = 1;
        } else if (startswith( argv[i], "--n_range=" )) {
            get_range( strchr( argv[i], '=' ) + 1, &start, &stop, &step );
        } else if (startswith( argv[i], "--m=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_M]) );
        } else if (startswith( argv[i], "--n=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_N]) );
            start = iparam[IPARAM_N];
            stop  = iparam[IPARAM_N];
            step  = 1;
        } else if (startswith( argv[i], "--nb=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_NB]) );
            iparam[IPARAM_MB] = iparam[IPARAM_NB];
        } else if (startswith( argv[i], "--nrhs=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_K]) );
        } else if (startswith( argv[i], "--k=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_K]) );
        } else if (startswith( argv[i], "--ib=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_IB]) );
        } else if (startswith( argv[i], "--ifmt=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_INPUTFMT]) );
        } else if (startswith( argv[i], "--ofmt=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_OUTPUTFMT]) );
        } else if (startswith( argv[i], "--thrdbypb=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_THRDNBR_SUBGRP]) );
        } else if (startswith( argv[i], "--niter=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &iparam[IPARAM_NITER] );
        } else if (startswith( argv[i], "--mx=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_MX]) );
        } else if (startswith( argv[i], "--nx=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_NX]) );
        } else if (startswith( argv[i], "--rhblk=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_RHBLK]) );
        } else if (startswith( argv[i], "--tntmode=" )) {
            int tmp;
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &tmp );
            iparam[IPARAM_TNTPIV_MODE] = (tmp == 1) ? PLASMA_TOURNAMENT_LU : PLASMA_TOURNAMENT_QR;
        } else if (startswith( argv[i], "--tntsize=" )) {
            sscanf( strchr( argv[i], '=' ) + 1, "%d", &(iparam[IPARAM_TNTPIV_SIZE]) );
        } else if (startswith( argv[i], "--inplace" )) {
            iparam[IPARAM_INPLACE] = PLASMA_INPLACE;
        } else if (startswith( argv[i], "--outplace" )) {
            iparam[IPARAM_INPLACE] = PLASMA_OUTOFPLACE;
        } else {
            fprintf( stderr, "Unknown option: %s\n", argv[i] );
        }
    }

    m  = iparam[IPARAM_M];
    mx = iparam[IPARAM_MX];
    nx = iparam[IPARAM_NX];

    /* Initialize Plasma */
	if ( iparam[IPARAM_RUNTIME] == PLASMA_QUARK )
		PLASMA_Runtime_Init(omp_get_num_threads(), PLASMA_QUARK);
	else
		PLASMA_Runtime_Init(omp_get_num_threads(), PLASMA_OMPSS);

    if ( iparam[IPARAM_SCHEDULER] ){
        PLASMA_Set(PLASMA_SCHEDULING_MODE, PLASMA_DYNAMIC_SCHEDULING );
        if ( iparam[IPARAM_RUNTIME] == PLASMA_QUARK )
            PLASMA_Set(PLASMA_RUNTIME_MODE, PLASMA_QUARK );
        else
            PLASMA_Set(PLASMA_RUNTIME_MODE, PLASMA_OMPSS );
    }
    else
        PLASMA_Set(PLASMA_SCHEDULING_MODE, PLASMA_STATIC_SCHEDULING );

    /* if ( !iparam[IPARAM_AUTOTUNING] ) { */
    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE,        iparam[IPARAM_NB] );
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, iparam[IPARAM_IB] );
    /* } else { */
    /*     PLASMA_Get(PLASMA_TILE_SIZE,        &iparam[IPARAM_NB] ); */
    /*     PLASMA_Get(PLASMA_INNER_BLOCK_SIZE, &iparam[IPARAM_IB] ); */
    /* } */

    /* Householder mode */
    if (iparam[IPARAM_RHBLK] < 1) {
        PLASMA_Set(PLASMA_HOUSEHOLDER_MODE, PLASMA_FLAT_HOUSEHOLDER);
    } else {
        PLASMA_Set(PLASMA_HOUSEHOLDER_MODE, PLASMA_TREE_HOUSEHOLDER);
        PLASMA_Set(PLASMA_HOUSEHOLDER_SIZE, iparam[IPARAM_RHBLK]);
    }

    /* Tournament pivoting */
    PLASMA_Set(PLASMA_TNTPIVOTING_MODE, iparam[IPARAM_TNTPIV_MODE]);
    PLASMA_Set(PLASMA_TNTPIVOTING_SIZE, iparam[IPARAM_TNTPIV_SIZE]);

    /* Layout conversion */
    PLASMA_Set(PLASMA_TRANSLATION_MODE, iparam[IPARAM_INPLACE]);

    print_header( argv[0], iparam);

    if (step < 1) step = 1;

    Test( -1, iparam ); /* print header */
    for (i = start; i <= stop; i += step)
    {
        if ( nx > 0 ) {
            iparam[IPARAM_M] = i;
            iparam[IPARAM_N] = max(1, i/nx);
        } else if ( mx > 0 ) {
            iparam[IPARAM_M] = max(1, i/mx);
            iparam[IPARAM_N] = i;
        } else {
            if ( m == -1 )
                iparam[IPARAM_M] = i;
            iparam[IPARAM_N] = i;
        }
        Test( iparam[IPARAM_N], iparam );
    }

    PLASMA_Finalize();

    /* if (gnuplot) { */
    /*         printf( "%s\n%s\n", */
    /*                 "e", */
    /*                 gnuplot > 1 ? "" : "pause 10" ); */
    /* } */

    return EXIT_SUCCESS;
}


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

#ifndef __GSL_MULTIFIT_H__
#define __GSL_MULTIFIT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t n; /* number of observations */
  size_t p; /* number of parameters */
  gsl_matrix * A;
  gsl_matrix * Q;
  gsl_matrix * QSI;
  gsl_vector * S;
  gsl_vector * t;
  gsl_vector * xt;
  gsl_vector * D;
}
gsl_multifit_linear_workspace;

GSL_FUN gsl_multifit_linear_workspace *
gsl_multifit_linear_alloc (size_t n, size_t p);

GSL_FUN void
gsl_multifit_linear_free (gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear (const gsl_matrix * X,
                     const gsl_vector * y,
                     gsl_vector * c,
                     gsl_matrix * cov,
                     double * chisq,
                     gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_svd (const gsl_matrix * X,
                         const gsl_vector * y,
                         double tol,
                         size_t * rank,
                         gsl_vector * c,
                         gsl_matrix * cov,
                         double *chisq,
                         gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double * chisq,
                      gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq,
                          gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_est (const gsl_vector * x,
                         const gsl_vector * c,
                         const gsl_matrix * cov, double *y, double *y_err);

GSL_FUN int
gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
                               const gsl_vector *c, gsl_vector *r);

__END_DECLS

#endif /* __GSL_MULTIFIT_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_histogram.h>
#include"mpi.h"

const double PI = 3.14159;

double rz(double red);

struct galaxy {
  double x, y, z, d, w;
};

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

  FILE * DD;
  FILE * PP;
  FILE * RR;

  const int Nr = 150;
  const int Nm = 100;
  const double RMIN = 0.0;
  const double RMAX = 150.0;
  const double RSTEP = 1.0;
  const double MSTEP = 0.01;

  const int maxNgal = 40000000;

  int dummy;

  char DDname[400];
  char PPname[400];
  char RRname[400];

  int i;
  int bit;
  int bittot;

  int *n;

  sprintf(DDname,argv[1]);
  sprintf(RRname,argv[2]);
  sprintf(PPname,argv[3]);

  struct galaxy * gald;
  struct galaxy * galr;

  double * DRcount;
  double * DRcount_in_place;
  gald = (struct galaxy*)malloc(maxNgal*sizeof(struct galaxy));
  galr = (struct galaxy*)malloc(maxNgal*sizeof(struct galaxy));
  DRcount = (double*)calloc(Nr*Nm,sizeof(double));
  DRcount_in_place = (double*)calloc(Nr*Nm,sizeof(double));
  struct galaxy * galpd, * galpr;
  galpd = gald;
  galpr = galr;

  DD = fopen(DDname,"r");
  RR = fopen(RRname,"r");

  double ra, dec, red, dist;
  double dd;
  int Nd = 0;
  double weight;
  double ddummy;
  int indexFKP;
  double weightFKP;
  double weight1, weight2, weight3, weight4;
  char line [2000];

  int N = 0;
  double sumW = 0;
  double sumW_in_place = 0;

  for (int i = 0; i < 0; i ++) {
  fgets (line, 2000, DD);
  printf("%s\n", line);
  }
  while (fscanf(DD,"%lf %lf %lf\n",&ra,&dec,&red) != EOF) {
    ra *= PI/180.0;
    dec *= PI/180.0;
    dist = rz(red);
    galpd->x = dist*cos(dec)*cos(ra);
    galpd->y = dist*cos(dec)*sin(ra);
    galpd->z = dist*sin(dec);
    galpd->d = dist;
//    galpd->w = weight1 * weight2;
    galpd++;
    Nd++;
  }
  fclose(DD);

  N = 0;
  for (int i = 0; i < 0; i ++) {
  fgets (line, 2000, RR);
  printf("%s\n", line);
  }
  while (fscanf(RR,"%le %le %le\n",&ra,&dec,&red) != EOF) {
    ra *= PI/180.0;
    dec *= PI/180.0;
    dist = rz(red);
    galpr->x = dist*cos(dec)*cos(ra);
    galpr->y = dist*cos(dec)*sin(ra);
    galpr->z = dist*sin(dec);
    galpr->d = dist;
//    galpr->w = weight1 * weight2;
    galpr++;
    N++;
  }
  fclose(RR);

  MPI_Init (&argc, &argv);
  int numproc;
  int totproc;
  MPI_Comm_rank (MPI_COMM_WORLD, &numproc);
  MPI_Comm_size (MPI_COMM_WORLD, &totproc);

  n = (int *) calloc(totproc + 1, sizeof(int));
  n[0] = 1;
  for (i = 1; i < totproc; i++) {
    n[i] = (int)floor(double(N)/totproc*i);
  }
  n[totproc] = N;
  bit = n[numproc];
  bittot = n[numproc + 1];

  sprintf(PPname, "%s.DR",PPname);
  double mu;
  double d1, d2, d3;
  double x1, y1, z1, x2, y2, z2, gd1, gd2, w1, w2;
  double r2, rat;
  double xb, yb, zb, db2;
  double rr;
  galpd = gald;
  int j;

  for (i = 0; i < Nd; i++) {
    x1 = galpd->x; y1 = galpd->y; z1 = galpd->z; gd1 = galpd->d; //w1 = galpd->w;
    galpr = galr + bit;
    for (j = bit; j < bittot; j++) {
       x2 = galpr->x; y2 = galpr->y; z2 = galpr->z; gd2 = galpr->d; //w2 = galpr->w;
       //sumW += w1*w2;
       sumW += 1.0;
       d1 = x1 - x2;
       d2 = y1 - y2;
       d3 = z1 - z2;
       r2 = d1*d1 + d2*d2 + d3*d3;
       if (r2 > RMAX*RMAX || r2 < RMIN*RMIN) {
	 galpr++;
         continue;
       }
       rat = gd1/gd2;
       xb = x1 + x2*rat;
       yb = y1 + y2*rat;
       zb = z1 + z2*rat;
       db2 = xb*xb + yb*yb + zb*zb;
       mu = fabs((xb*d1 + yb*d2 + zb*d3)/sqrt(r2)/sqrt(db2));
       rr = sqrt(r2);
       int binr = (int)((rr - RMIN)/RSTEP);
       int binm = (int)(mu/MSTEP);
       if (binr >= 0 && binm >= 0 && binr < Nr && binm < Nm){
         int ind = binr + Nr*binm;
//         DRcount[ind] += w1*w2;
         DRcount[ind] += 1.0;
       }
       galpr++;
    }
    galpd++;
  }
  free(gald);
  free(galr);
  free(n);

  MPI_Reduce(DRcount, DRcount_in_place, Nr*Nm, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&sumW, &sumW_in_place, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  printf("weighted number = %le %d\n", sumW, numproc);
  if (numproc == 0) {
    PP = fopen(PPname,"w");
    printf("%s\n", PPname);
    fprintf(PP,"# weighted number: %lf\n",sumW_in_place);
    fprintf(PP,"# RBINS: %d\n", Nr);
    fprintf(PP,"# MBINS: %d\n", Nm);
    int k, l;
    for (k = 0; k < Nm; k++) {
      for (l = 0; l < Nr; l++) {
      fprintf(PP,"%lf ",DRcount_in_place[k*Nr + l]);
      }
      fprintf(PP,"\n");
    }
    fclose(PP);
  }

  free(DRcount);
  free(DRcount_in_place);

  MPI_Finalize ();

  return 0;
}

double f (double x, void * p) {
  double Om = 0.310;
  double ff = 2997.92458/sqrt(Om*(1.0 + x)*(1.0 + x)*(1.0 + x) + 1.0 - Om);
  return ff;
}

double rz (double red) {
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  double result, error;
  gsl_function F;
  F.function = &f;
  gsl_integration_qags(&F, 0, red, 0, 1e-7, 1000, w, &result, &error);
  gsl_integration_workspace_free (w);
  return result;
}

ext: c

max_stars_repo_name: echaussidon/LSS

max_issues_repo_name: echaussidon/LSS

max_forks_repo_name: echaussidon/LSS


text:
/* -*- coding: utf-8; mode: c++; tab-width: 2; -*- -------------------------
   Interface/declaration file for GramSchmidt.cpp
------------------------------------------------------------------------- */

#ifndef GRAMSCHMIDT_H
#define GRAMSCHMIDT_H

#include <vector>
#include <numeric>
#include <algorithm>
#include <cblas.h>
#include <lapacke.h>
#include <iostream>

// TO DO: Improve performance with CBLAS/LAPACKE for large size x
// std::vector<double> v;
// double* a = &v[0];

// int GramSchmidt(size_t, size_t, void*, void*);
int GramSchmidt(size_t, size_t, double*, double*);

// int GramSchmidt(std::vector<std::vector<double>>, std::vector<std::vector<double>>);

#endif


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2015 Rokko Developers https://github.com/t-sakashita/rokko
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#ifndef ROKKO_PBLAS_INTERFACE_H
#define ROKKO_PBLAS_INTERFACE_H

#include <lapacke.h>
#undef I
#include <rokko/mangling.h>

/* PvCOPY */

#define PBLAS_PCOPY_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const int* N, const TYPE * X, const int* IX, const int* JX, const int* DESCX, const int* INCX, TYPE * Y, const int* IY, int* JY, const int* DESCY, const int* INCY);

#ifdef __cplusplus
extern "C" {
#endif

PBLAS_PCOPY_DECL(pscopy, PSCOPY, float);
PBLAS_PCOPY_DECL(pdcopy, PDCOPY, double);
PBLAS_PCOPY_DECL(pccopy, PCCOPY, lapack_complex_float);
PBLAS_PCOPY_DECL(pzcopy, PZCOPY, lapack_complex_double);

#ifdef __cplusplus
}
#endif

#undef PBLAS_PCOPY_DECL

/* PvDOT, PvDOTU, PvDOTC */

#define PBLAS_PDOT_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const int* N, TYPE * DOT, const TYPE * X, const int* IX, const int* JX, const int* DESCX, const int* INCX, const TYPE * Y, const int* IY, const int* JY, const int* DESCY, const int* INCY);

#ifdef __cplusplus
extern "C" {
#endif

PBLAS_PDOT_DECL(psdot, PSDOT, float);
PBLAS_PDOT_DECL(pddot, PDDOT, double);
PBLAS_PDOT_DECL(pcdotu, PCDOTU, lapack_complex_float);
PBLAS_PDOT_DECL(pzdotu, PZDOTU, lapack_complex_double);
PBLAS_PDOT_DECL(pcdotc, PCDOTC, lapack_complex_float);
PBLAS_PDOT_DECL(pzdotc, PZDOTC, lapack_complex_double);

#ifdef __cplusplus
}
#endif

#undef PBLAS_PDOT_DECL

/* PvGEMV */

#define PBLAS_PGEMV_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* TRANS, const int* M, const int* N, const TYPE * ALPHA, const TYPE * A, const int* IA, const int* JA, const int* DESCA, const TYPE * X, const int* IX, const int* JX, const int* DESCX, const int* INCX, const TYPE * BETA, TYPE * Y, const int* IY, const int* JY, const int* DESCY, const int* INCY);

#ifdef __cplusplus
extern "C" {
#endif

PBLAS_PGEMV_DECL(psgemv, PSGEMV, float);
PBLAS_PGEMV_DECL(pdgemv, PDGEMV, double);
PBLAS_PGEMV_DECL(pcgemv, PCGEMV, lapack_complex_float);
PBLAS_PGEMV_DECL(pzgemv, PZGEMV, lapack_complex_double);

#ifdef __cplusplus
}
#endif

#undef PBLAS_PGEMV_DECL

/* PvGEMM */

#define PBLAS_PGEMM_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* TRANSA, const char* TRANSB, const int* M, const int* N, const int* K, const TYPE * ALPHA, const TYPE * A, const int* IA, const int* JA, const int* DESCA, const TYPE * B, const int* IB, const int* JB, const int* DESCB, const TYPE * BETA, TYPE * C, const int* IC, const int* JC, const int* DESCC);

#ifdef __cplusplus
extern "C" {
#endif

PBLAS_PGEMM_DECL(psgemm, PSGEMM, float);
PBLAS_PGEMM_DECL(pdgemm, PDGEMM, double);
PBLAS_PGEMM_DECL(pcgemm, PCGEMM, lapack_complex_float);
PBLAS_PGEMM_DECL(pzgemm, PZGEMM, lapack_complex_double);

#ifdef __cplusplus
}
#endif

#undef PBLAS_PGEMM_DECL

#endif // ROKKO_PBLAS_INTERFACE_H

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


text:
#include <petsc.h>

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{

  PetscErrorCode ierr;

  PetscInitialize(&argc,&args,(char*)0,(char*)0);

  ierr = PetscPrintf(PETSC_COMM_WORLD," === Hello from Petsc! ===\n");CHKERRQ(ierr);

  ierr = PetscFinalize();

  return 0;
}


text:
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>
#include <gbpLib.h>
#include <gbpRNG.h>
#include <gbpMCMC.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_interp.h>

void read_MCMC_configuration(MCMC_info *MCMC, char *filename_output_dir, int chain) {
    SID_log("Reading the configuration of chain #%d from {%s}...", SID_LOG_OPEN, chain);

    // Set the names of the files we have to read
    char filename_chain_dir[SID_MAX_FILENAME_LENGTH];
    char filename_run[SID_MAX_FILENAME_LENGTH];
    char filename_chain_config[SID_MAX_FILENAME_LENGTH];
    sprintf(filename_chain_dir, "%s/chains/", filename_output_dir);
    sprintf(filename_run, "%s/run.dat", filename_output_dir);
    sprintf(filename_chain_config, "%s/chain_config_%06d.dat", filename_chain_dir, chain);

    // Read the run file to do some sanity checks
    int   n_chains;
    int   n_avg;
    int   flag_autocor_on;
    int   flag_no_map_write;
    int   n_P;
    FILE *fp;
    fp = fopen(filename_run, "r");
    SID_fread_verify(&n_chains, sizeof(int), 1, fp);
    SID_fread_verify(&n_avg, sizeof(int), 1, fp);
    SID_fread_verify(&flag_autocor_on, sizeof(int), 1, fp);
    SID_fread_verify(&flag_no_map_write, sizeof(int), 1, fp);
    SID_fread_verify(&n_P, sizeof(int), 1, fp);
    fclose(fp);
    if(n_P != MCMC->n_P)
        SID_exit_error("The number of parameters for the two runs does not match (ie. %d!=%d)", SID_ERROR_LOGIC, n_P,
                       MCMC->n_P);
    if(chain >= n_chains)
        SID_exit_error("Invalid chain number (ie. %d>%d)", SID_ERROR_LOGIC, chain, n_chains);

    // Read the configuration file
    int     n_iterations_file_total;
    int     n_iterations_file_burn;
    double *V;
    double  temperature;
    V  = (double *)SID_malloc(sizeof(double) * n_P * n_P);
    fp = fopen(filename_chain_config, "r");
    SID_fread_verify(&n_iterations_file_total, sizeof(int), 1, fp);
    SID_fread_verify(&n_iterations_file_burn, sizeof(int), 1, fp);
    SID_fread_verify(V, sizeof(double), n_P * n_P, fp);
    SID_fread_verify(&temperature, sizeof(double), 1, fp);
    fclose(fp);

    // Set the new state
    set_MCMC_covariance(MCMC, V);
    set_MCMC_temperature(MCMC, temperature);

    // Clean-up
    SID_free(SID_FARG V);
    SID_log("Done.", SID_LOG_CLOSE);
}

max_stars_repo_name: gbpoole/gbpCode

max_issues_repo_name: gbpoole/gbpCode

max_forks_repo_name: gbpoole/gbpCode


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gbpInterpolate.h>

void init_interpolate(double *x, double *y, size_t n, const gsl_interp_type *T, interp_info **interp) {
    size_t i;
    (*interp)    = (interp_info *)malloc(sizeof(interp_info));
    (*interp)->n = n;
    (*interp)->x = (double *)malloc(sizeof(double) * n);
    (*interp)->y = (double *)malloc(sizeof(double) * n);
    (*interp)->T = T;
    // If the x-array is not in ascending order, switch it
    if(x[0] < x[n - 1]) {
        for(i = 0; i < n; i++) {
            (*interp)->x[i] = x[i];
            (*interp)->y[i] = y[i];
        }
    } else {
        for(i = 0; i < n; i++) {
            (*interp)->x[i] = x[n - 1 - i];
            (*interp)->y[i] = y[n - 1 - i];
        }
    }
    (*interp)->accel  = gsl_interp_accel_alloc();
    (*interp)->interp = gsl_interp_alloc(T, n);
    gsl_interp_init((*interp)->interp, (const double *)(*interp)->x, (const double *)(*interp)->y, (*interp)->n);
}


text:
/**
 *
 * @file core_zlacpy.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Julien Langou
 * @author Henricus Bouwmeester
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 *  CORE_PLASMA_zlacpycopies all or part of a two-dimensional matrix A to another
 *  matrix B
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies the part of the matrix A to be copied to B.
 *            = PlasmaUpperLower: All the matrix A
 *            = PlasmaUpper: Upper triangular part
 *            = PlasmaLower: Lower triangular part
 *
 * @param[in] M
 *          The number of rows of the matrices A and B. M >= 0.
 *
 * @param[in] N
 *          The number of columns of the matrices A and B. N >= 0.
 *
 * @param[in] A
 *          The M-by-N matrix to copy.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[out] B
 *          The M-by-N copy of the matrix A.
 *          On exit, B = A ONLY in the locations specified by uplo.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDB >= max(1,M).
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zlacpy = PCORE_zlacpy
#define CORE_zlacpy PCORE_zlacpy
#endif
void CORE_zlacpy(PLASMA_enum uplo, int M, int N,
                 const PLASMA_Complex64_t *A, int LDA,
                 PLASMA_Complex64_t *B, int LDB)
{
    LAPACKE_zlacpy_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        M, N, A, LDA, B, LDB);
}