text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <string.h>
#include <mpi.h>
#include "parmt_mtsearch.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "parmt_utils.h"
#include "cps_mpe.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"

int parmt_obsSearch64f(const MPI_Comm globalComm,
                       const MPI_Comm obsComm, const MPI_Comm locComm,
                       const bool linObsComm, const bool linLocComm,
                       const struct parmtGeneralParms_struct parms,
                       struct parmtData_struct data,
                       struct localMT_struct mtloc,
                       double *__restrict__ phi)
{
    MPI_Win counterWin;
    const char *archiveName = NULL;
    double lagTime, *CeInv, *phiWork, *phiAll, *var, *varAll, *varOut, *wts,
           defaultWeight, wtNorm;
    int *lagsWork, *nlags, ierr, iobs, myid, nobsProcs,
        npmax, value, value0;
    bool ldefault, lwantLags;
    const double half = 0.5;
    const int master = 0;
    //------------------------------------------------------------------------//
    //
    // Have the root process determine lag information
    ierr = 0;
    MPI_Comm_rank(globalComm, &myid);
    nlags = memory_calloc32i(data.nobs);
    nobsProcs = 0;
    lagsWork = NULL;
    phiWork = NULL;
    phiAll = NULL;
    CeInv = NULL;
    var = NULL;
    varAll = NULL;
    varOut = NULL;
    wts = NULL;
    npmax = 0;
    lwantLags = false;
    if (myid == master)
    {
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            nlags[iobs] = 0;
            if (parms.lwantLags)
            {
                //lagTime = data.data[iobs].header.user0;
                //if (lagTime < zero){lagTime = parms.defaultMaxLagTime;}
                lagTime = parmt_utils_getLagTime(data.data[iobs],
                                                 parms.defaultMaxLagTime,
                                                 &ldefault);
                nlags[iobs] = (int)
                              (lagTime/data.data[iobs].header.delta + half);
            }
            lwantLags = false;
            if (nlags[iobs] > 0){lwantLags = true;}
            npmax = MAX(npmax, data.data[iobs].npts);
            if (parms.lwantLags)
            {
                printf("%s: Number of lags for waveform %d is %d\n",
                       __func__, iobs+1, nlags[iobs]);
            }
        }
    }
    MPI_Bcast(nlags, data.nobs, MPI_INT,    master, globalComm);
    MPI_Bcast(&npmax,        1, MPI_INT,    master, globalComm);
    MPI_Bcast(&lwantLags,    1, MPI_C_BOOL, master, globalComm);
    // Check for failures
    if (npmax < 1 || mtloc.nmtAll < 1 || parms.objFnType < 1 ||
        parms.objFnType > 3)
    {
        if (npmax < 1){printf("%s: No data points\n", __func__);}
        if (mtloc.nmtAll < 1){printf("%s: Empty mts\n", __func__);}
        if (parms.objFnType < 1 || parms.objFnType > 3)
        {
            printf("%s: Invalid objfn function\n", __func__);
        }
        return -1;
    }
    if (lwantLags)
    {
        if (mtloc.myid == master)
        {
            lagsWork = memory_calloc32i(data.nlocs*mtloc.nmtAll);
        }
        else
        {
            lagsWork = memory_calloc32i(1);
        }
    }
    // Set space
    var = memory_calloc64f(npmax);
    if (linObsComm){MPI_Comm_size(obsComm, &nobsProcs);}
    if (linObsComm)
    {
        if (nobsProcs > 1)
        {
            phiAll = memory_calloc64f(data.nlocs*mtloc.nmtAll);
        }
        else
        {
            if (phi == NULL)
            {
                printf("%s: phi not set on master - segfault\n", __func__);
                return -1;
            }
            phiAll = phi;
            array_zeros64f_work(data.nlocs*mtloc.nmtAll, phiAll);
        }
        phiWork = memory_calloc64f(data.nlocs*mtloc.nmtAll);
        varAll = memory_calloc64f(data.nobs*npmax);
        defaultWeight = 1.0; ///(double) data.nobs;
        wts = memory_calloc64f(data.nobs); //array_set64f(data.nobs, defaultWeight, &ierr);
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            wts[iobs] = parmt_utils_getWeight(data.data[iobs],
                                              defaultWeight,
                                              &ldefault);
        }
        // Normalize so weighted sum is unity
        wtNorm = 1.0/array_sum64f(data.nobs, wts, &ierr);
        // And rescale the weights
        for (iobs=0; iobs<data.nobs; iobs++)
        {
            wts[iobs] = wts[iobs]*wtNorm;
        }
    }
    else
    {
        phiWork = memory_calloc64f(1);
        phiAll = memory_calloc64f(1);
    }
    // Set up the MPI RMA window
    value = 0;
    value0 = 0;
    if (linObsComm)
    {
        MPE_setKeyval(MPI_KEYVAL_INVALID);
        MPE_counterCreate(obsComm, 1, &counterWin);
    }
    MPI_Barrier(globalComm);
    // Work-steal the observations to accomodate different sizes/lags
    while (true)
    {
        // Update the counter
        if (linObsComm){MPE_counterNextval(counterWin, 0, &value);}
        // Send the counter to the locations commuincator
        if (linLocComm){MPI_Bcast(&value, 1, MPI_INT, master, locComm);}
        // Trickle the value down to the moment tensor communicator
        MPI_Bcast(&value, 1, MPI_INT, master, mtloc.comm);
        if (value > data.nobs){break;}  // We are done
        if (value == value0){continue;} // Already did this
        value0 = value;   // Update the old value
        iobs = value - 1; // Map to C based numbering to indicate observation
        if (linObsComm)
        {
            printf("%s: Processing observation %d with %d points\n", __func__,
                   iobs+1, data.data[iobs].npts);
        }
        // Figure out the lags
        lwantLags = false;
        if (nlags[iobs] > 0){lwantLags = true;}
        if (parms.objFnType == 1)
        {
            ierr = parmt_locSearchL164f(locComm,
                                        iobs, parms.blockSize,
                                        nlags[iobs], lwantLags,
                                        parms.lrescale,
                                        mtloc, &data,
                                        CeInv, phiWork, var, lagsWork);
            if (ierr != 0)
            {
                printf("%s: L1 search error for observation %d\n",
                       __func__, iobs + 1);
                return -1;
            }
            // Put the variance into the global buffer
            if (linObsComm)
            {
                array_copy64f_work(npmax, var, &varAll[iobs*npmax]);
            }
        }
        else if (parms.objFnType == 2)
        {
            ierr = parmt_locSearchXC64f(locComm, iobs, parms.blockSize,
                                        nlags[iobs], lwantLags,
                                        mtloc, &data, phiWork, lagsWork);
            if (ierr != 0)
            {
                printf("%s: XC search error for observation %d\n",
                       __func__, iobs + 1);
                return -1;
            }
        }
        else if (parms.objFnType == 3)
        {
            printf("%s: need to write this function\n", __func__);
            return -1;
        }
        // Scale by number of observations and stack
        if (linObsComm)
        {
            //printf("%d %f\n", iobs, wts[iobs]);
            cblas_daxpy(data.nlocs*mtloc.nmtAll, wts[iobs], //1.0/(double) data.nobs,
                        phiWork, 1, phiAll, 1);
        }
    }
    // Reduce the objective function
    if (linObsComm)
    {
        if (nobsProcs > 1)
        {
            MPI_Reduce(phiAll, phi, data.nlocs*mtloc.nmtAll, MPI_DOUBLE,
                       MPI_SUM, master, obsComm);
        }
        MPI_Win_free(&counterWin);
    }
    MPI_Barrier(mtloc.comm);
    // Release resources
    if (linObsComm)
    {
        if (nobsProcs > 1){memory_free64f(&phiAll);}
        // Reduce the variance
        varOut = NULL;
        if (myid == master){varOut = memory_calloc64f(npmax*data.nobs);}
        MPI_Reduce(varAll, varOut, npmax*data.nobs, MPI_DOUBLE,
                   MPI_SUM, master, obsComm);
        // Have the master write it
        if (myid == master)
        {
            if (parms.programID == PARMT_ID)
            {
                archiveName = parms.parmtArchive;
            }
            else if (parms.programID == POLARMT_ID)
            {
                archiveName = parms.polarmtArchive;
            }
            for (iobs=0; iobs<data.nobs; iobs++)
            {
                parmt_io_writeVarianceForWaveform64f(
                   //parms.resultsDir, parms.projnm, parms.resultsFileSuffix,
                   archiveName,
                   iobs,
                   data.data[iobs].npts, &varOut[iobs*npmax]);
            }
            archiveName = NULL;
        }
        memory_free64f(&varAll);
        memory_free64f(&varOut);
    }
    else
    {
        memory_free64f(&phiAll);
    }
    memory_free64f(&phiWork);
    memory_free32i(&nlags);
    memory_free32i(&lagsWork);
    memory_free64f(&CeInv);
    memory_free64f(&wts);
    phiAll = NULL;
    MPI_Barrier(globalComm);
    return 0;
}


text:
/*	$Id$ */
/*
 * Copyright (c) 2014, 2015 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <stdint.h>
#include <stdlib.h>
#include <string.h>

#ifdef MAC_INTEGRATION
#include <gtkosxapplication.h>
#endif
#include <gtk/gtk.h>
#include <gdk/gdk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

#define DEFINE_GTK_INIT(_class, _type, _name) \
_class * \
win_init_ ##_name (GtkBuilder *b, const gchar *name) \
{ \
	g_assert(gtk_builder_get_object(b, name)); \
	g_assert(GTK_IS_##_type (gtk_builder_get_object(b, name))); \
	return(GTK_ ##_type (gtk_builder_get_object(b, name))); \
}

DEFINE_GTK_INIT(GtkAdjustment, ADJUSTMENT, adjustment)
DEFINE_GTK_INIT(GtkStatusbar, STATUSBAR, status)
DEFINE_GTK_INIT(GtkDrawingArea, DRAWING_AREA, draw)
DEFINE_GTK_INIT(GtkMenuBar, MENU_BAR, menubar)
DEFINE_GTK_INIT(GtkCheckMenuItem, CHECK_MENU_ITEM, menucheck)
DEFINE_GTK_INIT(GtkMenuItem, MENU_ITEM, menuitem)
DEFINE_GTK_INIT(GtkWindow, WINDOW, window)
DEFINE_GTK_INIT(GtkLabel, LABEL, label)
DEFINE_GTK_INIT(GtkBox, BOX, box)
DEFINE_GTK_INIT(GtkButton, BUTTON, button)
DEFINE_GTK_INIT(GtkRadioButton, RADIO_BUTTON, radio)
DEFINE_GTK_INIT(GtkNotebook, NOTEBOOK, notebook)
DEFINE_GTK_INIT(GtkToggleButton, TOGGLE_BUTTON, toggle)
DEFINE_GTK_INIT(GtkEntry, ENTRY, entry)
DEFINE_GTK_INIT(GtkFileChooser, FILE_CHOOSER, filechoose)

GtkBuilder *
builder_get(const gchar *name)
{
	gchar		*file = NULL;
#ifdef MAC_INTEGRATION
	gchar		*dir;
#endif
	GtkBuilder	*builder;

	/*
	 * Look up our `glade' file as follows.
	 * If we're in MAC_INTEGRATION, then intuit whether we're in a
	 * bundle and, if so, look up within the bundle.
	 * If we're not in a bundle, look in DATADIR.
	 */
