text:
/*
   Authors:
   G. Jungman
   J. Scott (james.scott@lexifi.com)

 */
/* Implement Niederreiter base 2 generator.
 * See:
 *   Bratley, Fox, Niederreiter, ACM Trans. Model. Comp. Sim. 2, 195 (1992)
 */
#include <config.h>
#include <gsl/gsl_qrng.h>
#include <gsl/gsl_rng.h>
#include <assert.h>

#define NIED2_BIT_COUNT 30
#define NIED2_NBITS (NIED2_BIT_COUNT+1)

/* Z_2 field operations */
#define NIED2_ADD(x,y) (((x)+(y))%2)
#define NIED2_MUL(x,y) (((x)*(y))%2)
#define NIED2_SUB(x,y) NIED2_ADD((x),(y))

void get_primitive_polynomials(int dimension, int *degree_table, int *primitive_polynomials);

static size_t nied2_state_size(unsigned int dimension);
static int nied2_init(void * state, unsigned int dimension);
static int nied2_get(void * state, unsigned int dimension, double * v);

static const gsl_qrng_type nied2_type =
{
  "niederreiter-base-2",
  0,
  nied2_state_size,
  nied2_init,
  nied2_get
};

const gsl_qrng_type * gsl_qrng_niederreiter_2 = &nied2_type;

typedef struct
{
  unsigned int sequence_count;
  int *cj[NIED2_NBITS];
  int *nextq;
} nied2_state_t;

static size_t nied2_state_size(unsigned int dimension)
{
  return sizeof(nied2_state_t) +             /* the struct */
    sizeof(int) * dimension * NIED2_NBITS +  /* cj */
    sizeof(int) * dimension;                 /* nextq */
}

/* Multiply polynomials over Z_2.
 * Notice use of a temporary vector,
 * side-stepping aliasing issues when
 * one of inputs is the same as the output
 * [especially important in the original fortran version, I guess].
 */
static void poly_multiply(
  const int pa[], int pa_degree,
  const int pb[], int pb_degree,
  int pc[], int  * pc_degree, int max_degree
  )
{
  int j, k;
  int *pt;
  int pt_degree;

  pt = (int *) malloc(sizeof(int) * (max_degree+1));
  pt_degree = pa_degree + pb_degree;

  assert(pt_degree <= max_degree);

  for(k=0; k<=pt_degree; k++) {
    int term = 0;
    for(j=0; j<=k; j++) {
      const int conv_term = NIED2_MUL(pa[k-j], pb[j]);
      term = NIED2_ADD(term, conv_term);
    }
    pt[k] = term;
  }

  for(k=0; k<=pt_degree; k++) {
    pc[k] = pt[k];
  }
  for(k=pt_degree+1; k<=max_degree; k++) {
    pc[k] = 0;
  }

  *pc_degree = pt_degree;
  free(pt);
}

/* Calculate the values of the constants V(J,R) as
 * described in BFN section 3.3.
 *
 *   px = appropriate irreducible polynomial for current dimension
 *   pb = polynomial defined in section 2.3 of BFN.
 * pb is modified
 */
static void calculate_v(
  const  int px[],  int px_degree,
   int pb[],  int * pb_degree,
  int v[],  int maxv, int max_degree, gsl_rng *rng
  )
{
  const int nonzero_element = 1;    /* nonzero element of Z_2  */

  /* The polynomial ph is px**(J-1), which is the value of B on arrival.
   * In section 3.3, the values of Hi are defined with a minus sign:
   * don't forget this if you use them later !
   */
  int *ph;
  /* int ph_degree = *pb_degree; */
  int bigm = *pb_degree;      /* m from section 3.3 */
  int m;                      /* m from section 2.3 */
  int r, k, kj;

  ph = (int *) malloc(sizeof(int) * (max_degree+1));

  for(k=0; k<=max_degree; k++) {
    ph[k] = pb[k];
  }

  /* Now multiply B by PX so B becomes PX**J.
   * In section 2.3, the values of Bi are defined with a minus sign :
   * don't forget this if you use them later !
   */
   poly_multiply(px, px_degree, pb, *pb_degree, pb, pb_degree, max_degree);
   m = *pb_degree;

  /* Now choose a value of Kj as defined in section 3.3.
   * We must have 0 <= Kj < E*J = M.
   * The limit condition on Kj does not seem very relevant
   * in this program.
   */
  /* Quoting from BFN: "Our program currently sets each K_q
   * equal to eq. This has the effect of setting all unrestricted
   * values of v to 1."
   * Actually, it sets them to the arbitrary chosen value.
   * Whatever.
   */
   kj = gsl_rng_uniform_int(rng, bigm+1);

  /* Now choose values of V in accordance with
   * the conditions in section 3.3.
   */
  for(r=0; r<kj; r++) {
    v[r] = 0;
  }
  v[kj] = 1;

  if(kj >= bigm) {
    for(r=kj+1; r<m; r++) {
      v[r] =  gsl_rng_uniform_int(rng, 2); /* arbitrary element */
    }
  }
  else {
    /* This block is never reached. */

    int term = NIED2_SUB(0, ph[kj]);

    for(r=kj+1; r<bigm; r++) {
      v[r] =  gsl_rng_uniform_int(rng, 2); /* arbitrary element */

      /* Check the condition of section 3.3,
       * remembering that the H's have the opposite sign.  [????????]
       */
      term = NIED2_SUB(term, NIED2_MUL(ph[r], v[r]));
    }

    /* Now v[bigm] != term. */
    v[bigm] = NIED2_ADD(nonzero_element, term);

    for(r=bigm+1; r<m; r++) {
      v[r] = gsl_rng_uniform_int(rng, 2); /* arbitrary element */
    }
  }

  /* Calculate the remaining V's using the recursion of section 2.3,
   * remembering that the B's have the opposite sign.
   */
  for(r=0; r<=maxv-m; r++) {
    int term = 0;
    for(k=0; k<m; k++) {
      term = NIED2_SUB(term, NIED2_MUL(pb[k], v[r+k]));
    }
    v[r+m] = term;
  }
  free(ph);
}

static int calculate_cj(nied2_state_t * ns,  int dimension)
{
  int ci[NIED2_NBITS][NIED2_NBITS];
  int *v;
  int r;
  int i_dim;
  int max_degree, maxv;
  int *pb, *px;

  gsl_rng *rng = gsl_rng_alloc(gsl_rng_default);

  int *primitive_polynomials;
  int *poly_degree;

  poly_degree = (int *) malloc(sizeof(int) * (1+dimension));
  primitive_polynomials = (int *) malloc(sizeof(int) * (1+dimension));

  if (poly_degree==0 || primitive_polynomials==0)
    GSL_ERROR_NULL ("allocation of degree table failed for niederreiter init", GSL_ENOMEM);

  /* Generate the primitive polynomials.  */
  get_primitive_polynomials(dimension, poly_degree+1, primitive_polynomials+1);

  poly_degree[0]=0;  poly_degree[1]=1;
  primitive_polynomials[0]=1;
  primitive_polynomials[1]=2;

  max_degree=NIED2_NBITS+poly_degree[dimension-1];
  maxv = max_degree;

  v = (int *) malloc(sizeof(int) * (1+maxv));
  pb = (int *) malloc(sizeof(int) * (1+max_degree));
  px = (int *) malloc(sizeof(int) * (1+max_degree));

  if (v==0 || pb==0 || px==0)
    GSL_ERROR_NULL ("allocation of degree table failed for niederreiter init", GSL_ENOMEM);

  for(i_dim=0; i_dim<dimension; i_dim++) {

    const int poly_index = i_dim + 1;
    int j, k;

    /* Niederreiter (page 56, after equation (7), defines two
     * variables Q and U.  We do not need Q explicitly, but we
     * do need U.
     */
    int u = 0;

    /* For each dimension, we need to calculate powers of an
     * appropriate irreducible polynomial, see Niederreiter
     * page 65, just below equation (19).
     * Copy the appropriate irreducible polynomial into PX,
     * and its degree into E.  Set polynomial B = PX ** 0 = 1.
     * M is the degree of B.  Subsequently B will hold higher
     * powers of PX.
     */
    int px_degree;
    int pb_degree;
    px_degree = poly_degree[poly_index];
    pb_degree = 0;

    for(k=0; k<=px_degree; k++) {
      px[k] = (primitive_polynomials[poly_index] >> k) & 1;
      pb[k] = 0;
    }

    for (;k<max_degree+1;k++) {
      px[k] = 0;
      pb[k] = 0;
    }

    pb[0] = 1;

    for(j=0; j<NIED2_NBITS; j++) {
      /* If U = 0, we need to set B to the next power of PX
       * and recalculate V.
       */
      if(u == 0) calculate_v(px, px_degree, pb, &pb_degree, v, maxv, max_degree, rng);

      /* Now C is obtained from V.  Niederreiter
       * obtains A from V (page 65, near the bottom), and then gets
       * C from A (page 56, equation (7)).  However this can be done
       * in one step.  Here CI(J,R) corresponds to
       * Niederreiter's C(I,J,R).
       */
      for(r=0; r<NIED2_NBITS; r++) {
        ci[r][j] = v[r+u];
      }

      /* Advance Niederreiter's state variables. */
      ++u;
      if(u == px_degree) u = 0;
    }

    /* The array CI now holds the values of C(I,J,R) for this value
     * of I.  We pack them into array CJ so that CJ(I,R) holds all
     * the values of C(I,J,R) for J from 1 to NBITS.
     */
    for(r=0; r<NIED2_NBITS; r++) {
      int term = 0;
      for(j=0; j<NIED2_NBITS; j++) {
        term = 2*term + ci[r][j];
      }
      ns->cj[r][i_dim] = term;
    }
  }
  free(primitive_polynomials);
  free(poly_degree);
  free(v);
  free(pb);
  free(px);
  gsl_rng_free(rng);
  return GSL_SUCCESS;
}

static int nied2_init(void * state, unsigned int dimension)
{
  nied2_state_t * n_state = (nied2_state_t *) state;
  unsigned int i_dim, i_bits;
  int ret;

  if(dimension < 1) return GSL_EINVAL;

  for (i_bits=0; i_bits<NIED2_NBITS; i_bits++)
    n_state->cj[i_bits] =
      (int *) ((char *) n_state +
	       sizeof(nied2_state_t) +
	       sizeof(int) * dimension * i_bits);

  n_state->nextq =
    (int *) ((char *) n_state +
	     sizeof(nied2_state_t) +
	     sizeof(int) * dimension * NIED2_NBITS);

  if ((ret = calculate_cj(n_state, dimension)) != GSL_SUCCESS)
    return ret;

  /* skip the '0th' state, which is 0 */
  for(i_dim=0; i_dim<dimension; i_dim++)
    n_state->nextq[i_dim] = n_state->cj[0][i_dim];

  n_state->sequence_count = 1;

  return GSL_SUCCESS;
}

static int nied2_get(void * state, unsigned int dimension, double * v)
{
  static const double recip = 1.0/(double)(1U << NIED2_NBITS); /* 2^(-nbits) */
  nied2_state_t * n_state = (nied2_state_t *) state;
  int r;
  int c;
  unsigned int i_dim;

  /* Load the result from the saved state. */
  for(i_dim=0; i_dim<dimension; i_dim++) {
    v[i_dim] = n_state->nextq[i_dim] * recip;
  }

  /* Find the position of the least-significant zero in sequence_count.
   * This is the bit that changes in the Gray-code representation as
   * the count is advanced.
   */
  r = 0;
  c = n_state->sequence_count;
  while(1) {
    if((c % 2) == 1) {
      ++r;
      c /= 2;
    }
    else break;
  }

  if(r >= NIED2_NBITS) return GSL_EFAILED; /* FIXME: better error code here */

  /* Calculate the next state. */
  for(i_dim=0; i_dim<dimension; i_dim++) {
    n_state->nextq[i_dim] ^= n_state->cj[r][i_dim];
  }

  n_state->sequence_count++;

  return GSL_SUCCESS;
}


text:
#ifndef TTT_AUDIOCACHE_H
#define TTT_AUDIOCACHE_H

#include <SDL_mixer.h>

#include <gsl/util>

#include <vector>
#include <memory>
#include <string_view>

namespace ttt
{
	// A class for caching sounds and music.
	class AudioCache final
	{
	public:
		using SizeType = gsl::index;

		Mix_Music *loadMusic(std::string_view path);
		Mix_Chunk *loadChunk(std::string_view path);

