text:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include <string.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_randist.h>

#include "allvars.h"
#include "proto.h"

int mcmc_sampling(char *fname_mcmc, double (* prob_fun)(double *theta))
{
  const int n_cov_update=10000;
  double prob, prob_new, scale, ratio, u;
  double theta[ntheta], theta_new[ntheta];
  double var_cov_mat[ntheta*ntheta], Pmatrix[ntheta*ntheta], Prandvec[ntheta], Py[ntheta];
  double *theta_mcmc;
  int istep, iaccpt, i, j, info, flag;

  FILE *fmcmc;
  fmcmc = fopen(fname_mcmc, "w");
  theta_mcmc = array_malloc(ntheta*n_cov_update);

/*
 * sampling in logrithm space.
 */
  memcpy(theta, theta_input, ntheta*sizeof(double));

  for(i=0; i<ntheta; i++)
  {
    var_cov_mat[i*ntheta+i] = sigma_input[i] * sigma_input[i];
    for(j=0; j<i; j++)
      var_cov_mat[i*ntheta+j] = var_cov_mat[j*ntheta + i] = 0.0;
  }
  memcpy(Pmatrix, var_cov_mat, ntheta*ntheta*sizeof(double));
  Chol_decomp_U(Pmatrix, ntheta, &info);

  scale = 1.0;

  prob = prob_fun(theta);
  istep = 0;
  iaccpt = 0;

  printf("step:%6d", 0);
  while(istep < nmcmc)
  {
    flag = 1;
    while(flag)
    {
      for(i=0; i<ntheta; i++)
      {
        Prandvec[i]  = gsl_ran_gaussian(gsl_r, 1.0);
      }
      multiply_matvec(Pmatrix, Prandvec, ntheta, Py);
      flag = 0;
      for(i=0; i<ntheta; i++)
      {
        if(theta_fixed[i]==1)
        {
          theta_new[i] = theta[i];
        }
        else
        {
          theta_new[i] = theta[i] + scale * Py[i];
          if(theta_new[i]<theta_range[i][0] || theta_new[i]>theta_range[i][1])
          {
            flag = 1;
            break;
          }
        }
      }
    }

    prob_new = prob_fun(theta_new);

    ratio = prob_new - prob;
    if(ratio > 0.0)
    {
      memcpy(theta, theta_new, ntheta*sizeof(double));
      prob = prob_new;
      iaccpt++;
    }
    else
    {
      u = gsl_rng_uniform(gsl_r);
      if(log(u) < ratio)
      {
        memcpy(theta, theta_new, ntheta*sizeof(double));
        prob = prob_new;
        iaccpt++;
      }
    }

    for(i=0; i<ntheta; i++)
      theta_mcmc[i*n_cov_update + istep%n_cov_update] = theta[i];
    if(istep%10==0)printf("\b\b\b\b\b\b%6d", istep);

    fprintf(fmcmc, "%d", istep);
    for(i=0;i<ntheta;i++)
    {
      fprintf(fmcmc,"\t%f", theta[i]);
    }
    fprintf(fmcmc, "\t%e\n", prob);

    istep++;

    if(istep%n_cov_update == 0)
    {
      get_cov_matrix_diag(theta_mcmc, n_cov_update, ntheta);
      memcpy(var_cov_mat, cov_matrix, ntheta*ntheta*sizeof(double));
      //display_mat(var_cov_mat, ntheta, ntheta);
      memcpy(Pmatrix, var_cov_mat, ntheta*ntheta*sizeof(double));
      Chol_decomp_U(Pmatrix, ntheta, &info);
    }
  }

  printf("\n");
  printf("Accept rate: %f\n", 1.0*iaccpt/istep);

  fclose(fmcmc);
  free(theta_mcmc);
  return 0;
}

ext: c

max_stars_repo_name: LiyrAstroph/MICA

max_issues_repo_name: LiyrAstroph/MICA

max_forks_repo_name: LiyrAstroph/MICA