#ifdef MAC_INTEGRATION
	g_object_new(GTKOSX_TYPE_APPLICATION, NULL);
	if (NULL != (dir = gtkosx_application_get_bundle_id())) {
		g_free(dir);
		dir = gtkosx_application_get_resource_path();
		file = g_strdup_printf
			("%s" G_DIR_SEPARATOR_S
			 "share" G_DIR_SEPARATOR_S
			 "bmigrate" G_DIR_SEPARATOR_S
			 "%s", dir, name);
		g_free(dir);
	}
#endif
	if (NULL == file)
		file = g_strdup_printf(DATADIR
			G_DIR_SEPARATOR_S "%s", name);

	builder = gtk_builder_new();
	assert(NULL != builder);

	/*
	 * This should be gtk_builder_new_from_file(), but we don't
	 * support that on older versions of GTK, so do it like this.
	 */
	if ( ! gtk_builder_add_from_file(builder, file, NULL)) {
		g_free(file);
		g_object_unref(G_OBJECT(builder));
		return(NULL);
	}

	g_free(file);
	return(builder);
}

ext: c

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_syr2k (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { -0.635f, 0.805f };
   int lda = 2;
   float B[] = { 0.773f, 0.375f };
   int ldb = 2;
   float C[] = { 0.616f };
   int ldc = 1;
   float C_expected[] = { 0.174988f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1622)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { -0.396f, -0.131f };
   int lda = 2;
   float B[] = { -0.603f, -0.288f };
   int ldb = 2;
   float C[] = { -0.434f };
   int ldc = 1;
   float C_expected[] = { -0.20931f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1623)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { -0.085f, -0.444f };
   int lda = 1;
   float B[] = { 0.936f, 0.752f };
   int ldb = 1;
   float C[] = { -0.64f };
   int ldc = 1;
   float C_expected[] = { 0.184069f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1624)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { 0.655f, 0.45f };
   int lda = 1;
   float B[] = { 0.16f, -0.747f };
   int ldb = 1;
   float C[] = { 0.576f };
   int ldc = 1;
   float C_expected[] = { 0.19641f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1625)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.0f;
   float A[] = { 0.259f, -0.334f };
   int lda = 1;
   float B[] = { -0.911f, -0.426f };
   int ldb = 1;
   float C[] = { 0.432f };
   int ldc = 1;
   float C_expected[] = { 0.056199f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1626)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.0f;
   float A[] = { -0.765f, 0.7f };
   int lda = 1;
   float B[] = { 0.487f, 0.768f };
   int ldb = 1;
   float C[] = { 0.836f };
   int ldc = 1;
   float C_expected[] = { -0.099027f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1627)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.0f;
   float A[] = { -0.584f, 0.056f };
   int lda = 2;
   float B[] = { 0.928f, -0.101f };
   int ldb = 2;
   float C[] = { -0.529f };
   int ldc = 1;
   float C_expected[] = { 0.328565f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1628)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha = -0.3f;
   float beta = 0.0f;
   float A[] = { 0.25f, 0.8f };
   int lda = 2;
   float B[] = { 0.489f, -0.642f };
   int ldb = 2;
   float C[] = { 0.322f };
   int ldc = 1;
   float C_expected[] = { 0.23481f };
   cblas_ssyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "ssyr2k(case 1629)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha = 0;
   double beta = 0;
   double A[] = { 0.591, 0.21 };
   int lda = 2;
   double B[] = { -0.718, -0.579 };
   int ldb = 2;
   double C[] = { -0.856 };
   int ldc = 1;
   double C_expected[] = { -0.0 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1630)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha = 0;
   double beta = 0;
   double A[] = { -0.971, -0.824 };
   int lda = 2;
   double B[] = { -0.227, 0.457 };
   int ldb = 2;
   double C[] = { 0.521 };
   int ldc = 1;
   double C_expected[] = { 0.0 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1631)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha = 0;
   double beta = 0;
   double A[] = { -0.274, 0.583 };
   int lda = 1;
   double B[] = { 0.668, -0.83 };
   int ldb = 1;
   double C[] = { 0.907 };
   int ldc = 1;
   double C_expected[] = { 0.0 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1632)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha = 0;
   double beta = 0;
   double A[] = { -0.512, -0.436 };
   int lda = 1;
   double B[] = { -0.443, -0.259 };
   int ldb = 1;
   double C[] = { -0.667 };
   int ldc = 1;
   double C_expected[] = { 0.0 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1633)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha = -0.3;
   double beta = 0;
   double A[] = { 0.741, -0.341 };
   int lda = 1;
   double B[] = { 0.743, -0.315 };
   int ldb = 1;
   double C[] = { -0.776 };
   int ldc = 1;
   double C_expected[] = { -0.3947868 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1634)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha = -0.3;
   double beta = 0;
   double A[] = { 0.03, 0.175 };
   int lda = 1;
   double B[] = { -0.832, 0.291 };
   int ldb = 1;
   double C[] = { 0.281 };
   int ldc = 1;
   double C_expected[] = { -0.015579 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1635)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha = -0.3;
   double beta = 0;
   double A[] = { 0.476, 0.496 };
   int lda = 2;
   double B[] = { -0.626, -0.159 };
   int ldb = 2;
   double C[] = { -0.964 };
   int ldc = 1;
   double C_expected[] = { 0.226104 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1636)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha = -0.3;
   double beta = 0;
   double A[] = { -0.489, 0.611 };
   int lda = 2;
   double B[] = { -0.285, -0.673 };
   int ldb = 2;
   double C[] = { -0.11 };
   int ldc = 1;
   double C_expected[] = { 0.1631028 };
   cblas_dsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dsyr2k(case 1637)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { 0.796f, 0.872f, -0.919f, 0.748f };
   int lda = 2;
   float B[] = { -0.945f, 0.915f, -0.252f, -0.276f };
   int ldb = 2;
   float C[] = { 0.07f, -0.957f };
   int ldc = 1;
   float C_expected[] = { -2.12843f, -0.054104f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1638) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1638) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { 0.984f, 0.526f, 0.284f, 0.806f };
   int lda = 2;
   float B[] = { -0.509f, -0.178f, 0.188f, -0.221f };
   int ldb = 2;
   float C[] = { -0.388f, 0.795f };
   int ldc = 1;
   float C_expected[] = { -0.43092f, -0.747044f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1639) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1639) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.16f, 0.628f, -0.06f, -0.645f };
   int lda = 1;
   float B[] = { 0.846f, 0.545f, 0.032f, 0.493f };
   int ldb = 1;
   float C[] = { -0.041f, -0.621f };
   int ldc = 1;
   float C_expected[] = { -0.26101f, 0.783636f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1640) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1640) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { -0.478f, -0.556f, 0.519f, 0.177f };
   int lda = 1;
   float B[] = { -0.946f, 0.423f, -0.859f, 0.736f };
   int ldb = 1;
   float C[] = { -0.54f, -0.035f };
   int ldc = 1;
   float C_expected[] = { 0.226066f, 1.05345f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1641) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1641) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { -0.582f, 0.09f, -0.176f, 0.784f };
   int lda = 1;
   float B[] = { 0.687f, -0.859f, 0.945f, 0.756f };
   int ldb = 1;
   float C[] = { -0.663f, -0.186f };
   int ldc = 1;
   float C_expected[] = { -0.663f, -0.186f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1642) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1642) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { 0.231f, -0.452f, -0.112f, -0.837f };
   int lda = 1;
   float B[] = { -0.258f, 0.464f, -0.224f, 0.893f };
   int ldb = 1;
   float C[] = { -0.448f, 0.046f };
   int ldc = 1;
   float C_expected[] = { -0.448f, 0.046f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1643) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1643) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { 0.115f, 0.178f, -0.193f, -0.491f };
   int lda = 2;
   float B[] = { 0.545f, -0.665f, 0.979f, -0.4f };
   int ldb = 2;
   float C[] = { 0.522f, 0.712f };
   int ldc = 1;
   float C_expected[] = { 0.522f, 0.712f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1644) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1644) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { -0.725f, -0.808f, -0.244f, 0.145f };
   int lda = 2;
   float B[] = { 0.447f, -0.413f, -0.226f, -0.585f };
   int ldb = 2;
   float C[] = { -0.531f, 0.227f };
   int ldc = 1;
   float C_expected[] = { -0.531f, 0.227f };
   cblas_csyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "csyr2k(case 1645) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "csyr2k(case 1645) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.337, -0.737, -0.993, 0.69 };
   int lda = 2;
   double B[] = { -0.39, -0.836, -0.32, 0.368 };
   int ldb = 2;
   double C[] = { 0.844, -0.763 };
   int ldc = 1;
   double C_expected[] = { 0.3494384, 0.5248712 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1646) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1646) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.386, -0.465, 0.719, -0.378 };
   int lda = 2;
   double B[] = { 0.099, -0.879, 0.864, 0.141 };
   int ldb = 2;
   double C[] = { -0.599, -0.47 };
   int ldc = 1;
   double C_expected[] = { 0.1664126, 0.5082238 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1647) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1647) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.914, 0.128, -0.262, -0.26 };
   int lda = 1;
   double B[] = { 0.431, 0.276, 0.75, 0.904 };
   int ldb = 1;
   double C[] = { 0.287, 0.537 };
   int ldc = 1;
   double C_expected[] = { -0.3532044, 0.0216788 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1648) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1648) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 111;
   int N = 1;
   int K = 2;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.618, 0.72, 0.392, -0.737 };
   int lda = 1;
   double B[] = { 0.783, 0.531, 0.375, 0.203 };
   int ldb = 1;
   double C[] = { 0.058, -0.116 };
   int ldc = 1;
   double C_expected[] = { -0.3837348, -0.2968344 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1649) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1649) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { -0.372, -0.735, -0.711, 0.051 };
   int lda = 1;
   double B[] = { 0.257, 0.097, 0.338, -0.484 };
   int ldb = 1;
   double C[] = { -0.142, -0.197 };
   int ldc = 1;
   double C_expected[] = { -0.414766, -0.676886 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1650) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1650) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { 0.1, -0.878, 0.28, -0.381 };
   int lda = 1;
   double B[] = { -0.208, 0.309, -0.276, 0.123 };
   int ldb = 1;
   double C[] = { 0.483, -0.541 };
   int ldc = 1;
   double C_expected[] = { -0.22324, -0.10083 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1651) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1651) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { -0.918, 0.515, -0.985, 0.067 };
   int lda = 2;
   double B[] = { -0.034, 0.408, 0.66, -0.945 };
   int ldb = 2;
   double C[] = { -0.063, -0.018 };
   int ldc = 1;
   double C_expected[] = { -1.228982, -1.549386 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1652) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1652) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int trans = 112;
   int N = 1;
   int K = 2;
   double alpha[2] = {0, 1};
   double beta[2] = {1, 0};
   double A[] = { 0.443, -0.009, -0.245, -0.008 };
   int lda = 2;
   double B[] = { 0.495, -0.239, -0.973, -0.032 };
   int ldb = 2;
   double C[] = { -0.85, -0.799 };
   int ldc = 1;
   double C_expected[] = { -0.660584, 0.111526 };
   cblas_zsyr2k(order, uplo, trans, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zsyr2k(case 1653) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zsyr2k(case 1653) imag");
     };
   };
  };

}