		void unloadMusic(SizeType index) noexcept;
		void unloadChunk(SizeType index) noexcept;

		void eraseMusic(SizeType index) noexcept;
		void eraseChunk(SizeType index) noexcept;

		void clearMusic() noexcept;
		void clearChunks() noexcept;
		void clear() noexcept;

		[[nodiscard]] Mix_Music *getMusic(SizeType index) const noexcept;
		[[nodiscard]] Mix_Chunk *getChunk(SizeType index) const noexcept;

	private:
		std::vector<std::unique_ptr<Mix_Music, decltype(&Mix_FreeMusic)>> music;
		std::vector<std::unique_ptr<Mix_Chunk, decltype(&Mix_FreeChunk)>> chunks;
	};

	inline void AudioCache::unloadMusic(SizeType index) noexcept { gsl::at(music, index).reset(); }
	inline void AudioCache::unloadChunk(SizeType index) noexcept { gsl::at(chunks, index).reset(); }

	inline void AudioCache::eraseMusic(SizeType index) noexcept { music.erase(music.cbegin() + index); }
	inline void AudioCache::eraseChunk(SizeType index) noexcept { chunks.erase(chunks.cbegin() + index); }

	inline void AudioCache::clearMusic() noexcept { music.clear(); }
	inline void AudioCache::clearChunks() noexcept { chunks.clear(); }

	inline Mix_Music *AudioCache::getMusic(SizeType index) const noexcept { return gsl::at(music, index).get(); }
	inline Mix_Chunk *AudioCache::getChunk(SizeType index) const noexcept { return gsl::at(chunks, index).get(); }
}

#endif

ext: h

max_stars_repo_name: itsArtem/TicTacToe

max_issues_repo_name: itsArtem/TicTacToe

max_forks_repo_name: itsArtem/TicTacToe


text:
/**
 *
 * @file qwrapper_dsygst.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:59 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

#undef COMPLEX
#define REAL

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dsygst(Quark *quark, Quark_Task_Flags *task_flags,
                       int itype, PLASMA_enum uplo, int n,
                       double *A, int lda,
                       double *B, int ldb,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    QUARK_Insert_Task(quark, CORE_dsygst_quark, task_flags,
        sizeof(int),                        &itype,      VALUE,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(double)*lda*n,    A,             INOUT,
        sizeof(int),                        &lda,       VALUE,
#ifdef COMPLEX
        sizeof(double)*ldb*n,    B,             INOUT,
#else
        sizeof(double)*ldb*n,    B,             INPUT,
#endif
        sizeof(int),                        &ldb,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dsygst_quark = PCORE_dsygst_quark
#define CORE_dsygst_quark PCORE_dsygst_quark
#endif
void CORE_dsygst_quark(Quark *quark)
{
    int itype;
    PLASMA_enum uplo;
    int n;
    double *A;
    int lda;
    double *B;
    int ldb;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_10(quark, itype, uplo, n, A, lda, B, ldb, sequence, request, iinfo);
    info = LAPACKE_dsygst_work(
        LAPACK_COL_MAJOR,
        itype,
        lapack_const(uplo),
        n, A, lda, B, ldb);
    if (sequence->status == PLASMA_SUCCESS && info != 0)
      plasma_sequence_flush(quark, sequence, request, iinfo+info);
}

ext: c


text:
/**
 * Set of functions to handle object slitless background subtraction
 *
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <unistd.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_multifit.h>

#include "fitsio.h"
#include "aXe_grism.h"
#include "aXe_utils.h"
#include "spce_sect.h"
#include "spce_fitting.h"
#include "spce_is_in.h"
#include "spc_back.h"

#define SQR(x) ((x)*(x))
#define MIN(x,y) (((x)<(y))?(x):(y))
#define MAX(x,y) (((x)>(y))?(x):(y))

/**
 * Function: is_pt_in_a_beam
 * The function checks whether a particular image pixel is part
 * of a beam or not. The pixel is checked against an array
 * of beams it could be part of.
 *
 * Parameters:
 * @param apoint  - the point to check
 * @param iids    - an array of is_in_descriptors for all the beams that
 *                  should be checked
 * @param tnbeams - number of beams in iids
 *
 * Returns:
 * @return 1/0   - fixed values
 */
int
is_pt_in_a_beam (const px_point * const apoint,
                 const is_in_descriptor * const iids, const int tnbeams)
{
  int i;

  for (i = 0; i < tnbeams; i++)
    {
      if (apoint->x < (iids + i)->mini)
        continue;
      if (apoint->x > (iids + i)->maxi)
        continue;
      if (apoint->y < (iids + i)->minj)
        continue;
      if (apoint->y > (iids + i)->maxj)
        continue;
      if (is_in (apoint->x, apoint->y, iids + i))
        {
          return 1;
        }
    }
  return 0;
}

/**
 * Function: get_window_points
 * The subroutine searches in a window around the tracepoint
 * for pixels which are suited for the background determination.
 * The extend of the window is specified in the beam structured,
 * and pixel which are part of any beam can not be used.
 *
 * Parameters:
 * @param  obs      - the observation
 * @param  bck_mask - the background mask
 * @param  actbeam  - the beam to search interp points for
 * @param  tr_point - the integer trace point
 *
 * Returns:
 * @return ret      - the vector with the row numbers
 */
/*
gsl_vector_int *
get_window_points(observation *obs, gsl_matrix *bck_mask, beam actbeam,
                  px_point tr_point)
{
  gsl_vector_int *tmp;
  gsl_vector_int *ret;

  int np_act=0;
  int l_act;
  int u_act;
  int ii;
  int ncols;
  int onbeam=0;

  double umax, lmax;

  // extract the lower and upper window extend
  umax = actbeam.backwindow.x;
  lmax = actbeam.backwindow.y;

  // get the number of columns
  ncols=obs->grism->size2;

  // allocate the vector
  tmp = gsl_vector_int_alloc((int)ceil(lmax) + (int)ceil(umax) + 1);

  // limit the starting point of the search
  // to values within the image dimension
  tr_point.y = MAX(0,tr_point.y);
  tr_point.y = MIN((int)obs->grism->size2,tr_point.y);

  // initialize the row numbers
  // to search up- and downwards
  l_act = tr_point.y -1;
  u_act = tr_point.y;

  if (!gsl_matrix_get(bck_mask,tr_point.x,tr_point.y))
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__,"The closest pixel to the",
                 "trace at %i, %i must be part of an object!",tr_point.x,
                 tr_point.y );

  // search downwards from the trace,
  // assume to start on the beeam
  onbeam=1;
  while (l_act > -1 && fabs((double)l_act - (double)tr_point.y) < lmax)
    {
      // if the pixel is not part of any beam and is not NAN
      if (gsl_matrix_get(bck_mask,tr_point.x,l_act) == 0 &&
          !isnan(gsl_matrix_get(obs->grism, tr_point.x, l_act)))
        {
          // store the pixel index
          gsl_vector_int_set(tmp, np_act, l_act);

          if (onbeam)
            onbeam = 0;

          // enhance the counter
          np_act++;
        }
      // if you are still on the beam
      else if (onbeam && gsl_matrix_get(bck_mask,tr_point.x,l_act))
        {
          // store the pixel index
          gsl_vector_int_set(tmp, np_act, l_act);

          // enhance the counter
          np_act++;
        }
      // decrease the search index
      l_act--;
    }

  // search upwards from the trace,
  // start always searching on the beam
  onbeam=1;
  while (u_act < ncols && fabs((double)u_act - (double)tr_point.y) < umax)
    {
      // if the pixel is not part of any beam and is not NAN
      if (gsl_matrix_get(bck_mask,tr_point.x,u_act) == 0 &&
          !isnan(gsl_matrix_get(obs->grism, tr_point.x, u_act)))
        {
          // store the pixel index
          gsl_vector_int_set(tmp, np_act, u_act);

          if (onbeam)
            onbeam = 0;

          // enhance the counter
          np_act++;
        }
      // if you are still on the beam
      else if (onbeam && gsl_matrix_get(bck_mask,tr_point.x,u_act))
        {
          // store the pixel index
          gsl_vector_int_set(tmp, np_act, u_act);

          // enhance the counter
          np_act++;
        }

      // enhance the search index
      u_act++;
    }

  // check whether the beam extends over
  // the image. Add an artificial start or
  // end point if necessary
  if (l_act < 0)
    gsl_vector_int_set(tmp, np_act++, -1);
  if (u_act > ncols-1)
    gsl_vector_int_set(tmp, np_act++, ncols);

  if (np_act >0)
    {
      // transfer the row numbers to a
      // vector of the right size
      ret = gsl_vector_int_alloc(np_act);
      for (ii=0; ii < np_act; ii++)
        gsl_vector_int_set(ret, ii, gsl_vector_int_get(tmp, ii));

      // sort the row numbers
      gsl_sort_vector_int(ret);

    }
  else
    {
      ret = NULL;
    }

  gsl_vector_int_free(tmp);

  //  return the result
  return ret;
}
*/

/**
 * Function: get_interp_points
 * The function searches for 2*n interpolation points above and
 * below the trace. The search is iteratively from the trace
 * in both directions. In case that the fram borders are met,
 * the respective direction is "closed". The subroutine
 * then tries to get more background pixels on the other
 * size to reach the desired number.
 *
 * Parameters:
 * @param  obs      - the observation
 * @param  bck_mask - the background mask
 * @param  np       - the desired number of interp. points
 * @param  tr_point - the integer trace point
 *
 * Returns:
 * @return ret      - the vector with the row numbers
 */
gsl_vector_int *
get_interp_points(observation *obs, gsl_matrix *bck_mask,
                  int np, px_point tr_point)
{
  gsl_vector_int *tmp;
  gsl_vector_int *ret;
  int np_act=0;
  //int y_low;
  //int y_upp;
  int l_space=1;
  int u_space=1;
  int l_act;
  int u_act;
  int l_np=0;
  int u_np=0;
  int ncols=obs->grism->size2;
  int ii;

  // allocate the vector
  tmp = gsl_vector_int_alloc(2*np+2);

  // limit the starting ppoint of the search
  // to values within the image dimension
  tr_point.y = MAX(0,tr_point.y);
  tr_point.y = MIN((int)obs->grism->size2,tr_point.y);

  // initialize the row numbers
  // to search up- and downwards
  l_act = tr_point.y -1;
  u_act = tr_point.y;

  // as long as interpolation points are missing
  // and one direction, either up or down,
  // is 'open', continue searching
  while (np_act < 2*np && (l_space || u_space))
    {

      // check whether the direction
      // downwards is still open
      if (l_space)
        {
          // move downward until you either find
          // and interp. point or the end of the frame
          while (l_act > -1 && (gsl_matrix_get(bck_mask,tr_point.x,l_act)
                                !=0 ||
                                isnan(gsl_matrix_get(obs->grism,
                                                     tr_point.x, l_act))))
            {
              l_act--;
            }

          // if you are at the end of the frame
          if (l_act < 0)
            {
              // close the direction downwards
              l_space=0;
            }
          else
            {
              // else store the interpolation point,
              // do the various increments
              gsl_vector_int_set(tmp, np_act, l_act);
              np_act++;
              l_np++;
              l_act--;
            }
        }

      // check whether the direction
      // upwards is still open
      if (u_space)
        {
          // move upward until you either find
          // and interp. point or the end of the frame
          while (u_act < ncols && (gsl_matrix_get(bck_mask,tr_point.x,u_act)
                                   !=0 ||
                                   isnan(gsl_matrix_get(obs->grism,
                                                        tr_point.x, u_act))))
            {
              u_act++;
            }

          // if you are at the end of the frame
          if (u_act >= ncols)
            {
              // close the direction upwards
              u_space=0;
            }
          else
            {
              // else store the interpolation point,
              // do the various increments
              gsl_vector_int_set(tmp, np_act, u_act);
              np_act++;
              u_np++;
              u_act++;
            }
        }
    }

  // check whether the beam extends over
  // the image. Add an artificial start or
  // end point if necessary
  if (!l_np)
    gsl_vector_int_set(tmp, np_act++, -1);
  if (!u_np)
    gsl_vector_int_set(tmp, np_act++, ncols);

  // transfer the row numbers to a
  // vector of the right size
  ret = gsl_vector_int_alloc(np_act);
  for (ii=0; ii < np_act; ii++)
    gsl_vector_int_set(ret, ii, gsl_vector_int_get(tmp, ii));

  // sort the row numbers
  gsl_sort_vector_int(ret);

  // release memory
  gsl_vector_int_free(tmp);

  //  return the result
  return ret;
}