text:
/*
 *  Copyright 2008-2014 NVIDIA Corporation
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

#include <cblas.h>

#include <cusp/complex.h>

#define CUSP_CBLAS_EXPAND_REAL_DEFS(FUNC_MACRO)                                             \
  FUNC_MACRO(float , float , s)                                                             \
  FUNC_MACRO(double, double, d)

#define CUSP_CBLAS_AMAX(T,V,name)                                                           \
  template<int dummy>                                                                       \
  int amax( const int n, const T* X, const int incX )                                       \
  {                                                                                \
    return cblas_i##name##amax(n, (const V*) X, incX);                                      \
  }

#define CUSP_CBLAS_ASUM(T,V,name)                                                           \
  template<int dummy>                                                                       \
  V asum( const int n, const T* X, const int incX )                                         \
  {                                                                                \
    return cblas_##name##asum(n, (const V*) X, incX);                                       \
  }

#define CUSP_CBLAS_AXPY(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void axpy( const int n, const T alpha, const T* X, const int incX, T* Y, const int incY ) \
  {                                                                                \
    cblas_##name##axpy(n, alpha, (const V*) X, incX, (V*) Y, incY);                         \
  }

#define CUSP_CBLAS_COPY(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void copy( const int n, const T* X, const int incX, T* Y, const int incY )                \
  {                                                                                \
    cblas_##name##copy(n, (const V*) X, incX, (V*) Y, incY);                                \
  }

#define CUSP_CBLAS_DOT(T,V,name)                                                            \
  template<int dummy>                                                                       \
  T dot( const int n, const T* X, const int incX, const T* Y, const int incY )              \
  {                                                                                \
    return cblas_##name##dot(n, (const V*) X, incX, (const V*) Y, incY);                    \
  }

#define CUSP_CBLAS_NRM2(T,V,name)                                                           \
  template<int dummy>                                                                       \
  V nrm2( const int n, const T* X, const int incX )                                         \
  {                                                                                \
    return cblas_##name##nrm2(n, (const V*) X, incX);                                       \
  }

#define CUSP_CBLAS_SCAL(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void scal( const int n, const T alpha, T* X, const int incX )                             \
  {                                                                                \
    cblas_##name##scal(n, alpha, (V*) X, incX);                                             \
  }

#define CUSP_CBLAS_SWAP(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void swap( const int n, T* X, const int incX, T* Y, const int incY )                      \
  {                                                                                \
    cblas_##name##swap(n, (V*) X, incX, (V*) Y, incY);                                      \
  }

#define CUSP_CBLAS_GEMV(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void gemv(  CBLAS_ORDER order,  CBLAS_TRANSPOSE trans,                                    \
              int m, int n, T alpha, const T* A, int lda,                                   \
              const T* x, int incx, T beta, T* y, int incy)                                 \
{                                                                                \
    cblas_##name##gemv(order, trans, m, n, alpha, A, lda, x, incx, beta, y, incy);          \
}

#define CUSP_CBLAS_GER(T,V,name)                                                            \
  template<int dummy>                                                                       \
  void ger(  CBLAS_ORDER order, int m, int n, T alpha, const T* x, int incx,                \
            const T* y, int incy, T* A, int lda)                                            \
{                                                                                \
    cblas_##name##ger(order, m, n, alpha,                                                   \
                      (const V*) x, incx, (const V*) y, incy,                               \
                      (V*) A, lda);                                                         \
}

#define CUSP_CBLAS_SYMV(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void symv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
             int n, T alpha, const T* A, int lda,                                           \
             const T* x, int incx, T beta, T* y, int incy)                                  \
{                                                                                \
    cblas_##name##symv(order, uplo, n, alpha, (const V*) A, lda,                            \
                       (const V*) x, incx, beta, (V*) y, incy);                             \
}

#define CUSP_CBLAS_SYR(T,V,name)                                                            \
  template<int dummy>                                                                       \
  void syr(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                           \
            int n, T alpha, const T* x, int incx, T* A, int lda)                            \
{                                                                                \
    cblas_##name##syr(order, uplo, n, alpha,                                                \
                      (const V*) x, incx, (V*) A, lda);                                     \
}

#define CUSP_CBLAS_TRMV(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void trmv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int n, const T* A, int lda, T* x, int incx)                                    \
{                                                                                \
    cblas_##name##trmv(order, uplo, trans, diag, n,                                         \
                       (const V*) A, lda, (V*) x, incx);                                    \
}

#define CUSP_CBLAS_TRSV(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void trsv(  CBLAS_ORDER order,  CBLAS_UPLO uplo,                                          \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int n, const T* A, int lda, T* x, int incx)                                    \
{                                                                                \
    cblas_##name##trsv(order, uplo, trans, diag, n,                                         \
                       (const V*) A, lda, (V*) x, incx);                                    \
}

#define CUSP_CBLAS_GEMM(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void gemm(  CBLAS_ORDER order,                                                            \
              CBLAS_TRANSPOSE transa,  CBLAS_TRANSPOSE transb,                              \
              int m, int n, int k, T alpha, const T* A, int lda,                            \
              const T* B, int ldb, T beta, T* C, int ldc)                                   \
{                                                                                \
    cblas_##name##gemm(order, transa, transb,                                               \
                       m, n, k, alpha, A, lda, B, ldb, beta, C, ldc);                       \
}

#define CUSP_CBLAS_SYMM(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void symm(  CBLAS_ORDER order,                                                            \
              CBLAS_SIDE side,  CBLAS_UPLO uplo,                                            \
             int m, int n, T alpha, const T* A, int lda,                                    \
             const T* B, int ldb, T beta, T* C, int ldc)                                    \
{                                                                                \
    cblas_##name##symm(order, side, uplo, m, n,                                             \
                       (V) alpha, (const V*) A, lda, (const V*) B, ldb,                     \
                       (V) beta, (V*) C, ldc);                                              \
}

#define CUSP_CBLAS_SYRK(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void syrk(  CBLAS_ORDER order,                                                            \
              CBLAS_UPLO uplo,  CBLAS_TRANSPOSE trans,                                      \
             int n, int k, T alpha, const T* A, int lda,                                    \
             T beta, T* C, int ldc)                                                         \
{                                                                                \
    cblas_##name##syrk(order, uplo, trans, n, k,                                            \
                       (V) alpha, (const V*) A, lda,                                        \
                       (V) beta, (V*) C, ldc);                                              \
}

#define CUSP_CBLAS_SYR2K(T,V,name)                                                          \
  template<int dummy>                                                                       \
  void syr2k(  CBLAS_ORDER order,                                                           \
               CBLAS_UPLO uplo,  CBLAS_TRANSPOSE trans,                                     \
              int n, int k, T& alpha, const T* A, int lda,                                  \
              const T* B, int ldb, T& beta, T* C, int ldc)                                  \
{                                                                                \
    cblas_##name##syr2k(order, uplo, trans, n, k,                                           \
                        alpha, (const V*) A, lda,                                           \
                        (const V*) B, ldb, beta, (V*) C, ldc);                              \
}

#define CUSP_CBLAS_TRMM(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void trmm(  CBLAS_ORDER order,                                                            \
              CBLAS_SIDE side,  CBLAS_UPLO uplo,                                            \
              CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                      \
             int m, int n, T alpha, const T* A, int lda,                                    \
             T* B, int ldb)                                                                 \
{                                                                                \
    cblas_##name##trmm(order, side, uplo, trans, diag, m, n,                                \
                       (V) alpha, (const V*) A, lda, (V*) B, ldb);                          \
}

#define CUSP_CBLAS_TRSM(T,V,name)                                                           \
  template<int dummy>                                                                       \
  void trsm( CBLAS_ORDER order,                                                             \
             CBLAS_SIDE side,  CBLAS_UPLO uplo,                                             \
             CBLAS_TRANSPOSE trans,  CBLAS_DIAG diag,                                       \
             int m, int n, T alpha, const T* A, int lda,                                    \
             T* B, int ldb)                                                                 \
{                                                                                \
    cblas_##name##trsm(order, side, uplo, trans, diag, m, n,                                \
                       (V) alpha, (const V*) A, lda, (V*) B, ldb);                          \
}

namespace cusp
{
namespace system
{
namespace cpp
{
namespace detail
{
namespace cblas
{

// LEVEL 1
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_AMAX);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_ASUM);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_AXPY);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_COPY);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_DOT);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_NRM2);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SCAL);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SWAP);

// LEVEL 2
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_GEMV);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_GER);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SYMV);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SYR);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_TRMV);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_TRSV);

// LEVEL 3
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_GEMM);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SYMM);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SYRK);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_SYR2K);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_TRMM);
CUSP_CBLAS_EXPAND_REAL_DEFS(CUSP_CBLAS_TRSM);

} // end namespace cblas
} // end namespace detail
} // end namespace cpp
} // end namespace system
} // end namespace cusp

ext: h

max_stars_repo_name: Raman-sh/cusplibrary


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

#include "tests.h"

void
test_syr (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha = 0.1f;
   float A[] = { -0.291f };
   float X[] = { 0.845f };
   int incX = -1;
   float A_expected[] = { -0.219597f };
   cblas_ssyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr(case 1402)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha = 0.1f;
   float A[] = { -0.291f };
   float X[] = { 0.845f };
   int incX = -1;
   float A_expected[] = { -0.219597f };
   cblas_ssyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr(case 1403)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha = 0.1f;
   float A[] = { -0.291f };
   float X[] = { 0.845f };
   int incX = -1;
   float A_expected[] = { -0.219597f };
   cblas_ssyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr(case 1404)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha = 0.1f;
   float A[] = { -0.291f };
   float X[] = { 0.845f };
   int incX = -1;
   float A_expected[] = { -0.219597f };
   cblas_ssyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], flteps, "ssyr(case 1405)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha = -0.3;
   double A[] = { -0.65 };
   double X[] = { -0.891 };
   int incX = -1;
   double A_expected[] = { -0.8881643 };
   cblas_dsyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr(case 1406)");
     }
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha = -0.3;
   double A[] = { -0.65 };
   double X[] = { -0.891 };
   int incX = -1;
   double A_expected[] = { -0.8881643 };
   cblas_dsyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr(case 1407)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha = -0.3;
   double A[] = { -0.65 };
   double X[] = { -0.891 };
   int incX = -1;
   double A_expected[] = { -0.8881643 };
   cblas_dsyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr(case 1408)");
     }
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha = -0.3;
   double A[] = { -0.65 };
   double X[] = { -0.891 };
   int incX = -1;
   double A_expected[] = { -0.8881643 };
   cblas_dsyr(order, uplo, N, alpha, X, incX, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[i], A_expected[i], dbleps, "dsyr(case 1409)");
     }
   };
  };

}


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

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#define M 107
#define N 239
#define M1 17
#define N1 23
#define MR 10
#define NR 5

int
main (void)
{
  double xr[MR + 1] =
  {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};

  double yr[NR + 1] =
  {90.0, 91.0, 92.0, 93.0, 94.0, 95.0};

  gsl_histogram2d *h, *h1, *g, *hr;
  size_t i, j, k;

  gsl_ieee_env_setup ();

  h = gsl_histogram2d_calloc (M, N);
  h1 = gsl_histogram2d_calloc (M, N);
  g = gsl_histogram2d_calloc (M, N);

  gsl_test (h->xrange == 0,
	    "gsl_histogram2d_calloc returns valid xrange pointer");
  gsl_test (h->yrange == 0,
	    "gsl_histogram2d_calloc returns valid yrange pointer");
  gsl_test (h->bin == 0,
	    "gsl_histogram2d_calloc returns valid bin pointer");

  gsl_test (h->nx != M, "gsl_histogram2d_calloc returns valid nx");
  gsl_test (h->ny != N, "gsl_histogram2d_calloc returns valid ny");

  hr = gsl_histogram2d_calloc_range (MR, NR, xr, yr);

  gsl_test (hr->xrange == 0,
	    "gsl_histogram2d_calloc_range returns valid xrange pointer");
  gsl_test (hr->yrange == 0,
	    "gsl_histogram2d_calloc_range returns valid yrange pointer");
  gsl_test (hr->bin == 0,
	    "gsl_histogram2d_calloc_range returns valid bin pointer");

  gsl_test (hr->nx != MR, "gsl_histogram2d_calloc_range returns valid nx");
  gsl_test (hr->ny != NR, "gsl_histogram2d_calloc_range returns valid ny");

  {
    int status = 0;
    for (i = 0; i <= MR; i++)
      {
	if (hr->xrange[i] != xr[i])
	  {
	    status = 1;
	  }
      };

    gsl_test (status, "gsl_histogram2d_calloc_range creates xrange correctly");
  }

  {
    int status = 0;
    for (i = 0; i <= NR; i++)
      {
	if (hr->yrange[i] != yr[i])
	  {
	    status = 1;
	  }
      };

    gsl_test (status, "gsl_histogram2d_calloc_range creates yrange correctly");
  }

  for (i = 0; i <= MR; i++)
    {
      hr->xrange[i] = 0.0;
    }

  for (i = 0; i <= NR; i++)
    {
      hr->yrange[i] = 0.0;
    }

  {
    int status = gsl_histogram2d_set_ranges (hr, xr, MR+1, yr, NR+1);

    for (i = 0; i <= MR; i++)
      {
	if (hr->xrange[i] != xr[i])
	  {
	    status = 1;
	  }
      };

    gsl_test (status, "gsl_histogram2d_set_ranges sets xrange correctly");
  }

  {
    int status = 0;
    for (i = 0; i <= NR; i++)
      {
	if (hr->yrange[i] != yr[i])
	  {
	    status = 1;
	  }
      };

    gsl_test (status, "gsl_histogram2d_set_ranges sets yrange correctly");
  }

  k = 0;
  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
	{
	  k++;
	  gsl_histogram2d_accumulate (h, (double) i, (double) j, (double) k);
	};
    }

  {
    int status = 0;
    k = 0;
    for (i = 0; i < M; i++)
      {
	for (j = 0; j < N; j++)
	  {
	    k++;
	    if (h->bin[i * N + j] != (double) k)
	      {
		status = 1;
	      }
	  }
      }
    gsl_test (status,
	      "gsl_histogram2d_accumulate writes into array correctly");
  }

  {
    int status = 0;
    k = 0;
    for (i = 0; i < M; i++)
      {
	for (j = 0; j < N; j++)
	  {
	    k++;
	    if (gsl_histogram2d_get (h, i, j) != (double) k)
	      status = 1;
	  };
      }
    gsl_test (status, "gsl_histogram2d_get reads from array correctly");
  }

  for (i = 0; i <= M; i++)
    {
      h1->xrange[i] = 100.0 + i;
    }

  for (i = 0; i <= N; i++)
    {
      h1->yrange[i] = 900.0 + i*i;
    }

  gsl_histogram2d_memcpy (h1, h);

  {
    int status = 0;
    for (i = 0; i <= M; i++)
      {
	if (h1->xrange[i] != h->xrange[i])
	  status = 1;
      };
    gsl_test (status, "gsl_histogram2d_memcpy copies bin xranges correctly");
  }

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
	if (h1->yrange[i] != h->yrange[i])
	  status = 1;
      };
    gsl_test (status, "gsl_histogram2d_memcpy copies bin yranges correctly");
  }

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            if (gsl_histogram2d_get (h1, i, j) != gsl_histogram2d_get (h, i, j))
              status = 1;
          }
      }
    gsl_test (status, "gsl_histogram2d_memcpy copies bin values correctly");
  }

  gsl_histogram2d_free (h1);

  h1 = gsl_histogram2d_clone (h);

  {
    int status = 0;
    for (i = 0; i <= M; i++)
      {
	if (h1->xrange[i] != h->xrange[i])
	  status = 1;
      };
    gsl_test (status, "gsl_histogram2d_clone copies bin xranges correctly");
  }

  {
    int status = 0;
    for (i = 0; i <= N; i++)
      {
	if (h1->yrange[i] != h->yrange[i])
	  status = 1;
      };
    gsl_test (status, "gsl_histogram2d_clone copies bin yranges correctly");
  }

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
        for (j = 0; j < N; j++)
          {
            if (gsl_histogram2d_get (h1, i, j) != gsl_histogram2d_get (h, i, j))
              status = 1;
          }
      }
    gsl_test (status, "gsl_histogram2d_clone copies bin values correctly");
  }

  gsl_histogram2d_reset (h);

  {
    int status = 0;

    for (i = 0; i < M * N; i++)
      {
	if (h->bin[i] != 0)
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_reset zeros array correctly");
  }

  gsl_histogram2d_free (h);
  h = gsl_histogram2d_calloc (M1, N1);

  {

    int status = 0;

    for (i = 0; i < M1; i++)
      {
	for (j = 0; j < N1; j++)
	  {
	    gsl_histogram2d_increment (h, (double) i, (double) j);

	    for (k = 0; k <= i * N1 + j; k++)
	      {
		if (h->bin[k] != 1)
		  {
		    status = 1;
		  }
	      }

	    for (k = i * N1 + j + 1; k < M1 * N1; k++)
	      {
		if (h->bin[k] != 0)
		  {
		    status = 1;
		  }
	      }
	  }
      }
    gsl_test (status, "gsl_histogram2d_increment works correctly");
  }

  gsl_histogram2d_free (h);
  h = gsl_histogram2d_calloc (M, N);

  {
    int status = 0;
    for (i = 0; i < M; i++)
      {
	double x0 = 0, x1 = 0;
	gsl_histogram2d_get_xrange (h, i, &x0, &x1);

	if (x0 != i || x1 != i + 1)
	  {
	    status = 1;
	  }
      }
    gsl_test (status,
      "gsl_histogram2d_get_xlowerlimit and xupperlimit works correctly");
  }

  {
    int status = 0;
    for (i = 0; i < N; i++)
      {
	double y0 = 0, y1 = 0;
	gsl_histogram2d_get_yrange (h, i, &y0, &y1);

	if (y0 != i || y1 != i + 1)
	  {
	    status = 1;
	  }
      }
    gsl_test (status, "gsl_histogram2d_get_ylowerlimit and yupperlimit works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_xmax (h) != M)
      status = 1;
    gsl_test (status, "gsl_histogram2d_xmax works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_xmin (h) != 0)
      status = 1;
    gsl_test (status, "gsl_histogram2d_xmin works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_nx (h) != M)
      status = 1;
    gsl_test (status, "gsl_histogram2d_nx works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ymax (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ymax works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ymin (h) != 0)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ymin works correctly");
  }

  {
    int status = 0;
    if (gsl_histogram2d_ny (h) != N)
      status = 1;
    gsl_test (status, "gsl_histogram2d_ny works correctly");
  }

  h->bin[3*N+2] = 123456.0;
  h->bin[4*N+3] = -654321;

  {
    double max = gsl_histogram2d_max_val (h);
    gsl_test (max != 123456.0, "gsl_histogram2d_max_val finds maximum value");
  }

  {
    double min = gsl_histogram2d_min_val (h);
    gsl_test (min != -654321.0, "gsl_histogram2d_min_val finds minimum value");
  }

  {
    size_t imax, jmax;
    gsl_histogram2d_max_bin (h, &imax, &jmax);
    gsl_test (imax != 3 || jmax != 2, "gsl_histogram2d_max_bin finds maximum value bin");
  }

  {
    size_t imin, jmin;
    gsl_histogram2d_min_bin (h, &imin, &jmin);
    gsl_test (imin != 4 || jmin != 3, "gsl_histogram2d_min_bin find minimum value bin");
  }

  for (i = 0; i < M*N; i++)
    {
      h->bin[i] = i + 27;
      g->bin[i] = (i + 27) * (i + 1);
    }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_add (h1, h);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != g->bin[i] + h->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_add works correctly");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_sub (h1, h);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != g->bin[i] - h->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_sub works correctly");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_mul (h1, h);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != g->bin[i] * h->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_mul works correctly");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_div (h1, h);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != g->bin[i] / h->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_div works correctly");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_scale (h1, 0.5);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != 0.5 * g->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_scale works correctly");
  }

  gsl_histogram2d_memcpy (h1, g);
  gsl_histogram2d_shift (h1, 0.25);

  {
    int status = 0;
    for (i = 0; i < M*N; i++)
      {
	if (h1->bin[i] != 0.25 + g->bin[i])
	  status = 1;
      }
    gsl_test (status, "gsl_histogram2d_shift works correctly");
  }

  gsl_histogram2d_free (h);	/* free whatever is in h */

  h = gsl_histogram2d_calloc_uniform (M1, N1, 0.0, 5.0, 0.0, 5.0);

  gsl_test (h->xrange == 0,
	    "gsl_histogram2d_calloc_uniform returns valid range pointer");
  gsl_test (h->yrange == 0,
	    "gsl_histogram2d_calloc_uniform returns valid range pointer");
  gsl_test (h->bin == 0,
	    "gsl_histogram2d_calloc_uniform returns valid bin pointer");
  gsl_test (h->nx != M1,
	    "gsl_histogram2d_calloc_uniform returns valid nx");
  gsl_test (h->ny != N1,
	    "gsl_histogram2d_calloc_uniform returns valid ny");

  gsl_histogram2d_accumulate (h, 0.0, 3.01, 1.0);
  gsl_histogram2d_accumulate (h, 0.1, 2.01, 2.0);
  gsl_histogram2d_accumulate (h, 0.2, 1.01, 3.0);
  gsl_histogram2d_accumulate (h, 0.3, 0.01, 4.0);

  {
    size_t i1, i2, i3, i4;
    size_t j1, j2, j3, j4;
    double expected;
    int status;
    status = gsl_histogram2d_find (h, 0.0, 3.01, &i1, &j1);
    status = gsl_histogram2d_find (h, 0.1, 2.01, &i2, &j2);
    status = gsl_histogram2d_find (h, 0.2, 1.01, &i3, &j3);
    status = gsl_histogram2d_find (h, 0.3, 0.01, &i4, &j4);

    for (i = 0; i < M1; i++)
      {
	for (j = 0; j < N1; j++)
	  {
	    if (i == i1 && j == j1)
	      {
		expected = 1.0;
	      }
	    else if (i == i2 && j == j2)
	      {
		expected = 2.0;
	      }
	    else if (i == i3 && j == j3)
	      {
		expected = 3.0;
	      }
	    else if (i == i4 && j == j4)
	      {
		expected = 4.0;
	      }
	    else
	      {
		expected = 0.0;
	      }

	    if (h->bin[i * N1 + j] != expected)
	      {
		status = 1;
	      }
	  }
      }
    gsl_test (status, "gsl_histogram2d_find works correctly");
  }

  {
    FILE *f = fopen ("test.txt", "w");
    gsl_histogram2d_fprintf (f, h, "%.19g", "%.19g");
    fclose (f);
  }

  {
    FILE *f = fopen ("test.txt", "r");
    gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
    int status = 0;

    gsl_histogram2d_fscanf (f, hh);

    for (i = 0; i <= M1; i++)
      {
	if (h->xrange[i] != hh->xrange[i])
	  {
	    printf ("xrange[%d] : %g orig vs %g\n",
		    (int)i, h->xrange[i], hh->xrange[i]);
	    status = 1;
	  }
      }

    for (j = 0; j <= N1; j++)
      {
	if (h->yrange[j] != hh->yrange[j])
	  {
	    printf ("yrange[%d] : %g orig vs %g\n",
		    (int)j, h->yrange[j], hh->yrange[j]);
	    status = 1;
	  }
      }

    for (i = 0; i < M1 * N1; i++)
      {
	if (h->bin[i] != hh->bin[i])
	  {
	    printf ("bin[%d] : %g orig vs %g\n",
		    (int)i, h->bin[i], hh->bin[i]);
	    status = 1;
	  }
      }

    gsl_test (status, "gsl_histogram2d_fprintf and fscanf work correctly");

    gsl_histogram2d_free (hh);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "wb");
    gsl_histogram2d_fwrite (f, h);
    fclose (f);
  }

  {
    FILE *f = fopen ("test.dat", "rb");
    gsl_histogram2d *hh = gsl_histogram2d_calloc (M1, N1);
    int status = 0;

    gsl_histogram2d_fread (f, hh);

    for (i = 0; i <= M1; i++)
      {
	if (h->xrange[i] != hh->xrange[i])
	  {
	    printf ("xrange[%d] : %g orig vs %g\n",
		    (int)i, h->xrange[i], hh->xrange[i]);
	    status = 1;
	  }
      }

    for (j = 0; j <= N1; j++)
      {
	if (h->yrange[j] != hh->yrange[j])
	  {
	    printf ("yrange[%d] : %g orig vs %g\n",
		    (int)j, h->yrange[j], hh->yrange[j]);
	    status = 1;
	  }
      }

    for (i = 0; i < M1 * N1; i++)
      {
	if (h->bin[i] != hh->bin[i])
	  {
	    printf ("bin[%d] : %g orig vs %g\n",
		    (int)i, h->bin[i], hh->bin[i]);
	    status = 1;
	  }
      }

    gsl_test (status, "gsl_histogram2d_fwrite and fread work correctly");

    gsl_histogram2d_free (hh);
    fclose (f);
  }

  gsl_histogram2d_free (h);
  gsl_histogram2d_free (h1);
  gsl_histogram2d_free (g);
  gsl_histogram2d_free (hr);

  exit (gsl_test_summary ());
}