text:
/* vector/gsl_vector_long.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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_VECTOR_LONG_H__
#define __GSL_VECTOR_LONG_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_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_long.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  long *data;
  gsl_block_long *block;
  int owner;
}
gsl_vector_long;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_view;

typedef _gsl_vector_long_view gsl_vector_long_view;

typedef struct
{
  gsl_vector_long vector;
} _gsl_vector_long_const_view;

typedef const _gsl_vector_long_const_view gsl_vector_long_const_view;

/* Allocation */

GSL_FUN gsl_vector_long *gsl_vector_long_alloc (const size_t n);
GSL_FUN gsl_vector_long *gsl_vector_long_calloc (const size_t n);

GSL_FUN gsl_vector_long *gsl_vector_long_alloc_from_block (gsl_block_long * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_long *gsl_vector_long_alloc_from_vector (gsl_vector_long * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_long_free (gsl_vector_long * v);

/* Views */

GSL_FUN _gsl_vector_long_view
gsl_vector_long_view_array (long *v, size_t n);

GSL_FUN _gsl_vector_long_view
gsl_vector_long_view_array_with_stride (long *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_long_const_view
gsl_vector_long_const_view_array (const long *v, size_t n);

GSL_FUN _gsl_vector_long_const_view
gsl_vector_long_const_view_array_with_stride (const long *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_long_view
gsl_vector_long_subvector (gsl_vector_long *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_long_view
gsl_vector_long_subvector_with_stride (gsl_vector_long *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_long_const_view
gsl_vector_long_const_subvector (const gsl_vector_long *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_long_const_view
gsl_vector_long_const_subvector_with_stride (const gsl_vector_long *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_long_set_zero (gsl_vector_long * v);
GSL_FUN void gsl_vector_long_set_all (gsl_vector_long * v, long x);
GSL_FUN int gsl_vector_long_set_basis (gsl_vector_long * v, size_t i);

GSL_FUN int gsl_vector_long_fread (FILE * stream, gsl_vector_long * v);
GSL_FUN int gsl_vector_long_fwrite (FILE * stream, const gsl_vector_long * v);
GSL_FUN int gsl_vector_long_fscanf (FILE * stream, gsl_vector_long * v);
GSL_FUN int gsl_vector_long_fprintf (FILE * stream, const gsl_vector_long * v,
                              const char *format);

GSL_FUN int gsl_vector_long_memcpy (gsl_vector_long * dest, const gsl_vector_long * src);

GSL_FUN int gsl_vector_long_reverse (gsl_vector_long * v);

GSL_FUN int gsl_vector_long_swap (gsl_vector_long * v, gsl_vector_long * w);
GSL_FUN int gsl_vector_long_swap_elements (gsl_vector_long * v, const size_t i, const size_t j);

GSL_FUN long gsl_vector_long_max (const gsl_vector_long * v);
GSL_FUN long gsl_vector_long_min (const gsl_vector_long * v);
GSL_FUN void gsl_vector_long_minmax (const gsl_vector_long * v, long * min_out, long * max_out);

GSL_FUN size_t gsl_vector_long_max_index (const gsl_vector_long * v);
GSL_FUN size_t gsl_vector_long_min_index (const gsl_vector_long * v);
GSL_FUN void gsl_vector_long_minmax_index (const gsl_vector_long * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_long_add (gsl_vector_long * a, const gsl_vector_long * b);
GSL_FUN int gsl_vector_long_sub (gsl_vector_long * a, const gsl_vector_long * b);
GSL_FUN int gsl_vector_long_mul (gsl_vector_long * a, const gsl_vector_long * b);
GSL_FUN int gsl_vector_long_div (gsl_vector_long * a, const gsl_vector_long * b);
GSL_FUN int gsl_vector_long_scale (gsl_vector_long * a, const long x);
GSL_FUN int gsl_vector_long_add_constant (gsl_vector_long * a, const double x);
GSL_FUN int gsl_vector_long_axpby (const long alpha, const gsl_vector_long * x, const long beta, gsl_vector_long * y);
GSL_FUN long gsl_vector_long_sum (const gsl_vector_long * a);

GSL_FUN int gsl_vector_long_equal (const gsl_vector_long * u,
                            const gsl_vector_long * v);

GSL_FUN int gsl_vector_long_isnull (const gsl_vector_long * v);
GSL_FUN int gsl_vector_long_ispos (const gsl_vector_long * v);
GSL_FUN int gsl_vector_long_isneg (const gsl_vector_long * v);
GSL_FUN int gsl_vector_long_isnonneg (const gsl_vector_long * v);

GSL_FUN INLINE_DECL long gsl_vector_long_get (const gsl_vector_long * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x);
GSL_FUN INLINE_DECL long * gsl_vector_long_ptr (gsl_vector_long * v, const size_t i);
GSL_FUN INLINE_DECL const long * gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
long
gsl_vector_long_get (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_long_set (gsl_vector_long * v, const size_t i, long x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
long *
gsl_vector_long_ptr (gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (long *) (v->data + i * v->stride);
}

INLINE_FUN
const long *
gsl_vector_long_const_ptr (const gsl_vector_long * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const long *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_LONG_H__ */

ext: h

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
#pragma once

#include <bgfx/bgfx.h>
#include <napi/napi.h>
#include <gsl/gsl>

namespace Babylon
{
    class VertexBuffer final
    {
    public:
        VertexBuffer(gsl::span<uint8_t> bytes, bool dynamic);
        ~VertexBuffer();

        void Dispose();

        void Update(Napi::Env env, gsl::span<uint8_t> bytes);
        void CreateHandle(const bgfx::VertexLayout& layout);
        void PromoteToFloats(bgfx::AttribType::Enum attribType, uint32_t numElements, uint32_t byteOffset, uint32_t byteStride);
        void Set(bgfx::Encoder* encoder, uint8_t stream, uint32_t startVertex, uint32_t numVertices, bgfx::VertexLayoutHandle layoutHandle);

    private:
        std::vector<uint8_t> m_bytes{};
        bool m_dynamic{};

        union
        {
            bgfx::VertexBufferHandle m_handle{bgfx::kInvalidHandle};
            bgfx::DynamicVertexBufferHandle m_dynamicHandle;
        };

        bool m_disposed{};
    };
}

max_stars_repo_name: darraghjburke/BabylonNative

max_issues_repo_name: darraghjburke/BabylonNative

max_forks_repo_name: darraghjburke/BabylonNative


text:
#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <cconfigspace.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>

#define NUM_DISTRIBS 2
#define NUM_SAMPLES 10000

void test_create_mixture_distribution() {
	ccs_distribution_t      distrib = NULL, distribs[NUM_DISTRIBS], distribs_ret[NUM_DISTRIBS];
	ccs_result_t            err = CCS_SUCCESS;
	int32_t                 refcount;
	ccs_object_type_t       otype;
	ccs_distribution_type_t dtype;
	ccs_numeric_type_t      data_type;
	ccs_interval_t          interval;
	const size_t            num_distribs = NUM_DISTRIBS;
	ccs_float_t             weights[NUM_DISTRIBS];
	size_t                  num_distribs_ret;
	ccs_float_t             weights_ret[NUM_DISTRIBS];
	const ccs_float_t       epsilon = 1e-15;

	for(size_t i = 0; i < num_distribs; i++) {
		weights[i] = (double)(i+1);
	}

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(-10.0),
		CCSF(0.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(0.0),
		CCSF(10.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs + 1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_mixture_distribution(
		num_distribs,
		distribs,
		weights,
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_object_get_type(distrib, &otype);
	assert( err == CCS_SUCCESS );
	assert( otype == CCS_DISTRIBUTION );

	err = ccs_distribution_get_type(distrib, &dtype);
	assert( err == CCS_SUCCESS );
	assert( dtype == CCS_MIXTURE );

	err = ccs_distribution_get_data_types(distrib, &data_type);
	assert( err == CCS_SUCCESS );
	assert( data_type == CCS_NUM_FLOAT );

        err = ccs_distribution_get_bounds(distrib, &interval);
	assert( err == CCS_SUCCESS );
	assert( interval.type == CCS_NUM_FLOAT );
	assert( interval.lower.f == -10.0 );
	assert( interval.lower_included == CCS_TRUE );
	assert( interval.upper.f == 10.0 );
	assert( interval.upper_included == CCS_FALSE );

	err = ccs_mixture_distribution_get_num_distributions(distrib, &num_distribs_ret);
	assert( err == CCS_SUCCESS );
	assert( num_distribs_ret == num_distribs );

	err = ccs_mixture_distribution_get_distributions(distrib, num_distribs, distribs_ret, &num_distribs_ret);
	assert( err == CCS_SUCCESS );
	assert( num_distribs_ret == num_distribs );
	for (size_t i = 0; i < num_distribs; i++)
		assert( distribs_ret[i] == distribs[i] );

	err = ccs_mixture_distribution_get_weights(distrib, num_distribs, weights_ret, &num_distribs_ret);
	assert( err == CCS_SUCCESS );
	assert( num_distribs_ret == num_distribs );

	ccs_float_t inv_sum = 2.0 / (num_distribs * (num_distribs + 1));
	for (size_t i = 0; i < num_distribs; i++)
		assert( weights_ret[i] <= weights[i] * inv_sum + epsilon &&
		        weights_ret[i] >= weights[i] * inv_sum - epsilon );

	err = ccs_object_get_refcount(distrib, &refcount);
	assert( err == CCS_SUCCESS );
	assert( refcount == 1 );

	for (size_t i = 0; i < num_distribs; i++) {
		err = ccs_release_object(distribs[i]);
		assert( err == CCS_SUCCESS );
	}
	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
}

void test_mixture_distribution() {
	ccs_distribution_t distrib = NULL, distribs[NUM_DISTRIBS];
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_distribs = NUM_DISTRIBS;
	ccs_float_t        weights[NUM_DISTRIBS];
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES];

	for(size_t i = 0; i < num_distribs; i++) {
		weights[i] = 1.0;
	}

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(-9.0),
		CCSF(1.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(1.0),
		CCSF(11.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs + 1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_mixture_distribution(
		num_distribs,
		distribs,
		weights,
		&distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	double mean = gsl_stats_mean((double*)samples, 1, num_samples);
	assert( mean < 1.0 + 0.1 );
	assert( mean > 1.0 - 0.1 );

	for (size_t i = 0; i < num_distribs; i++) {
		err = ccs_release_object(distribs[i]);
		assert( err == CCS_SUCCESS );
	}
	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

void test_mixture_distribution_strided_samples() {
	ccs_distribution_t t_distrib = NULL, distrib = NULL, distribs[NUM_DISTRIBS];
	ccs_distribution_t t_distribs[NUM_DISTRIBS];
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_t_distribs = NUM_DISTRIBS;
	const size_t       num_distribs = NUM_DISTRIBS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples[NUM_SAMPLES*(NUM_DISTRIBS+1)];
	ccs_float_t        weights[] = { 1.0, 1.0 };

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(-5.0),
		CCSF(5.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_INTEGER,
		CCSI(-5),
		CCSI(5),
		CCS_LINEAR,
		CCSI(0),
		distribs + 1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_multivariate_distribution(
		num_distribs,
		distribs,
		&distrib);
	assert( err == CCS_SUCCESS );

	t_distribs[0] = distrib;
	t_distribs[1] = distrib;
	err = ccs_create_mixture_distribution(
		num_t_distribs,
		t_distribs,
		weights,
		&t_distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_strided_samples(distrib, rng, num_samples, NUM_DISTRIBS + 1, samples);
	assert( err == CCS_SUCCESS );

	for (size_t i = 0; i < num_samples; i++) {
		assert( samples[(NUM_DISTRIBS + 1)*i].f >= -5.0 );
		assert( samples[(NUM_DISTRIBS + 1)*i].f <   5.0 );
		assert( samples[(NUM_DISTRIBS + 1)*i + 1].i >= -5 );
		assert( samples[(NUM_DISTRIBS + 1)*i + 1].i <   5 );
	}
	double mean = gsl_stats_mean((double*)samples, NUM_DISTRIBS + 1, num_samples);
	assert( mean < 0.0 + 0.1 );
	assert( mean > 0.0 - 0.1 );

	for (size_t i = 0; i < num_distribs; i++) {
		err = ccs_release_object(distribs[i]);
		assert( err == CCS_SUCCESS );
	}
	err = ccs_release_object(t_distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

void test_mixture_distribution_soa_samples() {
	ccs_distribution_t t_distrib = NULL, distrib = NULL, distribs[NUM_DISTRIBS];
	ccs_distribution_t t_distribs[NUM_DISTRIBS];
	ccs_rng_t          rng = NULL;
	ccs_result_t       err = CCS_SUCCESS;
	const size_t       num_t_distribs = NUM_DISTRIBS;
	const size_t       num_distribs = NUM_DISTRIBS;
	const size_t       num_samples = NUM_SAMPLES;
	ccs_numeric_t      samples1[NUM_SAMPLES];
	ccs_numeric_t      samples2[NUM_SAMPLES];
	ccs_numeric_t     *samples[] = { samples1, samples2 };
	ccs_float_t        weights[] = { 1.0, 1.0 };

	err = ccs_rng_create(&rng);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_FLOAT,
		CCSF(-5.0),
		CCSF(5.0),
		CCS_LINEAR,
		CCSF(0.0),
		distribs);
	assert( err == CCS_SUCCESS );

	err = ccs_create_uniform_distribution(
		CCS_NUM_INTEGER,
		CCSI(-5),
		CCSI(5),
		CCS_LINEAR,
		CCSI(0),
		distribs + 1);
	assert( err == CCS_SUCCESS );

	err = ccs_create_multivariate_distribution(
		num_distribs,
		distribs,
		&distrib);
	assert( err == CCS_SUCCESS );

	t_distribs[0] = distrib;
	t_distribs[1] = distrib;
	err = ccs_create_mixture_distribution(
		num_t_distribs,
		t_distribs,
		weights,
		&t_distrib);
	assert( err == CCS_SUCCESS );

	err = ccs_distribution_soa_samples(distrib, rng, num_samples, samples);
	assert( err == CCS_SUCCESS );

	for (size_t i = 0; i < num_samples; i++) {
		assert( samples1[i].f >= -5.0 );
		assert( samples1[i].f <   5.0 );
		assert( samples2[i].i >= -5 );
		assert( samples2[i].i <   5 );
	}
	double mean = gsl_stats_mean((double*)samples1, 1, num_samples);
	assert( mean < 0.0 + 0.1 );
	assert( mean > 0.0 - 0.1 );

	for (size_t i = 0; i < num_distribs; i++) {
		err = ccs_release_object(distribs[i]);
		assert( err == CCS_SUCCESS );
	}
	err = ccs_release_object(t_distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(distrib);
	assert( err == CCS_SUCCESS );
	err = ccs_release_object(rng);
	assert( err == CCS_SUCCESS );
}

int main() {
	ccs_init();
	test_create_mixture_distribution();
	test_mixture_distribution();
	test_mixture_distribution_strided_samples();
	test_mixture_distribution_soa_samples();
	ccs_fini();
	return 0;
}

max_stars_repo_name: deephyper/CCS

max_issues_repo_name: deephyper/CCS

max_forks_repo_name: deephyper/CCS


text:
#pragma once

#include <functional>
#include "ProgressBar.h"
#include <gsl_util.h>

template<typename TCollection>
void runParallel(
	std::function<void(typename TCollection::reference)> processElement,
	TCollection& collection,
	int maxThreadCount)
{
	if (maxThreadCount < 1) {
		throw std::invalid_argument(fmt::format("maxThreadCount cannot be {}.", maxThreadCount));
	}

	if (maxThreadCount == 1) {
		// Process synchronously
		for (auto& element : collection) {
			processElement(element);
		}
		return;
	}

	using future_type = std::future<void>;

	std::mutex mutex;
	int currentThreadCount = 0;
	std::condition_variable elementFinished;
	future_type finishedElement;

	// Before exiting, wait for all running tasks to finish, but don't re-throw exceptions.
	// This only applies if one task already failed with an exception.
	auto finishRunning = gsl::finally([&]{
		std::unique_lock<std::mutex> lock(mutex);
		elementFinished.wait(lock, [&] { return currentThreadCount == 0; });
	});

	// Asyncronously run all elements
	for (auto it = collection.begin(); it != collection.end(); ++it) {
		// This variable will later hold the future, but can be value-captured right now
		auto future = std::make_shared<future_type>();

		// Notifies that an element is done processing
		auto notifyElementDone = [&, future] {
			std::lock_guard<std::mutex> lock(mutex);
			finishedElement = std::move(*future);
			--currentThreadCount;
			elementFinished.notify_one();
		};

		// Processes the current element, then notifies
		auto wrapperFunction = [processElement, &element = *it, notifyElementDone]() {
			auto done = gsl::finally(notifyElementDone);
			processElement(element);
		};

		// Asynchronously process element
		{
			std::lock_guard<std::mutex> lock(mutex);
			*future = std::async(std::launch::async, wrapperFunction);
			++currentThreadCount;
		}

		// Wait for threads to finish, if necessary
		{
			std::unique_lock<std::mutex> lock(mutex);
			int targetThreadCount = it == collection.end() ? 0 : maxThreadCount - 1;
			while (currentThreadCount > targetThreadCount) {
				elementFinished.wait(lock);
				if (finishedElement.valid()) {
					// Re-throw any exception
					finishedElement.get();
					finishedElement = future_type();
				}
			}
		}
	}

}

template<typename TCollection>
void runParallel(
	std::function<void(typename TCollection::reference, ProgressSink&)> processElement,
	TCollection& collection,
	int maxThreadCount,
	ProgressSink& progressSink,
	std::function<double(const typename TCollection::reference)> getElementProgressWeight = [](typename TCollection::reference) { return 1.0; })
{
	// Create a collection of wrapper functions that take care of progress handling
	ProgressMerger progressMerger(progressSink);
	std::vector<std::function<void()>> functions;
	for (auto& element : collection) {
		auto& elementProgressSink = progressMerger.addSink(getElementProgressWeight(element));
		functions.push_back([&]() { processElement(element, elementProgressSink); });
	}

	// Run wrapper function
	runParallel([&](std::function<void()> function) { function(); }, functions, maxThreadCount);
}

inline int getProcessorCoreCount() {
	int coreCount = std::thread::hardware_concurrency();

	// If the number of cores cannot be determined, use a reasonable default
	return coreCount != 0 ? coreCount : 4;
}

max_stars_repo_name: meshonline/rhubarb-lip-sync

max_issues_repo_name: meshonline/rhubarb-lip-sync

max_forks_repo_name: meshonline/rhubarb-lip-sync


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

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

/*
gsl_movstat_alloc()
  Allocate a workspace for moving window statistics.
The window around sample x_i is defined as:

W_i^{H,J} = {x_{i-H},...,x_i,...x_{i+J}}

The total window size is:

K = H + J + 1

Inputs: K - total samples in window (H = J = K / 2)

Return: pointer to workspace

Notes:
1) If K is even, it is rounded up to the next odd
*/

gsl_movstat_workspace *
gsl_movstat_alloc(const size_t K)
{
  const size_t H = K / 2;
  return gsl_movstat_alloc_with_size(0, H, H);
}

/*
gsl_movstat_alloc2()
  Allocate a workspace for moving window statistics.
The window around sample x_i is defined as:

W_i^{H,J} = {x_{i-H},...,x_i,...x_{i+J}}

The total window size is:

K = H + J + 1

Inputs: H - number of samples before current sample
        J - number of samples after current sample

Return: pointer to workspace
*/

gsl_movstat_workspace *
gsl_movstat_alloc2(const size_t H, const size_t J)
{
  return gsl_movstat_alloc_with_size(0, H, J);
}

/*
gsl_movstat_alloc_with_size()
  Allocate a workspace for moving window statistics with predefined workspace
size for accumulators. The window around sample x_i is defined as:

W_i^{H,J} = {x_{i-H},...,x_i,...x_{i+J}}

The total window size is:

K = H + J + 1

Inputs: accum_state_size - state size for accumulator (set to zero to use default value)
        H                - number of samples before current sample
        J                - number of samples after current sample

Return: pointer to workspace
*/

gsl_movstat_workspace *
gsl_movstat_alloc_with_size(const size_t accum_state_size, const size_t H, const size_t J)
{
  gsl_movstat_workspace *w;
  size_t state_size = accum_state_size;

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

  w->H = H;
  w->J = J;
  w->K = H + J + 1;

  if (state_size == 0)
    {
      /*
       * determine maximum number of bytes needed for the various accumulators;
       * the accumulators will all share the same workspace
       */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_mad->size)(w->K));    /* MAD accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_mean->size)(w->K));   /* mean/variance/sd accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_min->size)(w->K));    /* min/max accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_sum->size)(w->K));    /* sum accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_median->size)(w->K)); /* median accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_Qn->size)(w->K));     /* Q_n accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_qqr->size)(w->K));    /* QQR accumulator */
      state_size = GSL_MAX(state_size, (gsl_movstat_accum_Sn->size)(w->K));     /* S_n accumulator */
    }

  w->state = malloc(state_size);
  if (w->state == 0)
    {
      gsl_movstat_free(w);
      GSL_ERROR_NULL ("failed to allocate space for accumulator state", GSL_ENOMEM);
    }

  w->work = malloc(w->K * sizeof(double));
  if (w->work == 0)
    {
      gsl_movstat_free(w);
      GSL_ERROR_NULL ("failed to allocate space for work", GSL_ENOMEM);
    }

  w->state_size = state_size;

  return w;
}

void
gsl_movstat_free(gsl_movstat_workspace * w)
{
  if (w->work)
    free(w->work);

  if (w->state)
    free(w->state);

  free(w);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <php.h>
#include <Zend/zend_interfaces.h>
#include <Zend/zend_exceptions.h>
#include <ext/spl/spl_iterators.h>
#include <ext/spl/spl_exceptions.h>
#if _MSC_VER
#include <complex.h>
#define lapack_complex_float _Fcomplex
#define lapack_complex_double _Dcomplex
#endif
#include <lapacke.h>
#include <Interop/Polite/Math/Matrix.h>

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "php_rindow_openblas.h"

#ifndef MAX
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#endif
#ifndef MIN
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#endif

static zend_object_handlers rindow_openblas_lapack_object_handlers;

// destractor
static void php_rindow_openblas_lapack_free_object(zend_object* object)
{
    zend_object_std_dtor(object);
}

// constructor
static zend_object* php_rindow_openblas_lapack_create_object(zend_class_entry* class_type) /* {{{ */
{
    zend_object* intern = NULL;

    intern = (zend_object*)ecalloc(1, sizeof(zend_object) + zend_object_properties_size(class_type));

    zend_object_std_init(intern, class_type);
    object_properties_init(intern, class_type);

    intern->handlers = &rindow_openblas_lapack_object_handlers;

    return intern;
} /* }}} */

/* Method Rindow\OpenBLAS\Lapack::
    public function gesvd(
        int $matrix_layout,
        string $jobu,
        string $jobvt,
        int $m,
        int $n,
        Buffer $A,  int $offsetA,  int $ldA,
        Buffer $S,  int $offsetS,
        Buffer $U,  int $offsetU,  int $ldU,
        Buffer $VT, int $offsetVT, int $ldVT,
        Buffer $SuperB,  int $offsetSuperB
     ) : void
 {{{ */
static PHP_METHOD(Lapack, gesvd)
{
/*
lapack_int LAPACKE_dgesvd( int matrix_layout, char jobu, char jobvt,
                           lapack_int m, lapack_int n, double* a,
                           lapack_int lda, double* s, double* u, lapack_int ldu,
                           double* vt, lapack_int ldvt, double* superb );
*/
    php_interop_polite_math_matrix_linear_buffer_t* bufferA;
    php_interop_polite_math_matrix_linear_buffer_t* bufferS;
    php_interop_polite_math_matrix_linear_buffer_t* bufferU;
    php_interop_polite_math_matrix_linear_buffer_t* bufferVT;
    php_interop_polite_math_matrix_linear_buffer_t* bufferSuperB;
    zend_long matrix_layout;
    zend_long jobu;
    zend_long jobvt;
    zend_long m;
    zend_long n;
    zval* objA=NULL;
    zend_long offsetA;
    zend_long ldA;
    zval* objS=NULL;
    zend_long offsetS;
    zval* objU=NULL;
    zend_long offsetU;
    zend_long ldU;
    zval* objVT=NULL;
    zend_long offsetVT;
    zend_long ldVT;
    zval* objSuperB=NULL;
    zend_long offsetSuperB;
    int info;

    ZEND_PARSE_PARAMETERS_START_EX(ZEND_PARSE_PARAMS_THROW, 18, 18)
        Z_PARAM_LONG(matrix_layout)
        Z_PARAM_LONG(jobu)
        Z_PARAM_LONG(jobvt)
        Z_PARAM_LONG(m)
        Z_PARAM_LONG(n)

        Z_PARAM_OBJECT(objA) // Interop\Polite\Math\Matrix\LinearBuffer
        Z_PARAM_LONG(offsetA)
        Z_PARAM_LONG(ldA)
        Z_PARAM_OBJECT(objS) // Interop\Polite\Math\Matrix\LinearBuffer
        Z_PARAM_LONG(offsetS)

        Z_PARAM_OBJECT(objU) // Interop\Polite\Math\Matrix\LinearBuffer
        Z_PARAM_LONG(offsetU)
        Z_PARAM_LONG(ldU)
        Z_PARAM_OBJECT(objVT) // Interop\Polite\Math\Matrix\LinearBuffer
        Z_PARAM_LONG(offsetVT)

        Z_PARAM_LONG(ldVT)
        Z_PARAM_OBJECT(objSuperB) // Interop\Polite\Math\Matrix\LinearBuffer
        Z_PARAM_LONG(offsetSuperB)
    ZEND_PARSE_PARAMETERS_END();

    if(php_rindow_openblas_assert_shape_parameter(
        PHP_RINDOW_OPENBLAS_ASSERT_M, m)) {
        return;
    }
    if(php_rindow_openblas_assert_shape_parameter(
        PHP_RINDOW_OPENBLAS_ASSERT_N, n)) {
        return;
    }
    if( offsetS < 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "offsetS must be greater than zero or equal", 0);
        return;
    }
    if( offsetU < 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "offsetU must be greater than zero or equal", 0);
        return;
    }
    if( ldU <= 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "ldU must be greater than zero", 0);
        return;
    }
    if( offsetVT < 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "offsetVT must be greater than zero or equal", 0);
        return;
    }
    if( ldVT <= 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "ldVT must be greater than zero", 0);
        return;
    }
    if( offsetSuperB < 0 ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "offsetVT must be greater than zero or equal", 0);
        return;
    }
    // Check Buffer A
    bufferA = Z_INTEROP_POLITE_MATH_MATRIX_LINEAR_BUFFER_OBJ_P(objA);
    if(php_rindow_openblas_assert_buffer_type(bufferA,"a")) {
        return;
    }
    if(php_rindow_openblas_assert_matrix_buffer_spec(
        PHP_RINDOW_OPENBLAS_ASSERT_A, bufferA,m,n,offsetA,ldA)) {
        return;
    }

    // Check Buffer S
    bufferS = Z_INTEROP_POLITE_MATH_MATRIX_LINEAR_BUFFER_OBJ_P(objS);
    if(php_rindow_openblas_assert_buffer_type(bufferS,"s")) {
        return;
    }
    if( offsetS+MIN(m,n) > bufferS->size) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "BufferS size is too small", 0);
        return;
    }

    // Check Buffer U
    bufferU = Z_INTEROP_POLITE_MATH_MATRIX_LINEAR_BUFFER_OBJ_P(objU);
    if(php_rindow_openblas_assert_buffer_type(bufferU,"u")) {
        return;
    }
    if( offsetU+m*ldU > bufferU->size) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "BufferU size is too small", 0);
        return;
    }

    // Check Buffer VT
    bufferVT = Z_INTEROP_POLITE_MATH_MATRIX_LINEAR_BUFFER_OBJ_P(objVT);
    if(php_rindow_openblas_assert_buffer_type(bufferVT,"vt")) {
        return;
    }
    if( offsetVT+ldVT*n > bufferVT->size) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "BufferVT size is too small", 0);
        return;
    }

    // Check Buffer SuperB
    bufferSuperB = Z_INTEROP_POLITE_MATH_MATRIX_LINEAR_BUFFER_OBJ_P(objSuperB);
    if(php_rindow_openblas_assert_buffer_type(bufferSuperB,"b")) {
        return;
    }
    if( offsetSuperB+MIN(m,n)-1 > bufferSuperB->size) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "bufferSuperB size is too small", 0);
        return;
    }

    // Check Buffer A and B and C
    if(bufferA->dtype!=bufferS->dtype ||
       bufferA->dtype!=bufferU->dtype ||
       bufferA->dtype!=bufferVT->dtype ||
       bufferA->dtype!=bufferSuperB->dtype
    ) {
        zend_throw_exception(spl_ce_InvalidArgumentException, "Unmatch data type", 0);
        return;
    }

    switch (bufferA->dtype) {
        case php_interop_polite_math_matrix_dtype_float32:
            info = LAPACKE_sgesvd(
                (int)matrix_layout,
                (char)jobu,
                (char)jobvt,
                (lapack_int)m,(lapack_int)n,
                &(((float *)bufferA->data)[offsetA]), (lapack_int)ldA,
                &(((float *)bufferS->data)[offsetS]),
                &(((float *)bufferU->data)[offsetU]), (lapack_int)ldU,
                &(((float *)bufferVT->data)[offsetVT]), (lapack_int)ldVT,
                &(((float *)bufferSuperB->data)[offsetSuperB])
            );
            break;
        case php_interop_polite_math_matrix_dtype_float64:
            info = LAPACKE_dgesvd(
                (int)matrix_layout,
                (char)jobu,
                (char)jobvt,
                (lapack_int)m,(lapack_int)n,
                &(((double *)bufferA->data)[offsetA]), (lapack_int)ldA,
                &(((double *)bufferS->data)[offsetS]),
                &(((double *)bufferU->data)[offsetU]), (lapack_int)ldU,
                &(((double *)bufferVT->data)[offsetVT]), (lapack_int)ldVT,
                &(((double *)bufferSuperB->data)[offsetSuperB])
            );
            break;
        default:
            zend_throw_exception(spl_ce_RuntimeException, "Unsupported data type.", 0);
            return;
    }
    if( info == LAPACK_WORK_MEMORY_ERROR ) {
        zend_throw_exception(spl_ce_RuntimeException, "Not enough memory to allocate work array.", info);
        return;
    } else if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
        zend_throw_exception(spl_ce_RuntimeException, "Not enough memory to transpose matrix.", info);
        return;
    } else if( info < 0 ) {
        zend_throw_exception(spl_ce_RuntimeException, "Wrong parameter.", info);
        return;
    }
}
/* }}} */

ZEND_BEGIN_ARG_INFO_EX(ai_Lapack_gesvd, 0, 0, 18)
    ZEND_ARG_INFO(0, matrix_layout)
    ZEND_ARG_INFO(0, jobu)
    ZEND_ARG_INFO(0, jobvt)
    ZEND_ARG_INFO(0, m)
    ZEND_ARG_INFO(0, n)

    ZEND_ARG_OBJ_INFO(0, objA, Interop\\Polite\\Math\\Matrix\\LinearBuffer, 0)
    ZEND_ARG_INFO(0, offsetA)
    ZEND_ARG_INFO(0, ldA)
    ZEND_ARG_OBJ_INFO(0, objS, Interop\\Polite\\Math\\Matrix\\LinearBuffer, 0)
    ZEND_ARG_INFO(0, offsetS)

    ZEND_ARG_OBJ_INFO(0, objU, Interop\\Polite\\Math\\Matrix\\LinearBuffer, 0)
    ZEND_ARG_INFO(0, offsetU)
    ZEND_ARG_INFO(0, ldU)
    ZEND_ARG_OBJ_INFO(0, objVT, Interop\\Polite\\Math\\Matrix\\LinearBuffer, 0)
    ZEND_ARG_INFO(0, offsetVT)

    ZEND_ARG_INFO(0, ldVT)
    ZEND_ARG_OBJ_INFO(0, objSuperB, Interop\\Polite\\Math\\Matrix\\LinearBuffer, 0)
    ZEND_ARG_INFO(0, offsetSuperB)
ZEND_END_ARG_INFO()

ZEND_BEGIN_ARG_INFO_EX(ai_Lapack_void, 0, 0, 0)
ZEND_END_ARG_INFO()

/* {{{ Rindow\OpenBLAS\Lapack function entries */
static zend_function_entry php_rindow_openblas_lapack_me[] = {
    /* clang-format off */
    PHP_ME(Lapack, gesvd,  ai_Lapack_gesvd,  ZEND_ACC_PUBLIC)
    PHP_FE_END
    /* clang-format on */
};
/* }}} */

/* Class Rindow\OpenBLAS\Lapack {{{ */
static zend_class_entry* rindow_openblas_lapack_ce;

void php_rindow_openblas_lapack_init_ce(INIT_FUNC_ARGS)
{
    zend_class_entry ce;

    INIT_NS_CLASS_ENTRY(ce, "Rindow\\OpenBLAS", "Lapack", php_rindow_openblas_lapack_me);
    rindow_openblas_lapack_ce = zend_register_internal_class(&ce);
    rindow_openblas_lapack_ce->create_object = php_rindow_openblas_lapack_create_object;

    memcpy(&rindow_openblas_lapack_object_handlers, zend_get_std_object_handlers(), sizeof(zend_object_handlers));
    rindow_openblas_lapack_object_handlers.offset    = 0;
    rindow_openblas_lapack_object_handlers.free_obj  = php_rindow_openblas_lapack_free_object;
    rindow_openblas_lapack_object_handlers.clone_obj = NULL;

    //zend_class_implements(rindow_openblas_lapack_ce, 2, spl_ce_ArrayAccess, spl_ce_Countable);
}
/* }}} */

max_stars_repo_name: yuichiis/rindow-openblas

max_issues_repo_name: yuichiis/rindow-openblas

max_forks_repo_name: yuichiis/rindow-openblas


text:
// Rob Farber
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <nlopt.h>

#include "myFunc.h"

void writeParam(char *filename, int nParam, double *x)
{
  FILE *fn=fopen(filename,"w");
  if(!fn) {
    fprintf(stderr,"Cannot open %s\n",filename);
    exit(1);
  }

  uint32_t n=nParam; // ensure size is uint32_t
  fwrite(&n,sizeof(uint32_t), 1, fn);
  for(int i=0; i < nParam; i++) {
    float tmp=x[i];
    fwrite(&tmp,sizeof(float), 1, fn);
  }
  fclose(fn);
}
int main(int argc, char* argv[])
{
  double begintime = getTime();

  nlopt_opt opt;
  userData_t uData = {0};

  memset(&uData, 0, sizeof(userData_t));

  if(argc < 3) {
    fprintf(stderr,"Use: datafile paramFile\n");
    return -1;
  }
  init(argv[1],&uData);
  printf("myFunc %s\n", desc);
  printf("nExamples %d\n", uData.nExamples);
  printf("Number Parameters %d\n", N_PARAM);

  double endtime = getTime();

  printf(" memset and read file cost time: %g\n", endtime-begintime);

  begintime = getTime();

  opt = nlopt_create(NLOPT_LN_PRAXIS, N_PARAM); // algorithm and dimensionality
  // NOTE: alternative optimization methods ...
  //opt = nlopt_create(NLOPT_LN_NEWUOA, N_PARAM);
  //opt = nlopt_create(NLOPT_LN_COBYLA, N_PARAM);
  //opt = nlopt_create(NLOPT_LN_BOBYQA, N_PARAM);
  //opt = nlopt_create(NLOPT_LN_AUGLAG, N_PARAM);

  nlopt_set_min_objective(opt, objFunc, (void*) &uData);
/*
#if defined(MAX_RUNTIME)
  fprintf(stderr,"Warning: performing a quick %d second test!\n", MAX_RUNTIME);
  nlopt_set_maxtime(opt, MAX_RUNTIME); // Use for running quick tests
#else
  fprintf(stderr,"MAX runtime %d seconds!\n", 120*60);
  nlopt_set_maxtime(opt, (120. * 60.)); // maximum runtime in seconds
#endif
*/
  double minf; /* the minimum objective value, upon return */

  __declspec(align(64)) double x[N_PARAM];
  for(int i=0; i < N_PARAM; i++) x[i] = 0.1*(rand()/(double)RAND_MAX);

  endtime= getTime();

  printf(" serial part costs time : %g \n ", endtime-begintime);

  begintime = getTime();

  int nExamples=uData.nExamples*EXAMPLE_SIZE;
  float* example=uData.example;
  float* param=uData.param;
  {
    double startTime=getTime();
    int ret=nlopt_optimize(opt, x, &minf);
    printf("Optimization Time %g\n",getTime()-startTime);

    if (ret < 0) {
      printf("nlopt failed! ret %d\n", ret);
    } else {
      printf("found minimum %0.10g ret %d\n", minf,ret);
    }
  }

  endtime= getTime();

  printf(" parallel part costs time: %g \n", endtime-begintime);

  writeParam(argv[2],N_PARAM, x);
  fini(&uData);
  nlopt_destroy(opt);

  return 0;
}

max_stars_repo_name: puckbee/PhiBench

max_issues_repo_name: puckbee/PhiBench

max_forks_repo_name: puckbee/PhiBench


text: #include <boost/mpl/aux_/begin_end_impl.hpp>

ext: hpp

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/**
 * Copyright Soramitsu Co., Ltd. All Rights Reserved.
 * SPDX-License-Identifier: Apache-2.0
 */

#ifndef IROHA_TEST_FRAMEWORK_EXECUTOR_ITF_HELPER_HPP
#define IROHA_TEST_FRAMEWORK_EXECUTOR_ITF_HELPER_HPP

#include <memory>

#include <boost/mpl/at.hpp>
#include <boost/mpl/contains.hpp>
#include <boost/mpl/map.hpp>
#include "common/result.hpp"
#include "interfaces/commands/command.hpp"
#include "interfaces/queries/query.hpp"

#include "interfaces/query_responses/account_asset_response.hpp"
#include "interfaces/query_responses/account_detail_response.hpp"
#include "interfaces/query_responses/account_response.hpp"
#include "interfaces/query_responses/asset_response.hpp"
#include "interfaces/query_responses/block_error_response.hpp"
#include "interfaces/query_responses/block_query_response.hpp"
#include "interfaces/query_responses/block_response.hpp"
#include "interfaces/query_responses/engine_receipts_response.hpp"
#include "interfaces/query_responses/error_query_response.hpp"
#include "interfaces/query_responses/pending_transactions_page_response.hpp"
#include "interfaces/query_responses/query_response.hpp"
#include "interfaces/query_responses/role_permissions.hpp"
#include "interfaces/query_responses/roles_response.hpp"
#include "interfaces/query_responses/signatories_response.hpp"
#include "interfaces/query_responses/transactions_page_response.hpp"
#include "interfaces/query_responses/transactions_response.hpp"

namespace iroha {
  namespace integration_framework {
    namespace detail {

      /**
       * A type map from specific query to the corresponding response.
       *
       * Is used only to get a default expected response type for some query.
       * The queries may also return different responses, and if you expect a
       * different response type, you should either use the interface of
       * conversion functions to specify it, or extract it from the response
       * variant manually.
       */
      typedef boost::mpl::map<
          boost::mpl::pair<shared_model::interface::GetAccount,
                           shared_model::interface::AccountResponse>,
          boost::mpl::pair<shared_model::interface::GetSignatories,
                           shared_model::interface::SignatoriesResponse>,
          boost::mpl::pair<shared_model::interface::GetAccountTransactions,
                           shared_model::interface::TransactionsResponse>,
          boost::mpl::pair<shared_model::interface::GetAccountAssetTransactions,
                           shared_model::interface::TransactionsResponse>,
          boost::mpl::pair<shared_model::interface::GetTransactions,
                           shared_model::interface::TransactionsResponse>,
          boost::mpl::pair<shared_model::interface::GetAccountAssets,
                           shared_model::interface::AccountAssetResponse>,
          boost::mpl::pair<shared_model::interface::GetAccountDetail,
                           shared_model::interface::AccountDetailResponse>,
          boost::mpl::pair<shared_model::interface::GetRoles,
                           shared_model::interface::RolesResponse>,
          boost::mpl::pair<shared_model::interface::GetRolePermissions,
                           shared_model::interface::RolePermissionsResponse>,
          boost::mpl::pair<shared_model::interface::GetAssetInfo,
                           shared_model::interface::AssetResponse>,
          boost::mpl::pair<
              shared_model::interface::GetPendingTransactions,
              shared_model::interface::PendingTransactionsPageResponse>,
          boost::mpl::pair<shared_model::interface::GetBlock,
                           shared_model::interface::BlockResponse>,
          boost::mpl::pair<shared_model::interface::GetEngineReceipts,
                           shared_model::interface::EngineReceiptsResponse>>
          SpecificQueryResponses;

      /// true for specific commands
      template <typename T>
      constexpr bool isSpecificCommand = boost::mpl::contains<
          shared_model::interface::Command::CommandVariantType::types,
          const std::decay_t<T> &>::type::value;

      /// true for specific queries
      template <typename T>
      constexpr bool isSpecificQuery = boost::mpl::contains<
          shared_model::interface::Query::QueryVariantType::types,
          const std::decay_t<T> &>::type::value;

      template <typename T>
      std::enable_if_t<isSpecificQuery<T>, const T &> getInterfaceQueryRef(
          const T &query) {
        return query;
      }

      template <typename T, typename SpecificQuery = typename T::ModelType>
      std::enable_if_t<isSpecificQuery<SpecificQuery>, const SpecificQuery &>
      getInterfaceQueryRef(const T &query) {
        return static_cast<const SpecificQuery &>(query);
      }

      template <typename T>
      using InterfaceQuery = decltype(getInterfaceQueryRef(std::declval<T>()));

      /// response for specific query
      template <typename SpecificQuery,
                typename = std::enable_if_t<isSpecificQuery<SpecificQuery>>>
      using GetSpecificQueryResponse =
          typename boost::mpl::at<SpecificQueryResponses,
                                  std::decay_t<SpecificQuery>>::type;

      /**
       * Try to extract the given specific response from a general query
       * response.
       * @tparam SpecificQueryResponse The specific response type to extract.
       * @param query_result The general query response.
       * @return The specific response if the general one has it, the general
       * response if not.
       */
      template <typename SpecificQueryResponse>
      iroha::expected::Result<const SpecificQueryResponse &,
                              iroha::ametsuchi::QueryExecutorResult &>
      convertToSpecificQueryResponse(
          iroha::ametsuchi::QueryExecutorResult &query_result) {
        if (auto specific_query_response =
                boost::strict_get<const SpecificQueryResponse &>(
                    &query_result->get())) {
          return iroha::expected::makeValue<const SpecificQueryResponse &>(
              *specific_query_response);
        }
        return iroha::expected::makeError<
            iroha::ametsuchi::QueryExecutorResult &>(query_result);
      }
    }  // namespace detail
  }    // namespace integration_framework
}  // namespace iroha

#endif  // IROHA_TEST_FRAMEWORK_EXECUTOR_ITF_HELPER_HPP

max_stars_repo_name: Insafin/iroha

max_issues_repo_name: Insafin/iroha

max_forks_repo_name: Insafin/iroha


text:
#ifndef PERSPECTIVE_PROJECTION_HPP
#define PERSPECTIVE_PROJECTION_HPP

#include <Eigen/Core>

#include "log.hpp"
#include "projection.hpp"

class PerspectiveProjection : public Projection {
  public:
    PerspectiveProjection(double left, double right, double bottom, double top, double near, double far);

  private:
    void generateTransformation();
};

#endif // PERSPECTIVE_PROJECTION_HPP

max_stars_repo_name: matheusportela/black-hole-ray-tracer

max_issues_repo_name: matheusportela/black-hole-ray-tracer

max_forks_repo_name: matheusportela/black-hole-ray-tracer


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_SIMD_SDK_SIMD_EXTENSIONS_X86_XOP_HPP_INCLUDED
#define BOOST_SIMD_SDK_SIMD_EXTENSIONS_X86_XOP_HPP_INCLUDED

#if defined(__XOP__)
#  ifndef BOOST_SIMD_HAS_XOP_SUPPORT
#    define BOOST_SIMD_HAS_XOP_SUPPORT
#  endif
#elif defined(BOOST_SIMD_HAS_XOP_SUPPORT) && !defined(_MSC_VER)
#  undef BOOST_SIMD_HAS_XOP_SUPPORT
#endif

#if defined(BOOST_SIMD_HAS_XOP_SUPPORT) && !defined(BOOST_SIMD_HAS_FMA4_SUPPORT)
#  define BOOST_SIMD_HAS_FMA4_SUPPORT
#endif

#if !defined(BOOST_SIMD_DETECTED) && defined(BOOST_SIMD_HAS_XOP_SUPPORT)

////////////////////////////////////////////////////////////////////////////////
// XOP extensions flags
////////////////////////////////////////////////////////////////////////////////
#define BOOST_SIMD_DETECTED
#define BOOST_SIMD_XOP
#define BOOST_SIMD_SSE_FAMILY
#define BOOST_SIMD_STRING             "XOP"
#define BOOST_SIMD_STRING_LIST        "SSE2 SSE3 SSE4A SSSE3 SSE4_1 SSE4_2 AVX FMA4 XOP"
#define BOOST_SIMD_BYTES              32
#define BOOST_SIMD_BITS               256
#define BOOST_SIMD_CARDINALS          (2)(4)(8)(16)(32)
#define BOOST_SIMD_TAG_SEQ            (::boost::simd::tag::avx_)(::boost::simd::tag::sse_)
#define BOOST_SIMD_DEFAULT_EXTENSION  ::boost::simd::tag::avx_
#define BOOST_SIMD_DEFAULT_SITE       ::boost::simd::tag::xop_

// XOP header not standardized
#ifdef _MSC_VER
#include <intrin.h>
#else
#include <x86intrin.h>
#include <xopintrin.h>
#endif

#include <boost/simd/sdk/simd/extensions/meta/sse.hpp>
#include <boost/simd/sdk/simd/extensions/meta/avx.hpp>

#endif
#endif


text:
#define BOOST_TEST_MODULE "test_omp_external_distance_interaction"

#ifdef BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>
#else
#include <boost/test/included/unit_test.hpp>
#endif

#include <test/util/utility.hpp>

#include <mjolnir/math/math.hpp>
#include <mjolnir/core/BoundaryCondition.hpp>
#include <mjolnir/core/SimulatorTraits.hpp>
#include <mjolnir/core/AxisAlignedPlane.hpp>
#include <mjolnir/forcefield/external/LennardJonesWallPotential.hpp>
#include <mjolnir/omp/System.hpp>
#include <mjolnir/omp/UnlimitedGridCellList.hpp>
#include <mjolnir/omp/ExternalDistanceInteraction.hpp>
#include <mjolnir/util/make_unique.hpp>

BOOST_AUTO_TEST_CASE(omp_ExternalDistacne_calc_force)
{
    namespace test = mjolnir::test;

    constexpr double tol = 1e-8;
    mjolnir::LoggerManager::set_default_logger("test_omp_external_distance_interaction.log");

    using seq_traits_type  = mjolnir::SimulatorTraits      <double, mjolnir::UnlimitedBoundary>;
    using omp_traits_type  = mjolnir::OpenMPSimulatorTraits<double, mjolnir::UnlimitedBoundary>;

    using real_type        = typename omp_traits_type::real_type;
    using coordinate_type  = typename omp_traits_type::coordinate_type;
    using boundary_type    = typename omp_traits_type::boundary_type;
    using topology_type    = mjolnir::Topology;

    using potential_type   = mjolnir::LennardJonesWallPotential<real_type>;
    using parameter_type   = typename potential_type::parameter_type;

    using omp_system_type      = mjolnir::System<omp_traits_type>;
    using omp_shape_type       = mjolnir::AxisAlignedPlane<omp_traits_type, mjolnir::PositiveZDirection<omp_traits_type>>;
    using omp_interaction_type = mjolnir::ExternalDistanceInteraction<omp_traits_type, potential_type, omp_shape_type>;

    using seq_system_type      = mjolnir::System<seq_traits_type>;
    using seq_shape_type       = mjolnir::AxisAlignedPlane<seq_traits_type, mjolnir::PositiveZDirection<seq_traits_type>>;
    using seq_interaction_type = mjolnir::ExternalDistanceInteraction<seq_traits_type, potential_type, seq_shape_type>;

    const int max_number_of_threads = omp_get_max_threads();
    BOOST_TEST_WARN(max_number_of_threads > 2);
    BOOST_TEST_MESSAGE("maximum number of threads = " << max_number_of_threads);

    const std::size_t N_particle = 64;
    std::vector<std::pair<std::size_t, parameter_type>> parameters(N_particle);
    for(std::size_t i=0; i<N_particle; ++i)
    {
        parameters.emplace_back(i, parameter_type(1.0, 1.0));
    }
    const potential_type potential(2.5, parameters);
    topology_type topol(N_particle);
    topol.construct_molecules();

    std::mt19937 rng(123456789);

    for(int num_thread=1; num_thread<=max_number_of_threads; ++num_thread)
    {
        omp_set_num_threads(num_thread);
        BOOST_TEST_MESSAGE("maximum number of threads = " << omp_get_max_threads());

        omp_system_type omp_sys(N_particle, boundary_type{});
        seq_system_type seq_sys(N_particle, boundary_type{});
        test::clear_everything(omp_sys);
        test::clear_everything(seq_sys);

        for(std::size_t i=0; i<omp_sys.size(); ++i)
        {
            const auto i_x = i % 4;
            const auto i_y = i / 4;
            const auto i_z = i / 16;

            omp_sys.position(i) = mjolnir::math::make_coordinate<coordinate_type>(i_x*2.0, i_y*2.0, i_z*2.0);
        }

        test::apply_random_perturbation(omp_sys, rng, 0.1);

        // init sequential one with the same coordinates
        for(std::size_t i=0; i<seq_sys.size(); ++i)
        {
            seq_sys.mass(i)     = omp_sys.mass(i);
            seq_sys.position(i) = omp_sys.position(i);
            seq_sys.velocity(i) = omp_sys.velocity(i);
            seq_sys.force(i)    = omp_sys.force(i);
            seq_sys.name(i)     = omp_sys.name(i);
            seq_sys.group(i)    = omp_sys.group(i);
        }

        omp_shape_type omp_xyplane(0.0);
        seq_shape_type seq_xyplane(0.0);

        omp_xyplane.initialize(omp_sys, potential);
        seq_xyplane.initialize(seq_sys, potential);

        omp_interaction_type omp_interaction(std::move(omp_xyplane), potential_type(potential));
        seq_interaction_type seq_interaction(std::move(seq_xyplane), potential_type(potential));

        omp_interaction.initialize(omp_sys);
        seq_interaction.initialize(seq_sys);

        test::check_force_consistency              (omp_sys, omp_interaction, seq_sys, seq_interaction, tol);
        test::check_force_and_energy_consistency   (omp_sys, omp_interaction, seq_sys, seq_interaction, tol);
        test::check_energy_consistency             (omp_sys, omp_interaction, seq_sys, seq_interaction, tol);
    }
}

ext: cpp

max_stars_repo_name: ToruNiina/Mjolnir

max_issues_repo_name: Mjolnir-MD/Mjolnir

max_forks_repo_name: Mjolnir-MD/Mjolnir


text:
// iorate::variant
//
// Copyright iorate 2017.
// 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)

#define BOOST_TEST_DYN_LINK

#include <boost/test/unit_test.hpp>

#include <string>
#include <type_traits>
#include <utility>
#include "common.hpp"

BOOST_AUTO_TEST_CASE(copy_assignment)
{
    {
        iorate::variant<int, double> v1 = 42;
        iorate::variant<int, double> v2 = 65;
        v1 = v2;
        BOOST_TEST(iorate::get<0>(v1) == 65);
    }
    {
        iorate::variant<int, double> v1 = 42;
        iorate::variant<int, double> v2 = 12.0;
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
    }
    {
        iorate::variant<int, std::string> v1 = 42;
        iorate::variant<int, std::string> v2 = "foo";
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1) == "foo");
    }
    {
        iorate::variant<int, susp_copy_move> v1 = 42;
        iorate::variant<int, susp_copy_move> v2(std::in_place_index<1>, 65);
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1).value == 65);
    }
    {
        iorate::variant<int, bad_copy_move> v1 = 42;
        iorate::variant<int, bad_copy_move> v2(std::in_place_index<1>);
        BOOST_CHECK_THROW(v1 = v2, int);
        BOOST_TEST(iorate::get<0>(v1) == 42);
    }
    {
        iorate::variant<susp_copy_move, double> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, double> v2 = 12.0;
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
        force_backup(v1, 42);
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
    }
    {
        iorate::variant<susp_copy_move, susp_copy_move> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, susp_copy_move> v2(std::in_place_index<1>, 65);
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1).value == 65);
        force_backup(v1, 42);
        v1 = v2;
        BOOST_TEST(iorate::get<1>(v1).value == 65);
    }
    {
        iorate::variant<susp_copy_move, bad_copy_move> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, bad_copy_move> v2(std::in_place_index<1>);
        BOOST_CHECK_THROW(v1 = v2, int);
        BOOST_TEST(iorate::get<0>(v1).value == 42);
        BOOST_CHECK_THROW(v1 = v2, int);
        BOOST_TEST(iorate::get<0>(v1).value == 42);
    }
}