/**
 * Function: compute_background
 * The function extracts possible background pixels for a given
 * beam using a specified interpolation functions.
 * After possibly rejecting cosmics, the background is
 * interpolated on the areas which are masked out. The interpolated
 * values are filled into a background structure. If kappa-sigma klipping
 * is applied, the clipped pixels are flagged in the dq-array
 * of the background image.
 *
 * Parameters:
 * @param  obs         - the object list
 * @param  actbeam     - the beam to compute the background for
 * @param  bck_mask    - the background mask
 * @param  fib         - the baground structure
 * @param  npoints     - the number of interpolation points
 * @param  interporder - the interpolation order
 * @param  niter_med   - the number of iterations around the median
 * @param  niter_fit   - the number of iterations around the fit
 * @param  kappa       - the kappa value
 */
void
compute_background (observation *obs, beam actbeam, gsl_matrix *bck_mask,
                    fullimg_background *fib, int npoints,
                    int interporder, const int niter_med,
                    const int niter_fit, const double kappa)
{
  px_point xborder;
  gsl_vector_int *yvec;
  trace_func *tracefun;
  px_point tpoint;
  double *ys, *fs, *ws, *yi;
  //double var;

  int i, ii;
  int j, n;
  int min_y, max_y;

  // define the beam and the trace function
  tracefun = actbeam.spec_trace;

  /* If this beam's ignore flag is set to 1 then do nothing */
  if (actbeam.ignore == 1)
    return;

  // determine the start and end point in x
  xborder = get_xrange(obs, actbeam);

  // Loop over all columns
  for (i = xborder.x; i < xborder.y; i++)
    {

      // determine the pixel closest to the trace
      tpoint.x = i;
      tpoint.y = (int)floor(tracefun->func((double)i-actbeam.refpoint.x,
                                           tracefun->data)
                            + actbeam.refpoint.y+0.5);

      // determine the interpolation points
      // around the trace
      yvec = get_interp_points(obs, bck_mask, npoints, tpoint);

      //-----------------------------------------------------------------
      //   some code for FORS2 MXU
      //      if (actbeam.backwindow.x < 0.0)
      //        yvec = get_interp_points(obs, bck_mask, npoints, tpoint);
      //      else
      //        yvec = get_window_points(obs, bck_mask, actbeam, tpoint);

      // give a warning and go to the next
      // column if there are no background points
      if (!yvec)
        {
          aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
                       "No backgound points could be found for beam %C. "
                       "at collumn %d %d",
                       BEAM(actbeam.ID), i, tpoint.y );
          continue;
        }

      // extract maximum and minimum 0f the y-values
      min_y = gsl_vector_int_get(yvec, 0);
      max_y = gsl_vector_int_get(yvec, yvec->size-1);
      //      if (actbeam.backwindow.x == 23.3 && actbeam.backwindow.y == 4.0)
      //        {
          //    fprintf(stdout, "%i <--> %i; ", min_y, max_y);
      //      if (min_y > tpoint.y - actbeam.width
      //          || max_y < tpoint.y + actbeam.width)
      //        fprintf(stdout, "%i %f <--> %f %i;  ", min_y, tpoint.y - actbeam.width,
      //                tpoint.y + actbeam.width, max_y);
              //        }

      // determine the size of the double vectors
      // to make the background determination
      // allocater the space and initialize
      // all values to default
      n = max_y - min_y + 1;
      ys = (double *) malloc (n * sizeof (double));
      fs = (double *) malloc (n * sizeof (double));
      ws = (double *) malloc (n * sizeof (double));
      yi = (double *) malloc (n * sizeof (double));
      for (ii = 0; ii < n; ii++)
        {
          ys[ii] = min_y + ii;
          fs[ii] = 0.0;
          ws[ii] = 0.0;
        }

      // transfer the values from the image column
      // to the double vectors, set the weight
      for (ii = 0; ii < (int)yvec->size; ii++)
        {
          // extract the row number
          j = gsl_vector_int_get (yvec, ii);

          // check whether the row is inside the imag
          // and whether the is no object on the pixel
          if ((j != -1) && (j != (int)obs->grism->size2)
              && !gsl_matrix_get(bck_mask, tpoint.x, j))
            {
              // set the values and the weight
              fs[j - min_y] = gsl_matrix_get (obs->grism, tpoint.x, j);
              ws[j - min_y] = gsl_matrix_get (obs->pixerrs, tpoint.x, j);
            }
        }

      // iterate on the background points to
      // reject e.g. cosmics
      if (niter_med > 0 || niter_fit > 0)
        {
          // iterate on the median
          for (j=0; j < niter_med; j++)
            kappa_sigma_clipp(ys, fs,ws, n,kappa, obs, tpoint.x);

          // iterate on the fit
          if (niter_fit > 0)
            comp_kappasigma_interp( ys, fs, ws, n, interporder,
                                    niter_fit, kappa, obs, tpoint.x);
        }

      // do the final background determination
      comp_vector_interp( ys, fs, ws, yi, n, interporder, 1);

      // copy the intepolated values
      // to the background matrix
      for (j = 0; j < n; j++)
        {
          if ((ys[j] < 0) || (ys[j] >= obs->grism->size2))
            continue;
          gsl_matrix_set (fib->bck, tpoint.x, (int) floor (ys[j]), fs[j]);
          gsl_matrix_set (fib->err, tpoint.x, (int) floor (ys[j]), ws[j]);
        }

      // release memory
      free (ys);
      ys = NULL;
      free (fs);
      fs = NULL;
      free (yi);
      yi = NULL;
      free (ws);
      ws = NULL;

      // release memory
      gsl_vector_int_free(yvec);
    }
}

/**
 * Function: get_xrange
 * The subroutine determines the extend of a beam in x-direction.
 * The minimum and maximum value in x of pixels which are part
 * of the particular beam are determined and returned.
 *
 * Parameters:
 * @param  obs     - the object list
 * @param  actbeam - the beam to determine the extent for
 *
 * Returns:
 * @return ret     - Min/Max values of the beam in x
 */
px_point
get_xrange(observation *obs, beam actbeam)
{

  px_point ret;

  // Find the object starting and ending column
  // for the beam of interest actbeam.corners
  ret.x = MIN (MIN (MIN (actbeam.corners[0].x, actbeam.corners[1].x),
                    actbeam.corners[2].x), actbeam.corners[3].x);
  ret.y  = MAX (MAX (MAX (actbeam.corners[0].x, actbeam.corners[1].x),
                     actbeam.corners[2].x), actbeam.corners[3].x);

  // limit the start and end column
  // to the image size
  ret.x = MAX(ret.x, 0);
  ret.y = MIN(ret.y, (int)obs->grism->size1);

  // warn if the beam is completely outside
  // the grism image
  if (ret.x > (int)obs->grism->size1 || ret.y < 0)
    {
      aXe_message (aXe_M_WARN4, __FILE__, __LINE__,
                   "Object  is not in the image start_i:%d end_i:%d",
                   ret.x, ret.y);
    }

  // return the result
  return ret;
}

/**
 * Function: comp_kappasigma_interp
 * The function performs a kappa-sigma clipping rejection
 * on a set of data given in vectors for indipendent, dipendent
 * and weight values. The differences to base the clipping on
 * is "true value minus background value", where background
 * value is determined using the identical algorithm as the
 * final background determination.
 *
 * Parameters:
 * @param xs     - array for independent values
 * @param ys     - array for dependent values
 * @param ws     - weight array
 * @param n      - number of pixels
 * @param interp - number indicating the interpolation scheme
 * @param niter  - the number of iterations
 * @param kappa  - the kappa value for rejection
 * @param obs    - the object list
 * @param colnum - the column number
 */
void
comp_kappasigma_interp(const double *const xs, double *const ys,
                       double *const ws, const int n,
                       const int interp, const int niter, const double kappa,
                       observation *obs, int colnum)
{
  double *ys_tmp;
  double *yi_tmp;
  double *y_diff;

  double stdev;

  int *iindex;

  int i, m, j;

  // allocate temporary vectors
  ys_tmp = (double *) malloc (n * sizeof (double));
  yi_tmp = (double *) malloc (n * sizeof (double));
  y_diff = (double *) malloc (n * sizeof (double));
  iindex = (int *) malloc (n * sizeof (int));

  // transfer the dependent value
  // to the tmp vector
  for (i=0; i < n; i++)
    ys_tmp[i] = ys[i];

  // do niter times
  for (j=0; j < niter; j++)
    {

      // make the background determination
      comp_vector_interp(xs, ys_tmp, ws, yi_tmp, n, interp, 0);

      // calculate for all background
      // pixels the differences between
      // the original and background value
      m=0;
      for (i=0; i < n; i++)
        {
          if (ws[i] != 0.0)
            {
              y_diff[m] = ys[i] - yi_tmp[i];
              iindex[m] = i;
              m++;
            }
        }

      // compute the standard deviation
      // on the differences
      stdev = gsl_stats_sd (y_diff, 1, m);

      // do the clipping
      for (i=0; i < m; i++)
        {
          // check for pixels to exclude
          if (fabs(y_diff[i]) > kappa*stdev)
            {
              // set its weight to 0.0
              ws[iindex[i]] = 0.0;
              // mark the pixel in the dq-array
              if ( xs[iindex[i]] >= 0 &&  xs[iindex[i]] < obs->dq->size2)
                gsl_matrix_set(obs->dq, colnum, xs[iindex[i]], DQ_KAPPA_SIGMA);
            }
        }
    }

  for (i=0; i < n; i++)
    ys[i] = ys_tmp[i];

  // free memory
  free(ys_tmp);
  free(yi_tmp);
  free(y_diff);
  free(iindex);
}

/**
 * Function: kappa_sigma_clipp
 * The subroutine performs one kappa-sigma step on the data
 * given in various vectors for independent, dipendent and
 * weight data. The difference to apply the clipping on
 * is "original value <minus> median of the data set".
 *
 * Parameters:
 * @param xs     - array for independent values
 * @param ys     - array for dependent values
 * @param ws     - weight array
 * @param n      - number of pixels
 * @param kappa  - the kappa value for rejection
 * @param obs    - the object list
 * @param colnum - the column number
 */
void
kappa_sigma_clipp(const double *const xs, double *const ys, double *const ws,
                  const int n, const double kappa, observation *obs,
                  int colnum)
{
  double *ys_tmp, *ys_med;

  int *iindex;

  int ii, m=0, npixel=0;

  double median=0.0, stdev=0.0;

  // allocate memory for temporay vectors
  ys_tmp = (double *) malloc (n * sizeof (double));
  ys_med = (double *) malloc (n * sizeof (double));
  iindex = (int *) malloc (n * sizeof (int));

  // store the background values
  // in the temporary vectors
  for (ii=0; ii < n; ii++)
    {
      if (ws[ii] != 0.0)
        {
          ys_tmp[m] = ys[ii];
          ys_med[m] = ys[ii];
          iindex[m] = ii;
          m++;
        }
    }

  // derive median and standard deviation
  gsl_sort (ys_tmp, 1, m);
  median = gsl_stats_median_from_sorted_data(ys_tmp, 1, m);
  stdev = gsl_stats_sd_m (ys_med, 1, m, median);

  // apply the clipping
  npixel=0;
  for (ii=0; ii < m; ii++)
    {
      // check whether the pixel should
      // be clipped
      if (fabs(ys_med[ii]-median) > kappa*stdev)
        {
          // set the weight of a clipped pixel to 0.0
          ws[iindex[ii]] = 0.0;
          // mark the clipped pixel in the dq-array
          gsl_matrix_set(obs->dq, colnum, xs[iindex[ii]], DQ_KAPPA_SIGMA);
        }
      else
        {
          npixel++;
        }

    }

  // release memory
  free(ys_med);
  free(ys_tmp);
  free(iindex);
}

/**
 * Function: comp_vector_interp
 * The function passes the interpolation data to the
 * desired interpolating function.
 *
 * Parameters:
 * @param xs     - double vector containing the x values
 * @param ys     - double vector containing the y values
 * @param ws     - double vector containing the weights associated with ys
 * @param n      - number of points in xs,ys, and ws (must be greater than m!)
 * @param interp - desired interpolation type
 * @param final  - indicates a final interpolation
 */