text:
// Implementation of interface described in float_blob.h.

#include <errno.h>
#include <fcntl.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>

#include <glib.h>
#include <glib/gstdio.h>

#include <gsl/gsl_math.h>

#include "float_blob.h"
#include "float_image.h"
#include "utilities.h"

FloatBlob *
float_blob_new (gssize size_x, gssize size_y, const char *scratch_dir,
		const char *name)
{
  g_assert (size_x > 0 && size_y > 0);
  g_assert (g_file_test (scratch_dir, G_FILE_TEST_IS_DIR));
  g_assert (g_access (scratch_dir, R_OK | W_OK | X_OK) == 0);

  FloatBlob *self = g_new (FloatBlob, 1);

  self->size_x = size_x;
  self->size_y = size_y;

  self->swapped_bytes = FALSE;

  if ( name != NULL ) {
    self->file = ensure_slash_terminated (g_string_new (scratch_dir));
    g_string_append (self->file, name);
  }
  else {
    self->file = make_unique_tmp_file_name (scratch_dir, "float_blob_tmp_");
    g_assert (!g_file_test (self->file->str, G_FILE_TEST_EXISTS));
  }
  self->is_file_owner = TRUE;

  self->fd = open (self->file->str, O_RDWR | O_CREAT | O_TRUNC,
		   S_IRUSR | S_IWUSR);
  g_assert (self->fd != -1);

  // Out of carefulness we didn't truncate, but we want to make sure
  // we really got a new file.
  g_assert (file_size (self->fd) == 0);

  // Seek to last byte of file.
  off_t last_byte_offset
    = (off_t) self->size_x * (off_t) self->size_y * sizeof (float) - 1;
  off_t file_offset = lseek (self->fd, last_byte_offset, SEEK_SET);
  g_assert (file_offset == last_byte_offset);

  // Write a byte at the end of the file.
  size_t bytes_to_write = 1;
  g_assert (sizeof (unsigned char) == 1);
  unsigned char good_old_zero = 0x0;
  ssize_t bytes_written = write (self->fd, &good_old_zero, bytes_to_write);
  g_assert (bytes_written == bytes_to_write);

  // We should now have reserved all the space needed in this file.
  g_assert (file_size (self->fd) == last_byte_offset + 1);

  self->reference_count = 1;

  return self;
}

FloatBlob *
float_blob_new_using (gssize size_x, gssize size_y, const char *file,
		      gboolean swapped_bytes)
{
  g_assert (size_x > 0 && size_y > 0);
  g_assert (g_file_test (file, G_FILE_TEST_IS_REGULAR));
  g_assert (g_access (file, R_OK | W_OK) == 0);

  FloatBlob *self = g_new (FloatBlob, 1);

  self->fd = open (file, O_RDWR);
  g_assert (self->fd != -1);

  g_assert (file_size (self->fd) >= (off_t) size_x * size_y * sizeof (float));

  self->size_x = size_x;
  self->size_y = size_y;

  self->swapped_bytes = swapped_bytes;

  self->file = g_string_new (file);

  self->is_file_owner = FALSE;

  self->reference_count = 1;

  return self;
}

static void
float_blob_region_transfer (FloatBlob *self, gssize start_x, gssize start_y,
			    gssize w, gssize h, float *buffer,
			    gboolean is_get_transfer)
{
  g_assert (start_x >= 0 && start_x < self->size_x);
  g_assert (start_y >= 0 && start_y < self->size_y);
  g_assert (w > 0 && start_x + w <= self->size_x);
  g_assert (h > 0 && start_y + h <= self->size_y);

  ssize_t (*transfer)(int fd, void *buf, size_t count)
    = is_get_transfer ? read : ((ssize_t (*)(int, void *, size_t)) write);

  size_t jj;
  for ( jj = start_y ; jj < start_y + h ; jj++ ) {
    off_t file_offset
      = lseek (self->fd,
	       sizeof (float) * ((off_t) jj * self->size_x + start_x),
	       SEEK_SET);
    g_assert (file_offset != (off_t) -1);
    size_t transfer_size = w * sizeof (float);
    ssize_t transfer_result
      = transfer (self->fd, buffer + (jj - start_y) * w, transfer_size);
    if ( transfer_result == -1 ) {
      g_error ("I/O error: %s\n", strerror (errno));
    }
    g_assert (transfer_result == transfer_size);
  }
}

void
float_blob_get_region (FloatBlob *self, gssize start_x, gssize start_y,
		       gssize w, gssize h, float *buffer)
{
  float_blob_region_transfer (self, start_x, start_y, w, h, buffer, TRUE);

  if ( self->swapped_bytes ) {
    swap_array_bytes_32 (buffer, w * h);
  }
}

void
float_blob_set_region (FloatBlob *self, gssize start_x, gssize start_y,
		       gssize w, gssize h, float *buffer)
{
  if ( self->swapped_bytes ) {
    swap_array_bytes_32 (buffer, w * h);
  }

  float_blob_region_transfer (self, start_x, start_y, w, h, buffer, FALSE);

  if ( self->swapped_bytes ) {
    swap_array_bytes_32 (buffer, w * h);
  }
}

GString *
float_blob_steal_file (FloatBlob *self)
{
  g_assert (self->is_file_owner == TRUE);

  self->is_file_owner = FALSE;

  return g_string_new (self->file->str);

}

void
float_blob_export_as_jpeg (FloatBlob *self, const char *file, double mask)
{
  // Self as float image.
  FloatImage *safi = float_image_new (self->size_x, self->size_y);

  float *row_buffer = g_new (float, self->size_x);

  // We are about to use swap_32 type function, so we have a quick
  // attack of paranoia and check sizeof float...
  g_assert (sizeof (float) == 4);

  // Seek to the beginning of the file.
  off_t resultant_offset = lseek (self->fd, 0, SEEK_SET);
  if ( resultant_offset == (off_t) -1 ) {
    g_error ("lseek failed: %s\n", strerror (errno));
  }
  g_assert (resultant_offset == 0);

  size_t ii, jj;
  for ( jj = 0 ; jj < self->size_y ; jj++ ) {

    // Read one row of data from the file underlying self.
    ssize_t bytes_to_read = self->size_x * sizeof (float);
    ssize_t read_count = read (self->fd, row_buffer, bytes_to_read);
    g_assert (read_count == bytes_to_read);

    // Swap row bytes if necessary.
    if ( self->swapped_bytes ) {
      swap_array_bytes_32 (row_buffer, self->size_x);
    }

    // Write that row to the FloatImage.
    for ( ii = 0 ; ii < self->size_x ; ii++ ) {
      float_image_set_pixel (safi, ii, jj, row_buffer[ii]);
    }
  }

  g_free (row_buffer);

  float_image_export_as_jpeg (safi, file, GSL_MAX (safi->size_x, safi->size_y),
			      mask);

  float_image_unref (safi);
}

