text:
/////////////////////////////////////////////////////////////////////////////////
//
//  Solution of linear systems involved in the Levenberg - Marquardt
//  minimization algorithm
//  Copyright (C) 2004  Manolis Lourakis (lourakis at ics forth gr)
//  Institute of Computer Science, Foundation for Research & Technology - Hellas
//  Heraklion, Crete, Greece.
//
//  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.
//
/////////////////////////////////////////////////////////////////////////////////

/* Solvers for the linear systems Ax=b. Solvers should NOT modify their A & B arguments! */

#ifndef LM_REAL // not included by Axb.c
#error This file should not be compiled directly!
#endif

#ifdef LINSOLVERS_RETAIN_MEMORY
#define __STATIC__ static
#else
#define __STATIC__ // empty
#endif /* LINSOLVERS_RETAIN_MEMORY */

#ifdef HAVE_LAPACK

/* prototypes of LAPACK routines */

#define GEQRF LM_MK_LAPACK_NAME(geqrf)
#define ORGQR LM_MK_LAPACK_NAME(orgqr)
#define TRTRS LM_MK_LAPACK_NAME(trtrs)
#define POTF2 LM_MK_LAPACK_NAME(potf2)
#define POTRF LM_MK_LAPACK_NAME(potrf)
#define POTRS LM_MK_LAPACK_NAME(potrs)
#define GETRF LM_MK_LAPACK_NAME(getrf)
#define GETRS LM_MK_LAPACK_NAME(getrs)
#define GESVD LM_MK_LAPACK_NAME(gesvd)
#define GESDD LM_MK_LAPACK_NAME(gesdd)
#define SYTRF LM_MK_LAPACK_NAME(sytrf)
#define SYTRS LM_MK_LAPACK_NAME(sytrs)
#define PLASMA_POSV LM_CAT_(PLASMA_, LM_ADD_PREFIX(posv))

#ifdef __cplusplus
extern "C" {
#endif
/* QR decomposition */
extern int GEQRF(int *m, int *n, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info);
extern int ORGQR(int *m, int *n, int *k, LM_REAL *a, int *lda, LM_REAL *tau, LM_REAL *work, int *lwork, int *info);

/* solution of triangular systems */
extern int TRTRS(char *uplo, char *trans, char *diag, int *n, int *nrhs, LM_REAL *a, int *lda, LM_REAL *b, int *ldb, int *info);

/* Cholesky decomposition and systems solution */
extern int POTF2(char *uplo, int *n, LM_REAL *a, int *lda, int *info);
extern int POTRF(char *uplo, int *n, LM_REAL *a, int *lda, int *info); /* block version of dpotf2 */
extern int POTRS(char *uplo, int *n, int *nrhs, LM_REAL *a, int *lda, LM_REAL *b, int *ldb, int *info);

/* LU decomposition and systems solution */
extern int GETRF(int *m, int *n, LM_REAL *a, int *lda, int *ipiv, int *info);
extern int GETRS(char *trans, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);

/* Singular Value Decomposition (SVD) */
extern int GESVD(char *jobu, char *jobvt, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu,
                   LM_REAL *vt, int *ldvt, LM_REAL *work, int *lwork, int *info);

/* lapack 3.0 new SVD routine, faster than xgesvd().
 * In case that your version of LAPACK does not include them, use the above two older routines
 */
extern int GESDD(char *jobz, int *m, int *n, LM_REAL *a, int *lda, LM_REAL *s, LM_REAL *u, int *ldu, LM_REAL *vt, int *ldvt,
                   LM_REAL *work, int *lwork, int *iwork, int *info);

/* LDLt/UDUt factorization and systems solution */
extern int SYTRF(char *uplo, int *n, LM_REAL *a, int *lda, int *ipiv, LM_REAL *work, int *lwork, int *info);
extern int SYTRS(char *uplo, int *n, int *nrhs, LM_REAL *a, int *lda, int *ipiv, LM_REAL *b, int *ldb, int *info);
#ifdef __cplusplus
}
#endif

/* precision-specific definitions */
#define AX_EQ_B_QR LM_ADD_PREFIX(Ax_eq_b_QR)
#define AX_EQ_B_QRLS LM_ADD_PREFIX(Ax_eq_b_QRLS)
#define AX_EQ_B_CHOL LM_ADD_PREFIX(Ax_eq_b_Chol)
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU)
#define AX_EQ_B_SVD LM_ADD_PREFIX(Ax_eq_b_SVD)
#define AX_EQ_B_BK LM_ADD_PREFIX(Ax_eq_b_BK)
#define AX_EQ_B_PLASMA_CHOL LM_ADD_PREFIX(Ax_eq_b_PLASMA_Chol)

/*
 * This function returns the solution of Ax = b
 *
 * The function is based on QR decomposition with explicit computation of Q:
 * If A=Q R with Q orthogonal and R upper triangular, the linear system becomes
 * Q R x = b or R x = Q^T b.
 * The last equation can be solved directly.
 *
 * A is mxm, b is mx1
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_QR(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;

static int nb=0; /* no __STATIC__ decl. here! */

LM_REAL *a, *tau, *r, *work;
int a_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
      return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

    /* calculate required memory size */
    a_sz=m*m;
    tau_sz=m;
    r_sz=m*m; /* only the upper triangular part really needed */
    if(!nb){
      LM_REAL tmp;

      worksz=-1; // workspace query; optimal size is returned in tmp
      GEQRF((int *)&m, (int *)&m, NULL, (int *)&m, NULL, (LM_REAL *)&tmp, (int *)&worksz, (int *)&info);
      nb=((int)tmp)/m; // optimal worksize is m*nb
    }
    worksz=nb*m;
    tot_sz=a_sz + tau_sz + r_sz + worksz;

#ifdef LINSOLVERS_RETAIN_MEMORY
    if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
      if(buf) free(buf); /* free previously allocated memory */

      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n");
        exit(1);
      }
    }
#else
      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_QR) "() failed!\n");
        exit(1);
      }
#endif /* LINSOLVERS_RETAIN_MEMORY */

    a=buf;
    tau=a+a_sz;
    r=tau+tau_sz;
    work=r+r_sz;

  /* store A (column major!) into a */
	for(i=0; i<m; i++)
		for(j=0; j<m; j++)
			a[i+j*m]=A[i*m+j];

  /* QR decomposition of A */
  GEQRF((int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GEQRF) " in ", AX_EQ_B_QR) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* R is stored in the upper triangular part of a; copy it in r so that ORGQR() below won't destroy it */
  memcpy(r, a, r_sz*sizeof(LM_REAL));

  /* compute Q using the elementary reflectors computed by the above decomposition */
  ORGQR((int *)&m, (int *)&m, (int *)&m, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", ORGQR) " in ", AX_EQ_B_QR) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("Unknown LAPACK error (%d) in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* Q is now in a; compute Q^T b in x */
  for(i=0; i<m; i++){
    for(j=0, sum=0.0; j<m; j++)
      sum+=a[i*m+j]*B[j];
    x[i]=sum;
  }

  /* solve the linear system R x = Q^t b */
  TRTRS("U", "N", "N", (int *)&m, (int *)&nrhs, r, (int *)&m, x, (int *)&m, &info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QR) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QR) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

/*
 * This function returns the solution of min_x ||Ax - b||
 *
 * || . || is the second order (i.e. L2) norm. This is a least squares technique that
 * is based on QR decomposition:
 * If A=Q R with Q orthogonal and R upper triangular, the normal equations become
 * (A^T A) x = A^T b  or (R^T Q^T Q R) x = A^T b or (R^T R) x = A^T b.
 * This amounts to solving R^T y = A^T b for y and then R x = y for x
 * Note that Q does not need to be explicitly computed
 *
 * A is mxn, b is mx1
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_QRLS(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m, int n)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;

static int nb=0; /* no __STATIC__ decl. here! */

LM_REAL *a, *tau, *r, *work;
int a_sz, tau_sz, r_sz, tot_sz;
register int i, j;
int info, worksz, nrhs=1;
register LM_REAL sum;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
      return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

    if(m<n){
		  printf(RCAT("Normal equations require that the number of rows is greater than number of columns in ", AX_EQ_B_QRLS) "() [%d x %d]! -- try transposing\n", m, n);
		  exit(1);
	  }

    /* calculate required memory size */
    a_sz=m*n;
    tau_sz=n;
    r_sz=n*n;
    if(!nb){
      LM_REAL tmp;

      worksz=-1; // workspace query; optimal size is returned in tmp
      GEQRF((int *)&m, (int *)&m, NULL, (int *)&m, NULL, (LM_REAL *)&tmp, (int *)&worksz, (int *)&info);
      nb=((int)tmp)/m; // optimal worksize is m*nb
    }
    worksz=nb*m;
    tot_sz=a_sz + tau_sz + r_sz + worksz;

#ifdef LINSOLVERS_RETAIN_MEMORY
    if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
      if(buf) free(buf); /* free previously allocated memory */

      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");
        exit(1);
      }
    }
#else
      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_QRLS) "() failed!\n");
        exit(1);
      }
#endif /* LINSOLVERS_RETAIN_MEMORY */

    a=buf;
    tau=a+a_sz;
    r=tau+tau_sz;
    work=r+r_sz;

  /* store A (column major!) into a */
	for(i=0; i<m; i++)
		for(j=0; j<n; j++)
			a[i+j*m]=A[i*n+j];

  /* compute A^T b in x */
  for(i=0; i<n; i++){
    for(j=0, sum=0.0; j<m; j++)
      sum+=A[j*n+i]*B[j];
    x[i]=sum;
  }

  /* QR decomposition of A */
  GEQRF((int *)&m, (int *)&n, a, (int *)&m, tau, work, (int *)&worksz, (int *)&info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GEQRF) " in ", AX_EQ_B_QRLS) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT(RCAT("Unknown LAPACK error %d for ", GEQRF) " in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* R is stored in the upper triangular part of a. Note that a is mxn while r nxn */
  for(j=0; j<n; j++){
    for(i=0; i<=j; i++)
      r[i+j*n]=a[i+j*m];

    /* lower part is zero */
    for(i=j+1; i<n; i++)
      r[i+j*n]=0.0;
  }

  /* solve the linear system R^T y = A^t b */
  TRTRS("U", "T", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, x, (int *)&n, &info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QRLS) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* solve the linear system R x = y */
  TRTRS("U", "N", "N", (int *)&n, (int *)&nrhs, r, (int *)&n, x, (int *)&n, &info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_QRLS) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_QRLS) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

/*
 * This function returns the solution of Ax=b
 *
 * The function assumes that A is symmetric & postive definite and employs
 * the Cholesky decomposition:
 * If A=L L^T with L lower triangular, the system to be solved becomes
 * (L L^T) x = b
 * This amounts to solving L y = b for y and then L^T x = y for x
 *
 * A is mxm, b is mx1
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_CHOL(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;

LM_REAL *a;
int a_sz, tot_sz;
int info, nrhs=1;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
      return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

    /* calculate required memory size */
    a_sz=m*m;
    tot_sz=a_sz;

#ifdef LINSOLVERS_RETAIN_MEMORY
    if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
      if(buf) free(buf); /* free previously allocated memory */

      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");
        exit(1);
      }
    }
#else
      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_CHOL) "() failed!\n");
        exit(1);
      }