void
comp_vector_interp(const double *const xs, double *const ys,
                   double *const ws, double *const yi, const int n,
                   const int interp, const int final)
{
  /* Median the background */
  if (interp == -1)
    {
      comp_vector_median(xs, ys, ws, yi, n, final);
    }
  /* Straight average of the background */
  else if (interp == 0)
    {
      comp_vector_average(xs, ys, ws, yi, n, final);
    }
  /* Linear interpolation of the background */
  else if (interp == 1)
    {
      comp_vector_linear(xs, ys, ws, yi, n, final);
    }
  /* n(>1) order interpolation of the background */
  else if (interp > 1)
    {
      comp_vector_polyN (interp + 1, xs, ys, ws, yi, n, final);
    }
  else
    {
      aXe_message(aXe_M_FATAL, __FILE__, __LINE__,
                  "Do not know what to do with interpolation: %i %s\n", interp);
    }
}

void
compute_global_background (object **oblist, const int obj_index,
                           gsl_matrix *bck_mask, fullimg_background * fib,
                           int interporder)
{
  int i, j, n;
  double *ys, *fs, *ws, *yi;
  //double *ws0;
  observation *grism = oblist[obj_index]->grism_obs;
  //long ma;
  //double var;

  for (i = 0; i < (int)grism->grism->size1; i++)
    {
      /* Loop over the columns of interest */

      n = grism->grism->size2;

      ys = (double *) malloc (n * sizeof (double));
      fs = (double *) malloc (n * sizeof (double));
      ws = (double *) malloc (n * sizeof (double));
      yi = (double *) malloc (n * sizeof (double));

      for (j = 0; j < n; j++)
        {
          ys[j] = j;
          fs[j] = 0.0;
          ws[j] = 0.0;
          if ((gsl_matrix_get(bck_mask,i,j)==0) &&
              (!(isnan(gsl_matrix_get(grism->grism, i,j)))))
            {
              fs[j]=gsl_matrix_get (grism->grism,i,j);
              ws[j]=gsl_matrix_get (grism->grism,i,j);
            }
        }

      /* Median the background */
      if (interporder == -1)
        {
          double *tmp, med, std;
          int nn = 0;
          for (j = 0; j < n; j++)
            {
              if (ws[j] != 0.0)
                {
                  nn++;
                }
            }
          tmp = malloc (nn * sizeof (double));
          nn = 0;
          for (j = 0; j < n; j++)
            {
              if (ws[j] != 0.0)
                {
                  tmp[nn] = fs[j];
                  nn++;
                }
            }
          gsl_sort (tmp, 1, nn);
          med = gsl_stats_median_from_sorted_data (tmp, 1, nn);
          std = gsl_stats_sd (tmp, 1, nn);
          //fprintf(stderr,"med: %g\n",med);
          for (j = 0; j < n; j++)
            {
              if (ws[j] == 0.0)
                {
                  fs[j] = med;
                  ws[j] = std;
                }
              //              else
              //                {
              //                  ws[j] = 1.0/sqrt(ws[j]);
              //                }
            }
          free (tmp);
          tmp = NULL;
        }

      /* Straight average of the background */
      if (interporder == 0)
        {
          double *tmp, sum = 0.0, avg = 0.0, std = 0.0;
          int nn = 0;

          nn = 0;
          for (j = 0; j < n; j++)
            {
              if (ws[j] != 0.0)
                {
                  nn++;
                }
            }
          tmp = malloc (nn * sizeof (double));
          nn = 0;
          for (j = 0; j < n; j++)
            {
              if (ws[j] != 0.0)
                {
                  sum += fs[j];
                  tmp[nn] = fs[j];
                  nn++;
                }
            }
          if (nn > 0)
            avg = sum / nn;
          std = gsl_stats_sd (tmp, 1, nn);
          for (j = 0; j < n; j++)
            {
              if (ws[j] == 0.0)
                {
                  fs[j] = avg;
                  ws[j] = std;
                }
              //              else
              //                {
              //                  ws[j] = 1.0/sqrt(ws[j]);
              //                }
            }
          free (tmp);
          tmp = NULL;
        }

      /* Linear interpolation of the background */
      if (interporder == 1)
        {
          comp_vector_linear (ys, fs, ws, yi, n, 1);
          //      fit_vector_linear_t (ys, fs, ws, n);
        }

      /* n(>1) order interpolation of the background */
      if (interporder > 1)
        {
          comp_vector_polyN (interporder + 1, ys,fs, ws, yi, n, 1);
      //          fit_vector_poly_N_t (interporder + 1, ys, fs, ws, n);
        }

      for (j = 0; j < n; j++)
        {
          if ((ys[j] < 0) || (ys[j] >= grism->grism->size2))
            continue;
            gsl_matrix_set (fib->bck, i, j,fs[j]);
            gsl_matrix_set (fib->err, i, j,ws[j]);
        }
      free (ys);
      ys = NULL;
      free (fs);
      fs = NULL;
      free (ws);
      ws = NULL;
    }
}

/**
 * Function: fullimg_background_function
 * Returns the background level for the point x, y if a complete
 * gsl_matrix of the background is available.  This function is
 * exposed to the outside through a pointer in the observation
 * structure.
 *
 * Parameters:
 * @param x ad nauseam
 * @param y ad nauseam
 * @param val pointer to a double to leave the background value in
 * @param err pointer to a double to leave the absolute error of val in
 * @param pars a gsl_matrix containing the background
 */
void
fullimg_background_function (const int x, const int y, PIXEL_T * const val,
                             PIXEL_T * const err,
                             const background * const back)
{
  const fullimg_background *const fib = back->pars;

  *val = gsl_matrix_get (fib->bck, (int) rint (x), (int) rint (y));
  if (fib->err)
    {
      *err = gsl_matrix_get (fib->err, (int) rint (x), (int) rint (y));
    }
  else
    {
      *err = 0;
    }
}

/**
 * Function: compute_fullimg_background
 * Computes a background image. All beams with ignore=1 are
 * completely ignored. Anything else is not.
 *
 * Parameters:
 * @param obs         - a pointer to the observation structure to fill out
 * @param oblist      - a list of all objects of the observation
 * @param npoints     - number of points to examine when fitting the background
 * @param interporder - order of the polynomial to fit to the background
 * @param niter_med   - order of the polynomial to fit to the background
 * @param niter_fit   - order of the polynomial to fit to the background
 * @param kappa       - order of the polynomial to fit to the background
 *
 * Returns:
 * @return background - an allocated background structure
 *
 */
background *
compute_fullimg_background (observation *obs, object **oblist,
                            int npoints, int interporder, const int niter_med,
                            const int niter_fit, const double kappa,
                            int nor_flag, const int sm_length, const double fwhm)
{
  fullimg_background *fib;
  background *backg;
  gsl_matrix *bck_mask;
  int i, j;
  //object *const *obp;

  // allocate space for the backgrounds
  fib   = (fullimg_background *)malloc (sizeof (fullimg_background));
  fib->bck = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);
  gsl_matrix_set_all (fib->bck, 0.);
  fib->err = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);
  gsl_matrix_set_all (fib->err, 0.);

  if (nor_flag)
    {
      for (i = 0; i < (int)fib->bck->size1; i++)
        {
          for (j = 0; j < (int)fib->bck->size2; j++)
            {
              gsl_matrix_set (fib->bck, i, j,
                              gsl_matrix_get (obs->grism, i, j));
              gsl_matrix_set (fib->err, i, j,
                              gsl_matrix_get (obs->pixerrs, i, j));
            }
        }
    }

  // allocate memory
  backg = (background *)malloc (sizeof (background));

  // create the mask image
  bck_mask = aperture_mask(obs,oblist);

  // in case that pixels may get dq values,
  // make sure that there will be a dq-array
  //  if (niter_med > 0 || niter_fit > 0)
  //    {
  if (!obs->dq)
    obs->dq = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);

  // initialize the dq-array
  //      gsl_matrix_set_all (obs->dq, 0.0);
  //    }

  if (oblist != NULL)
    {
      // Now compute background for each beam, one after the other
      // go over the whole object list
      i=0;
      while (oblist[i] != NULL) {
        // go over each beam
        for (j = 0; j < oblist[i]->nbeams; j++)
          {
            // check for beams to be neglected
            if (oblist[i]->beams[j].ignore == 1)
              {
                continue;
              }
            else
              {
                // start the background interpolation
                // for a specific beam
                fprintf(stdout,"Computing background of BEAM %d%c.",
                        oblist[i]->ID,BEAM(oblist[i]->beams[j].ID));
                compute_background(obs, oblist[i]->beams[j], bck_mask,
                                   fib, npoints, interporder, niter_med,
                                   niter_fit, kappa);
                fprintf(stdout," Done.\n");
              }
          }
        // increment the counter
        i++;
      }
    }
  // replace all NAN's with values 0.0
  for (i = 0; i < (int)fib->bck->size1; i++)
    {
      for (j = 0; j < (int)fib->bck->size2; j++)
        {
          if(isnan(gsl_matrix_get (fib->bck, i, j)))
            gsl_matrix_set (fib->bck, i, j,0.0);
          if(isnan(gsl_matrix_get (fib->err, i, j)))
            gsl_matrix_set (fib->err, i, j,0.0);
        }
    }

  if (sm_length && fwhm)
    gsmooth_background (bck_mask, sm_length, fwhm, fib);

  // put together the result
  backg->pars = fib;
  backg->bck_func = &fullimg_background_function;

  // free space
  gsl_matrix_free(bck_mask);

  // return the result
  return backg;
}

/**
 * Function: compute_backsub_mask
 * Computes a mask image for the background subtraction. All pixels
 * covered by at least one beam are set to a vale -100000 to distinguish
 * them from pixels which are part of the background.
 *
 * Parameters:
 * @param obs a pointer to the observation structure to fill out
 * @param oblist a list of all objects of the observation
 * @param npoints number of points to examine when fitting the background
 * @param interporder order of the polynomial to fit to the background
 *
 * Returns:
 * @return an allocated background structure
 *
 */
background *
compute_backsub_mask (observation *obs, object **oblist)
{
  fullimg_background *fib;// = malloc (sizeof (fullimg_background));
  background *backg; // = malloc (sizeof (background));
  gsl_matrix *bck_mask;
  int i, j;
  //object *const *obp;
  //int tnbeams;

  // make the mask image
  bck_mask = aperture_mask(obs,oblist);

  // allocate memory for the return structure
  // and fill in some dummy values
  fib = malloc (sizeof (fullimg_background));
  fib->bck = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);
  gsl_matrix_set_all (fib->bck, 0.);
  fib->err = NULL;

  // allocate memory for the return structure
  backg = malloc (sizeof (background));

  // transfer the pixel values from the grism image
  // to the mask image
  for (i = 0; i < (int)fib->bck->size1; i++)
    {
      for (j = 0; j < (int)fib->bck->size2; j++)
        {
          gsl_matrix_set (fib->bck, i, j,
                          gsl_matrix_get (obs->grism, i, j));
        }
    }

  // set pixels occupied by beams
  // to the value -10000000
  for (i = 0; i < (int)fib->bck->size1; i++)
    {
      for (j = 0; j < (int)fib->bck->size2; j++)
        {

          if (gsl_matrix_get (bck_mask, i, j))
            gsl_matrix_set (fib->bck, i, j, -1000000.0);
        }
    }

  // release the dq-structure,
  // otherwise it is saved to
  // the mask image
  if (obs->dq != NULL)
    {
      gsl_matrix_free (obs->dq);
      obs->dq=NULL;
    }

  // release memory
  gsl_matrix_free(bck_mask);

  // assemble the return structure
  backg->pars = fib;
  backg->bck_func = &fullimg_background_function;

  // return the result
  return backg;
}

/**
 * Function: compute_fullimg_global_background
 * The subroutine computes the background image based on all pixels
 * in a column which are not part of an object.
 *
 * Parameters:
 * @param obs         - a pointer to the observation structure to fill out
 * @param oblist      - a list of all objects of the observation
 * @param interporder - order of the polynomial to fit to the background
 *
 * Returns:
 * @return background - an allocated background structure
 *
 */