FloatBlob *
float_blob_ref (FloatBlob *self)
{
  g_assert (self->reference_count > 0);

  self->reference_count++;

  return self;
}

void
float_blob_unref (FloatBlob *self)
{
  g_assert (self->reference_count > 0);

  self->reference_count--;

  if ( self->reference_count == 0 ) {
    int return_code = close (self->fd);
    g_assert (return_code == 0);

    if ( self->is_file_owner ) {
      return_code = unlink (self->file->str);
      if ( return_code != 0 ) {
	g_error ("unlink of float_blob file '%s' in function %s, (file "
		 __FILE__ ", line %d) failed: %s", self->file->str, __func__,
		 __LINE__, strerror (errno));
      }
      g_assert (return_code == 0);
    }

    g_string_free (self->file, TRUE);

    g_free (self);

    self = NULL;
  }
}

max_stars_repo_name: glshort/MapReady

max_issues_repo_name: glshort/MapReady

max_forks_repo_name: glshort/MapReady


text:
/*
 * Copyright 2021 The DAPHNE Consortium
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef SRC_RUNTIME_LOCAL_KERNELS_SYRK_H
#define SRC_RUNTIME_LOCAL_KERNELS_SYRK_H

#include <runtime/local/context/DaphneContext.h>
#include <runtime/local/datastructures/CSRMatrix.h>
#include <runtime/local/datastructures/DataObjectFactory.h>
#include <runtime/local/datastructures/DenseMatrix.h>

#include <cblas.h>

// ****************************************************************************
// Struct for partial template specialization
// ****************************************************************************

template<class DTRes, class DTArg>
struct Syrk {
    static void apply(DTRes *& res, const DTArg * arg, DCTX(ctx)) = delete;
};

// ****************************************************************************
// Convenience function
// ****************************************************************************

template<class DTRes, class DTArg>
void syrk(DTRes *& res, const DTArg * arg, DCTX(ctx)) {
    Syrk<DTRes, DTArg>::apply(res, arg, ctx);
}

// ****************************************************************************
// (Partial) template specializations for different data/value types
// ****************************************************************************

// ----------------------------------------------------------------------------
// DenseMatrix <- DenseMatrix
// ----------------------------------------------------------------------------

template<>
struct Syrk<DenseMatrix<double>, DenseMatrix<double>> {
    static void apply(DenseMatrix<double> *& res, const DenseMatrix<double> * arg, DCTX(ctx)) {
        const size_t numRows = arg->getNumRows();
        const size_t numCols = arg->getNumCols();

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<double>>(numCols, numCols, false);

        cblas_dsyrk(CblasRowMajor,
            CblasUpper,
            CblasTrans,
            numCols,
            numRows,
            1.0,
            arg->getValues(),
            arg->getRowSkip(),
            0.0,
            res->getValues(),
            res->getRowSkip());
        for (auto r = 0u; r < numCols; ++r) {
            for (auto c = r + 1; c < numCols; ++c) {
                res->set(c, r, res->get(r, c));
            }
        }
    }
};

template<>
struct Syrk<DenseMatrix<float>, DenseMatrix<float>> {
    static void apply(DenseMatrix<float> *& res, const DenseMatrix<float> * arg, DCTX(ctx)) {
        const size_t numRows = arg->getNumRows();
        const size_t numCols = arg->getNumCols();

        if(res == nullptr)
            res = DataObjectFactory::create<DenseMatrix<float>>(numCols, numCols, false);

        cblas_ssyrk(CblasRowMajor,
            CblasUpper,
            CblasTrans,
            numCols,
            numRows,
            1.0,
            arg->getValues(),
            arg->getRowSkip(),
            0.0,
            res->getValues(),
            res->getRowSkip());
        for (auto r = 0u; r < numCols; ++r) {
            for (auto c = r + 1; c < numCols; ++c) {
                res->set(c, r, res->get(r, c));
            }
        }
    }
};

// ----------------------------------------------------------------------------
// CSRMatrix <- CSRMatrix
// ----------------------------------------------------------------------------

template<typename VT>
struct Syrk<CSRMatrix<VT>, CSRMatrix<VT>> {
    static void apply(CSRMatrix<VT> *& res, const CSRMatrix<VT> * arg, DCTX(ctx)) {
        const size_t numRows = arg->getNumRows();
        const size_t numCols = arg->getNumCols();

        if(res == nullptr)
            res = DataObjectFactory::create<CSRMatrix<VT>>(numCols, numRows, arg->getNumNonZeros(), false);
        assert(false && "TODO: Syrk for Sparse");
    }
};

#endif //SRC_RUNTIME_LOCAL_KERNELS_SYRK_H

max_stars_repo_name: daphne-eu/daphne

max_issues_repo_name: daphne-eu/daphne

max_forks_repo_name: daphne-eu/daphne


text:
#include <iostream>
#include <vector>
#include <memory>
#include <cstdlib>
#include <cstddef>
#ifdef ENABLE_GSL
#include <gsl/gsl_integration.h>
#include <gsl/gsl_deriv.h>
#include <cmath>

#endif

#ifndef PI
#define PI  3.141593
#endif

#ifndef Grav
#define Grav 4.7788e-20  // Newton's over c^2 in Mpc / Msun
#endif

#ifndef lightspeed
#define lightspeed 2.99792458e5  // km / s
#endif

#ifndef error_message
#define error_message
#define ERROR_MESSAGE() std::cerr << "ERROR: file: " << __FILE__ << " line: " << __LINE__ << std::endl;
#endif

#define CRITD 2.49783e18    /* critical density / Ho^2 solar/Mpc */
#define CRITD2 2.7752543e11 /* critical density / h^2 M_sun/Mpc^3 */

#ifndef cosmo_declare

enum CosmoParamSet {WMAP5yr,Millennium,Planck1yr,Planck15,Planck18,BigMultiDark,none};

std::ostream &operator<<(std::ostream &os, CosmoParamSet p);

/**
 *
 * \brief The cosmology and all the functions required to calculated quantities based on the cosmology.
 *
 * This class is used to store the cosmological parameters, calculate cosmological distances, calculate
 * the power spectrum of density fluctuations and the mass function of halos.
 *
 * As set now, there are no baryon acoustic oscillations in the power spectrum, but this can be changed
 * at the expense of not including neutrinos.
 */
class COSMOLOGY{
public:

  COSMOLOGY(CosmoParamSet cosmo_p = WMAP5yr);
  /// if justdistances== true
  COSMOLOGY(double omegam,double omegal,double h, double w = -1,double wa = 0
            ,bool justdistances=false /// if true the internals needed calculate structure formation will not be calculated making constructuion faster
  );
  COSMOLOGY(const COSMOLOGY &cosmo);
  ~COSMOLOGY();

  COSMOLOGY& operator=(const COSMOLOGY &cosmo);

  // returns the parameter set if any
  CosmoParamSet ParamSet(){return cosmo_set;}

  void PrintCosmology(short physical = 0) const;
  std::string print(short physical= 0) const;

  // Lengths
  double rcurve() const;
  //double emptyDist(double zo,double z);
  double coorDist(double zo,double z) const;
  double coorDist(double z) const {return coorDist(0,z);}
  double radDist(double zo,double z) const;
  double radDist(double z) const {return radDist(0,z);}
  double angDist(double zo,double z) const;
  double angDist(double z) const {return angDist(0,z);}
  double lumDist(double zo,double z) const;
  double lumDist(double z) const {return lumDist(0,z);}

  /// derivative of coordinate distance with respect to Omo in flat cosmology.  Useful for Fisher matrix calculations
  double d_coorDist_dOmo(double zo,double z) const;
  /// derivative of coordinate distance with respect to w.  Useful for Fisher matrix calculations
  double d_coorDist_dw(double zo,double z) const;
  /// derivative of coordinate distance with respect to w1.  Useful for Fisher matrix calculations
  double d_coorDist_dw1(double zo,double z) const;

  double DRradius(double zo,double z,double pfrac);
  double DRradius2(double zo,double z);

  double invCoorDist(double d) const;
  double invRadDist(double d) const;
  double invComovingDist(double d) const;

  double scalefactor(double rad) const;
  double Omegam(double z) const;
  double rho_crit_comoving(double z) const;
  //double drdz_empty(double x);
  double drdz(double x) const;
  //double adrdz(double x) const;
  double drdz_dark(double x) const;
  inline double dark_factor(double x) const;
  //double adrdz_dark(double x) const;

  /// the critical over density
  double delta_c() const;

  double DeltaVir(double z,int caseunit=0) const;
  double Deltao(double m) const;
  double time(double z) const;
  double nonlinMass(double z) const;

  // Stuff having to do with the power spectrum
  double power_normalize(double sigma8);
  /**
   * \brief Linear power spectrum P(k,z)/a^2
   */
  double power_linear(double k,double z);
  double Dgrowth(double z) const;

  /**
   * \brief  powerCDM.c calculates the nonlinear P(k,z)/a(r)^2
   *
   * The method of Peacock & Dodds 1996 is used to convert the linear
   * power spectrum into the nonlinear one.
   * This could be updated to a more recent nonlinear power spectrum
   */
  double powerCDMz(double k,double z);
  double psdfdm(double z,double m,int caseunit=0);
  double halo_bias (double m, double z, int t=0);
  double stdfdm(double z,double m,int caseunit=0);
  double powerlawdfdm(double z,double m,double alpha,int caseunit=0);
  double haloNumberDensity(double m,double z,double a, int t,double alpha = 0.0);

  /** \brief Dark matter correlation function
   *
   *  This integrates powerCDMz() to get the correlation function as a function of comoving radius.
   *  Care should be taken that the right range of k is integrated over if the radius is very small or large.
   */
  double CorrelationFunction(double radius,double redshift
                      ,double k_max = 100,double k_min = 1.0e-3);

  struct CorrFunctorType{
    CorrFunctorType(COSMOLOGY *cosmo,double radius,double redshift)
    : cosmology(cosmo),z(redshift),r(radius)
    {
      norm = 0.5/PI/PI/(1+z)/(1+z);
    };

    COSMOLOGY *cosmology;
    double z;
    double r;

    double norm;

    double operator () (double k) {

      double rk = r*k;
      double jo = sin(rk)/rk;
      if(rk < 1.0e-3) jo = 1.0;

      return norm*jo*k*k*cosmology->powerCDMz(k,z);
    }
  };

  double haloNumberDensityOnSky (double m,double z1,double z2,int t,double alpha = 0.0);
  double haloNumberInBufferedCone (double mass ,double z1,double z2,double fov,double buffer ,int type ,double alpha=0.0);
  double haloMassInBufferedCone (double mass ,double z1,double z2,double fov,double buffer ,int type ,double alpha=0.0);

  double TopHatVariance(double m) const;
  double TopHatVarianceR(double R,double z);
  double TopHatVarianceM(double M,double z);
  //double gradius(double R,double rd);

  double getTimefromDeltaC(double dc);
  double getZfromDeltaC(double dc);

  /// accesser functions:

  /// Hubble paremters in units of 100 km/s/Mpc, renormalizes P(k) to keep sig8 fixed
  void sethubble(double ht){ h = ht; TFmdm_set_cosm();
    power_normalize(sig8); 	calc_interp_dist(); cosmo_set = none;}
  double gethubble() const {return h;}
  /// Hubble parameter in 1/Mpc units
  double getHubble() const {return 100*h/lightspeed;}

  /// Primordial spectral index, renormalizes P(k) to keep sig8 fixed
  void setindex(double nn){ n = nn;  power_normalize(sig8); cosmo_set = none;}
  double getindex() const { return n;}

  /// Omega matter, renormalizes P(k) to keep sig8 fixed
  void setOmega_matter(double Omega_matter,bool FLAT = false){Omo = Omega_matter; if(FLAT) Oml = 1-Omo ; TFmdm_set_cosm(); power_normalize(sig8); calc_interp_dist(); cosmo_set = none;}
  double getOmega_matter() const {return Omo;}