#endif /* LINSOLVERS_RETAIN_MEMORY */

  a=buf;

  /* store A into a and B into x. A is assumed symmetric,
   * hence no transposition is needed
   */
  memcpy(a, A, a_sz*sizeof(LM_REAL));
  memcpy(x, B, m*sizeof(LM_REAL));

  /* Cholesky decomposition of A */
  //POTF2("L", (int *)&m, a, (int *)&m, (int *)&info);
  POTRF("L", (int *)&m, a, (int *)&m, (int *)&info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTF2) "/", POTRF) " in ",
                      AX_EQ_B_CHOL) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT(RCAT(RCAT("LAPACK error: the leading minor of order %d is not positive definite,\nthe factorization could not be completed for ", POTF2) "/", POTRF) " in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* solve using the computed Cholesky in one lapack call */
  POTRS("L", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
  if(info<0){
    printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", POTRS) " in ", AX_EQ_B_CHOL) "()\n", -info);
    exit(1);
  }

#if 0
  /* alternative: solve the linear system L y = b ... */
  TRTRS("L", "N", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) " in ", AX_EQ_B_CHOL) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  /* ... solve the linear system L^T x = y */
  TRTRS("L", "T", "N", (int *)&m, (int *)&nrhs, a, (int *)&m, x, (int *)&m, &info);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", TRTRS) "in ", AX_EQ_B_CHOL) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: the %d-th diagonal element of A is zero (singular matrix) in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }
#endif /* 0 */

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

#ifdef HAVE_PLASMA

/* Linear algebra using PLASMA parallel library for multicore CPUs.
 * http://icl.cs.utk.edu/plasma/
 *
 * WARNING: BLAS multithreading should be disabled, e.g. setenv MKL_NUM_THREADS 1
 */

#ifndef _LM_PLASMA_MISC_
/* avoid multiple inclusion of helper code */
#define _LM_PLASMA_MISC_

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <plasma_tmg.h>
#include <core_blas.h>

/* programmatically determine the number of cores on the current machine */
#ifdef _WIN32
#include <windows.h>
#elif __linux
#include <unistd.h>
#endif
static int getnbcores()
{
#ifdef _WIN32
  SYSTEM_INFO sysinfo;
  GetSystemInfo(&sysinfo);
  return sysinfo.dwNumberOfProcessors;
#elif __linux
  return sysconf(_SC_NPROCESSORS_ONLN);
#else // unknown system
  return 2<<1; // will be halved by right shift below
#endif
}

static int PLASMA_ncores=-(getnbcores()>>1); // >0 if PLASMA initialized, <0 otherwise

/* user-specified number of cores */
void levmar_PLASMA_setnbcores(int cores)
{
  PLASMA_ncores=(cores>0)? -cores : ((cores)? cores : -2);
}
#endif /* _LM_PLASMA_MISC_ */

/*
 * This function returns the solution of Ax=b
 *
 * The function assumes that A is symmetric & positive definite and employs the
 * Cholesky decomposition implemented by PLASMA for homogeneous multicore processors.
 *
 * A is mxm, b is mx1
 *
 * The function returns 0 in case of error, 1 if successfull
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_PLASMA_CHOL(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;

LM_REAL *a;
int a_sz, tot_sz;
int info, nrhs=1;

    if(A==NULL){
#ifdef LINSOLVERS_RETAIN_MEMORY
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;
#endif /* LINSOLVERS_RETAIN_MEMORY */

      PLASMA_Finalize();
      PLASMA_ncores=-PLASMA_ncores;

      return 1;
    }

    /* calculate required memory size */
    a_sz=m*m;
    tot_sz=a_sz;

#ifdef LINSOLVERS_RETAIN_MEMORY
    if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
      if(buf) free(buf); /* free previously allocated memory */

      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_PLASMA_CHOL) "() failed!\n");
        exit(1);
      }
    }
#else
    buf_sz=tot_sz;
    buf=(LM_REAL *)malloc(buf_sz*sizeof(LM_REAL));
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_PLASMA_CHOL) "() failed!\n");
      exit(1);
    }
#endif /* LINSOLVERS_RETAIN_MEMORY */

    a=buf;

    /* store A into a and B into x; A is assumed to be symmetric,
     * hence no transposition is needed
     */
    memcpy(a, A, a_sz*sizeof(LM_REAL));
    memcpy(x, B, m*sizeof(LM_REAL));

  /* initialize PLASMA */
  if(PLASMA_ncores<0){
    PLASMA_ncores=-PLASMA_ncores;
    PLASMA_Init(PLASMA_ncores);
    printf(RCAT("\n", AX_EQ_B_PLASMA_CHOL) "(): PLASMA is running on %d cores.\n\n", PLASMA_ncores);
  }

  /* Solve the linear system */
  info=PLASMA_POSV(PlasmaLower, m, 1, a, m, x, m);
  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", PLASMA_POSV) " in ",
                      AX_EQ_B_PLASMA_CHOL) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT(RCAT("LAPACK error: the leading minor of order %d is not positive definite,\n"
                                "the factorization could not be completed for ", PLASMA_POSV) " in ", AX_EQ_B_CHOL) "()\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif
      return 0;
    }
  }

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}
#endif /* HAVE_PLASMA */

/*
 * This function returns the solution of Ax = b
 *
 * The function employs LU decomposition:
 * If A=L U with L lower and U upper triangular, then the original system
 * amounts to solving
 * L y = b, U x = y
 *
 * A is mxm, b is mx1
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;

int a_sz, ipiv_sz, tot_sz;
register int i, j;
int info, *ipiv, nrhs=1;
LM_REAL *a;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
      return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

    /* calculate required memory size */
    ipiv_sz=m;
    a_sz=m*m;
    tot_sz=a_sz*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */

#ifdef LINSOLVERS_RETAIN_MEMORY
    if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
      if(buf) free(buf); /* free previously allocated memory */

      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz);
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
        exit(1);
      }
    }
#else
      buf_sz=tot_sz;
      buf=(LM_REAL *)malloc(buf_sz);
      if(!buf){
        printf(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
        exit(1);
      }
#endif /* LINSOLVERS_RETAIN_MEMORY */

    a=buf;
    ipiv=(int *)(a+a_sz);

    /* store A (column major!) into a and B into x */
	  for(i=0; i<m; i++){
		  for(j=0; j<m; j++)
        a[i+j*m]=A[i*m+j];

      x[i]=B[i];
    }

  /* LU decomposition for A */
	GETRF((int *)&m, (int *)&m, a, (int *)&m, ipiv, (int *)&info);
	if(info!=0){
		if(info<0){
      printf(RCAT(RCAT("argument %d of ", GETRF) " illegal in ", AX_EQ_B_LU) "()\n", -info);
			exit(1);
		}
		else{
      printf(RCAT(RCAT("singular matrix A for ", GETRF) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

			return 0;
		}
	}

  /* solve the system with the computed LU */
  GETRS("N", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, x, (int *)&m, (int *)&info);
	if(info!=0){
		if(info<0){
			printf(RCAT(RCAT("argument %d of ", GETRS) " illegal in ", AX_EQ_B_LU) "()\n", -info);
			exit(1);
		}
		else{
			printf(RCAT(RCAT("unknown error for ", GETRS) " in ", AX_EQ_B_LU) "()\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

			return 0;
		}
	}

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

/*
 * This function returns the solution of Ax = b
 *
 * The function is based on SVD decomposition:
 * If A=U D V^T with U, V orthogonal and D diagonal, the linear system becomes
 * (U D V^T) x = b or x=V D^{-1} U^T b
 * Note that V D^{-1} U^T is the pseudoinverse A^+
 *
 * A is mxm, b is mx1.
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_SVD(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0;
static LM_REAL eps=LM_CNST(-1.0);

register int i, j;
LM_REAL *a, *u, *s, *vt, *work;
int a_sz, u_sz, s_sz, vt_sz, tot_sz;
LM_REAL thresh, one_over_denom;
register LM_REAL sum;
int info, rank, worksz, *iwork, iworksz;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
      return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

  /* calculate required memory size */
#if 1 /* use optimal size */
  worksz=-1; // workspace query. Keep in mind that GESDD requires more memory than GESVD
  /* note that optimal work size is returned in thresh */
  GESVD("A", "A", (int *)&m, (int *)&m, NULL, (int *)&m, NULL, NULL, (int *)&m, NULL, (int *)&m, (LM_REAL *)&thresh, (int *)&worksz, &info);
  //GESDD("A", (int *)&m, (int *)&m, NULL, (int *)&m, NULL, NULL, (int *)&m, NULL, (int *)&m, (LM_REAL *)&thresh, (int *)&worksz, NULL, &info);
  worksz=(int)thresh;
#else /* use minimum size */
  worksz=5*m; // min worksize for GESVD
  //worksz=m*(7*m+4); // min worksize for GESDD
#endif
  iworksz=8*m;
  a_sz=m*m;
  u_sz=m*m; s_sz=m; vt_sz=m*m;

  tot_sz=(a_sz + u_sz + s_sz + vt_sz + worksz)*sizeof(LM_REAL) + iworksz*sizeof(int); /* should be arranged in that order for proper doubles alignment */

#ifdef LINSOLVERS_RETAIN_MEMORY
  if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
    if(buf) free(buf); /* free previously allocated memory */

    buf_sz=tot_sz;
    buf=(LM_REAL *)malloc(buf_sz);
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");
      exit(1);
    }
  }
#else
    buf_sz=tot_sz;
    buf=(LM_REAL *)malloc(buf_sz);
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_SVD) "() failed!\n");
      exit(1);
    }
