text:
#include <cblas.h>
#include <lapacke.h>

#include "math.h"
#include "minimal_opencv.h"
#include "stdbool.h"
#include "stdio.h"
#include "string.h"

#include <limits.h>
#include <stdarg.h>

#include "linmath.h"

#ifdef _WIN32
#define SURVIVE_LOCAL_ONLY
#include <malloc.h>
#define alloca _alloca
#else
#define SURVIVE_LOCAL_ONLY __attribute__((visibility("hidden")))
#endif

SURVIVE_LOCAL_ONLY int cvRound(float f) { return roundf(f); }
#define CV_Error(code, msg) assert(0 && msg); // cv::error( code, msg, CV_Func, __FILE__, __LINE__ )

const int DECOMP_SVD = 1;
const int DECOMP_LU = 2;

void print_mat(const CvMat *M);

static size_t mat_size_bytes(const CvMat *mat) { return (size_t)CV_ELEM_SIZE(mat->type) * mat->cols * mat->rows; }

SURVIVE_LOCAL_ONLY void cvCopy(const CvMat *srcarr, CvMat *dstarr, const CvMat *mask) {
	assert(mask == 0 && "This isn't implemented yet");
	assert(srcarr->rows == dstarr->rows);
	assert(srcarr->cols == dstarr->cols);
	assert(dstarr->type == srcarr->type);
	memcpy(dstarr->data.db, srcarr->data.db, mat_size_bytes(srcarr));
}

SURVIVE_LOCAL_ONLY void cvGEMM(const CvMat *src1, const CvMat *src2, double alpha, const CvMat *src3, double beta,
							   CvMat *dst, int tABC) {

	int rows1 = (tABC & CV_GEMM_A_T) ? src1->cols : src1->rows;
	int cols1 = (tABC & CV_GEMM_A_T) ? src1->rows : src1->cols;

	int rows2 = (tABC & CV_GEMM_B_T) ? src2->cols : src2->rows;
	int cols2 = (tABC & CV_GEMM_B_T) ? src2->rows : src2->cols;

	if (src3) {
		int rows3 = (tABC & CV_GEMM_C_T) ? src3->cols : src3->rows;
		int cols3 = (tABC & CV_GEMM_C_T) ? src3->rows : src3->cols;
		assert(rows3 == dst->rows);
		assert(cols3 == dst->cols);
	}

	// assert(src3 == 0 || beta != 0);
	assert(cols1 == rows2);
    assert(rows1 == dst->rows);
    assert(cols2 == dst->cols);

	lapack_int lda = src1->cols;
	lapack_int ldb = src2->cols;

	if (src3)
		cvCopy(src3, dst, 0);
	else
		beta = 0;

	assert(dst->data.db != src1->data.db);
	assert(dst->data.db != src2->data.db);

	cblas_dgemm(CblasRowMajor, (tABC & CV_GEMM_A_T) ? CblasTrans : CblasNoTrans,
				(tABC & CV_GEMM_B_T) ? CblasTrans : CblasNoTrans, dst->rows, dst->cols, cols1, alpha, src1->data.db,
				lda, src2->data.db, ldb, beta, dst->data.db, dst->cols);
}

SURVIVE_LOCAL_ONLY void cvMulTransposed(const CvMat *src, CvMat *dst, int order, const CvMat *delta, double scale) {
	lapack_int rows = src->rows;
	lapack_int cols = src->cols;

	lapack_int drows = dst->rows;
	assert(drows == cols);
	assert(order == 1 ? (dst->cols == src->cols) : (dst->cols == src->rows));
	assert(delta == 0 && "This isn't implemented yet");
	double beta = 0;

	bool isAT = order == 1;
	bool isBT = !isAT;

	lapack_int dstCols = dst->cols;

	cblas_dgemm(CblasRowMajor, isAT ? CblasTrans : CblasNoTrans, isBT ? CblasTrans : CblasNoTrans, cols, dstCols, rows,
				scale, src->data.db, cols, src->data.db, cols, beta, dst->data.db, dstCols);
}

SURVIVE_LOCAL_ONLY void *cvAlloc(size_t size) { return malloc(size); }

static void icvCheckHuge(CvMat *arr) {
	if ((int64_t)arr->step * arr->rows > INT_MAX)
		arr->type &= ~CV_MAT_CONT_FLAG;
}

SURVIVE_LOCAL_ONLY CvMat *cvCreateMatHeader(int rows, int cols, int type) {
	type = CV_MAT_TYPE(type);

	assert(!(rows < 0 || cols < 0));

	int min_step = CV_ELEM_SIZE(type);
	assert(!(min_step <= 0));
	min_step *= cols;

	CvMat *arr = (CvMat *)cvAlloc(sizeof(*arr));

	arr->step = min_step;
	arr->type = CV_MAT_MAGIC_VAL | type | CV_MAT_CONT_FLAG;
	arr->rows = rows;
	arr->cols = cols;
	arr->data.ptr = 0;
	arr->refcount = 0;
	arr->hdr_refcount = 1;

	icvCheckHuge(arr);
	return arr;
}

/* the alignment of all the allocated buffers */
#define CV_MALLOC_ALIGN 16

/* IEEE754 constants and macros */
#define CV_TOGGLE_FLT(x) ((x) ^ ((int)(x) < 0 ? 0x7fffffff : 0))
#define CV_TOGGLE_DBL(x) ((x) ^ ((int64)(x) < 0 ? CV_BIG_INT(0x7fffffffffffffff) : 0))

#define CV_DbgAssert assert

static inline void *cvAlignPtr(const void *ptr, int align) {
	CV_DbgAssert((align & (align - 1)) == 0);
	return (void *)(((size_t)ptr + align - 1) & ~(size_t)(align - 1));
}

SURVIVE_LOCAL_ONLY void cvCreateData(CvMat *arr) {
	if (CV_IS_MAT_HDR_Z(arr)) {
		size_t step, total_size;
		CvMat *mat = (CvMat *)arr;
		step = mat->step;

		if (mat->rows == 0 || mat->cols == 0)
			return;

		if (mat->data.ptr != 0)
			CV_Error(CV_StsError, "Data is already allocated");

		if (step == 0)
			step = CV_ELEM_SIZE(mat->type) * mat->cols;

		int64_t _total_size = (int64_t)step * mat->rows + sizeof(int) + CV_MALLOC_ALIGN;
		total_size = (size_t)_total_size;
		if (_total_size != (int64_t)total_size)
			CV_Error(CV_StsNoMem, "Too big buffer is allocated");
		mat->refcount = (int *)cvAlloc((size_t)total_size);
		mat->data.ptr = (uchar *)cvAlignPtr(mat->refcount + 1, CV_MALLOC_ALIGN);
		*mat->refcount = 1;
	} else if (CV_IS_MATND_HDR(arr)) {
		assert("There is no support for ND types");
	} else
		CV_Error(CV_StsBadArg, "unrecognized or unsupported array type");
}

SURVIVE_LOCAL_ONLY CvMat *cvCreateMat(int height, int width, int type) {
	CvMat *arr = cvCreateMatHeader(height, width, type);
	cvCreateData(arr);

	return arr;
}

SURVIVE_LOCAL_ONLY double cvInvert(const CvMat *srcarr, CvMat *dstarr, int method) {
	lapack_int inf;
	lapack_int rows = srcarr->rows;
	lapack_int cols = srcarr->cols;
	lapack_int lda = srcarr->cols;

	cvCopy(srcarr, dstarr, 0);
	double *a = dstarr->data.db;

#ifdef DEBUG_PRINT
	printf("a: \n");
	print_mat(srcarr);
#endif
	if (method == DECOMP_LU) {
		lapack_int *ipiv = malloc(sizeof(lapack_int) * MIN(srcarr->rows, srcarr->cols));
		inf = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, rows, cols, a, lda, ipiv);
		assert(inf == 0);

		inf = LAPACKE_dgetri(LAPACK_ROW_MAJOR, rows, a, lda, ipiv);
		assert(inf >= 0);
		if (inf > 0) {
			printf("Warning: Singular matrix: \n");
			// print_mat(srcarr);
		}

		free(ipiv);

	} else if (method == DECOMP_SVD) {
		// TODO: There is no way this needs this many allocations,
		// but in my defense I was very tired when I wrote this code
		CvMat *w = cvCreateMat(1, MIN(dstarr->rows, dstarr->cols), dstarr->type);
		CvMat *u = cvCreateMat(dstarr->cols, dstarr->cols, dstarr->type);
		CvMat *v = cvCreateMat(dstarr->rows, dstarr->rows, dstarr->type);
		CvMat *um = cvCreateMat(w->cols, w->cols, w->type);

		cvSVD(dstarr, w, u, v, 0);

		cvSetZero(um);
		for (int i = 0; i < w->cols; i++) {
			cvmSet(um, i, i, 1. / w->data.db[i]);
		}

		CvMat *tmp = cvCreateMat(dstarr->cols, dstarr->rows, dstarr->type);
		cvGEMM(v, um, 1, 0, 0, tmp, CV_GEMM_A_T);
		cvGEMM(tmp, u, 1, 0, 0, dstarr, CV_GEMM_B_T);

		cvReleaseMat(&tmp);
		cvReleaseMat(&w);
		cvReleaseMat(&u);
		cvReleaseMat(&v);
		cvReleaseMat(&um);
	}
	return 0;
}

SURVIVE_LOCAL_ONLY CvMat *cvCloneMat(const CvMat *mat) {
	CvMat *rtn = cvCreateMat(mat->rows, mat->cols, mat->type);
	cvCopy(mat, rtn, 0);
	return rtn;
}