background *
compute_fullimg_global_background(observation *obs, object **oblist,
                                  int interporder, const int sm_length, const double fwhm)
{
  fullimg_background *fib;
  background *backg;
  gsl_matrix *bck_mask;
  int i, j;
  object *const *obp;
  int tnbeams;

  // allocate memory
  fib = (fullimg_background *)malloc (sizeof (fullimg_background));
  backg = (background *)malloc (sizeof (background));

  // allocate memory
  bck_mask = aperture_mask(obs,oblist);
  fib->bck = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);
  fib->err = gsl_matrix_alloc (obs->grism->size1, obs->grism->size2);

  // initialize the new arrays
  gsl_matrix_set_all (fib->bck, 0.);
  gsl_matrix_set_all (fib->err, 0.);

  /* Count beams in observation */
  if (oblist!=NULL)
    {
      tnbeams = 0;
      for (obp = oblist; *obp; obp++)
        {
          for (i = 0; i < (*obp)->nbeams; i++)
            {
              // if (((*obp)->beams[i]).ignore!=1)
              tnbeams++; /* Count ALL the beams */
            }
        }
    }

  /* Now compute global background using first beam grism info */
  if (oblist != NULL) {
    compute_global_background (oblist, 0, bck_mask,
                               fib, interporder);
  }

  // replace all NAN's with values 0.0
  for (i = 0; i < (int)fib->bck->size1; i++)
    {
      for (j = 0; j < (int)fib->bck->size2; j++)
        {
          if(isnan(gsl_matrix_get (fib->bck, i, j)))
            {
              gsl_matrix_set (fib->bck, i, j,0.0);
            }
          if(isnan(gsl_matrix_get (fib->err, i, j)))
            {
              gsl_matrix_set (fib->err, i, j,0.0);
            }
        }
    }

  // make a Gaussian smoothing
  // if requested
  if (sm_length && fwhm)
    gsmooth_background (bck_mask, sm_length, fwhm, fib);

  // compose the result structure
  backg->pars = fib;
  backg->bck_func = &fullimg_background_function;

  // return the result
  return backg;
}

/**
 * Function: free_fullimg_background
 * Frees a background structure allocated using make_fullimg_background
 *
 * Parameters:
 * @param backg - a pointer to the background structure to free
 */
void
free_fullimg_background (background * backg)
{
  fullimg_background *fib = backg->pars;

  gsl_matrix_free (fib->bck);
  if (fib->err)
    {
      gsl_matrix_free (fib->err);
    }
  free (backg);
  backg = NULL;
}

/**
 * Function: aperture_mask
 * This functions returns an image mask where pixels that are within
 * an aperture are set to the number of aperture they appear.
 * Skips any beams which are set to be ignored (ignore=1)
 *
 * Parameters:
 * @param  obs    - a pointer to an observation structure
 * @param  oblist - a pointer to an object list
 *
 * Returns:
 * @return bck   - a pointer to a gsl_matrix
 */
gsl_matrix *
aperture_mask (observation * const obs, object **oblist)
{

  int x, y,i=0,j;

  double xrel, yrel, width;

  is_in_descriptor iid;
  px_point ll, ur;
  gsl_matrix *bck;
  sectionfun sf;

  // create the image,
  // set all values to zero
  bck = gsl_matrix_alloc(obs->grism->size1,obs->grism->size2);
  gsl_matrix_set_all(bck,0.0);

  // Return a zero image when there is no
  // beam in the list
  if (oblist==NULL)
    return bck;

  // go over each object
  while(oblist[i]!=NULL)
    {
      // go over each beam
      for(j=0;j<oblist[i]->nbeams;j++)
        {
          // continue if the beam is to be ignored
          if ((oblist[i]->beams[j]).ignore==1)
            continue;

          // check whther the trace is second order or
          // even higher
          if ((oblist[i]->beams[j]).spec_trace->type > 1)
            {
              // create the section function for second order
              // traces such as FORSII
              fill_in_sectionfun (&sf, (oblist[i]->beams[j]).orient,
                                  &(oblist[i]->beams[j]));
            }
          else
            {
              // create the descriptor to set up the quadrangle routines
              fill_is_in_descriptor (&iid, (oblist[i]->beams[j]).corners);
            }

          // check for the corners of the quadrangle.
          // go over each point in x and y.
          quad_to_bbox ((oblist[i]->beams[j]).corners, &ll, &ur);
          for (x = ll.x; x <= ur.x; x++)
            {
              for (y = ll.y; y <= ur.y; y++)
                {

                  // neglect the pixel if it is a]outside the
                  // image area
                  if ((x < 0) || (y < 0) || (x >= (int)obs->grism->size1)
                      || (y >= (int)obs->grism->size2))
                    continue;

                  // check the order of the trace polynomial
                  if ((oblist[i]->beams[j]).spec_trace->type > 1)
                    {

                      // determine the trace distance and
                      // enhance the image value if the pixel
                      // satisfies the distance criterium
                      xrel = x-(oblist[i]->beams[j]).refpoint.x;
                      yrel = y-(oblist[i]->beams[j]).refpoint.y;
                      width = (oblist[i]->beams[j]).width+0.5;
                      if (tracedist_criteria(xrel, yrel, &sf, (oblist[i]->beams[j]).spec_trace, width))
                        gsl_matrix_set (bck,x,y,gsl_matrix_get (bck,x,y)+1);
                    }
                  else
                    {
                      //                      if (oblist[i]->ID == 11 && x == 59)
                      //                        fprintf(stdout, "xx: %i, yy: %i: %i\n", x, y, is_in (x, y, &iid));

                      // check whether the pixel is inside the quadrangle
                      // which defines the beam are, and enhance the image
                      // value if yes.
                      if (is_in (x, y, &iid))
                        gsl_matrix_set (bck,x,y,gsl_matrix_get (bck,x,y)+1);
                    }
                }
            }
        }

      // enhance the object counter
      i++;
    }

  // return the resulting image
  return bck;
}

/**
 * Function: background_to_FITSimage
 * Function to write the data and error content of a backgound
 * observation pars component (if it is of type fullimg_background only)
 * into the main HDU (data) and first extension (error) of a FITS file.
 * A GQ array can be appended if an observation with a non NULL DQ is
 * passed to this function
 *
 * Parameters:
 * @param filename - name of the image
 * @param bck      - filled background structure
 * @param obs      - the observation the background is based upon
 */
void
background_to_FITSimage (char filename[], background * bck, observation *obs)
{
  fitsfile *output;
  long naxes[2];
  int f_status = 0;
  PIXEL_T *storage, *dp;
  int x, y;
  fullimg_background *pars;
  int hdunum,hdutype;

  pars = (fullimg_background *) bck->pars;

  // Open the file for creating/appending
  create_FITSimage (filename, 1);
  fits_open_file (&output, filename, READWRITE, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "gsl_to_FITSimage: " "Could not open file: %s",
                   filename);
    }

  // count the number of extentions
  fits_get_num_hdus (output, &hdunum, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "gsl_to_FITSimage: Could not get"
                   " number of HDU from: %s", filename);
    }

  // Move to last HDU
  fits_movabs_hdu (output, hdunum, &hdutype, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "gsl_to_FITSimage: Could not mov"
                   " to HDU number %d in file: %s", hdunum, filename);
    }
  /* Get current HDU number */
  fits_get_hdu_num (output, &hdunum);

  // Deal with the SCI part of the background
  naxes[0] = pars->bck->size1;
  naxes[1] = pars->bck->size2;

  // Allocate storage room
  if (!(storage = malloc (naxes[0] * naxes[1] * sizeof (double))))
    aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory");

  // Populate the storage array
  dp = storage;
  for (y = 0; y < naxes[1]; y++)
    {
      for (x = 0; x < naxes[0]; x++)
        {
          *dp = gsl_matrix_get (pars->bck, x, y);
          dp++;
        }
    }
  /* create HDU extname */
  fits_create_img (output, -32, 2, naxes, &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "gsl_to_FITSimage: Could create SCI HDU in file: %s", filename);
    }
  fits_write_img (output, TFLOAT, 1, naxes[0] * naxes[1], storage,
                  &f_status);
  if (f_status)
    {
      ffrprt (stderr, f_status);
      aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
                   "gsl_to_FITSimage: Could write SCI HDU in file: %s", filename);
    }

  /* write the HDU EXTNAME */
  {
    char comment[FLEN_COMMENT];
    char str[FLEN_KEYWORD];
    strcpy (str, "SCI");
    strcpy (comment, "Extension name");
    fits_write_key_str (output, "EXTNAME", str, comment, &f_status);
  }
  free (storage);
  storage = NULL;

  if (pars->err)
    {
      /* Deal with the error part of the background */
      naxes[0] = pars->err->size1;
      naxes[1] = pars->err->size2;
      /* Allocate storage room */
      if (!(storage = malloc (naxes[0] * naxes[1] * sizeof (double))))
        {
          aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory");
        }
      /* Populate the storage array */
      dp = storage;
      for (y = 0; y < naxes[1]; y++)
        {
          for (x = 0; x < naxes[0]; x++)
            {
              *dp = gsl_matrix_get (pars->err, x, y);
              dp++;
            }
        }
      /* create HDU extname */
      /* Get current HDU number */
      fits_get_hdu_num (output, &hdunum);

      fits_create_img (output, -32, 2, naxes, &f_status);
      fits_write_img (output, TFLOAT, 1, naxes[0] * naxes[1], storage,
                      &f_status);

      /* Get current HDU number */
      fits_get_hdu_num (output, &hdunum);
      /* write the HDU EXTNAME */
      {
        char comment[FLEN_COMMENT];
        char str[FLEN_KEYWORD];
        strcpy (str, "ERR");
        strcpy (comment, "Extension name");
        fits_write_key_str (output, "EXTNAME", str, comment, &f_status);
      }
      free (storage);
      storage = NULL;
    }

  /* Deal with the DQ part of the background */
  if (obs->dq != NULL)
    {
      naxes[0] = obs->dq->size1;
      naxes[1] = obs->dq->size2;
      /* Allocate storage room */
      if (!(storage = malloc (naxes[0] * naxes[1] * sizeof (double))))
        {
          aXe_message (aXe_M_FATAL, __FILE__, __LINE__, "Out of memory");
        }
      /* Populate the storage array */
      dp = storage;
      for (y = 0; y < naxes[1]; y++)
        {
          for (x = 0; x < naxes[0]; x++)
            {
              *dp = gsl_matrix_get (obs->dq, x, y);
              dp++;
            }
        }
      /* create HDU extname */
      fits_create_img (output, 16, 2, naxes, &f_status);
      fits_write_img (output, TFLOAT, 1, naxes[0] * naxes[1], storage,
                      &f_status);
      /* write the HDU EXTNAME */
      {
        char comment[FLEN_COMMENT];
        char str[FLEN_KEYWORD];
        strcpy (str, "DQ");
        strcpy (comment, "Extension name");
        fits_write_key_str (output, "EXTNAME", str, comment, &f_status);
      }
      free (storage);
      storage = NULL;
    }

  /* close file */
  fits_close_file (output, &f_status);
}

/**
 * Function: gsmooth_background
 * Smooth all interpolated pixel in the background using a Gaussian
 * function. The smoothing is done exclusively towards the x-values.
 * The smoothing should help to reduce the noise from the limited number
 * of background pixels.
 *
 * Parameters:
 * @param  bck_mask     - the background mask
 * @param  smoot_length - number of pixels on either side to use for smoothing
 * @param  fwhm         - fwhm of the Gaussian
 * @param  fib          - the background structure
 */
void
gsmooth_background (const gsl_matrix *bck_mask, const int smooth_length,
                    const double fwhm, fullimg_background *fib)
{
  double efactor;

  int ix, iy;

  gsl_vector *pixvalues;
  gsl_vector *weights;
  gsl_vector *pmask;

  gsl_matrix *new_bck;

  // allocate memory for the vectors
  pixvalues = gsl_vector_alloc(2 * smooth_length + 1);
  weights   = gsl_vector_alloc(2 * smooth_length + 1);
  pmask     = gsl_vector_alloc(2 * smooth_length + 1);

  // allocate the new background
  // and initialize it
  new_bck = gsl_matrix_alloc(bck_mask->size1, bck_mask->size2);
  gsl_matrix_set_all(new_bck, 0.0);

  // prepare the Gaussian
  efactor = compute_efactor(fwhm);

  // fill the weights
  for (ix=0; ix < (int)weights->size; ix++)
    gsl_vector_set(weights, ix, compute_gvalue((double)ix - (double)smooth_length, efactor));

  // go over all rows
  for (iy=0; iy < (int)fib->bck->size2; iy++)
    {
      // go over all columns
      for (ix=0; ix < (int)fib->bck->size1; ix++)
        {
          // check whether the pixel IS part of a beam
          if (!(gsl_matrix_get(bck_mask, ix, iy) != 0.0
                && gsl_matrix_get(fib->bck, ix, iy) !=0.0))
            {
              // transfer the old background value
              gsl_matrix_set(new_bck, ix, iy, gsl_matrix_get(fib->bck, ix, iy));
            }
          else
            {
              // fill the pixel and mask values
              fill_pixvalues(bck_mask, smooth_length, fib, ix, iy, pixvalues, pmask);

              // fill in the weighted mean
              gsl_matrix_set(new_bck, ix, iy, get_weighted_mean(pixvalues, weights, pmask));
            }
        }
    }

  // release the allocated dspace
  gsl_vector_free(pixvalues);
  gsl_vector_free(weights);
  gsl_vector_free(pmask);

  // release the memory of the
  // old background
  gsl_matrix_free(fib->bck);

  // transfer the new background
  // to the background structure
  fib->bck = new_bck;
}

