text:
// Copyright (c) 2005 - 2010 Marc de Kamps
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
//
//    * 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 the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
//      without specific prior written permission.
//
// 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 THE COPYRIGHT OWNER OR CONTRIBUTORS 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.
//
//      If you use this software in work leading to a scientific publication, you should cite
//      the 'currently valid reference', which can be found at http://miind.sourceforge.net
#ifndef _CODE_LIBS_NUMTOOLSLIB_NODE_GSLOBJECTS_INCLUDE_GUARD
#define _CODE_LIBS_NUMTOOLSLIB_NODE_GSLOBJECTS_INCLUDE_GUARD

#include <boost/shared_ptr.hpp>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>
#include "Precision.h"

using boost::shared_ptr;

namespace NumtoolsLib {
  typedef unsigned int Number;

	typedef int (*Function)   (double t, const double y[], double dydt[], void * params);
	typedef int (*Derivative) (double t, const double y[], double * dfdy, double dfdt[], void * params);

	//! auxilliary struct to clean up some of the code of AbstractDVIntegrator
	struct GSLObjects{

		const gsl_odeiv_step_type*		_p_step_algorithm; //!< pointer to function, so no ownership issues
		gsl_odeiv_step*					_p_step;
		gsl_odeiv_control*				_p_controller;
		gsl_odeiv_evolve*				_p_evolver;
		gsl_odeiv_system				_system;

		GSLObjects
		(
			const gsl_odeiv_step_type*,	//<! choice for solver type
			Number,						//<! state size
			const Precision&			//<! absolute and relative precision conform GSL documentation
		);

		~GSLObjects();

	};
}

#endif // include guard

max_stars_repo_name: dekamps/miind

max_issues_repo_name: dekamps/miind

max_forks_repo_name: dekamps/miind


text:
/* linalg/cholesky_band.c
 *
 * Copyright (C) 2018, 2019, 2020 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 <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>

static double cholesky_band_norm1(const gsl_matrix * A);
static int cholesky_band_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);

/*
gsl_linalg_cholesky_band_decomp()
  Cholesky decomposition of a square symmetric positive definite banded
matrix

Inputs: A - matrix in banded format, N-by-ndiag where N is the size of
            the matrix and ndiag is the number of nonzero diagonals.

Notes:
1) The main diagonal of the Cholesky factor is stored in the first column
of A; the first subdiagonal in the second column and so on.

2) If ndiag > 1, the 1-norm of A is stored in A(N,ndiag) on output

3) At each diagonal element, the matrix is factored as

A(j:end,j:end) = [ A11 A21^T ] = [ alpha  0 ] [ alpha v^T ]
                 [ A21 A22   ]   [   v    L ] [   0   L^T ]

where:

alpha = sqrt(A(j,j))
v = A(j+1:end, j) / alpha
A22 = L L^T + v v^T

So we start at A(1,1) and work right. Pseudo-code is:

loop j = 1, ..., N
  alpha = sqrt(A(j,j))
  A(j+1:end, j) := A(j+1:end, j) / alpha     (DSCAL)
  A(j+1:end, j+1:end) -= v v^T               (DSYR)

Due to the banded structure, v has at most p non-zero elements, where
p is the lower bandwidth
*/