SURVIVE_LOCAL_ONLY int cvSolve(const CvMat *Aarr, const CvMat *xarr, CvMat *Barr, int method) {
	lapack_int inf;
	lapack_int arows = Aarr->rows;
	lapack_int acols = Aarr->cols;
	lapack_int xcols = xarr->cols;
	lapack_int xrows = xarr->rows;
	lapack_int lda = acols; // Aarr->step / sizeof(double);
	lapack_int type = CV_MAT_TYPE(Aarr->type);

	if (method == DECOMP_LU) {
		assert(Aarr->cols == Barr->rows);
		assert(xarr->rows == Aarr->rows);
		assert(Barr->cols == xarr->cols);
		assert(type == CV_MAT_TYPE(Barr->type) && (type == CV_32F || type == CV_64F));

		cvCopy(xarr, Barr, 0);
		CvMat *a_ws = cvCloneMat(Aarr);

		lapack_int brows = Barr->rows;
		lapack_int bcols = Barr->cols;
		lapack_int ldb = bcols; // Barr->step / sizeof(double);

		lapack_int *ipiv = malloc(sizeof(lapack_int) * MIN(Aarr->rows, Aarr->cols));

		inf = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, arows, acols, a_ws->data.db, lda, ipiv);
		assert(inf >= 0);
		if (inf > 0) {
			printf("Warning: Singular matrix: \n");
			// print_mat(a_ws);
		}

#ifdef DEBUG_PRINT
		printf("Solve A * x = B:\n");
		print_mat(a_ws);
		print_mat(Barr);
#endif

		inf =
			LAPACKE_dgetrs(LAPACK_ROW_MAJOR, CblasNoTrans, arows, bcols, a_ws->data.db, lda, ipiv, Barr->data.db, ldb);
		assert(inf == 0);

		free(ipiv);
		cvReleaseMat(&a_ws);
	} else if (method == DECOMP_SVD) {

#ifdef DEBUG_PRINT
		printf("Solve |b - A * x|:\n");
		print_mat(Aarr);
		print_mat(xarr);
#endif
		bool xLargerThanB = xarr->rows > acols;
		CvMat *xCpy = 0;
		if (xLargerThanB) {
			xCpy = cvCloneMat(xarr);
		} else {
			xCpy = Barr;
			memcpy(Barr->data.db, xarr->data.db, mat_size_bytes(xarr));
		}

		CvMat *aCpy = cvCloneMat(Aarr);

		double *S = malloc(sizeof(double) * MIN(arows, acols));
		double rcond = -1;
		lapack_int *rank = malloc(sizeof(lapack_int) * MIN(arows, acols));
		lapack_int inf = LAPACKE_dgelss(LAPACK_ROW_MAJOR, arows, acols, xcols, aCpy->data.db, acols, xCpy->data.db,
										xcols, S, rcond, rank);
		free(rank);
		free(S);

		assert(Barr->rows == acols);
		assert(Barr->cols == xCpy->cols);

		if (xLargerThanB) {
			xCpy->rows = acols;
			cvCopy(xCpy, Barr, 0);
			cvReleaseMat(&xCpy);
		}

		cvReleaseMat(&aCpy);
#ifdef DEBUG_PRINT
		print_mat(Barr);
#endif
		assert(inf == 0);
	}
	return 0;
}

SURVIVE_LOCAL_ONLY void cvTranspose(const CvMat *M, CvMat *dst) {
	bool inPlace = M == dst || M->data.db == dst->data.db;
	double *src = M->data.db;

	CvMat *tmp = 0;
	if (inPlace) {
		tmp = cvCloneMat(dst);
		src = tmp->data.db;
	} else {
	  assert(M->rows == dst->cols);
	  assert(M->cols == dst->rows);
	}

	for (unsigned i = 0; i < M->rows; i++) {
		for (unsigned j = 0; j < M->cols; j++) {
			dst->data.db[j * M->rows + i] = src[i * M->cols + j];
		}
	}

	if (inPlace) {
		cvReleaseMat(&tmp);
	}
}

SURVIVE_LOCAL_ONLY void cvSVD(CvMat *aarr, CvMat *warr, CvMat *uarr, CvMat *varr, int flags) {
	char jobu = 'A';
	char jobvt = 'A';

	lapack_int inf;

	if ((flags & CV_SVD_MODIFY_A) == 0) {
		aarr = cvCloneMat(aarr);
	}

	if (uarr == 0)
		jobu = 'N';
	if (varr == 0)
		jobvt = 'N';

	double *pw, *pu, *pv;
	lapack_int arows = aarr->rows, acols = aarr->cols;

	pw = warr ? warr->data.db : (double *)alloca(sizeof(double) * arows * acols);
	pu = uarr ? uarr->data.db : (double *)alloca(sizeof(double) * arows * arows);
	pv = varr ? varr->data.db : (double *)alloca(sizeof(double) * acols * acols);

	lapack_int ulda = uarr ? uarr->cols : acols;
	lapack_int plda = varr ? varr->cols : acols;

	double *superb = malloc(sizeof(double) * MIN(arows, acols));
	inf = LAPACKE_dgesvd(LAPACK_ROW_MAJOR, jobu, jobvt, arows, acols, aarr->data.db, acols, pw, pu, ulda, pv, plda,
						 superb);

	free(superb);

	switch (inf) {
	case -6:
		assert(false && "matrix has NaNs");
		break;
	case 0:
		break;
	default:
		assert(inf == 0);
	}

	if (uarr && (flags & CV_SVD_U_T)) {
		cvTranspose(uarr, uarr);
	}

	if (varr && (flags & CV_SVD_V_T) == 0) {
		cvTranspose(varr, varr);
	}

	if ((flags & CV_SVD_MODIFY_A) == 0) {
		cvReleaseMat(&aarr);
	}
}

SURVIVE_LOCAL_ONLY void cvSetZero(CvMat *arr) {
	for (int i = 0; i < arr->rows; i++)
		for (int j = 0; j < arr->cols; j++)
			arr->data.db[i * arr->cols + j] = 0;
}
SURVIVE_LOCAL_ONLY void cvSetIdentity(CvMat *arr) {
	for (int i = 0; i < arr->rows; i++)
		for (int j = 0; j < arr->cols; j++)
			arr->data.db[i * arr->cols + j] = i == j;
}

SURVIVE_LOCAL_ONLY void cvReleaseMat(CvMat **mat) {
	assert(*(*mat)->refcount == 1);
	free((*mat)->refcount);
	free(*mat);
	*mat = 0;
}

SURVIVE_LOCAL_ONLY double cvDet(const CvMat *M) {
	assert(M->rows == M->cols);
	assert(M->rows <= 3 && "cvDet unimplemented for matrices >3");
	assert(CV_64F == CV_MAT_TYPE(M->type) && "cvDet unimplemented for float");
	double *m = M->data.db;

	switch (M->rows) {
	case 1:
		return m[0];
	case 2: {
		return m[0] * m[3] - m[1] * m[2];
	}
	case 3: {
		double m00 = m[0], m01 = m[1], m02 = m[2], m10 = m[3], m11 = m[4], m12 = m[5], m20 = m[6], m21 = m[7],
			   m22 = m[8];

		return m00 * (m11 * m22 - m12 * m21) - m01 * (m10 * m22 - m12 * m20) + m02 * (m10 * m21 - m11 * m20);
	}
	default:
		abort();
	}
}

max_issues_repo_name: ChristophHaag/libsurvive

max_forks_repo_name: ChristophHaag/libsurvive


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* The sum of N samples from an exponential distribution gives an
   Erlang distribution

   p(x) dx = x^(n-1) exp (-x/a) / ((n-1)!a^n) dx

   for x = 0 ... +infty */

double
gsl_ran_erlang (const gsl_rng * r, const double a, const double n)
{
  return gsl_ran_gamma (r, n, a);
}

double
gsl_ran_erlang_pdf (const double x, const double a, const double n)
{
  if (x <= 0)
    {
      return 0 ;
    }
  else
    {
      double p;
      double lngamma = gsl_sf_lngamma (n);

      p = exp ((n - 1) * log (x/a) - x/a - lngamma) / a;
      return p;
    }
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// --------------------------------------------------------------------------
//                   OpenMS -- Open-Source Mass Spectrometry
// --------------------------------------------------------------------------
// Copyright The OpenMS Team -- Eberhard Karls University Tuebingen,
// ETH Zurich, and Freie Universitaet Berlin 2002-2013.
//
// This software is released under a three-clause BSD license:
//  * Redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//  * Redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//  * Neither the name of any author or any participating institution
//    may be used to endorse or promote products derived from this software
//    without specific prior written permission.
// For a full list of authors, refer to the file AUTHORS.
// --------------------------------------------------------------------------
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL ANY OF THE AUTHORS OR THE CONTRIBUTING
// INSTITUTIONS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
// OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
// WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
// OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
// ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// --------------------------------------------------------------------------
// $Maintainer: Erhan Kenar $
// $Authors: $
// --------------------------------------------------------------------------

#ifndef OPENMS_DATASTRUCTURES_MATRIX_H
#define OPENMS_DATASTRUCTURES_MATRIX_H

#include <OpenMS/CONCEPT/Macros.h>

#include <cmath> // pow()
#include <iomanip>
#include <vector>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>

namespace OpenMS
{

  /**
    @brief A two-dimensional matrix.  Similar to std::vector, but uses a binary
    operator(,) for element access.

    Think of it as a random access container.  You can also generate gray
    scale images.  This data structure is not designed to be used for linear algebra,
    but rather a simple two-dimensional array.

    The following member functions of the base class std::vector<ValueType>
    can also be used:

begin

end

rbegin

rend

front

back

assign

empty

size

capacity

max_size

         @ingroup Datastructures
  */
  template <typename Value>
  class Matrix :
    protected std::vector<Value>
  {
protected:
    typedef std::vector<Value> Base;

public:

    ///@name STL compliance type definitions
    //@{
    typedef Base container_type;

    typedef typename Base::difference_type difference_type;
    typedef typename Base::size_type size_type;

    typedef typename Base::const_iterator const_iterator;
    typedef typename Base::const_reverse_iterator const_reverse_iterator;
    typedef typename Base::iterator iterator;
    typedef typename Base::reverse_iterator reverse_iterator;

    typedef typename Base::const_reference const_reference;
    typedef typename Base::pointer pointer;
    typedef typename Base::reference reference;
    typedef typename Base::value_type value_type;

    typedef typename Base::allocator_type allocator_type;
    //@}

    ///@name OpenMS compliance type definitions
    //@{
    typedef Base ContainerType;
    typedef difference_type DifferenceType;
    typedef size_type SizeType;

    typedef const_iterator ConstIterator;
    typedef const_reverse_iterator ConstReverseIterator;
    typedef iterator Iterator;
    typedef reverse_iterator ReverseIterator;

    typedef const_reference ConstReference;
    typedef pointer Pointer;
    typedef reference Reference;
    typedef value_type ValueType;

    typedef allocator_type AllocatorType;
    //@}

    ///@name Constructors, assignment, and destructor
    //@{
    Matrix() :
      Base(),
      rows_(0),
      cols_(0)
    {}

    Matrix(const SizeType rows, const SizeType cols, ValueType value = ValueType()) :
      Base(rows * cols, value),
      rows_(rows),
      cols_(cols)
    {}

    Matrix(const Matrix& source) :
      Base(source),
      rows_(source.rows_),
      cols_(source.cols_)
    {}

    Matrix& operator=(const Matrix& rhs)
    {
      Base::operator=(rhs);
      rows_ = rhs.rows_;
      cols_ = rhs.cols_;
      return *this;
    }

    ~Matrix() {}
    //@}

    ///@name Accessors
    //@{
    const_reference operator()(size_type const i, size_type const j) const
    {
      return getValue(i, j);
    }

    reference operator()(size_type const i, size_type const j)
    {
      return getValue(i, j);
    }

    const_reference getValue(size_type const i, size_type const j) const
    {
      return Base::operator[](index(i, j));
    }

    reference getValue(size_type const i, size_type const j)
    {
      return Base::operator[](index(i, j));
    }

    void setValue(size_type const i, size_type const j, value_type value)
    {
      Base::operator[](index(i, j)) = value;
    }

    /// Return the i-th row of the matrix as a vector.
    container_type row(size_type const i) const
    {
#ifdef OPENMS_DEBUG
      if (i >= rows_) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, i, rows_);
#endif
      container_type values(cols_);
      for (size_type j = 0; j < cols_; j++)
      {
        values[j] = Base::operator[](index(i, j));
      }
      return values;
    }

    /// Return the i-th column of the matrix as a vector.
    container_type col(size_type const i) const
    {
#ifdef OPENMS_DEBUG
      if (i >= cols_) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, i, cols_);
#endif
      container_type values(rows_);
      for (size_type j = 0; j < rows_; j++)
      {
        values[j] = Base::operator[](index(j, i));
      }
      return values;
    }

    //@}

    /**
      @name Pure access declarations

      These make begin(), end() etc. from container_type accessible.
    */
    //@{
public:

    using Base::begin;
    using Base::end;
    using Base::rbegin;
    using Base::rend;

    using Base::front;
    using Base::back;
    using Base::assign;

    using Base::empty;
    using Base::size;

    using Base::capacity;
    using Base::max_size;

    //@}

    void clear()
    {
      Base::clear();
      rows_ = 0;
      cols_ = 0;
    }

    void resize(size_type i, size_type j, value_type value = value_type())
    {
      rows_ = i;
      cols_ = j;
      Base::resize(rows_ * cols_, value);
    }

    void resize(std::pair<Size, Size> const& size_pair, value_type value = value_type())
    {
      rows_ = size_pair.first;
      cols_ = size_pair.second;
      Base::resize(rows_ * cols_, value);
    }

    /// Number of rows
    SizeType rows() const
    {
      return rows_;
    }

    /// Number of columns
    SizeType cols() const
    {
      return cols_;
    }

    std::pair<Size, Size> sizePair() const
    {
      return std::pair<Size, Size>(rows_, cols_);
    }

    /**
      @brief Calculate the index into the underlying vector from row and column.
      Note that Matrix uses the (row,column) lexicographic ordering for indexing.
    */
    SizeType const index(SizeType row, SizeType col) const
    {
#ifdef OPENMS_DEBUG
      if (row >= rows_) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, row, rows_);
      if (col >= cols_) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, col, cols_);