#endif /* LINSOLVERS_RETAIN_MEMORY */

  a=buf;
  u=a+a_sz;
  s=u+u_sz;
  vt=s+s_sz;
  work=vt+vt_sz;
  iwork=(int *)(work+worksz);

  /* store A (column major!) into a */
  for(i=0; i<m; i++)
    for(j=0; j<m; j++)
      a[i+j*m]=A[i*m+j];

  /* SVD decomposition of A */
  GESVD("A", "A", (int *)&m, (int *)&m, a, (int *)&m, s, u, (int *)&m, vt, (int *)&m, work, (int *)&worksz, &info);
  //GESDD("A", (int *)&m, (int *)&m, a, (int *)&m, s, u, (int *)&m, vt, (int *)&m, work, (int *)&worksz, iwork, &info);

  /* error treatment */
  if(info!=0){
    if(info<0){
      printf(RCAT(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", GESVD), "/" GESDD) " in ", AX_EQ_B_SVD) "()\n", -info);
      exit(1);
    }
    else{
      printf(RCAT("LAPACK error: dgesdd (dbdsdc)/dgesvd (dbdsqr) failed to converge in ", AX_EQ_B_SVD) "() [info=%d]\n", info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

      return 0;
    }
  }

  if(eps<0.0){
    LM_REAL aux;

    /* compute machine epsilon */
    for(eps=LM_CNST(1.0); aux=eps+LM_CNST(1.0), aux-LM_CNST(1.0)>0.0; eps*=LM_CNST(0.5))
                                          ;
    eps*=LM_CNST(2.0);
  }

  /* compute the pseudoinverse in a */
	for(i=0; i<a_sz; i++) a[i]=0.0; /* initialize to zero */
  for(rank=0, thresh=eps*s[0]; rank<m && s[rank]>thresh; rank++){
    one_over_denom=LM_CNST(1.0)/s[rank];

    for(j=0; j<m; j++)
      for(i=0; i<m; i++)
        a[i*m+j]+=vt[rank+i*m]*u[j+rank*m]*one_over_denom;
  }

	/* compute A^+ b in x */
	for(i=0; i<m; i++){
	  for(j=0, sum=0.0; j<m; j++)
      sum+=a[i*m+j]*B[j];
    x[i]=sum;
  }

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

/*
 * This function returns the solution of Ax = b for a real symmetric matrix A
 *
 * The function is based on LDLT factorization with the pivoting
 * strategy of Bunch and Kaufman:
 * A is factored as L*D*L^T where L is lower triangular and
 * D symmetric and block diagonal (aka spectral decomposition,
 * Banachiewicz factorization, modified Cholesky factorization)
 *
 * A is mxm, b is mx1.
 *
 * The function returns 0 in case of error, 1 if successfull
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_BK(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ LM_REAL *buf=NULL;
__STATIC__ int buf_sz=0, nb=0;

LM_REAL *a, *work;
int a_sz, ipiv_sz, work_sz, tot_sz;
int info, *ipiv, nrhs=1;

  if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
  {
    if(buf) free(buf);
    buf=NULL;
    buf_sz=0;

    return 1;
  }
#else
  return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

  /* calculate required memory size */
  ipiv_sz=m;
  a_sz=m*m;
  if(!nb){
    LM_REAL tmp;

    work_sz=-1; // workspace query; optimal size is returned in tmp
    SYTRF("L", (int *)&m, NULL, (int *)&m, NULL, (LM_REAL *)&tmp, (int *)&work_sz, (int *)&info);
    nb=((int)tmp)/m; // optimal worksize is m*nb
  }
  work_sz=(nb!=-1)? nb*m : 1;
  tot_sz=(a_sz + work_sz)*sizeof(LM_REAL) + ipiv_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */

#ifdef LINSOLVERS_RETAIN_MEMORY
  if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
    if(buf) free(buf); /* free previously allocated memory */

    buf_sz=tot_sz;
    buf=(LM_REAL *)malloc(buf_sz);
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n");
      exit(1);
    }
  }
#else
  buf_sz=tot_sz;
  buf=(LM_REAL *)malloc(buf_sz);
  if(!buf){
    printf(RCAT("memory allocation in ", AX_EQ_B_BK) "() failed!\n");
    exit(1);
  }
#endif /* LINSOLVERS_RETAIN_MEMORY */

  a=buf;
  work=a+a_sz;
  ipiv=(int *)(work+work_sz);

  /* store A into a and B into x; A is assumed to be symmetric, hence
   * the column and row major order representations are the same
   */
  memcpy(a, A, a_sz*sizeof(LM_REAL));
  memcpy(x, B, m*sizeof(LM_REAL));

  /* LDLt factorization for A */
	SYTRF("L", (int *)&m, a, (int *)&m, ipiv, work, (int *)&work_sz, (int *)&info);
	if(info!=0){
		if(info<0){
      printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", SYTRF) " in ", AX_EQ_B_BK) "()\n", -info);
			exit(1);
		}
		else{
      printf(RCAT(RCAT("LAPACK error: singular block diagonal matrix D for", SYTRF) " in ", AX_EQ_B_BK)"() [D(%d, %d) is zero]\n", info, info);
#ifndef LINSOLVERS_RETAIN_MEMORY
      free(buf);
#endif

			return 0;
		}
	}

  /* solve the system with the computed factorization */
  SYTRS("L", (int *)&m, (int *)&nrhs, a, (int *)&m, ipiv, x, (int *)&m, (int *)&info);
  if(info<0){
    printf(RCAT(RCAT("LAPACK error: illegal value for argument %d of ", SYTRS) " in ", AX_EQ_B_BK) "()\n", -info);
    exit(1);
	}

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

	return 1;
}

/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_QR
#undef AX_EQ_B_QRLS
#undef AX_EQ_B_CHOL
#undef AX_EQ_B_LU
#undef AX_EQ_B_SVD
#undef AX_EQ_B_BK
#undef AX_EQ_B_PLASMA_CHOL

#undef GEQRF
#undef ORGQR
#undef TRTRS
#undef POTF2
#undef POTRF
#undef POTRS
#undef GETRF
#undef GETRS
#undef GESVD
#undef GESDD
#undef SYTRF
#undef SYTRS
#undef PLASMA_POSV

#else // no LAPACK

/* precision-specific definitions */
#define AX_EQ_B_LU LM_ADD_PREFIX(Ax_eq_b_LU_noLapack)

/*
 * This function returns the solution of Ax = b
 *
 * The function employs LU decomposition followed by forward/back substitution (see
 * also the LAPACK-based LU solver above)
 *
 * A is mxm, b is mx1
 *
 * The function returns 0 in case of error, 1 if successful
 *
 * This function is often called repetitively to solve problems of identical
 * dimensions. To avoid repetitive malloc's and free's, allocated memory is
 * retained between calls and free'd-malloc'ed when not of the appropriate size.
 * A call with NULL as the first argument forces this memory to be released.
 */
int AX_EQ_B_LU(LM_REAL *A, LM_REAL *B, LM_REAL *x, int m)
{
__STATIC__ void *buf=NULL;
__STATIC__ int buf_sz=0;

register int i, j, k;
int *idx, maxi=-1, idx_sz, a_sz, work_sz, tot_sz;
LM_REAL *a, *work, max, sum, tmp;

    if(!A)
#ifdef LINSOLVERS_RETAIN_MEMORY
    {
      if(buf) free(buf);
      buf=NULL;
      buf_sz=0;

      return 1;
    }
#else
    return 1; /* NOP */
#endif /* LINSOLVERS_RETAIN_MEMORY */

  /* calculate required memory size */
  idx_sz=m;
  a_sz=m*m;
  work_sz=m;
  tot_sz=(a_sz+work_sz)*sizeof(LM_REAL) + idx_sz*sizeof(int); /* should be arranged in that order for proper doubles alignment */

#ifdef LINSOLVERS_RETAIN_MEMORY
  if(tot_sz>buf_sz){ /* insufficient memory, allocate a "big" memory chunk at once */
    if(buf) free(buf); /* free previously allocated memory */

    buf_sz=tot_sz;
    buf=(void *)malloc(tot_sz);
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
      exit(1);
    }
  }
#else
    buf_sz=tot_sz;
    buf=(void *)malloc(tot_sz);
    if(!buf){
      printf(RCAT("memory allocation in ", AX_EQ_B_LU) "() failed!\n");
      exit(1);
    }
#endif /* LINSOLVERS_RETAIN_MEMORY */

  a=buf;
  work=a+a_sz;
  idx=(int *)(work+work_sz);

  /* avoid destroying A, B by copying them to a, x resp. */
  memcpy(a, A, a_sz*sizeof(LM_REAL));
  memcpy(x, B, m*sizeof(LM_REAL));

  /* compute the LU decomposition of a row permutation of matrix a; the permutation itself is saved in idx[] */
	for(i=0; i<m; ++i){
		max=0.0;
		for(j=0; j<m; ++j)
			if((tmp=FABS(a[i*m+j]))>max)
        max=tmp;
		  if(max==0.0){
        printf(RCAT("Singular matrix A in ", AX_EQ_B_LU) "()!\n");
#ifndef LINSOLVERS_RETAIN_MEMORY
        free(buf);
#endif

        return 0;
      }
		  work[i]=LM_CNST(1.0)/max;
	}

	for(j=0; j<m; ++j){
		for(i=0; i<j; ++i){
			sum=a[i*m+j];
			for(k=0; k<i; ++k)
        sum-=a[i*m+k]*a[k*m+j];
			a[i*m+j]=sum;
		}
		max=0.0;
		for(i=j; i<m; ++i){
			sum=a[i*m+j];
			for(k=0; k<j; ++k)
        sum-=a[i*m+k]*a[k*m+j];
			a[i*m+j]=sum;
			if((tmp=work[i]*FABS(sum))>=max){
				max=tmp;
				maxi=i;
			}
		}
		if(j!=maxi){
			for(k=0; k<m; ++k){
				tmp=a[maxi*m+k];
				a[maxi*m+k]=a[j*m+k];
				a[j*m+k]=tmp;
			}
			work[maxi]=work[j];
		}
		idx[j]=maxi;
		if(a[j*m+j]==0.0)
      a[j*m+j]=LM_REAL_EPSILON;
		if(j!=m-1){
			tmp=LM_CNST(1.0)/(a[j*m+j]);
			for(i=j+1; i<m; ++i)
        a[i*m+j]*=tmp;
		}
	}

  /* The decomposition has now replaced a. Solve the linear system using
   * forward and back substitution
   */
	for(i=k=0; i<m; ++i){
		j=idx[i];
		sum=x[j];
		x[j]=x[i];
		if(k!=0)
			for(j=k-1; j<i; ++j)
        sum-=a[i*m+j]*x[j];
		else
      if(sum!=0.0)
			  k=i+1;
		x[i]=sum;
	}

	for(i=m-1; i>=0; --i){
		sum=x[i];
		for(j=i+1; j<m; ++j)
      sum-=a[i*m+j]*x[j];
		x[i]=sum/a[i*m+i];
	}

#ifndef LINSOLVERS_RETAIN_MEMORY
  free(buf);
#endif

  return 1;
}

/* undefine all. IT MUST REMAIN IN THIS POSITION IN FILE */
#undef AX_EQ_B_LU

#endif /* HAVE_LAPACK */

max_stars_repo_name: jongwonjlee/CLAPACK

max_issues_repo_name: jongwonjlee/CLAPACK

max_forks_repo_name: jongwonjlee/CLAPACK


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

#include <cblas.h>

#define SWAP(x, y, T) do { T SWAP = x; x = y; y = SWAP; } while (0)

//The Struct that contains the k nearest neighbors of m queries
typedef struct knnresult{
  int    * nidx;    //!< Indices (0-based) of nearest neighbors [m-by-k]
  double * ndist;   //!< Distance of nearest neighbors          [m-by-k]
  int      m;       //!< Number of query points                 [scalar]
  int      k;       //!< Number of nearest neighbors            [scalar]
} knnresult;

//A sorting alogorithm for the k first items of a list
void k_select(int * nidx,
            double * ndist,
            int k,
            int n)
{
    for(int i=0; i<k; i++){
        int minidx = i;
        double min = ndist[i];
        for(int j=i+1; j<n; j++){
            if(min > ndist[j]){
                min = ndist[j];
                minidx = nidx[j];
                SWAP(ndist[i],ndist[j], double);
                SWAP(nidx[i], nidx[j], int);
            }
        }
    }
}

int partition(int * nidx,
            double * ndist,
            int left,
            int right,
            int pivot)
{

    double val = ndist[pivot];

    SWAP(ndist[right], ndist[pivot], double);
    SWAP(nidx[right], nidx[pivot], int);

    int storeIdx = left;

    for(int i=left; i<right; i++){
        if(ndist[i] < val){
            SWAP(ndist[storeIdx], ndist[i], double);
            SWAP(nidx[storeIdx], nidx[i], int);
            storeIdx++;

        }
    }
    SWAP(ndist[storeIdx], ndist[right], double);
    SWAP(nidx[storeIdx], nidx[right], int);
    return storeIdx;
}