BOOST_AUTO_TEST_CASE(move_assignment)
{
    {
        static_assert(std::is_nothrow_move_assignable_v<iorate::variant<int, double>>);
        iorate::variant<int, double> v1 = 42;
        iorate::variant<int, double> v2 = 65;
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<0>(v1) == 65);
    }
    {
        iorate::variant<int, double> v1 = 42;
        iorate::variant<int, double> v2 = 12.0;
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
    }
    {
        static_assert(std::is_nothrow_move_assignable_v<iorate::variant<int, std::string>>);
        iorate::variant<int, std::string> v1 = 42;
        iorate::variant<int, std::string> v2 = "foo";
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1) == "foo");
    }
    {
        static_assert(!std::is_nothrow_move_assignable_v<iorate::variant<int, susp_copy_move>>);
        iorate::variant<int, susp_copy_move> v1 = 42;
        iorate::variant<int, susp_copy_move> v2(std::in_place_index<1>, 65);
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1).value == 65);
        BOOST_TEST(iorate::get<1>(v2).value == 0);
    }
    {
        static_assert(!std::is_nothrow_move_assignable_v<iorate::variant<int, bad_copy_move>>);
        iorate::variant<int, bad_copy_move> v1 = 42;
        iorate::variant<int, bad_copy_move> v2(std::in_place_index<1>);
        BOOST_CHECK_THROW(v1 = std::move(v2), int);
        BOOST_TEST(iorate::get<0>(v1) == 42);
    }
    {
        static_assert(!std::is_nothrow_move_assignable_v<iorate::variant<susp_copy_move, double>>);
        iorate::variant<susp_copy_move, double> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, double> v2 = 12.0;
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
        force_backup(v1, 42);
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1) == 12.0);
    }
    {
        static_assert(!std::is_nothrow_move_assignable_v<iorate::variant<susp_copy_move, susp_copy_move>>);
        iorate::variant<susp_copy_move, susp_copy_move> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, susp_copy_move> v2(std::in_place_index<1>, 65);
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1).value == 65);
        BOOST_TEST(iorate::get<1>(v2).value == 0);
        force_backup(v1, 42);
        v2.emplace<1>(65);
        v1 = std::move(v2);
        BOOST_TEST(iorate::get<1>(v1).value == 65);
        BOOST_TEST(iorate::get<1>(v2).value == 0);
    }
    {
        static_assert(!std::is_nothrow_move_assignable_v<iorate::variant<susp_copy_move, bad_copy_move>>);
        iorate::variant<susp_copy_move, bad_copy_move> v1(std::in_place_index<0>, 42);
        iorate::variant<susp_copy_move, bad_copy_move> v2(std::in_place_index<1>);
        BOOST_CHECK_THROW(v1 = std::move(v2), int);
        BOOST_TEST(iorate::get<0>(v1).value == 42);
        BOOST_CHECK_THROW(v1 = std::move(v2), int);
        BOOST_TEST(iorate::get<0>(v1).value == 42);
    }
}