#endif
      return row * cols_ + col;
    }

    /**
      @brief Calculate the row and column from an index into the underlying vector.
      Note that Matrix uses the (row,column) lexicographic ordering for indexing.
    */
    std::pair<Size, Size> const indexPair(Size index) const
    {
#ifdef OPENMS_DEBUG
      if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, index, size() - 1);
#endif
      return std::pair<SizeType, SizeType>(index / cols_, index % cols_);
    }

    /**
      @brief Calculate the column from an index into the underlying vector.
      Note that Matrix uses the (row,column) lexicographic ordering for indexing.
    */
    SizeType colIndex(SizeType index) const
    {
#ifdef OPENMS_DEBUG
      if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, index, size() - 1);
#endif
      return index % cols_;
    }

    /**
      @brief Calculate the row from an index into the underlying vector.
      Note that Matrix uses the (row,column) lexicographic ordering for indexing.
    */
    SizeType rowIndex(SizeType index) const
    {
#ifdef OPENMS_DEBUG
      if (index >= size()) throw Exception::IndexOverflow(__FILE__, __LINE__, __PRETTY_FUNCTION__, index, size() - 1);
#endif

      return index / cols_;
    }

    /**
      @brief Equality comparator.

      If matrices have different row or column numbers, throws a precondition exception.
    */
    bool operator==(Matrix const& rhs) const
    {
      OPENMS_PRECONDITION(cols_ == rhs.cols_,
                          "Matrices have different row sizes.");
      OPENMS_PRECONDITION(rows_ == rhs.rows_,
                          "Matrices have different column sizes.");
      return static_cast<typename Matrix<Value>::Base const&>(*this) == static_cast<typename Matrix<Value>::Base const&>(rhs);
    }

    /**
      @brief Less-than comparator.  Comparison is done lexicographically: first by row, then by column.

      If matrices have different row or column numbers, throws a precondition exception.
    */
    bool operator<(Matrix const& rhs) const
    {
      OPENMS_PRECONDITION(cols_ == rhs.cols_,
                          "Matrices have different row sizes.");
      OPENMS_PRECONDITION(rows_ == rhs.rows_,
                          "Matrices have different column sizes.");
      return static_cast<typename Matrix<Value>::Base const&>(*this) < static_cast<typename Matrix<Value>::Base const&>(rhs);
    }

    /// set matrix to 2D arrays values
    template <int ROWS, int COLS>
    void setMatrix(const ValueType matrix[ROWS][COLS])
    {
      resize(ROWS, COLS);
      for (SizeType i = 0; i < this->rows_; ++i)
      {
        for (SizeType j = 0; j < this->cols_; ++j)
        {
          setValue(i, j, matrix[i][j]);
        }
      }
    }

    /**
      @brief create gsl_matrix*

      allocate and return an equivalent GSL matrix
      @note Works only for Matrix<double> and Matrix<float>
      @note Clean up the gsl_matrix using gsl_matrix_free (gsl_matrix * m)
    */
    gsl_matrix* toGslMatrix()
    {
      gsl_matrix* m_ptr = gsl_matrix_alloc(rows_, cols_);

      for (size_type i = 0; i < this->rows_; ++i)
      {
        for (size_type j = 0; j < this->cols_; ++j)
        {
          gsl_matrix_set(m_ptr, i, j, (double) (*this)(i, j));
        }
      }

      return m_ptr;
    }

protected:

    ///@name Data members
    //@{
    /// Number of rows (height of a column)
    SizeType rows_;
    /// Number of columns (width of a row)
    SizeType cols_;
    //@}

  }; // class Matrix

  // template<> OPENMS_DLLAPI gsl_matrix* Matrix<double>::toGslMatrix();
  // template<> OPENMS_DLLAPI gsl_matrix* Matrix<float>::toGslMatrix();

  /**
    @brief Print the contents to a stream.

    @relatesalso Matrix
  */
  template <typename Value>
  std::ostream& operator<<(std::ostream& os, const Matrix<Value>& matrix)
  {
    typedef typename Matrix<Value>::size_type size_type;
    for (size_type i = 0; i < matrix.rows(); ++i)
    {
      for (size_type j = 0; j < matrix.cols(); ++j)
      {
        os << std::setprecision(6) << std::setw(6) << matrix(i, j) << ' ';
      }
      os << std::endl;
    }
    return os;
  }

} // namespace OpenMS

#endif // OPENMS_DATASTRUCTURES_MATRIX_H

ext: h

max_stars_repo_name: open-ms/all-svn-branches

max_issues_repo_name: open-ms/all-svn-branches

max_forks_repo_name: open-ms/all-svn-branches


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_statistics.h>
#include <assert.h>
#include <float.h>

extern "C" {
#include <quadmath.h>
}

#ifndef N
#define N 32
#endif

#define K N/2

#ifndef S
#define S K/2
#endif

#ifndef SEED
#define SEED 1
#endif

#ifndef FB
#define FB 64
#endif

#if(FB==64)
#define FT double
#elif(FB==80)
#define FT long double
#endif