/**
 * Function: get_weighted_mean
 * The function computes the weighted mean of values stored in a value vector
 * and a weight vector. As mask vector marks values pixels not to be considered.
 * Using a separate mask vectorhas the advantage the weights can be kept
 * constant and do not have to be re-calculated in repeated runs.
 *
 * Parameters:
 * @param pixvalues - vector with pixel values
 * @param weights   - vector with weights
 * @param pmask     - mask vector
 *
 * Returns:
 * @return sum/www  - the weighted mean
 */
double
get_weighted_mean(const gsl_vector *pixvalues, const gsl_vector *weights,
                  const gsl_vector *pmask)
{
  int index;

  // initialize the total
  // sum and weight
  double sum=0.0;
  double www=0.0;

  for (index=0; index < (int)pixvalues->size; index ++)
    {

      if (gsl_vector_get(pmask, index))
        {
          // enhance the total sum
          sum += gsl_vector_get(pixvalues, index) * gsl_vector_get(weights, index);

          // enhance the total weight
          www += gsl_vector_get(weights, index);
        }
    }

  // return the total sum,
  // divided by the total weight
  return sum / www;
}

/**
 * Function: fill_pixvalues
 * The function provides the essential information for Gaussian smoothing
 * for a single pixel. It fills a vector with the values of all pixels
 * within the smoothing length. Not interpolated pixels are excluded.
 * A mask vector provides the information on which position is filled
 * with pixel values.
 *
 * Parameters:
 * @param  bck_mask     - the background mask
 * @param  smoot_length - number of pixels on either side to use for smoothing
 * @param  fib          - the background structure
 * @param  bck_mask     - the background mask
 * @param  ix           - x-value of central pixel
 * @param  iy           - y-value of central pixel
 * @param  pixvalues    - vector for pixel values
 * @param  pmask        - vector for pixel mask
 */
void
fill_pixvalues(const gsl_matrix *bck_mask, const int smooth_length,
               const fullimg_background *fib, const int ix, const int iy,
               gsl_vector *pixvalues, gsl_vector *pmask)
{
  int iact;
  int index;

  // initialize the pixel values
  // and the mask
  gsl_vector_set_all(pixvalues, 0.0);
  gsl_vector_set_all(pmask, 0.0);

  // iterate over the x-direction
  index=0;
  for (iact=ix-smooth_length; iact<=ix+smooth_length; iact++)
    {
      // check if you are within the chip
      if (iact > -1 && iact < (int)bck_mask->size1)
        {
          // check whether the pixel was interpolated
          // and whether the background is non-zero
          if (gsl_matrix_get(bck_mask, iact, iy) != 0.0
               && gsl_matrix_get(fib->bck, iact, iy) != 0.0)
            {
              // if yes, get the pixel value and set the mask
              gsl_vector_set(pixvalues, index, gsl_matrix_get(fib->bck, iact, iy));
              gsl_vector_set(pmask,     index, 1.0);
            }
        }
      // enhance the vector index
      index++;
    }
}

/**
 * Function: compute_gvalue
 * The function computes the values of a Gauss function
 * [exp(factor * xdiff^2)]. NO normalization factor is applied.
 *
 * Parameters:
 * @param xdiff   - the value [x-x_0]
 * @param efactor - the factor for the exponent
 *
 * Returns:
 * @return value - the Gaussian value
 */
double
compute_gvalue(const double xdiff, const double efactor)
{
  double value;

  // just compose the exp-function
  value = exp(efactor * xdiff * xdiff);

  // return the value
  return value;
}

/**
 * Function: compute_efactor
 * The function computes the appropriate factor of the Gaussian for any
 * FWHM given as input. This speeds up any later computation of the
 * Gaussian.
 *
 * Parameters:
 * @param fwhm - the input FWHM
 *
 * Returns:
 * @return expfactor - the factor for the Gaussian
 */
double
compute_efactor(const double fwhm)
{
  double sigma;
  double expfactor;

  // get the sigma value [sig = fwhm / (2 sqrt(2 ln(2)))]
  sigma = fwhm / 2.3458;

  // compute the factor for the Gaussian
  expfactor =  -1.0 / (2.0* sigma * sigma);

  // return the factor
  return expfactor;
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
#pragma once

#include "Common.h"
#include "PEFile.h"
#include "DelayLoadedModule.h"

#include <gsl/gsl>

namespace winlogo::pe_parser {

namespace details {
/**
 * This is an input iterator for iterating over the delay-loaded modules of the delay-load
 * directory.
 */
class DelayLoadDirectoryIterator {
public:
    DelayLoadDirectoryIterator(PEFile file,
                               gsl::span<IMAGE_DELAYLOAD_DESCRIPTOR>::iterator iterator);

    bool operator==(const DelayLoadDirectoryIterator& other) const noexcept;
    bool operator!=(const DelayLoadDirectoryIterator& other) const noexcept;

    DelayLoadedModule operator*() const noexcept;
    DelayLoadDirectoryIterator& operator++() noexcept;

private:
    PEFile m_peFile;
    gsl::span<IMAGE_DELAYLOAD_DESCRIPTOR>::iterator m_current;
};

}  // namespace details

/**
 * This class represents an iterable view of a PE file's delay-loaded directory.
 */
class DelayLoadDirectory {
public:
    /**
     * Initialize the delay-load directory from its PE file.
     */
    explicit DelayLoadDirectory(PEFile peFile);

    /**
     * Get the raw import descriptors.
     */
    gsl::span<IMAGE_DELAYLOAD_DESCRIPTOR> rawDelayLoadDescriptors() const noexcept;

    // Iterators for iterating over the imported modules.

    details::DelayLoadDirectoryIterator begin() const noexcept;
    details::DelayLoadDirectoryIterator end() const noexcept;

private:
    /// The owning PE file.
    PEFile m_peFile;
    gsl::span<IMAGE_DELAYLOAD_DESCRIPTOR> m_delayLoadDescriptors;
};

}  // namespace winlogo::pe_parser


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

#ifndef SIM_MATH_LINALG_H_
#define SIM_MATH_LINALG_H_

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

// We use GSL style conventions here to maintain consistency with the
// GSL library.  See GSL docs for comment.

#define GSL_VECTOR_STACK_ALLOC(v_, size_)                               \
  double v_##_data[size_];                                              \
  gsl_block v_##_block = {size_, v_##_data};                            \
  gsl_vector v_##_vector = {size_, 1, v_##_block.data, &v_##_block, 0}; \
  gsl_vector *v_ = &v_##_vector;

#define GSL_MATRIX_STACK_ALLOC(m_, size1_, size2_)                \
  double m_##_data[(size1_) * (size2_)];                          \
  gsl_block m_##_block = {(size1_) * (size2_), m_##_data};        \
  gsl_matrix m_##_matrix = {size1_,          size2_,      size2_, \
                            m_##_block.data, &m_##_block, 0};     \
  gsl_matrix *m_ = &m_##_matrix;

double gsl_matrix_trace(const gsl_matrix *m);
const gsl_vector *gsl_matrix_diag(const gsl_matrix *m, gsl_vector *d);
void gsl_matrix_disp(const gsl_matrix *m);
void gsl_vector_disp(const gsl_vector *v);
const gsl_vector *gsl_vector_saturate(const gsl_vector *x, double low,
                                      double high, gsl_vector *y);
double gsl_vector_norm_bound(const gsl_vector *x, double low);
double gsl_vector_dot(const gsl_vector *x, const gsl_vector *y);

#endif  // SIM_MATH_LINALG_H_


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

/***************************************************************************//**
 *
 * Version of zgemv for tile storage, to avoid dependency problem when
 * computations are done within the tile. alpha and beta are passed as
 * pointers so they can depend on runtime values.
 *
 * @param[in] Alock
 *          Pointer to tile owning submatrix A.
 *
 * @param[in] xlock
 *          Pointer to tile owning subvector x.
 *
 * @param[in] ylock
 *          Pointer to tile owning subvector y.
 *
 **/
void QUARK_CORE_cgemv_tile(Quark *quark, Quark_Task_Flags *task_flags,
                           PLASMA_enum trans,
                           int m, int n,
                           const PLASMA_Complex32_t *alpha, const PLASMA_Complex32_t *A, int lda,
                                                            const PLASMA_Complex32_t *x, int incx,
                           const PLASMA_Complex32_t *beta,        PLASMA_Complex32_t *y, int incy,
                           const PLASMA_Complex32_t *Alock,
                           const PLASMA_Complex32_t *xlock,
                           const PLASMA_Complex32_t *ylock)
{
    /* Quick return. Bad things happen if sizeof(...)*m*n is zero in QUARK_Insert_Task */
    if ( m == 0 || n == 0 )
        return;

    DAG_SET_PROPERTIES("gemv", "lightslateblue");
    QUARK_Insert_Task(quark, CORE_cgemv_tile_quark, task_flags,
        sizeof(PLASMA_enum),             &trans,  VALUE,
        sizeof(int),                     &m,      VALUE,
        sizeof(int),                     &n,      VALUE,
        sizeof(PLASMA_Complex32_t),      alpha,           INPUT,
        sizeof(PLASMA_Complex32_t)*m*n,  A,               NODEP,          /* input; see Alock */
        sizeof(int),                     &lda,    VALUE,
        sizeof(PLASMA_Complex32_t)*n,    x,               NODEP,          /* input; see xlock */
        sizeof(int),                     &incx,   VALUE,
        sizeof(PLASMA_Complex32_t),      beta,            INPUT,
        sizeof(PLASMA_Complex32_t)*m,    y,                       NODEP,  /* inout; see ylock */
        sizeof(int),                     &incy,   VALUE,
        sizeof(PLASMA_Complex32_t)*m*n,  Alock,           INPUT,
        sizeof(PLASMA_Complex32_t)*n,    xlock,           INPUT,
        sizeof(PLASMA_Complex32_t)*m,    ylock,                   INOUT,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cgemv_tile_quark = PCORE_cgemv_tile_quark
#define CORE_cgemv_tile_quark PCORE_cgemv_tile_quark
#endif
void CORE_cgemv_tile_quark(Quark *quark)
{
    PLASMA_enum trans;
    int m, n, lda, incx, incy;
    const PLASMA_Complex32_t *alpha, *beta;
    const PLASMA_Complex32_t *A, *x;
    PLASMA_Complex32_t *y;

    quark_unpack_args_11( quark, trans, m, n, alpha, A, lda, x, incx, beta, y, incy );
    cblas_cgemv(
        CblasColMajor,
        (CBLAS_TRANSPOSE)trans,
        m, n,
        CBLAS_SADDR(*alpha), A, lda,
                             x, incx,
        CBLAS_SADDR(*beta),  y, incy );
}


text:
// Copyright 2012 Jesse Windle - jesse.windle@gmail.com

// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.

// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.

// You should have received a copy of the GNU General Public License
// along with this program.  If not, see
// <http://www.gnu.org/licenses/>.

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

  This class wraps GSL's random number generator and random
  distribution functions into a class.  We use the Mersenne Twister
  for random number generation since it has a large period, which is
  what we want for MCMC simulation.

  When compiling include -lgsl -lcblas -llapack .

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

#ifndef __BASICRNG__
#define __BASICRNG__

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_sf.h>

#include <iostream>
#include <fstream>
#include <ctime>
#include <string>
#include <cmath>

using std::string;
using std::ofstream;
using std::ifstream;

//////////////////////////////////////////////////////////////////////
			      // RNG //
//////////////////////////////////////////////////////////////////////

class BasicRNG {

 protected:

  gsl_rng * r;

 public:

  // Constructors and destructors.
  BasicRNG(unsigned long seed);

  virtual ~BasicRNG()
    { gsl_rng_free (r); }

  // Assignment=
  BasicRNG& operator=(const BasicRNG& rng);

  // Read / Write / Set
  bool read (const string& filename);
  bool write(const string& filename);
  void set(unsigned long seed);

  // Get rng -- be careful.  Needed for other random variates.
  gsl_rng* getrng() { return r; }

  // Random variates.
  double unif  ();                             // Uniform
  double expon_mean(double mean);     // Exponential
  double expon_rate(double rate);                  // Exponential
  double chisq (double df);                    // Chisq
  double norm  (double sd);                    // Normal
  double norm  (double mean , double sd);      // Normal
  double gamma_scale (double shape, double scale); // Gamma_Scale
  double gamma_rate  (double shape, double rate);  // Gamma_Rate
  double igamma(double shape, double scale);   // Inv-Gamma
  double flat  (double a=0  , double b=1  );   // Flat
  double beta  (double a=1.0, double b=1.0);   // Beta

  int bern  (double p);                     // Bernoulli

  // CDF
  static double p_norm (double x, int use_log=0);
  static double p_gamma_rate(double x, double shape, double rate, int use_log=0);

  // Density
  static double d_beta(double x, double a, double b);

  // Utility
  static double Gamma (double x, int use_log=0);

}; // BasicRNG

#endif

////////////////////////////////////////////////////////////////////////////////
				 // APPENDIX //
////////////////////////////////////////////////////////////////////////////////

// If you make everything inline within the same translation unit then that
// function will not be callable from anohter translation unit.  You can see
// that the function is missing by using the nm command.

max_stars_repo_name: TeoGiane/SPMIX

max_issues_repo_name: TeoGiane/SPMIX

max_forks_repo_name: TeoGiane/SPMIX


text:
#include <pygsl/solver.h>
#include <gsl/gsl_roots.h>

const char  * filename = __FILE__;
PyObject *module = NULL;

static const char root_f_type_name[] = "F-RootSolver";
static const char root_fdf_type_name[] = "FdF-RootSolver";
static const char root_f_root_doc[] = "Get the value of F";
static const char root_set_f_doc[] = "";
static const char root_set_fdf_doc[] = "";
static const char root_fdf_root_doc[] = "";
static const char root_x_lower_doc  [] = "Get the lower bound of x";
static const char root_x_upper_doc  [] = "Get the upper bound of x";

static PyObject*
PyGSL_root_f_root(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_root_fsolver_root);
}

static PyObject*
PyGSL_root_fdf_root(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_root_fdfsolver_root);
}

static PyObject*
PyGSL_root_x_lower(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_root_fsolver_x_lower);
}

static PyObject*
PyGSL_root_x_upper(PyGSL_solver *self, PyObject *args)
{
     return PyGSL_solver_ret_double(self, args, (double_m_t) gsl_root_fsolver_x_upper);
}

static PyObject*
PyGSL_root_solver_test_interval(PyGSL_solver * self, PyObject *args)
{
     double epsabs, epsrel;
     gsl_root_fsolver *s = (gsl_root_fsolver *) self->solver;
     if(!PyArg_ParseTuple(args, "dd", &epsabs, &epsrel))
	  return NULL;
     return PyInt_FromLong(gsl_root_test_interval(s->x_lower, s->x_upper, epsabs, epsrel));
}

static PyObject*
PyGSL_root_set_f(PyGSL_solver *self, PyObject *args, PyObject *kw)
{
     return PyGSL_solver_set_f(self, args, kw, (void *)gsl_root_fsolver_set, 0);
}

static PyObject*
PyGSL_root_set_fdf(PyGSL_solver *self, PyObject *args, PyObject *kw)
{
     return PyGSL_solver_set_f(self, args, kw,  (void *)gsl_root_fdfsolver_set, 1);
}

static PyMethodDef PyGSL_root_fmethods[] = {
     {"root",      (PyCFunction)PyGSL_root_f_root,     METH_NOARGS, (char *)root_f_root_doc},
     {"x_lower",   (PyCFunction)PyGSL_root_x_lower,    METH_NOARGS, (char *)root_x_lower_doc},
     {"x_upper",   (PyCFunction)PyGSL_root_x_upper,    METH_NOARGS, (char *)root_x_upper_doc},
     {"set",       (PyCFunction)PyGSL_root_set_f,      METH_VARARGS|METH_KEYWORDS, (char *)root_set_f_doc},
     {"test_interval",(PyCFunction)PyGSL_root_solver_test_interval, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}           /* sentinel */
};
static PyMethodDef PyGSL_root_fdfmethods[] = {
     {"root",      (PyCFunction)PyGSL_root_fdf_root,     METH_NOARGS, (char *)root_fdf_root_doc},
     {"set",       (PyCFunction)PyGSL_root_set_fdf,      METH_VARARGS|METH_KEYWORDS, (char *)root_set_fdf_doc},
     {NULL, NULL, 0, NULL}           /* sentinel */
};

const struct _SolverStatic
root_solver_f = {{(void_m_t) gsl_root_fsolver_free,
		  /* gsl_multimin_fminimizer_restart */  (void_m_t) NULL,
		  (name_m_t) gsl_root_fsolver_name,
		  (int_m_t) gsl_root_fsolver_iterate},
		 1, PyGSL_root_fmethods, root_f_type_name},
root_solver_fdf = {{(void_m_t) gsl_root_fdfsolver_free,
		    /* gsl_multimin_fminimizer_restart */  (void_m_t) NULL,
		    (name_m_t) gsl_root_fdfsolver_name,
		     (int_m_t) gsl_root_fdfsolver_iterate},
		   3, PyGSL_root_fdfmethods, root_fdf_type_name};

static PyObject*
PyGSL_root_f_init(PyObject *self, PyObject *args,
		 const gsl_root_fsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_root_fsolver_alloc,
			      &root_solver_f};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 0);
     FUNC_MESS_END();
     return tmp;
}

static PyObject*
PyGSL_root_fdf_init(PyObject *self, PyObject *args,
		 const gsl_root_fdfsolver_type * type)
{

     PyObject *tmp=NULL;
     solver_alloc_struct s = {type, (void_an_t) gsl_root_fdfsolver_alloc,
			      &root_solver_fdf};
     FUNC_MESS_BEGIN();
     tmp = PyGSL_solver_dn_init(self, args, &s, 0);
     FUNC_MESS_END();
     return tmp;
}