BOOST_AUTO_TEST_CASE(assignment)
{
    {
        iorate::variant<int, double> v = 42;
        static_assert(noexcept(v = 65));
        v = 65;
        BOOST_TEST(iorate::get<0>(v) == 65);
    }
    {
        iorate::variant<int, double> v = 42;
        static_assert(noexcept(v = 12.0));
        v = 12.0;
        BOOST_TEST(iorate::get<1>(v) == 12.0);
    }
    {
        iorate::variant<int, std::string> v = 42;
        static_assert(!noexcept(v = "foo"));
        v = "foo";
        BOOST_TEST(iorate::get<1>(v) == "foo");
    }
    {
        iorate::variant<int, bad_copy_move> v = 42;
        static_assert(!noexcept(v = bad_copy_move()));
        BOOST_CHECK_THROW(v = bad_copy_move(), int);
        BOOST_TEST(iorate::get<0>(v) == 42);
    }
    {
        iorate::variant<susp_copy_move, double> v(std::in_place_index<0>, 42);
        static_assert(noexcept(v = 12.0));
        v = 12.0;
        BOOST_TEST(iorate::get<1>(v) == 12.0);
        force_backup(v, 42);
        v = 12.0;
        BOOST_TEST(iorate::get<1>(v) == 12.0);
    }
    {
        iorate::variant<susp_copy_move, std::string> v(std::in_place_index<0>, 42);
        static_assert(!noexcept(v = "foo"));
        v = "foo";
        BOOST_TEST(iorate::get<1>(v) == "foo");
        force_backup(v, 42);
        v = "foo";
        BOOST_TEST(iorate::get<1>(v) == "foo");
    }
    {
        iorate::variant<susp_copy_move, bad_copy_move> v(std::in_place_index<0>, 42);
        static_assert(!noexcept(v = bad_copy_move()));
        BOOST_CHECK_THROW(v = bad_copy_move(), int);
        BOOST_TEST(iorate::get<0>(v).value == 42);
        BOOST_CHECK_THROW(v = bad_copy_move(), int);
        BOOST_TEST(iorate::get<0>(v).value == 42);
    }
}

max_stars_repo_name: iorate/variant

max_issues_repo_name: iorate/variant

max_forks_repo_name: iorate/variant