int
main(int argc, char *argv[])
{
  assert(argc == 4);

  int i, j;
  char *cinname = argv[1];
  char *sinname = argv[2];
  char *outname = argv[3];
  FILE *cinfile = fopen(cinname, "r");
  FILE *infile = fopen(sinname, "r");
  FILE *outfile = fopen(outname, "w");

  srand(SEED);

  int index[K];
  for (i=0; i<K; i++)
    index[i]=i;
  for (i=0; i<K-1; i++){
    j = rand()%(K-1);
    int tmp = index[i];
    index[i] = index[j];
    index[j] = tmp;
  }

  FT data[K], w[K];
  FT wtss;

  // read concrete input file
  float cinput[N]; i=0;
  while(i<N && fscanf(cinfile, "%f", &cinput[i])>0){
      i++;
  }
  assert(i>=N);
  fclose(cinfile);

  for (i = 0; i < K-S; i++){
    w[index[i]]=cinput[index[i]];
  }

  for (i = 0 ; i <K ; i++){
    data[i] = cinput[i+K];
  }

  // read symbolic input data
  float sinput[S];
  for (i = 0 ; i < S ; i++) {
    __float128 in_data;
    fread(&in_data, sizeof(__float128), 1, infile);
    sinput[i] = (float) in_data;
  }
  fclose(infile);

  for (i = 0 ; i<S; i++){
    w[index[K-1-i]] = sinput[i] / FLT_MAX * 100;
  }

  // library call
  #if(FB==64)
    wtss     = gsl_stats_wtss(w, 1, data, 1, K);
  #elif(FB==80)
    wtss     = gsl_stats_long_double_wtss(w, 1, data, 1, K);
  #else
    exit(-1);
  #endif

  __float128 out_data;
  out_data = (__float128) wtss;
  fwrite(&out_data, sizeof(__float128), 1, outfile);
  fclose(outfile);

  return 0;
}

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \addtogroup ioda_cxx_attribute
 *
 * @{
 * \file Attribute_Creator.h
 * \brief Flywheel creation of ioda::Attribute. Used by ioda::Has_Attributes and
 * ioda::VariableCreationParameters.
 */

#include <gsl/gsl-lite.hpp>
#include <iostream>
#include <memory>
#include <string>
#include <vector>

#include "ioda/Attributes/Attribute.h"
#include "ioda/Attributes/Has_Attributes.h"
#include "ioda/Misc/Dimensions.h"
#include "ioda/Types/Type.h"
#include "ioda/defs.h"

namespace ioda {
namespace detail {
/// \brief Flywheel creation of ioda::Attribute.
/// \ingroup ioda_cxx_attribute
class IODA_DL Attribute_Creator_Base {
protected:
  std::string name_;

public:
  virtual ~Attribute_Creator_Base();
  virtual void apply(Has_Attributes& obj) const = 0;
  Attribute_Creator_Base(const std::string& name);
};
}  // namespace detail

/// \brief Flywheel creation of ioda::Attribute.
/// \ingroup ioda_cxx_attribute
template <class DataType>
class Attribute_Creator : public detail::Attribute_Creator_Base {
private:
  ::std::vector<Dimensions_t> dimensions_;
  ::std::vector<DataType> data_;

public:
  virtual ~Attribute_Creator() {}
  void apply(Has_Attributes& obj) const override { obj.add<DataType>(name_, data_, dimensions_); }

  /// \details Totally taking advantage of vector constructors.
  template <class DataInput, class DimensionInput>
  Attribute_Creator(const std::string& name, DataInput data, DimensionInput dimensions)
      : Attribute_Creator_Base(name), data_(data.begin(), data.end()), dimensions_(dimensions) {}

  template <class DimensionInput>
  Attribute_Creator(const std::string& name, DimensionInput dimensions)
      : Attribute_Creator_Base(name), dimensions_(dimensions) {}

  template <class DataInput_junk_param = DataType>
  void write(const ::gsl::span<const DataType>& data) {
    data_ = ::std::vector<DataType>(data.begin(), data.end());
  }
};

/// \brief Flywheel creation of ioda::Attribute objects.
/// \ingroup ioda_cxx_attribute
/// This is needed because you might want to make the same Attribute in multiple places.
class IODA_DL Attribute_Creator_Store : public detail::CanAddAttributes<Attribute_Creator_Store> {
  std::vector<std::shared_ptr<detail::Attribute_Creator_Base>> atts_;

public:
  Attribute_Creator_Store();
  virtual ~Attribute_Creator_Store();

  void apply(Has_Attributes& obj) const;

  /// @name Convenience functions for adding attributes
  /// @{
  /// \see CanAddAttributes
  /// \see CanAddAttributes::add

  template <class DataType>
  Attribute_Creator_Store& create(const std::string& attrname, ::gsl::span<const DataType> data,
                                  ::std::initializer_list<Dimensions_t> dimensions) {
    atts_.push_back(std::make_shared<Attribute_Creator<DataType>>(attrname, data, dimensions));
    return *this;
  }

  template <class DataType>
  Attribute_Creator_Store& create(const std::string& attrname,
                                  ::std::initializer_list<DataType> data,
                                  ::std::initializer_list<Dimensions_t> dimensions) {
    atts_.push_back(std::make_shared<Attribute_Creator<DataType>>(attrname, data, dimensions));
    return *this;
  }

  template <class DataType2>
  struct AttWrapper {
    std::shared_ptr<Attribute_Creator<DataType2>> inner;
    AttWrapper(std::shared_ptr<Attribute_Creator<DataType2>> d) : inner(d) {}
    template <class DataInput_junk_param = DataType2>
    void write(const ::gsl::span<const DataType2>& data) {
      inner->write(data);
    }
  };

  template <class DataType>
  AttWrapper<DataType> create(const std::string& attrname, ::std::vector<Dimensions_t> dimensions) {
    auto res = std::make_shared<Attribute_Creator<DataType>>(attrname, dimensions);
    atts_.push_back(res);
    return AttWrapper<DataType>{res};
  }

  /// @}
};

}  // namespace ioda

/// @} // End Doxygen block

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


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

/**
  This is a lagged Fibonacci generator which supposedly excellent
  statistical properties (I do not concur)

  I got it from the net and translated into C.

* ======================================================================
* NIST Guide to Available Math Software.
* Fullsource for module UNI from package CMLIB.
* Retrieved from CAMSUN on Tue Oct  8 14:04:10 1996.
* ======================================================================

C***BEGIN PROLOGUE  UNI
C***DATE WRITTEN   810915
C***REVISION DATE  830805
C***CATEGORY NO.  L6A21
C***KEYWORDS  RANDOM NUMBERS, UNIFORM RANDOM NUMBERS
C***AUTHOR    BLUE, JAMES, SCIENTIFIC COMPUTING DIVISION, NBS
C             KAHANER, DAVID, SCIENTIFIC COMPUTING DIVISION, NBS
C             MARSAGLIA, GEORGE, COMPUTER SCIENCE DEPT., WASH STATE UNIV
C
C***PURPOSE  THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON [0,1
C             AND CAN BE USED ON ANY COMPUTER WITH WHICH ALLOWS INTEGERS
C             AT LEAST AS LARGE AS 32767.
C***DESCRIPTION
C
C       THIS ROUTINE GENERATES QUASI UNIFORM RANDOM NUMBERS ON THE INTER
C       [0,1).  IT CAN BE USED WITH ANY COMPUTER WHICH ALLOWS
C       INTEGERS AT LEAST AS LARGE AS 32767.
C
C
C   USE
C       FIRST TIME....
C                   Z = UNI(JD)
C                     HERE JD IS ANY  N O N - Z E R O  INTEGER.
C                     THIS CAUSES INITIALIZATION OF THE PROGRAM
C                     AND THE FIRST RANDOM NUMBER TO BE RETURNED AS Z.
C       SUBSEQUENT TIMES...
C                   Z = UNI(0)
C                     CAUSES THE NEXT RANDOM NUMBER TO BE RETURNED AS Z.
C
C
C..................................................................
C   NOTE: USERS WHO WISH TO TRANSPORT THIS PROGRAM FROM ONE COMPUTER
C         TO ANOTHER SHOULD READ THE FOLLOWING INFORMATION.....
C
C   MACHINE DEPENDENCIES...
C      MDIG = A LOWER BOUND ON THE NUMBER OF BINARY DIGITS AVAILABLE
C              FOR REPRESENTING INTEGERS, INCLUDING THE SIGN BIT.
C              THIS VALUE MUST BE AT LEAST 16, BUT MAY BE INCREASED
C              IN LINE WITH REMARK A BELOW.
C
C   REMARKS...
C     A. THIS PROGRAM CAN BE USED IN TWO WAYS:
C        (1) TO OBTAIN REPEATABLE RESULTS ON DIFFERENT COMPUTERS,
C            SET 'MDIG' TO THE SMALLEST OF ITS VALUES ON EACH, OR,
C        (2) TO ALLOW THE LONGEST SEQUENCE OF RANDOM NUMBERS TO BE
C            GENERATED WITHOUT CYCLING (REPEATING) SET 'MDIG' TO THE
C            LARGEST POSSIBLE VALUE.
C     B. THE SEQUENCE OF NUMBERS GENERATED DEPENDS ON THE INITIAL
C          INPUT 'JD' AS WELL AS THE VALUE OF 'MDIG'.
C          IF MDIG=16 ONE SHOULD FIND THAT
   Editors Note: set the seed using 152 in order to get uni(305)
   -jt
C            THE FIRST EVALUATION
C              Z=UNI(305) GIVES Z=.027832881...
C            THE SECOND EVALUATION
C              Z=UNI(0) GIVES   Z=.56102176...
C            THE THIRD EVALUATION
C              Z=UNI(0) GIVES   Z=.41456343...
C            THE THOUSANDTH EVALUATION
C              Z=UNI(0) GIVES   Z=.19797357...
C
C***REFERENCES  MARSAGLIA G., "COMMENTS ON THE PERFECT UNIFORM RANDOM
C                 NUMBER GENERATOR", UNPUBLISHED NOTES, WASH S. U.
C***ROUTINES CALLED  I1MACH,XERROR
C***END PROLOGUE  UNI

  **/

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

static inline unsigned long int uni_get (void *vstate);
static double uni_get_double (void *vstate);
static void uni_set (void *state, unsigned long int s);

static const unsigned int MDIG = 16;    /* Machine digits in int */
static const unsigned int m1 = 32767;   /* 2^(MDIG-1) - 1 */
static const unsigned int m2 = 256;     /* 2^(MDIG/2) */

typedef struct
  {
    int i, j;
    unsigned long m[17];
  }
uni_state_t;

static inline unsigned long
uni_get (void *vstate)
{
  uni_state_t *state = (uni_state_t *) vstate;
  const int i = state->i;
  const int j = state->j;

  /* important k not be unsigned */
  long k = state->m[i] - state->m[j];

  if (k < 0)
    k += m1;
  state->m[j] = k;

  if (i == 0)
    {
      state->i = 16;
    }
  else
    {
      (state->i)--;
    }

  if (j == 0)
    {
      state->j = 16;
    }
  else
    {
      (state->j)--;
    }

  return k;
}

static double
uni_get_double (void *vstate)
{
  return uni_get (vstate) / 32767.0 ;
}

static void
uni_set (void *vstate, unsigned long int s)
{
  unsigned int i, seed, k0, k1, j0, j1;

  uni_state_t *state = (uni_state_t *) vstate;

  /* For this routine, the seeding is very elaborate! */
  /* A flaw in this approach is that seeds 1,2 give exactly the
     same random number sequence!  */

  s = 2 * s + 1;        /* enforce seed be odd */
  seed = (s < m1 ? s : m1);     /* seed should be less than m1 */

  k0 = 9069 % m2;
  k1 = 9069 / m2;
  j0 = seed % m2;
  j1 = seed / m2;

  for (i = 0; i < 17; ++i)
    {
      seed = j0 * k0;
      j1 = (seed / m2 + j0 * k1 + j1 * k0) % (m2 / 2);
      j0 = seed % m2;
      state->m[i] = j0 + m2 * j1;
    }
  state->i = 4;
  state->j = 16;

  return;
}

static const gsl_rng_type uni_type =
{"uni",                         /* name */
 32766,                         /* RAND_MAX */
 0,                             /* RAND_MIN */
 sizeof (uni_state_t),
 &uni_set,
 &uni_get,
 &uni_get_double};

const gsl_rng_type *gsl_rng_uni = &uni_type;


text:
#pragma once

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

#include "defines.h"
#include "matrix_vector_ops.h"

#define NUM_STATES 1
#define NUM_INPUTS 0

// A-stability test of RK method:
// y' = ky
int testDynamics(gsl_vector *dy, double t, const gsl_vector *y, const gsl_vector *u);

int setupTestDynamics();
int teardownTestDynamics();

gsl_vector *testFeedback(gsl_vector *yd, gsl_vector *y);

int setupTestFeedback();
int teardownTestFeedback();

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim


text:
#if MKL_FOUND
#include <mkl.h>
#elif BLAS_FOUND
  #if WASM_COMPATIBLE_BLAS
    #include "3rd_party/onnxjs/src/wasm-ops/gemm.h"
  #else
    #include <cblas.h>
  #endif // WASM_COMPATIBLE_BLAS
#endif

inline void sgemm(bool transA,
                  bool transB,
                  int rows_a,
                  int rows_b,
                  int width,
                  float alpha,
                  float* a,
                  int lda,
                  float* b,
                  int ldb,
                  float beta,
                  float* c,
                  int ldc) {
#if BLAS_FOUND
    #if WASM_COMPATIBLE_BLAS
        gemm_f32_imp(transA, transB, rows_a, rows_b, width, alpha, a, b, beta, c);
    #else
        cblas_sgemm(CblasRowMajor,
                    transA ? CblasTrans : CblasNoTrans,
                    transB ? CblasTrans : CblasNoTrans,
                    rows_a,
                    rows_b,
                    width,
                    alpha,
                    a,
                    lda,
                    b,
                    ldb,
                    beta,
                    c,
                    ldc);
    #endif
#else
    transA; transB; rows_a; rows_b; width; alpha; a; lda; b; ldb; beta; c; ldc; // make compiler happy
    ABORT("Marian must be compiled with a BLAS library");
#endif
}

max_stars_repo_name: felipesantosk/marian-dev

max_issues_repo_name: felipesantosk/marian-dev

max_forks_repo_name: felipesantosk/marian-dev


text:
//#include <gsl/gsl>

#include "../Problem15/ipv4Base.h"

class ipv4Extend : public ipv4Base
{
public:
    ipv4Extend() : ipv4Base()
    {
    };

    ipv4Extend( std::string   aInput ) : ipv4Base( aInput )
    {
    };

    ~ipv4Extend()
    {
    };

    ipv4Extend& operator= ( const ipv4Extend & aOther )
    {
        if ( this != & aOther )
        {
            mIPAddress = aOther.mIPAddress;
        }

        return *this;
    }

    //ipv4Extend& operator= ( ipv4Extend && aOther )
    //{
    //    if ( this != &aOther )
    //    {
    //    }

    //    return *this;
    //}

    ipv4Extend& operator++()
    {
        increase();

        return *this;
    }

    ipv4Extend operator++( int )
    {
        ipv4Extend sTemp( *this );
        operator++();
        return sTemp;
    }

    friend bool operator== ( const ipv4Extend   & aLhs,
                             const ipv4Extend   & aRhs )
    {
        return (aLhs.mIPAddress == aRhs.mIPAddress);
    }

    friend bool operator!= ( const ipv4Extend   & aLhs,
                             const ipv4Extend   & aRhs )
    {
        return !( aLhs == aRhs );
    }

    friend bool operator < ( const ipv4Extend   & aLhs,
                             const ipv4Extend   & aRhs )
    {
        bool    sIsRigthBigger = false;

        for ( int i = 0; i < 4; i++ )
        {
            if ( aLhs.mIPAddress[i] < aRhs.mIPAddress[i] )
            {
                sIsRigthBigger = true;
                break;
            }
        }

        return sIsRigthBigger;
    }

    friend bool operator > ( const ipv4Extend   & aLhs,
                             const ipv4Extend   & aRhs )
    {
        return aRhs < aLhs;
    }

    friend bool operator <= ( const ipv4Extend   & aLhs,
                              const ipv4Extend   & aRhs )
    {
        return !( aLhs > aRhs );
    }

    friend bool operator >= ( const ipv4Extend   & aLhs,
                              const ipv4Extend   & aRhs )
    {
        return !( aLhs < aRhs );
    }

private:

    void increase()
    {
        bool    sIsRoundUp = false;

        mIPAddress[3]++;

        for ( auto sIt = mIPAddress.rbegin();
              sIt != mIPAddress.rend();
              sIt++ )
        {
            //std::cout << *sIt << "\n";
            if ( sIsRoundUp == true )
            {
                *sIt = *sIt + 1;
                //*sIt++;
            }

            //std::cout << *sIt << "\n";
            if ( (*sIt) > 255 )
            {
                sIsRoundUp = true;
                *sIt = 0;
            }
            else
            {
                sIsRoundUp = false;
            }
        }

        if ( sIsRoundUp == true )
        {
            mIPAddress[3]++;
        }
    }

    void decrease()
    {
        bool sIsRoundDown = false;
    }
};

max_stars_repo_name: akaseon/ModernCppChallengeStudy

max_issues_repo_name: akaseon/ModernCppChallengeStudy

max_forks_repo_name: akaseon/ModernCppChallengeStudy


text:
#ifndef tabulate_H_INCLUDED
#define tabulate_H_INCLUDED

//#define DEBUG

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include "wf.h"
#include "glasma.h"
#include "jet.h"
#include "bfkl.h"
#include <gsl/gsl_errno.h>

#define END_MSG 0
#define PIPE_MSG 1
#define RET_MSG 2

int myid;
//number of processors
int np;
//each point in phi is done on an individual slave node
//the code is setup such that the number of points in phi is the number of processors - 1

MPI_Status Status;

int master_main(int, char**) ;
int slave_main(int, char**) ;

char tag[256];

//*
//*
// This section of code initializes and reads in the
// following global variables from file
//*
//*
#define X_FIELDS \
   X(int, wfTAG, "%d") \
   X(int, A1, "%d") \
   X(int, A2, "%d") \
   X(double, rts, "%lf") \

#include "ReadParameters.cxx"
//

#endif

max_stars_repo_name: kdusling/mpc

max_issues_repo_name: kdusling/mpc

max_forks_repo_name: kdusling/mpc


text:
/**
 *
 * @file qwrapper_dsygst.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.8.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated d Fri Apr  1 11:02:42 2016
 *
 **/
#include <lapacke.h>
#include <al4san.h>
#include <runtime/al4san_quark.h>
#include <runtime/al4san_starpu.h>
//#include <runtime/al4san_parsec.h>
#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 **/

AL4SAN_TASK_CPU(dsygst, dsygst_cpu_func)
void EIG_AL4SAN_CORE_dsygst(AL4SAN_option_t *options,
                int itype, int uplo, int n,
                AL4SAN_desc_t *A, int Am, int An, int lda,
                AL4SAN_desc_t *B, int Bm, int Bn, int ldb,
                int iinfo)
{

    AL4SAN_Insert_Task(AL4SAN_TASK(dsygst), (AL4SAN_option_t*)options,
        AL4SAN_VALUE,                       &itype,                      sizeof(int),
        AL4SAN_VALUE,                       &uplo,                      sizeof(int),
        AL4SAN_VALUE,                       &n,                      sizeof(int),
        AL4SAN_INOUT,                       AL4SAN_ADDR(A, double, Am, An),   AL4SAN_DEP,
        AL4SAN_VALUE,                       &lda,                      sizeof(int),
        AL4SAN_INPUT,                       AL4SAN_ADDR(B, double, Bm, Bn),   AL4SAN_DEP,
        AL4SAN_VALUE,                       &ldb,                      sizeof(int),
        AL4SAN_VALUE,                       &iinfo,                      sizeof(int),
       ARG_END);
}

/***************************************************************************//**
 *
 **/
void dsygst_cpu_func(AL4SAN_arg_list *al4san_arg)
{
        int itype;
        int uplo;
        int n;
        double *A;
        int lda;
        double *B;
        int ldb;
        int iinfo;

        int info;

        AL4SAN_Unpack_Arg(al4san_arg, &itype, &uplo, &n, &A, &lda, &B, &ldb, &iinfo);

        char u;
         if(uplo==Al4sanUpper)
            u='U';
         else
           u='L';

        info = LAPACKE_dsygst_work(
                        LAPACK_COL_MAJOR,
                        itype,
                        u,
                        n, A, lda, B, ldb);
}


text:
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
// reserved. See files LICENSE and NOTICE for details.
//
// This file is part of CEED, a collection of benchmarks, miniapps, software
// libraries and APIs for efficient high-order finite element and spectral
// element discretizations for exascale applications. For more information and
// source code availability see http://github.com/ceed.
//
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
// a collaborative effort of two U.S. Department of Energy organizations (Office
// of Science and the National Nuclear Security Administration) responsible for
// the planning and preparation of a capable exascale ecosystem, including
// software, applications, hardware, advanced system engineering and early
// testbed platforms, in support of the nation's exascale computing imperative.

//                        libCEED + PETSc Example: CEED BPs
//
// This example demonstrates a simple usage of libCEED with PETSc to solve the
// CEED BP benchmark problems, see http://ceed.exascaleproject.org/bps.
//
// The code uses higher level communication protocols in DMPlex.
//
// Build with:
//
//     make bps [PETSC_DIR=</path/to/petsc>] [CEED_DIR=</path/to/libceed>]
//
// Sample runs:
//
//     ./bps -problem bp1 -degree 3
//     ./bps -problem bp2 -degree 3
//     ./bps -problem bp3 -degree 3
//     ./bps -problem bp4 -degree 3
//     ./bps -problem bp5 -degree 3 -ceed /cpu/self
//     ./bps -problem bp6 -degree 3 -ceed /gpu/cuda
//
//TESTARGS -ceed {ceed_resource} -test -problem bp5 -degree 3 -ksp_max_it_clip 15,15

/// @file
/// CEED BPs example using PETSc with DMPlex
/// See bpsraw.c for a "raw" implementation using a structured grid.
const char help[] = "Solve CEED BPs using PETSc with DMPlex\n";

#include <stdbool.h>
#include <string.h>
#include <ceed.h>
#include <petsc.h>
#include <petscdmplex.h>
#include <petscksp.h>
#include <petscsys.h>

#include "bps.h"
#include "include/bpsproblemdata.h"
#include "include/petscmacros.h"
#include "include/petscutils.h"
#include "include/matops.h"
#include "include/structs.h"
#include "include/libceedsetup.h"

#if PETSC_VERSION_LT(3,12,0)
#ifdef PETSC_HAVE_CUDA
#include <petsccuda.h>
// Note: With PETSc prior to version 3.12.0, providing the source path to
//       include 'cublas_v2.h' will be needed to use 'petsccuda.h'.
#endif
#endif

// -----------------------------------------------------------------------------
// Utilities
// -----------------------------------------------------------------------------

// Utility function, compute three factors of an integer
static void Split3(PetscInt size, PetscInt m[3], bool reverse) {
  for (PetscInt d=0, size_left=size; d<3; d++) {
    PetscInt try = (PetscInt)PetscCeilReal(PetscPowReal(size_left, 1./(3 - d)));
    while (try * (size_left / try) != size_left) try++;
    m[reverse ? 2-d : d] = try;
    size_left /= try;
  }
}

static int Max3(const PetscInt a[3]) {
  return PetscMax(a[0], PetscMax(a[1], a[2]));
}

static int Min3(const PetscInt a[3]) {
  return PetscMin(a[0], PetscMin(a[1], a[2]));
}

// -----------------------------------------------------------------------------
// Parameter structure for running problems
// -----------------------------------------------------------------------------
typedef struct RunParams_ *RunParams;
struct RunParams_ {
  MPI_Comm comm;
  PetscBool test_mode, read_mesh, user_l_nodes, write_solution;
  char *filename, *hostname;
  PetscInt local_nodes, degree, q_extra, dim, num_comp_u, *mesh_elem;
  PetscInt ksp_max_it_clip[2];
  PetscMPIInt ranks_per_node;
  BPType bp_choice;
  PetscLogStage solve_stage;
};

// -----------------------------------------------------------------------------
// Main body of program, called in a loop for performance benchmarking purposes
// -----------------------------------------------------------------------------
static PetscErrorCode RunWithDM(RunParams rp, DM dm,
                                const char *ceed_resource) {
  PetscErrorCode ierr;
  double my_rt_start, my_rt, rt_min, rt_max;
  PetscInt xl_size, l_size, g_size;
  PetscScalar *r;
  Vec X, X_loc, rhs, rhs_loc;
  Mat mat_O;
  KSP ksp;
  UserO user_O;
  Ceed ceed;
  CeedData ceed_data;
  CeedQFunction qf_error;
  CeedOperator op_error;
  CeedVector rhs_ceed, target;
  VecType vec_type;
  PetscMemType mem_type;

  PetscFunctionBeginUser;
  // Set up libCEED
  CeedInit(ceed_resource, &ceed);
  CeedMemType mem_type_backend;
  CeedGetPreferredMemType(ceed, &mem_type_backend);

  ierr = DMGetVecType(dm, &vec_type); CHKERRQ(ierr);
  if (!vec_type) { // Not yet set by user -dm_vec_type
    switch (mem_type_backend) {
    case CEED_MEM_HOST: vec_type = VECSTANDARD; break;
    case CEED_MEM_DEVICE: {
      const char *resolved;
      CeedGetResource(ceed, &resolved);
      if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
      else if (strstr(resolved, "/gpu/hip/occa"))
        vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678
      else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
      else vec_type = VECSTANDARD;
    }
    }
    ierr = DMSetVecType(dm, vec_type); CHKERRQ(ierr);
  }

  // Create global and local solution vectors
  ierr = DMCreateGlobalVector(dm, &X); CHKERRQ(ierr);
  ierr = VecGetLocalSize(X, &l_size); CHKERRQ(ierr);
  ierr = VecGetSize(X, &g_size); CHKERRQ(ierr);
  ierr = DMCreateLocalVector(dm, &X_loc); CHKERRQ(ierr);
  ierr = VecGetSize(X_loc, &xl_size); CHKERRQ(ierr);
  ierr = VecDuplicate(X, &rhs); CHKERRQ(ierr);

  // Operator
  ierr = PetscMalloc1(1, &user_O); CHKERRQ(ierr);
  ierr = MatCreateShell(rp->comm, l_size, l_size, g_size, g_size,
                        user_O, &mat_O); CHKERRQ(ierr);
  ierr = MatShellSetOperation(mat_O, MATOP_MULT,
                              (void(*)(void))MatMult_Ceed); CHKERRQ(ierr);
  ierr = MatShellSetOperation(mat_O, MATOP_GET_DIAGONAL,
                              (void(*)(void))MatGetDiag); CHKERRQ(ierr);
  ierr = MatShellSetVecType(mat_O, vec_type); CHKERRQ(ierr);

  // Print summary
  if (!rp->test_mode) {
    PetscInt P = rp->degree + 1, Q = P + rp->q_extra;

    const char *used_resource;
    CeedGetResource(ceed, &used_resource);

    VecType vec_type;
    ierr = VecGetType(X, &vec_type); CHKERRQ(ierr);

    PetscInt c_start, c_end;
    ierr = DMPlexGetHeightStratum(dm, 0, &c_start, &c_end); CHKERRQ(ierr);
    PetscMPIInt comm_size;
    ierr = MPI_Comm_size(rp->comm, &comm_size); CHKERRQ(ierr);
    ierr = PetscPrintf(rp->comm,
                       "\n-- CEED Benchmark Problem %d -- libCEED + PETSc --\n"
                       "  MPI:\n"
                       "    Hostname                           : %s\n"
                       "    Total ranks                        : %d\n"
                       "    Ranks per compute node             : %d\n"
                       "  PETSc:\n"
                       "    PETSc Vec Type                     : %s\n"
                       "  libCEED:\n"
                       "    libCEED Backend                    : %s\n"
                       "    libCEED Backend MemType            : %s\n"
                       "  Mesh:\n"
                       "    Number of 1D Basis Nodes (P)       : %d\n"
                       "    Number of 1D Quadrature Points (Q) : %d\n"
                       "    Global nodes                       : %D\n"
                       "    Local Elements                     : %D\n"
                       "    Owned nodes                        : %D\n"
                       "    DoF per node                       : %D\n",
                       rp->bp_choice+1, rp->hostname, comm_size,
                       rp->ranks_per_node, vec_type, used_resource,
                       CeedMemTypes[mem_type_backend],
                       P, Q, g_size/rp->num_comp_u, c_end - c_start, l_size/rp->num_comp_u,
                       rp->num_comp_u);
    CHKERRQ(ierr);
  }

  // Create RHS vector
  ierr = VecDuplicate(X_loc, &rhs_loc); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs_loc); CHKERRQ(ierr);
  ierr = VecGetArrayAndMemType(rhs_loc, &r, &mem_type); CHKERRQ(ierr);
  CeedVectorCreate(ceed, xl_size, &rhs_ceed);
  CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r);

  ierr = PetscMalloc1(1, &ceed_data); CHKERRQ(ierr);
  ierr = SetupLibceedByDegree(dm, ceed, rp->degree, rp->dim, rp->q_extra,
                              rp->dim, rp->num_comp_u, g_size, xl_size, bp_options[rp->bp_choice],
                              ceed_data, true, rhs_ceed, &target); CHKERRQ(ierr);

  // Gather RHS
  CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL);
  ierr = VecRestoreArrayAndMemType(rhs_loc, &r); CHKERRQ(ierr);
  ierr = VecZeroEntries(rhs); CHKERRQ(ierr);
  ierr = DMLocalToGlobal(dm, rhs_loc, ADD_VALUES, rhs); CHKERRQ(ierr);
  CeedVectorDestroy(&rhs_ceed);

  // Create the error QFunction
  CeedQFunctionCreateInterior(ceed, 1, bp_options[rp->bp_choice].error,
                              bp_options[rp->bp_choice].error_loc, &qf_error);
  CeedQFunctionAddInput(qf_error, "u", rp->num_comp_u, CEED_EVAL_INTERP);
  CeedQFunctionAddInput(qf_error, "true_soln", rp->num_comp_u, CEED_EVAL_NONE);
  CeedQFunctionAddOutput(qf_error, "error", rp->num_comp_u, CEED_EVAL_NONE);

  // Create the error operator
  CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
                     &op_error);
  CeedOperatorSetField(op_error, "u", ceed_data->elem_restr_u,
                       ceed_data->basis_u, CEED_VECTOR_ACTIVE);
  CeedOperatorSetField(op_error, "true_soln", ceed_data->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, target);
  CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_u_i,
                       CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);

  // Set up Mat
  user_O->comm = rp->comm;
  user_O->dm = dm;
  user_O->X_loc = X_loc;
  ierr = VecDuplicate(X_loc, &user_O->Y_loc); CHKERRQ(ierr);
  user_O->x_ceed = ceed_data->x_ceed;
  user_O->y_ceed = ceed_data->y_ceed;
  user_O->op = ceed_data->op_apply;
  user_O->ceed = ceed;

  ierr = KSPCreate(rp->comm, &ksp); CHKERRQ(ierr);
  {
    PC pc;
    ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);
    if (rp->bp_choice == CEED_BP1 || rp->bp_choice == CEED_BP2) {
      ierr = PCSetType(pc, PCJACOBI); CHKERRQ(ierr);
      ierr = PCJacobiSetType(pc, PC_JACOBI_ROWSUM); CHKERRQ(ierr);
    } else {
      ierr = PCSetType(pc, PCNONE); CHKERRQ(ierr);
    }
    ierr = KSPSetType(ksp, KSPCG); CHKERRQ(ierr);
    ierr = KSPSetNormType(ksp, KSP_NORM_NATURAL); CHKERRQ(ierr);
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT,
                            PETSC_DEFAULT); CHKERRQ(ierr);
  }
  ierr = KSPSetOperators(ksp, mat_O, mat_O); CHKERRQ(ierr);

  // First run's performance log is not considered for benchmarking purposes
  ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT, 1);
  CHKERRQ(ierr);
  my_rt_start = MPI_Wtime();
  ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr);
  my_rt = MPI_Wtime() - my_rt_start;
  ierr = MPI_Allreduce(MPI_IN_PLACE, &my_rt, 1, MPI_DOUBLE, MPI_MIN, rp->comm);
  CHKERRQ(ierr);
  // Set maxits based on first iteration timing
  if (my_rt > 0.02) {
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT,
                            rp->ksp_max_it_clip[0]);
    CHKERRQ(ierr);
  } else {
    ierr = KSPSetTolerances(ksp, 1e-10, PETSC_DEFAULT, PETSC_DEFAULT,
                            rp->ksp_max_it_clip[1]);
    CHKERRQ(ierr);
  }
  ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);

  // Timed solve
  ierr = VecZeroEntries(X); CHKERRQ(ierr);
  ierr = PetscBarrier((PetscObject)ksp); CHKERRQ(ierr);

  // -- Performance logging
  ierr = PetscLogStagePush(rp->solve_stage); CHKERRQ(ierr);

  // -- Solve
  my_rt_start = MPI_Wtime();
  ierr = KSPSolve(ksp, rhs, X); CHKERRQ(ierr);
  my_rt = MPI_Wtime() - my_rt_start;

  // -- Performance logging
  ierr = PetscLogStagePop();

  // Output results
  {
    KSPType ksp_type;
    KSPConvergedReason reason;
    PetscReal rnorm;
    PetscInt its;
    ierr = KSPGetType(ksp, &ksp_type); CHKERRQ(ierr);
    ierr = KSPGetConvergedReason(ksp, &reason); CHKERRQ(ierr);
    ierr = KSPGetIterationNumber(ksp, &its); CHKERRQ(ierr);
    ierr = KSPGetResidualNorm(ksp, &rnorm); CHKERRQ(ierr);
    if (!rp->test_mode || reason < 0 || rnorm > 1e-8) {
      ierr = PetscPrintf(rp->comm,
                         "  KSP:\n"
                         "    KSP Type                           : %s\n"
                         "    KSP Convergence                    : %s\n"
                         "    Total KSP Iterations               : %D\n"
                         "    Final rnorm                        : %e\n",
                         ksp_type, KSPConvergedReasons[reason], its,
                         (double)rnorm); CHKERRQ(ierr);
    }
    if (!rp->test_mode) {
      ierr = PetscPrintf(rp->comm,"  Performance:\n"); CHKERRQ(ierr);
    }
    {
      PetscReal max_error;
      ierr = ComputeErrorMax(user_O, op_error, X, target, &max_error);
      CHKERRQ(ierr);
      PetscReal tol = 5e-2;
      if (!rp->test_mode || max_error > tol) {
        ierr = MPI_Allreduce(&my_rt, &rt_min, 1, MPI_DOUBLE, MPI_MIN, rp->comm);
        CHKERRQ(ierr);
        ierr = MPI_Allreduce(&my_rt, &rt_max, 1, MPI_DOUBLE, MPI_MAX, rp->comm);
        CHKERRQ(ierr);
        ierr = PetscPrintf(rp->comm,
                           "    Pointwise Error (max)              : %e\n"
                           "    CG Solve Time                      : %g (%g) sec\n",
                           (double)max_error, rt_max, rt_min); CHKERRQ(ierr);
      }
    }
    if (!rp->test_mode) {
      ierr = PetscPrintf(rp->comm,
                         "    DoFs/Sec in CG                     : %g (%g) million\n",
                         1e-6*g_size*its/rt_max,
                         1e-6*g_size*its/rt_min); CHKERRQ(ierr);
    }
  }

  if (rp->write_solution) {
    PetscViewer vtk_viewer_soln;

    ierr = PetscViewerCreate(rp->comm, &vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerSetType(vtk_viewer_soln, PETSCVIEWERVTK); CHKERRQ(ierr);
    ierr = PetscViewerFileSetName(vtk_viewer_soln, "solution.vtu"); CHKERRQ(ierr);
    ierr = VecView(X, vtk_viewer_soln); CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&vtk_viewer_soln); CHKERRQ(ierr);
  }

  // Cleanup
  ierr = VecDestroy(&X); CHKERRQ(ierr);
  ierr = VecDestroy(&X_loc); CHKERRQ(ierr);
  ierr = VecDestroy(&user_O->Y_loc); CHKERRQ(ierr);
  ierr = MatDestroy(&mat_O); CHKERRQ(ierr);
  ierr = PetscFree(user_O); CHKERRQ(ierr);
  ierr = CeedDataDestroy(0, ceed_data); CHKERRQ(ierr);

  ierr = VecDestroy(&rhs); CHKERRQ(ierr);
  ierr = VecDestroy(&rhs_loc); CHKERRQ(ierr);
  ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
  CeedVectorDestroy(&target);
  CeedQFunctionDestroy(&qf_error);
  CeedOperatorDestroy(&op_error);
  CeedDestroy(&ceed);
  PetscFunctionReturn(0);
}