void quickselect(int * nidx,
            double * ndist,
            int left,
            int right,
            int k)
{
    if(left >= right) return;

    int pivot = left + rand() % (right - left + 1);
    pivot = partition(nidx, ndist, left, right, pivot);

    if(k == pivot){
        return;
    }else if(k< pivot){
        quickselect(nidx, ndist, left, pivot - 1, k);
    }else{
        quickselect(nidx, ndist, pivot + 1, right, k);
    }

}

//Calculates a column of the Euclidean distance matrix D
double *calc_Dcol(  double *X,
                    double *Y,
                    int n,
                    int d,
                    int row)
{
    double *Dcol= (double *)malloc(n * sizeof(double));

    for(int i=0; i<n; i++) Dcol[i] = 0;

    for(int i=0; i<n; i++){
        for(int k=0; k<d; k++){
            Dcol[i] += X[k + i*d]*X[k + i*d] + Y[k + row*d]*Y[k + row*d] - 2*X[k + i*d]*Y[k + row*d];
        }
    }

    return Dcol;
}

//Calculates the Euclidean distance matrix D
double *calc_D( double *X,
                double *Y,
                int n,
                int d,
                int m)
{
    double *D= (double *)malloc(n * m * sizeof(double));

    double xsum[n],
           ysum[m];
    for(int i=0; i<n; i++) xsum[i] = cblas_ddot(d, X + i*d, 1, X + i*d, 1);
    for(int i=0; i<m; i++) ysum[i] = cblas_ddot(d, Y + i*d, 1, Y + i*d, 1);

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, n, m, d, -2, X, d, Y, d, 0, D, m);

    for(int i=0; i<n; i++){
        for(int j=0; j<m; j++){
            D[j + m*i] += xsum[i] + ysum[j];
        }
    }

    return D;
}

void calc_part_knn( double *D,
                    int n,
                    int m,
                    int k,
                    int col,
                    int *final_nidx,
                    double *final_ndist)
{

    int * nidx = (int *)malloc(n * sizeof(int));
    double * ndist = (double *)malloc(n * sizeof(double));

    for(int i=0; i<n; i++){
        nidx[i] = i;
        ndist[i] = D[i * m + col];
    }

    quickselect(nidx, ndist, 0, n-1, k);

    memcpy(final_nidx, nidx, k * sizeof(int));
    memcpy(final_ndist, ndist, k * sizeof(double));

    free(nidx);
    free(ndist);
}

//Finds for each point in a query set Y the k nearest neighbors
//in the corpus set X
knnresult kNN(  double *X,
                double *Y,
                int n,
                int m,
                int d,
                int k)
{

    knnresult knn;
    knn.m = m;
    knn.k = k;
    knn.nidx = (int *)malloc(m * k * sizeof(int));
    knn.ndist = (double *)malloc(m * k * sizeof(double));

    int max_m = 1e3;
    if(m > max_m){

        int iteration_m[m/max_m];
        for(int i=0; i<m/max_m; i++) iteration_m[i] = max_m;
        iteration_m[m/max_m - 1] = max_m + m%max_m;

        int displ = 0;
        for(int i=0; i<m/max_m; i++){
            double *D = calc_D(X, Y + displ * d, n, d, iteration_m[i]);

            for(int j=0; j<iteration_m[i]; j++) calc_part_knn(D, n, iteration_m[i], k, j, knn.nidx + (j + displ) * k, knn.ndist + (j + displ) * k);

            free(D);
            displ += iteration_m[i];
        }

    }else{

        double *D = calc_D(X, Y, n, d, m);

        for(int j=0; j<m; j++) calc_part_knn(D, n, m, k, j, knn.nidx + j * k, knn.ndist + j * k);

        free(D);

    }

    return knn;
}

knnresult alt_kNN(  double *X,
                    double *Y,
                    int n,
                    int m,
                    int d,
                    int k)
{

    knnresult knn;
    knn.m = m;
    knn.k = k;
    knn.nidx = (int *)malloc(m * k * sizeof(int));
    knn.ndist = (double *)malloc(m * k * sizeof(double));

    for(int j=0; j<m; j++){

        double * ndist = calc_Dcol(X, Y, n, d, j);

        int * nidx = (int *)malloc(n * sizeof(int));
        for(int i=0; i<n; i++) nidx[i] = i;

        quickselect(nidx, ndist, 0, n-1, k);

        memcpy(knn.nidx + j * k, nidx, k * sizeof(int));
        memcpy(knn.ndist + j * k, ndist, k * sizeof(double));

        free(nidx);
        free(ndist);
    }

    return knn;
}

ext: h


text:
#ifndef OPENMC_TALLIES_FILTER_CELL_H
#define OPENMC_TALLIES_FILTER_CELL_H

#include <cstdint>
#include <unordered_map>

#include <gsl/gsl>

#include "openmc/tallies/filter.h"
#include "openmc/vector.h"

namespace openmc {

//==============================================================================
//! Specifies which geometric cells tally events reside in.
//==============================================================================

class CellFilter : public Filter
{
public:
  //----------------------------------------------------------------------------
  // Constructors, destructors

  ~CellFilter() = default;

  //----------------------------------------------------------------------------
  // Methods

  std::string type() const override {return "cell";}

  void from_xml(pugi::xml_node node) override;

  void get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match)
  const override;

  void to_statepoint(hid_t filter_group) const override;

  std::string text_label(int bin) const override;

  //----------------------------------------------------------------------------
  // Accessors

  const vector<int32_t>& cells() const { return cells_; }

  void set_cells(gsl::span<int32_t> cells);

protected:
  //----------------------------------------------------------------------------
  // Data members

  //! The indices of the cells binned by this filter.
  vector<int32_t> cells_;

  //! A map from cell indices to filter bin indices.
  std::unordered_map<int32_t, int> map_;
};

} // namespace openmc
#endif // OPENMC_TALLIES_FILTER_CELL_H

max_stars_repo_name: cjwyett/openmc

max_issues_repo_name: cjwyett/openmc

max_forks_repo_name: cjwyett/openmc


text:
#ifndef _SPC_SPC_H

#define _SPC_SPC_H

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

#include <gsl/gsl_nan.h>
#include <gsl/gsl_sys.h>

#include "aXe_grism.h"
#include "aXe_utils.h"
#include "fitsio.h"

#define SPCTOL 1.0e-5

// interpolation type used for the
// response function
#define RESP_FUNC_INTERP_TYPE gsl_interp_linear

extern spectrum *
allocate_spectrum (const int numbin);

extern void
free_spectrum (spectrum * const spec);

extern void
fprintf_spectrum (FILE * output, const spectrum * const sp);

extern spectrum *
subtract_spectra (spectrum * a, spectrum * b);

extern int
get_ID_index_to_SPC (char filename[], char ID[]);

extern void
add_ID_index_to_SPC (char filename[], char ID[], int hdunum);

extern int
add_ID_to_SPC (char filename[], int N, char ID[]);

extern int
find_ID_in_SPC (char filename[], char ID[]);

extern void
create_SPC (char filename[], int overwrite);

extern int
get_SPC_colnum (fitsfile * input, char colname[]);

extern void
add_spectra_to_SPC (char filename[], spectrum * obj_spec,
		    spectrum * bck_spec, spectrum * sobj_spec,
		    int aperID, int beamID);

extern void
add_data_to_SPC (spectrum * spec, char countcolname[],
		 char errorcolname[], char weightcolname[], char ID[],
		 char filename[], int hdunum, long N);

extern spectrum *
trim_spectrum (spectrum * spc);

extern spectrum *
empty_counts_spectrum_copy (spectrum * a);

extern void
add_ID_to_SPC_opened (fitsfile *output, int N, char ID[]);

extern void
add_spectra_to_SPC_opened (fitsfile *input, spectrum * obj_spec,
			   spectrum * bck_spec, spectrum * sobj_spec,
			   int aperID, int beamID);

extern void
add_data_to_SPC_opened (spectrum * spec, char countcolname[],
			char errorcolname[], char weightcolname[],
			char ID[], fitsfile *input, long N);

extern fitsfile *
create_SPC_opened (char filename[], int overwrite);

#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <galpy_potentials.h>
// MovingObjectPotential
// 3 arguments: amp, t0, tf
void constrain_range3(double * d) {
  // Constrains index to be within interpolation range
  if (*d < 0) *d = 0.0;
  if (*d > 1) *d = 1.0;
}

void cosinerule(double R1, double R2, double Rdist, double * cosphi, double *sinphi) {
// R1 is the distane of the LMC from the MW
// R2 is the distance of the Sat from the MW
// Rdist is the estimated distance of the sat from the LMC
  *cosphi = (R1*R1 - R2*R2 - Rdist*Rdist)/(2*R2*Rdist);
  *sinphi = pow(1 - (*cosphi)*(*cosphi), 0.5);
}