int
gsl_linalg_cholesky_band_decomp(gsl_matrix * A)
{
  const size_t N = A->size1;     /* size of matrix */
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */

  if (ndiag > N)
    {
      GSL_ERROR ("invalid matrix dimensions", GSL_EBADLEN);
    }
  else
    {
      const size_t p = ndiag - 1; /* lower bandwidth */
      const int kld = (int) GSL_MAX(1, p);
      size_t j;

      if (ndiag > 1)
        {
          /*
           * calculate 1-norm of A and store in lower right of matrix, which is not accessed
           * by rest of routine. gsl_linalg_cholesky_band_rcond() will use this later. If
           * A is diagonal, there is no empty slot to store the 1-norm, so the rcond routine
           * will have to reconstruct it from the Cholesky factor.
           */
          double Anorm = cholesky_band_norm1(A);
          gsl_matrix_set(A, N - 1, p, Anorm);
        }

      for (j = 0; j < N; ++j)
        {
          double ajj = gsl_matrix_get(A, j, 0);
          size_t lenv;

          if (ajj <= 0.0)
            {
              GSL_ERROR("matrix is not positive definite", GSL_EDOM);
            }

          ajj = sqrt(ajj);
          gsl_matrix_set(A, j, 0, ajj);

          /* number of elements in v, which will normally be p, unless we
           * are in lower right corner of matrix */
          lenv = GSL_MIN(p, N - j - 1);

          if (lenv > 0)
            {
              gsl_vector_view v = gsl_matrix_subrow(A, j, 1, lenv);
              gsl_matrix_view m = gsl_matrix_submatrix(A, j + 1, 0, lenv, lenv);

              gsl_blas_dscal(1.0 / ajj, &v.vector);

              m.matrix.tda = kld;
              gsl_blas_dsyr(CblasUpper, -1.0, &v.vector, &m.matrix);
            }
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_band_solve (const gsl_matrix * LLT,
                                const gsl_vector * b,
                                gsl_vector * x)
{
  if (LLT->size1 != b->size)
    {
      GSL_ERROR ("matrix size must match b size", GSL_EBADLEN);
    }
  else if (LLT->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      int status;

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

      status = gsl_linalg_cholesky_band_svx(LLT, x);

      return status;
    }
}

int
gsl_linalg_cholesky_band_svx (const gsl_matrix * LLT, gsl_vector * x)
{
  if (LLT->size1 != x->size)
    {
      GSL_ERROR ("matrix size must match solution size", GSL_EBADLEN);
    }
  else
    {
      /* solve for c using forward-substitution, L c = b */
      cblas_dtbsv(CblasColMajor, CblasLower, CblasNoTrans, CblasNonUnit,
                  (int) LLT->size1, (int) (LLT->size2 - 1), LLT->data, LLT->tda,
                  x->data, x->stride);

      /* perform back-substitution, L^T x = c */
      cblas_dtbsv(CblasColMajor, CblasLower, CblasTrans, CblasNonUnit,
                  (int) LLT->size1, (int) (LLT->size2 - 1), LLT->data, LLT->tda,
                  x->data, x->stride);

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_band_solvem (const gsl_matrix * LLT,
                                 const gsl_matrix * B,
                                 gsl_matrix * X)
{
  if (LLT->size1 != B->size1)
    {
      GSL_ERROR ("LLT size1 must match B size1", GSL_EBADLEN);
    }
  else if (LLT->size1 != X->size1)
    {
      GSL_ERROR ("LLT size1 must match solution size1", GSL_EBADLEN);
    }
  else if (B->size2 != X->size2)
    {
      GSL_ERROR ("B size2 must match X size2", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* copy X <- B */
      gsl_matrix_memcpy (X, B);

      status = gsl_linalg_cholesky_band_svxm(LLT, X);

      return status;
    }
}

int
gsl_linalg_cholesky_band_svxm (const gsl_matrix * LLT, gsl_matrix * X)
{
  if (LLT->size1 != X->size1)
    {
      GSL_ERROR ("LLT size1 must match solution size1", GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t nrhs = X->size2;
      size_t j;

      for (j = 0; j < nrhs; ++j)
        {
          gsl_vector_view xj = gsl_matrix_column(X, j);

          status = gsl_linalg_cholesky_band_svx (LLT, &xj.vector);
          if (status)
            return status;
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_band_invert (const gsl_matrix * LLT, gsl_matrix * Ainv)
{
  if (Ainv->size1 != Ainv->size2)
    {
      GSL_ERROR("Ainv must be square", GSL_ENOTSQR);
    }
  else if (LLT->size1 != Ainv->size1)
    {
      GSL_ERROR("cholesky matrix has different dimensions from Ainv", GSL_EBADLEN);
    }
  else
    {
      int status;

      /* unpack Cholesky factor into lower triangle of Ainv */
      status = gsl_linalg_cholesky_band_unpack(LLT, Ainv);
      if (status)
        return status;

      /* call the standard Cholesky inversion routine */
      status = gsl_linalg_cholesky_invert(Ainv);
      if (status)
        return status;

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_band_unpack (const gsl_matrix * LLT, gsl_matrix * L)
{
  const size_t N = LLT->size1;

  if (N != L->size1)
    {
      GSL_ERROR("L matrix does not match LLT dimensions", GSL_EBADLEN);
    }
  else if (L->size1 != L->size2)
    {
      GSL_ERROR("L matrix is not square", GSL_ENOTSQR);
    }
  else
    {
      const size_t p = LLT->size2 - 1; /* lower bandwidth */
      size_t i;

      for (i = 0; i < p + 1; ++i)
        {
          gsl_vector_const_view v = gsl_matrix_const_subcolumn(LLT, i, 0, N - i);
          gsl_vector_view w = gsl_matrix_subdiagonal(L, i);

          gsl_vector_memcpy(&w.vector, &v.vector);
        }

      /* zero out remaining subdiagonals */
      for (i = p + 1; i < N; ++i)
        {
          gsl_vector_view w = gsl_matrix_subdiagonal(L, i);
          gsl_vector_set_zero(&w.vector);
        }

      return GSL_SUCCESS;
    }
}

int
gsl_linalg_cholesky_band_rcond (const gsl_matrix * LLT, double * rcond, gsl_vector * work)
{
  const size_t N = LLT->size1;

  if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      int status;
      const size_t ndiag = LLT->size2;
      double Anorm;    /* ||A||_1 */
      double Ainvnorm; /* ||A^{-1}||_1 */

      if (ndiag == 1)
        {
          /* diagonal matrix, compute 1-norm since it has not been stored */
          gsl_vector_const_view v = gsl_matrix_const_column(LLT, 0);
          Anorm = gsl_vector_max(&v.vector);
          Anorm = Anorm * Anorm;
        }
      else
        {
          /* 1-norm is stored in A(N, ndiag) by gsl_linalg_cholesky_band_decomp() */
          Anorm = gsl_matrix_get(LLT, N - 1, ndiag - 1);
        }

      *rcond = 0.0;

      /* return if matrix is singular */
      if (Anorm == 0.0)
        return GSL_SUCCESS;

      status = gsl_linalg_invnorm1(N, cholesky_band_Ainv, (void *) LLT, &Ainvnorm, work);
      if (status)
        return status;

      if (Ainvnorm != 0.0)
        *rcond = (1.0 / Anorm) / Ainvnorm;

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_cholesky_band_scale()
  This function computes scale factors diag(S), such that

diag(S) A diag(S)

has a condition number within a factor N of the matrix
with the smallest condition number over all possible
diagonal scalings. See Corollary 7.6 of:

N. J. Higham, Accuracy and Stability of Numerical Algorithms (2nd Edition),
SIAM, 2002.

Inputs: A - symmetric positive definite matrix
        S - (output) scale factors, S_i = 1 / sqrt(A_ii)
*/

int
gsl_linalg_cholesky_band_scale(const gsl_matrix * A, gsl_vector * S)
{
  const size_t N = A->size1;     /* size of matrix */
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */

  if (ndiag > N)
    {
      GSL_ERROR ("invalid matrix dimensions", GSL_EBADLEN);
    }
  else if (N != S->size)
    {
      GSL_ERROR("S must have length N", GSL_EBADLEN);
    }
  else
    {
      size_t i;

      /* compute S_i = 1/sqrt(A_{ii}) */
      for (i = 0; i < N; ++i)
        {
          double Aii = gsl_matrix_get(A, i, 0);

          if (Aii <= 0.0)
            gsl_vector_set(S, i, 1.0); /* matrix not positive definite */
          else
            gsl_vector_set(S, i, 1.0 / sqrt(Aii));
        }

      return GSL_SUCCESS;
    }
}

/*
gsl_linalg_cholesky_band_scale_apply()
  This function applies scale transformation to A:

A <- diag(S) A diag(S)

Inputs: A - (input/output)
            on input, symmetric positive definite matrix in banded format
            on output, diag(S) * A * diag(S) in banded format
        S - (input) scale factors
*/

int
gsl_linalg_cholesky_band_scale_apply(gsl_matrix * A, const gsl_vector * S)
{
  const size_t N = A->size1;     /* size of matrix */
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */

  if (ndiag > N)
    {
      GSL_ERROR ("invalid matrix dimensions", GSL_EBADLEN);
    }
  else if (N != S->size)
    {
      GSL_ERROR("S must have length N", GSL_EBADLEN);
    }
  else
    {
      size_t i, j;

      for (j = 0; j < N; ++j)
        {
          double sj = gsl_vector_get(S, j);

          for (i = j; i < GSL_MIN(N, j + ndiag); ++i)
            {
              double si = gsl_vector_get(S, i);
              double * ptr = gsl_matrix_ptr(A, j, i - j);
              *ptr *= sj * si;
            }
        }

      return GSL_SUCCESS;
    }
}

/* compute 1-norm of symmetric banded matrix */
static double
cholesky_band_norm1(const gsl_matrix * A)
{
  const size_t N = A->size1;
  const size_t ndiag = A->size2; /* number of diagonals in band, including main diagonal */
  double value;

  if (ndiag == 1)
    {
      /* diagonal matrix */
      gsl_vector_const_view v = gsl_matrix_const_column(A, 0);
      CBLAS_INDEX_t idx = gsl_blas_idamax(&v.vector);
      value = gsl_vector_get(&v.vector, idx);
    }
  else
    {
      size_t j;

      value = 0.0;
      for (j = 0; j < N; ++j)
        {
          size_t ncol = GSL_MIN(ndiag, N - j); /* number of elements in column j below and including main diagonal */
          gsl_vector_const_view v = gsl_matrix_const_subrow(A, j, 0, ncol);
          double sum = gsl_blas_dasum(&v.vector);
          size_t k, l;

          /* sum now contains the absolute sum of elements below and including main diagonal for column j; we
           * have to add the symmetric elements above the diagonal */
          k = j;
          l = 1;
          while (k > 0 && l < ndiag)
            {
              double Akl = gsl_matrix_get(A, --k, l++);
              sum += fabs(Akl);
            }

          value = GSL_MAX(value, sum);
        }
    }

  return value;
}

/* x := A^{-1} x = A^{-t} x, A = L L^T */
static int
cholesky_band_Ainv(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  gsl_matrix * LLT = (gsl_matrix * ) params;

  (void) TransA; /* unused parameter warning */

  /* compute x := L^{-1} x */
  cblas_dtbsv(CblasColMajor, CblasLower, CblasNoTrans, CblasNonUnit,
              (int) LLT->size1, (int) (LLT->size2 - 1), LLT->data, LLT->tda,
              x->data, x->stride);

  /* compute x := L^{-T} x */
  cblas_dtbsv(CblasColMajor, CblasLower, CblasTrans, CblasNonUnit,
              (int) LLT->size1, (int) (LLT->size2 - 1), LLT->data, LLT->tda,
              x->data, x->stride);

  return GSL_SUCCESS;
}

ext: c

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once

#include <gsl/span>

#include <cstddef>
#include <type_traits>

namespace imageview {
namespace detail {

template <class T, typename Enable = void>
class HasColorTypeTypedef : public std::false_type {};

template <class T>
class HasColorTypeTypedef<T, std::void_t<typename T::color_type>> : public std::true_type {};

template <class T, typename Enable = void>
class HasKBytesPerPixelConstant : public std::false_type {};

template <class T>
class HasKBytesPerPixelConstant<T, std::enable_if_t<std::is_integral_v<decltype(T::kBytesPerPixel)> &&
                                                    std::is_const_v<decltype(T::kBytesPerPixel)>>>
    : public std::true_type {};

template <class T, class Enable = void>
class HasRead : public std::false_type {};

template <class T>
class HasRead<T, std::enable_if_t<std::is_same_v<typename T::color_type,
                                                 decltype(std::declval<const T&>().read(
                                                     std::declval<gsl::span<const std::byte, T::kBytesPerPixel>>()))>>>
    : public std::true_type {};

template <class T, class Enable = void>
class HasWrite : public std::false_type {};

template <class T>
class HasWrite<T, std::enable_if_t<std::is_same_v<void, decltype(std::declval<const T&>().write(
                                                            std::declval<const typename T::color_type&>(),
                                                            std::declval<gsl::span<std::byte, T::kBytesPerPixel>>()))>>>
    : public std::true_type {};

}  // namespace detail

template <class T>
class IsPixelFormat : public std::conjunction<
                          // Has a member typedef 'color_type'.
                          detail::HasColorTypeTypedef<T>,
                          // Has an integral static member constant 'kBytesPerPixel'.
                          detail::HasKBytesPerPixelConstant<T>,
                          // Has a member function read() with the signature equivalent to
                          //   color_type read(gsl::span<const std::byte, kBytesPerPixel>) const;
                          detail::HasRead<T>,
                          // Has a member function write() with the signature equivalent to
                          //   void write(const color_type&, gsl::span<std::byte, kBytesPerPixel>) const;
                          detail::HasWrite<T>> {};

}  // namespace imageview

ext: h

max_stars_repo_name: alexanderbelous/imageview

max_issues_repo_name: alexanderbelous/imageview

max_forks_repo_name: alexanderbelous/imageview


text:
#ifndef KGS_EIGENVALUE_H
#define KGS_EIGENVALUE_H

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

class Eigenvalue {
    protected:
        const int m, n; ///< Dimensions of matrix

    public:
        Eigenvalue(gsl_matrix* matrix1, gsl_matrix* matrix2, gsl_matrix* matrix3);
        gsl_matrix * const Jacobianmatrix;  //TODO: Make private
        gsl_matrix * const Hessiancartesian;       //TODO: Make private
        gsl_matrix * const Hessiantorsionangle;       //TODO: Make private
        gsl_matrix * const Massmatrix;
        gsl_vector_complex * const eigenvaluealpha;
        gsl_vector * const eigenvaluebeta;

        gsl_matrix_complex * const eigenvector;
        gsl_matrix* times(gsl_matrix* matrix1, gsl_matrix* matrix2) const;

        virtual ~Eigenvalue();

        gsl_matrix* getHessiantorsionangle() const;

        void setSingularvalue() const;

        gsl_vector* getSingularvalue() const;

        gsl_vector * const singularvalue;
    };

#endif //KGS_EIGENVALUE_H


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
#include <mpi.h>
#include <StGermain/libStGermain/src/StGermain.h>
#include <petsc.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

const Type StGermain_Type = "StGermain";

#include "StGermain_Tools.h"

StgData* StgInit( int argc, char* argv[] ) {
   StgData* data = (StgData*) malloc(sizeof(StgData));
   *data = (StgData){.comm = NULL, .rank=-1, .nProcs=-1, .dictionary=NULL, .argcCpy=0, .argvCpy=NULL };

   //lets copy all this data for safety
   data->argcCpy = argc;
   data->argvCpy = (char**) malloc((argc+1)*sizeof(char*));
   int ii;
   for(ii = 0; ii<argc; ii++){
      data->argvCpy[ii] = (char*)malloc(strlen(argv[ii])+1);
      strcpy(data->argvCpy[ii],argv[ii]);
   }
   data->argvCpy[argc] = NULL;  //add sentinel

   int is_inited;
   MPI_Initialized( &is_inited );
   if(!is_inited)
      MPI_Init( &argc, &argv );
   MPI_Comm_dup( MPI_COMM_WORLD, &data->comm );
   MPI_Comm_size( data->comm, &data->nProcs );
   MPI_Comm_rank( data->comm, &data->rank );
   StGermain_Init( &(data->argcCpy), &(data->argvCpy) );

   return data;
}

int StgFinalise(StgData* data){
   /* Close off everything */
//   StGermain_Finalise();   // avoid this finalise, as py interpreter has its own ideas about when it will destroy objects.
   PetscFinalize();
//   MPI_Finalize();

   /* free up these guys created earlier */
   int ii;
   for(ii = 0; ii<data->argcCpy; ii++)
   	  free(data->argvCpy[ii]);
   free(data->argvCpy);
   free(data);
   return 0; /* success */
}

void StgAbort(StgData* data){
     MPI_Abort( data->comm, EXIT_FAILURE );
}


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code for the implementation of the Fourier domain response for LISA-like detectors.
 *
 */

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"
#include "LISAgeometry.h"
#include "LISAFDresponse.h"

/***************************************/
/********* Core functions **************/

/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle - here simplified version for just the y_21 observable */
/* Older version of th FD response: two stages, first for the orbital delay (Bessel phase) and second for the constellation delay/modulation */
int LISASimFDResponse21(
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,  /* Input/Output: list of modes in Frequency-domain amplitude and phase form as produced by the ROM, and output after FD response processing */
  const double inclination,                                   /* Inclination of the source */
  const double lambda,                                        /* First angle for the position in the sky */
  const double beta,                                          /* Second angle for the position in the sky */
  const double psi)                                           /* Polarization angle */
{
  /* Computing the complicated trigonometric coefficients */
  //clock_t begsetcoeffs = clock();
  SetCoeffsG(lambda, beta, psi);
  //clock_t endsetcoeffs = clock();
  //printf("Set Coeffs time: %g s\n", (double)(endsetcoeffs - begsetcoeffs) / CLOCKS_PER_SEC);

  /* Main loop over the modes - goes through all the modes present */
  ListmodesCAmpPhaseFrequencySeries* listelement = *list;
  while(listelement) {

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;
    int len = (int) freq->size;
    double f, tf, bphi;
    double complex camp;

    /* Computing the Ylm combined factors for plus and cross for this mode */
    /* Capital Phi is set to 0 by convention */
    double complex Yfactorplus;
    double complex Yfactorcross;
    if (!(l%2)) {
      Yfactorplus = 1/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }
    else {
      Yfactorplus = 1/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }

    /* First step of the processing: orbital delay */
    /* Initializing spline for the phase */
    gsl_spline* spline_phi = gsl_spline_alloc(gsl_interp_cspline, len);
    gsl_interp_accel* accel_phi = gsl_interp_accel_alloc();
    gsl_spline_init(spline_phi, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(phase, 0), len);
    /* Vector keeping track of the Bessel Phase correction - to be used in the next step */
    gsl_vector* besselphi = gsl_vector_alloc(len);

    /* Loop over the frequencies - computing the correction due to the orbital delay */
    for(int j=0; j<len; j++) {
      f = gsl_vector_get(freq, j);
      tf =  gsl_spline_eval_deriv(spline_phi, f, accel_phi)/(2*PI);
      bphi = -2*PI*f*variant->OrbitR/C_SI*cos(beta) * cos( variant->OrbitOmega*tf + variant->OrbitPhi0 - lambda );
      camp = gsl_vector_get(amp_real, j) * cexp(I*bphi); /* Amplitude is real before applying this first delay */
      gsl_vector_set(amp_real, j, creal(camp));
      gsl_vector_set(amp_imag, j, cimag(camp));
      gsl_vector_set(besselphi, j, bphi);
    }

    /* Second step of the processing: constellation delay/modulation */
    /* Initializing spline for the bessel phase */
    gsl_spline* spline_besselphi = gsl_spline_alloc(gsl_interp_cspline, len);
    gsl_interp_accel* accel_besselphi = gsl_interp_accel_alloc();
    gsl_spline_init(spline_besselphi, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(besselphi, 0), len);
    /* Loop over the frequencies - computing the correction due to the orbital delay */
    for(int j=0; j<len; j++) {
      f = gsl_vector_get(freq, j);
      tf =  (gsl_spline_eval_deriv(spline_phi, f, accel_phi) + gsl_spline_eval_deriv(spline_besselphi, f, accel_besselphi))/(2*PI);
      camp = G21mode(variant, f, tf, Yfactorplus, Yfactorcross) * (gsl_vector_get(amp_real, j) + I * gsl_vector_get(amp_imag, j));
      /**/
      gsl_vector_set(amp_real, j, creal(camp));
      gsl_vector_set(amp_imag, j, cimag(camp));
    }

    listelement = listelement->next;

    /* Clean up */
    gsl_spline_free(spline_phi);
    gsl_interp_accel_free(accel_phi);
    gsl_vector_free(besselphi);
    gsl_spline_free(spline_besselphi);
    gsl_interp_accel_free(accel_besselphi);
  }

  return SUCCESS;
}

//WARNING: tRef is ignored for now in the response - i.e. set to 0
/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle */
int LISASimFDResponsey12(
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **listy12,   /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the y12 observable */
  const double torb,                                       /* Reference orbital time - tf as read from the hlm gives t-tinj, this arg allows to pass tinj to the response */
  const double lambda,                                        /* First angle for the position in the sky */
  const double beta,                                          /* Second angle for the position in the sky */
  const double inclination,                                   /* Inclination of the source */
  const double psi,                                           /* Polarization angle */
  const int tagfrozenLISA,                                    /* Tag to treat LISA as frozen at its torb configuration  */
  const ResponseApproxtag responseapprox)                     /* Tag to select possible low-f approximation level in FD response */
{
  /* Computing the complicated trigonometric coefficients */
  //clock_t begsetcoeffs = clock();
  SetCoeffsG(lambda, beta, psi);
  //clock_t endsetcoeffs = clock();
  //printf("Set Coeffs time: %g s\n", (double)(endsetcoeffs - begsetcoeffs) / CLOCKS_PER_SEC);

  /* Main loop over the modes - goes through all the modes present, stopping when encountering NULL */
  ListmodesCAmpPhaseFrequencySeries* listelement = *list;
  while(listelement) {

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;
    int len = (int) freq->size;
    double f, tf, tforb;
    double complex g12mode = 0.;
    double complex g21mode = 0.;
    double complex g23mode = 0.;
    double complex g32mode = 0.;
    double complex g31mode = 0.;
    double complex g13mode = 0.;
    double complex camp_y;

    /* Computing the Ylm combined factors for plus and cross for this mode */
    /* Capital Phi is set to 0 by convention */
    double complex Yfactorplus;
    double complex Yfactorcross;
    if (!(l%2)) {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }
    else {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }

    /* Initializing spline for the phase */
    gsl_spline* spline_phi = gsl_spline_alloc(gsl_interp_cspline, len);
    gsl_interp_accel* accel_phi = gsl_interp_accel_alloc();
    gsl_spline_init(spline_phi, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(phase, 0), len);

    /* Orbital delay, constellation delay/modulation */
    /* Initializing frequency series structure for this mode, for each of the TDI observables */
    CAmpPhaseFrequencySeries *modefreqseries = NULL;
    CAmpPhaseFrequencySeries_Init(&modefreqseries, len);
    gsl_vector* freq_y = modefreqseries->freq;
    gsl_vector* amp_real_y = modefreqseries->amp_real;
    gsl_vector* amp_imag_y = modefreqseries->amp_imag;
    gsl_vector* phase_y = modefreqseries->phase;
    /* Loop over the frequencies - computing the correction due to the constellation delay/modulation */
    //clock_t tbegcontesllation = clock();
    //double timingcumulativeGABmode = 0;
    for(int j=0; j<len; j++) {
      f = gsl_vector_get(freq, j);
      tf = (gsl_spline_eval_deriv(spline_phi, f, accel_phi))/(2*PI);
      /* tf read from hlm is t-tinj - here convert to orbital time */
      if(!(tagfrozenLISA)) {
        tforb = tf + torb;
      } else {
        tforb = torb;
      }
      //clock_t tbegGAB = clock();
	  EvaluateGABmode(variant, &g12mode, &g21mode, &g23mode, &g32mode, &g31mode, &g13mode, f, tforb, Yfactorplus, Yfactorcross, 1, responseapprox); /* does include the R-delay term */
      //clock_t tendGAB = clock();
      //timingcumulativeGABmode += (double) (tendGAB-tbegGAB) /CLOCKS_PER_SEC;
      /**/
      camp_y = (gsl_vector_get(amp_real, j) + I * gsl_vector_get(amp_imag, j)) * g12mode;
      /**/
      gsl_vector_set(amp_real_y, j, creal(camp_y));
      gsl_vector_set(amp_imag_y, j, cimag(camp_y));
    }
    //clock_t tendcontesllation = clock();
  //printf("Set constellation time: %g s\n", (double)(tendcontesllation - tbegcontesllation) / CLOCKS_PER_SEC);
  //printf("GAB cumulated time: %g s\n", timingcumulativeGABmode);

    /* Copying the vectors of frequencies and phases */
    gsl_vector_memcpy(freq_y, freq);
    gsl_vector_memcpy(phase_y, phase);

    /* Append the modes to the ouput list-of-modes structures */
    *listy12 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*listy12, modefreqseries, l, m);

    /* Going to the next mode in the list */
    listelement = listelement->next;

    /* Clean up */
    gsl_spline_free(spline_phi);
    gsl_interp_accel_free(accel_phi);
  }

  return SUCCESS;
}

//WARNING: tRef is ignored for now in the response - i.e. set to 0
/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LISA response, for given values of the inclination, position in the sky and polarization angle */
int LISASimFDResponseTDI3Chan(
  int tagtRefatLISA,                                          /* 0 to measure Tref from SSB arrival, 1 at LISA guiding center */
  LISAconstellation *variant,                                 /* Provides specifics on the variant of LISA */
  struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI1,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 1 */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI2,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 2 */
  struct tagListmodesCAmpPhaseFrequencySeries **listTDI3,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the TDI channel 3 */
  const double torb,                                       /* Reference orbital time - tf as read from the hlm gives t-tinj, this arg allows to pass tinj to the response */
  const double lambda,                                     /* First angle for the position in the sky */
  const double beta,                                       /* Second angle for the position in the sky */
  const double inclination,                                /* Inclination of the source */
  const double psi,                                        /* Polarization angle */
  const double m1,                                         /* m1 in solar masses - used for resampling */
  const double m2,                                         /* m2 in solar masses - used for resampling */
  const double maxf,                                       /* Maximal frequency to consider - used to ignore hard-to-resolve response at f>1Hz - NOTE: for now, no recomputation of the boundary, so when not resampling can lose a bit of support between the last frequency point covered and maxf */
  const TDItag tditag,                                     /* Selector for the set of TDI observables */
  const int tagfrozenLISA,                                 /* Tag to treat LISA as frozen at its torb configuration  */
  const ResponseApproxtag responseapprox)                  /* Tag to select possible low-f approximation level in FD response */
{
  /* Computing the complicated trigonometric coefficients */
  //clock_t begsetcoeffs = clock();
  SetCoeffsG(lambda, beta, psi);
  //clock_t endsetcoeffs = clock();
  //printf("Set Coeffs time: %g s\n", (double)(endsetcoeffs - begsetcoeffs) / CLOCKS_PER_SEC);

  /* Chirp mass for resampling */
  double mchirp = Mchirpofm1m2(m1, m2);

  /* Main loop over the modes - goes through all the modes present, stopping when encountering NULL */
  ListmodesCAmpPhaseFrequencySeries* listelement = *list;
  while(listelement) {

    /* Definitions: l,m, frequency series and length */
    int l = listelement->l;
    int m = listelement->m;
    CAmpPhaseFrequencySeries* freqseries = listelement->freqseries;
    gsl_vector* freq = freqseries->freq;
    gsl_vector* amp_real = freqseries->amp_real;
    gsl_vector* amp_imag = freqseries->amp_imag;
    gsl_vector* phase = freqseries->phase;
    int len = (int) freq->size;

    //
    //printf("len: %d\n", len);
    //
    //if(l==2&&m==1) printf("no resampling 21\n");
    //      if(l==2&&m==1) {for(int i=0; i<len; i++) printf("%d, %g, %g, %g, %g\n", i, gsl_vector_get(freq, i), gsl_vector_get(amp_real, i), gsl_vector_get(amp_imag, i), gsl_vector_get(phase, i));};

    double f, tf, tforb;
    double complex g21mode = 0.;
    double complex g12mode = 0.;
    double complex g32mode = 0.;
    double complex g23mode = 0.;
    double complex g13mode = 0.;
    double complex g31mode = 0.;
    double complex camp;
    double complex camp1;
    double complex camp2;
    double complex camp3;
    double complex factor1 = 0.;
    double complex factor2 = 0.;
    double complex factor3 = 0.;

    /* Computing the Ylm combined factors for plus and cross for this mode */
    /* Capital Phi is set to 0 by convention */
    double complex Yfactorplus;
    double complex Yfactorcross;
    if (!(l%2)) {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }
    else {
      Yfactorplus = 1./2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) - conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
      Yfactorcross = I/2 * (SpinWeightedSphericalHarmonic(inclination, 0., -2, l, m) + conj(SpinWeightedSphericalHarmonic(inclination, 0., -2, l, -m)));
    }

    /* Initializing spline for the phase - will be used to compute tf */
    gsl_spline* spline_phi = gsl_spline_alloc(gsl_interp_cspline, len);
    gsl_interp_accel* accel_phi = gsl_interp_accel_alloc();
    gsl_spline_init(spline_phi, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(phase, 0), len);

    /* Resampling at high f to achieve a deltaf of at most 0.002 Hz */
    /* Resample linearly at this deltaf when this threshold is reached */
    /* NOTE: Assumes input frequencies are logarithmic (except maybe first interval) to evaluate when to resample */
    gsl_vector* freqrhigh = NULL;
    SetMaxdeltafResampledFrequencies(&freqrhigh, freq, maxf, 0.002); /* Use 0.002Hz as a default maximal deltaf */

    /* Resampling at low f to achieve a deltat of at most 2 weeks */
    /* Resample linearly in time until this threshold is reached */
    /* NOTE: uses Newtonian estimates - also requires the chirp mass as extra information */
    /* It would also be possible to use t(f) (more accurate), but we would also need to interpolate f(t) */
    /* Accuracy in time of this resampling is not critical, estimate should be enough */
    gsl_vector* freqr = NULL;
    SetMaxdeltatResampledFrequencies(&freqr, freqrhigh, 1./24, mchirp, m); /* Use half a month as a default maximal deltaft */

    /* Evaluate resampled waveform */
    CAmpPhaseFrequencySeries* freqseriesr = NULL;
    CAmpPhaseFrequencySeries_Resample(&freqseriesr, freqseries, freqr);
    gsl_vector* freq_resample = freqseriesr->freq;
    gsl_vector* amp_real_resample = freqseriesr->amp_real;
    gsl_vector* amp_imag_resample = freqseriesr->amp_imag;
    gsl_vector* phase_resample = freqseriesr->phase;
    int len_resample = (int) freq_resample->size;

    // /* Determine frequencies to use */
    // /* Because we will need the interpolation on the Re/Im amplitude to resolve the structure of the L-response at high fequencies, we add points beyond 0.01 Hz, with a linear sampling */
    // /* WARNING : It seemed 600 points between 0.1 and 3 Hz (deltaf=0.005) should give interpolation errors below 1e-4 - considering a simple sin(2 pi f L) */
    // /* But first test with TDIA show the sampling should be reduced 10-fold - possibly large increase in cost */
    // /* We use here deltaf=0.0005 until the cause is better understood */
    // /* NOTE: the structure in the response due to the R-delay gives much larger interpolation errors - here we assume the R-delay term is now treated as a phase */
    // int resampled = 0; /* Keeps track of wether or not we resampled and allocated new resources we need to free */
    // double maxfsignal = gsl_vector_get(freq, len-1);
    // double fHigh = fmin(maxf, maxfsignal);
    // /* BEWARE : as Mathematica tests show, this is way too pessimistic - normally deltaf=0.002Hz sould work */
    // /* To be investigated */
    // double fHigh_log_samp = 0.002;
    // double deltaflineartarget = 0.00002; //1e-5 better, but slower
    // int ifmax = len-1; /* last index to be covered in original sampling overall */
    // while((gsl_vector_get(freq, ifmax)>fHigh) && ifmax>0) ifmax--;
    // int imaxlogsampling = ifmax; /* last index to be covered with original sampling */
    // while((gsl_vector_get(freq, imaxlogsampling)>fHigh_log_samp) && imaxlogsampling>0) imaxlogsampling--;
    // gsl_vector* freq_resample = NULL; /*  */
    // gsl_vector* amp_real_resample = NULL;
    // gsl_vector* amp_imag_resample = NULL;
    // gsl_vector* phase_resample = NULL;
    // int len_resample;
    //
    // if((fHigh>fHigh_log_samp) && ((fHigh - gsl_vector_get(freq, imaxlogsampling))/deltaflineartarget)>ifmax-imaxlogsampling) { /* condition to check if the linear sampling will add points - if not, do nothing */
    //
    //   resampled = 1;
    //   /* Number of pts in resampled part */
    //   int nbfreqlinear = ceil((fHigh - gsl_vector_get(freq, imaxlogsampling))/deltaflineartarget);
    //   double deltaflinear = (fHigh - gsl_vector_get(freq, imaxlogsampling))/(nbfreqlinear + 1);
    //   /* Initialize new vectors */
    //   len_resample = imaxlogsampling + 1 + nbfreqlinear;
    //   freq_resample = gsl_vector_alloc(len_resample);
    //   amp_real_resample = gsl_vector_alloc(len_resample);
    //   amp_imag_resample = gsl_vector_alloc(len_resample);
    //   phase_resample = gsl_vector_alloc(len_resample);
    //   /* Build interpolation for original amp_real, amp_imag and phase */
    //   /* NOTE: we could use spline_phi here, written this way for clarity */
    //   gsl_spline* spline_amp_real = gsl_spline_alloc(gsl_interp_cspline, len);
    //   gsl_spline* spline_amp_imag = gsl_spline_alloc(gsl_interp_cspline, len);
    //   gsl_spline* spline_phase = gsl_spline_alloc(gsl_interp_cspline, len);
    //   gsl_interp_accel* accel_amp_real = gsl_interp_accel_alloc();
    //   gsl_interp_accel* accel_amp_imag = gsl_interp_accel_alloc();
    //   gsl_interp_accel* accel_phase = gsl_interp_accel_alloc();
    //   gsl_spline_init(spline_amp_real, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(amp_real, 0), len);
    //   gsl_spline_init(spline_amp_imag, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(amp_imag, 0), len);
    //   gsl_spline_init(spline_phase, gsl_vector_const_ptr(freq, 0), gsl_vector_const_ptr(phase, 0), len);
    //   /* Set resampled frequencies and values */
    //   for(int j=0; j<=imaxlogsampling; j++) {
    //     gsl_vector_set(freq_resample, j, gsl_vector_get(freq, j));
    //     gsl_vector_set(amp_real_resample, j, gsl_vector_get(amp_real, j));
    //     gsl_vector_set(amp_imag_resample, j, gsl_vector_get(amp_imag, j));
    //     gsl_vector_set(phase_resample, j, gsl_vector_get(phase, j));
    //   }
    //   double fimax = gsl_vector_get(freq, imaxlogsampling);
    //   for(int j=imaxlogsampling+1; j<len_resample; j++) {
    //     f = fimax + (j-imaxlogsampling) * deltaflinear;
    //     gsl_vector_set(freq_resample, j, f);
    //     gsl_vector_set(amp_real_resample, j, gsl_spline_eval(spline_amp_real, f, accel_amp_real));
    //     gsl_vector_set(amp_imag_resample, j, gsl_spline_eval(spline_amp_imag, f, accel_amp_imag));
    //     gsl_vector_set(phase_resample, j, gsl_spline_eval(spline_phase, f, accel_phase));
    //   }
    //   /* Free interpolation functions */
    //   gsl_spline_free(spline_amp_real);
    //   gsl_spline_free(spline_amp_imag);
    //   gsl_spline_free(spline_phase);
    //   gsl_interp_accel_free(accel_amp_real);
    //   gsl_interp_accel_free(accel_amp_imag);
    //   gsl_interp_accel_free(accel_phase);
    // }
    // else { /* If no resampling, use the values we had as input */
    //   resampled = 0;
    //   len_resample = imaxlogsampling + 1; /* If maxf < maxfsignal, we will cut the signal above maxf by simply adjusting the range of indices included (NOTE: without recomputing the exact boundary, so some support is lost due to discretization) */
    //   freq_resample = freq;
    //   amp_real_resample = amp_real;
    //   amp_imag_resample = amp_imag;
    //   phase_resample = phase;
    // }

    /* Initializing frequency series structure for this mode, for each of the TDI observables */
    CAmpPhaseFrequencySeries *modefreqseries1 = NULL;
    CAmpPhaseFrequencySeries *modefreqseries2 = NULL;
    CAmpPhaseFrequencySeries *modefreqseries3 = NULL;
    CAmpPhaseFrequencySeries_Init(&modefreqseries1, len_resample);
    CAmpPhaseFrequencySeries_Init(&modefreqseries2, len_resample);
    CAmpPhaseFrequencySeries_Init(&modefreqseries3, len_resample);
    gsl_vector* freq1 = modefreqseries1->freq;
    gsl_vector* amp_real1 = modefreqseries1->amp_real;
    gsl_vector* amp_imag1 = modefreqseries1->amp_imag;
    gsl_vector* phase1 = modefreqseries1->phase;
    gsl_vector* freq2 = modefreqseries2->freq;
    gsl_vector* amp_real2 = modefreqseries2->amp_real;
    gsl_vector* amp_imag2 = modefreqseries2->amp_imag;
    gsl_vector* phase2 = modefreqseries2->phase;
    gsl_vector* freq3 = modefreqseries3->freq;
    gsl_vector* amp_real3 = modefreqseries3->amp_real;
    gsl_vector* amp_imag3 = modefreqseries3->amp_imag;
    gsl_vector* phase3 = modefreqseries3->phase;

    /* Loop over the frequencies */
    //clock_t tbegcontesllation = clock();
    //double timingcumulativeGABmode = 0;
    for(int j=0; j<len_resample; j++) {
      f = gsl_vector_get(freq_resample, j);
      tf = (gsl_spline_eval_deriv(spline_phi, f, accel_phi))/(2*PI);
      /* tf read from hlm is t-tinj - here convert to orbital time using torb - ignore tf if orbit is frozen */
      if(!(tagfrozenLISA)) {
        tforb = tf + torb;
      } else {
        tforb = torb;
      }
      //clock_t tbegGAB = clock();
      EvaluateGABmode(variant, &g12mode, &g21mode, &g23mode, &g32mode, &g31mode, &g13mode, f, tforb, Yfactorplus, Yfactorcross, 0, responseapprox); /* does not include the R-delay term */
      //clock_t tendGAB = clock();
      //timingcumulativeGABmode += (double) (tendGAB-tbegGAB) /CLOCKS_PER_SEC;
      /**/
      EvaluateTDIfactor3Chan(variant, &factor1, &factor2, &factor3, g12mode, g21mode, g23mode, g32mode, g31mode, g13mode, f, tditag, responseapprox);
      double complex amphtilde = gsl_vector_get(amp_real_resample, j) + I * gsl_vector_get(amp_imag_resample, j);
      camp1 = factor1 * amphtilde;
      camp2 = factor2 * amphtilde;
      camp3 = factor3 * amphtilde;
      /* Phase term due to the R-delay, including correction to first order */
      double phase=variant->OrbitOmega*tforb + variant->OrbitPhi0 - lambda;
      double OrbitRoC=variant->OrbitR/C_SI;
      //double phaseRdelay = -2*PI*R_SI/C_SI*f*cos(beta)*cos(Omega_SI*tf - lambda) * (1 + R_SI/C_SI*cos(beta)*Omega_SI*sin(Omega_SI*tf - lambda));
      double phaseRdelay;
      if(tagtRefatLISA==0){//delay so that tinj=torb refers to arrival time at SSB
        phaseRdelay = -2*PI*OrbitRoC*f*cos(beta)*cos(phase) * (1 + cos(beta)*OrbitRoC*variant->OrbitOmega*sin(phase));
      } else {
        //In this version we delay so that tinj=torb is relative to LISAcenter arrival time, so there is no orbital delay when tf = tinj
        //If the original delay realized td=t+d(t), now we want td=t+d(t)-d(t0), that we we change d(t) -> d(t) + d(t0)
        //Then with the approximation d(td)=d(t)(1-ddot(t)),
        double phase0 = variant->OrbitOmega*torb + variant->OrbitPhi0 - lambda;
        phaseRdelay = -2*PI*OrbitRoC*f*cos(beta)*( cos(phase)-cos(phase0) ) * (1 + cos(beta)*OrbitRoC*variant->OrbitOmega*sin(phase));
      }
      if(responseapprox==lowf) { /* In the full low-f approximation, ignore this delay term */
        phaseRdelay = 0.;
      }
      double phasewithRdelay = gsl_vector_get(phase_resample, j) + phaseRdelay;

      /**/
      gsl_vector_set(amp_real1, j, creal(camp1));
      gsl_vector_set(amp_imag1, j, cimag(camp1));
      gsl_vector_set(amp_real2, j, creal(camp2));
      gsl_vector_set(amp_imag2, j, cimag(camp2));
      gsl_vector_set(amp_real3, j, creal(camp3));
      gsl_vector_set(amp_imag3, j, cimag(camp3));
      gsl_vector_set(phase1, j, phasewithRdelay);
      gsl_vector_set(phase2, j, phasewithRdelay);
      gsl_vector_set(phase3, j, phasewithRdelay);
    }
    //clock_t tendconstellation = clock();
    //printf("Set constellation time: %g s\n", (double)(tendconstellation - tbegconstellation) / CLOCKS_PER_SEC);
    //printf("GAB cumulated time: %g s\n", timingcumulativeGABmode);

    /* Copying the vectors of frequencies - we have to allow for the case where it has been shortened */
    gsl_vector_view freq_resample_subview = gsl_vector_subvector(freq_resample, 0, len_resample);
    gsl_vector_memcpy(freq1, &freq_resample_subview.vector);
    gsl_vector_memcpy(freq2, &freq_resample_subview.vector);
    gsl_vector_memcpy(freq3, &freq_resample_subview.vector);

    /* Append the modes to the ouput list-of-modes structures */
    *listTDI1 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*listTDI1, modefreqseries1, l, m);
    *listTDI2 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*listTDI2, modefreqseries2, l, m);
    *listTDI3 = ListmodesCAmpPhaseFrequencySeries_AddModeNoCopy(*listTDI3, modefreqseries3, l, m);

    /* Going to the next mode in the list */
    listelement = listelement->next;

    /* Clean up */
    gsl_spline_free(spline_phi);
    gsl_interp_accel_free(accel_phi);
    gsl_vector_free(freqrhigh);
    gsl_vector_free(freqr);
    CAmpPhaseFrequencySeries_Cleanup(freqseriesr);
    // /* If we used resampling, then we need to free the additional resources that were allocated */
    // if(resampled) {
    //   gsl_vector_free(freq_resample);
    //   gsl_vector_free(amp_real_resample);
    //   gsl_vector_free(amp_imag_resample);
    //   gsl_vector_free(phase_resample);
    // }
  }

  return SUCCESS;
}

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
// @file im2col_conv.c
//
//  \date Created on: Sep 30, 2017
//  \author Gopalakrishna Hegde
//
//   Description:
//
//
//

#include <stdbool.h>
#include <cblas.h>
#include "common_types.h"

//From Berkeley Vision's Caffe
//Refer to Caffe's license : https://github.com/BVLC/caffe/blob/master/LICENSE
static inline bool is_a_ge_zero_and_a_lt_b(int a, int b) {
  return (unsigned int)a < (unsigned int)(b);
}

void Im2Col(const float *data_im, const int channels,
    const int height, const int width, const int kernel_h, const int kernel_w,
    const int pad_h, const int pad_w,
    const int stride_h, const int stride_w,
    float *data_col) {
  const int output_h = (height + 2 * pad_h -
    ((kernel_h - 1) + 1)) / stride_h + 1;
  const int output_w = (width + 2 * pad_w -
    ((kernel_w - 1) + 1)) / stride_w + 1;
  const int channel_size = height * width;
  for (int channel = channels; channel--; data_im += channel_size) {
    for (int kernel_row = 0; kernel_row < kernel_h; kernel_row++) {
      for (int kernel_col = 0; kernel_col < kernel_w; kernel_col++) {
        int input_row = -pad_h + kernel_row;
        for (int output_rows = output_h; output_rows; output_rows--) {
          if (!is_a_ge_zero_and_a_lt_b(input_row, height)) {
            for (int output_cols = output_w; output_cols; output_cols--) {
              *(data_col++) = 0;
            }
          } else {
            int input_col = -pad_w + kernel_col;
            for (int output_col = output_w; output_col; output_col--) {
              if (is_a_ge_zero_and_a_lt_b(input_col, width)) {
                *(data_col++) = data_im[input_row * width + input_col];
              } else {
                *(data_col++) = 0;
              }
              input_col += stride_w;
            }
          }
          input_row += stride_h;
        }
      }
    }
  }
}

void Im2ColConvLayer(const float *input, const float *weight,
    const float *bias, float *scratchpad, const TensorDim in_dim,
    const TensorDim out_dim, const int ker_size, const int group,
    const int pad, const int stride, const int bias_en, float *output) {

  int C = in_dim.c;
  int H = in_dim.h;
  int W = in_dim.w;
  int in_ch_per_group = C / group;
  int out_ch_per_group = out_dim.c / group;
  float alpha = 1;
  float beta = 0;

  for (int b = 0; b < in_dim.n; b++) {
    int in_offset = b * C * H * W;
    Im2Col(input + in_offset, C, H, W, ker_size, ker_size, pad, pad, stride,
           stride, scratchpad);
    for (int g = 0; g < group; g++) {
      int weight_offset = g * ker_size * ker_size * in_ch_per_group * out_ch_per_group;
      int data_offset = g * (ker_size * ker_size * in_ch_per_group) *
          (out_dim.h * out_dim.w);
      int out_offset = (b * out_dim.c + g * out_ch_per_group) *
          out_dim.h * out_dim.w;

      int m = out_ch_per_group;
      int k = in_ch_per_group * ker_size * ker_size;
      int n = out_dim.h * out_dim.w;
      cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                 m, n, k, alpha, weight + weight_offset,
                 k, scratchpad + data_offset, n, beta, output + out_offset, n);
    }
  }

}

max_stars_repo_name: gplhegde/convolution-flavors

max_issues_repo_name: chayitw/convolution-flavors

max_forks_repo_name: chayitw/convolution-flavors


text:
//
// Created by pedram pakseresht on 2/18/21.
//

// #include <petsc.h>

#include <printf.h>
#include <stdlib.h>

void setNumber(double *zxy);

typedef struct{

    double x;
    double z;
    int g[3];

}Pstruct;

void printStruct(Pstruct str) {

    printf("print struct %f %f %d %d %d \n", str.x, str.z, str.g[0], str.g[1], str.g[2]);
//    printf("print x %lf \n", str.z );
//    printf("print x %lf \n", str.g[1] );

};

void setValues(){

}

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

    Pstruct a;
    a.x = 7.0;
    a.z = 8.0;
    a.g[0] = 1;
    a.g[1] = 2;
    a.g[2] = 3;

    Pstruct b;
    b.x = 5.0;
    b.z = 4.0;
    b.g[0] = 3;
    b.g[1] = 1;
    b.g[2] = 2;

    printStruct(a);
    printStruct(b);

//    printf("x of b is %lf \n",b.x);
//printf("size of struc %lu \n", sizeof(Pstruct));
//printf("size of int %lu", sizeof(int [4]));

}

max_stars_repo_name: pakserep/ablateClient

max_issues_repo_name: pakserep/ablateClient

max_forks_repo_name: pakserep/ablateClient


text:
/*

Excited States software: KGS
Contributors: See CONTRIBUTORS.txt
Contact: kgs-contact@simtk.org

Copyright (C) 2009-2017 Stanford University

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

This entire text, including the above copyright notice and this permission notice
shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
IN THE SOFTWARE.

*/
//
// Created by Dominik Budday on 15.03.16.
//

#ifndef KGS_FASTCLASHAVOIDINGMOVE_H
#define KGS_FASTCLASHAVOIDINGMOVE_H

#include <vector>

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

#include "metrics/Metric.h"
#include "metrics/RMSD.h"
#include "moves/Move.h"
#include "core/Configuration.h"

class FastClashAvoidingMove : public Move
{
 public:
  FastClashAvoidingMove( double maxRotation,
                         int trialSteps,
                         const std::string& atomTypes,
                         bool projectConstraints
  );

 protected:
  Configuration* performMove(Configuration* current, gsl_vector* gradient);

 private:
  Configuration* projectOnClashNullspace(Configuration *conf,
                                      gsl_vector *gradient,
                                      std::set<std::pair<Atom *, Atom *> > &collisions);

  gsl_matrix* computeClashAvoidingJacobian( Configuration* conf,
                                            std::map<int,int>& dofMap,
                                            std::set< std::pair<Atom*,Atom*> >& collisions);

  /** Return a map that associates cycle-dofs and constrained dofs with a general dofs. */
  std::map<int,int> collectConstrainedDofMap(Configuration* conf, std::set< std::pair<Atom*,Atom*> >& allCollisions);
  const int m_trialSteps;
  const bool m_projectConstraints;
  const std::string m_collisionCheckAtomTypes;
};

#endif //KGS_FASTCLASHAVOIDINGMOVE_H


text:
#pragma once

#include "halley_gl.h"
#include <gsl/gsl>

namespace Halley
{
	class GLBuffer
	{
	public:
		GLBuffer();
		~GLBuffer();

		void bind();
		void bindToTarget(GLuint index);
		void init(GLenum target, GLenum usage = GL_STREAM_DRAW);
		void setData(gsl::span<const gsl::byte> data);
		size_t getSize() const;

	private:
		GLenum target = 0;
		GLenum usage = 0;
        GLuint name = 0;
		size_t capacity = 0;
		size_t size = 0;
	};
}

max_stars_repo_name: sunhay/halley

max_issues_repo_name: sunhay/halley

max_forks_repo_name: sunhay/halley


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

#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_legendre.h>

#include "fftlog.h"

#include "ccl.h"
#include "ccl_params.h"

/*--------ROUTINE: taper_cl ------
TASK:n Apply cosine tapering to Cls to reduce aliasing
INPUT: number of ell bins for Cl, ell vector, C_ell vector, limits for tapering
       e.g., ell_limits=[low_ell_limit_lower,low_ell_limit_upper,high_ell_limit_lower,high_ell_limit_upper]
*/
static int taper_cl(int n_ell,double *ell,double *cl, double *ell_limits)
{

  for(int i=0;i<n_ell;i++) {
    if(ell[i]<ell_limits[0] || ell[i]>ell_limits[3]) {
      cl[i]=0;//ell outside desirable range
      continue;
    }
    if(ell[i]>=ell_limits[1] && ell[i]<=ell_limits[2])
      continue;//ell within good ell range

    if(ell[i]<ell_limits[1])//tapering low ell
      cl[i]*=cos((ell[i]-ell_limits[1])/(ell_limits[1]-ell_limits[0])*M_PI/2.);

    if(ell[i]>ell_limits[2])//tapering high ell
      cl[i]*=cos((ell[i]-ell_limits[2])/(ell_limits[3]-ell_limits[2])*M_PI/2.);
  }

  return 0;
}

/*--------ROUTINE: ccl_tracer_corr_fftlog ------
TASK: For a given tracer, get the correlation function
      Following function takes a function to calculate angular cl as well.
      By default above function will call it using ccl_angular_cl
INPUT: type of tracer, number of theta values to evaluate = NL, theta vector
 */
static void ccl_tracer_corr_fftlog(ccl_cosmology *cosmo,
				   int n_ell,double *ell,double *cls,
				   int n_theta,double *theta,double *wtheta,
				   int corr_type,int do_taper_cl,double *taper_cl_limits,
				   int *status)
{
  int i;
  double *l_arr,*cl_arr,*th_arr,*wth_arr;

  l_arr=ccl_log_spacing(ccl_splines->ELL_MIN_CORR,ccl_splines->ELL_MAX_CORR,ccl_splines->N_ELL_CORR);
  if(l_arr==NULL) {
    *status=CCL_ERROR_LINSPACE;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  cl_arr=malloc(ccl_splines->N_ELL_CORR*sizeof(double));
  if(cl_arr==NULL) {
    free(l_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  //Interpolate input Cl into array needed for FFTLog
  SplPar *cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
  if(cl_spl==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  double cl_tilt,l_edge,cl_edge;
  l_edge=ell[n_ell-1];
  if((cls[n_ell-1]*cls[n_ell-2]<0) || (cls[n_ell-2]==0)) {
    cl_tilt=0;
    cl_edge=0;
  }
  else {
    cl_tilt=log(cls[n_ell-1]/cls[n_ell-2])/log(ell[n_ell-1]/ell[n_ell-2]);
    cl_edge=cls[n_ell-1];
  }
  for(i=0;i<ccl_splines->N_ELL_CORR;i++) {
    if(l_arr[i]>=l_edge)
      cl_arr[i]=cl_edge*pow(l_arr[i]/l_edge,cl_tilt);
    else
      cl_arr[i]=ccl_spline_eval(l_arr[i],cl_spl);
  }
  ccl_spline_free(cl_spl);

  if (do_taper_cl)
    taper_cl(ccl_splines->N_ELL_CORR,l_arr,cl_arr,taper_cl_limits);

  th_arr=malloc(sizeof(double)*ccl_splines->N_ELL_CORR);
  if(th_arr==NULL) {
    free(l_arr);
    free(cl_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }
  wth_arr=(double *)malloc(sizeof(double)*ccl_splines->N_ELL_CORR);
  if(wth_arr==NULL) {
    free(l_arr); free(cl_arr); free(th_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_fftlog ran out of memory\n");
    return;
  }

  for(i=0;i<ccl_splines->N_ELL_CORR;i++)
    th_arr[i]=0;
  //Although set here to 0, theta is modified by FFTlog to obtain the correlation at ~1/l

  int i_bessel=0;
  if(corr_type==CCL_CORR_GG) i_bessel=0;
  if(corr_type==CCL_CORR_GL) i_bessel=2;
  if(corr_type==CCL_CORR_LP) i_bessel=0;
  if(corr_type==CCL_CORR_LM) i_bessel=4;
  fftlog_ComputeXi2D(i_bessel,ccl_splines->N_ELL_CORR,l_arr,cl_arr,th_arr,wth_arr);

  // Interpolate to output values of theta
  SplPar *wth_spl=ccl_spline_init(ccl_splines->N_ELL_CORR,th_arr,wth_arr,wth_arr[0],0);
  for(i=0;i<n_theta;i++)
    wtheta[i]=ccl_spline_eval(theta[i]*M_PI/180.,wth_spl);
  ccl_spline_free(wth_spl);

  free(l_arr); free(cl_arr);
  free(th_arr); free(wth_arr);

  return;
}

typedef struct {
  int nell;
  double ell0;
  double ellf;
  double cl0;
  double clf;
  int extrapol_0;
  int extrapol_f;
  double tilt0;
  double tiltf;
  SplPar *cl_spl;
  int i_bessel;
  double th;
} corr_int_par;

static double corr_bessel_integrand(double l,void *params)
{
  double cl,jbes;
  corr_int_par *p=(corr_int_par *)params;
  double x=l*p->th;

  if(l<p->ell0) {
    if(p->extrapol_0)
      cl=p->cl0*pow(l/p->ell0,p->tilt0);
    else
      cl=0;
  }
  else if(l>p->ellf) {
    if(p->extrapol_f)
      cl=p->clf*pow(l/p->ellf,p->tiltf);
    else
      cl=0;
  }
  else
    cl=ccl_spline_eval(l,p->cl_spl);

  jbes=gsl_sf_bessel_Jn(p->i_bessel,x);

  return l*jbes*cl;
}

static void ccl_tracer_corr_bessel(ccl_cosmology *cosmo,
				   int n_ell,double *ell,double *cls,
				   int n_theta,double *theta,double *wtheta,
				   int corr_type,int *status)
{
  corr_int_par *cp=malloc(sizeof(corr_int_par));
  if(cp==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_bessel ran out of memory\n");
    return;
  }

  cp->nell=n_ell;
  cp->ell0=ell[0];
  cp->ellf=ell[n_ell-1];
  cp->cl0=cls[0];
  cp->clf=cls[n_ell-1];
  cp->cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
  if(cp->cl_spl==NULL) {
    free(cp);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_bessel ran out of memory\n");
    return;
  }
  switch(corr_type) {
  case CCL_CORR_GG :
    cp->i_bessel=0;
    break;
  case CCL_CORR_GL :
    cp->i_bessel=2;
    break;
  case CCL_CORR_LP :
    cp->i_bessel=0;
    break;
  case CCL_CORR_LM :
    cp->i_bessel=4;
    break;
  }

  if(cls[0]*cls[1]<=0)
    cp->extrapol_0=0;
  else {
    cp->extrapol_0=1;
    cp->tilt0=log10(cls[1]/cls[0])/log10(ell[1]/ell[0]);
  }

  if(cls[n_ell-2]*cls[n_ell-1]<=0)
    cp->extrapol_f=0;
  else {
    cp->extrapol_f=1;
    cp->tiltf=log10(cls[n_ell-1]/cls[n_ell-2])/log10(ell[n_ell-1]/ell[n_ell-2]);
  }

  int ith, gslstatus;
  double result,eresult;
  gsl_function F;
  gsl_integration_workspace *w=gsl_integration_workspace_alloc(ccl_gsl->N_ITERATION);
  for(ith=0;ith<n_theta;ith++) {
    cp->th=theta[ith]*M_PI/180;
    F.function=&corr_bessel_integrand;
    F.params=cp;
    //TODO: Split into intervals between first bessel zeros before integrating
    //This will help both speed and accuracy of the integral.
    gslstatus = gsl_integration_qag(&F, 0, ccl_splines->ELL_MAX_CORR, 0,
                                    ccl_gsl->INTEGRATION_EPSREL, ccl_gsl->N_ITERATION,
                                    ccl_gsl->INTEGRATION_GAUSS_KRONROD_POINTS,
                                    w, &result, &eresult);
    if(gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_correlation.c: ccl_tracer_corr_bessel():");
      *status |= gslstatus;
    }
    wtheta[ith]=result/(2*M_PI);
  }
  gsl_integration_workspace_free(w);
  ccl_spline_free(cp->cl_spl);
  free(cp);
}

/*--------ROUTINE: ccl_compute_legendre_polynomial ------
TASK: Compute input factor for ccl_tracer_corr_legendre
INPUT: tracer 1, tracer 2, i_bessel, theta array, n_theta, L_max, output Pl_theta
 */
static void ccl_compute_legendre_polynomial(int corr_type,double theta,int ell_max,double *Pl_theta)
{
  int i,j;
  double k=0;
  double cth=cos(theta*M_PI/180);

  //Initialize Pl_theta
  for (j=0;j<=ell_max;j++)
      Pl_theta[j]=0.;

  if(corr_type==CCL_CORR_GG) {
    gsl_sf_legendre_Pl_array(ell_max,cth,Pl_theta);
    for (j=0;j<=ell_max;j++)
      Pl_theta[j]*=(2*j+1);
  }
  else if(corr_type==CCL_CORR_GL) {
    for (j=2;j<=ell_max;j++) {//https://arxiv.org/pdf/1007.4809.pdf
      Pl_theta[j]=gsl_sf_legendre_Plm(j,2,cth);
      Pl_theta[j]*=(2*j+1.)/((j+0.)*(j+1.));
    }
  }
}

/*--------ROUTINE: ccl_tracer_corr_legendre ------
TASK: Compute correlation function via Legendre polynomials
INPUT: cosmology, number of theta bins, theta array, tracer 1, tracer 2, i_bessel, boolean
       for tapering, vector of tapering limits, correlation vector, angular_cl function.
 */
static void ccl_tracer_corr_legendre(ccl_cosmology *cosmo,
				     int n_ell,double *ell,double *cls,
				     int n_theta,double *theta,double *wtheta,
				     int corr_type,int do_taper_cl,double *taper_cl_limits,
				     int *status)
{
  int i;
  double *l_arr,*cl_arr,*Pl_theta;
  SplPar *cl_spl;

  if(corr_type==CCL_CORR_LM || corr_type==CCL_CORR_LP){
    *status=CCL_ERROR_NOT_IMPLEMENTED;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: CCL does not support full-sky xi+- calcuations.\nhttps://arxiv.org/abs/1702.05301 indicates flat-sky to be sufficient.\n");
  }

  if(*status==0) {
    l_arr=malloc(((int)(ccl_splines->ELL_MAX_CORR)+1)*sizeof(double));
    if(l_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    cl_arr=malloc(((int)(ccl_splines->ELL_MAX_CORR)+1)*sizeof(double));
    if(cl_arr==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    //Interpolate input Cl into
    cl_spl=ccl_spline_init(n_ell,ell,cls,cls[0],0);
    if(cl_spl==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    double cl_tilt,l_edge,cl_edge;
    l_edge=ell[n_ell-1];
    if((cls[n_ell-1]*cls[n_ell-2]<0) || (cls[n_ell-2]==0)) {
      cl_tilt=0;
      cl_edge=0;
    }
    else {
      cl_tilt=log(cls[n_ell-1]/cls[n_ell-2])/log(ell[n_ell-1]/ell[n_ell-2]);
      cl_edge=cls[n_ell-1];
    }
    for(i=0;i<=(int)(ccl_splines->ELL_MAX_CORR);i++) {
      double l=(double)i;
      l_arr[i]=l;
      if(l>=l_edge)
	cl_arr[i]=cl_edge*pow(l/l_edge,cl_tilt);
      else
	cl_arr[i]=ccl_spline_eval(l,cl_spl);
    }
    ccl_spline_free(cl_spl);

    if (do_taper_cl)
      *status=taper_cl((int)(ccl_splines->ELL_MAX_CORR)+1,l_arr,cl_arr,taper_cl_limits);
  }

  if(*status==0) {
    Pl_theta=malloc(sizeof(double)*((int)(ccl_splines->ELL_MAX_CORR)+1));
    if(Pl_theta==NULL) {
      *status=CCL_ERROR_MEMORY;
      ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_tracer_corr_legendre ran out of memory\n");
    }
  }

  if(*status==0) {
    for (int i=0;i<n_theta;i++) {
      wtheta[i]=0;
      ccl_compute_legendre_polynomial(corr_type,theta[i],(int)(ccl_splines->ELL_MAX_CORR),Pl_theta);
      for(int i_L=1;i_L<(int)(ccl_splines->ELL_MAX_CORR);i_L+=1)
	wtheta[i]+=cl_arr[i_L]*Pl_theta[i_L];
      wtheta[i]/=(M_PI*4);
    }
  }

  free(Pl_theta);
  free(l_arr);
  free(cl_arr);
}

/*--------ROUTINE: ccl_tracer_corr ------
TASK: For a given tracer, get the correlation function. Do so by running
      ccl_angular_cls. If you already have Cls calculated, go to the next
      function to pass them directly.
INPUT: cosmology, number of theta values to evaluate = NL, theta vector,
       tracer 1, tracer 2, i_bessel, key for tapering, limits of tapering
       correlation function.
 */
void ccl_correlation(ccl_cosmology *cosmo,
		     int n_ell,double *ell,double *cls,
		     int n_theta,double *theta,double *wtheta,
		     int corr_type,int do_taper_cl,double *taper_cl_limits,int flag_method,
		     int *status)
{
  switch(flag_method) {
  case CCL_CORR_FFTLOG :
    ccl_tracer_corr_fftlog(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
			   do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_LGNDRE :
    ccl_tracer_corr_legendre(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,
			     do_taper_cl,taper_cl_limits,status);
    break;
  case CCL_CORR_BESSEL :
    ccl_tracer_corr_bessel(cosmo,n_ell,ell,cls,n_theta,theta,wtheta,corr_type,status);
    break;
  default :
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation. Unknown algorithm\n");
  }

  ccl_check_status(cosmo,status);
}

/*--------ROUTINE: ccl_correlation_3d ------
TASK: Calculate the 3d-correlation function. Do so by using FFTLog.

INPUT: cosmology, scale factor a,
       number of r values, r values,
       key for tapering, limits of tapering

Correlation function result will be in array xi
 */

void ccl_correlation_3d(ccl_cosmology *cosmo, double a,
			int n_r,double *r,double *xi,
			int do_taper_pk,double *taper_pk_limits,
			int *status)
{
  int i,N_ARR;
  double *k_arr,*pk_arr,*r_arr,*xi_arr;

  //number of data points for k and pk array
  N_ARR=(int)(ccl_splines->N_K_3DCOR*log10(ccl_splines->K_MAX/ccl_splines->K_MIN));

  k_arr=ccl_log_spacing(ccl_splines->K_MIN,ccl_splines->K_MAX,N_ARR);
  if(k_arr==NULL) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  pk_arr=malloc(N_ARR*sizeof(double));
  if(pk_arr==NULL) {
    free(k_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  for (i=0; i<N_ARR; i++)
    pk_arr[i] = ccl_nonlin_matter_power(cosmo, k_arr[i], a, status);

  if (do_taper_pk)
    taper_cl(N_ARR,k_arr,pk_arr,taper_pk_limits);

  r_arr=malloc(sizeof(double)*N_ARR);
  if(r_arr==NULL) {
    free(k_arr);
    free(pk_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }
  xi_arr=malloc(sizeof(double)*N_ARR);
  if(xi_arr==NULL) {
    free(k_arr); free(pk_arr); free(r_arr);
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_correlation.c: ccl_correlation_3d ran out of memory\n");
    return;
  }

  for(i=0;i<N_ARR;i++)
    r_arr[i]=0;

  pk2xi(N_ARR,k_arr,pk_arr,r_arr,xi_arr);

  // Interpolate to output values of r
  SplPar *xi_spl=ccl_spline_init(N_ARR,r_arr,xi_arr,xi_arr[0],0);
  for(i=0;i<n_r;i++)
    xi[i]=ccl_spline_eval(r[i],xi_spl);
  ccl_spline_free(xi_spl);

  free(k_arr); free(pk_arr);
  free(r_arr); free(xi_arr);

  ccl_check_status(cosmo,status);

  return;
}

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
#ifndef GENERIC_API_NAMESPACE
	#error "Wrong usage of this header"
#endif

#ifndef GENERIC_API_ROUTINES_NAMESPACE
	#error "Wrong usage of this header"
#endif

#define ROUTINES_NAMESPACE namespace GENERIC_API_ROUTINES_NAMESPACE
#define BLAS_NAMESPACE namespace GENERIC_API_NAMESPACE

#ifndef GENERIC_API_DEFINE

namespace cl
{
	namespace routines
	{
		ROUTINES_NAMESPACE
		{
			template<MathDomain md>
			static void Add(MemoryBuffer&, const MemoryBuffer&, const MemoryBuffer&, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void AddEqual(MemoryBuffer&, const MemoryBuffer&, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void AddEqualMatrix(MemoryTile&, const MemoryTile&, const MatrixOperation, const MatrixOperation, const double, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void Scale(MemoryBuffer&, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ScaleColumns(MemoryTile&, const MemoryBuffer&)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ElementwiseProduct(MemoryBuffer&, const MemoryBuffer&, const MemoryBuffer&, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void SubMultiply(MemoryTile&, const MemoryTile&, const MemoryTile&, const unsigned, const unsigned, const unsigned, const MatrixOperation, const MatrixOperation, const double, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void BatchedMultiply(MemoryCube&, const MemoryCube&, const MemoryCube&, const unsigned, const unsigned, const MatrixOperation, const MatrixOperation, const double, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void Dot(MemoryBuffer&, const MemoryTile&, const MemoryBuffer&, const MatrixOperation, const double, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void KroneckerProduct(MemoryTile&, const MemoryBuffer&, const MemoryBuffer&, const double)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void Solve(const MemoryTile&, MemoryTile&, const MatrixOperation, const LinearSystemSolverType)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ArgAbsMin(int&, const MemoryBuffer&)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ColumnWiseArgAbsMin(MemoryBuffer&, const MemoryTile&)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ArgAbsMax(int&, const MemoryBuffer&)
			{
				throw NotImplementedException();
			}

			template<MathDomain md>
			static void ColumnWiseArgAbsMax(MemoryBuffer&, const MemoryTile&)
			{
				throw NotImplementedException();
			}

			// norm = ||x||_2
			template<MathDomain md>
			static void EuclideanNorm(double&, const MemoryBuffer&)
			{
				throw NotImplementedException();
			}
		}
	}	 // namespace routines
}	 // namespace cl

#else

	#include <cmath>
	#include <complex>

BLAS_NAMESPACE
{
	#include <cblas.h>
	#include <lapacke.h>
}

namespace cl
{
	namespace routines
	{
		ROUTINES_NAMESPACE
		{
			static constexpr GENERIC_API_NAMESPACE::CBLAS_ORDER columnMajorLayout = { GENERIC_API_NAMESPACE::CBLAS_ORDER::CblasColMajor };
			static constexpr std::array<GENERIC_API_NAMESPACE::CBLAS_TRANSPOSE, 2> operationsEnum = { GENERIC_API_NAMESPACE::CBLAS_TRANSPOSE::CblasNoTrans, GENERIC_API_NAMESPACE::CBLAS_TRANSPOSE::CblasTrans };

			template<MathDomain md>
			static void Add(MemoryBuffer & z, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha);

			template<>
			inline void Add<MathDomain::Float>(MemoryBuffer & z, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha)
			{
				Copy<MathDomain::Float>(z, y);
				GENERIC_API_NAMESPACE::cblas_saxpy(static_cast<int>(z.size), static_cast<float>(alpha), reinterpret_cast<float*>(x.pointer), 1, reinterpret_cast<float*>(z.pointer), 1);
			}
			template<>
			inline void Add<MathDomain::Double>(MemoryBuffer & z, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha)
			{
				Copy<MathDomain::Double>(z, y);
				GENERIC_API_NAMESPACE::cblas_daxpy(static_cast<int>(z.size), alpha, reinterpret_cast<double*>(x.pointer), 1, reinterpret_cast<double*>(z.pointer), 1);
			}

			template<MathDomain md>
			static void AddEqual(MemoryBuffer & z, const MemoryBuffer& x, const double alpha);

			template<>
			inline void AddEqual<MathDomain::Float>(MemoryBuffer & z, const MemoryBuffer& x, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_saxpy(static_cast<int>(z.size), static_cast<float>(alpha), reinterpret_cast<float*>(x.pointer), 1, reinterpret_cast<float*>(z.pointer), 1);
			}
			template<>
			inline void AddEqual<MathDomain::Double>(MemoryBuffer & z, const MemoryBuffer& x, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_daxpy(static_cast<int>(z.size), alpha, reinterpret_cast<double*>(x.pointer), 1, reinterpret_cast<double*>(z.pointer), 1);
			}

			template<MathDomain md>
			static void AddEqualMatrix(MemoryTile & A, const MemoryTile& B, const MatrixOperation aOperation, const MatrixOperation bOperation, const double alpha, const double beta);

			template<>
			inline void AddEqualMatrix<MathDomain::Float>(MemoryTile & A, const MemoryTile& B, const MatrixOperation aOperation, const MatrixOperation bOperation, const double alpha, const double beta)
			{
				if (aOperation != MatrixOperation::None)
					throw NotImplementedException();
				if (bOperation != MatrixOperation::None)
					throw NotImplementedException();

				GENERIC_API_NAMESPACE::cblas_sgeadd(columnMajorLayout, static_cast<int>(A.nRows), static_cast<int>(A.nCols), static_cast<float>(beta), reinterpret_cast<float*>(B.pointer), static_cast<int>(B.leadingDimension), static_cast<float>(alpha), reinterpret_cast<float*>(A.pointer), static_cast<int>(A.leadingDimension));
			}
			template<>
			inline void AddEqualMatrix<MathDomain::Double>(MemoryTile & A, const MemoryTile& B, const MatrixOperation aOperation, const MatrixOperation bOperation, const double alpha, const double beta)
			{
				if (aOperation != MatrixOperation::None)
					throw NotImplementedException();
				if (bOperation != MatrixOperation::None)
					throw NotImplementedException();

				GENERIC_API_NAMESPACE::cblas_dgeadd(columnMajorLayout, static_cast<int>(A.nRows), static_cast<int>(A.nCols), beta, reinterpret_cast<double*>(B.pointer), static_cast<int>(B.leadingDimension), alpha, reinterpret_cast<double*>(A.pointer), static_cast<int>(A.leadingDimension));
			}

			template<MathDomain md>
			static void Scale(MemoryBuffer & z, const double alpha);

			template<>
			inline void Scale<MathDomain::Float>(MemoryBuffer & z, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_sscal(static_cast<int>(z.size), static_cast<float>(alpha), reinterpret_cast<float*>(z.pointer), 1);
			}
			template<>
			inline void Scale<MathDomain::Double>(MemoryBuffer & z, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_dscal(static_cast<int>(z.size), alpha, reinterpret_cast<double*>(z.pointer), 1);
			}

			template<MathDomain md>
			static void ScaleColumns(MemoryTile & z, const MemoryBuffer& alpha);

			template<>
			inline void ScaleColumns<MathDomain::Float>(MemoryTile & z, const MemoryBuffer& alpha)
			{
				for (size_t i = 0; i < z.nCols; ++i)
					GENERIC_API_NAMESPACE::cblas_sscal(static_cast<int>(z.nRows), *reinterpret_cast<float*>(alpha.pointer + i * alpha.ElementarySize()), reinterpret_cast<float*>(z.pointer + i * z.nRows * z.ElementarySize()), 1);
			}
			template<>
			inline void ScaleColumns<MathDomain::Double>(MemoryTile & z, const MemoryBuffer& alpha)
			{
				for (size_t i = 0; i < z.nCols; ++i)
					GENERIC_API_NAMESPACE::cblas_dscal(static_cast<int>(z.nRows), *reinterpret_cast<double*>(alpha.pointer + i * alpha.ElementarySize()), reinterpret_cast<double*>(z.pointer + i * z.nRows * z.ElementarySize()), 1);
			}

			template<MathDomain md>
			static void SubMultiply(MemoryTile & A, const MemoryTile& B, const MemoryTile& C, const unsigned nRowsB, const unsigned nColsB, const unsigned nColsC, const MatrixOperation bOperation, const MatrixOperation cOperation, const double alpha, const double beta);

			template<>
			inline void SubMultiply<MathDomain::Float>(MemoryTile & A, const MemoryTile& B, const MemoryTile& C, const unsigned nRowsB, const unsigned nColsB, const unsigned nColsC, const MatrixOperation bOperation, const MatrixOperation cOperation, const double alpha, const double beta)
			{
				GENERIC_API_NAMESPACE::cblas_sgemm(columnMajorLayout, operationsEnum[static_cast<unsigned>(bOperation)], operationsEnum[static_cast<unsigned>(cOperation)], static_cast<int>(nRowsB), static_cast<int>(nColsC), static_cast<int>(nColsB), static_cast<float>(alpha), reinterpret_cast<float*>(B.pointer), static_cast<int>(B.leadingDimension), reinterpret_cast<float*>(C.pointer), static_cast<int>(C.leadingDimension), static_cast<float>(beta), reinterpret_cast<float*>(A.pointer), static_cast<int>(A.leadingDimension));
			}

			template<>
			inline void SubMultiply<MathDomain::Double>(MemoryTile & A, const MemoryTile& B, const MemoryTile& C, const unsigned nRowsB, const unsigned nColsB, const unsigned nColsC, const MatrixOperation bOperation, const MatrixOperation cOperation, const double alpha, const double beta)
			{
				GENERIC_API_NAMESPACE::cblas_dgemm(columnMajorLayout, operationsEnum[static_cast<unsigned>(bOperation)], operationsEnum[static_cast<unsigned>(cOperation)], static_cast<int>(nRowsB), static_cast<int>(nColsC), static_cast<int>(nColsB), alpha, reinterpret_cast<double*>(B.pointer), static_cast<int>(B.leadingDimension), reinterpret_cast<double*>(C.pointer), static_cast<int>(C.leadingDimension), beta, reinterpret_cast<double*>(A.pointer), static_cast<int>(A.leadingDimension));
			}

			template<MathDomain md>
			static void Dot(MemoryBuffer & y, const MemoryTile& A, const MemoryBuffer& x, const MatrixOperation aOperation, const double alpha = 1.0, const double beta = 0.0);

			template<>
			inline void Dot<MathDomain::Float>(MemoryBuffer & y, const MemoryTile& A, const MemoryBuffer& x, const MatrixOperation aOperation, const double alpha, const double beta)
			{
				GENERIC_API_NAMESPACE::cblas_sgemv(columnMajorLayout, operationsEnum[static_cast<unsigned>(aOperation)], static_cast<int>(A.nRows), static_cast<int>(A.nCols), static_cast<float>(alpha), reinterpret_cast<float*>(A.pointer), static_cast<int>(A.leadingDimension), reinterpret_cast<float*>(x.pointer), 1, static_cast<float>(beta), reinterpret_cast<float*>(y.pointer), 1);
			}
			template<>
			inline void Dot<MathDomain::Double>(MemoryBuffer & y, const MemoryTile& A, const MemoryBuffer& x, const MatrixOperation aOperation, const double alpha, const double beta)
			{
				GENERIC_API_NAMESPACE::cblas_dgemv(columnMajorLayout, operationsEnum[static_cast<unsigned>(aOperation)], static_cast<int>(A.nRows), static_cast<int>(A.nCols), alpha, reinterpret_cast<double*>(A.pointer), static_cast<int>(A.leadingDimension), reinterpret_cast<double*>(x.pointer), 1, beta, reinterpret_cast<double*>(y.pointer), 1);
			}

			template<MathDomain md>
			static void KroneckerProduct(MemoryTile & A, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha);

			template<>
			inline void KroneckerProduct<MathDomain::Float>(MemoryTile & A, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_sger(columnMajorLayout, static_cast<int>(x.size), static_cast<int>(y.size), static_cast<float>(alpha), reinterpret_cast<float*>(x.pointer), 1, reinterpret_cast<float*>(y.pointer), 1, reinterpret_cast<float*>(A.pointer), static_cast<int>(A.nRows));
			}
			template<>
			inline void KroneckerProduct<MathDomain::Double>(MemoryTile & A, const MemoryBuffer& x, const MemoryBuffer& y, const double alpha)
			{
				GENERIC_API_NAMESPACE::cblas_dger(columnMajorLayout, static_cast<int>(x.size), static_cast<int>(y.size), alpha, reinterpret_cast<double*>(x.pointer), 1, reinterpret_cast<double*>(y.pointer), 1, reinterpret_cast<double*>(A.pointer), static_cast<int>(A.nRows));
			}

			template<MathDomain md>
			static void Solve(const MemoryTile& A, MemoryTile& B, const MatrixOperation aOperation, const LinearSystemSolverType solver);

			template<>
			inline void Solve<MathDomain::Float>(const MemoryTile& A, MemoryTile& B, const MatrixOperation aOperation, const LinearSystemSolverType solver)
			{
				// Need to copy A, as it will be overwritten by its factorization
				MemoryTile aCopy(A);
				Alloc(aCopy);
				Copy<MathDomain::Float>(aCopy, A);

				const auto nra = static_cast<int>(A.nRows);
				const auto lda = static_cast<int>(A.leadingDimension);
				const auto ncb = static_cast<int>(B.nCols);
				const auto ldb = static_cast<int>(B.leadingDimension);

				// Initializes auxliary value for solver
				int info = 0;

				switch (solver)
				{
					case LinearSystemSolverType::Lu:
					{
						// allocate memory for pivoting
						MemoryBuffer pivot(0, A.nRows, A.memorySpace, MathDomain::Int);
						Alloc(pivot);

						// Factorize A (and overwrite it with L)
						info = GENERIC_API_NAMESPACE::LAPACKE_sgetrf(static_cast<int>(columnMajorLayout), nra, nra, reinterpret_cast<float*>(aCopy.pointer), lda, reinterpret_cast<int*>(pivot.pointer));
						if (info != 0)
							throw OpenBlasException(__func__);

						// Solve factorized system
						info = GENERIC_API_NAMESPACE::LAPACKE_sgetrs(static_cast<int>(columnMajorLayout), openBlasOperation[static_cast<unsigned>(aOperation)], nra, ncb, reinterpret_cast<float*>(aCopy.pointer), lda, reinterpret_cast<int*>(pivot.pointer), reinterpret_cast<float*>(B.pointer), ldb);
						if (info != 0)
							throw OpenBlasException(__func__);

						// free memory
						Free(pivot);

						break;
					}
					case LinearSystemSolverType::Qr:
					{
						// allocate memory for tau
						MemoryBuffer tau(0, A.nRows, A.memorySpace, MathDomain::Float);
						Alloc(tau);

						// A = Q * R
						/* int matrix_layout, lapack_int m, lapack_int n,
								   float* a, lapack_int lda, float* tau */
						info = GENERIC_API_NAMESPACE::LAPACKE_sgeqrf(static_cast<int>(columnMajorLayout), nra, nra, reinterpret_cast<float*>(aCopy.pointer), lda, reinterpret_cast<float*>(tau.pointer));
						if (info != 0)
							throw OpenBlasException(__func__);

						// B = Q^T * B
						/*
						 *  int matrix_layout, char side, char trans,
								   lapack_int m, lapack_int n, lapack_int k,
								   const float* a, lapack_int lda, const float* tau,
								   float* c, lapack_int ldc */
						info = GENERIC_API_NAMESPACE::LAPACKE_sormqr(static_cast<int>(columnMajorLayout), 'L', 'T', nra, nra, ncb, reinterpret_cast<float*>(aCopy.pointer), lda, reinterpret_cast<float*>(tau.pointer), reinterpret_cast<float*>(B.pointer), ldb);
						if (info != 0)
							throw OpenBlasException(__func__);

						// Solve (x = R \ (Q^T * B))
						GENERIC_API_NAMESPACE::cblas_strsm(columnMajorLayout, GENERIC_API_NAMESPACE::CBLAS_SIDE::CblasLeft, GENERIC_API_NAMESPACE::CBLAS_UPLO::CblasUpper, operationsEnum[static_cast<unsigned>(aOperation)], GENERIC_API_NAMESPACE::CBLAS_DIAG::CblasNonUnit, nra, nra, 1.0, reinterpret_cast<float*>(aCopy.pointer), lda, reinterpret_cast<float*>(B.pointer), ldb);

						// free memory
						Free(tau);

						break;
					}
					default:
						throw NotImplementedException();
				}

				Free(aCopy);
			}

			template<>
			inline void Solve<MathDomain::Double>(const MemoryTile& A, MemoryTile& B, const MatrixOperation aOperation, const LinearSystemSolverType solver)
			{
				// Need to copy A, as it will be overwritten by its factorization
				MemoryTile aCopy(A);
				Alloc(aCopy);
				Copy<MathDomain::Double>(aCopy, A);

				const auto nra = static_cast<int>(A.nRows);
				const auto lda = static_cast<int>(A.leadingDimension);
				const auto ncb = static_cast<int>(B.nCols);
				const auto ldb = static_cast<int>(B.leadingDimension);

				// Initializes auxliary value for solver
				int info = 0;

				switch (solver)
				{
					case LinearSystemSolverType::Lu:
					{
						// allocate memory for pivoting
						MemoryBuffer pivot(0, A.nRows, A.memorySpace, MathDomain::Int);
						Alloc(pivot);

						// Factorize A (and overwrite it with L)
						info = GENERIC_API_NAMESPACE::LAPACKE_dgetrf(static_cast<int>(columnMajorLayout), nra, nra, reinterpret_cast<double*>(aCopy.pointer), lda, reinterpret_cast<int*>(pivot.pointer));
						if (info != 0)
							throw OpenBlasException(__func__);

						// Solve factorized system
						info = GENERIC_API_NAMESPACE::LAPACKE_dgetrs(static_cast<int>(columnMajorLayout), openBlasOperation[static_cast<unsigned>(aOperation)], nra, ncb, reinterpret_cast<double*>(aCopy.pointer), lda, reinterpret_cast<int*>(pivot.pointer), reinterpret_cast<double*>(B.pointer), ldb);
						if (info != 0)
							throw OpenBlasException(__func__);

						// free memory
						Free(pivot);

						break;
					}
					case LinearSystemSolverType::Qr:
					{
						// allocate memory for tau
						MemoryBuffer tau(0, A.nRows, A.memorySpace, MathDomain::Double);
						Alloc(tau);

						// A = Q * R
						/* int matrix_layout, lapack_int m, lapack_int n,
								   double* a, lapack_int lda, double* tau */
						info = GENERIC_API_NAMESPACE::LAPACKE_dgeqrf(static_cast<int>(columnMajorLayout), nra, nra, reinterpret_cast<double*>(aCopy.pointer), lda, reinterpret_cast<double*>(tau.pointer));
						if (info != 0)
							throw OpenBlasException(__func__);

						// B = Q^T * B
						/*
						 *  int matrix_layout, char side, char trans,
								   lapack_int m, lapack_int n, lapack_int k,
								   const double* a, lapack_int lda, const double* tau,
								   double* c, lapack_int ldc */
						info = GENERIC_API_NAMESPACE::LAPACKE_dormqr(static_cast<int>(columnMajorLayout), 'L', 'T', nra, nra, ncb, reinterpret_cast<double*>(aCopy.pointer), lda, reinterpret_cast<double*>(tau.pointer), reinterpret_cast<double*>(B.pointer), ldb);
						if (info != 0)
							throw OpenBlasException(__func__);

						// Solve (x = R \ (Q^T * B))
						cblas_dtrsm(columnMajorLayout, GENERIC_API_NAMESPACE::CBLAS_SIDE::CblasLeft, GENERIC_API_NAMESPACE::CBLAS_UPLO::CblasUpper, operationsEnum[static_cast<unsigned>(aOperation)], GENERIC_API_NAMESPACE::CBLAS_DIAG::CblasNonUnit, nra, nra, 1.0, reinterpret_cast<double*>(aCopy.pointer), lda, reinterpret_cast<double*>(B.pointer), ldb);

						// free memory
						Free(tau);

						break;
					}
					default:
						throw NotImplementedException();
				}

				Free(aCopy);
			}

			template<MathDomain md>
			static void ArgAbsMin(int& argMin, const MemoryBuffer& x);

			template<>
			inline void ArgAbsMin<MathDomain::Float>(int&, const MemoryBuffer&)
			{
				// TODO: apparently there's isamax but not isamin??
				// argMin = static_cast<int>(cblas_isamin(static_cast<int>(x.size), reinterpret_cast<float*>(x.pointer), 1));
			}
			template<>
			inline void ArgAbsMin<MathDomain::Double>(int&, const MemoryBuffer&)
			{
				// TODO: apparently there's isamax but not isamin??
				// argMin = static_cast<int>(cblas_idamin(static_cast<int>(x.size), reinterpret_cast<double*>(x.pointer), 1));
			}

			template<MathDomain md>
			static void ColumnWiseArgAbsMin(MemoryBuffer & argMin, const MemoryTile& A);

			template<>
			inline void ColumnWiseArgAbsMin<MathDomain::Float>(MemoryBuffer&, const MemoryTile&)
			{
				// TODO
				//		auto* argMinPtr = reinterpret_cast<int*>(argMin.pointer);
				//		// 1 added for compatibility
				//		for (size_t j = 0; j < A.nCols; ++j)
				//			argMinPtr[j] = 1 + static_cast<int>(cblas_isamin(static_cast<int>(A.nRows), reinterpret_cast<float*>(A.pointer + j * A.nRows * A.ElementarySize()), 1));
			}
			template<>
			inline void ColumnWiseArgAbsMin<MathDomain::Double>(MemoryBuffer&, const MemoryTile&)
			{
				// TODO
				//		// 1 added for compatibility
				//		auto* argMinPtr = reinterpret_cast<int*>(argMin.pointer);
				//		for (size_t j = 0; j < A.nCols; ++j)
				//			argMinPtr[j] = 1 + static_cast<int>(cblas_idamin(static_cast<int>(A.nRows), reinterpret_cast<double*>(A.pointer + j * A.nRows * A.ElementarySize()), 1));
			}

			template<MathDomain md>
			static void ArgAbsMax(int& argMax, const MemoryBuffer& x);

			template<>
			inline void ArgAbsMax<MathDomain::Float>(int& argMax, const MemoryBuffer& x)
			{
				argMax = static_cast<int>(GENERIC_API_NAMESPACE::cblas_isamax(static_cast<int>(x.size), reinterpret_cast<float*>(x.pointer), 1));
			}
			template<>
			inline void ArgAbsMax<MathDomain::Double>(int& argMax, const MemoryBuffer& x)
			{
				argMax = static_cast<int>(GENERIC_API_NAMESPACE::cblas_idamax(static_cast<int>(x.size), reinterpret_cast<double*>(x.pointer), 1));
			}

			template<MathDomain md>
			static void ColumnWiseArgAbsMax(MemoryBuffer & argMax, const MemoryTile& A);

			template<>
			inline void ColumnWiseArgAbsMax<MathDomain::Float>(MemoryBuffer & argMax, const MemoryTile& A)
			{
				// 1 added for compatibility
				auto* argMaxPtr = reinterpret_cast<int*>(argMax.pointer);
				for (size_t j = 0; j < A.nCols; ++j)
					argMaxPtr[j] = 1 + static_cast<int>(GENERIC_API_NAMESPACE::cblas_isamax(static_cast<int>(A.nRows), reinterpret_cast<float*>(A.pointer + j * A.nRows * A.ElementarySize()), 1));
			}
			template<>
			inline void ColumnWiseArgAbsMax<MathDomain::Double>(MemoryBuffer & argMax, const MemoryTile& A)
			{
				// 1 added for compatibility
				auto* argMaxPtr = reinterpret_cast<int*>(argMax.pointer);
				for (size_t j = 0; j < A.nCols; ++j)
					argMaxPtr[j] = 1 + static_cast<int>(GENERIC_API_NAMESPACE::cblas_idamax(static_cast<int>(A.nRows), reinterpret_cast<double*>(A.pointer + j * A.nRows * A.ElementarySize()), 1));
			}

			// norm = ||x||_2
			template<MathDomain md>
			static void EuclideanNorm(double& norm, const MemoryBuffer& z);

			template<>
			inline void EuclideanNorm<MathDomain::Float>(double& norm, const MemoryBuffer& z)
			{
				norm = static_cast<double>(GENERIC_API_NAMESPACE::cblas_snrm2(static_cast<int>(z.size), reinterpret_cast<float*>(z.pointer), 1));
			}
			template<>
			inline void EuclideanNorm<MathDomain::Double>(double& norm, const MemoryBuffer& z)
			{
				norm = GENERIC_API_NAMESPACE::cblas_dnrm2(static_cast<int>(z.size), reinterpret_cast<double*>(z.pointer), 1);
			}
		}
	}	 // namespace routines
}	 // namespace cl

#endif

#undef GENERIC_API_DEFINE
#undef GENERIC_API_NAMESPACE
#undef GENERIC_API_ROUTINES_NAMESPACE
#undef ROUTINES_NAMESPACE
#undef BLAS_NAMESPACE

max_stars_repo_name: pmontalb/CudaLight

max_issues_repo_name: pmontalb/CudaLight

max_forks_repo_name: pmontalb/CudaLight