static PetscErrorCode Run(RunParams rp, PetscInt num_resources,
                          char *const *ceed_resources, PetscInt num_bp_choices,
                          const BPType *bp_choices) {
  PetscInt ierr;
  DM dm;

  PetscFunctionBeginUser;
  // Setup DM
  if (rp->read_mesh) {
    ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD, rp->filename, PETSC_TRUE, &dm);
    CHKERRQ(ierr);
  } else {
    if (rp->user_l_nodes) {
      // Find a nicely composite number of elements no less than global nodes
      PetscMPIInt size;
      ierr = MPI_Comm_size(rp->comm, &size); CHKERRQ(ierr);
      for (PetscInt g_elem =
             PetscMax(1, size * rp->local_nodes / PetscPowInt(rp->degree, rp->dim));
           ;
           g_elem++) {
        Split3(g_elem, rp->mesh_elem, true);
        if (Max3(rp->mesh_elem) / Min3(rp->mesh_elem) <= 2) break;
      }
    }
    ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, rp->dim, PETSC_FALSE,
                               rp->mesh_elem,
                               NULL, NULL, NULL, PETSC_TRUE, &dm); CHKERRQ(ierr);
  }

  {
    DM dm_dist = NULL;
    PetscPartitioner part;

    ierr = DMPlexGetPartitioner(dm, &part); CHKERRQ(ierr);
    ierr = PetscPartitionerSetFromOptions(part); CHKERRQ(ierr);
    ierr = DMPlexDistribute(dm, 0, NULL, &dm_dist); CHKERRQ(ierr);
    if (dm_dist) {
      ierr = DMDestroy(&dm); CHKERRQ(ierr);
      dm  = dm_dist;
    }
  }
  // Disable default VECSTANDARD *after* distribution (which creates a Vec)
  ierr = DMSetVecType(dm, NULL); CHKERRQ(ierr);

  for (PetscInt b = 0; b < num_bp_choices; b++) {
    DM dm_deg;
    VecType vec_type;
    PetscInt q_extra = rp->q_extra;
    rp->bp_choice = bp_choices[b];
    rp->num_comp_u = bp_options[rp->bp_choice].num_comp_u;
    rp->q_extra = q_extra < 0 ? bp_options[rp->bp_choice].q_extra : q_extra;
    ierr = DMClone(dm, &dm_deg); CHKERRQ(ierr);
    ierr = DMGetVecType(dm, &vec_type); CHKERRQ(ierr);
    ierr = DMSetVecType(dm_deg, vec_type); CHKERRQ(ierr);
    // Create DM
    PetscInt dim;
    ierr = DMGetDimension(dm_deg, &dim); CHKERRQ(ierr);
    ierr = SetupDMByDegree(dm_deg, rp->degree, rp->num_comp_u, dim,
                           bp_options[rp->bp_choice].enforce_bc,
                           bp_options[rp->bp_choice].bc_func); CHKERRQ(ierr);
    for (PetscInt r = 0; r < num_resources; r++) {
      ierr = RunWithDM(rp, dm_deg, ceed_resources[r]); CHKERRQ(ierr);
    }
    ierr = DMDestroy(&dm_deg); CHKERRQ(ierr);
    rp->q_extra = q_extra;
  }

  ierr = DMDestroy(&dm); CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

int main(int argc, char **argv) {
  PetscInt ierr, comm_size;
  RunParams rp;
  MPI_Comm comm;
  char filename[PETSC_MAX_PATH_LEN];
  char *ceed_resources[30];
  PetscInt num_ceed_resources = 30;
  char hostname[PETSC_MAX_PATH_LEN];

  PetscInt dim = 3, mesh_elem[3] = {3, 3, 3};
  PetscInt num_degrees = 30, degree[30] = {}, num_local_nodes = 2,
                                          local_nodes[2] = {};
  PetscMPIInt ranks_per_node;
  PetscBool degree_set;
  BPType bp_choices[10];
  PetscInt num_bp_choices = 10;

  // Initialize PETSc
  ierr = PetscInitialize(&argc, &argv, NULL, help);
  if (ierr) return ierr;
  comm = PETSC_COMM_WORLD;
  ierr = MPI_Comm_size(comm, &comm_size);
  if (ierr != MPI_SUCCESS) return ierr;
  #if defined(PETSC_HAVE_MPI_PROCESS_SHARED_MEMORY)
  {
    MPI_Comm splitcomm;
    ierr = MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL,
                               &splitcomm);
    CHKERRQ(ierr);
    ierr = MPI_Comm_size(splitcomm, &ranks_per_node); CHKERRQ(ierr);
    ierr = MPI_Comm_free(&splitcomm); CHKERRQ(ierr);
  }
  #else
  ranks_per_node = -1; // Unknown
  #endif

  // Setup all parameters needed in Run()
  ierr = PetscMalloc1(1, &rp); CHKERRQ(ierr);
  rp->comm = comm;

  // Read command line options
  ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL);
  CHKERRQ(ierr);
  {
    PetscBool set;
    ierr = PetscOptionsEnumArray("-problem", "CEED benchmark problem to solve",
                                 NULL,
                                 bp_types, (PetscEnum *)bp_choices, &num_bp_choices, &set);
    CHKERRQ(ierr);
    if (!set) {
      bp_choices[0] = CEED_BP1;
      num_bp_choices = 1;
    }
  }
  rp->test_mode = PETSC_FALSE;
  ierr = PetscOptionsBool("-test",
                          "Testing mode (do not print unless error is large)",
                          NULL, rp->test_mode, &rp->test_mode, NULL); CHKERRQ(ierr);
  rp->write_solution = PETSC_FALSE;
  ierr = PetscOptionsBool("-write_solution", "Write solution for visualization",
                          NULL, rp->write_solution, &rp->write_solution, NULL);
  CHKERRQ(ierr);
  degree[0] = rp->test_mode ? 3 : 2;
  ierr = PetscOptionsIntArray("-degree",
                              "Polynomial degree of tensor product basis", NULL,
                              degree, &num_degrees, &degree_set); CHKERRQ(ierr);
  if (!degree_set)
    num_degrees = 1;
  rp->q_extra = PETSC_DECIDE;
  ierr = PetscOptionsInt("-q_extra",
                         "Number of extra quadrature points (-1 for auto)", NULL,
                         rp->q_extra, &rp->q_extra, NULL); CHKERRQ(ierr);
  {
    PetscBool set;
    ierr = PetscOptionsStringArray("-ceed",
                                   "CEED resource specifier (comma-separated list)", NULL,
                                   ceed_resources, &num_ceed_resources, &set); CHKERRQ(ierr);
    if (!set) {
      ierr = PetscStrallocpy( "/cpu/self", &ceed_resources[0]); CHKERRQ(ierr);
      num_ceed_resources = 1;
    }
  }
  ierr = PetscGetHostName(hostname, sizeof hostname); CHKERRQ(ierr);
  ierr = PetscOptionsString("-hostname", "Hostname for output", NULL, hostname,
                            hostname, sizeof(hostname), NULL); CHKERRQ(ierr);
  rp->read_mesh = PETSC_FALSE;
  ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL, filename,
                            filename, sizeof(filename), &rp->read_mesh);
  CHKERRQ(ierr);
  rp->filename = filename;
  if (!rp->read_mesh) {
    PetscInt tmp = dim;
    ierr = PetscOptionsIntArray("-cells", "Number of cells per dimension", NULL,
                                mesh_elem, &tmp, NULL); CHKERRQ(ierr);
  }
  local_nodes[0] = 1000;
  ierr = PetscOptionsIntArray("-local_nodes",
                              "Target number of locally owned nodes per "
                              "process (single value or min,max)",
                              NULL, local_nodes, &num_local_nodes, &rp->user_l_nodes);
  CHKERRQ(ierr);
  if (num_local_nodes < 2)
    local_nodes[1] = 2 * local_nodes[0];
  {
    PetscInt two = 2;
    rp->ksp_max_it_clip[0] = 5;
    rp->ksp_max_it_clip[1] = 20;
    ierr = PetscOptionsIntArray("-ksp_max_it_clip",
                                "Min and max number of iterations to use during benchmarking",
                                NULL, rp->ksp_max_it_clip, &two, NULL); CHKERRQ(ierr);
  }
  if (!degree_set) {
    PetscInt max_degree = 8;
    ierr = PetscOptionsInt("-max_degree",
                           "Range of degrees [1, max_degree] to run with",
                           NULL, max_degree, &max_degree, NULL);
    CHKERRQ(ierr);
    for (PetscInt i = 0; i < max_degree; i++)
      degree[i] = i + 1;
    num_degrees = max_degree;
  }
  {
    PetscBool flg;
    PetscInt p = ranks_per_node;
    ierr = PetscOptionsInt("-p", "Number of MPI ranks per node", NULL,
                           p, &p, &flg);
    CHKERRQ(ierr);
    if (flg) ranks_per_node = p;
  }

  ierr = PetscOptionsEnd();
  CHKERRQ(ierr);

  // Register PETSc logging stage
  ierr = PetscLogStageRegister("Solve Stage", &rp->solve_stage);
  CHKERRQ(ierr);

  rp->hostname = hostname;
  rp->dim = dim;
  rp->mesh_elem = mesh_elem;
  rp->ranks_per_node = ranks_per_node;

  for (PetscInt d = 0; d < num_degrees; d++) {
    PetscInt deg = degree[d];
    for (PetscInt n = local_nodes[0]; n < local_nodes[1]; n *= 2) {
      rp->degree = deg;
      rp->local_nodes = n;
      ierr = Run(rp, num_ceed_resources, ceed_resources,
                 num_bp_choices, bp_choices); CHKERRQ(ierr);
    }
  }
  // Clear memory
  ierr = PetscFree(rp); CHKERRQ(ierr);
  for (PetscInt i=0; i<num_ceed_resources; i++) {
    ierr = PetscFree(ceed_resources[i]); CHKERRQ(ierr);
  }
  return PetscFinalize();
}

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
/**
 *
 * @file core_sgeqp3_init.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @generated s Tue Jan  7 11:44:49 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  CORE_sgeqp3_init initializes jpvt to [1, ..., n].
 *  Uses 1-based indexing for Fortran compatability.
 *
 *******************************************************************************
 *
 *  @param[in] n
 *          Size of vector.
 *
 *  @param[in,out] jpvt
 *          Vector of size n.
 *          On exit, jpvt[i] = i+1.
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_sgeqp3_init = PCORE_sgeqp3_init
#define CORE_sgeqp3_init PCORE_sgeqp3_init
#endif
void CORE_sgeqp3_init( int n, int *jpvt )
{
    int j;
    for( j = 0; j < n; ++j ) {
        jpvt[j] = j+1;
    }
}


text:
/* -*- Mode: C++; tab-width: 4; c-basic-offset: 4; indent-tabs-mode: nil -*- */
/*
 *     Copyright 2018 Couchbase, Inc.
 *
 *   Licensed under the Apache License, Version 2.0 (the "License");
 *   you may not use this file except in compliance with the License.
 *   You may obtain a copy of the License at
 *
 *       http://www.apache.org/licenses/LICENSE-2.0
 *
 *   Unless required by applicable law or agreed to in writing, software
 *   distributed under the License is distributed on an "AS IS" BASIS,
 *   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *   See the License for the specific language governing permissions and
 *   limitations under the License.
 */
#pragma once

#include <gsl/gsl-lite.hpp>
#include <nlohmann/json.hpp>
#include <platform/dirutils.h>
#include <vector>

class Event;

/**
 * The Module class represents the configuration for a single module.
 * See ../README.md for more information.
 */
class Module {
public:
    Module() = delete;
    Module(const nlohmann::json& json,
           const std::string& srcRoot,
           const std::string& objRoot);

    void createHeaderFile();

    /**
     * The name of the module
     */
    std::string name;
    /**
     * The lowest identifier for the audit events in this module. All
     * audit descriptor defined for this module MUST be within the range
     * [start, start + max_events_per_module]
     */
    int64_t start;
    /**
     * The name of the file containing the audit descriptors for this
     * module.
     */
    std::string file;
    /**
     * The JSON data describing the audit descriptors for this module
     */
    nlohmann::json json;
    /**
     * Is this module enterprise only?
     */
    bool enterprise = false;

    /**
     * If present this is the name of a C headerfile to generate with
     * #defines for all audit identifiers for the module.
     */
    std::string header;

    /**
     * A list of all of the events defined for this module
     */
    std::vector<std::unique_ptr<Event>> events;

protected:
    /**
     * Add the event to the list of events for the module
     *
     * @param event the event to add
     * @throws std::invalid_argument if the event is outside the legal range
     *                               for the module
     */
    void addEvent(std::unique_ptr<Event> event);

    /// Parse the event descriptor file and add all of the events into
    /// the list of events
    void parseEventDescriptorFile();
};


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#define COMPEARTH_PRIVATE_UPDOWN_ARGSORT3 1
#define COMPEARTH_PRIVATE_UPDOWN_ABS_ARGSORT3 1
#include "compearth.h"
#ifdef COMPEARTH_USE_MKL
#ifdef __clang__
#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wreserved-id-macro"
#pragma clang diagnostic ignored "-Wstrict-prototypes"
#endif
#include <mkl_lapacke.h>
#ifdef __clang__
#pragma clang diagnostic pop
#endif
#else
#include <lapacke.h>
#endif

/*
static int argsort3_upDown(const double *__restrict__ x,
                           const bool lascend,
                           int *__restrict__ iperm);
static int argsort3(const double *__restrict__ x, int *__restrict__ iperm);
*/

#define LWORK 102 //(NB+2)*N hence NB=32; must be atleast 3*N-1 which is 8

/*!
 * @brief Decomposes a set of moment tensors into eigenvalues + basis.
 *
 * @param[in] nmt    Number of moment tensors.
 * @param[in] M      6 x n moment tensors with some unspecified basis
 *                   (e.g., up-south-east) packed
 *                    M = [M11 M22 M33 M12 M13 M23].
 * @param[in] isort  Sorting of eigenvalues. \n
 *                   If 1 then the eigenvalues are sorted highest to lowest. \n
 *                   If 2 then the eigenvalues are sorted lowest to highest. \n
 *                   If 3 then the eigenvalues are sorted by absolute value
 *                   from highest to lowest. \n
 *                   If 4 then the eigenvalues are sorted by absolute value
 *                   from lowest to highest.
 *
 * @param[out] lam   3 x n set of eigenvalues.
 * @param[out] U     3 x 3 x n set of bases.  The first 3 x 3 matrix is in
 *                   column major order.  This has the same basis as M.
 *
 * @author Carl Tape.  Converted to C by Ben Baker.
 *
 * @copyright MIT
 *
 * @bug The
 */
int compearth_CMTdecom(const int nmt, const double *__restrict__ M,
                       const int isort,
                       double *__restrict__ lam,
                       double *__restrict__ U)
{
    double Lams[3], Mx[9], Ut[9], work[LWORK];
    int perm[3], c, i, ierr, info, r;
    // Error checks
    if (nmt < 1 || M == NULL || lam == NULL || U == NULL)
    {
        if (nmt < 1){fprintf(stderr, "%s: Error no mts\n", __func__);}
        if (M == NULL){fprintf(stderr, "%s: Error M is NULL\n", __func__);}
        if (lam == NULL){fprintf(stderr, "%s: Error lam is NULL\n", __func__);}
        if (U == NULL){fprintf(stderr, "%s: Error U is NULL\n", __func__);}
        return -1;
    }
    if (isort < 1 || isort > 4)
    {
        fprintf(stderr, "%s: Error isort=%d is invalid\n", __func__, isort);
        return -1;
    }
    // Decompose the moment tensors
    for (i=0; i<nmt; i++)
    {
        // Create the symmetric moment tensor matrix from the 6x1 vector
        compearth_Mvec2Mmat(1, M, 1, Mx);
        // Compute eigenvalues in ascending order
        info = LAPACKE_dsyev_work(LAPACK_COL_MAJOR, 'V', 'U', 3, Mx, 3,
                                  Lams, work, LWORK);
        if (info != 0)
        {
            fprintf(stderr, "%s: Error computing eigenvalues\n", __func__);
            return -1;
        }
        // Descending
        if (isort == 1)
        {
            argsort3_upDown(Lams, false, perm);
        }
        else if (isort == 2)
        {
            argsort3_upDown(Lams, true, perm);
        }
        // Descending on absolute value
        else if (isort == 3)
        {
            argsort3_absUpDown(Lams, false, perm);
        }
        else //if (isort == 4)
        {
            argsort3_absUpDown(Lams, true,  perm);
        }
        // And permute and save
        for (r=0; r<3; r++)
        {
            lam[3*i+r] = Lams[perm[r]];
        }
        // Permute columns
        for (c=0; c<3; c++)
        {
            for (r=0; r<3; r++)
            {
                Ut[3*c+r] = Mx[3*perm[c]+r];
            }
        }
        // Ensure the determinant is > 0
        ierr = compearth_Udetcheck(1, Ut, &U[9*i]);
        if (ierr != 0)
        {
            fprintf(stderr, "%s: Error checking determinant\n", __func__);
            return -1;
        }
    }
    return 0;
}

/*
static int argsort3_upDown(const double *__restrict__ x,
                           const bool lascend,
                           int *__restrict__ iperm)
{
    int ipermt[3], ierr;
    if (lascend)
    {
        ierr = argsort3(x, iperm);
    }
    else
    {
        ierr = argsort3(x, ipermt);
        iperm[0] = ipermt[2];
        iperm[1] = ipermt[1];
        iperm[2] = ipermt[0];
    }
    return ierr;
}

static int argsort3(const double *__restrict__ x, int *__restrict__ iperm)
{
    int i, temp;
    const int a = 0;
    const int b = 1;
    const int c = 2;
    // Copy
    iperm[a] = a;
    iperm[b] = b;
    iperm[c] = c;
    if (x[iperm[a]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[a];
        iperm[a] = temp;
    }
    if (x[iperm[a]] > x[iperm[b]])
    {
        temp = iperm[b];
        iperm[b] = iperm[a];
        iperm[a] = temp;
    }
    //Now the smallest element is the first one. Just check the 2-nd and 3-rd
    if (x[iperm[b]] > x[iperm[c]])
    {
        temp = iperm[c];
        iperm[c] = iperm[b];
        iperm[b] = temp;
    }
    // Verify
    for (i=1; i<3; i++)
    {
        if (x[iperm[i-1]] > x[iperm[i]])
        {
            fprintf(stderr, "%s: Failed to sort numbers in ascending order\n",
                     __func__);
            return -1;
        }
    }
    return 0;
}
*/

max_stars_repo_name: OUCyf/mtbeach

max_issues_repo_name: carltape/mtbeach

max_forks_repo_name: carltape/mtbeach