#define AROOT_F(name)                                                  \
static PyObject* PyGSL_root_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_root_f_init(self, args,  gsl_root_fsolver_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}
#define AROOT_FDF(name)                                                  \
static PyObject* PyGSL_root_init_ ## name (PyObject *self, PyObject *args)\
{                                                                             \
     PyObject *tmp = NULL;                                                    \
     FUNC_MESS_BEGIN();                                                       \
     tmp = PyGSL_root_fdf_init(self, args,  gsl_root_fdfsolver_ ## name); \
     if (tmp == NULL){                                                        \
	  PyGSL_add_traceback(module, __FILE__, __FUNCTION__, __LINE__); \
     }                                                                        \
     FUNC_MESS_END();                                                         \
     return tmp;                                                              \
}
AROOT_F(bisection)
AROOT_F(falsepos)
AROOT_F(brent)
AROOT_FDF(newton)
AROOT_FDF(secant)
AROOT_FDF(steffenson)

static PyObject*
PyGSL_root_test_delta(PyObject * self, PyObject *args)
{
     double x_lower, x_upper, epsabs, epsrel;
     if(!PyArg_ParseTuple(args, "dddd", &x_lower, &x_upper, &epsabs, &epsrel))
	  return NULL;
     return PyInt_FromLong(gsl_root_test_delta(x_lower, x_upper, epsabs, epsrel));
}

static PyObject*
PyGSL_root_test_interval(PyObject * self, PyObject *args)
{
     double x_lower, x_upper, epsabs, epsrel;
     if(!PyArg_ParseTuple(args, "dddd", &x_lower, &x_upper, &epsabs, &epsrel))
	  return NULL;
     return PyInt_FromLong(gsl_root_test_interval(x_lower, x_upper, epsabs, epsrel));
}

static const char PyGSL_roots_module_doc [] = "XXX Missing ";
static PyMethodDef mMethods[] = {
     /* solver */
     {"bisection",  PyGSL_root_init_bisection, METH_NOARGS, NULL},
     {"falsepos",  PyGSL_root_init_falsepos, METH_NOARGS, NULL},
     {"brent",  PyGSL_root_init_brent, METH_NOARGS, NULL},
     {"newton",  PyGSL_root_init_newton, METH_NOARGS, NULL},
     {"secant",  PyGSL_root_init_secant, METH_NOARGS, NULL},
     {"steffenson",  PyGSL_root_init_steffenson, METH_NOARGS, NULL},
     /* functions */
     {"test_delta",  PyGSL_root_test_delta, METH_VARARGS, NULL},
     {"test_interval",  PyGSL_root_test_interval, METH_VARARGS, NULL},
     {NULL, NULL, 0, NULL}
};

void
initroots(void)
{
     PyObject* m, *dict, *item;
     FUNC_MESS_BEGIN();

     m=Py_InitModule("roots", mMethods);
     module = m;
     assert(m);
     dict = PyModule_GetDict(m);
     if(!dict)
	  goto fail;

     init_pygsl()
     import_pygsl_solver();
     assert(PyGSL_API);

     if (!(item = PyString_FromString((char*)PyGSL_roots_module_doc))){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not generate module doc string!");
	  goto fail;
     }

     if (PyDict_SetItemString(dict, "__doc__", item) != 0){
	  PyErr_SetString(PyExc_ImportError,
			  "I could not init doc string!");
	  goto fail;
     }

     FUNC_MESS_END();

 fail:
     FUNC_MESS("FAIL");
     return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/******************************************************************************

  (c) 2017 - 2019 Scientific Computation Research Center,
      Rensselaer Polytechnic Institute. All rights reserved.

  This work is open source software, licensed under the terms of the
  BSD license as described in the LICENSE file in the top-level directory.

*******************************************************************************/
#ifndef MSI_PETSC_H
#define MSI_PETSC_H
#include "msi.h"
#include "msi_solver.h"
#include <pumi.h>
#include <petsc.h>
#include <map>
#include <vector>
// Added for the synchronization function
// TODO : make it so these are not necessary
#include "apfFieldData.h"
#include "apfNew.h"
#include "apfNumbering.h"
#include "apfNumberingClass.h"
#include "apfShape.h"
int copyField2PetscVec(pField f, Vec& petscVec);
int copyPetscVec2Field(Vec& petscVec, pField f, bool);
void printMemStat( );
// NOTE: all field related interaction is done through msi api rather than apf
class msi_matrix
{
  public:
  msi_matrix(pField f);
  virtual ~msi_matrix( );
  virtual int initialize( ) = 0;  // create a matrix and solver object
  int destroy( );                 // delete a matrix and solver object
  int set_value(msi_int row,
                msi_int col,
                int operation,
                double real_val,
                double imag_val);  // insertion/addition with global numbering
  // values use column-wise, size * size block
  int add_values(msi_int rsize,
                 msi_int* rows,
                 msi_int csize,
                 msi_int* columns,
                 double* values);
  int get_values(std::vector<msi_int>& rows,
                 std::vector<msi_int>& n_columns,
                 std::vector<msi_int>& columns,
                 std::vector<double>& values);
  void set_status(int s) { mat_status = s; }
  int get_status( ) { return mat_status; }
  int get_scalar_type( )
  {
#ifdef MSI_COMPLEX
    return 1;
#else
    return 0;
#endif
  }
  pField get_field( ) { return field; }
  int write(const char* file_name);
  virtual int get_type( ) const = 0;
  virtual int assemble( ) = 0;
  virtual int setupMat( ) = 0;
  virtual int preAllocate( ) = 0;
  virtual int flushAssembly( );
  int printInfo( );
  // PETSc data structures
  Mat* A;
  protected:
  int setupSeqMat( );
  int setupParaMat( );
  int preAllocateSeqMat( );
  int preAllocateParaMat( );
  int mat_status;
  pField field;  // the field that provide numbering
};
class matrix_mult : public msi_matrix
{
  public:
  matrix_mult(pField f) : msi_matrix(f), localMat(1) { initialize( ); }
  virtual int initialize( );
  void set_mat_local(bool flag) { localMat = flag; }
  int is_mat_local( ) { return localMat; }
  int multiply(pField in_f, pField out_f, bool sync = true);
  virtual int get_type( ) const { return 0; }  // MSI_MULTIPLY; }
  virtual int assemble( );
  virtual int setupMat( );
  virtual int preAllocate( );
  private:
  bool localMat;
};
class matrix_solve : public msi_matrix
{
  public:
  matrix_solve(pField f);
  virtual int initialize( );
  virtual ~matrix_solve( );
  int solve(pField rhs, pField sol, bool sync = true);
  int set_bc(msi_int row);
  int set_row(msi_int row, msi_int numVals, msi_int* colums, double* vals);
  int add_blockvalues(msi_int rbsize,
                      msi_int* rows,
                      msi_int cbsize,
                      msi_int* columns,
                      double* values);
  virtual int get_type( ) const { return 1; }
  virtual int assemble( );
  virtual int setupMat( );
  virtual int preAllocate( );
  int iterNum;
  private:
  int setUpRemoteAStruct( );
  int setUpRemoteAStructParaMat( );
  int setKspType( );
  int kspSet;
  KSP* ksp;
  Mat remoteA;
  std::set<int> remotePidOwned;
  std::map<int, std::map<int, int> >
    remoteNodeRow;  // <pid, <locnode>, numAdj >
  std::map<int, int> remoteNodeRowSize;
};
#endif

max_stars_repo_name: SCOREC/msi

max_issues_repo_name: SCOREC/msi

max_forks_repo_name: SCOREC/msi


text:
#include <sys/time.h>
#include <cblas.h>
#include "kadbg.h"
#include "matrix.h"
#include "kalas.h"

#define MAX_BUF_SIZE 128*1024*1024
#define NUM_OF_STEPS 1024

void test_sgemm( KalasState *state)
{
	int rowA, colA, rowB, colB, rowC, colC;
	Matrix *matA = NULL, *matB = NULL, *matC1 = NULL, *matC2 = NULL;
	float alpha = 1.0, beta = 0.0, e;
  struct timeval t;
  double ts, te;

	for ( rowA = NUM_OF_STEPS; rowA <= MAX_BUF_SIZE / NUM_OF_STEPS
					/ sizeof(float);
				rowA += NUM_OF_STEPS) {
		for ( colA = NUM_OF_STEPS; rowA * colA <= MAX_BUF_SIZE / sizeof(float);
					colA += NUM_OF_STEPS) {

			rowB = colA;
			rowC = rowA;
			for ( colB = NUM_OF_STEPS; rowB * colB <= MAX_BUF_SIZE / sizeof(float);
						colB += NUM_OF_STEPS) {
				colC = colB;

				if ( rowC * colC > MAX_BUF_SIZE / sizeof(float)) break;

				if ( IS_FAILED( ( matA = matrixNew( rowA, colA, colA,
																						sizeof(float), true)) != NULL))
					break;
				if ( IS_FAILED( ( matB = matrixNew( rowB, colB, colB,
																						sizeof(float), true)) != NULL))
					break;
				if ( IS_FAILED( ( matC1 = matrixNew( rowC, colC, colC,
																						 sizeof(float), true)) != NULL))
					break;
				if ( IS_FAILED( ( matC2 = matrixNew( rowC, colC, colC,
																						 sizeof(float), true)) != NULL))
					break;

				printf( "(%d,%d)(%d,%d)=(%d,%d) ",
								rowA, colA, rowB, colB, rowC, colC);

				gettimeofday( &t, NULL);
				ts = t.tv_sec + (double)(t.tv_usec) / 1000000.0;

				cblas_sgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
										 matC1->row, matC1->col, matA->col,
										 alpha, matA->elm, matA->col,
										 matB->elm, matB->col, beta,
										 matC1->elm, matC1->col);

				gettimeofday( &t, NULL);
				te = t.tv_sec + (double)(t.tv_usec) / 1000000.0;
				printf( "%f s %f GFlops ", (te - ts),
								(double)rowA * (double)colA * (double)colB * 2.0
								/ ( te - ts) / 1000000000.0);

				gettimeofday( &t, NULL);
				ts = t.tv_sec + (double)(t.tv_usec) / 1000000.0;

				kalasSgemm( state, clblasRowMajor, clblasNoTrans, clblasNoTrans,
										matC2->row, matC2->col, matA->col,
										alpha, matA->elm, matA->ld,
										matB->elm, matB->ld, beta,
										matC2->elm, matC2->ld);

				gettimeofday( &t, NULL);
				te = t.tv_sec + (double)(t.tv_usec) / 1000000.0;
				printf( "%f s %f GFlops ", (te - ts),
								(double)rowA * (double)colA * (double)colB * 2.0
								/ ( te - ts) / 1000000000.0);

				printf( "%e err", matrixCalcDiff( matC1, matC2));
				putchar( '\n');
				matrixDelete( matA); matA = NULL;
				matrixDelete( matB); matB = NULL;
				matrixDelete( matC1); matC1 = NULL;
				matrixDelete( matC2); matC2 = NULL;
			}
		}
  }
}

void test_dgemm( KalasState *state)
{
	int rowA, colA, rowB, colB, rowC, colC;
	Matrix *matA = NULL, *matB = NULL, *matC1 = NULL, *matC2 = NULL;
	double alpha = 1.0, beta = 0.0, e;
  struct timeval t;
  double ts, te;

	for ( rowA = NUM_OF_STEPS; rowA <= MAX_BUF_SIZE / NUM_OF_STEPS
					/ sizeof(double);
				rowA += NUM_OF_STEPS) {
		for ( colA = NUM_OF_STEPS; rowA * colA <= MAX_BUF_SIZE / sizeof(double);
					colA += NUM_OF_STEPS) {

			rowB = colA;
			rowC = rowA;
			for ( colB = NUM_OF_STEPS; rowB * colB <= MAX_BUF_SIZE / sizeof(double);
						colB += NUM_OF_STEPS) {
				colC = colB;

				if ( rowC * colC > MAX_BUF_SIZE / sizeof(double)) break;

				if ( IS_FAILED( ( matA = matrixNew( rowA, colA, colA,
																						sizeof(double), true)) != NULL))
					break;
				if ( IS_FAILED( ( matB = matrixNew( rowB, colB, colB,
																						sizeof(double), true)) != NULL))
					break;
				if ( IS_FAILED( ( matC1 = matrixNew( rowC, colC, colC,
																						 sizeof(double), true)) != NULL))
					break;
				if ( IS_FAILED( ( matC2 = matrixNew( rowC, colC, colC,
																						 sizeof(double), true)) != NULL))
					break;

				printf( "(%d,%d)(%d,%d)=(%d,%d) ",
								rowA, colA, rowB, colB, rowC, colC);

				gettimeofday( &t, NULL);
				ts = t.tv_sec + (double)(t.tv_usec) / 1000000.0;

				cblas_dgemm( CblasRowMajor, CblasNoTrans, CblasNoTrans,
										 matC1->row, matC1->col, matA->col,
										 alpha, matA->elm, matA->col,
										 matB->elm, matB->col, beta,
										 matC1->elm, matC1->col);

				gettimeofday( &t, NULL);
				te = t.tv_sec + (double)(t.tv_usec) / 1000000.0;
				printf( "%f s %f GFlops ", (te - ts),
								(double)rowA * (double)colA * (double)colB * 2.0
								/ ( te - ts) / 1000000000.0);

				gettimeofday( &t, NULL);
				ts = t.tv_sec + (double)(t.tv_usec) / 1000000.0;

				kalasDgemm( state, clblasRowMajor, clblasNoTrans, clblasNoTrans,
										matC2->row, matC2->col, matA->col,
										alpha, matA->elm, matA->ld,
										matB->elm, matB->ld, beta,
										matC2->elm, matC2->ld);

				gettimeofday( &t, NULL);
				te = t.tv_sec + (double)(t.tv_usec) / 1000000.0;
				printf( "%f s %f GFlops ", (te - ts),
								(double)rowA * (double)colA * (double)colB * 2.0
								/ ( te - ts) / 1000000000.0);

				printf( "%e err", matrixCalcDiff( matC1, matC2));
				putchar( '\n');
				matrixDelete( matA); matA = NULL;
				matrixDelete( matB); matB = NULL;
				matrixDelete( matC1); matC1 = NULL;
				matrixDelete( matC2); matC2 = NULL;
			}
		}
  }
}

int main( void)
{
	KalasState *state;

	state = kalasStateNew();

	test_sgemm( state);
	test_dgemm( state);

	kalasStateDelete( state);
	return 0;
}


text:
/// \file linalg.h
/// \brief Wrapper for including CBLAS and LAPACKE

#ifndef LINALG_H
#define LINALG_H

#define LAPACK_DISABLE_NAN_CHECK

#ifdef USE_MKL

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

#else

// generic include
#include <cblas.h>
#include <lapacke.h>

#endif

#endif


text:
/* this program is in the public domain */

#include "bench-user.h"
#include <math.h>

#ifdef HAVE_GSL_GSL_VERSION_H
#include <gsl/gsl_version.h>
#endif
static const char *mkvers(void)
{
#ifdef HAVE_GSL_GSL_VERSION_H
     return gsl_version;
#else
     return "unknown"
#endif
}

BEGIN_BENCH_DOC
BENCH_DOC("name", NAME)
BENCH_DOC("package", "GNU Scientific Library (GSL)")
BENCH_DOC("author", "Brian Gough")
BENCH_DOC("year", "2002")
BENCH_DOC("url", "http://sources.redhat.com/gsl")
BENCH_DOC("url-was-valid-on", "Thu Jul 26 07:48:45 EDT 2001")
BENCH_DOC("copyright",
  "Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough\n"
  "\n"
  "This program is free software; you can redistribute it and/or modify\n"
  "it under the terms of the GNU General Public License as published by\n"
  "the Free Software Foundation; either version 2 of the License, or (at\n"
  "your option) any later version.\n")
BENCH_DOC("language", "C")
BENCH_DOCF("version", mkvers)
END_BENCH_DOC

#define CAT0(a, b) a ## b
#define CAT(a, b) CAT0(a, b)

#ifdef BENCHFFT_SINGLE
#include <gsl/gsl_fft_complex_float.h>
#include <gsl/gsl_fft_real_float.h>
#include <gsl/gsl_fft_halfcomplex_float.h>
#define C(x) CAT(gsl_fft_complex_float_, P(x))
#define R(x) CAT(gsl_fft_real_float_radix2_, x)
#define H(x) CAT(gsl_fft_halfcomplex_float_radix2_, x)
#else
#include <gsl/gsl_fft_complex.h>
#include <gsl/gsl_fft_real.h>
#include <gsl/gsl_fft_halfcomplex.h>
#define C(x) CAT(gsl_fft_complex_, P(x))
#define R(x) CAT(gsl_fft_real_radix2_, x)
#define H(x) CAT(gsl_fft_halfcomplex_radix2_, x)
#endif

int can_do(struct problem *p)
{
     return (p->rank == 1 &&
#ifdef DO_COMPLEX
	     p->kind == PROBLEM_COMPLEX
#else
	     p->kind == PROBLEM_REAL
#endif
	     && problem_power_of_two(p, 1));
}

void copy_h2c(struct problem *p, bench_complex *out)
{
     copy_h2c_1d_halfcomplex(p, out, -1.0);
}

void copy_c2h(struct problem *p, bench_complex *in)
{
     copy_c2h_1d_halfcomplex(p, in, -1.0);
}

void setup(struct problem *p)
{
     BENCH_ASSERT(can_do(p));
     /* nothing to do */
}

void doit(int iter, struct problem *p)
{
     int i;
     int n = p->n[0];
     void *in = p->in;

     if (p->kind == PROBLEM_COMPLEX) {
#ifdef DO_COMPLEX
	  if (p->sign == -1) {
	    for (i = 0; i < iter; ++i)
	      C(forward)(in, 1, n);
	  } else {
	    for (i = 0; i < iter; ++i)
	      C(backward)(in, 1, n);
	  }
#endif
     } else {
#ifdef DO_REAL
	  if (p->sign == -1) {
	    for (i = 0; i < iter; ++i)
	      R(transform)(in, 1, n);
	  } else {
	    for (i = 0; i < iter; ++i)
	      H(transform)(in, 1, n);
	  }
#endif
     }
}

void done(struct problem *p)
{
     UNUSED(p);
}

max_stars_repo_name: mreineck/benchfft

max_issues_repo_name: mreineck/benchfft

max_forks_repo_name: mreineck/benchfft