double MovingObjectDissipativeRforce(double R,double z, double phi,
				   double t,
				   struct potentialArg * potentialArgs,
				   double vR,double vT,
				   double vz){
  double amp,t0,tf,d_ind, Rdist, Rorb, RF, TF;
  double x,y,vx,vy, obj_x,obj_y,obj_z,obj_vx,obj_vy,obj_vz;
  double costheta, sintheta, phi2;
  double vxd, vyd, vzd, vR1, vT1;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  constrain_range3(&d_ind);

  // convert to cartesian coodinates
  x= R*cos(phi);
  y= R*sin(phi);
  vx=  vR*cos(phi) - vT*sin(phi);
  vy=  vR*sin(phi) + vT*cos(phi);

  // Interpolate x, y, z, vx, vy, vz
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  obj_vx= gsl_spline_eval(*(potentialArgs->spline1d+3),d_ind,
			 *(potentialArgs->acc1d+3));
  obj_vy= gsl_spline_eval(*(potentialArgs->spline1d+4),d_ind,
			 *(potentialArgs->acc1d+4));
  obj_vz= gsl_spline_eval(*(potentialArgs->spline1d+5),d_ind,
			 *(potentialArgs->acc1d+5));

  // convert velocity to reference frame of LMC
  vxd = vx - obj_vx;
  vyd = vy - obj_vy;
  vzd = vz - obj_vz;

  // determine angles for transformation
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  Rorb = pow(obj_x*obj_x + obj_y*obj_y, 0.5);
  phi2 = atan2(y-obj_y, x-obj_x);
  cosinerule(Rorb,R,Rdist,&costheta,&sintheta);

  // convert velocity to cyclindrical coordinates in reference frame of LMC.
  vR1=    cos(phi2)*vxd    + sin(phi2)*vyd;
  vT1=   -sin(phi2)*vxd    + cos(phi2)*vyd;

  // Calculate Radial and transverse forces
  RF= calcRforce(Rdist,(z-obj_z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg,vR1,vT1,vzd);
  TF= calcPhiforce(Rdist,(z-obj_z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg,vR1,vT1,vzd);

  // Transform force back into the refernce frame of the MW
  return  -amp*(RF*costheta  +TF*sintheta);
}

double MovingObjectDissipativezforce(double R,double z,double phi,
				      double t,
				      struct potentialArg * potentialArgs,
				      double vR,double vT,
				      double vz){
  double amp,t0,tf,d_ind, Rdist, Rorb;
  double x,y,vx,vy, obj_x,obj_y,obj_z,obj_vx,obj_vy,obj_vz;
  double costheta, sintheta, phi2;
  double vxd, vyd, vzd, vR1, vT1;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  constrain_range3(&d_ind);

  // convert to cartesian coodinates
  x= R*cos(phi);
  y= R*sin(phi);
  vx=  vR*cos(phi) - vT*sin(phi);
  vy=  vR*sin(phi) + vT*cos(phi);

  // Interpolate x, y, z, vx, vy, vz
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  obj_vx= gsl_spline_eval(*(potentialArgs->spline1d+3),d_ind,
			 *(potentialArgs->acc1d+3));
  obj_vy= gsl_spline_eval(*(potentialArgs->spline1d+4),d_ind,
			 *(potentialArgs->acc1d+4));
  obj_vz= gsl_spline_eval(*(potentialArgs->spline1d+5),d_ind,
			 *(potentialArgs->acc1d+5));

  // convert velocity to reference frame of LMC
  vxd = vx - obj_vx;
  vyd = vy - obj_vy;
  vzd = vz - obj_vz;

  // determine angles for transformation
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  Rorb = pow(obj_x*obj_x + obj_y*obj_y, 0.5);
  phi2 = atan2(y-obj_y, x-obj_x);
  cosinerule(Rorb,R,Rdist,&costheta,&sintheta);

  // convert velocity to cyclindrical coordinates in reference frame of LMC.
  vR1=    cos(phi2)*vxd    + sin(phi2)*vyd;
  vT1=   -sin(phi2)*vxd    + cos(phi2)*vyd;

  // Calculate z force
  return amp * calczforce(Rdist,(z-obj_z),phi,t,potentialArgs->nwrapped,
			   potentialArgs->wrappedPotentialArg,vR1,vT1,vzd);
}

double MovingObjectDissipativephiforce(double R,double z,double phi,
					double t,
					struct potentialArg * potentialArgs,
					double vR,double vT,
					double vz){
  double amp,t0,tf,d_ind, Rdist, Rorb, RF, TF;
  double x,y,vx,vy, obj_x,obj_y,obj_z,obj_vx,obj_vy,obj_vz;
  double costheta, sintheta, phi2;
  double vxd, vyd, vzd, vR1, vT1;
  double * args= potentialArgs->args;
  //Get args
  amp= *args;
  t0= *(args+1);
  tf= *(args+2);
  d_ind= (t-t0)/(tf-t0);
  constrain_range3(&d_ind);

  // convert to cartesian coodinates
  x= R*cos(phi);
  y= R*sin(phi);
  vx=  vR*cos(phi) - vT*sin(phi);
  vy=  vR*sin(phi) + vT*cos(phi);

  // Interpolate x, y, z, vx, vy, vz
  obj_x= gsl_spline_eval(*potentialArgs->spline1d,d_ind,*potentialArgs->acc1d);
  obj_y= gsl_spline_eval(*(potentialArgs->spline1d+1),d_ind,
			 *(potentialArgs->acc1d+1));
  obj_z= gsl_spline_eval(*(potentialArgs->spline1d+2),d_ind,
			 *(potentialArgs->acc1d+2));
  obj_vx= gsl_spline_eval(*(potentialArgs->spline1d+3),d_ind,
			 *(potentialArgs->acc1d+3));
  obj_vy= gsl_spline_eval(*(potentialArgs->spline1d+4),d_ind,
			 *(potentialArgs->acc1d+4));
  obj_vz= gsl_spline_eval(*(potentialArgs->spline1d+5),d_ind,
			 *(potentialArgs->acc1d+5));

  // convert velocity to reference frame of LMC
  vxd = vx - obj_vx;
  vyd = vy - obj_vy;
  vzd = vz - obj_vz;

  // determine angles for transformation
  Rdist= pow(pow(x-obj_x, 2)+pow(y-obj_y, 2), 0.5);
  Rorb = pow(obj_x*obj_x + obj_y*obj_y, 0.5);
  phi2 = atan2(y-obj_y, x-obj_x);
  cosinerule(Rorb,R,Rdist,&costheta,&sintheta);

  // convert velocity to cyclindrical coordinates in reference frame of LMC.
  vR1=    cos(phi2)*vxd    + sin(phi2)*vyd;
  vT1=   -sin(phi2)*vxd    + cos(phi2)*vyd;

  // Calculate Radial and transverse forces
  RF= calcRforce(Rdist,(z-obj_z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg,vR1,vT1,vzd);
  TF= calcPhiforce(Rdist,(z-obj_z),phi,t,potentialArgs->nwrapped,
		 potentialArgs->wrappedPotentialArg,vR1,vT1,vzd);

  // Transform force back into the refernce frame of the MW
  return  -amp*(-RF*sintheta  +TF*costheta);

}

ext: c

max_stars_repo_name: radsouza/galpy

max_issues_repo_name: radsouza/galpy

max_forks_repo_name: radsouza/galpy


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

#ifndef __GSL_VECTOR_FLOAT_H__
#define __GSL_VECTOR_FLOAT_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_float.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  float *data;
  gsl_block_float *block;
  int owner;
}
gsl_vector_float;

typedef struct
{
  gsl_vector_float vector;
} _gsl_vector_float_view;

typedef _gsl_vector_float_view gsl_vector_float_view;

typedef struct
{
  gsl_vector_float vector;
} _gsl_vector_float_const_view;

typedef const _gsl_vector_float_const_view gsl_vector_float_const_view;

/* Allocation */

gsl_vector_float *gsl_vector_float_alloc (const size_t n);
gsl_vector_float *gsl_vector_float_calloc (const size_t n);

gsl_vector_float *gsl_vector_float_alloc_from_block (gsl_block_float * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

gsl_vector_float *gsl_vector_float_alloc_from_vector (gsl_vector_float * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

void gsl_vector_float_free (gsl_vector_float * v);

/* Views */

_gsl_vector_float_view
gsl_vector_float_view_array (float *v, size_t n);

_gsl_vector_float_view
gsl_vector_float_view_array_with_stride (float *base,
                                         size_t stride,
                                         size_t n);

_gsl_vector_float_const_view
gsl_vector_float_const_view_array (const float *v, size_t n);

_gsl_vector_float_const_view
gsl_vector_float_const_view_array_with_stride (const float *base,
                                               size_t stride,
                                               size_t n);

_gsl_vector_float_view
gsl_vector_float_subvector (gsl_vector_float *v,
                            size_t i,
                            size_t n);

_gsl_vector_float_view
gsl_vector_float_subvector_with_stride (gsl_vector_float *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

_gsl_vector_float_const_view
gsl_vector_float_const_subvector (const gsl_vector_float *v,
                                  size_t i,
                                  size_t n);

_gsl_vector_float_const_view
gsl_vector_float_const_subvector_with_stride (const gsl_vector_float *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

float gsl_vector_float_get (const gsl_vector_float * v, const size_t i);
void gsl_vector_float_set (gsl_vector_float * v, const size_t i, float x);

float *gsl_vector_float_ptr (gsl_vector_float * v, const size_t i);
const float *gsl_vector_float_const_ptr (const gsl_vector_float * v, const size_t i);

void gsl_vector_float_set_zero (gsl_vector_float * v);
void gsl_vector_float_set_all (gsl_vector_float * v, float x);
int gsl_vector_float_set_basis (gsl_vector_float * v, size_t i);

int gsl_vector_float_fread (FILE * stream, gsl_vector_float * v);
int gsl_vector_float_fwrite (FILE * stream, const gsl_vector_float * v);
int gsl_vector_float_fscanf (FILE * stream, gsl_vector_float * v);
int gsl_vector_float_fprintf (FILE * stream, const gsl_vector_float * v,
                              const char *format);

int gsl_vector_float_memcpy (gsl_vector_float * dest, const gsl_vector_float * src);

int gsl_vector_float_reverse (gsl_vector_float * v);

int gsl_vector_float_swap (gsl_vector_float * v, gsl_vector_float * w);
int gsl_vector_float_swap_elements (gsl_vector_float * v, const size_t i, const size_t j);

float gsl_vector_float_max (const gsl_vector_float * v);
float gsl_vector_float_min (const gsl_vector_float * v);
void gsl_vector_float_minmax (const gsl_vector_float * v, float * min_out, float * max_out);

size_t gsl_vector_float_max_index (const gsl_vector_float * v);
size_t gsl_vector_float_min_index (const gsl_vector_float * v);
void gsl_vector_float_minmax_index (const gsl_vector_float * v, size_t * imin, size_t * imax);

int gsl_vector_float_add (gsl_vector_float * a, const gsl_vector_float * b);
int gsl_vector_float_sub (gsl_vector_float * a, const gsl_vector_float * b);
int gsl_vector_float_mul (gsl_vector_float * a, const gsl_vector_float * b);
int gsl_vector_float_div (gsl_vector_float * a, const gsl_vector_float * b);
int gsl_vector_float_scale (gsl_vector_float * a, const double x);
int gsl_vector_float_add_constant (gsl_vector_float * a, const double x);

int gsl_vector_float_isnull (const gsl_vector_float * v);
int gsl_vector_float_ispos (const gsl_vector_float * v);
int gsl_vector_float_isneg (const gsl_vector_float * v);

#ifdef HAVE_INLINE

extern inline
float
gsl_vector_float_get (const gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

extern inline
void
gsl_vector_float_set (gsl_vector_float * v, const size_t i, float x)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

extern inline
float *
gsl_vector_float_ptr (gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (float *) (v->data + i * v->stride);
}

extern inline
const float *
gsl_vector_float_const_ptr (const gsl_vector_float * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (i >= v->size)
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const float *) (v->data + i * v->stride);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_FLOAT_H__ */

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


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

#undef REAL
#define COMPLEX

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex32_t
 *
 *  CORE_clarfy applies an elementary reflector, or Householder matrix, H,
 *  to a N-by-N hermitian matrix C, from both the left and the right.
 *
 *  H is represented in the form
 *
 *     H = I - tau * v * v'
 *
 *  where  tau  is a scalar and  v  is a vector.
 *
 *  If tau is zero, then H is taken to be the unit matrix.
 *
 *******************************************************************************
 *
 * @param[in] N
 *          The number of rows and columns of the matrix C.  N >= 0.
 *
 * @param[in,out] A
 *          COMPLEX*16 array, dimension (LDA, N)
 *          On entry, the Hermetian matrix A.
 *          On exit, A is overwritten by H * A * H'.
 *
 * @param[in] LDA
 *         The leading dimension of the array A.  LDA >= max(1,N).
 *
 * @param[in] V
 *          The vector V that contains the Householder reflectors.
 *
 * @param[in] TAU
 *          The value tau.
 *
 * @param[out] WORK
 *          Workspace.
 *
 ******************************************************************************/
void
CORE_clarfy(int N,
            PLASMA_Complex32_t *A, int LDA,
            const PLASMA_Complex32_t *V,
            const PLASMA_Complex32_t *TAU,
            PLASMA_Complex32_t *WORK)
{
    static PLASMA_Complex32_t zzero =  0.0;
    static PLASMA_Complex32_t zmone = -1.0;

    int j;
    PLASMA_Complex32_t dtmp;

    /* Compute dtmp = X'*V */
    /* X = AVtau */
    cblas_chemv(CblasColMajor, CblasLower,
                N, CBLAS_SADDR(*TAU), A, LDA,
                V, 1, CBLAS_SADDR(zzero), WORK, 1);

    /* cblas_cdotc_sub(N, WORK, 1, V, 1, &dtmp);*/
    dtmp = 0.;
    for (j = 0; j < N ; j++)
        dtmp = dtmp + conjf(WORK[j]) * V[j];

    /* Compute 1/2 X'*V*t = 1/2*dtmp*tau  */
    dtmp = -dtmp * 0.5 * (*TAU);

   /* Compute W=X-1/2VX'Vt = X - dtmp*V */
    cblas_caxpy(N, CBLAS_SADDR(dtmp),
                V, 1, WORK, 1);

    /*
     * Performs the symmetric rank 2 operation
     *    A := alpha*x*y' + alpha*y*x' + A
     */
    cblas_cher2(CblasColMajor, CblasLower, N,
                CBLAS_SADDR(zmone), WORK, 1,
                                    V,    1,
                                    A,    LDA);

    return;
}
#undef COMPLEX


text:
#include <stdio.h>
#include <stdlib.h>
#include <gsl/vector/gsl_vector.h>
#include <gsl/gsl_math.h>

#define BASE_DOUBLE
#include <gsl/templates_on.h>
#include <gsl/vector/minmax_source.c>
#include <gsl/templates_off.h>
#undef BASE_DOUBLE

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
/*
 * SPDX-License-Identifier:  BSD-3-Clause
 *
 * example_10-StructOfArrays-CellLinkedList-OuterLoop-SymmetricalLoadBalancing.c :
 *      Example of SPH Density Calculation using
 *      fast neighbor search the main density loop via
 *      Cell Linked List method, Struct of Arrays (SoA)
 *      data layout, OpenMP parallelization at the
 *      cell-pair level, SIMD directives in the kernel
 *      and in the inner-most loop. It also implements
 *      symmetrical load balancing by moving the parallelism
 *      from iterating over cells to iterate over
 *      unique pairs of cell (i.e. i<j) and recyling the
 *      intermediary calculations.
 *
 * (C) Copyright 2021 José Hugo Elsas
 * Author: José Hugo Elsas <jhelsas@gmail.com>
 *
 * Command Line Options:
 *   -runs  <int>   : Set the number of repetitions (runs) for
 *                      calculating the density. The value of
 *                      the density is based on the last
 *                      iteration.
 *                    Default value: 1
 *   -run_seed <int>: Flag to set an alternative seed use for
 *                      for the PRNG. Instead of feeding seed
 *                      to the PRNG directly, it feeds
 *                      seed + iteration, as to generate different
 *                      configurations for each iteration.
 *                    Default value: 0 - (possible 0/1)
 *   -seed     <int>: Set the seed to use for the SPH particles
 *                      uniform position generation in the box
 *                    Default value: 123123123
 *
 *   -N        <int>: Set the number of SPH particles to be used
 *                    Default value: 1e5 = 100,000
 *   -h      <float>: Set the value of the smoothing kernel
 *                      parameter h, which corresponds to half
 *                      of the support of the kernel.
 *                    Default value: 0.05
 *
 *   -Nx       <int>: Set the number of Cells in the X direction
 *                    Default value: 10
 *   -Ny       <int>: Set the number of Cells in the Y direction
 *                    Default value: 10
 *   -Nz       <int>: Set the number of Cells in the Z direction
 *                    Default value: 10
 *
 *   -Xmin   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Y direction for
 *                    the Cell Linked List box
 *                    Default value: 0.0
 *   -Ymin   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 0.0
 *
 *   -Xmax   <float>: Set the lower bound in the X direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Ymax   <float>: Set the lower bound in the Y direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 *   -Zmax   <float>: Set the lower bound in the Z direction for
 *                      the Cell Linked List box
 *                    Default value: 1.0
 */

#include <math.h>
#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <limits.h>
#include <unistd.h>
#include <stdbool.h>
#include <sys/time.h>
#include <inttypes.h>

#include <omp.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_heapsort.h>

#include "sph_data_types.h"
#include "sph_linked_list.h"
#include "sph_utils.h"

#ifndef M_PI
#define M_PI (3.14159265358979323846)
#endif

#define COMPUTE_BLOCKS 5

#define dbg false

int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr,int64_t *temp_hash, linkedListBox *box, SPHparticle *lsph, double *times);

int compute_density_3d_symmetrical_load_ballance(int N, double h, SPHparticle *lsph, linkedListBox *box);

int compute_density_3d_chunk_symmetrical(int64_t node_begin, int64_t node_end,
                                         int64_t nb_begin, int64_t nb_end,double h,
                                         double* restrict x, double* restrict y,
                                         double* restrict z, double* restrict nu,
                                         double* restrict rhoi, double* restrict rhoj);

//void quickSort_int64_t2(int64_t *arr, int64_t low, int64_t high);

void quickSort_int64_t_2(int64_t *arr, int64_t low, int64_t high);

void mergeSort(int64_t arr[], int64_t l, int64_t r);

double w_bspline_3d_constant(double h);

#pragma omp declare simd
double w_bspline_3d_simd(double q);

#define ip_swap(a,b) (b)^=((a)^=((b)^=(a)))

KHASH_MAP_INIT_INT64(2, int64_t)

void counting_sort(int64_t N, linkedListBox *box, SPHparticle *lsph,void *swap_arr,int64_t *temp_hash);

int main(int argc, char **argv){
  bool run_seed = false;       // By default the behavior is is to use the same seed
  int runs = 1,err;            // it only runs once
  long int seed = 123123123;   // The default seed is 123123123
  int64_t N = 100000;          // The default number of particles is N = 1e5 = 100,000
  double h=0.05;               // The default kernel smoothing length is h = 0.05
  linkedListBox *box;          // Uninitialized Box containing the cells for the cell linked list method
  SPHparticle *lsph;           // Uninitialized array of SPH particles

  box = (linkedListBox*)malloc(1*sizeof(linkedListBox)); // Create a box representing the entire 3d domain

  // allow for command line customization of the run
  arg_parse(argc,argv,&N,&h,&seed,&runs,&run_seed,box);  // Parse the command line options
                                                         // line arguments and override default values

  err = SPHparticle_SoA_malloc(N,&lsph);
  if(err)
    fprintf(stderr,"error in SPHparticle_SoA_malloc\n");

  void *swap_arr = malloc(N*sizeof(double));
  int64_t *temp_hash = (int64_t*)malloc(2*N*sizeof(int64_t));
  double times[runs*COMPUTE_BLOCKS];

  for(int run=0;run<runs;run+=1)
    main_loop(run,run_seed,N,h,seed,swap_arr,temp_hash,box,lsph,times);

  bool is_cll = true;
  const char *prefix = "ex11,cll,SoA,outer,simd,symmLB,quicker";
  print_time_stats(prefix,is_cll,N,h,seed,runs,lsph,box,times);
  print_sph_particles_density(prefix,is_cll,N,h,seed,runs,lsph,box);

  SPHparticleSOA_safe_free(N,&lsph);
  safe_free_box(box);
  free(swap_arr);
  free(temp_hash);

  return 0;
}

/*
 *  Function main_loop:
 *    Runs the main loop of the program, including the particle array generation,
 *    density calculation and the timings annotations.
 *
 *    Arguments:
 *       run <int>            : index (or value) or the present iteration
 *       run_seed <bool>      : boolean defining whether to use run index for seed or not
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       seed <long int>      : seed for GSL PRNG generator to generate particle positions
 *       box  <linkedListBox> : Box of linked list cells, encapsulating the 3d domain
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *       times <double>       : Array to store the computation timings to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 *       times <double>       : Times is updated by reference
 */

int main_loop(int run, bool run_seed, int64_t N, double h, long int seed,
              void *swap_arr, int64_t *temp_hash, linkedListBox *box, SPHparticle *lsph, double *times)
{
  int err;

  if(run_seed)
    err = gen_unif_rdn_pos_box(N,seed+run,box,lsph);
  else
    err = gen_unif_rdn_pos_box(N,seed,box,lsph);

  if(err)
    fprintf(stderr,"error in gen_unif_rdn_pos\n");

  // ------------------------------------------------------ //

  double t0,t1,t2,t3,t4,t5;

  t0 = omp_get_wtime();

  err = compute_hash_MC3D(N,lsph,box);                               // Compute Morton Z 3D hash based on the
  if(err)                                                            // cell index for each of the X, Y and Z
    fprintf(stderr,"error in compute_hash_MC3D\n");                  // directions, in which a given particle reside

  t1 = omp_get_wtime();

  counting_sort(N,box,lsph,swap_arr,temp_hash);

  t2 = omp_get_wtime();

  t3 = omp_get_wtime();

  t4 = omp_get_wtime();

  err = compute_density_3d_symmetrical_load_ballance(N,h,lsph,box);  // Compute the density of the particles based
  if(err)                                                            // on the cell linked list method for fast
    fprintf(stderr,"error in compute_density_3d_load_ballanced\n");  // neighbor search

  // --------------------------------------------------------------- //

  t5 = omp_get_wtime();

  kh_clear(0, box->hbegin);
  kh_clear(1, box->hend);

  times[COMPUTE_BLOCKS*run+0] = t1-t0;                               // Time for compute morton Z 3d hash
  times[COMPUTE_BLOCKS*run+1] = t2-t1;                               // Time for sorting the particles' hashes
  times[COMPUTE_BLOCKS*run+2] = t3-t2;                               // Time for reordering all other arrays accordingly
  times[COMPUTE_BLOCKS*run+3] = t4-t3;                               // Time for setting up the interval hash tables
  times[COMPUTE_BLOCKS*run+4] = t5-t4;                               // Time for computing the SPH particle densities

  return 0;
}

/*
 *  Function compute_density_3d_symmetrical_load_ballance:
 *    Computes the SPH density from the particles using cell linked list with
 *    vectorization at the compute_density_3d_chunk level, but the parallelization
 *    done at the level of the outer-most loop of the compute_density_3d_cll_outerOmp
 *    function, not at the chunk level.
 *
 *    The parallelization is done at the level of unique cell pair instead of cells,
 *    with the indexes for the cell pairs pre-computed before parallelization.
 *
 *    Arguments:
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       lsph <SPHparticle>   : Array (pointer) of SPH particles to be updated
 *    Returns:
 *       0                    : error code returned
 *       lsph <SPHparticle>   : SPH particle array is updated in the rho field by reference
 */
int compute_density_3d_symmetrical_load_ballance(int N, double h, SPHparticle *lsph, linkedListBox *box){
  int64_t *node_begin,*node_end,*nb_begin,*nb_end;                        // Define the arrays for cell boundaries
  int64_t max_cell_pair_count = 0;                                        // and the number of cell pairs
  const double kernel_constant = w_bspline_3d_constant(h);

  max_cell_pair_count = count_box_pairs(box);                             // compute the maximum number of cell pairs

  node_begin = (int64_t*)malloc(max_cell_pair_count*sizeof(int64_t));     // allocate space for node_begin
  node_end   = (int64_t*)malloc(max_cell_pair_count*sizeof(int64_t));     // allocate space for node_end
  nb_begin   = (int64_t*)malloc(max_cell_pair_count*sizeof(int64_t));     // allocate space for nb_begin
  nb_end     = (int64_t*)malloc(max_cell_pair_count*sizeof(int64_t));     // allocate space for nb_end

  max_cell_pair_count = setup_unique_box_pairs(box,                       // set the values for cell pairs
                                               node_begin,node_end,
                                               nb_begin,nb_end);

  memset(lsph->rho,(int)0,N*sizeof(double));                              // Pre-initialize the density to zero

                                                                          // Parallelism was moved
                                                                          // to the level of unique pairs
  #pragma omp parallel for schedule(dynamic,5) proc_bind(master)          // Execute in parallel
  for(size_t i=0;i<max_cell_pair_count;i+=1){                             // over the unique pairs' array
    double local_rhoi[node_end[i] - node_begin[i]];                       // partial density array for node indexs
    double local_rhoj[  nb_end[i] -   nb_begin[i]];                       // partial density array for nb   indexs

    memset(local_rhoi,(int)0,(node_end[i]-node_begin[i])*sizeof(double)); // initialize node partial density to zero
    memset(local_rhoj,(int)0,    (nb_end[i]-nb_begin[i])*sizeof(double)); // initialize nb partial density to zero

    compute_density_3d_chunk_symmetrical(node_begin[i],node_end[i],       // Compute the density contribution
                                         nb_begin[i],nb_end[i],h,         // from this particular cell pair
                                         lsph->x,lsph->y,lsph->z,         // for both node and nb partial density
                                         lsph->nu,local_rhoi,
                                         local_rhoj);

    // merging the results can result in race conditions, therefore needs to be serialized
    #pragma omp critical                                                  // this serializes this code section
    {

      for(size_t ii=node_begin[i];ii<node_end[i];ii+=1){                  // iterate over the node_ cell
        lsph->rho[ii] += kernel_constant*local_rhoi[ii-node_begin[i]];    // add the partial density contribution
      }

      if(node_begin[i] != nb_begin[i])                                    // if sender and receiver are different
        for(size_t ii=nb_begin[i];ii<nb_end[i];ii+=1){                    // iterate over the nb_ cell
          lsph->rho[ii] += kernel_constant*local_rhoj[ii-nb_begin[i]];    // add the partial density contribution
        }
    }
  }

  free(node_begin);
  free(node_end);
  free(nb_begin);
  free(nb_end);

  return 0;
}

/*
 *  Function compute_density_3d_chunk_symmetrical:
 *    Computes the SPH density contribution to both the node_ cell and the nb_ cell.
 *    Vectorization in the inner-most loop, but no parallelization.
 *    The density contribution is symmetrical.
 *
 *    Arguments:
 *       N <int>              : Number of SPH particles to be used in the run
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *       x       <double*>    : Array of particles' X positions
 *       y       <double*>    : Array of particles' Y positions
 *       z       <double*>    : Array of particles' Z positions
 *       nu      <double*>    : Array of particles' density weights (i.e. masses)
 *    Returns:
 *       0                    : error code returned
 *       rho       <double*>  : Array of particles' densities
 */
int compute_density_3d_chunk_symmetrical(int64_t node_begin, int64_t node_end,
                                         int64_t nb_begin, int64_t nb_end,double h,
                                         double* restrict x, double* restrict y,
                                         double* restrict z, double* restrict nu,
                                         double* restrict rhoi, double* restrict rhoj){
  const double inv_h = 1./h;

  for(int64_t ii=node_begin;ii<node_end;ii+=1){ // Iterate over the ii index of the chunk
    double xii = x[ii];                         // Load the X component of the ii particle position
    double yii = y[ii];                         // Load the Y component of the ii particle position
    double zii = z[ii];                         // Load the Z component of the ii particle position

    #pragma omp simd                            // Hint at the compiler to vectorize the inner most loop
    for(int64_t jj=nb_begin;jj<nb_end;jj+=1){   // Iterate over the each other particle in jj loop
      double q = 0.;                            // Initialize the distance

      double xij = xii-x[jj];                   // Load and subtract jj particle's X position component
      double yij = yii-y[jj];                   // Load and subtract jj particle's Y position component
      double zij = zii-z[jj];                   // Load and subtract jj particle's Z position component

      q += xij*xij;                             // Add the jj contribution to the ii distance in X
      q += yij*yij;                             // Add the jj contribution to the ii distance in Y
      q += zij*zij;                             // Add the jj contribution to the ii distance in Z

      q = sqrt(q)*inv_h;                        // Sqrt to compute the normalized distance, measured in h

      double wij = w_bspline_3d_simd(q);        // compute the smoothing kernel separately for re-use

      rhoi[ii-node_begin] += nu[jj]*wij;        // add the jj contribution to ii density
      rhoj[jj-nb_begin]   += nu[ii]*wij;        // add the ii contribution to jj density
    }
  }

  return 0;
}

/*
 *  Function w_bspline_3d_constant:
 *    Returns the 3d normalization constant for the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       h <double>           : Smoothing Length for the Smoothing Kernel w_bspline
 *    Returns:
 *       3d bspline normalization density <double>
 */
double w_bspline_3d_constant(double h){
  return 3./(2.*M_PI*h*h*h);  // 3d normalization value for the b-spline kernel
}

/*
 *  Function w_bspline_3d_simd:
 *    Returns the un-normalized value of the cubic b-spline SPH smoothing kernel
 *
 *    Arguments:
 *       q <double>           : Distance between particles normalized by the smoothing length h
 *    Returns:
 *       wq <double>          : Unnormalized value of the kernel
 *
 *    Observation:
 *       Why not else if(q<2.)?
 *       Because if you use "else if", the compiler refuses to vectorize,
 *       This results in a large slowdown, as of 2.5x slower for example_04
 */
#pragma omp declare simd
double w_bspline_3d_simd(double q){
  double wq=0;
  double wq1 = (0.6666666666666666 - q*q + 0.5*q*q*q);             // The first polynomial of the spline
  double wq2 = 0.16666666666666666*(2.-q)*(2.-q)*(2.-q);           // The second polynomial of the spline

  if(q<2.)                                                         // If the distance is below 2
    wq = wq2;                                                      // Use the 2nd polynomial for the spline

  if(q<1.)                                                         // If the distance is below 1
    wq = wq1;                                                      // Use the 1st polynomial for the spline

  return wq;                                                       // return which ever value corresponds to the distance
}

int cmp_int64_t(const void *p,const void *q){
  int64_t *data1,*data2;
  data1 = (int64_t*)p;
  data2 = (int64_t*)q;

  if(*data1 < *data2)         // data[0] is the hash value,
    return -1;
  else if(*data1 == *data2)   // in the unsorted array
    return 0;
  else
    return 1;
}

#define swap_loop(N,lsph,temp_swap,member,type) for(int64_t i=0;i<(N);i+=1)                            \
                                                  (temp_swap)[i] = (lsph)->member[(lsph)->hash[2*i+1]];\
                                                memcpy((lsph)->member,temp_swap,(N)*sizeof(type))

/*
 *  Function counting_sort:
 *
 *
 *    Arguments:
 *     N            <int64_t>    :
 *     box      <linkedListBox*> :
 *     lsph       <SPHparticle*> :
 *     swap_arr     <double*>    :
 *     temp_hash   <int64_t*>    :
 *    Returns:
 *     box      <linkedListBox*> :
 *
 *    Observation:
 *       Why not else if(q<2.)?
 *       Because if you use "else if", the compiler refuses to vectorize,
 *       This results in a large slowdown, as of 2.5x slower for example_04
 */
void counting_sort(int64_t N, linkedListBox *box, SPHparticle *lsph,void *swap_arr,int64_t *temp_hash){

  double t0,t1,t2,t3,t4,t5,t6;
  double t31,t32,t33,t34;

  t0 = omp_get_wtime();

  for(int64_t i=0;i<N;i+=1){
    int ret;
    khiter_t k = kh_put(0, box->hbegin, lsph->hash[2*i+0], &ret);
    if(ret==1)
      kh_value(box->hbegin, k) = 1;
    else if(ret==0){
      int64_t val = kh_value(box->hbegin, k);
      kh_value(box->hbegin, k) = val + 1;
    }
    else{
      printf("error counting_sort init");
      return;
    }
  }

  t1 = omp_get_wtime();

  unsigned int dict_size = kh_size(box->hbegin);
  int64_t hash[dict_size];
  int64_t prefix[dict_size];
  int idx = 0;

  for (khiter_t k = kh_begin(box->hbegin); k != kh_end(box->hbegin); ++k){
    if (kh_exist(box->hbegin, k)){
      hash[idx] = kh_key(box->hbegin,k);
      idx+=1;
    }
  }

  t2 = omp_get_wtime();

  qsort(hash,dict_size,sizeof(int64_t),cmp_int64_t);

  t3 = omp_get_wtime();

  prefix[0] = 0;
  for(int i=1;i<dict_size;i+=1){
    prefix[i] = prefix[i-1];

    khiter_t k = kh_get(0, box->hbegin, hash[i-1]);
    prefix[i] += kh_value(box->hbegin, k);
  }

  t31 = omp_get_wtime();

  for(int i = 0; i< dict_size;i+=1){
    khiter_t kp = kh_get(0, box->hbegin, hash[i]);
    kh_value(box->hbegin, kp) = prefix[i];
  }

  t32 = omp_get_wtime();

  for (khiter_t k = kh_begin(box->hbegin); k != kh_end(box->hbegin); ++k){
    if (kh_exist(box->hbegin, k)){
      int ret;
      int64_t hash = kh_key(box->hbegin,k);
      khiter_t ke = kh_put(1, box->hend, hash, &ret);
      kh_value(box->hend, ke) = kh_value(box->hbegin, k);
    }
  }

  t33 = omp_get_wtime();

  for(int64_t i=0;i<N;i+=1){
    khiter_t kp = kh_get(1, box->hend, lsph->hash[2*i]);
    if (kh_exist(box->hend, kp)){
      lsph->hash[2*i+1] = kh_value(box->hend, kp);
      kh_value(box->hend, kp) += 1;
    }
    else{
      printf("there is an issue\n");
    }
  }

  t34 = omp_get_wtime();

  t4 = omp_get_wtime();

  for(int64_t i=0;i<N;i+=1){
    int64_t idx = lsph->hash[2*i+1];
    temp_hash[2*idx+0] = lsph->hash[2*i+0];
    temp_hash[2*idx+1] = i;
  }

  for(int64_t i=0;i<N;i+=1){
    lsph->hash[2*i+0] = temp_hash[2*i+0];
    lsph->hash[2*i+1] = temp_hash[2*i+1];
  }

  t4 = omp_get_wtime();

  int64_t *int64_temp_swap = (int64_t *)swap_arr;
  swap_loop(N,lsph,int64_temp_swap,id ,int64_t);
  double *double_temp_swap = (double *)swap_arr;
  swap_loop(N,lsph,double_temp_swap,nu ,double);
  swap_loop(N,lsph,double_temp_swap,rho,double);
  swap_loop(N,lsph,double_temp_swap,x  ,double);
  swap_loop(N,lsph,double_temp_swap,y  ,double);
  swap_loop(N,lsph,double_temp_swap,z  ,double);
  swap_loop(N,lsph,double_temp_swap,ux ,double);
  swap_loop(N,lsph,double_temp_swap,uy ,double);
  swap_loop(N,lsph,double_temp_swap,uz ,double);
  swap_loop(N,lsph,double_temp_swap,Fx ,double);
  swap_loop(N,lsph,double_temp_swap,Fy ,double);
  swap_loop(N,lsph,double_temp_swap,Fz ,double);

  t5 = omp_get_wtime();

  for(int64_t i=0;i<N;i+=1)
    lsph->hash[i] = lsph->hash[2*i];

  t6 = omp_get_wtime();

  if(dbg){
    printf("1: %lf \n",t1-t0);
    printf("2: %lf \n",t2-t1);
    printf("3: %lf \n",t3-t2);
    printf("4: %lf \n",t4-t3);
    printf("  31: %lf\n",t31-t3);
    printf("  32: %lf\n",t32-t31);
    printf("  33: %lf\n",t33-t32);
    printf("  34: %lf\n",t34-t33);
    printf("  35: %lf\n",t4-t34);
    printf("5: %lf \n",t5-t4);
    printf("6: %lf \n",t6-t5);
  }

  return;
}

max_stars_repo_name: jhelsas/sphalerite

max_issues_repo_name: jhelsas/sphalerite

max_forks_repo_name: jhelsas/sphalerite


text:
////////////////////////////////////////////////////////////
//
// Copyright (c) 2018 Jan Filipowicz, Filip Turobos
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of this software and associated documentation files (the "Software"), to deal
// in the Software without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of the Software, and to permit persons to whom the Software is
// furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included in all
// copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
// SOFTWARE.
//
////////////////////////////////////////////////////////////

#ifndef GENETIC_ALGORITHM_LIBRARY_ELITIST_SELECTION_H
#define GENETIC_ALGORITHM_LIBRARY_ELITIST_SELECTION_H

#include <algorithm>
#include <cstddef>
#include <functional>
#include <vector>
#include <gsl/gsl_assert>

template<class Compare = std::less<>>
class elitist_selection {
public:
	explicit elitist_selection(const Compare& comp = Compare()) noexcept(noexcept(Compare(comp)));
	template<class Specimen>
	void operator()(std::vector<Specimen>& specimens, std::size_t n) const;
private:
	Compare comparator;
};

template<class Compare>
inline elitist_selection<Compare>::elitist_selection(const Compare& comp) noexcept(noexcept(Compare(comp)))
	: comparator(comp) {}

template<class Compare>
template<class Specimen>
inline void elitist_selection<Compare>::operator()(std::vector<Specimen>& specimens, std::size_t n) const {
	Expects(specimens.size() >= n);
	std::nth_element(specimens.begin(), specimens.begin() + n, specimens.end(), [this](const Specimen& lhs, const Specimen& rhs) {
		return comparator(rhs.rating(), lhs.rating());
	});
	specimens.resize(n);
}

#endif

max_stars_repo_name: SirEmentaler/Eugenics-Wars

max_issues_repo_name: SirEmentaler/Eugenics-Wars

max_forks_repo_name: SirEmentaler/Eugenics-Wars


text:
#include <cblas.h>
#include "util.h"
#include "conv.h"

extern float *scratch;

convLayer convLayer_init(int Sm, int Sn, int p)
{
     convLayer cl; CONVL_INIT(cl);
     cl.Sm=Sm; cl.Sn=Sn; cl.p=p;
     return cl;
}

void convLayer_free(convLayer *cl)
{
     ftens_free(&cl->W);   ftens_free(&cl->b);
     ftens_free(&cl->dW);  ftens_free(&cl->db);
     ftens_free(&cl->out); ftens_free(&cl->in);
}

void convLayer_print_shape(convLayer *cl)
{
     printf("conv: D=%d M=%d N=%d L=%d Sm=%d Sn=%d p=%d\n",
            cl->D, cl->M, cl->N, cl->L, cl->Sm, cl->Sn, cl->p);
}

void convLayer_set(ftens *W, convLayer *cl)
{
     int D=W->D, M=W->M, N=W->N, L=W->L;
     ftens_free(&cl->W);
     cl->D=D; cl->M=M; cl->N=N; cl->L=L;
     cl->W = ftens_copy(W);
}

void convLayer_copy_input(ftens *t, convLayer *cl)
{
     if (!cl->in.data)
          cl->in=ftens_init(t->D, t->M, t->N, t->L);
     memcpy(cl->in.data, t->data, t->bytes);
}

ftens convLayer_pad_input(ftens *t, float *scr,
                          int *M, int *N, int p)
{
     ftens tp; const int D=t->D, L=t->L;
     *M=PAD(*M, p); *N=PAD(*N, p);
     if (!scratch) tp=ftens_copy_pad(t, p);
     else {
          tp = ftens_from_ptr(D, *M, *N, L, scr);
          ftens_pad(t, &tp, p);
          scr += (*M)*(*N)*L*D;
     }

     return tp;
}

void convLayer_forward(ftens *t, convLayer *cl, int save)
{
     float *scr = scratch; ftens tp, tmp;
     int D=t->D,  Ms=t->M, Ns=t->N, Ls=t->L;
     int F=cl->D, W=cl->M, H=cl->N, L=cl->L;
     int p=cl->p, Sy=cl->Sm, Sx=cl->Sn;
     ASSERT(t->L == cl->L, "err: conv shape\n");

     if (save)      convLayer_copy_input(t, cl);
     if (p)    tp = convLayer_pad_input(t, scr, &Ms, &Ns, p);

     // lower
     const int Md = OUT_LEN(Ms, H, Sy);
     const int Nd = OUT_LEN(Ns, W, Sx);
     const int Ld =  W*H*L;
     if (!scratch) tmp=ftens_init(D, Md, Nd, Ld);
     else          tmp=ftens_from_ptr(D, Md, Nd, Ld, scr);

     ftens_lower(p ? &tp : t, &tmp, W, H, Sx, Sy);

     // mat mul
     if (!cl->out.data) cl->out=ftens_init(D, Md, Nd, F);
     int M=Md*Nd, N=F, K=cl->W.MNL;
     cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
                 M, N, K, 1, tmp.data, K, cl->W.data, K,
                 0, cl->out.data, N);

     if (!scratch)      ftens_free(&tmp);
     if (!scratch && p) ftens_free(&tp);
}

void convLayer_backward(ftens *dout, convLayer *cl)
{
     exit(-2);
}

void convLayer_update(convLayer *cl)
{
     exit(-3);
}

max_stars_repo_name: pppoe/espresso

max_issues_repo_name: pppoe/espresso

max_forks_repo_name: pppoe/espresso


text:
/*  Copyright (c) 2014, J.M. Hernandez-Lobato, M.W. Hoffman, Z. Ghahramani
 This function is from the code for the paper
 Hernández-Lobato J. M., Hoffman M. W. and Ghahramani Z.
 Predictive Entropy Search for Efficient Global Optimization of Black-box
 Functions, In NIPS, 2014.
 https://bitbucket.org/jmh233/codepesnips2014
 */
#include "mex.h"
#include <gsl/gsl_math.h>
#include <gsl/gsl_linalg.h>

/* Function that computes the inverse of a matrix using its Cholesky decomposition */

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {

	#define L_matlab prhs[0]
	#define ret_matlab plhs[0]

	int i, j, n, m;
	double *L, *ret;
	gsl_matrix *result;

	n = mxGetN(L_matlab);
	m = mxGetM(L_matlab);

	/* We ask for memory to store the matrix in the gnu library format */

	result = gsl_matrix_alloc(n, m);

	/* We copy the matrix from the matlab format to the gnu library format */

	L = mxGetPr(L_matlab);
	for (i = 0 ; i < n ; i++)
		for (j = 0 ; j < m ; j++)
			gsl_matrix_set(result, i, j, L[ i + n * j ]);

	gsl_linalg_cholesky_decomp(result);

	/* We obtain the inverse of the matrix */

	gsl_linalg_cholesky_invert(result);

	/* We ask for memory to return the solution */

	ret_matlab = mxCreateDoubleMatrix(n, m, mxREAL);

	/* We copy the solution from the gnu library representation to the matlab representation */

	ret = mxGetPr(ret_matlab);
	for (i = 0 ; i < n ; i++)
		for (j = 0 ; j < m ; j++)
			ret[ i + j * n ] = gsl_matrix_get(result, i, j);

	/* We are done */

	gsl_matrix_free(result);

	return;
}

max_stars_repo_name: pubino/Max-value-Entropy-Search


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

#ifndef __GSL_STATISTICS_INT_H__
#define __GSL_STATISTICS_INT_H__

#include <stddef.h>
#include <gsl/gsl_types.h>

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

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_int_mean (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_variance (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_sd (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_variance_with_fixed_mean (const int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_int_sd_with_fixed_mean (const int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_int_absdev (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_skew (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_kurtosis (const int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_int_lag1_autocorrelation (const int data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_int_covariance (const int data1[], const size_t stride1,const int data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_int_variance_m (const int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_int_sd_m (const int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_int_absdev_m (const int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_int_skew_m_sd (const int data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_int_kurtosis_m_sd (const int data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_int_lag1_autocorrelation_m (const int data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_int_covariance_m (const int data1[], const size_t stride1,const int data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

GSL_EXPORT double gsl_stats_int_pvariance (const int data1[], const size_t stride1, const size_t n1, const int data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_int_ttest (const int data1[], const size_t stride1, const size_t n1, const int data2[], const size_t stride2, const size_t n2);

GSL_EXPORT int gsl_stats_int_max (const int data[], const size_t stride, const size_t n);
GSL_EXPORT int gsl_stats_int_min (const int data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_int_minmax (int * min, int * max, const int data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_int_max_index (const int data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_int_min_index (const int data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_int_minmax_index (size_t * min_index, size_t * max_index, const int data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_int_median_from_sorted_data (const int sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_int_quantile_from_sorted_data (const int sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_INT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#pragma once

#include <gsl/gsl>

#include <vector>
#include <memory>

struct AVCodecContext;
struct AVPacket;
struct AVFrame;

namespace video_streamer
{

class FrameEncoder
{
public:
    explicit FrameEncoder(gsl::not_null<AVCodecContext*> codecContext);

    std::vector<unsigned char> encode(gsl::not_null<const AVFrame*> frame);
private:
    gsl::not_null<AVCodecContext*> m_codecContext;
    gsl::not_null<AVPacket*> m_packet;
};

std::shared_ptr<FrameEncoder> createFrameEncoder(
    const std::string& encoderName, gsl::not_null<const AVCodecContext*> codecContext);

} // video_streamer

max_stars_repo_name: nicledomaS/VideoStream

max_issues_repo_name: nicledomaS/VideoStream

max_forks_repo_name: nicledomaS/VideoStream


text:
#ifndef BASIS_H
#define BASIS_H

#include <algorithm>
#include <cstddef>
#include <cmath>
#include <map>
#include <memory>
#include <mutex>
#include <vector>

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_wavelet.h>

#include "constants.h"
#include "vector.h"

using std::array;
using std::shared_ptr;
using std::unique_ptr;
using std::vector;
using std::map;
using std::pair;

using BasisMap = Vector[]; // length
using ConvMap = Vector[]; // length * length

void projForward(gsl_vector *f);
void projBackward(gsl_vector *u);

void projForward(gsl_vector *f, gsl_vector *u);
void projBackward(gsl_vector *u, gsl_vector *f);

const Vector& getBasis(size_t k);
const Vector& getConvoluted(size_t p, size_t f);

double convolutedBasis(size_t p, size_t f, size_t mu);

void normalize(gsl_vector *f);

extern unique_ptr<BasisMap> storedBases;
extern unique_ptr<ConvMap> storedConvs;

#endif // WAVELETS_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse