text:
/*
 * SpectrumTransformer.h
 *
 * Created on: Jul 30, 2015
 *     Author: dpayne
 */

#ifndef _SPECTRUM_TRANSFORMER_H
#define _SPECTRUM_TRANSFORMER_H

#include <fftw.h>
#include "Transformer/GenericTransformer.h"
#include "Writer/NcursesWriter.h"
#include "Domain/Settings.h"

namespace vis
{

class SpectrumTransformer : public GenericTransformer
{
  public:
    explicit SpectrumTransformer(const Settings *const settings);

    ~SpectrumTransformer() override;

    void execute_mono(pcm_stereo_sample *buffer,
                      vis::NcursesWriter *writer) override;
    void execute_stereo(pcm_stereo_sample *buffer,
                        vis::NcursesWriter *writer) override;

  private:
    void execute(pcm_stereo_sample *buffer, vis::NcursesWriter *writer,
                 const bool is_stereo);
    const Settings *const m_settings;

    /** --- BEGIN MEMBER VARIABLES --- */

    /** --- BEGIN --- fft calculations vars **/
    size_t m_fftw_results;

    double *m_fftw_input_left;
    double *m_fftw_input_right;

    fftw_complex *m_fftw_output_left;
    fftw_complex *m_fftw_output_right;

    fftw_plan m_fftw_plan_left;
    fftw_plan m_fftw_plan_right;
    /** --- END --- fft calculations vars **/

    /** --- BEGIN --- frequency cutoff calculations vars **/
    std::vector<uint32_t> m_low_cutoff_frequencies;

    std::vector<uint32_t> m_high_cutoff_frequencies;

    std::vector<double> m_frequency_constants_per_bin;

    int32_t m_previous_win_width; // Used to determine if freq cutoffs need to
                                  // be re-calculated
    /** --- END --- frequency cutoff calculations vars **/

    uint64_t m_silent_runs; // Used to determine if the transformer should sleep
                            // and wait for input

    // Holds the current bar heights after processing, this is done as a member
    // function to avoid memory allocations on every run
    std::vector<double> m_bars_left;
    std::vector<double> m_bars_right;

    // falloff vectors hold the previous runs values, this is used to to apply
    // falloff effect for the current run
    std::vector<double> m_bars_falloff_left;
    std::vector<double> m_bars_falloff_right;

    std::vector<double> m_previous_max_heights;

    // Used by monstercat smoothing to apply weights to a bar's height as
    // determined by it's frequency range
    // Note: this is only re-computed when screen width changes
    std::vector<double> m_monstercat_smoothing_weights;

    // Pre-compute colors calculations to avoid duplicate work
    // Note: this is only re-computed when screen height changes
    std::vector<vis::ColorIndex> m_precomputed_colors;

    /** --- END MEMBER VARIABLES --- */

    /** --- BEGIN MEMBER FUNCTIONS --- */

    /**
     * Copies the channel given by "channel_mode" to the fftw_input buffer
     */
    bool prepare_fft_input(pcm_stereo_sample *buffer, uint32_t sample_size,
                           double *fftw_input, ChannelMode channel_mode);

    /**
     * Populates "bars" and "bars_falloff" with the bar heights to be displayed
     */
    virtual void create_spectrum_bars(fftw_complex *fftw_output,
                                      const size_t fftw_results,
                                      const int32_t win_height,
                                      const int32_t win_width,
                                      const uint32_t number_of_bars,
                                      std::vector<double> &bars,
                                      std::vector<double> &bars_falloff);

    void
    generate_bars(std::vector<double> &bars, const uint32_t number_of_bars,
                  const fftw_complex *fftw_output, const size_t fftw_results,
                  const std::vector<uint32_t> &low_cutoff_frequencies,
                  const std::vector<uint32_t> &high_cutoff_frequencies) const;

    void recalculate_cutoff_frequencies(
        uint32_t number_of_bars, std::vector<uint32_t> *low_cutoff_frequencies,
        std::vector<uint32_t> *high_cutoff_frequencies,
        std::vector<double> *freqconst_per_bin);

    /**
     * Applies the smoothing operations based on the settings in m_settings to
     * "bars"
     */
    void smooth_bars(std::vector<double> &bars);

    /**
     * Applies the falloff effect based on the settings in m_settings to "bars".
     * The old falloff from the previous run should be in "falloff_bars". The
     * new falloff values will be updated in-place in "falloff_bars"
     */
    std::vector<double> apply_falloff(const std::vector<double> &bars,
                                      std::vector<double> &falloff_bars) const;

    /**
     * Calculates the moving average and the standard deviation for a given
     * range of elements.
     *
     * "new_value" is added to "old_values".
     *
     * If "old_values.size()" >  "max_number_of_elements" the first element of
     * old_values is erased.
     */
    void calculate_moving_average_and_std_dev(
        const double new_value, const size_t max_number_of_elements,
        std::vector<double> &old_values, double *moving_average,
        double *std_dev) const;

    /**
     * A long term and short term running average are kept of the max bar
     * heights for each frame. If the short term running average is very
     * different than the long term running average then the scaling size
     * should be reset to better suit the current music.
     *
     * This happens commonly if a there is a new song that is a lot quieter
     * or louder than the previous song.
     */
    void maybe_reset_scaling_window(const double current_max_height,
                                    const size_t max_number_of_elements,
                                    std::vector<double> *values,
                                    double *moving_average, double *std_dev);

    /**
     * Renders the spectrum bars to the screen.
     */
    virtual void draw_bars(const std::vector<double> &bars,
                           const std::vector<double> &bars_falloff,
                           int32_t win_height, const bool flipped,
                           const std::wstring &bar_row_msg,
                           vis::NcursesWriter *writer);

    /**
     * Scaling the given vector of points "bars" to a fit a screen with a window
     * height of "height".
     */
    virtual void scale_bars(std::vector<double> &bars, const int32_t height);

    /**
     * Creates the to be used for every section of the bar to be printed. For
     * example if the bar width is set to two, and the bar character to '#'.
     * Then the bar row msg would be "##";
     */
    std::wstring create_bar_row_msg(const wchar_t character,
                                    uint32_t bar_width);

    /**
     * Savitzky-Golay smoothng. This type of smoothing is usually much faster
     * than monstercat.
     *
     *  https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
     */
    void sgs_smoothing(std::vector<double> &bars);

    /**
     * This type of smoothing is inspired by monstercat
     * (https://www.youtube.com/user/MonstercatMedia). The code was largely
     * taken from cava git@github.com:karlstav/cava.git
     */
    void monstercat_smoothing(std::vector<double> &bars);

    /** --- END MEMBER FUNCTIONS --- */
};
}

#endif

max_stars_repo_name: kevinarefunny/cli-visualizer

max_issues_repo_name: kevinarefunny/cli-visualizer

max_forks_repo_name: kevinarefunny/cli-visualizer


text:
/* ensure once-only inclusion. */
#ifndef __VFL_BLAS_H__
#define __VFL_BLAS_H__

/* include the matrix and vector headers. */
#include <vfl/util/matrix.h>
#include <vfl/util/vector.h>

/* determine whether or not atlas is used. */
#ifdef __VFL_USE_ATLAS
/* include the atlas blas header. */
#include <cblas.h>
#include <clapack.h>
#else
/* define values for cblas enumerations. */
#define CblasNoTrans 111
#define CblasTrans   112
#define CblasUpper   121
#define CblasLower   122
#endif

/* BlasTranspose: enumeration of all possible ways to transpose
 * a matrix (or not) during calculations.
 */
typedef enum {
  BLAS_NO_TRANS = CblasNoTrans,
  BLAS_TRANS    = CblasTrans
}
BlasTranspose;

/* BlasTriangle: enumeration of all possible ways to access
 * a triangular matrix.
 */
typedef enum {
  BLAS_UPPER = CblasUpper,
  BLAS_LOWER = CblasLower
}
BlasTriangle;

/* function declarations (util/blas.c): */

double blas_dasum (const Vector *x);

double blas_dnrm2 (const Vector *x);

double blas_ddot (const Vector *x, const Vector *y);

void blas_daxpy (double alpha, const Vector *x, Vector *y);

void blas_dscal (double alpha, Vector *y);

/* --- */

void blas_dgemv (BlasTranspose trans, double alpha, const Matrix *A,
                 const Vector *x, double beta, Vector *y);

void blas_dtrmv (BlasTranspose trans, const Matrix *L,
                 const Vector *x, Vector *y);

void blas_dtrsv (BlasTriangle tri, const Matrix *L, Vector *x);

#endif /* !__VFL_BLAS_H__ */

ext: h

max_stars_repo_name: geekysuavo/vfl

max_issues_repo_name: geekysuavo/vfl

max_forks_repo_name: geekysuavo/vfl


text:
// @(#)root/matrix:$Name:  $:$Id: TMatrixTSparse.h,v 1.5 2006/05/17 06:22:06 brun Exp $
// Authors: Fons Rademakers, Eddy Offermann   Feb 2004

/*************************************************************************
 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

//#ifndef ROOT_TMatrixTSparse
//#define ROOT_TMatrixTSparse

//#ifndef ROOT_TMatrixTBase
#include "TMatrixTBase.cpp" //KEVIN .... program wont' link correctedly when using templates unless i include cxx file instead of .h
#include "Rtypes.h"   //KEVIN ADDED THIS
//#endif
//#ifndef ROOT_TMatrixTUtils
//#include "TMatrixTUtils.h"
//#endif

//#ifdef CBLAS
//#include <vecLib/vBLAS.h>
//#include <cblas.h>
//#endif

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TMatrixTSparse                                                       //
//                                                                      //
// Template class of a general sparse matrix in the Harwell-Boeing      //
// format                                                               //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

//template<class Element> class TMatrixT;

template<class Element> class TMatrixTSparse : public TMatrixTBase<Element> {

protected:

   Int_t   *fRowIndex;  //[fNrowIndex] row index
   Int_t   *fColIndex;  //[fNelems]    column index
   Element *fElements;  //[fNelems]

   void Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,
                 Int_t init = 0,Int_t nr_nonzeros = 0);

  // Elementary constructors
   //void AMultB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
   //             const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,b); AMultBt(a,bt,constr); }
   //void AMultB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0) {
   //             const TMatrixTSparse<Element> bsp = b;
   //             const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,bsp); AMultBt(a,bt,constr); }
   //void AMultB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
   //             const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,b); AMultBt(a,bt,constr); }

   //void AMultBt(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
   //void AMultBt(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);

   //void APlusB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
   //void APlusB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
   //void APlusB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) { APlusB(b,a,constr); }

   //void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
   //void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
   //void AMinusB(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);

public:

   enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kAtA }; //IN ROSA'S CODE THESE TWO LINES ARE IN TMATRIXDBASE.H
   enum EMatrixCreatorsOp2 { kMult,kMultTranspose,kPlus,kMinus };

   void AMultBt(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=1); //CHANGED CONSTR=0 TO CONSTR = 1 FROM ROSA'S CODE
   //IN ADDITION THIS WAS ORIGINALLY PRIVATE... MOVED IT TO PUBLIC FOR TESTTMATRIXDSPARSE.CPP

   TMatrixTSparse() { fElements = 0; fRowIndex = 0; fColIndex = 0; }
   TMatrixTSparse(Int_t nrows,Int_t ncols);
   //TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
   TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
                  Int_t *row, Int_t *col,Element *data);
   //TMatrixTSparse(const TMatrixTSparse<Element> &another);
   //TMatrixTSparse(const TMatrixT<Element>       &another);

   //TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype);
   //TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);

   //virtual ~TMatrixTSparse() { Clear(); }

   virtual const Element *GetMatrixArray  () const;
   virtual       Element *GetMatrixArray  ();
   virtual const Int_t    *GetRowIndexArray() const;
   virtual       Int_t    *GetRowIndexArray();
   virtual const Int_t    *GetColIndexArray() const;
   virtual       Int_t    *GetColIndexArray();

   //virtual TMatrixTBase<Element>   &SetRowIndexArray(Int_t *data) { memmove(fRowIndex,data,(this->fNrows+1)*sizeof(Int_t)); return *this; }
   //virtual TMatrixTBase<Element>   &SetColIndexArray(Int_t *data) { memmove(fColIndex,data,this->fNelems*sizeof(Int_t)); return *this; }

   //virtual void                     GetMatrix2Array (Element *data,Option_t *option="") const;
   //virtual TMatrixTBase<Element>   &SetMatrixArray  (const Element *data,Option_t * /*option*/="")
   //                                                 { memcpy(fElements,data,this->fNelems*sizeof(Element)); return *this; }

   //KEVIN: SetMatrixArray used to be virtual..... gives me vtable errors.... don't know what it is though... googled it
   TMatrixTBase<Element>   &SetMatrixArray  (Int_t nr_nonzeros,Int_t *irow,Int_t *icol,Element *data);
           TMatrixTSparse<Element> &SetSparseIndex  (Int_t nelem_new);
	   //        TMatrixTSparse<Element> &SetSparseIndex  (const TMatrixTBase<Element> &another);
           //TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b);
	   //virtual TMatrixTBase<Element>   &InsertRow       (Int_t row,Int_t col,const Element *v,Int_t n=-1);
	   //virtual void                     ExtractRow      (Int_t row,Int_t col,      Element *v,Int_t n=-1) const;

	   //virtual TMatrixTBase<Element>   &ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros=-1);
	   //  virtual TMatrixTBase<Element>   &ResizeTo(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros=-1);
	   //inline  TMatrixTBase<Element>   &ResizeTo(const TMatrixTSparse<Element> &m) {return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),
	   //                                                                                m.GetColUpb(),m.GetNoElements()); }

//virtual void Clear(Option_t * /*option*/ ="") { if (this->fIsOwner) {
//                                                   if (fElements) delete [] fElements; fElements = 0;
//                                                     if (fRowIndex) delete [] fRowIndex; fRowIndex = 0;
//                                                     if (fColIndex) delete [] fColIndex; fColIndex = 0;
//                                                   }
//                                                   this->fNelems    = 0;
//                                                   this->fNrowIndex = 0;
//                                                 }

//           TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
//                                           Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
//           TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
//                                           Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
//           TMatrixTSparse<Element> &Use   (TMatrixTSparse<Element> &a);
//   virtual TMatrixTBase<Element>   &GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
//                                            TMatrixTBase<Element> &target,Option_t *option="S") const;
//           TMatrixTSparse<Element>  GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
//   virtual TMatrixTBase<Element>   &SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source);

//   virtual Bool_t IsSymmetric() const { return (*this == TMatrixTSparse<Element>(kTransposed,*this)); }
	   TMatrixTSparse<Element> &Transpose (const TMatrixTSparse<Element> &source);
	   inline TMatrixTSparse<Element> &T () { return this->Transpose(*this); }

//   inline void Mult(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b) { AMultB(a,b,0); }

//   virtual TMatrixTBase<Element> &Zero       ();
//   virtual TMatrixTBase<Element> &UnitMatrix ();

//   virtual Element RowNorm () const;
//   virtual Element ColNorm () const;
      virtual Int_t   NonZeros() const { return this->fNelems; }

//   virtual TMatrixTBase<Element> &NormByDiag(const TVectorT<Element> &/*v*/,Option_t * /*option*/)
//                                              { MayNotUse("NormByDiag"); return *this; }

   // Either access a_ij as a(i,j)
//   inline       Element                          operator()(Int_t rown,Int_t coln) const;
//                Element                         &operator()(Int_t rown,Int_t coln);

   // or as a[i][j]
//   inline const TMatrixTSparseRow_const<Element> operator[](Int_t rown) const { return TMatrixTSparseRow_const<Element>(*this,rown); }
//   inline       TMatrixTSparseRow      <Element> operator[](Int_t rown)       { return TMatrixTSparseRow      <Element>(*this,rown); }

//   TMatrixTSparse<Element> &operator=(const TMatrixT<Element>       &source);
//   TMatrixTSparse<Element> &operator=(const TMatrixTSparse<Element> &source);

//   TMatrixTSparse<Element> &operator= (Element val);
//   TMatrixTSparse<Element> &operator-=(Element val);
//   TMatrixTSparse<Element> &operator+=(Element val);
//   TMatrixTSparse<Element> &operator*=(Element val);

//   TMatrixTSparse<Element> &operator+=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                if (this == &source) APlusB (tmp,tmp,1);
//                                                                                else                 APlusB (tmp,source,1); return *this; }
//   TMatrixTSparse<Element> &operator+=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                APlusB(tmp,source,1); return *this; }
//   TMatrixTSparse<Element> &operator-=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                if (this == &source) AMinusB (tmp,tmp,1);
//                                                                                else                 AMinusB(tmp,source,1); return *this; }
//   TMatrixTSparse<Element> &operator-=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                AMinusB(tmp,source,1); return *this; }
//   TMatrixTSparse<Element> &operator*=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                if (this == &source) AMultB (tmp,tmp,1);
//                                                                                else                 AMultB (tmp,source,1); return *this; }
//   TMatrixTSparse<Element> &operator*=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this);
//                                                                                AMultB(tmp,source,1); return *this; }

//   virtual TMatrixTBase  <Element> &Randomize  (Element alpha,Element beta,Double_t &seed);
//   virtual TMatrixTSparse<Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);

//   ClassDef(TMatrixTSparse,3) // Template of Sparse Matrix class
};

template <class Element> inline const Element *TMatrixTSparse<Element>::GetMatrixArray  () const { return fElements; }
template <class Element> inline       Element *TMatrixTSparse<Element>::GetMatrixArray  ()       { return fElements; }
template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetRowIndexArray() const { return fRowIndex; }
template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetRowIndexArray()       { return fRowIndex; }
template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetColIndexArray() const { return fColIndex; }
template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetColIndexArray()       { return fColIndex; }

//template <class Element>
//inline TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
//                                                                Int_t *pRowIndex,Int_t *pColIndex,Element *pData)
//                                                                  { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
//template <class Element>
//inline TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (TMatrixTSparse<Element> &a)
//                                                                  { R__ASSERT(a.IsValid());
//                                                                     return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
//                                                                                a.GetNoElements(),a.GetRowIndexArray(),
//                                                                                a.GetColIndexArray(),a.GetMatrixArray()); }
//template <class Element>
//inline TMatrixTSparse<Element>  TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
//                                                                Option_t *option) const
//                                                                  {
//                                                                    TMatrixTSparse<Element> tmp;
//                                                                    this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
//                                                                    return tmp;
//                                                                  }

//template <class Element>
//inline Element TMatrixTSparse<Element>::operator()(Int_t rown,Int_t coln) const
//{
//   R__ASSERT(this->IsValid());
//   if (this->fNrowIndex > 0 && this->fRowIndex[this->fNrowIndex-1] == 0) {
//      Error("operator=()(Int_t,Int_t) const","row/col indices are not set");
//      printf("fNrowIndex = %d fRowIndex[fNrowIndex-1] = %d\n",this->fNrowIndex,this->fRowIndex[this->fNrowIndex-1]);
//      return 0.0;
//   }
//   const Int_t arown = rown-this->fRowLwb;
//   const Int_t acoln = coln-this->fColLwb;
//   R__ASSERT(arown < this->fNrows && arown >= 0);
//   R__ASSERT(acoln < this->fNcols && acoln >= 0);
//   const Int_t sIndex = fRowIndex[arown];
//   const Int_t eIndex = fRowIndex[arown+1];
//   const Int_t index = TMath::BinarySearch(eIndex-sIndex,fColIndex+sIndex,acoln)+sIndex;
//   if (index < sIndex || fColIndex[index] != acoln) return 0.0;
//   else                                             return fElements[index];
//}

//template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
//template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
//template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
//template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source ,      Element                  val    );
//template <class Element> TMatrixTSparse<Element>  operator+ (      Element                  val    ,const TMatrixTSparse<Element> &source );
//template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
//template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
//template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
//template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source ,      Element                  val    );
//template <class Element> TMatrixTSparse<Element>  operator- (      Element                  val    ,const TMatrixTSparse<Element> &source );
//template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
//template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
//template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
///template <class Element> TMatrixTSparse<Element>  operator* (      Element                  val    ,const TMatrixTSparse<Element> &source );
//template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source,       Element                  val    );

//template <class Element> TMatrixTSparse<Element> &Add        (TMatrixTSparse<Element> &target,      Element                   scalar,
//                                                              const TMatrixTSparse<Element> &source);
//template <class Element> TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
//template <class Element> TMatrixTSparse<Element> &ElementDiv (TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);

//template <class Element> Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose=0);

//#endif

max_stars_repo_name: sharugupta/OpenUH

max_issues_repo_name: sharugupta/OpenUH

max_forks_repo_name: sharugupta/OpenUH


text:
/*  Standard C Libraries  */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "interpolation.h"
#include "float.h"

#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

typedef struct {
    int start_index;
    int end_index;
    int n;
    double* x;
    double* y;
} segment;

typedef struct {
    int n;
    segment * seg;
} segments;

void free_segments(segments segs)
{
    for(int i=0; i<segs.n; i++)
    {
        free(segs.seg[i].x);
        free(segs.seg[i].y);
    }
//    free(segs.seg);
}

/***
 * What is segment:
 * mask                    | 00000000000000000111111111111100000000000000000 |
 *                         |       signal    |  no signal  |   signal        |
 * signal data             | ~~|~~~~~~~|~~~~~|_____________|~~~~~|~~~~~~~|~~ |
 * params                  |   | edges | cut |             | cut | edges |   |
 * stored for intrpolation |   |-------|                         |-------|   |
 * will interpolate        |           |-------------------------|           |
 * start_index             |                 X                               |
 * end_index               |                               X                 |
 * start_data              |   X                                             |
 * start_data_to           |           X                                     |
 * end_data                |                                     X           |
 ***/

segments add_new_segment(segments segs, int start_index, int end_index, int edges, int cut, vector data)
{
    int start_data = start_index - edges - cut;
    if(start_data < 0) start_data = 0;

    int start_data_to = start_index - cut;
    if(start_data_to < 0) start_data_to = 0;

    int end_data = end_index + cut;
    if(end_data > data.x) end_data = data.x - 1;

    int new_size = segs.n + 1;
    segs.seg = realloc(segs.seg, sizeof(segment)*new_size);

    int data_size = edges * 2;

    segs.seg[segs.n].start_index = start_data_to;
    segs.seg[segs.n].end_index = end_data;
    segs.seg[segs.n].n = data_size;
    segs.seg[segs.n].x = calloc(data_size, sizeof(double));
    segs.seg[segs.n].y = calloc(data_size, sizeof(double));

    for(int i=0; i<edges; i++)
    {
        int x = i + start_data;
        segs.seg[segs.n].x[i] = x;
        segs.seg[segs.n].y[i] = getv(data, x);

        x = edges + i + end_data;
        segs.seg[segs.n].x[i+edges] = x;
        segs.seg[segs.n].y[i+edges] = getv(data, x);
    }

    segs.n = new_size;
    return segs;
}

segments vector_split_by_mask(vector data, vector mask, int edges, int cut)
{
    segments segs;
    segs.n = 0;
    segs.seg = calloc(1, sizeof(segment));

    int start_index = -1;
    int end_index = -1;

    for(int i=0; i<data.x; i++)
    {
        if(start_index == -1 && getv(mask, i) == 1.0)
        {
            start_index = i;
            continue;
        }
        if(start_index != -1 && getv(mask, i) != 1.0)
        {
            end_index = i;
            segs = add_new_segment(segs, start_index, end_index, edges, cut, data);

            start_index = -1;
            end_index = -1;
            continue;
        }
    }

//    if(start_index != -1)
//    {
//        end_index = data.x - 1;
//        segs = add_new_segment(segs, start_index, end_index, edges, 0, data);
//    }

    return segs;
}

vector vector_interpolate(vector data, vector mask, int cut, const gsl_interp_type * T)
{
    int edges = gsl_interp_type_min_size(T);
    segments segs = vector_split_by_mask(data, mask, edges, cut);

    vector result = copyv(data);

    for(int i=0; i<segs.n; i++)
    {
        gsl_interp_accel *acc = gsl_interp_accel_alloc ();
        gsl_spline *liner = gsl_spline_alloc(T, segs.seg[i].n);
        gsl_spline_init (liner, segs.seg[i].x, segs.seg[i].y, segs.seg[i].n);

        for(int j=segs.seg[i].start_index; j<segs.seg[i].end_index; j++)
        {
            setv(result, j, gsl_spline_eval (liner, j, acc));
        }

        gsl_spline_free (liner);
        gsl_interp_accel_free (acc);
    }

    free_segments(segs);

    return result;
}

vector vector_interpolate_by_mask(vector data, vector mask, int cut, int type)
{
    vector result;
    switch (type) {
    case 1: // Linear interpolation
        result = vector_interpolate(data, mask, cut, gsl_interp_linear);
        break;
    case 2: // Cubic spline (тatural boundary conditions)
        result = vector_interpolate(data, mask, cut, gsl_interp_cspline);
        break;
    case 3: // Cubic spline (periodic boundary conditions)
        result = vector_interpolate(data, mask, cut, gsl_interp_cspline_periodic);
        break;
    case 4: // Non-rounded Akima spline (natural boundary conditions)
        result = vector_interpolate(data, mask, cut, gsl_interp_akima);
        break;
    case 5: // Non-rounded Akima spline (periodic boundary conditions)
        result = vector_interpolate(data, mask, cut, gsl_interp_akima_periodic);
        break;
    case 6: // Steffen’s method
        result = vector_interpolate(data, mask, cut, gsl_interp_steffen);
        break;
    default:
        result = copyv(data);
        break;
    }
    return result;
}

void interpolate_part(vector * data, int start, int end, const gsl_interp_type * T)
{
    int edges = gsl_interp_type_min_size(T);

    int start_index = start - edges;
    if (start_index < 0) start_index = 0;
    int end_index = end + edges;
    if (end_index >= data->x) end_index = data->x;

    int n = (start - start_index) + (end_index - end) + 1;
    int index = 0;
    vector x = zerov(n);
    vector y = zerov(n);

    for(int i=start_index; i<=start; i++)
    {
//        printf("i %i\n", i);
//        fprintf(stdout, "val %f\n", data->v[i]);
        setv(y, index, data->v[i]);
        setv(x, index, i);
        index++;
    }

    for(int i=end; i<end_index; i++)
    {
//        printf("i %i\n", i);
//        fprintf(stdout, "val %f\n", data->v[i]);
        setv(y, index, getv(*data, i));
        setv(x, index, i);
        index++;
    }

//    for (int i=start_index; i<end_index; i++)
//    {
//        if((i<=start_index || i>=end_index)&&(index<n)&&(data->v[i]!=0||i==start_index||i==end_index))
//        {
//            setv(y, index, getv(data, i));
//            setv(x, index, i);
//            index++;
//        }
//    }

    n = index;

    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    gsl_spline *liner = gsl_spline_alloc(T, n);
    gsl_spline_init (liner, x.v, y.v, n);

    for(int j=(start); j<end; j++)
    {
        data->v[j] = gsl_spline_eval (liner, j, acc);
    }

    gsl_spline_free (liner);
    gsl_interp_accel_free (acc);

    freev(x);
    freev(y);
}

void vector_interpolate_part(vector * data, int start, int end, int type)
{
    if(end - start > 1)
        switch (type) {
        case 1: // Linear interpolation
            interpolate_part(data, start, end, gsl_interp_linear);
            break;
        case 2: // Cubic spline (тatural boundary conditions)
            interpolate_part(data, start, end, gsl_interp_cspline);
            break;
        case 3: // Cubic spline (periodic boundary conditions)
            interpolate_part(data, start, end, gsl_interp_cspline_periodic);
            break;
        case 4: // Non-rounded Akima spline (natural boundary conditions)
            interpolate_part(data, start, end, gsl_interp_akima);
            break;
        case 5: // Non-rounded Akima spline (periodic boundary conditions)
            interpolate_part(data, start, end, gsl_interp_akima_periodic);
            break;
        case 6: // Steffen’s method
             interpolate_part(data, start, end, gsl_interp_steffen);
            break;
        }
}

ext: c

max_stars_repo_name: zhitko/sptk-analyzer

max_issues_repo_name: zhitko/sptk-analyzer

max_forks_repo_name: zhitko/inton-trainer


text:
#include <cblas.h>

#include "tasks.h"
#include "utils.h"

void trsm_task_par_reconfigure(int nth)
{
    // empty
}

void trsm_task_par_finalize(void)
{
    // empty
}

void trsm_task_par(void *ptr, int nth, int me)
{
    struct trsm_task_arg *arg = (struct trsm_task_arg*) ptr;

    int n       = arg->n;
    int m       = arg->m;
    double *A21 = arg->A21;
    double *A11 = arg->A11;
    int ldA     = arg->ldA;

    // Compute nominal block size.
    int blksz = iceil(m, nth);

    // Determine my share of the rows of A21.
    int my_first_row = blksz * me;
    int my_num_rows  = min(blksz, m - my_first_row);

    // Compute A21 := A21 * inv(A11'), using my block of A21.
    if (my_num_rows > 0) {
        cblas_dtrsm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit,
                    my_num_rows, n,
                    1.0, A11,                ldA,
                         A21 + my_first_row, ldA);
    }

}

max_stars_repo_name: NLAFET/pcp-runtime

max_issues_repo_name: NLAFET/pcp-runtime

max_forks_repo_name: NLAFET/pcp-runtime


text:
#ifndef problems_h
#define problems_h

#include <ceed.h>
#include <petsc.h>
#include "../problems/cl-problems.h"
#include "../problems/neo-hookean.h"
#include "../problems/mooney-rivlin.h"

// Physics options
#define SOLIDS_PROBLEM_REGISTER(list, name, fname, physics) \
  ierr = PetscFunctionListAdd(&list->setupPhysics, name,                          \
                              PhysicsContext_ ## physics); CHKERRQ(ierr);         \
  ierr = PetscFunctionListAdd(&list->setupSmootherPhysics, name,                  \
                              PhysicsSmootherContext_ ## physics); CHKERRQ(ierr); \
  ierr = PetscFunctionListAdd(&list->setupLibceedFineLevel, name,                 \
                              SetupLibceedFineLevel_ ## fname); CHKERRQ(ierr);    \
  ierr = PetscFunctionListAdd(&list->setupLibceedLevel, name,                     \
                              SetupLibceedLevel_ ## fname); CHKERRQ(ierr);        \

typedef struct ProblemFunctions_ *ProblemFunctions;
struct ProblemFunctions_ {
  PetscFunctionList setupPhysics, setupSmootherPhysics, setupLibceedFineLevel,
                    setupLibceedLevel;
};

PetscErrorCode RegisterProblems(ProblemFunctions problem_functions);

#define SOLIDS_PROBLEM(name) \
  PetscErrorCode SetupLibceedFineLevel_ ## name (DM dm, DM dm_energy,           \
    DM dm_diagnostic, Ceed ceed, AppCtx app_ctx, CeedQFunctionContext phys_ctx, \
    PetscInt fine_level, PetscInt num_comp_u, PetscInt U_g_size,                \
    PetscInt U_loc_size, CeedVector force_ceed, CeedVector neumann_ceed,        \
    CeedData *data);                                                            \
  PetscErrorCode SetupLibceedLevel_ ## name (DM dm, Ceed ceed,                  \
    AppCtx app_ctx, PetscInt level, PetscInt num_comp_u, PetscInt U_g_size,     \
    PetscInt u_loc_size, CeedVector fine_mult, CeedData *data);                 \

SOLIDS_PROBLEM(ElasLinear);
SOLIDS_PROBLEM(ElasSSNH);
SOLIDS_PROBLEM(ElasFSCurrentNH1);
SOLIDS_PROBLEM(ElasFSCurrentNH2);
SOLIDS_PROBLEM(ElasFSInitialNH1);
SOLIDS_PROBLEM(ElasFSInitialNH2);
SOLIDS_PROBLEM(ElasFSInitialMR1);

#endif //problems_h

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


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

void
cblas_dsyr2k (const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
	      const enum CBLAS_TRANSPOSE Trans, const int N, const int K,
	      const double alpha, const double *A, const int lda,
	      const double *B, const int ldb, const double beta, double *C,
	      const int ldc)
{
#define BASE double
#include "source_syr2k_r.h"
#undef BASE
}


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_cblas.h>
#include "cblas.h"
#include "error_cblas_l2.h"

void
cblas_sspmv (const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
             const int N, const float alpha, const float *Ap, const float *X,
             const int incX, const float beta, float *Y, const int incY)
{
#define BASE float
#include "source_spmv.h"
#undef BASE
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/**
 *
 * @file zcposv.c
 *
 *  PLASMA computational routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Emmanuel Agullo
 * @date 2010-11-15
 * @precisions mixed zc -> ds
 *
 **/
#include <stdlib.h>
#include <math.h>
#include <lapacke.h>
#include "common.h"

#define PLASMA_zlag2c(_descA, _descSB)                \
  plasma_parallel_call_4(plasma_pzlag2c,              \
                         PLASMA_desc,      (_descA),  \
                         PLASMA_desc,      (_descSB), \
                         PLASMA_sequence*, sequence,  \
                         PLASMA_request*,  request)

#define PLASMA_clag2z(_descSA, _descB)                \
  plasma_parallel_call_4(plasma_pclag2z,              \
                         PLASMA_desc,      (_descSA), \
                         PLASMA_desc,      (_descB),  \
                         PLASMA_sequence*, sequence,  \
                         PLASMA_request*,  request)

#define PLASMA_zlange(_norm, _descA, _result, _work)   \
  _result = 0;                                         \
  plasma_parallel_call_6(plasma_pzlange,               \
                         PLASMA_enum,      (_norm),    \
                         PLASMA_desc,      (_descA),   \
                         double*,          (_work),    \
                         double*,          &(_result), \
                         PLASMA_sequence*, sequence,   \
                         PLASMA_request*,  request);

#define PLASMA_zlanhe(_norm, _uplo, _descA, _result, _work) \
  _result = 0;                                              \
  plasma_parallel_call_7(plasma_pzlanhe,                    \
                         PLASMA_enum,      (_norm),         \
                         PLASMA_enum,      (_uplo),         \
                         PLASMA_desc,      (_descA),        \
                         double*,          (_work),         \
                         double*,          &(_result),      \
                         PLASMA_sequence*, sequence,        \
                         PLASMA_request*,  request);

#define PLASMA_zlacpy(_descA, _descB)                        \
  plasma_parallel_call_5(plasma_pzlacpy,                     \
                         PLASMA_enum,      PlasmaUpperLower, \
                         PLASMA_desc,      (_descA),         \
                         PLASMA_desc,      (_descB),         \
                         PLASMA_sequence*, sequence,         \
                         PLASMA_request*,  request)

#define PLASMA_zgeadd(_alpha, _descA, _descB)           \
  plasma_parallel_call_5(plasma_pzgeadd,                \
                         PLASMA_Complex64_t, (_alpha),  \
                         PLASMA_desc,        (_descA),  \
                         PLASMA_desc,        (_descB),  \
                         PLASMA_sequence*,   sequence,  \
                         PLASMA_request*,    request)

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t
 *
 *  PLASMA_zcposv - Computes the solution to a system of linear equations A * X = B,
 *  where A is an N-by-N symmetric positive definite (or Hermitian positive definite
 *  in the complex case) matrix and X and B are N-by-NRHS matrices.
 *  The Cholesky decomposition is used to factor A as
 *
 *    A = U**H * U, if uplo = PlasmaUpper, or
 *    A = L * L**H, if uplo =  PlasmaLower,
 *
 *  where U is an upper triangular matrix and  L is a lower triangular matrix.
 *  The factored form of A is then used to solve the system of equations A * X = B.
 *
 *  PLASMA_zcposv first attempts to factorize the matrix in COMPLEX and use this
 *  factorization within an iterative refinement procedure to produce a
 *  solution with COMPLEX*16 normwise backward error quality (see below).
 *  If the approach fails the method switches to a COMPLEX*16
 *  factorization and solve.
 *
 *  The iterative refinement is not going to be a winning strategy if
 *  the ratio COMPLEX performance over COMPLEX*16 performance is too
 *  small. A reasonable strategy should take the number of right-hand
 *  sides and the size of the matrix into account. This might be done
 *  with a call to ILAENV in the future. Up to now, we always try
 *  iterative refinement.
 *
 *  The iterative refinement process is stopped if ITER > ITERMAX or
 *  for all the RHS we have: RNRM < N*XNRM*ANRM*EPS*BWDMAX
 *  where:
 *
 *  - ITER is the number of the current iteration in the iterative refinement process
 *  - RNRM is the infinity-norm of the residual
 *  - XNRM is the infinity-norm of the solution
 *  - ANRM is the infinity-operator-norm of the matrix A
 *  - EPS is the machine epsilon returned by DLAMCH('Epsilon').
 *
 *  Actually, in its current state (PLASMA 2.1.0), the test is slightly relaxed.
 *
 *  The values ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in] N
 *          The number of linear equations, i.e., the order of the matrix A. N >= 0.
 *
 * @param[in] NRHS
 *          The number of right hand sides, i.e., the number of columns of the matrix B.
 *          NRHS >= 0.
 *
 * @param[in] A
 *          The N-by-N symmetric positive definite (or Hermitian) coefficient matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular part of A
 *          contains the upper triangular part of the matrix A, and the strictly lower triangular
 *          part of A is not referenced.
 *          If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
 *          triangular part of the matrix A, and the strictly upper triangular part of A is not
 *          referenced.
 *          This matrix is not modified.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,N).
 *
 * @param[in] B
 *          The N-by-NRHS matrix of right hand side matrix B.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDB >= max(1,N).
 *
 * @param[out] X
 *          If return value = 0, the N-by-NRHS solution matrix X.
 *
 * @param[in] LDX
 *          The leading dimension of the array B. LDX >= max(1,N).
 *
 * @param[out] ITER
 *          The number of the current iteration in the iterative refinement process
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *          \retval >0 if i, the leading minor of order i of A is not positive definite, so the
 *               factorization could not be completed, and the solution has not been computed.
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcposv_Tile
 * @sa PLASMA_zcposv_Tile_Async
 * @sa PLASMA_dsposv
 * @sa PLASMA_zposv
 *
 ******************************************************************************/
int PLASMA_zcposv(PLASMA_enum uplo, int N, int NRHS,
                  PLASMA_Complex64_t *A, int LDA,
                  PLASMA_Complex64_t *B, int LDB,
                  PLASMA_Complex64_t *X, int LDX, int *ITER)
{
    int NB;
    int status;
    PLASMA_desc  descA;
    PLASMA_desc  descB;
    PLASMA_desc  descX;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zcposv", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }

    *ITER = 0;

    /* Check input arguments */
    if (uplo != PlasmaUpper && uplo != PlasmaLower) {
        plasma_error("PLASMA_zcposv", "illegal value of uplo");
        return -1;
    }
    if (N < 0) {
        plasma_error("PLASMA_zcposv", "illegal value of N");
        return -2;
    }
    if (NRHS < 0) {
        plasma_error("PLASMA_zcposv", "illegal value of NRHS");
        return -3;
    }
    if (LDA < max(1, N)) {
        plasma_error("PLASMA_zcposv", "illegal value of LDA");
        return -5;
    }
    if (LDB < max(1, N)) {
        plasma_error("PLASMA_zcposv", "illegal value of LDB");
        return -7;
    }
    if (LDX < max(1, N)) {
        plasma_error("PLASMA_zcposv", "illegal value of LDX");
        return -9;
    }
    /* Quick return - currently NOT equivalent to LAPACK's
     * LAPACK does not have such check for ZCPOSV */
    if (min(N, NRHS) == 0)
        return PLASMA_SUCCESS;

    /* Tune NB depending on M, N & NRHS; Set NBNBSIZE */
    status = plasma_tune(PLASMA_FUNC_ZCPOSV, N, N, NRHS);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcposv", "plasma_tune() failed");
        return status;
    }

    NB = PLASMA_NB;

    plasma_sequence_create(plasma, &sequence);

    /* DOUBLE PRECISION INITIALIZATION */
    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_zooplap2tile( descA, A, NB, NB, LDA, N,    0, 0, N, N,    sequence, &request,
                             plasma_desc_mat_free(&(descA)) );
        plasma_zooplap2tile( descB, B, NB, NB, LDB, NRHS, 0, 0, N, NRHS, sequence, &request,
                             plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descB)) );
        plasma_zdesc_alloc(  descX, NB, NB, N, NRHS, 0, 0, N, NRHS, plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descB)); plasma_desc_mat_free(&(descX)) );
    } else {
        plasma_ziplap2tile( descA, A, NB, NB, LDA, N,    0, 0, N, N,
                            sequence, &request);
        plasma_ziplap2tile( descB, B, NB, NB, LDB, NRHS, 0, 0, N, NRHS,
                            sequence, &request);

        descX = plasma_desc_init(
            PlasmaComplexDouble, NB, NB, (NB*NB),
            LDX, NRHS, 0, 0, N, NRHS);
        descX.mat = X;
    }

    /* Call the native interface */
    status = PLASMA_zcposv_Tile_Async(uplo, &descA, &descB, &descX, ITER, sequence, &request);

    if (status == PLASMA_SUCCESS) {
        if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
            plasma_zooptile2lap( descX, X, NB, NB, LDX, NRHS,  sequence, &request);
            plasma_dynamic_sync();
            plasma_desc_mat_free(&descA);
            plasma_desc_mat_free(&descB);
            plasma_desc_mat_free(&descX);
        } else {
            plasma_ziptile2lap( descA, A, NB, NB, LDA, N,     sequence, &request);
            plasma_ziptile2lap( descB, B, NB, NB, LDB, NRHS,  sequence, &request);
            plasma_ziptile2lap( descX, X, NB, NB, LDX, NRHS,  sequence, &request);
            plasma_dynamic_sync();
        }
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile
 *
 *  PLASMA_zcposv_Tile - Solves a symmetric positive definite or Hermitian positive definite
 *  system of linear equations using the Cholesky factorization and mixed-precision iterative refinement.
 *  Tile equivalent of PLASMA_zcposv().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] uplo
 *          Specifies whether the matrix A is upper triangular or lower triangular:
 *          = PlasmaUpper: Upper triangle of A is stored;
 *          = PlasmaLower: Lower triangle of A is stored.
 *
 * @param[in,out] A
 *          On entry, the N-by-N symmetric positive definite (or Hermitian) coefficient matrix A.
 *          If uplo = PlasmaUpper, the leading N-by-N upper triangular part of A
 *          contains the upper triangular part of the matrix A, and the strictly lower triangular
 *          part of A is not referenced.
 *          If UPLO = 'L', the leading N-by-N lower triangular part of A contains the lower
 *          triangular part of the matrix A, and the strictly upper triangular part of A is not
 *          referenced.
 *          - If the iterative refinement converged, A is not modified;
 *          - otherwise, it falled backed to double precision solution,
 *
 * @param[in] B
 *          On entry, the N-by-NRHS matrix of right hand side matrix B.
 *
 * @param[out] X
 *          On exit, if return value = 0, the N-by-NRHS solution matrix X.
 *
 * @param[out] ITER
 *          The number of the current iteration in the iterative refinement process
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval >0 if i, the leading minor of order i of A is not positive definite, so the
 *               factorization could not be completed, and the solution has not been computed.
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcposv
 * @sa PLASMA_zcposv_Tile_Async
 * @sa PLASMA_dsposv_Tile
 * @sa PLASMA_zposv_Tile
 *
 ******************************************************************************/
int PLASMA_zcposv_Tile(PLASMA_enum uplo, PLASMA_desc *A, PLASMA_desc *B,
                       PLASMA_desc *X, int *ITER)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zcposv_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    status = PLASMA_zcposv_Tile_Async(uplo, A, B, X, ITER, sequence, &request);
    if (status != PLASMA_SUCCESS)
        return status;
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile_Async
 *
 *  PLASMA_zcposv_Tile_Async - Solves a symmetric positive definite or Hermitian
 *  positive definite system of linear equations using the Cholesky factorization
 *  and mixed-precision iterative refinement.
 *  Non-blocking equivalent of PLASMA_zcposv_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcposv
 * @sa PLASMA_zcposv_Tile
 * @sa PLASMA_dsposv_Tile_Async
 * @sa PLASMA_zposv_Tile_Async
 *
 ******************************************************************************/
int PLASMA_zcposv_Tile_Async(PLASMA_enum uplo, PLASMA_desc *A, PLASMA_desc *B,
                             PLASMA_desc *X, int *ITER,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    int N, NB;
    PLASMA_desc descA;
    PLASMA_desc descB;
    PLASMA_desc descX;
    plasma_context_t *plasma;
    double *work;
    PLASMA_desc descR, descSA, descSX;

    const int itermax = 30;
    const double bwdmax = 1.0;
    const PLASMA_Complex64_t negone = -1.0;
    const PLASMA_Complex64_t one = 1.0;
    int iiter;
    double Anorm, cte, eps, Rnorm, Xnorm;
    *ITER=0;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zcposv_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_zcposv_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_zcposv_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcposv_Tile_Async", "invalid first descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descA = *A;
    }
    if (plasma_desc_check(B) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcposv_Tile_Async", "invalid second descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descB = *B;
    }
    if (plasma_desc_check(X) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcposv_Tile_Async", "invalid third descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descX = *X;
    }
    /* Check input arguments */
    if (descA.nb != descA.mb || descB.nb != descB.mb || descX.nb != descX.mb) {
        plasma_error("PLASMA_zcposv_Tile_Async", "only square tiles supported");
        return PLASMA_ERR_ILLEGAL_VALUE;
    }
    if (uplo != PlasmaUpper && uplo != PlasmaLower) {
        plasma_error("PLASMA_zcposv_Tile_Async", "illegal value of uplo");
        return -1;
    }
    /* Quick return - currently NOT equivalent to LAPACK's
     * LAPACK does not have such check for DPOSV */

/*
    if (min(N, NRHS) == 0)
        return PLASMA_SUCCESS;
*/

    /* Set N, NRHS */
    N  = descA.m;
    NB = descA.nb;

    work = (double *)plasma_shared_alloc(plasma, PLASMA_SIZE, PlasmaRealDouble);
    if (work == NULL) {
        plasma_error("PLASMA_zcposv_Tile_Async", "plasma_shared_alloc() failed");
        plasma_shared_free(plasma, work);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }

    plasma_zdesc_alloc( descR,  NB, NB, descB.m, descB.n, 0, 0, descB.m, descB.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR) );
    plasma_cdesc_alloc( descSA, NB, NB, descA.m, descA.n, 0, 0, descA.m, descA.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR); plasma_desc_mat_free(&descSA) );
    plasma_cdesc_alloc( descSX, NB, NB, descX.m, descX.n, 0, 0, descX.m, descX.n, plasma_shared_free( plasma, work ); plasma_desc_mat_free(&descR); plasma_desc_mat_free(&descSA); plasma_desc_mat_free(&descSX) );

    /* Compute some constants */
    PLASMA_zlanhe(PlasmaInfNorm, uplo, descA, Anorm, work);
    eps = LAPACKE_dlamch_work('e');

    /* Convert B from double precision to single precision and store
       the result in SX. */
    PLASMA_zlag2c(descB, descSX);
    if (sequence->status != PLASMA_SUCCESS)
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Convert A from double precision to single precision and store
       the result in SA. */
    PLASMA_zlag2c(descA, descSA);
    if (sequence->status != PLASMA_SUCCESS)
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Compute the Cholesky factorization of SA */
    plasma_parallel_call_4(plasma_pcpotrf,
        PLASMA_enum, uplo,
        PLASMA_desc, descSA,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Solve the system SA*SX = SB */
    /* Forward substitution */
    plasma_parallel_call_9(plasma_pctrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, uplo,
        PLASMA_enum, uplo == PlasmaUpper ? PlasmaConjTrans : PlasmaNoTrans,
        PLASMA_enum, PlasmaNonUnit,
        PLASMA_Complex32_t, 1.0,
        PLASMA_desc, descSA,
        PLASMA_desc, descSX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Backward substitution */
    plasma_parallel_call_9(plasma_pctrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, uplo,
        PLASMA_enum, uplo == PlasmaUpper ? PlasmaNoTrans : PlasmaConjTrans,
        PLASMA_enum, PlasmaNonUnit,
        PLASMA_Complex32_t, 1.0,
        PLASMA_desc, descSA,
        PLASMA_desc, descSX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Convert SX back to double precision */
    PLASMA_clag2z(descSX, descX);

    /* Compute R = B - AX. */
    PLASMA_zlacpy(descB,descR);
    plasma_parallel_call_9(plasma_pzhemm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, uplo,
        PLASMA_Complex64_t, negone,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_Complex64_t, one,
        PLASMA_desc, descR,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Check whether the NRHS normwise backward error satisfies the
       stopping criterion. If yes return. Note that ITER=0 (already set). */
    PLASMA_zlange(PlasmaInfNorm, descX, Xnorm, work);
    PLASMA_zlange(PlasmaInfNorm, descR, Rnorm, work);

    /* Wait for the end of Anorm, Xnorm and Bnorm computations */
    plasma_dynamic_sync();

    cte = Anorm*eps*((double) N)*bwdmax;
    if (Rnorm < Xnorm * cte){
        /* The NRHS normwise backward errors satisfy the
           stopping criterion. We are good to exit. */
        plasma_desc_mat_free(&descSA);
        plasma_desc_mat_free(&descSX);
        plasma_desc_mat_free(&descR);
        plasma_shared_free(plasma, work);
        return PLASMA_SUCCESS;
    }

    /* Iterative refinement */
    for (iiter = 0; iiter < itermax; iiter++){

        /* Convert R from double precision to single precision
           and store the result in SX. */
        PLASMA_zlag2c(descR, descSX);

        /* Solve the system SA*SX = SR */
        /* Forward substitution */
        plasma_parallel_call_9(plasma_pctrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, uplo,
            PLASMA_enum, uplo == PlasmaUpper ? PlasmaConjTrans : PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex32_t, 1.0,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        /* Backward substitution */
        plasma_parallel_call_9(plasma_pctrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, uplo,
            PLASMA_enum, uplo == PlasmaUpper ? PlasmaNoTrans : PlasmaConjTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex32_t, 1.0,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        /* Convert SX back to double precision and update the current
           iterate. */
        PLASMA_clag2z(descSX, descR);
        PLASMA_zgeadd(one, descR, descX);

        /* Compute R = B - AX. */
        PLASMA_zlacpy(descB,descR);
        plasma_parallel_call_9(plasma_pzhemm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, uplo,
            PLASMA_Complex64_t, negone,
            PLASMA_desc, descA,
            PLASMA_desc, descX,
            PLASMA_Complex64_t, one,
            PLASMA_desc, descR,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        /* Check whether the NRHS normwise backward errors satisfy the
           stopping criterion. If yes, set ITER=IITER>0 and return. */
        PLASMA_zlange(PlasmaInfNorm, descX, Xnorm, work);
        PLASMA_zlange(PlasmaInfNorm, descR, Rnorm, work);

        /* Wait for the end of Xnorm and Bnorm computations */
        plasma_dynamic_sync();

        if (Rnorm < Xnorm * cte){
            /* The NRHS normwise backward errors satisfy the
               stopping criterion. We are good to exit. */
            *ITER = iiter;

            plasma_desc_mat_free(&descSA);
            plasma_desc_mat_free(&descSX);
            plasma_desc_mat_free(&descR);
            plasma_shared_free(plasma, work);
            return PLASMA_SUCCESS;
        }
    }

    /* We have performed ITER=itermax iterations and never satisified
       the stopping criterion, set up the ITER flag accordingly and
       follow up on double precision routine. */
    *ITER = -itermax - 1;

    plasma_desc_mat_free(&descSA);
    plasma_desc_mat_free(&descSX);
    plasma_desc_mat_free(&descR);
    plasma_shared_free(plasma, work);

    /* Single-precision iterative refinement failed to converge to a
       satisfactory solution, so we resort to double precision. */

    plasma_parallel_call_4(plasma_pzpotrf,
        PLASMA_enum, uplo,
        PLASMA_desc, descA,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    PLASMA_zlacpy(descB,descX);

    plasma_parallel_call_9(plasma_pztrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, uplo,
        PLASMA_enum, uplo == PlasmaUpper ? PlasmaConjTrans : PlasmaNoTrans,
        PLASMA_enum, PlasmaNonUnit,
        PLASMA_Complex64_t, 1.0,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_parallel_call_9(plasma_pztrsm,
        PLASMA_enum, PlasmaLeft,
        PLASMA_enum, uplo,
        PLASMA_enum, uplo == PlasmaUpper ? PlasmaNoTrans : PlasmaConjTrans,
        PLASMA_enum, PlasmaNonUnit,
        PLASMA_Complex64_t, 1.0,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    return PLASMA_SUCCESS;
}


text:
// stdafx.h : include file for standard system include files,
// or project specific include files that are used frequently, but
// are changed infrequently
//

#pragma once
#pragma warning (disable : 4514 4710 4711) // inlining
#pragma warning (disable: 5045) // TODO Spectre

// #define WIN32_LEAN_AND_MEAN             // Exclude rarely-used stuff from Windows headers
// Windows Header Files:
#pragma warning (push, 1)
#pragma warning (disable : 5039 4668 4996 26814 4548 4355 4917 4702 26400 4987 4820 4365 4623 4625 4626 5026 5027 4571 4774 26412 26461 26426 26432 26447 26472 26446 26473 26440 26429 26496 26472 26482 26486 26487 26434)
#include <gsl.h>
#pragma warning (pop) // unbalanced push in span.h
#include <windows.h>
#include <Shlobj.h>
#include <Strsafe.h>

#include <algorithm>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <memory>
#include <string>
#include <vector>
#include <array>
#include <ranges>
#include <span>

#include "Resource.h"
#pragma warning (pop)

max_stars_repo_name: SteveGilham/drophash

max_issues_repo_name: SteveGilham/drophash

max_forks_repo_name: SteveGilham/drophash


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

/*
 * This module contains routines for tracking minimum/maximum values of a
 * moving fixed-sized window. It is based on the algorithm of:
 *
 * [1] Daniel Lemire, Streaming Maximum-Minimum Filter Using No More than Three Comparisons per Element,
 *     Nordic Journal of Computing, Volume 13, Number 4, pages 328-339, 2006
 *
 * Also available as a preprint here: https://arxiv.org/abs/cs/0610046
 */

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

typedef double mmacc_type_t;
typedef mmacc_type_t ringbuf_type_t;

#include "deque.c"
#include "ringbuf.c"

typedef struct
{
  size_t n;             /* window size */
  size_t k;             /* number of samples in current window */
  mmacc_type_t xprev;   /* previous sample added to window */
  ringbuf *rbuf;        /* ring buffer storing current window, size n */
  deque *minque;        /* double-ended queue of min values (L) */
  deque *maxque;        /* double-ended queue of max values (U) */
} mmacc_state_t;

static size_t mmacc_size(const size_t n);
static int mmacc_init(const size_t n, void * vstate);
static int mmacc_insert(const mmacc_type_t x, void * vstate);
static int mmacc_delete(void * vstate);
static int mmacc_min(void * params, mmacc_type_t * result, const void * vstate);
static int mmacc_max(void * params, mmacc_type_t * result, const void * vstate);
static int mmacc_minmax(void * params, mmacc_type_t * result, const void * vstate);

static size_t
mmacc_size(const size_t n)
{
  size_t size = 0;

  size += sizeof(mmacc_state_t);
  size += ringbuf_size(n);          /* rbuf */
  size += 2 * deque_size(n + 1);    /* minque/maxque */

  return size;
}

static int
mmacc_init(const size_t n, void * vstate)
{
  mmacc_state_t * state = (mmacc_state_t *) vstate;

  state->n = n;
  state->k = 0;
  state->xprev = 0.0;

  state->rbuf = (ringbuf *) ((unsigned char *) vstate + sizeof(mmacc_state_t));
  state->minque = (deque *) ((unsigned char *) state->rbuf + ringbuf_size(n));
  state->maxque = (deque *) ((unsigned char *) state->minque + deque_size(n + 1));

  ringbuf_init(n, state->rbuf);
  deque_init(n + 1, state->minque);
  deque_init(n + 1, state->maxque);

  return GSL_SUCCESS;
}

static int
mmacc_insert(const mmacc_type_t x, void * vstate)
{
  mmacc_state_t * state = (mmacc_state_t *) vstate;
  int head, tail;

  if (state->k == 0)
    {
      /* first sample */
      ringbuf_insert(x, state->rbuf);
      head = state->rbuf->head;
      deque_push_back(head, state->maxque);
      deque_push_back(head, state->minque);
    }
  else
    {
      if (x > state->xprev)
        {
          deque_pop_back(state->maxque);

          while (!deque_is_empty(state->maxque))
            {
              if (x <= state->rbuf->array[deque_peek_back(state->maxque)])
                break;

              deque_pop_back(state->maxque);
            }
        }
      else
        {
          deque_pop_back(state->minque);

          while (!deque_is_empty(state->minque))
            {
              if (x >= state->rbuf->array[deque_peek_back(state->minque)])
                break;

              deque_pop_back(state->minque);
            }
        }

      /* store new sample into ring buffer */
      tail = state->rbuf->tail;
      ringbuf_insert(x, state->rbuf);
      head = state->rbuf->head;

      deque_push_back(head, state->maxque);
      deque_push_back(head, state->minque);

      if (state->k == state->n)
        {
          /*
           * window is full - check if oldest window element is a global minimum/maximum
           * of current window - if so pop it from U/L queues;
           * the check head != tail ensures there is more than 1 element in the
           * queue, do not pop if queue has only 1 element, since this element would
           * be the newest sample
           */
          if (state->maxque->head != state->maxque->tail && tail == deque_peek_front(state->maxque))
            deque_pop_front(state->maxque);
          else if (state->minque->head != state->minque->tail && tail == deque_peek_front(state->minque))
            deque_pop_front(state->minque);
        }
    }

  if (state->k < state->n)
    ++(state->k);

  state->xprev = x;

  return GSL_SUCCESS;
}

static int
mmacc_delete(void * vstate)
{
  mmacc_state_t * state = (mmacc_state_t *) vstate;

  if (state->k > 0)
    {
      /*
       * check if oldest window element is a global minimum/maximum; if so
       * pop it from U/L queues
       */
      if (state->rbuf->tail == deque_peek_front(state->maxque))
        deque_pop_front(state->maxque);
      else if (state->rbuf->tail == deque_peek_front(state->minque))
        deque_pop_front(state->minque);

      /* remove oldest element from ring buffer */
      ringbuf_pop_back(state->rbuf);

      --(state->k);
    }

  return GSL_SUCCESS;
}

static int
mmacc_min(void * params, mmacc_type_t * result, const void * vstate)
{
  const mmacc_state_t * state = (const mmacc_state_t *) vstate;

  (void) params;

  if (state->k == 0)
    {
      GSL_ERROR ("no samples yet added to workspace", GSL_EINVAL);
    }
  else
    {
      *result = state->rbuf->array[deque_peek_front(state->minque)];
      return GSL_SUCCESS;
    }
}

static int
mmacc_max(void * params, mmacc_type_t * result, const void * vstate)
{
  const mmacc_state_t * state = (const mmacc_state_t *) vstate;

  (void) params;

  if (state->k == 0)
    {
      GSL_ERROR ("no samples yet added to workspace", GSL_EINVAL);
    }
  else
    {
      *result = state->rbuf->array[deque_peek_front(state->maxque)];
      return GSL_SUCCESS;
    }
}

static int
mmacc_minmax(void * params, mmacc_type_t * result, const void * vstate)
{
  const mmacc_state_t * state = (const mmacc_state_t *) vstate;

  (void) params;

  if (state->k == 0)
    {
      GSL_ERROR ("no samples yet added to workspace", GSL_EINVAL);
    }
  else
    {
      result[0] = state->rbuf->array[deque_peek_front(state->minque)];
      result[1] = state->rbuf->array[deque_peek_front(state->maxque)];
      return GSL_SUCCESS;
    }
}

static const gsl_movstat_accum min_accum_type =
{
  mmacc_size,
  mmacc_init,
  mmacc_insert,
  mmacc_delete,
  mmacc_min
};

const gsl_movstat_accum *gsl_movstat_accum_min = &min_accum_type;

static const gsl_movstat_accum max_accum_type =
{
  mmacc_size,
  mmacc_init,
  mmacc_insert,
  mmacc_delete,
  mmacc_max
};

const gsl_movstat_accum *gsl_movstat_accum_max = &max_accum_type;

static const gsl_movstat_accum minmax_accum_type =
{
  mmacc_size,
  mmacc_init,
  mmacc_insert,
  mmacc_delete,
  mmacc_minmax
};

const gsl_movstat_accum *gsl_movstat_accum_minmax = &minmax_accum_type;

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/*
 * bspline.c: Routines for calculating values of B-splines and their
 *  derivatives, as well as efficient computation of the values of
 *  N-dimensional B-spline surfaces.
 */

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

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

#include "photospline/bspline.h"

/*
 * Compute the value of the ith nth-order basis spline of a set
 * defined by knots at the point x.
 *
 * This is implemented using the De Boor algorithm, as outlined on
 * Wikipedia.
 */

double
bspline(const double *knots, double x, int i, int n)
{
	double result;

	if (n == 0) {
		/*
		 * Special case the 0th order case, where B-Splines
		 * are constant functions from one knot to the next.
		 */

		if (x >= knots[i] && x < knots[i+1])
			return 1.0;
		else
			return 0.0;
	}

	result = (x - knots[i])*bspline(knots, x, i, n-1) /
	    (knots[i+n] - knots[i]);
	result += (knots[i+n+1] - x)*bspline(knots, x, i+1, n-1) /
	    (knots[i+n+1] - knots[i+1]);

	return result;
}

/*
 * A brain-dead reimplementation of de Boor's BSPLVB, which generates
 * the values of the non-zero B-splines at x from the bottom up without
 * unnecessarily recalculating terms.
 *
 * NB: for bsplvb_simple(), bspline_nonzero(), and bspline_deriv_nonzero(),
 * `left' must be the index of the nearest fully-supported knot
 * span for splines of order n, i.e. n <= left <= nknots-n-2. For bsplvb(),
 * `left' must be the index of the nonzero 0-th order spline, i.e.
 * knots[left] <= x < knots[left+1].
 *
 * See Chapter X in:
 *
 * Carl de Boor. A Practical Guide to Splines, volume 27 of Applied
 *     Mathematical Sciences. Springer-Verlag, 1978.
 */

void
bsplvb_simple(const double *knots, const unsigned nknots,
    double x, int left, int degree, float *restrict biatx)
{
	assert(degree>0);
	int i, j;
	double saved, term;
	double delta_l[degree], delta_r[degree];

	biatx[0] = 1.0;

	/*
	 * Handle the (rare) cases where x is outside the full
	 * support of the spline surface.
	 */
	if (left == degree-1)
		while (left >= 0 && x < knots[left])
			left--;
	else if (left == nknots-degree-1)
		while (left < nknots-1 && x > knots[left+1])
			left++;

	/*
	 * NB: if left < degree-1 or left > nknots-degree-1,
	 * the following loop will dereference addresses ouside
	 * of knots[0:nknots]. While terms involving invalid knot
	 * indices will be discarded, it is important that `knots'
	 * have (maxdegree-1)*sizeof(double) bytes of padding
	 * before and after its valid range to prevent segfaults
	 * (see parsefitstable()).
	 */
	for (j = 0; j < degree-1; j++) {
		delta_r[j] = knots[left+j+1] - x;
		delta_l[j] = x - knots[left-j];

		saved = 0.0;

		for (i = 0; i < j+1; i++) {
			term = biatx[i] / (delta_r[i] + delta_l[j-i]);
			biatx[i] = saved + delta_r[i]*term;
			saved = delta_l[j-i]*term;
		}

		biatx[j+1] = saved;
	}

	/*
	 * If left < (spline order), only the first (left+1)
	 * splines are valid; the remainder are utter nonsense.
	 */
	if ((i = degree-1-left) > 0) {
		for (j = 0; j < left+1; j++)
			biatx[j] = biatx[j+i]; /* Move valid splines over. */
		for ( ; j < degree; j++)
			biatx[j] = 0.0; /* The rest are zero by construction. */
	} else if ((i = left+degree+1-nknots) > 0) {
		for (j = degree-1; j > i-1; j--)
			biatx[j] = biatx[j-i];
		for ( ; j >= 0; j--)
			biatx[j] = 0.0;
	}
}

void
bsplvb(const double *knots, const double x, const int left, const int jlow,
    const int jhigh, float *restrict biatx,
    double *restrict delta_l, double *restrict delta_r)
{
	int i, j;
	double saved, term;

	if (jlow == 0)
		biatx[0] = 1.0;

	for (j = jlow; j < jhigh-1; j++) {
		delta_r[j] = knots[left+j+1] - x;
		delta_l[j] = x - knots[left-j];

		saved = 0.0;

		for (i = 0; i < j+1; i++) {
			term = biatx[i] / (delta_r[i] + delta_l[j-i]);
			biatx[i] = saved + delta_r[i]*term;
			saved = delta_l[j-i]*term;
		}

		biatx[j+1] = saved;
	}
}

void
bspline_deriv_nonzero(const double *knots, const unsigned nknots,
    const double x, int left, const int n, float *restrict biatx)
{
	/* NB: it might be tempting to use unsigned integers *left* and *n* here,
	   but indices into the knot vector may be negative (up to -order) before
	   the first fully-supported knot. */
	assert(n>0);
	int i, j;
	double temp, a;
	double delta_l[n], delta_r[n];

	/* Special case for constant splines */
	if (n == 0)
		return;

	/*
	 * Handle the (rare) cases where x is outside the full
	 * support of the spline surface.
	 */
	if (left == n)
		while (left >= 0 && x < knots[left])
			left--;
	else if (left == nknots-n-2)
		while (left < nknots-1 && x > knots[left+1])
			left++;

	/* Get the non-zero n-1th order B-splines at x */
	bsplvb(knots, x, left, 0 /* jlow */, n /* jhigh */,
	    biatx, delta_l, delta_r);

	/*
	 * Now, form the derivatives of the nth order B-splines from
	 * linear combinations of the lower-order splines.
	 */

	/*
	 * On the last supported segment of the ith nth order spline,
	 * only the i+1th n-1th order spline is nonzero.
	 */
	temp = biatx[0];
	biatx[0] =  - n*temp / ((knots[left+1] - knots[left+1-n]));

	/* On the middle segments, both the ith and i+1th splines contribute. */
	for (i = 1; i < n; i++) {
		a = n*temp/((knots[left+i] - knots[left+i-n]));
		temp = biatx[i];
		biatx[i] = a - n*temp/(knots[left+i+1] - knots[left+i+1-n]);
	}
	/*
	 * On the first supported segment of the i+nth nth order spline,
	 * only the ith n-1th order spline is nonzero.
	 */
	biatx[n] = n*temp/((knots[left+n] - knots[left]));

	/* Rearrange for partially-supported points. */
	if ((i = n-left) > 0) {
		for (j = 0; j < left+1; j++)
			biatx[j] = biatx[j+i]; /* Move valid splines over. */
		for ( ; j < n+1; j++)
			biatx[j] = 0.0; /* The rest are zero by construction. */
	} else if ((i = left+n+2-nknots) > 0) {
		for (j = n; j > i-1; j--)
			biatx[j] = biatx[j-i];
		for ( ; j >= 0; j--)
			biatx[j] = 0.0;
	}

}

double
bspline_deriv(const double *knots, double x, int i, int n, unsigned order)
{
	double result;

	if (n == 0) {
		/*
		 * Special case the 0th order case, where B-Splines
		 * are constant functions from one knot to the next.
		 */

		return 0.0;
	}

	if (order <= 1) {
		result = n * bspline(knots, x, i, n-1) / (knots[i+n] - knots[i]);
		result -= n * bspline(knots, x, i+1, n-1) / (knots[i+n+1] - knots[i+1]);
	} else {
		result = n * bspline_deriv(knots, x, i, n-1, order-1) / (knots[i+n] - knots[i]);
		result -= n * bspline_deriv(knots, x, i+1, n-1, order-1) / (knots[i+n+1] - knots[i+1]);
	}
	return result;
}

double
bspline_deriv_2(const double *knots, double x, int i, int n)
{
	double result;

	if (n <= 1) {
		/*
		 * Special case the 1st order case, where B-Splines
		 * are linear functions from one knot to the next.
		 */

		return 0.0;
	}

	result = bspline(knots, x, i, n-2) /
	    ((knots[i+n] - knots[i])*(knots[i+n-1] - knots[i]));
	result -= bspline(knots, x, i+1, n-2) *
	    (1./(knots[i+n] - knots[i]) + 1./(knots[i+n+1] - knots[i+1])) /
	    (knots[i+n] - knots[i+1]);
	result += bspline(knots, x, i+2, n-2) /
	    ((knots[i+n+1] - knots[i+1])*(knots[i+n+1] - knots[i+2]));

	result *= n*(n-1);

	return result;
}

/*
 * Evaluates the results of a full spline basis given a set of knots,
 * a position, an order, and a central spline for the position (or -1).
 * The central spline should be the index of the 0th order basis spline
 * that is non-zero at the position x.
 */

double
splineeval(const double *knots, const double *weights, int nknots, double x, int order,
    int center)
{
	double work = 0.0;
	int i;

	if (center < 0) {
		/* XXX: should be a binary search */
		for (center = 0; center+1 < nknots; center++) {
			if (x > knots[center] && x < knots[center+1])
				break;
		}

		if (center+1 >= nknots)
			return 0.0;
	}

	i = center - order;
	if (i < 0)
		i = 0;

	while (i < nknots-order-1 && i <= center) {
		work += weights[i]*bspline(knots, x, i, order);
		i++;
	}

	return work;
}

int
tablesearchcenters(const struct splinetable *table, const double *x, int *centers)
{
	int i, min, max;

	for (i = 0; i < table->ndim; i++) {

		/* Ensure we are actually inside the table. */
		if (x[i] <= table->knots[i][0] ||
		    x[i] > table->knots[i][table->nknots[i]-1])
			return (-1);

		/*
		 * If we're only a few knots in, take the center to be
		 * the nearest fully-supported knot.
		 */
		if (x[i] < table->knots[i][table->order[i]]) {
			centers[i] = table->order[i];
			continue;
		} else if (x[i] >= table->knots[i][table->naxes[i]]) {
			centers[i] = table->naxes[i]-1;
			continue;
		}

		min = table->order[i];
		max = table->nknots[i]-2;
		do {
			centers[i] = (max+min)/2;

			if (x[i] < table->knots[i][centers[i]])
				max = centers[i]-1;
			else
				min = centers[i]+1;
		} while (x[i] < table->knots[i][centers[i]] ||
		    x[i] >= table->knots[i][centers[i]+1]);

		/*
		 * B-splines are defined on a half-open interval. For the
		 * last point of the interval, move center one point to the
		 * left to get the limit of the sum without evaluating
		 * absent basis functions.
		 */
		if (centers[i] == table->naxes[i])
			centers[i]--;
	}

	return (0);
}

static int
maxorder(int *order, int ndim)
{
	int i, max = 0;

	for (i = 0; i < ndim; i++)
		if (order[i] > max)
			max = order[i];

	return (max);
}

/*
 * The N-Dimensional tensor product basis version of splineeval.
 * Evaluates the results of a full spline basis given a set of knots,
 * a position, an order, and a central spline for the position (or -1).
 * The central spline should be the index of the 0th order basis spline
 * that is non-zero at the position x.
 *
 * x is the vector at which we will evaluate the space
 */

static double
ndsplineeval_core(const struct splinetable *table, const int *centers, int maxdegree,
    float localbasis[table->ndim][maxdegree])
{
	assert(table->ndim>0);
	int i, j, n, tablepos;
	float result;
	float basis_tree[table->ndim+1];
	int nchunks;
	int decomposedposition[table->ndim];

	tablepos = 0;
	for (n = 0; n < table->ndim; n++) {
		decomposedposition[n] = 0;
		tablepos += (centers[n] - table->order[n])*table->strides[n];
	}

	basis_tree[0] = 1;
	for (n = 0; n < table->ndim; n++)
		basis_tree[n+1] = basis_tree[n]*localbasis[n][0];
	nchunks = 1;
	for (n = 0; n < table->ndim - 1; n++)
		nchunks *= (table->order[n] + 1);

	result = 0;
	n = 0;
	while (1) {
		for (i = 0; __builtin_expect(i < table->order[table->ndim-1] +
		    1, 1); i++) {

			result += basis_tree[table->ndim-1]*
			    localbasis[table->ndim-1][i]*
			    table->coefficients[tablepos + i];
		}

		if (__builtin_expect(++n == nchunks, 0))
			break;

		tablepos += table->strides[table->ndim-2];
		decomposedposition[table->ndim-2]++;

		/* Carry to higher dimensions */
		for (i = table->ndim-2;
		    decomposedposition[i] > table->order[i]; i--) {
			decomposedposition[i-1]++;
			tablepos += (table->strides[i-1]
			    - decomposedposition[i]*table->strides[i]);
			decomposedposition[i] = 0;
		}
		for (j = i; __builtin_expect(j < table->ndim-1, 1); j++)
			basis_tree[j+1] = basis_tree[j]*
			    localbasis[j][decomposedposition[j]];
	}

	return result;
}

/* This function returns bspline coefficients along a given dimension, fixing the values+coefficients for the
other dimensions. Used to obtain a 1-d spline representation that can be easily further convolved. */

void
ndsplineeval_slice_coeffs(const struct splinetable *table, const double *x, const int *centers, double *results,int slice_dimension, int derivative, int area_norm )
{
	assert(table->ndim>0);
	int n;
	int maxdegree = maxorder(table->order, table->ndim) + 1;
	float localbasis[table->ndim][maxdegree];

	if(slice_dimension==table->ndim-1)
	{
		printf("ERROR!!! slice dimension cannot be ndim-1 in this implementation!");
	}

	for (n = 0; n < table->ndim; n++) {

			bsplvb_simple(table->knots[n], table->nknots[n],
			    x[n], centers[n], table->order[n] + 1,
			    localbasis[n]);
	}

	int i, j, tablepos, slice_dimension_stride, num_coeffs;
	//double result;
	num_coeffs=table->nknots[slice_dimension]-table->order[slice_dimension]-1;

	memset(results, 0, sizeof(double)*num_coeffs);
	// temp_result is stored to teomporarily store coefficients for derivative calculation
	double temp_result[num_coeffs];

	float basis_tree[table->ndim+1]; // the last basis_tree dimension is unused ..
	int nchunks;
	int decomposedposition[table->ndim];

	//int derivative_correction[table->ndim];

	nchunks = 1;
	for (n = 0; n < table->ndim - 1; n++)
	{
		//derivative_correction[n]=0;

		if(n==slice_dimension)
		{
			continue;
		}
		nchunks *= (table->order[n] + 1);

	}

	// initialize overall table position
	tablepos = 0;
	for (n = 0; n < table->ndim; n++) {
		decomposedposition[n]=0;
		if(n!=slice_dimension)
		{

			tablepos += (centers[n] - table->order[n])*table->strides[n];
		}
		else
		{
			slice_dimension_stride=table->strides[n];
		}
	}

	// initialize local basis
	basis_tree[0] = 1;
	for (n = 0; n < table->ndim; n++)
	{
		if(n==slice_dimension)
		{
			basis_tree[n+1] = basis_tree[n];
			continue;

		}
		basis_tree[n+1] = basis_tree[n]*localbasis[n][0];
	}

	double temp_base_eval=0.0;
	n=0;
	while (1) {

		for (i = 0; __builtin_expect(i < table->order[table->ndim-1] +
		    1, 1); i++) {

			// in contrast to ndsplineeval_core, save a value for each coefficient, separated by the slice dimension stride
			temp_base_eval=basis_tree[table->ndim-1]*localbasis[table->ndim-1][i];
			for(int nc=0; __builtin_expect(nc<num_coeffs,1) ;nc++)
			{
				results[nc] += temp_base_eval*table->coefficients[tablepos + i + nc*slice_dimension_stride];
			}

		}

		if (__builtin_expect(++n == nchunks, 0))
			break;

		// special case when the slicing dimension is ndim-2 .. skip over this dimension immediately by pushing tablepos forward by (order+1)*strides
		if(slice_dimension==table->ndim-2)
		{
			tablepos += table->strides[table->ndim-2]*(table->order[table->ndim-2]+1);
			decomposedposition[table->ndim-2]+=(table->order[table->ndim-2]+1);
		}
		else
		{
			// otherwise just add one stride
			tablepos += table->strides[table->ndim-2];
			decomposedposition[table->ndim-2]++;
		}

		// Carry to higher dimensions
		for (i = table->ndim-2;
		    decomposedposition[i] > table->order[i]; i--) {

			if(i==slice_dimension+1)
			{
				//printf("i=slicedimsino +1 .. \n");
				decomposedposition[i-2]++;
				tablepos += (table->strides[i-2]
			    - decomposedposition[i]*table->strides[i]);
				decomposedposition[i] = 0;
				// add one extra -1, since we want to skip the slicing dimension
				i=i-1;
			}
			else
			{
				decomposedposition[i-1]++;
				tablepos += (table->strides[i-1]
			    - decomposedposition[i]*table->strides[i]);
				decomposedposition[i] = 0;
			}

		}

		// stacks the tree basis up .. never include the dimension of interest, ie.e the slice dimension
		for (j = i; __builtin_expect(j < table->ndim-1, 1); j++)
		 {
			//printf("last loop index .. %d\n", j);
			if(j==slice_dimension)
			{
				basis_tree[j+1] = basis_tree[j];
			}
			else
			{
				basis_tree[j+1] = basis_tree[j]*
			    localbasis[j][decomposedposition[j]];
			}
		}

	}

	for(int nc=0; __builtin_expect(nc<num_coeffs,1) ;nc++)
	{

		if(derivative>0)
		{
			double y_diff;
			double x_diff;
			int deriv_order=table->order[slice_dimension];

			if(nc==0)
			{
				// first one is special
				y_diff=deriv_order*results[0];
			}
			else
			{
				// first form derivatives
				// requires two coefficient results usually ... so only start once we know >=2 coefficients
				// since table->order is really the degree, the new derivative order is just the degree of the original spline

				y_diff=deriv_order*(results[nc]-results[nc-1]);

			}
			x_diff=table->knots[slice_dimension][deriv_order+nc]-table->knots[slice_dimension][nc];

			temp_result[nc]=y_diff/((double)x_diff);

			if(area_norm)
			{
				double norm_factor= ((double)deriv_order)/(   table->knots[slice_dimension][deriv_order+nc] - table->knots[slice_dimension][nc]);
				temp_result[nc]/=norm_factor;
			}

			// check again for area normalization
		}
		else
		{
			if(area_norm)
			{
				int real_order=table->order[slice_dimension]+1;
				double norm_factor= ((double)(real_order))/ (   table->knots[slice_dimension][real_order+nc] - table->knots[slice_dimension][nc]);
				results[nc]/=norm_factor;
			}

			//printf("SLICE: abs res: %d %.10f\n", nc, results[nc]);
		}

	}

	// calculated the derivate .. copy temp_result into result
	if(derivative>0)
	{
		memcpy(results, temp_result, sizeof(temp_result));
	}

}

double
ndsplineeval(const struct splinetable *table, const double *x, const int *centers,
    int derivatives)
{
	assert(table->ndim>0);
	int n;
	int maxdegree = maxorder(table->order, table->ndim) + 1;
	float localbasis[table->ndim][maxdegree];

	for (n = 0; n < table->ndim; n++) {
		if (derivatives & (1 << n)) {
			bspline_deriv_nonzero(table->knots[n],
			    table->nknots[n], x[n], centers[n],
			    table->order[n], localbasis[n]);
		} else {
			bsplvb_simple(table->knots[n], table->nknots[n],
			    x[n], centers[n], table->order[n] + 1,
			    localbasis[n]);
		}
	}

	return ndsplineeval_core(table, centers, maxdegree, localbasis);
}

double
ndsplineeval_deriv(const struct splinetable *table, const double *x,
    const int *centers, const unsigned *derivatives)
{

	assert(table->ndim>0);
	int i, n;
	int maxdegree = maxorder(table->order, table->ndim) + 1;
	float localbasis[table->ndim][maxdegree];

	for (n = 0; n < table->ndim; n++) {
		if (derivatives == NULL || derivatives[n] == 0) {
			bsplvb_simple(table->knots[n], table->nknots[n],
			    x[n], centers[n], table->order[n] + 1,
			    localbasis[n]);
		} else if (derivatives[n] == 1) {
			bspline_deriv_nonzero(table->knots[n],
			    table->nknots[n], x[n], centers[n],
			    table->order[n], localbasis[n]);
		} else {
			for (i = 0; i <= table->order[n]; i++)
				localbasis[n][i] = bspline_deriv(
				    table->knots[n], x[n],
				    centers[n] - table->order[n] + i,
				    table->order[n], derivatives[n]);
		}
	}

	return ndsplineeval_core(table, centers, maxdegree, localbasis);
}

double
ndsplineeval_linalg(const struct splinetable *table, const double *x,
    const int *centers, int derivatives)
{
	assert(table->ndim>0);
	int totalcoeff, n;
	int coeffstrides[table->ndim];
	gsl_matrix_float *basis1, *basis2, *basis_elem;

	assert(table->ndim > 0);
	coeffstrides[table->ndim - 1] = totalcoeff = 1;
        for (n = table->ndim-1; n >= 0; n--) {
                totalcoeff *= (table->order[n] + 1);
                if (n > 0)
                        coeffstrides[n-1] = totalcoeff;
        }

	float basis1_data[totalcoeff], basis2_data[totalcoeff],
	    elem_data[maxorder(table->order, table->ndim) + 1];
	gsl_matrix_float b1, b2, be;
	basis1 = &b1; basis2 = &b2; basis_elem = &be;
	basis1->data = basis1_data;
	basis2->data = basis2_data;
	basis_elem->data = elem_data;

	/*
	 * Form outer product basis1 = basis2 x basis_elem, filling basis_elem
	 * every time with the non-zero basis functions on each axis and
	 * swapping basis1 and basis2 via tmp_basis.
	 */
	basis2->size1 = 1;
	basis2->size2 = 1;
	basis2->data[0] = 1.0;

	for (n = table->ndim-1; n >= 0; n--) {
		gsl_matrix_float *tmp_basis;
		if (derivatives & (1 << n)) {
			bspline_deriv_nonzero(table->knots[n],
			    table->nknots[n], x[n], centers[n],
			    table->order[n], basis_elem->data);
		} else {
			bsplvb_simple(table->knots[n], table->nknots[n],
			    x[n], centers[n], table->order[n] + 1,
			    basis_elem->data);
		}

		basis_elem->size1 = table->order[n] + 1;
		basis_elem->size2 = 1;

		basis1->size2 = basis2->size2;
		basis1->size1 = basis_elem->size1;
		cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
		    basis1->size1, basis1->size2, basis_elem->size2, 1.0,
		    basis_elem->data, 1, basis2->data, basis2->size2, 0,
		    basis1->data, basis1->size2);
		basis1->size2 = basis1->size1 * basis1->size2;
		basis1->size1 = 1;

		tmp_basis = basis1;
		basis1 = basis2;
		basis2 = tmp_basis;
	}

	/* Now basis1 is free, so fill it with the spline coefficients */
	int i, tablepos;
	int decomposedposition[table->ndim];
	tablepos = 0;
	for (n = 0; n < table->ndim; n++) {
		decomposedposition[n] = 0;
		tablepos += (centers[n] - table->order[n])*table->strides[n];
	}

	for (i = 0; i < table->order[table->ndim-1] + 1; i++)
		basis1->data[i] = table->coefficients[tablepos + i];
	for (n = 1; n < coeffstrides[0] /* number of chunks */; n++) {
		tablepos += table->strides[table->ndim-2];
		decomposedposition[table->ndim-2]++;
		/* Carry to higher dimensions */
		for (i = table->ndim-2; __builtin_expect(i > 0 &&
		    decomposedposition[i] > table->order[i], 0); i--) {
			decomposedposition[i-1]++;
			tablepos += (table->strides[i-1]
			    - decomposedposition[i]*table->strides[i]);
			decomposedposition[i] = 0;
		}

		for (i = 0; i < table->order[table->ndim-1] + 1; i++)
			basis1->data[n*(table->order[table->ndim-1] + 1) + i] =
			    table->coefficients[tablepos + i];
	}

	/* Take the dot product */
	__builtin_prefetch(basis1->data);
	__builtin_prefetch(basis2->data);
	return cblas_sdot(totalcoeff, basis1->data, 1, basis2->data, 1);
}