  /// Omega lambda, renormalizes P(k) to keep sig8 fixed
  void setOmega_lambda(double Omega_lambda,bool FLAT = false){Oml = Omega_lambda;  if(FLAT) Oml = 1-Omo ; TFmdm_set_cosm(); power_normalize(sig8); calc_interp_dist(); cosmo_set = none;}
  double getOmega_lambda() const {return Oml;}

  /// Omega baryon, renormalizes P(k) to keep sig8 fixed
  void setOmega_baryon(double Omega_baryon){Omb = Omega_baryon; TFmdm_set_cosm(); power_normalize(sig8); calc_interp_dist(); cosmo_set = none;
  }
  double getOmega_baryon() const {return Omb;}

  /// Omega neutrino, renormalizes P(k) to keep sig8 fixed
  void setOmega_neutrino(double Omega_neutrino){Omnu = Omega_neutrino; TFmdm_set_cosm(); power_normalize(sig8); 	calc_interp_dist(); cosmo_set = none;
  }
  double getOmega_neutrino() const {return Omnu;}

  /// Number of neutrino species, renormalizes P(k) to keep sig8 fixed
  void setNneutrino(double Nneutrino){Nnu = Nneutrino; TFmdm_set_cosm(); power_normalize(sig8); cosmo_set = none;}
  double getNneutrino() const {return Nnu;}

  /// Dark energy equation of state parameter p/rho = w + w_1 (1+z)
  void setW(double w){ww = w; 	calc_interp_dist(); cosmo_set = none;}
  double getW() const {return ww;}
  void setW1(double w){ww1 = w; 	calc_interp_dist(); cosmo_set = none;}
  double getW1() const {return ww1;}

  /// Running of primordial spectral index, P(k)_primordial \propto pow(k/h,n+dndlnk*log(k)), renormalizes P(k) to keep sig8 fixed
  void setdndlnk(double w){dndlnk = w; power_normalize(sig8); cosmo_set = none;}
  double getdndlnk() const {return dndlnk;}

  /// Alternative to w for dark energy/ alt. grav. structure evolution
  void setgamma(double gamm){gamma = gamm; cosmo_set = none;}
  double getgamma() const {return gamma;}
  // If physical = 1 all Omega are Omega*h^2 otherwise they have the usual definitions.
  // 2 gamma parameterization is used for dark energy
  void setDEtype(short tt){darkenergy = tt; cosmo_set = none;}
  short getDEtype() const {return darkenergy;}

  void setSigma8(double my_sig8){power_normalize(my_sig8); cosmo_set = none;}
  double getSigma8() const {return sig8;}

  void dzdangDist(double D,double z[],double dzdD[]);

  double totalMassDensityinHalos(int t,double alpha,double m_min,double z,double z1,double z2);

  /// set interpolation range and number of points
  void setInterpolation(double z_interp, std::size_t n_interp);

  /// The lensing critical density in Msun / Mpc^2
  double SigmaCrit(double zlens,double zsource) const;

protected:
  void SetConcordenceCosmology(CosmoParamSet cosmo_p);

  CosmoParamSet cosmo_set;

  // structure rappers to make integration thread safe
  struct drdz_struct{
    drdz_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.drdz(x);}
    COSMOLOGY const &cos;
  };
  struct drdzdark_struct{
    drdzdark_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.drdz_dark(x);}
    COSMOLOGY const &cos;
  };
  struct adrdz_struct{
    adrdz_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.drdz(x)/x;}
    COSMOLOGY const &cos;
  };
  struct adrdzdark_struct{
    adrdzdark_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.drdz_dark(x)/x;}
    COSMOLOGY const &cos;
  };

  struct normL_struct{
    normL_struct(COSMOLOGY &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.normL(x);}
    COSMOLOGY &cos;
  };

  double ddrdzdOmo(double x) const;
  double ddrdzdw(double x) const;
  double ddrdzdw1(double x) const;

  struct ddrdzdOmo_struct{
    ddrdzdOmo_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.ddrdzdOmo(x);}
    COSMOLOGY const &cos;
  };
  struct ddrdzdw_struct{
    ddrdzdw_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.ddrdzdw(x);}
    COSMOLOGY const &cos;
  };
  struct ddrdzdw1_struct{
    ddrdzdw1_struct(COSMOLOGY const &cosmo):cos(cosmo){};
    double operator()(double x) const {return cos.ddrdzdw1(x);}
    COSMOLOGY const &cos;
  };

  bool init_structure_functions;

  /// Hubble paremters in units of 100 km/s/Mpc
  double h;
  /// Primordial spectral index
  double n;
  /// Omega matter, dark matter + baryons
  double Omo;
  /// Omega lambda
  double Oml;
  /// Omega baryon
  double Omb;
  /// Omega neutrino
  double Omnu;
  /// Number of neutrino species
  double Nnu;
  /// Dark energy equation of state parameter p/rho = ww + ww_1 (1+z)
  double ww;
  /// Dark energy equation of state parameter p/rho = ww + ww_1 (1+z)
  double ww1;
  /// Running of primordial spectral index, P(k)_primordial \propto pow(k/h,n+dndlnk*log(k))
  double dndlnk;

  /// Alternative to w for dark energy/ alt. grav. structure evolution
  double gamma;
  // If physical = 1 all Omega are Omega*h^2 otherwise they have the usual definitions.
  //short physical;
  // 2 gamma parameterization is used for dark energy
  short darkenergy;

  /* table for growth parameter */
  //std::auto_ptr<double> *aa;
  //std::auto_ptr<double> *growth;
  //int Ntable;

  double A;
  double sig8;  /* do not access these normalization outside */

  double Rtophat;
  double ztmp;

  double powerEH(double k,double z);
  double powerEHv2(double k);
  double powerloc(double k,double z);
  double npow(double k);
  double normL(double lgk);
  double gradius(double R,double rd) const;

  int nbin;  // number of points when binning to interpolate
  std:: vector<double> vDeltaCz,vlz,vt;
  double DpropDz(double z);
  double dsigdM(double m);
  double timeEarly(double a) const;
  double dNdz(double z);

  typedef double (COSMOLOGY::*pt2MemFunc)(double) const;
  typedef double (COSMOLOGY::*pt2MemFuncNonConst)(double);

  double nintegrateDcos(pt2MemFunc func,double a,double b,double tols) const;
  double trapzdDcoslocal(pt2MemFunc func, double a, double b, int n, double *s2) const;
  void polintD(double xa[], double ya[], int n, double x, double *y, double *dy) const;
  double nintegrateDcos(pt2MemFuncNonConst func,double a,double b,double tols);
  double trapzdDcoslocal(pt2MemFuncNonConst func, double a, double b, int n, double *s2);
  double dfridrDcos(pt2MemFunc func, double x, double h, double *err);
  double f4(double u) const;

  int ni;
  std::vector<float> xf, wf;

  // The parameters used in Eisenstein & Hu power spectrum
  /* in powerEH.c */
  short TFmdm_set_cosm_change_z(double redshift);
  short TFmdm_set_cosm();
  double TFmdm_onek_mpc(double kk);
  double TFmdm_onek_hmpc(double kk);

  double f_baryon;	// Baryon fraction
  double f_bnu;		// Baryon + Massive Neutrino fraction
  double f_cb;		// Baryon + CDM fraction
  double f_cdm;		// CDM fraction
  double f_hdm;		// Massive Neutrino fraction

  double alpha_gamma;	 // sqrt(alpha_nu)
  double alpha_nu;	     // The small-scale suppression
  double beta_c;		 // The correction to the log in the small-scale
  double num_degen_hdm;	 // Number of degenerate massive neutrino species
  double growth_k0;	     // D_1(z) -- the growth function as k->0
  double growth_to_z0;	 // D_1(z)/D_1(0) -- the growth relative to z=0
  double k_equality;	 // The comoving wave number of the horizon at equality
  double obhh;		     // Omega_baryon * hubble^2
  double omega_curv;	 // = 1 - omega_matter - omega_lambda
  double omega_lambda_z; // Omega_lambda at the given redshift
  double omega_matter_z; // Omega_matter at the given redshift
  double omhh;		     // Omega_matter * hubble^2
  double onhh;		     // Omega_hdm * hubble^2
  double p_c;		     // The correction to the exponent before drag epoch
  double p_cb;		     // The correction to the exponent after drag epoch
  double sound_horizon_fit;  // The sound horizon at the drag epoch
  double theta_cmb;	     // The temperature of the CMB, in units of 2.7 K
  double y_drag;		 // Ratio of z_equality to z_drag
  double z_drag;		 // Redshift of the drag epoch
  double z_equality;	 // Redshift of matter-radiation equality

  void setinternals();
  // in powerEHv2.c
  void TFset_parameters(double omega0hh, double f_baryon, double Tcmb);
  double TFfit_onek(double k, double *tf_baryon, double *tf_cdm);

  double R_drag;		// Photon-baryon ratio at drag epoch
  double R_equality;	// Photon-baryon ratio at equality epoch
  double sound_horizon;	// Sound horizon at drag epoch, in Mpc
  double k_silk;		// Silk damping scale, in Mpc^-1
  double alpha_c;	    // CDM suppression
  double alpha_b;	    // Baryon suppression
  double beta_b;		// Baryon envelope shift
  double beta_node;	// Sound horizon shift
  double k_peak;		// Fit to wavenumber of first peak, in Mpc^-1

  // temporary variables for doing interations
  int tmp_type;
  double tmp_alpha;
  double tmp_mass;
  double tmp_a;

  // interpolation of functions
  double z_interp;
  std::size_t n_interp;
  void calc_interp_dist();
  double interp(const std::vector<double>& table, double z) const;
  double invert(const std::vector<double>& table, double f_z) const;
  std::vector<double> redshift_interp;
  std::vector<double> coorDist_interp;
  std::vector<double> radDist_interp;
};

typedef COSMOLOGY *CosmoHndl;
/**
 *  \brief Class for calculating properties of NFW halo profile.
 *
 *   This class does not take into affect the cosmological correlation between concentration and mass.
 *   For this see the HALOCalculator class.
 */
class NFW_Utility {
public:
  NFW_Utility(){return;}
  ~NFW_Utility(){};

  // methods for NFW profile
  double NFW_V200(double M200,double R200);
  double NFW_Vmax(double cons,double M200,double R200);
  double NFW_Vr(double x,double cons,double M200,double R200);
  double NFW_M(double x,double cons,double M200);
  double NFW_deltac(double cons);
  double NFW_Concentration(double Vmax,double M200,double R200);
  double NFW_rho(double cons,double x);

  void match_nfw(float Vmax,float R_half,float mass,float *cons,float *Rsize);
  float Rsize(float cons,float Vmax,float mass);
  float g_func(float x);

private:
  float Vmax,R_half,mass;       /// Mass (solar masses)

  typedef float (NFW_Utility::*MemFunc)(float);
  float zbrentD(MemFunc func, float a,float b,float tols);
  float nfwfunc(float cons);
  float funcforconcentration(float cons);
};

/// wrapper functions for gsl integration / differentiation
double drdz_wrapper(double x, void *params);
double drdz_dark_wrapper(double x, void *params);
double Deltao_wrapper(double m, void *params);

#define cosmo_declare
#endif

/*** in cosmo.c ***/
int cosmo_compare(CosmoHndl cos1,CosmoHndl cos2);
void cosmo_copy(CosmoHndl cos1,CosmoHndl cos2);
void ders(double z,double Da[],double dDdz[]);
void dir(double r,double a[],double dadr[]);
double arctanh(double x);
double fmini(double a,double b);
double fmaxi(double a,double b);
double **dmatrixcos(long nrl, long nrh, long ncl, long nch);
void free_dmatrixcos(double **m, long nrl, long nrh, long ncl, long nch);

/* in nfw.c */

/* in powerEHv2.c */
void TFset_parameters(double omega0hh, double f_baryon, double Tcmb);
double TFfit_onek(double k, double *tf_baryon, double *tf_cdm);

max_stars_repo_name: glenco/CosmoLib

max_issues_repo_name: glenco/CosmoLib

max_forks_repo_name: glenco/CosmoLib


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

/*!
 * @brief Computes a synthetic seismogram corresponding for the
 *        given observations, location (depth), and moment tensor.
 *
 * @param[in] iobs      Observations number.
 * @param[in] ilocOpt   Location index.
 * @param[in] imtOpt    Moment tensor number.
 * @param[in] ldm       Leading dimension of mts.  This must be at leats 6.
 * @param[in] mts       Matrix of moement tensors with dimension [nmt x ldm].
 *
 * @param[in,out] data  On input contains the data.
 * @param[in,out] data  On exit contains the corresponding synthetic.
 *
 * @result 0 indicates success.
 *
 */
int parmt_computeSynthetic(const int iobs,
                           const int ilocOpt,
                           const int imtOpt,
                           const int ldm,
                           const double *__restrict__ mts,
                           struct parmtData_struct *data)
{
    double *G;
    int k, npts;
    k = iobs*data->nlocs + ilocOpt;
    fprintf(stdout, "%s: Processing observation=%d, location=%d, mt=%d, k=%d\n",
            __func__, iobs, ilocOpt, imtOpt, k);
    npts = data->sacGxx[k].header.npts;
    G = memory_calloc64f(6*npts);
    sacio_copyHeader(data->sacGxx[k].header, &data->est[iobs].header);
    cblas_dcopy(data->sacGxx[k].header.npts,
                data->sacGxx[k].data, 1, &G[0*npts], 1);
    cblas_dcopy(data->sacGyy[k].header.npts,
                data->sacGyy[k].data, 1, &G[1*npts], 1);
    cblas_dcopy(data->sacGzz[k].header.npts,
                data->sacGzz[k].data, 1, &G[2*npts], 1);
    cblas_dcopy(data->sacGxy[k].header.npts,
                data->sacGxy[k].data, 1, &G[3*npts], 1);
    cblas_dcopy(data->sacGxz[k].header.npts,
                data->sacGxz[k].data, 1, &G[4*npts], 1);
    cblas_dcopy(data->sacGyz[k].header.npts,
                data->sacGyz[k].data, 1, &G[5*npts], 1);
    data->est[iobs].data = memory_calloc64f(npts);
    data->est[iobs].npts = npts;
    strcpy(data->est[iobs].header.kcmpnm, data->data[iobs].header.kcmpnm);
    cblas_dgemv(CblasColMajor, CblasNoTrans,
                npts, 6, 1.0, G, npts,
                &mts[ldm*imtOpt], 1, 0.0, data->est[iobs].data, 1);
    memory_free64f(&G);
    return 0;
}


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_statistics.h>

/*
gsl_movstat_Sn()
  Calculate moving S_n statistic for input vector

Inputs: endtype - how to handle end points
        x       - input vector, size n
        xscale  - (output) vector of "S_n" statistics, size n
                  xscale_i = (S_n)_i for i-th window:
        w       - workspace
*/

int
gsl_movstat_Sn(const gsl_movstat_end_t endtype, const gsl_vector * x,
               gsl_vector * xscale, gsl_movstat_workspace * w)
{
  int status = gsl_movstat_apply_accum(endtype, x, gsl_movstat_accum_Sn, NULL, xscale, NULL, w);
  return status;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#pragma once

#define QT_IS_BROKEN
#define QAPPLICATION_CLASS QApplication // for SingleApplication

// Windows
#include <Windows.h>
#include <SetupAPI.h>

#include <winioctl.h>
#include <BluetoothAPIs.h>
#include <Bthioctl.h>

#include <Dbt.h>
#include <Hidsdi.h>
#include <shellapi.h>
#include <Xinput.h>

// STL
#include <array>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <deque>
#include <functional>
#include <iomanip>
#include <list>
#include <map>
#include <memory>
#include <mutex>
#include <optional>
#include <sstream>
#include <string>
#include <thread>
#include <unordered_map>

// Qt
#include <QDialog>
#include <QtWidgets/QApplication>
#include <QtWidgets/QMainWindow>
#include <QtWidgets>

#include <singleapplication.h>

// GSL
#include <gsl/span>

// better-enums
#include <enum.h>

// fmt
#include <fmt/format.h>

// libhid
#include <hid_handle.h>
#include <hid_instance.h>
#include <hid_util.h>

// ViGEm
#include <ViGEm/Client.h>
#include <ViGEm/Common.h>
#include <ViGEm/km/BusShared.h>

#include "MapCache.h"
#include "average.h"
#include "AxisOptions.h"
#include "Bluetooth.h"
#include "busenum.h"
#include "ConnectionType.h"
#include "DeviceIdleOptions.h"
#include "DeviceProfile.h"
#include "DeviceProfileCache.h"
#include "DevicePropertiesDialog.h"
#include "DeviceSettings.h"
#include "DeviceSettingsCommon.h"
#include "Ds4AutoLightColor.h"
#include "Ds4Color.h"
#include "Ds4Device.h"
#include "Ds4DeviceManager.h"
#include "Ds4Input.h"
#include "Ds4InputData.h"
#include "Ds4ItemModel.h"
#include "Ds4LightOptions.h"
#include "Ds4Output.h"
#include "Ds4TouchRegion.h"
#include "enums.h"
#include "Event.h"
#include "gmath.h"
#include "InputMap.h"
#include "InputSimulator.h"
#include "JsonData.h"
#include "KeyboardSimulator.h"
#include "Latency.h"
#include "lock.h"
#include "Logger.h"
#include "MainWindow.h"
#include "MouseSimulator.h"
#include "pathutil.h"
#include "Pressable.h"
#include "ProfileEditorDialog.h"
#include "program.h"
#include "Settings.h"
#include "Stopwatch.h"
#include "stringutil.h"
#include "Trackball.h"
#include "Vector2.h"
#include "Vector3.h"
#include "XInputGamepad.h"


text:
/**
 *
 * @file
 *
 *  PLASMA is a software package provided by:
 *  University of Tennessee, US,
 *  University of Manchester, UK.
 *
 **/
#ifndef ICL_CORE_LAPACK_H
#define ICL_CORE_LAPACK_H

#ifdef PLASMA_WITH_MKL
    #include <mkl_cblas.h>
    #include <mkl_lapacke.h>

    // MKL LAPACKE doesn't provide LAPACK_GLOBAL macro, so define it here.
    // MKL provides all 3 name manglings (foo, foo_, FOO); pick foo_.
    #ifndef LAPACK_GLOBAL
    #define LAPACK_GLOBAL(lcname,UCNAME)  lcname##_
    #endif
#else
    #include <cblas.h>
    #include <lapacke.h>

    // Intel mkl_cblas.h does: typedef enum {...} CBLAS_ORDER;
    // Netlib    cblas.h does: enum CBLAS_ORDER {...};
    // OpenBLAS  cblas.h does: typedef enum CBLAS_ORDER {...} CBLAS_ORDER;
    // We use (CBLAS_ORDER), so add these typedefs for Netlib.
    #ifndef OPENBLAS_VERSION
    typedef enum CBLAS_ORDER CBLAS_ORDER;
    typedef enum CBLAS_TRANSPOSE CBLAS_TRANSPOSE;
    typedef enum CBLAS_UPLO CBLAS_UPLO;
    typedef enum CBLAS_DIAG CBLAS_DIAG;
    typedef enum CBLAS_SIDE CBLAS_SIDE;
    #endif
#endif

#include "core_lapack_s.h"
#include "core_lapack_d.h"
#include "core_lapack_c.h"
#include "core_lapack_z.h"

#endif // ICL_CORE_LAPACK_H

max_stars_repo_name: NLAFET/plasma

max_issues_repo_name: NLAFET/plasma

max_forks_repo_name: NLAFET/plasma


text:
#pragma once

#include "iobserver.h"
#include "irequest.h"
#include <gsl/pointers>
#include <atomic>
#include <vector>
#include <Windows.h>

namespace died
{
	namespace observer
	{
		namespace callback
		{
			unsigned WINAPI start_thread_proc(LPVOID);
			VOID CALLBACK terminate_proc(__in ULONG_PTR);
			VOID CALLBACK add_directory_proc(__in ULONG_PTR);
		}
	}

	class observer_impl final : public iobserver
	{
	public:
		observer_impl(idirectory_watcher*);

		// diable copy
		observer_impl(observer_impl const&) = delete;
		observer_impl& operator=(observer_impl const&) = delete;

		// interface
	private:
		unsigned int do_inc_request() final;
		unsigned int do_dec_request() final;
		idirectory_watcher* do_get_watcher() const final;

	private:
		void run();
		bool terminated() const;
		bool empty_request() const;
		bool add_directory(irequest* pBlock);
		void request_termination();

		friend unsigned WINAPI observer::callback::start_thread_proc(LPVOID);
		friend VOID CALLBACK observer::callback::terminate_proc(__in ULONG_PTR);
		friend VOID CALLBACK observer::callback::add_directory_proc(__in ULONG_PTR);

	private:
		gsl::not_null<idirectory_watcher*> mDirWatcher;
		std::atomic_bool mTerminated{};
		std::atomic_uint mOutstandingRequests{};
		std::vector<gsl::not_null<irequest*>> mBlocks;
	};
}

max_stars_repo_name: pvthuyet/file-watcher-demo

max_issues_repo_name: pvthuyet/file-watcher-demo

max_forks_repo_name: pvthuyet/file-watcher-demo


text:
#include <stdio.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>

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 <string_view>
#include <gsl/span>
#include <spirv_cross.hpp>

namespace glslang
{
    class TShader;
}

namespace Babylon
{
    class ShaderCompiler
    {
    public:
        ShaderCompiler();
        ~ShaderCompiler();

        struct ShaderInfo
        {
            std::unique_ptr<const spirv_cross::Compiler> Compiler;
            gsl::span<uint8_t> Bytes;
        };

        void Compile(std::string_view vertexSource, std::string_view fragmentSource, std::function<void(ShaderInfo, ShaderInfo)> onCompiled);

    protected:
        // Invert dFdy operands similar to bgfx_shader.sh
        // https://github.com/bkaradzic/bgfx/blob/7be225bf490bb1cd231cfb4abf7e617bf35b59cb/src/bgfx_shader.sh#L44-L45
        // https://github.com/bkaradzic/bgfx/blob/7be225bf490bb1cd231cfb4abf7e617bf35b59cb/src/bgfx_shader.sh#L62-L65
        static void InvertYDerivativeOperands(glslang::TShader& shader);
    };
}


text:
#ifndef GUARD_QDWF_H
#define GUARD_QDWF_H

#include "Numerics/Spline3D/Config.h"
#include "Numerics/Spline3D/uGrid1D.h"
#include "Numerics/Spline3D/CubicSpline.h"
#include "Optimize/VarList.h"
#include <gsl/gsl_sf_airy.h>
#include <math.h>
#include <iostream>

using namespace std;

class QDwf
{

  std::vector<double> X;
  std::vector<double> sigx;
  std::vector<double> a;

  double sigy;//var
  double Y0;//var

  double z0;//var
  double Ez;//var

  double VB;
  double zB;
  double kappa;
  double EL;
  double psiB;

  double onethird;

  double fac1,fac2;

  uGrid1D* m_grid;
  CubicSpline* m_spline;

  gsl_mode_t mode;

  int funcz(double z)
  {
    double denom = gsl_sf_airy_Ai(z,mode);
    double numer = fac2*gsl_sf_airy_Ai_deriv(z,mode);
    double xsign = numer-denom*kappa;
    if(xsign < 0 )
    {
      return -1;
    }
    else
    {
      return 1;
    }
  }

  inline double zprime(double z)
  {
    return fac2*(z-z0-EL/Ez);
  }

public:

  typedef double value_type;

  QDwf()
  {
    VB = 24.535;
    zB = 9.1854;
    //    z0 = 10.6858;
    //    Ez = 2.2;
    Ez = 2.04866;
    z0 = 10.8049;
    //    sigxr = 0.08;
    //    xr = 80.72;
    sigy = 0.107;
    Y0 = 25.3;
    /// set the x parameters
    int nxp = 1;
    a.resize(nxp);
    sigx.resize(nxp);
    X.resize(nxp);
    a[0] = 1.00000;
    sigx[0] = 0.08;
    X[0] = 80.628;
    //    a[1] = 1.00e-4; sigx[1] = 0.08;  X[1] = 62;
    /*
    a[0] = 1.00000; sigx[0] = 0.08; X[0] = 80.628;
    a[1] = 0.08260; sigx[1] = 0.08; X[1] = 76.000;
    a[2] = 2.17e-4; sigx[2] = 0.04; X[2] = 67.000;
    */
    /// optimised parameters
    /*
    a[0] = 1.00000; sigx[0] = 0.092403; X[0] = 80.569;
    a[1] = 0.071703; sigx[1] = 0.18445; X[1] = 76.064;
    a[2] = 2.86e-3; sigx[2] = 0.052224; X[2] = 67.000;
    sigy = 0.11112;
    Y0 = 25.167;
    */
    onethird = 1.0/3.0;
    reset();
    cout << "Airy Solved: " << EL << '\t' << kappa << '\t'
         << psiB << endl;
    /*
    /// create and initialise the uniform one-dimensional uGrid1D
    if( !m_grid ) m_grid = new uGrid1D;
    double zstart = zB; double zend = 16.00; int npts = 1000;
    m_grid->init(zstart,zend,npts);

    /// create and initialise the CubicSpline
    if(!m_spline) m_spline = new CubicSpline(m_grid);
    for(int i = 0; i < npts; i++){
      double z = m_grid->m_coord[i]; double zz = zprime(z);
      //(*m_spline)(i) = gsl_sf_airy_Ai(zz,mode);//evaluate(z);
      m_spline->F(i)[0] = gsl_sf_airy_Ai(zz,mode);
      m_spline->F(i)[1] = fac2*gsl_sf_airy_Ai_deriv(zz,mode);
    }
    */
  }

  template<class T1>
  void put(xmlNodePtr cur, VarRegistry<T1>& vlist)
  {
    //    vlist.add("C_sigxr",&sigxr,1);
    //    vlist.add("C_xr",&xr,1);
    vlist.add("C_X",&X[0],1);
    vlist.add("C_sigx",&sigx[0],1);
    vlist.add("C_a",&a[0],1);
    vlist.add("C_sigy",&sigy,1);
    vlist.add("C_Y0",&Y0,1);
    vlist.add("C_Ez",&Ez,1);
    vlist.add("C_z0",&z0,1);
  }

  void reset();

  inline void set_point(const posvec_t& r) {}

  inline double evaluate(const posvec_t& r,
                         posvec_t& gradf,
                         double& lapf)
  {
    int nsize = sigx.size();
    /// \f$\Psi_{x}\f$ sum of Gaussians
    double psix = 0.0,gfx=0.0,lfx=0.0;
    for(int i = 0; i < nsize; i++)
    {
      double x = r[0] - X[i];
      double sx = sigx[i] * x;
      double sxx = sx * x;
      double phix = a[i] * exp(-sxx);
      psix += phix;
      gfx += -2 * sx * phix;
      lfx += -2 * sigx[i] * ( 1.0 - 2.0 * sxx ) * phix;
    }
    /// \f$\Psi_{y}\f$ single Gaussian
    double y = r[1] - Y0;
    double sy = sigy * y;
    double syy = sy * y;
    double psiy = exp(-syy);
    double gfy = - 2 * sy * psiy;
    double lfy = - 2 * sigy * ( 1.0 - 2.0 * syy ) * psiy;
    /// \f$\Psi_{z}\f$ AiryAi with exponential tail
    double z = r[2];
    double psiz,gfz,lfz;
    if( z <= zB )
    {
      psiz = psiB * exp ( kappa * ( z - zB ) );
      gfz = kappa * psiz;
      lfz = kappa * gfz;
    }
    else
    {
      psiz = m_spline->evaluate(z,gfz,lfz);
    }
    double psi = psix * psiy * psiz;
    gradf[0] = gfx * psiy * psiz;
    gradf[1] = psix * gfy * psiz;
    gradf[2] = psix * gfy * psiz;
    lapf = lfx * psiy * psiz + psix * lfy * psiz + psix * psiy * lfz;
    return psi;
  }

};
#endif

max_stars_repo_name: FeiyangJin/hclib

max_issues_repo_name: FeiyangJin/hclib

max_forks_repo_name: FeiyangJin/hclib


text:
/* $Id$   */
/*--------------------------------------------------------------------*/
/*;  Copyright (C) 2003-2018                                          */
/*;  Associated Universities, Inc. Washington DC, USA.                */
/*;                                                                   */
/*;  This program is free software; you can redistribute it and/or    */
/*;  modify it under the terms of the GNU General Public License as   */
/*;  published by the Free Software Foundation; either version 2 of   */
/*;  the License, or (at your option) any later version.              */
/*;                                                                   */
/*;  This program is distributed in the hope that it will be useful,  */
/*;  but WITHOUT ANY WARRANTY; without even the implied warranty of   */
/*;  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    */
/*;  GNU General Public License for more details.                     */
/*;                                                                   */
/*;  You should have received a copy of the GNU General Public        */
/*;  License along with this program; if not, write to the Free       */
/*;  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,     */
/*;  MA 02139, USA.                                                   */
/*;                                                                   */
/*;Correspondence about this software should be addressed as follows: */
/*;         Internet email: bcotton@nrao.edu.                         */
/*;         Postal address: William Cotton                            */
/*;                         National Radio Astronomy Observatory      */
/*;                         520 Edgemont Road                         */
/*;                         Charlottesville, VA 22903-2475 USA        */
/*--------------------------------------------------------------------*/
#ifndef OBITFARRAY_H
#define OBITFARRAY_H

#include "Obit.h"
#include "ObitErr.h"
#include "ObitInfoList.h"
#include "ObitThread.h"
#if HAVE_GSL==1  /* GSL stuff */
#include <gsl/gsl_randist.h>
#endif /* HAVE_GSL */

/*-------- Obit: Merx mollis mortibus nuper ------------------*/
/**
 * \file ObitFArray.h
 * ObitFArray numeric array class definition.
 *
 * This class is derived from the #Obit class.
 * Related functions are in the
 * \link ObitFArrayUtil.h ObitFArrayUtil
 * \endlink module.
 *
 * This class is for creating and manipulating a Array as a memory resident
 * multidimensional rectangular array of floats.
 * Elements are stored in order of the increasing axis order (the reverse of the
 * usual c definition).
 * Except as noted, magic value blanking is supported (OBIT_MAGIC).
 *
 * \section ObitFArrayaccess Creators and Destructors
 * An ObitFArray will usually be created using ObitFArrayCreate which allows
 * specifying a name for the object as well as dimensionality of the array.
 *
 * A copy of a pointer to an ObitFArray should always be made using the
 * #ObitFArrayRef function which updates the reference count in the object.
 * Then whenever freeing an ObitFArray or changing a pointer, the function
 * #ObitFArrayUnref will decrement the reference count and destroy the object
 * when the reference count hits 0.
 * There is no explicit destructor.
 */

/*--------------Class definitions-------------------------------------*/
/** ObitFArray Class structure. */
typedef struct {
#include "ObitFArrayDef.h"   /* this class definition */
} ObitFArray;

/*----------------- Macroes ---------------------------*/
/**
 * Macro to unreference (and possibly destroy) an ObitFArray
 * returns a ObitFArray*.
 * in = object to unreference
 */
#define ObitFArrayUnref(in) ObitUnref (in)

/**
 * Macro to reference (update reference count) an ObitFArray.
 * returns a ObitFArray*.
 * in = object to reference
 */
#define ObitFArrayRef(in) ObitRef (in)

/**
 * Macro to determine if an object is the member of this or a
 * derived class.
 * Returns TRUE if a member, else FALSE
 * in = object to reference
 */
#define ObitFArrayIsA(in) ObitIsA (in, ObitFArrayGetClass())

/** Maximum ObitFArray number of dimensions */
#ifndef MAXFARRAYDIM
#define MAXFARRAYDIM 10
#endif

/*---------------Public functions---------------------------*/
/** Public: Class initializer. */
void ObitFArrayClassInit (void);

/** Public: Default Constructor. */
ObitFArray* newObitFArray (gchar* name);

/** Public: Create/initialize ObitFArray structures */
ObitFArray* ObitFArrayCreate (gchar* name, olong ndim, olong *naxis);
/** Typedef for definition of class pointer structure */
typedef void (*ObitFArrayCreateFP) (gchar* name, olong ndim, olong *naxis);

/** Public: ClassInfo pointer */
gconstpointer ObitFArrayGetClass (void);

/** Public: Copy (deep) constructor. */
ObitFArray* ObitFArrayCopy  (ObitFArray *in, ObitFArray *out, ObitErr *err);

/** Public: Copy structure. */
void ObitFArrayClone (ObitFArray *in, ObitFArray *out, ObitErr *err);

/** Public: Copy Subarray constructor. */
ObitFArray* ObitFArraySubArr  (ObitFArray *in, olong *blc, olong *trc,
			       ObitErr *err);
typedef ObitFArray* (*ObitFArraySubArrFP) (ObitFArray *in, olong *blc, olong *trc,
					   ObitErr *err);

/** Public: Transpose constructor. */
ObitFArray* ObitFArrayTranspose  (ObitFArray *in, olong *order, ObitErr *err);
typedef ObitFArray* (*ObitFArrayTransposeFP) (ObitFArray *in, olong *order,
					      ObitErr *err);

/** Public: Are two FArrays of compatable geometry. */
gboolean ObitFArrayIsCompatable  (ObitFArray *in1, ObitFArray *in2);
typedef gboolean (*ObitFArrayIsCompatableFP) (ObitFArray *in1, ObitFArray *in2);

/** Public: Reallocate/initialize ObitFArray structures */
ObitFArray* ObitFArrayRealloc (ObitFArray* in, olong ndim, olong *naxis);
typedef void (*ObitFArrayReallocFP) (ObitFArray* in, olong ndim, olong *naxis);

/** Public: return pointer to a specified element */
ofloat* ObitFArrayIndex (ObitFArray* in, olong *pos);
typedef ofloat* (*ObitFArrayIndexFP) (ObitFArray* in, olong *pos);

/** Public: Find Maximum value in an ObitFArray */
ofloat ObitFArrayMax (ObitFArray* in, olong *pos);
typedef ofloat (*ObitFArrayMaxFP) (ObitFArray* in, olong *pos);

/** Public: Find Maximum abs value in an ObitFArray */
ofloat ObitFArrayMaxAbs (ObitFArray* in, olong *pos);
typedef ofloat (*ObitFArrayMaxAbsFP) (ObitFArray* in, olong *pos);

/** Public: Find Minimum value in an ObitFArray */
ofloat ObitFArrayMin (ObitFArray* in, olong *pos);
typedef ofloat (*ObitFArrayMinFP) (ObitFArray* in, olong *pos);

/** Public: Replace blanks in an ObitFArray */
void ObitFArrayDeblank (ObitFArray* in, ofloat scalar);
typedef void (*ObitFArrayDeblankFP) (ObitFArray* in, ofloat scalar);

/** Public: RMS of pixel distribution from histogram */
ofloat ObitFArrayRMS (ObitFArray* in);
typedef ofloat (*ObitFArrayRMSFP) (ObitFArray* in);

/** Public: RMS of pixel distribution. */
ofloat ObitFArrayRawRMS (ObitFArray* in);
typedef ofloat (*ObitFArrayRawRMSFP) (ObitFArray* in);

/** Public: RMS of pixel about zero. */
ofloat ObitFArrayRMS0 (ObitFArray* in);
typedef ofloat (*ObitFArrayRMS0FP) (ObitFArray* in);

/** Public: Function for combining planes of multi frequency images */
ofloat ObitFArrayRMS0ST (ObitFArray* in);
typedef ofloat (*ObitFArrayRMS0STFP) (ObitFArray* in);

/** Public: RMS of pixel in potentially quantized image. */
ofloat ObitFArrayRMSQuant (ObitFArray* in);
typedef ofloat (*ObitFArrayRMSQuantFP) (ObitFArray* in);

/** Public: Determine quantization and offset in an image */
void ObitFArrayQuant (ObitFArray* in, ofloat *quant, ofloat *zero);
typedef void (*ObitFArrayQuantFP) (ObitFArray* in, ofloat *quant, ofloat *zero);

/** Public: Mode of pixel distribution. */
ofloat ObitFArrayMode (ObitFArray* in);
typedef ofloat (*ObitFArrayModeFP) (ObitFArray* in);

/** Public: Mean of pixel distribution. */
ofloat ObitFArrayMean (ObitFArray* in);
typedef ofloat (*ObitFArrayMeanFP) (ObitFArray* in);

/** Public: fill elements of an FArray */
void ObitFArrayFill (ObitFArray* in, ofloat scalar);
typedef void (*ObitFArrayFillFP) (ObitFArray* in, ofloat scalar);

/** Public: negate elements of an FArray */
void ObitFArrayNeg (ObitFArray* in);
typedef void (*ObitFArrayNegFP) (ObitFArray* in);

/** Public: taks absolute value of  elements of an FArray */
void ObitFArrayAbs (ObitFArray* in);
typedef void (*ObitFArrayAbsFP) (ObitFArray* in);

/** Public: sine of  elements of an FArray */
void ObitFArraySin (ObitFArray* in);
typedef void (*ObitFArraySinFP) (ObitFArray* in);

/** Public: cosine of  elements of an FArray */
void ObitFArrayCos (ObitFArray* in);
typedef void (*ObitFArrayCosFP) (ObitFArray* in);

/** Public: sine/cosine of  elements of an FArray */
void ObitFArraySinCos (ObitFArray* in, ObitFArray* outS, ObitFArray* outC);
typedef void (*ObitFArraySinCosFP) (ObitFArray* in, ObitFArray* outS, ObitFArray* outC);

/** Public: square root of  elements of an FArray */
void ObitFArraySqrt (ObitFArray* in);
typedef void (*ObitFArraySqrtFP) (ObitFArray* in);

/** Public: sum elements of an FArray */
ofloat ObitFArraySum (ObitFArray* in);
typedef ofloat (*ObitFArraySumFP) (ObitFArray* in);

/** Public: number of valid elements in an FArray */
olong ObitFArrayCount (ObitFArray* in);
typedef olong (*ObitFArrayCountFP) (ObitFArray* in);

/** Public: Add a scalar to elements of an FArray */
void ObitFArraySAdd (ObitFArray* in, ofloat scalar);
typedef void (*ObitFArraySAddFP) (ObitFArray* in, ofloat scalar);

/** Public: Multiply elements of an FArray by a scalar*/
void ObitFArraySMul (ObitFArray* in, ofloat scalar);
typedef void (*ObitFArraySMulFP) (ObitFArray* in, ofloat scalar);

/** Public: Divide elements of an FArray into a scalar*/
void ObitFArraySDiv (ObitFArray* in, ofloat scalar);
typedef void (*ObitFArraySDivFP) (ObitFArray* in, ofloat scalar);

/** Public: Clip elements of an FArray outside of a given range */
void ObitFArrayClip (ObitFArray* in, ofloat minVal, ofloat maxVal, ofloat newVal);
typedef void (*ObitFArrayClipFP) (ObitFArray* in, ofloat minVal, ofloat maxVal,
				  ofloat newVal);

/** Public: Clip elements of an FArray inside of a given range */
void ObitFArrayInClip (ObitFArray* in, ofloat minVal, ofloat maxVal, ofloat newVal);
typedef void (*ObitFArrayInClipFP) (ObitFArray* in, ofloat minVal, ofloat maxVal,
				    ofloat newVal);

/** Public: Blank elements of an array where another is blanked */
void ObitFArrayBlank (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayBlankFP) (ObitFArray* in1, ObitFArray* in2,
				  ObitFArray* out);

/** Public: Get larger elements of two FArrays */
void ObitFArrayMaxArr (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayMaxArrFP) (ObitFArray* in1, ObitFArray* in2,
				    ObitFArray* out);

/** Public: Get lesser elements of two FArrays */
void ObitFArrayMinArr (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayMinArrFP) (ObitFArray* in1, ObitFArray* in2,
				    ObitFArray* out);

/** Public: Get more extreme elements of two FArrays */
void ObitFArrayExtArr (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayExtArrFP) (ObitFArray* in1, ObitFArray* in2,
				    ObitFArray* out);

/** Public: Sum nonblanked elements of two FArrays */
void ObitFArraySumArr (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArraySumArrFP) (ObitFArray* in1, ObitFArray* in2,
				    ObitFArray* out);

/** Public: Average nonblanked elements of two FArrays */
void ObitFArrayAvgArr (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayAvgArrFP) (ObitFArray* in1, ObitFArray* in2,
				    ObitFArray* out);

/** Public: Add elements of two FArrays */
void ObitFArrayAdd (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayAddFP) (ObitFArray* in1, ObitFArray* in2,
				  ObitFArray* out);

/** Public: Abs Add elements of two FArrays */
void ObitFArrayAddAbs (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayAddAbsFP) (ObitFArray* in1, ObitFArray* in2,
			  	  ObitFArray* out);

/** Public: Subtract elements of two FArrays */
void ObitFArraySub (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArraySubFP) (ObitFArray* in1, ObitFArray* in2,
				  ObitFArray* out);
/** Public: copy elements from one FArray to another */
void ObitFArrayCopyData (ObitFArray* in, ObitFArray* out);
typedef void (*ObitFArrayCopyDataFP) (ObitFArray* in1, ObitFArray* out);

/** Public: Give the elements of one array the sign of the other */
void ObitFArraySign (ObitFArray* in1, ObitFArray* in2);
typedef void (*ObitFArraySignFP) (ObitFArray* in1, ObitFArray* in2);

/** Public: Multiply elements of two FArrays */
void ObitFArrayMul (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayMulFP) (ObitFArray* in1, ObitFArray* in2,
				  ObitFArray* out);

/** Public: Divide elements of two FArrays */
void ObitFArrayDiv (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void (*ObitFArrayDivFP) (ObitFArray* in1, ObitFArray* in2,
				  ObitFArray* out);

/** Public: Divide elements of two FArrays with clipping*/
void ObitFArrayDivClip (ObitFArray* in1, ObitFArray* in2, ofloat minVal, ObitFArray* out);
typedef void (*ObitFArrayDivClipFP) (ObitFArray* in1, ObitFArray* in2,
				     ofloat minVal, ObitFArray* out);

/** Public: "Dot" product to two arrays */
ofloat ObitFArrayDot (ObitFArray* in1, ObitFArray* in2);
typedef ofloat (*ObitFArrayDotFP) (ObitFArray* in1, ObitFArray* in2);

/** Public: Multiply a 2D array by a Col vector * Row vector */
void ObitFArrayMulColRow (ObitFArray* in, ObitFArray* row, ObitFArray* col,
			  ObitFArray* out);
typedef void (*ObitFArrayMulColRowFP) (ObitFArray* in, ObitFArray* row,
				       ObitFArray* col, ObitFArray* out);

/** Public: Convert a 1D "center at edges" array to proper order */
void ObitFArray1DCenter (ObitFArray* in);
typedef void (*ObitFArray1DCenterFP) (ObitFArray* in);

/** Public: Convert a 2D "center at edges" array to proper order */
void ObitFArray2DCenter (ObitFArray* in);
typedef void (*ObitFArray2DCenterFP) (ObitFArray* in);

/** Public: inplace invert a symmetric 2D array */
void ObitFArray2DSymInv (ObitFArray* in, olong *ierr);
typedef void (*ObitFArray2DSymInvFP) (ObitFArray* in, olong *ierr);

/** Public: Make 2-D Circular Gaussian in FArray */
void ObitFArray2DCGauss (ObitFArray* in, olong Cen[2], ofloat FWHM);
typedef void (*ObitFArray2DCGaussFP) (ObitFArray* in, olong Cen[2], ofloat FWHM);

/** Public: Make 2-D Eliptical Gaussian in FArray */
void ObitFArray2DEGauss (ObitFArray* in, ofloat amp, ofloat Cen[2], ofloat GauMod[3]);
typedef void (*ObitFArray2DEGaussFP) (ObitFArray* in, ofloat amp, ofloat Cen[2],
				      ofloat GauMod[3] );

/** Public: Shift and Add scaled array */
void ObitFArrayShiftAdd (ObitFArray* in1, olong *pos1,
			 ObitFArray* in2, olong *pos2,
			 ofloat scalar, ObitFArray* out);
typedef void
(*ObitFArrayShiftAddFP) (ObitFArray* in1, olong *pos1,
			 ObitFArray* in2, olong *pos2,
			 ofloat scalar, ObitFArray* out);

/** Public: Shift and Add scaled array, no threading */
void ObitFArrayShiftAddNT (ObitFArray* in1, olong *pos1,
			   ObitFArray* in2, olong *pos2,
			   ofloat scalar, ObitFArray* out);
/** Public: Zero pad an array */
void  ObitFArrayPad (ObitFArray* in, ObitFArray* out, ofloat factor);
typedef void (*ObitFArrayPadFP) (ObitFArray* in, ObitFArray* out,
				 ofloat factor);

/** Public: Convolve a list of Gaussians onto an FArray */
void  ObitFArrayConvGaus (ObitFArray* in, ObitFArray* list, olong ncomp,
			  ofloat gauss[3]);
typedef void (*ObitFArrayConvGausFP) (ObitFArray* in, ObitFArray* list,
				      olong ncomp, ofloat gauss[3]);

/** Public: Select elements in an FArray by increment */
void  ObitFArraySelInc (ObitFArray* in, ObitFArray* out, olong *blc, olong *trc,
			olong* inc, ObitErr *err);
typedef void (*ObitFArraySelIncFP) (ObitFArray* in, ObitFArray* out,
				    olong *blc, olong *trc, olong* inc, ObitErr *err);

/** Public: return histogram of elements in an FArray */
ObitFArray*  ObitFArrayHisto (ObitFArray* in, olong n, ofloat min, ofloat max);
typedef ObitFArray*  (*ObitFArrayHistoFP) (ObitFArray* in, olong n, ofloat min, ofloat max);

/** Public: exponentiate elements in an FArray */
void ObitFArrayExp (ObitFArray* in, ObitFArray* out);
typedef void  (*ObitFArrayExpFP) (ObitFArray* in, ObitFArray* out);

/** Public: natural log of elements in an FArray */
void ObitFArrayLog (ObitFArray* in, ObitFArray* out);
typedef void  (*ObitFArrayLogFP) (ObitFArray* in, ObitFArray* out);

/** Public: natural log of elements in an FArray */
void ObitFArrayPow (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);
typedef void  (*ObitFArrayPowFP) (ObitFArray* in1, ObitFArray* in2, ObitFArray* out);

/** Public: Gaussian distributed random numbers */
ofloat ObitFArrayRandom (ofloat mean, ofloat sigma);
typedef ofloat  (*ObitFArrayRandomFP) (ofloat mean, ofloat sigma);

/** Public: Fill with Gaussian distributed random numbers */
void ObitFArrayRandomFill (ObitFArray* in, ofloat mean, ofloat sigma);
typedef void  (*ObitFArrayRandomFillFP) (ObitFArray* in, ofloat mean, ofloat sigma);

/*----------- ClassInfo Structure -----------------------------------*/
/**
 * ClassInfo Structure.
 * Contains class name, a pointer to any parent class
 * (NULL if none) and function pointers.
 */
typedef struct  {
#include "ObitFArrayClassDef.h"
} ObitFArrayClassInfo;

#endif /* OBITFARRAY_H */

max_stars_repo_name: sarrvesh/Obit

max_issues_repo_name: sarrvesh/Obit

max_forks_repo_name: sarrvesh/Obit