text:
/* matrix/gsl_matrix_ulong.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_MATRIX_ULONG_H__
#define __GSL_MATRIX_ULONG_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_ulong.h>
#include <gsl/gsl_blas_types.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  unsigned long * data;
  gsl_block_ulong * block;
  int owner;
} gsl_matrix_ulong;

typedef struct
{
  gsl_matrix_ulong matrix;
} _gsl_matrix_ulong_view;

typedef _gsl_matrix_ulong_view gsl_matrix_ulong_view;

typedef struct
{
  gsl_matrix_ulong matrix;
} _gsl_matrix_ulong_const_view;

typedef const _gsl_matrix_ulong_const_view gsl_matrix_ulong_const_view;

/* Allocation */

GSL_FUN gsl_matrix_ulong *
gsl_matrix_ulong_alloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ulong *
gsl_matrix_ulong_calloc (const size_t n1, const size_t n2);

GSL_FUN gsl_matrix_ulong *
gsl_matrix_ulong_alloc_from_block (gsl_block_ulong * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_FUN gsl_matrix_ulong *
gsl_matrix_ulong_alloc_from_matrix (gsl_matrix_ulong * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN gsl_vector_ulong *
gsl_vector_ulong_alloc_row_from_matrix (gsl_matrix_ulong * m,
                                        const size_t i);

GSL_FUN gsl_vector_ulong *
gsl_vector_ulong_alloc_col_from_matrix (gsl_matrix_ulong * m,
                                        const size_t j);

GSL_FUN void gsl_matrix_ulong_free (gsl_matrix_ulong * m);

/* Views */

GSL_FUN _gsl_matrix_ulong_view
gsl_matrix_ulong_submatrix (gsl_matrix_ulong * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_row (gsl_matrix_ulong * m, const size_t i);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_column (gsl_matrix_ulong * m, const size_t j);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_diagonal (gsl_matrix_ulong * m);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_subdiagonal (gsl_matrix_ulong * m, const size_t k);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_superdiagonal (gsl_matrix_ulong * m, const size_t k);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_subrow (gsl_matrix_ulong * m, const size_t i,
                         const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_matrix_ulong_subcolumn (gsl_matrix_ulong * m, const size_t j,
                            const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ulong_view
gsl_matrix_ulong_view_array (unsigned long * base,
                             const size_t n1,
                             const size_t n2);

GSL_FUN _gsl_matrix_ulong_view
gsl_matrix_ulong_view_array_with_tda (unsigned long * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_FUN _gsl_matrix_ulong_view
gsl_matrix_ulong_view_vector (gsl_vector_ulong * v,
                              const size_t n1,
                              const size_t n2);

GSL_FUN _gsl_matrix_ulong_view
gsl_matrix_ulong_view_vector_with_tda (gsl_vector_ulong * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_FUN _gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_submatrix (const gsl_matrix_ulong * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_row (const gsl_matrix_ulong * m,
                            const size_t i);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_column (const gsl_matrix_ulong * m,
                               const size_t j);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_diagonal (const gsl_matrix_ulong * m);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_subdiagonal (const gsl_matrix_ulong * m,
                                    const size_t k);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_superdiagonal (const gsl_matrix_ulong * m,
                                      const size_t k);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_subrow (const gsl_matrix_ulong * m, const size_t i,
                               const size_t offset, const size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_matrix_ulong_const_subcolumn (const gsl_matrix_ulong * m, const size_t j,
                                  const size_t offset, const size_t n);

GSL_FUN _gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_array (const unsigned long * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_FUN _gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_array_with_tda (const unsigned long * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_FUN _gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_vector (const gsl_vector_ulong * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_FUN _gsl_matrix_ulong_const_view
gsl_matrix_ulong_const_view_vector_with_tda (const gsl_vector_ulong * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_FUN void gsl_matrix_ulong_set_zero (gsl_matrix_ulong * m);
GSL_FUN void gsl_matrix_ulong_set_identity (gsl_matrix_ulong * m);
GSL_FUN void gsl_matrix_ulong_set_all (gsl_matrix_ulong * m, unsigned long x);

GSL_FUN int gsl_matrix_ulong_fread (FILE * stream, gsl_matrix_ulong * m) ;
GSL_FUN int gsl_matrix_ulong_fwrite (FILE * stream, const gsl_matrix_ulong * m) ;
GSL_FUN int gsl_matrix_ulong_fscanf (FILE * stream, gsl_matrix_ulong * m);
GSL_FUN int gsl_matrix_ulong_fprintf (FILE * stream, const gsl_matrix_ulong * m, const char * format);

GSL_FUN int gsl_matrix_ulong_memcpy(gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);
GSL_FUN int gsl_matrix_ulong_swap(gsl_matrix_ulong * m1, gsl_matrix_ulong * m2);
GSL_FUN int gsl_matrix_ulong_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);

GSL_FUN int gsl_matrix_ulong_swap_rows(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ulong_swap_columns(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ulong_swap_rowcol(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_FUN int gsl_matrix_ulong_transpose (gsl_matrix_ulong * m);
GSL_FUN int gsl_matrix_ulong_transpose_memcpy (gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);
GSL_FUN int gsl_matrix_ulong_transpose_tricpy (CBLAS_UPLO_t Uplo_src, CBLAS_DIAG_t Diag, gsl_matrix_ulong * dest, const gsl_matrix_ulong * src);

GSL_FUN unsigned long gsl_matrix_ulong_max (const gsl_matrix_ulong * m);
GSL_FUN unsigned long gsl_matrix_ulong_min (const gsl_matrix_ulong * m);
GSL_FUN void gsl_matrix_ulong_minmax (const gsl_matrix_ulong * m, unsigned long * min_out, unsigned long * max_out);

GSL_FUN void gsl_matrix_ulong_max_index (const gsl_matrix_ulong * m, size_t * imax, size_t *jmax);
GSL_FUN void gsl_matrix_ulong_min_index (const gsl_matrix_ulong * m, size_t * imin, size_t *jmin);
GSL_FUN void gsl_matrix_ulong_minmax_index (const gsl_matrix_ulong * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_FUN int gsl_matrix_ulong_equal (const gsl_matrix_ulong * a, const gsl_matrix_ulong * b);

GSL_FUN int gsl_matrix_ulong_isnull (const gsl_matrix_ulong * m);
GSL_FUN int gsl_matrix_ulong_ispos (const gsl_matrix_ulong * m);
GSL_FUN int gsl_matrix_ulong_isneg (const gsl_matrix_ulong * m);
GSL_FUN int gsl_matrix_ulong_isnonneg (const gsl_matrix_ulong * m);

GSL_FUN unsigned long gsl_matrix_ulong_norm1 (const gsl_matrix_ulong * m);

GSL_FUN int gsl_matrix_ulong_add (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_FUN int gsl_matrix_ulong_sub (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_FUN int gsl_matrix_ulong_mul_elements (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_FUN int gsl_matrix_ulong_div_elements (gsl_matrix_ulong * a, const gsl_matrix_ulong * b);
GSL_FUN int gsl_matrix_ulong_scale (gsl_matrix_ulong * a, const double x);
GSL_FUN int gsl_matrix_ulong_scale_rows (gsl_matrix_ulong * a, const gsl_vector_ulong * x);
GSL_FUN int gsl_matrix_ulong_scale_columns (gsl_matrix_ulong * a, const gsl_vector_ulong * x);
GSL_FUN int gsl_matrix_ulong_add_constant (gsl_matrix_ulong * a, const double x);
GSL_FUN int gsl_matrix_ulong_add_diagonal (gsl_matrix_ulong * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_FUN int gsl_matrix_ulong_get_row(gsl_vector_ulong * v, const gsl_matrix_ulong * m, const size_t i);
GSL_FUN int gsl_matrix_ulong_get_col(gsl_vector_ulong * v, const gsl_matrix_ulong * m, const size_t j);
GSL_FUN int gsl_matrix_ulong_set_row(gsl_matrix_ulong * m, const size_t i, const gsl_vector_ulong * v);
GSL_FUN int gsl_matrix_ulong_set_col(gsl_matrix_ulong * m, const size_t j, const gsl_vector_ulong * v);
/***********************************************************************/

/* inline functions if you are using GCC */

GSL_FUN INLINE_DECL unsigned long   gsl_matrix_ulong_get(const gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL void    gsl_matrix_ulong_set(gsl_matrix_ulong * m, const size_t i, const size_t j, const unsigned long x);
GSL_FUN INLINE_DECL unsigned long * gsl_matrix_ulong_ptr(gsl_matrix_ulong * m, const size_t i, const size_t j);
GSL_FUN INLINE_DECL const unsigned long * gsl_matrix_ulong_const_ptr(const gsl_matrix_ulong * m, const size_t i, const size_t j);

#ifdef HAVE_INLINE
INLINE_FUN
unsigned long
gsl_matrix_ulong_get(const gsl_matrix_ulong * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
        }
    }
#endif
  return m->data[i * m->tda + j] ;
}

INLINE_FUN
void
gsl_matrix_ulong_set(gsl_matrix_ulong * m, const size_t i, const size_t j, const unsigned long x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  m->data[i * m->tda + j] = x ;
}

INLINE_FUN
unsigned long *
gsl_matrix_ulong_ptr(gsl_matrix_ulong * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (unsigned long *) (m->data + (i * m->tda + j)) ;
}

INLINE_FUN
const unsigned long *
gsl_matrix_ulong_const_ptr(const gsl_matrix_ulong * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(1))
    {
      if (i >= m->size1)
        {
          GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
        }
      else if (j >= m->size2)
        {
          GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
        }
    }
#endif
  return (const unsigned long *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_ULONG_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
// @(#)root/matrix:$Id$
// 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

#include "TMatrixTBase.h"
#include "TMatrixTUtils.h"

#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 TMatrixTSparse<Element> &b,Int_t constr=0);
   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 };
   enum EMatrixCreatorsOp2 { kMult,kMultTranspose,kPlus,kMinus };

   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);
   TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT      <Element> &b);
   TMatrixTSparse(const TMatrixT      <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; }

           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);
           TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixT      <Element> &a,const TMatrixTSparse<Element> &b);
           TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixT      <Element> &b)
                                              { return SetSparseIndexAB(b,a); }

   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; }
   virtual TMatrixTBase<Element>   &SetMatrixArray  (Int_t nr_nonzeros,Int_t *irow,Int_t *icol,Element *data);
   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);
   const   TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
                                           const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
                                            { return (const TMatrixTSparse<Element>&)
                                                     ((const_cast<TMatrixTSparse<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb,nr_nonzeros,
                                                                                const_cast<Int_t *>(pRowIndex),
                                                                                const_cast<Int_t *>(pColIndex),
                                                                                const_cast<Element *>(pData))); }
           TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
                                           Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
   const   TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
                                           const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const;
           TMatrixTSparse<Element> &Use   (TMatrixTSparse<Element> &a);
   const   TMatrixTSparse<Element> &Use   (const TMatrixTSparse<Element> &a) const;

   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)
   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); Clear();
                                                                                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); Clear();
                                                                                APlusB(tmp,source,1); return *this; }
   TMatrixTSparse<Element> &operator-=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
                                                                                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); Clear();
                                                                                AMinusB(tmp,source,1); return *this; }
   TMatrixTSparse<Element> &operator*=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
                                                                                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); Clear();
                                                                                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
};

#ifndef __CINT__
// When building with -fmodules, it instantiates all pending instantiations,
// instead of delaying them until the end of the translation unit.
// We 'got away with' probably because the use and the definition of the
// explicit specialization do not occur in the same TU.
//
// In case we are building with -fmodules, we need to forward declare the
// specialization in order to compile the dictionary G__Matrix.cxx.
template <> TClass *TMatrixTSparse<double>::Class();
#endif // __CINT__

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 const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
                                                                      const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
                                                                        { 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 const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (const TMatrixTSparse<Element> &a) const
                                                                        { 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> 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

ext: h


text:
/**
 * \file FIRFilter.h
 */

#ifndef ATK_EQ_FIRFILTER_H
#define ATK_EQ_FIRFILTER_H

#include <ATK/config.h>
#include <ATK/EQ/config.h>

#include <gsl/gsl>

#include <cassert>
#include <cstdint>
#include <vector>

namespace ATK
{
  /**
   * FIR filter template class
   */
  template<class Coefficients >
  class ATK_EQ_EXPORT FIRFilter final : public Coefficients
  {
  protected:
    using Parent = Coefficients;
    using typename Parent::DataType;
    using typename Parent::AlignedScalarVector;

    using Parent::converted_inputs;
    using Parent::outputs;
    using Parent::coefficients_in;
    using Parent::input_sampling_rate;
    using Parent::output_sampling_rate;
    using Parent::nb_input_ports;
    using Parent::nb_output_ports;

    using Parent::in_order;
    using Parent::input_delay;
    using Parent::setup;

  public:
    /*!
     * @brief Constructor
     * @param nb_channels is the number of input and output channels
     */
    explicit FIRFilter(gsl::index nb_channels = 1)
      :Parent(nb_channels)
    {
    }

    /// Move constructor
    FIRFilter(FIRFilter&& other)
    :Parent(std::move(other))
    {
    }

    void setup() override
    {
      Parent::setup();
      input_delay = in_order;
    }

    void process_impl(gsl::index size) const final
    {
      assert(input_sampling_rate == output_sampling_rate);
      assert(nb_input_ports == nb_output_ports);
      assert(coefficients_in.data());

      const auto* ATK_RESTRICT coefficients_in_ptr = coefficients_in.data();

      for (gsl::index channel = 0; channel < nb_input_ports; ++channel)
      {
        const DataType* ATK_RESTRICT input = converted_inputs[channel] - static_cast<int64_t>(in_order);
        DataType* ATK_RESTRICT output = outputs[channel];

        for (gsl::index i = 0; i < size; ++i)
        {
          output[i] = 0;
        }

        for (gsl::index j = 0; j < in_order + 1; ++j)
        {
          for (gsl::index i = 0; i < size; ++i)
          {
            output[i] += coefficients_in_ptr[j] * input[i + j];
          }
        }
      }

    }

    /// Returns the vector of internal coefficients for the MA section
    const AlignedScalarVector& get_coefficients_in() const
    {
      return coefficients_in;
    }
  };

}

#endif

max_stars_repo_name: AudioTK/AudioTK

max_issues_repo_name: AudioTK/AudioTK

max_forks_repo_name: AudioTK/AudioTK


text:
#ifndef __MATH_GSL_MATRIX__
#define __MATH_GSL_MATRIX__

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

class GslMatrixItem {
 public:
  GslMatrixItem(gsl_matrix* ptr, size_t index1, size_t index2) :
    ptr_(ptr),
    index1_(index1),
    index2_(index2) { }

  operator const double() {
    return gsl_matrix_get(ptr_, index1_, index2_);
  }

  double operator =(const double v) {
    gsl_matrix_set(ptr_, index1_, index2_, v);
    return v;
  }

  double operator +=(const double v) {
    double old_v = gsl_matrix_get(ptr_, index1_, index2_);
    gsl_matrix_set(ptr_, index1_, index2_, v + old_v);
    return v + old_v;
  }
 private:
  gsl_matrix* ptr_;
  size_t index1_;
  size_t index2_;
};

class GslMatrixBase {
 public:
  GslMatrixBase& operator=(const double v) {
    if (v == 0.0) {
      SetZero();
    } else {
      SetAll(v);
    }
    return *this;
  }

  GslMatrixItem operator()(const size_t index1, const size_t index2) const {
    assert(ptr_ != NULL);
    return GslMatrixItem(ptr_, index1, index2);
  }

  void SetZero() {
    assert(ptr_ != NULL);
    gsl_matrix_set_zero(ptr_);
  }

  void SetAll(const double v) {
    assert(ptr_ != NULL);
    gsl_matrix_set_all(ptr_, v);
  }

  void Reset(gsl_matrix* val) {
    if(ptr_ != NULL) {
      gsl_matrix_free(ptr_);
    }
    ptr_ = val;
  }

  int Fprintf(FILE* stream, const char* format) const {
    assert(ptr_ != NULL);
    return gsl_matrix_fprintf(stream, ptr_, format);
  }

  int Fscanf(FILE* stream) {
    assert(ptr_ != NULL);
    return gsl_matrix_fscanf(stream, ptr_);
  }

  void Set(const int i, const int j, double val) {
    gsl_matrix_set(ptr_, i, j, val);
  }

  /*
    double operator()(const int nCol, const int nRow) {
    return gsl_matrix_get(ptr_, nCol, nRow);
    }
  */

  int size1() const {
    return ptr_->size1;
  }

  int size2() const {
    return ptr_->size2;
  }

  double Trace() const {
    double val = 0;
    assert(ptr_ != NULL);
    assert(ptr_->size1 == ptr_->size2);
    for (size_t ii = 0; ii < ptr_->size1; ++ii) {
      val += gsl_matrix_get(ptr_, ii, ii);
    }
    return val;
  }

  double Sum() const {
    double val = 0;
    assert(ptr_ != NULL);
    for (size_t ii = 0; ii < ptr_->size1; ++ii) {
      for (size_t jj = 0; jj < ptr_->size2; ++jj) {
	val += gsl_matrix_get(ptr_, ii, jj);
      }
    }
    return val;
  }

  /*
   * Apply the transpose of this matrix to a vector x and store the result.

   int TransMul(const GslVector& x, GslVector& res, double scale = 0.0) {
   return gsl_blas_dgemv(CblasTrans, 1.0, ptr_, x.ptr(), scale, res.ptr());
   }

   int Mul(const GslVector& x, GslVector& res, double scale = 0.0) {
   return gsl_blas_dgemv(CblasNoTrans, 1.0, ptr_, x.ptr(), scale, res.ptr());
   }
  */

  const gsl_matrix* ptr() const { return ptr_; }
  gsl_matrix* mutable_ptr() { return ptr_; }

 protected:
  GslMatrixBase() : ptr_(NULL) {
  }
  gsl_matrix* ptr_;

 private:
  GslMatrixBase(const GslMatrixBase&) { }
};

class GslMatrix : public GslMatrixBase {
 public:
  GslMatrix(const size_t size1, const size_t size2) : GslMatrixBase() {
    Allocate(size1, size2);
  }

  void Allocate(const size_t size1, const size_t size2) {
    assert(ptr_ == NULL);
    ptr_ = gsl_matrix_alloc(size1, size2);
  }

  GslMatrix() : GslMatrixBase() {
  }

  GslMatrix(gsl_matrix* val) : GslMatrixBase() {
    ptr_ = val;
  }

  ~GslMatrix() {
    if(ptr_ != NULL) {
      gsl_matrix_free(ptr_);
    }
  }

  GslMatrixBase& operator=(const double v) {
    GslMatrixBase::operator=(v);
    return *this;
  }
 private:
  GslMatrix(const GslMatrix&) { }
};

class GslSubmatrix : public GslMatrixBase {
 public:
 GslSubmatrix(GslMatrixBase& matrix, size_t k1, size_t k2, size_t n1, size_t n2) :
  view_(gsl_matrix_submatrix(matrix.mutable_ptr(), k1, k2, n1, n2)) {
     ptr_ = &view_.matrix;
  }

 GslSubmatrix(gsl_matrix* matrix, size_t k1, size_t k2, size_t n1, size_t n2) :
  view_(gsl_matrix_submatrix(matrix, k1, k2, n1, n2)) {
     ptr_ = &view_.matrix;
  }

  GslMatrixBase& operator=(const double v) {
    GslMatrixBase::operator=(v);
    return *this;
  }
 private:
  gsl_matrix_view view_;
  GslSubmatrix(const GslSubmatrix&) { }
};

#endif  // __MATH_GSL_MATRIX__

max_stars_repo_name: boomsbloom/dtm-fmri

max_issues_repo_name: boomsbloom/dtm-fmri

max_forks_repo_name: boomsbloom/dtm-fmri


text:
///
/// @file This file contains the Schur reduction experiment.
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <starneig_test_config.h>
#include <starneig/configuration.h>
#include "experiment.h"
#include "solvers.h"
#include "../common/common.h"
#include "../common/threads.h"
#include "../common/parse.h"
#include "../common/init.h"
#include "../common/checks.h"
#include "../common/hooks.h"
#include "../common/local_pencil.h"
#ifdef STARNEIG_ENABLE_MPI
#include "../common/starneig_pencil.h"
#endif
#include "../common/io.h"
#include "../common/crawler.h"
#include "../common/complex_distr.h"
#include "../hessenberg/solvers.h"
#include <starneig/starneig.h>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#ifdef GSL_FOUND
#include <gsl/gsl_randist.h>
#endif

#define DEFLATE         0x1
#define SET_TO_INF      0x2

static int deflate_and_place_infinities_crawler(
    int offset, int size, int m, int n, int count, size_t *lds,
    void **ptrs, void *arg)
{
    int const *sub = arg;

    double *A = ptrs[0];
    size_t ldA = lds[0];

    double *B = NULL; size_t ldB = 0;
    if (ptrs[1] != NULL) {
        B = ptrs[1];
        ldB = lds[1];
    }

    for (int i = 1; i < size; i++) {
        if (sub[offset+i] & DEFLATE)
            A[(i-1)*ldA+i] = 0.0;
        if (B != NULL && sub[offset+i] & SET_TO_INF)
            B[i*ldB+i] = 0.0;
    }

    return size;
}

///
/// @brief Modifies a matrix pencil such that it is decoupled to multiple
/// independent subproblems and adds
///
/// @param[in] cuts
///         The desired number of deflations.
///
/// @param[in] infinities
///         The desired number of infinities.
///
/// @param[in,out] pencil
///         The matrix pencil.
///
static void deflate_and_place_infinities(
    int cuts, int infinities, pencil_t pencil)
{
    int n = GENERIC_MATRIX_N(pencil->mat_a);
    int *sub = malloc(n*sizeof(int));
    memset(sub, 0, n*sizeof(int));

    cuts = MIN(cuts, n-1);
    for (int i = 0; i < cuts; i++) {
        int p = prand() % (n-1) + 1;
        while (sub[p] & DEFLATE)
            p = prand() % (n-1) + 1;
        sub[p] |= DEFLATE;
    }

    for (int i = 0; i < infinities; i++) {
        int p = prand() % (n-1) + 1;
        while (sub[p] & SET_TO_INF)
            p = prand() % (n-1) + 1;
        sub[p] |= SET_TO_INF;
    }

    crawl_matrices(CRAWLER_RW, CRAWLER_DIAG_WINDOW,
        &deflate_and_place_infinities_crawler, sub, 0,
        pencil->mat_a, pencil->mat_b, NULL);

    free(sub);

    free_matrix_descr(pencil->mat_ca);
    pencil->mat_ca = NULL;

    free_matrix_descr(pencil->mat_cb);
    pencil->mat_cb = NULL;

    fill_pencil(pencil);
}

////////////////////////////////////////////////////////////////////////////////

static void random_initializer_print_usage(int argc, char * const *argv)
{
    printf(
        "  --n (num) -- Problem dimension\n"
        "  --generalized -- Generalized problem\n"
        "  --decouple (num) -- Decouple the problem\n"
        "  --set-to-inf (num) -- Place infinities\n"
    );

    init_helper_print_usage("", INIT_HELPER_ALL, argc, argv);
}

static void random_initializer_print_args(int argc, char * const *argv)
{
    printf(" --n %d", read_int("--n", argc, argv, NULL, -1));
    if (read_opt("--generalized", argc, argv, NULL))
        printf(" --generalized");

    printf(" --decouple %d", read_int("--decouple", argc, argv, NULL, 0));
    printf(" --set-to-inf %d", read_int("--set-to-inf", argc, argv, NULL, 0));

    init_helper_print_args("", INIT_HELPER_ALL, argc, argv);
}

static int random_initializer_check_args(
    int argc, char * const *argv, int *argr)
{
    if (read_int("--n", argc, argv, argr, -1) < 1)
        return 1;

    read_opt("--generalized", argc, argv, argr);

    if (read_int("--decouple", argc, argv, argr, 0) < 0)
        return 1;

    if (read_int("--set-to-inf", argc, argv, argr, 0) < 0)
        return 1;

    return init_helper_check_args("", INIT_HELPER_ALL, argc, argv, argr);
}

static struct hook_data_env* random_initializer_init(
    hook_data_format_t format, int argc, char * const *argv)
{
    printf("INIT...\n");

    int n = read_int("--n", argc, argv, NULL, -1);
    int generalized = read_opt("--generalized", argc, argv, NULL);
    int decouple = read_int("--decouple", argc, argv, NULL, 0);
    int set_to_inf = read_int("--set-to-inf", argc, argv, NULL, 0);

    init_helper_t helper = init_helper_init_hook(
        "", format, n, n, PREC_DOUBLE | NUM_REAL, argc, argv);

    struct hook_data_env *env = malloc(sizeof(struct hook_data_env));
    env->format = format;
    env->copy_data = (hook_data_env_copy_t) copy_pencil;
    env->free_data = (hook_data_env_free_t) free_pencil;
    pencil_t data = env->data = init_pencil();

    data->mat_a = generate_random_hessenberg(n, n, helper);
    data->mat_q = generate_random_householder(n, helper);

    if (generalized) {
        data->mat_b = generate_random_uptriag(n, n, helper);
        data->mat_z = generate_random_householder(n, helper);
    }

    if (0 < decouple || 0 < set_to_inf)
        deflate_and_place_infinities(decouple, set_to_inf, data);

    init_helper_free(helper);

    return env;
}

static const struct hook_initializer_t random_initializer = {
    .name = "random",
    .desc = "Generates a random upper Hessenberg matrix",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_MPI
        HOOK_DATA_FORMAT_PENCIL_STARNEIG,
#endif
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &random_initializer_print_usage,
    .print_args = &random_initializer_print_args,
    .check_args = &random_initializer_check_args,
    .init = &random_initializer_init
};

////////////////////////////////////////////////////////////////////////////////

static void known_initializer_print_usage(int argc, char * const *argv)
{
    printf(
        "  --n (num) -- Problem dimension\n"
        "  --generalized -- Generalized problem\n"
        "  --complex-distr (complex distribution) -- 2-by-2 block "
        "distribution module\n"
    );

    init_helper_print_usage("", INIT_HELPER_ALL, argc, argv);
}

static void known_initializer_print_args(int argc, char * const *argv)
{
    printf(" --n %d", read_int("--n", argc, argv, NULL, -1));

    int generalized = read_opt("--generalized", argc, argv, NULL);

    if (generalized)
        printf(" --generalized");

    struct complex_distr const *complex_distr =
        read_complex_distr("--complex-distr", argc, argv, NULL);

    printf(" --complex-distr %s", complex_distr->name);

    if (complex_distr->print_args != NULL)
        complex_distr->print_args(argc, argv);

    init_helper_print_args("", INIT_HELPER_ALL, argc, argv);
}

static int known_initializer_check_args(
    int argc, char * const *argv, int *argr)
{
    if (read_int("--n", argc, argv, argr, -1) < 1)
        return 1;

    read_opt("--generalized", argc, argv, argr);

    struct complex_distr const *complex_distr =
        read_complex_distr("--complex-distr", argc, argv, argr);
    if (complex_distr == NULL) {
        fprintf(stderr, "Invalid 2-by-2 block distribution module.\n");
        return -1;
    }

    if (complex_distr->check_args != NULL) {
        int ret = complex_distr->check_args(argc, argv, argr);
        if (ret)
            return ret;
    }

    return init_helper_check_args("", INIT_HELPER_ALL, argc, argv, argr);
}

static struct hook_data_env* known_initializer_init(
    hook_data_format_t format, int argc, char * const *argv)
{
    printf("INIT...\n");

    int n = read_int("--n", argc, argv, NULL, -1);

    int generalized = read_opt("--generalized", argc, argv, NULL);

    struct complex_distr const *complex_distr =
        read_complex_distr("--complex-distr", argc, argv, NULL);

    init_helper_t helper = init_helper_init_hook(
        "", format, n, n, PREC_DOUBLE | NUM_REAL, argc, argv);

    struct hook_data_env *env = malloc(sizeof(struct hook_data_env));
    env->format = format;
    env->copy_data = (hook_data_env_copy_t) copy_pencil;
    env->free_data = (hook_data_env_free_t) free_pencil;
    pencil_t pencil = env->data = init_pencil();

    double *real, *imag, *beta;
    init_supplementary_known_eigenvalues(n, &real, &imag, &beta, &pencil->supp);

    // generate (generalized) Schur form and multiply with Householder
    // reflectors from both sides

    if (generalized) {
        matrix_t mat_s = generate_random_uptriag(n, n, helper);
        matrix_t mat_t = generate_identity(n, n, helper);

        complex_distr->init(argc, argv, mat_s, mat_t);
        extract_eigenvalues(mat_s, mat_t, real, imag, beta);

        matrix_t mat_q = generate_random_householder(n, helper);
        matrix_t mat_z = generate_random_householder(n, helper);

        mul_QAZT(mat_q, mat_s, mat_z, &pencil->mat_a);
        mul_QAZT(mat_q, mat_t, mat_z, &pencil->mat_b);

        free_matrix_descr(mat_s);
        free_matrix_descr(mat_t);
        free_matrix_descr(mat_q);
        free_matrix_descr(mat_z);
    }
    else {
        matrix_t mat_s = generate_random_uptriag(n, n, helper);

        complex_distr->init(argc, argv, mat_s, NULL);
        extract_eigenvalues(mat_s, NULL, real, imag, beta);

        matrix_t mat_q = generate_random_householder(n, helper);

        mul_QAZT(mat_q, mat_s, mat_q, &pencil->mat_a);

        free_matrix_descr(mat_s);
        free_matrix_descr(mat_q);
    }

    // reduce the dense matrix (pencil) to Hessenberg(-triangular) form

    pencil->mat_q = generate_identity(n, n, helper);
    pencil->mat_ca = copy_matrix_descr(pencil->mat_a);
    if (generalized) {
        pencil->mat_z = generate_identity(n, n, helper);
        pencil->mat_cb = copy_matrix_descr(pencil->mat_b);
    }

#ifdef STARNEIG_ENABLE_BLACS
    if (format == HOOK_DATA_FORMAT_PENCIL_BLACS) {
        starneig_node_init(threads_get_workers(), STARNEIG_USE_ALL,
            STARNEIG_HINT_DM | STARNEIG_FXT_DISABLE);

        if (generalized) {
            starneig_GEP_DM_HessenbergTriangular(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_b),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q),
                STARNEIG_MATRIX_HANDLE(pencil->mat_z));
        }
        else {
            starneig_SEP_DM_Hessenberg(
                STARNEIG_MATRIX_HANDLE(pencil->mat_a),
                STARNEIG_MATRIX_HANDLE(pencil->mat_q));
        }

        starneig_node_finalize();
    }
#endif

    if (format == HOOK_DATA_FORMAT_PENCIL_LOCAL) {
        starneig_node_init(threads_get_workers(), STARNEIG_USE_ALL,
            STARNEIG_HINT_SM | STARNEIG_FXT_DISABLE);

        if (generalized) {
            starneig_GEP_SM_HessenbergTriangular(LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_b), LOCAL_MATRIX_LD(pencil->mat_b),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q),
                LOCAL_MATRIX_PTR(pencil->mat_z), LOCAL_MATRIX_LD(pencil->mat_z)
            );
        }

        else {
            starneig_SEP_SM_Hessenberg(LOCAL_MATRIX_N(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_a), LOCAL_MATRIX_LD(pencil->mat_a),
                LOCAL_MATRIX_PTR(pencil->mat_q), LOCAL_MATRIX_LD(pencil->mat_q));
        }

        starneig_node_finalize();
    }

    init_helper_free(helper);

    return env;
}

static const struct hook_initializer_t known_initializer = {
    .name = "known",
    .desc = "Generates an upper Hessenberg matrix with known eigenvalues",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_MPI
        HOOK_DATA_FORMAT_PENCIL_STARNEIG,
#endif
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &known_initializer_print_usage,
    .print_args = &known_initializer_print_args,
    .check_args = &known_initializer_check_args,
    .init = &known_initializer_init
};

////////////////////////////////////////////////////////////////////////////////

#ifdef GSL_FOUND

static void hessrand_initializer_print_usage(int argc, char * const *argv)
{
    printf(
        "  --n (num) -- Problem dimension\n"
        "  --generalized -- Generalized problem\n"
        "  --decouple (num) -- Decouple the problem\n"
        "  --set-to-inf (num) -- Place infinities\n"
    );

    init_helper_print_usage("", INIT_HELPER_ALL, argc, argv);
}

static void hessrand_initializer_print_args(int argc, char * const *argv)
{
    printf(" --n %d", read_int("--n", argc, argv, NULL, -1));
    if (read_opt("--generalized", argc, argv, NULL))
        printf(" --generalized");

    printf(" --decouple %d", read_int("--decouple", argc, argv, NULL, 0));
    printf(" --set-to-inf %d", read_int("--set-to-inf", argc, argv, NULL, 0));

    init_helper_print_args("", INIT_HELPER_ALL, argc, argv);
}

static int hessrand_initializer_check_args(
    int argc, char * const *argv, int *argr)
{
    if (read_int("--n", argc, argv, argr, -1) < 1)
        return 1;

    read_opt("--generalized", argc, argv, argr);

    if (read_int("--decouple", argc, argv, argr, 0) < 0)
        return 1;

    if (read_int("--set-to-inf", argc, argv, argr, 0) < 0)
        return 1;

    return init_helper_check_args("", INIT_HELPER_ALL, argc, argv, argr);
}

static int hessrand_crawler(
    int offset, int width, int m, int n, int count, size_t *lds,
    void **ptrs, void *arg)
{
    gsl_rng *r = arg;

    {
        double *A = ptrs[0];
        size_t ldA = lds[0];

        for (int i = 0; i < width; i++) {
            for (int j = 0; j < offset+i+1; j++)
                A[i*ldA+j] = gsl_ran_gaussian(r, 1.0);
            if (offset+i+1 < m)
                A[i*ldA+offset+i+1] = sqrt(gsl_ran_chisq(r, n-offset-i-1));
            for (int j = offset+i+2; j < m; j++)
                A[i*ldA+j] = 0.0;
        }
    }

    if (1 < count) {
        double *B = ptrs[1];
        size_t ldB = lds[1];

        for (int i = 0; i < width; i++) {
            for (int j = 0; j < offset+i; j++)
                B[i*ldB+j] = gsl_ran_gaussian(r, 1.0);
            B[i*ldB+offset+i] = sqrt(gsl_ran_chisq(r, offset+i));
            for (int j = offset+i+1; j < m; j++)
                B[i*ldB+j] = 0.0;
        }

        if (offset == 0)
            B[0] = sqrt(gsl_ran_chisq(r, n));
    }

    return width;
}

static struct hook_data_env* hessrand_initializer_init(
    hook_data_format_t format, int argc, char * const *argv)
{
    printf("INIT...\n");

    int n = read_int("--n", argc, argv, NULL, -1);
    int generalized = read_opt("--generalized", argc, argv, NULL);
    int decouple = read_int("--decouple", argc, argv, NULL, 0);
    int set_to_inf = read_int("--set-to-inf", argc, argv, NULL, 0);

    init_helper_t helper = init_helper_init_hook(
        "", format, n, n, PREC_DOUBLE | NUM_REAL, argc, argv);

    struct hook_data_env *env = malloc(sizeof(struct hook_data_env));
    env->format = format;
    env->copy_data = (hook_data_env_copy_t) copy_pencil;
    env->free_data = (hook_data_env_free_t) free_pencil;
    pencil_t data = env->data = init_pencil();

    data->mat_a = init_matrix(n, n, helper);
    if (generalized)
        data->mat_b = init_matrix(n, n, helper);

    gsl_rng_env_setup();

    const gsl_rng_type *T = gsl_rng_default;
    gsl_rng *r = gsl_rng_alloc(T);
    gsl_rng_set(r, prand());

    crawl_matrices(
        CRAWLER_W, CRAWLER_PANEL, &hessrand_crawler, r, 0,
        data->mat_a, data->mat_b, NULL);

    gsl_rng_free(r);

    data->mat_q = generate_random_householder(n, helper);
    if (generalized)
        data->mat_z = generate_random_householder(n, helper);

    if (0 < decouple || 0 < set_to_inf)
        deflate_and_place_infinities(decouple, set_to_inf, data);

    init_helper_free(helper);

    return env;
}

static const struct hook_initializer_t hessrand_initializer = {
    .name = "hessrand",
    .desc = "Generates a hessrand upper Hessenberg matrix",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_MPI
        HOOK_DATA_FORMAT_PENCIL_STARNEIG,
#endif
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &hessrand_initializer_print_usage,
    .print_args = &hessrand_initializer_print_args,
    .check_args = &hessrand_initializer_check_args,
    .init = &hessrand_initializer_init
};

#endif

////////////////////////////////////////////////////////////////////////////////

static void generic_initializer_print_usage(int argc, char * const *argv)
{
    printf(
        "  --left-input (mtx filename) -- Left-hand side matrix input "
        "file name\n"
        "  --right-input (mtx filename) -- Right-hand side matrix input "
        "file name\n"
        "  --input-begin (num) -- First matrix row/column to be read\n"
        "  --input-end (num) -- Last matrix row/column to be read + 1\n"
        "  --n (num) -- Problem dimension\n"
        "  --generalized -- Generate a generalized problem\n"
        "  --decouple (num) -- Decouple the problem\n"
        "  --set-to-inf (num) -- Place infinities\n"
    );
    init_helper_print_usage("", INIT_HELPER_BLACS_PENCIL, argc, argv);
}

static void generic_initializer_print_args(int argc, char * const *argv)
{
    char const *left_input = read_str("--left-input", argc, argv, NULL, NULL);
    char const *right_input = read_str("--right-input", argc, argv, NULL, NULL);

    if (left_input != NULL) {
        printf(" --left-input %s", left_input);
        if (right_input)
            printf(" --right-input %s", right_input);
        int input_begin = read_int("--input-begin", argc, argv, NULL, -1);
        int input_end = read_int("--input-end", argc, argv, NULL, -1);
        if (0 <= input_begin && 0 <= input_end)
            printf(" --input-begin %d --input-end %d", input_begin, input_end);
    }
    else {
        printf(" --n %d", read_int("--n", argc, argv, NULL, -1));
        if (read_opt("--generalized", argc, argv, NULL))
            printf(" --generalized");
    }

    printf(" --decouple %d", read_int("--decouple", argc, argv, NULL, 0));
    printf(" --set-to-inf %d", read_int("--set-to-inf", argc, argv, NULL, 0));

    init_helper_print_args("", INIT_HELPER_BLACS_PENCIL, argc, argv);
}

static int generic_initializer_check_args(
    int argc, char * const *argv, int *argr)
{
    char const *left_input = read_str("--left-input", argc, argv, argr, NULL);
    char const *right_input = read_str("--right-input", argc, argv, argr, NULL);
    if (left_input != NULL) {
        if (access(left_input, R_OK) != 0) {
            fprintf(stderr, "Left-hand side input file does not exists.\n");
            return 1;
        }
        int input_begin = read_int("--input-begin", argc, argv, argr, -1);
        int input_end = read_int("--input-end", argc, argv, argr, -1);
        if (input_begin < 0 && input_end < input_begin)
            return 1;
    }
    else if (right_input != NULL) {
        fprintf(stderr, "Left-hand side input filename is missing.\n");
        return 1;
    }
    else {
        if (read_int("--n", argc, argv, argr, -1) < 1)
            return 1;
        read_opt("--generalized", argc, argv, argr);
    }

    if (read_int("--decouple", argc, argv, argr, 0) < 0)
        return 1;

    if (read_int("--set-to-inf", argc, argv, argr, 0) < 0)
        return 1;

    return init_helper_check_args(
        "", INIT_HELPER_BLACS_PENCIL, argc, argv, argr);
}

static struct hook_data_env* lapack_initializer_init(
    hook_data_format_t format, int argc, char * const *argv)
{
    printf("INIT...\n");

    int n = read_int("--n", argc, argv, NULL, -1);
    char const *left_input = read_str("--left-input", argc, argv, NULL, NULL);
    char const *right_input = read_str("--right-input", argc, argv, NULL, NULL);
    int input_begin = read_int("--input-begin", argc, argv, NULL, -1);
    int input_end = read_int("--input-end", argc, argv, NULL, -1);
    int generalized = read_opt("--generalized", argc, argv, NULL);
    int decouple = read_int("--decouple", argc, argv, NULL, 0);
    int set_to_inf = read_int("--set-to-inf", argc, argv, NULL, 0);

    if (left_input != NULL) {
        int m;
        read_mtx_dimensions_from_file(left_input, &m, &n);
    }

    if (input_begin == -1)
        input_begin = 0;
    if (input_end == -1)
        input_end = n;

    init_helper_t helper = init_helper_init_hook(
        "", format, n, n, PREC_DOUBLE | NUM_REAL, argc, argv);

    struct hook_data_env *env = malloc(sizeof(struct hook_data_env));
    env->format = HOOK_DATA_FORMAT_PENCIL_LOCAL;
    env->copy_data = (hook_data_env_copy_t) copy_pencil;
    env->free_data = (hook_data_env_free_t) free_pencil;
    pencil_t data = env->data = data = init_pencil();

    // initialize A

    if (left_input != NULL)
        data->mat_a = read_mtx_sub_matrix_from_file(
            input_begin, input_end, left_input, helper);
    else
        data->mat_a = generate_random_fullpos(n, n, helper);

    // initialize B

    if (right_input != NULL) {
        generalized = 1;
        data->mat_b = read_mtx_sub_matrix_from_file(
            input_begin, input_end, right_input, helper);
    }
    else if (generalized) {
        data->mat_b = generate_random_fullpos(n, n, helper);
    }

    // initialize Q and Z

    data->mat_q = generate_identity(n, n, helper);
    if (generalized)
        data->mat_z = generate_identity(n, n, helper);

    // prepare for decoupling

    if (decouple < 1) {
        data->mat_ca = copy_matrix_descr(data->mat_a);
        data->mat_cb = copy_matrix_descr(data->mat_b);
    }

    // reduce

    hook_solver_state_t state =
        hessenberg_lapack_solver.prepare(argc, argv, env);
    hessenberg_lapack_solver.run(state);
    hessenberg_lapack_solver.finalize(state, env);

    if (0 < decouple || 0 < set_to_inf)
        deflate_and_place_infinities(decouple, set_to_inf, data);

    init_helper_free(helper);

    printf("INIT FINISHED.\n");

    return env;
}

static struct hook_data_env* starpu_initializer_init(
    hook_data_format_t format, int argc, char * const *argv)
{
    printf("INIT...\n");

    int n = read_int("--n", argc, argv, NULL, -1);
    char const *left_input = read_str("--left-input", argc, argv, NULL, NULL);
    char const *right_input = read_str("--right-input", argc, argv, NULL, NULL);
    int input_begin = read_int("--input-begin", argc, argv, NULL, -1);
    int input_end = read_int("--input-end", argc, argv, NULL, -1);
    int generalized = read_opt("--generalized", argc, argv, NULL);
    int decouple = read_int("--decouple", argc, argv, NULL, 0);
    int set_to_inf = read_int("--set-to-inf", argc, argv, NULL, 0);

    if (left_input != NULL) {
        int m;
        read_mtx_dimensions_from_file(left_input, &m, &n);
    }

    init_helper_t helper =
        init_helper_init_hook("", format, n, n, PREC_DOUBLE | NUM_REAL, argc, argv);

    struct hook_data_env *env = malloc(sizeof(struct hook_data_env));
    env->format = format;
    env->copy_data = (hook_data_env_copy_t) copy_pencil;
    env->free_data = (hook_data_env_free_t) free_pencil;
    pencil_t data = env->data = data = init_pencil();

    // initialize A

    if (left_input != NULL) {
        if (0 <= input_begin)
            data->mat_a = read_mtx_sub_matrix_from_file(
                input_begin, input_end, left_input, helper);
        else
            data->mat_a = read_mtx_matrix_from_file(left_input, helper);
    }
    else {
        data->mat_a = generate_random_fullpos(n, n, helper);
    }

    // initialize B

    if (right_input != NULL) {
        generalized = 1;
        if (0 <= input_begin)
            data->mat_b = read_mtx_sub_matrix_from_file(
                input_begin, input_end, right_input, helper);
        else
            data->mat_b = read_mtx_matrix_from_file(right_input, helper);
    }
    else if (generalized) {
        data->mat_b = generate_random_fullpos(n, n, helper);
    }

    // initialize Q and Z

    data->mat_q = generate_identity(n, n, helper);
    if (generalized)
        data->mat_z = generate_identity(n, n, helper);

    // prepare for decoupling

    if (decouple < 1 && set_to_inf < 1) {
        data->mat_ca = copy_matrix_descr(data->mat_a);
        data->mat_cb = copy_matrix_descr(data->mat_b);
    }

    // reduce

    if (format == HOOK_DATA_FORMAT_PENCIL_LOCAL) {

        starneig_node_init(threads_get_workers(), STARNEIG_USE_ALL,
            STARNEIG_HINT_SM | STARNEIG_FXT_DISABLE);

        if (generalized)
            starneig_GEP_SM_HessenbergTriangular(LOCAL_MATRIX_N(data->mat_a),
                LOCAL_MATRIX_PTR(data->mat_a), LOCAL_MATRIX_LD(data->mat_a),
                LOCAL_MATRIX_PTR(data->mat_b), LOCAL_MATRIX_LD(data->mat_b),
                LOCAL_MATRIX_PTR(data->mat_q), LOCAL_MATRIX_LD(data->mat_q),
                LOCAL_MATRIX_PTR(data->mat_z), LOCAL_MATRIX_LD(data->mat_z));
        else
            starneig_SEP_SM_Hessenberg(LOCAL_MATRIX_N(data->mat_a),
                LOCAL_MATRIX_PTR(data->mat_a), LOCAL_MATRIX_LD(data->mat_a),
                LOCAL_MATRIX_PTR(data->mat_q), LOCAL_MATRIX_LD(data->mat_q));

        starneig_node_finalize();

    }

#ifdef STARNEIG_ENABLE_BLACS
    if (format == HOOK_DATA_FORMAT_PENCIL_BLACS) {

        starneig_node_init(threads_get_workers(), STARNEIG_USE_ALL,
            STARNEIG_HINT_DM | STARNEIG_FXT_DISABLE);

        if (generalized)
            starneig_GEP_DM_HessenbergTriangular(
                STARNEIG_MATRIX_HANDLE(data->mat_a),
                STARNEIG_MATRIX_HANDLE(data->mat_b),
                STARNEIG_MATRIX_HANDLE(data->mat_q),
                STARNEIG_MATRIX_HANDLE(data->mat_z));
        else
            starneig_SEP_DM_Hessenberg(
                STARNEIG_MATRIX_HANDLE(data->mat_a),
                STARNEIG_MATRIX_HANDLE(data->mat_q));

        starneig_node_finalize();
    }
#endif

    if (0 < decouple || 0 < set_to_inf)
        deflate_and_place_infinities(decouple, set_to_inf, data);

    init_helper_free(helper);

    printf("INIT FINISHED.\n");

    return env;
}

static const struct hook_initializer_t lapack_initializer = {
    .name = "lapack",
    .desc =
        "Reduces a matrix pencil to upper Hessenberg / Hessenberg-triangular "
        "form using a LAPACK algorithm",
    .formats = (hook_data_format_t[]) { HOOK_DATA_FORMAT_PENCIL_LOCAL, 0 },
    .print_usage = &generic_initializer_print_usage,
    .print_args = &generic_initializer_print_args,
    .check_args = &generic_initializer_check_args,
    .init = &lapack_initializer_init
};

static const struct hook_initializer_t starpu_initializer = {
    .name = "starneig",
    .desc =
        "Reduces a matrix pencil to upper Hessenberg / Hessenberg-triangular "
        "form using a StarPU based parallel algorithm",
    .formats = (hook_data_format_t[]) {
        HOOK_DATA_FORMAT_PENCIL_LOCAL,
#ifdef STARNEIG_ENABLE_BLACS
        HOOK_DATA_FORMAT_PENCIL_BLACS,
#endif
        0 },
    .print_usage = &generic_initializer_print_usage,
    .print_args = &generic_initializer_print_args,
    .check_args = &generic_initializer_check_args,
    .init = &starpu_initializer_init
};

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////

static void print_usage(int argc, char * const *argv)
{
    print_avail_complex_distr();
    print_opt_complex_distr();
}

const struct hook_experiment_descr schur_experiment = {
    .print_usage = &print_usage,
    .initializers = (struct hook_initializer_t const *[])
    {
#ifdef GSL_FOUND
        &hessrand_initializer,
#endif
        &starpu_initializer,
        &lapack_initializer,
        &random_initializer,
        &known_initializer,
        &mtx_initializer,
        &raw_initializer,
        0
    },
    .supplementers = (struct hook_supplementer_t const *[])
    {
        0
    },
    .solvers = (struct hook_solver const *[])
    {
        &schur_starpu_solver,
        &schur_starpu_simple_solver,
        &schur_lapack_solver,
#ifdef PDLAHQR_FOUND
        &schur_pdlahqr_solver,
#endif
#ifdef PDHSEQR_FOUND
        &schur_pdhseqr_solver,
#endif
#ifdef CUSTOM_PDHSEQR
        &schur_custom_pdhseqr_solver,
#endif
#ifdef PDHGEQZ_FOUND
        &schur_pdhgeqz_solver,
#endif
        0
    },
    .hook_descrs = (struct hook_descr_t const *[])
    {
        &default_schur_test_descr,
        &default_eigenvalues_descr,
        &default_known_eigenvalues_descr,
        &default_analysis_descr,
        &default_residual_test_descr,
        &default_print_pencil_descr,
        &default_print_input_pencil_descr,
        &default_store_raw_pencil_descr,
        &default_store_raw_input_pencil_descr,
        0
    }
};

ext: c

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
/* specfunc/test_sf.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman */

#ifndef TEST_SF_H
#define TEST_SF_H

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_machine.h>
#include <gsl/gsl_sf_result.h>

double test_sf_frac_diff(double x1, double x2);
int test_sf_check_result(char * message_buff, gsl_sf_result r, double val, double tol);
int test_sf_check_val(char * message_buff, double rval, double val, double tol);
int test_sf_check_return(char * message_buff, int val_return, int expected_return);
int test_sf_check_result_relax(char * message_buff, gsl_sf_result r, double val, double tol);

/* Include an overall test factor to allow for differences between
   compilers, otherwise there are too many bug reports on the released
   versions.  Turn this value down to 1.0 for development purposes */

#ifndef TEST_FACTOR
#ifdef RELEASED
#define TEST_FACTOR 100.0
#else
#define TEST_FACTOR 1.0
#endif
#endif

#ifndef TEST_SIGMA
#ifdef RELEASED
#define TEST_SIGMA 1.5
#else
#define TEST_SIGMA 1.0
#endif
#endif

#define TEST_TOL0  (2.0*GSL_DBL_EPSILON)
#define TEST_TOL1  (16.0*GSL_DBL_EPSILON)
#define TEST_TOL2  (256.0*GSL_DBL_EPSILON)
#define TEST_TOL3  (2048.0*GSL_DBL_EPSILON)
#define TEST_TOL4  (16384.0*GSL_DBL_EPSILON)
#define TEST_TOL5  (131072.0*GSL_DBL_EPSILON)
#define TEST_TOL6  (1048576.0*GSL_DBL_EPSILON)
#define TEST_SQRT_TOL0 (2.0*GSL_SQRT_DBL_EPSILON)
#define TEST_SNGL  (1.0e-06)

#define TEST_SF_INCONS  1
#define TEST_SF_ERRNEG  2
#define TEST_SF_TOLBAD  4
#define TEST_SF_RETBAD  8
#define TEST_SF_ERRBAD  16
#define TEST_SF_ERRBIG  32
#define TEST_SF_EXPBAD  64

int test_sf (gsl_sf_result r, double val_in, double tol, int status, int expect_return, const char * desc);
int test_sf_e10 (gsl_sf_result_e10 r, double val_in, int e10_in, double tol, int status, int expect_return, const char * desc);
int test_sf_val (double val, double val_in, double tol, const char * desc);
int test_sf_rlx (gsl_sf_result r, double val_in, double tol, int status, int expect_return, const char * desc);
int test_sf_2 (gsl_sf_result r1, double val1, double tol1, gsl_sf_result r2, double val2, double tol2, int status, int expect_return, const char * desc);
int test_sf_sgn (gsl_sf_result r, double sgn, double val_in, double tol, double expect_sgn, int status, int expect_return, const char * desc);
int test_sf_return (int status, int expect_return, const char * desc);

#define TEST_SF(stat, func, args, val_in, tol, expect_return) { int status = func args; stat += test_sf(r, val_in, tol, status, expect_return, #func #args); }

#define TEST_SF_E10(stat, func, args, val_in, e10_in, tol, expect_return) { int status = func args; stat += test_sf_e10(re, val_in, e10_in, tol, status, expect_return, #func #args); }

#define TEST_SF_VAL(stat, func, args, val_in, tol) { double val = func args; stat += test_sf_val(val, val_in, tol, #func #args); }

#define TEST_SF_RLX(stat, func, args, val_in, tol, expect_return) { int status = func args; stat += test_sf_rlx(r, val_in, tol, status, expect_return, #func #args); }

#define TEST_SF_2(stat, func, args, val1, tol1, val2, tol2, expect_return) { int status = func args; stat += test_sf_2(r1, val1, tol1, r2, val2, tol2, status, expect_return, #func #args); }

#define TEST_SF_SGN(stat, func, args, val_in, tol, expect_sgn, expect_return) { int status = func args; stat += test_sf_sgn(r, sgn, val_in, tol, expect_sgn, status, expect_return, #func #args); }

#define TEST_SF_THETA(stat, func, args, val_in, tol) { int status; theta=args; status = func (&theta);  stat += test_sf_val(theta, val_in, tol, #func #args); }

#define TEST_SF_RETURN(stat, func, args, expect_return) { int status = func args; stat += test_sf_return(status, expect_return, #func #args); }

int test_airy(void);
int test_bessel(void);
int test_coulomb(void);
int test_dilog(void);
int test_gamma(void);
int test_mathieu(void);
int test_hermite(void);
int test_hyperg(void);
int test_legendre(void);
int test_sincos_pi(void);

#endif /* !TEST_SF_H */

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <gsl/gsl_ntuple.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

struct data
{
  double x;
  double y;
  double z;
};

int
main (void)
{
  const gsl_rng_type * T;
  gsl_rng * r;

  struct data ntuple_row;
  int i;

  gsl_ntuple *ntuple
    = gsl_ntuple_create ("test.dat", &ntuple_row,
                         sizeof (ntuple_row));

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for (i = 0; i < 10000; i++)
    {
      ntuple_row.x = gsl_ran_ugaussian (r);
      ntuple_row.y = gsl_ran_ugaussian (r);
      ntuple_row.z = gsl_ran_ugaussian (r);

      gsl_ntuple_write (ntuple);
    }

  gsl_ntuple_close (ntuple);
  gsl_rng_free (r);

  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>

/* This is an implementation of Martin Luescher's second generation
   double-precision (48-bit) version of the RANLUX generator.

   Thanks to Martin Luescher for providing information on this
   generator.

*/

static inline unsigned long int ranlxd_get (void *vstate);
static double ranlxd_get_double (void *vstate);
static void ranlxd_set_lux (void *state, unsigned long int s, unsigned int luxury);
static void ranlxd1_set (void *state, unsigned long int s);
static void ranlxd2_set (void *state, unsigned long int s);

static const int next[12] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0};

static const double one_bit = 1.0 / 281474976710656.0;	/* 1/2^48 */

#define RANLUX_STEP(x1,x2,i1,i2,i3)      \
          x1=xdbl[i1] - xdbl[i2];        \
          if (x2 < 0)                    \
          {                              \
            x1-=one_bit;                 \
            x2+=1;                       \
          }                              \
         xdbl[i3]=x2

typedef struct
  {
    double xdbl[12];
    double carry;
    unsigned int ir;
    unsigned int jr;
    unsigned int ir_old;
    unsigned int pr;
  }
ranlxd_state_t;

static inline void increment_state (ranlxd_state_t * state);

static inline void
increment_state (ranlxd_state_t * state)
{
  int k, kmax;
  double y1, y2, y3;

  double *xdbl = state->xdbl;
  double carry = state->carry;
  unsigned int ir = state->ir;
  unsigned int jr = state->jr;

  for (k = 0; ir > 0; ++k)
    {
      y1 = xdbl[jr] - xdbl[ir];
      y2 = y1 - carry;
      if (y2 < 0)
	{
	  carry = one_bit;
	  y2 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[ir] = y2;
      ir = next[ir];
      jr = next[jr];
    }

  kmax = state->pr - 12;

  for (; k <= kmax; k += 12)
    {
      y1 = xdbl[7] - xdbl[0];
      y1 -= carry;

      RANLUX_STEP (y2, y1, 8, 1, 0);
      RANLUX_STEP (y3, y2, 9, 2, 1);
      RANLUX_STEP (y1, y3, 10, 3, 2);
      RANLUX_STEP (y2, y1, 11, 4, 3);
      RANLUX_STEP (y3, y2, 0, 5, 4);
      RANLUX_STEP (y1, y3, 1, 6, 5);
      RANLUX_STEP (y2, y1, 2, 7, 6);
      RANLUX_STEP (y3, y2, 3, 8, 7);
      RANLUX_STEP (y1, y3, 4, 9, 8);
      RANLUX_STEP (y2, y1, 5, 10, 9);
      RANLUX_STEP (y3, y2, 6, 11, 10);

      if (y3 < 0)
	{
	  carry = one_bit;
	  y3 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[11] = y3;
    }

  kmax = state->pr;

  for (; k < kmax; ++k)
    {
      y1 = xdbl[jr] - xdbl[ir];
      y2 = y1 - carry;
      if (y2 < 0)
	{
	  carry = one_bit;
	  y2 += 1;
	}
      else
	{
	  carry = 0;
	}
      xdbl[ir] = y2;
      ir = next[ir];
      jr = next[jr];
    }
  state->ir = ir;
  state->ir_old = ir;
  state->jr = jr;
  state->carry = carry;
}

static inline unsigned long int
ranlxd_get (void *vstate)
{
  return ranlxd_get_double (vstate) * 4294967296.0;	/* 2^32 */
}

static double
ranlxd_get_double (void *vstate)
{
  ranlxd_state_t *state = (ranlxd_state_t *) vstate;

  int ir = state->ir;

  state->ir = next[ir];

  if (state->ir == state->ir_old)
    increment_state (state);

  return state->xdbl[state->ir];
}

static void
ranlxd_set_lux (void *vstate, unsigned long int s, unsigned int luxury)
{
  ranlxd_state_t *state = (ranlxd_state_t *) vstate;

  int ibit, jbit, i, k, l, xbit[31];
  double x, y;

  long int seed;

  if (s == 0)
    s = 1;			/* default seed is 1 */

  seed = s;

  i = seed & 0xFFFFFFFFUL;

  for (k = 0; k < 31; ++k)
    {
      xbit[k] = i % 2;
      i /= 2;
    }

  ibit = 0;
  jbit = 18;

  for (k = 0; k < 12; ++k)
    {
      x = 0;

      for (l = 1; l <= 48; ++l)
	{
	  y = (double) ((xbit[ibit] + 1) % 2);
	  x += x + y;
	  xbit[ibit] = (xbit[ibit] + xbit[jbit]) % 2;
	  ibit = (ibit + 1) % 31;
	  jbit = (jbit + 1) % 31;
	}
      state->xdbl[k] = one_bit * x;
    }

  state->carry = 0;
  state->ir = 11;
  state->jr = 7;
  state->ir_old = 0;
  state->pr = luxury;
}

static void
ranlxd1_set (void *vstate, unsigned long int s)
{
  ranlxd_set_lux (vstate, s, 202);
}

static void
ranlxd2_set (void *vstate, unsigned long int s)
{
  ranlxd_set_lux (vstate, s, 397);
}

static const gsl_rng_type ranlxd1_type =
{"ranlxd1",			/* name */
 0xffffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlxd_state_t),
 &ranlxd1_set,
 &ranlxd_get,
 &ranlxd_get_double};

static const gsl_rng_type ranlxd2_type =
{"ranlxd2",			/* name */
 0xffffffffUL,			/* RAND_MAX */
 0,				/* RAND_MIN */
 sizeof (ranlxd_state_t),
 &ranlxd2_set,
 &ranlxd_get,
 &ranlxd_get_double};

const gsl_rng_type *gsl_rng_ranlxd1 = &ranlxd1_type;
const gsl_rng_type *gsl_rng_ranlxd2 = &ranlxd2_type;


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

#include <config.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

#include "initialise.c"
#include "append.c"
#include "qpsrt.c"
#include "util.c"
#include "qc25s.c"

int
gsl_integration_qaws (gsl_function * f,
                      const double a, const double b,
                      gsl_integration_qaws_table * t,
                      const double epsabs, const double epsrel,
                      const size_t limit,
                      gsl_integration_workspace * workspace,
                      double *result, double *abserr)
{
  double area, errsum;
  double result0, abserr0;
  double tolerance;
  size_t iteration = 0;
  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;

  /* Initialize results */

  initialise (workspace, a, b);

  *result = 0;
  *abserr = 0;

  if (limit > workspace->limit)
    {
      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
    }

  if (b <= a)
    {
      GSL_ERROR ("limits must form an ascending sequence, a < b", GSL_EINVAL) ;
    }

  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
    {
      GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
                 GSL_EBADTOL);
    }

  /* perform the first integration */

  {
    double area1, area2;
    double error1, error2;
    int err_reliable1, err_reliable2;
    double a1 = a;
    double b1 = 0.5 * (a + b);
    double a2 = b1;
    double b2 = b;

    qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
    qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);

    if (error1 > error2)
      {
        append_interval (workspace, a1, b1, area1, error1);
        append_interval (workspace, a2, b2, area2, error2);
      }
    else
      {
        append_interval (workspace, a2, b2, area2, error2);
        append_interval (workspace, a1, b1, area1, error1);
      }

    result0 = area1 + area2;
    abserr0 = error1 + error2;
  }

  /* Test on accuracy */

  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));

  /* Test on accuracy, use 0.01 relative error as an extra safety
     margin on the first iteration (ignored for subsequent iterations) */

  if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0))
    {
      *result = result0;
      *abserr = abserr0;

      return GSL_SUCCESS;
    }
  else if (limit == 1)
    {
      *result = result0;
      *abserr = abserr0;

      GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
    }

  area = result0;
  errsum = abserr0;

  iteration = 2;

  do
    {
      double a1, b1, a2, b2;
      double a_i, b_i, r_i, e_i;
      double area1 = 0, area2 = 0, area12 = 0;
      double error1 = 0, error2 = 0, error12 = 0;
      int err_reliable1, err_reliable2;

      /* Bisect the subinterval with the largest error estimate */

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      a1 = a_i;
      b1 = 0.5 * (a_i + b_i);
      a2 = b1;
      b2 = b_i;

      qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
      qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);

      area12 = area1 + area2;
      error12 = error1 + error2;

      errsum += (error12 - e_i);
      area += area12 - r_i;

      if (err_reliable1 && err_reliable2)
        {
          double delta = r_i - area12;

          if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
            {
              roundoff_type1++;
            }
          if (iteration >= 10 && error12 > e_i)
            {
              roundoff_type2++;
            }
        }

      tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));

      if (errsum > tolerance)
        {
          if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
            {
              error_type = 2;   /* round off error */
            }

          /* set error flag in the case of bad integrand behaviour at
             a point of the integration range */

          if (subinterval_too_small (a1, a2, b2))
            {
              error_type = 3;
            }
        }

      update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      iteration++;

    }
  while (iteration < limit && !error_type && errsum > tolerance);

  *result = sum_results (workspace);
  *abserr = errsum;

  if (errsum <= tolerance)
    {
      return GSL_SUCCESS;
    }
  else if (error_type == 2)
    {
      GSL_ERROR ("roundoff error prevents tolerance from being achieved",
                 GSL_EROUND);
    }
  else if (error_type == 3)
    {
      GSL_ERROR ("bad integrand behavior found in the integration interval",
                 GSL_ESING);
    }
  else if (iteration == limit)
    {
      GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
    }
  else
    {
      GSL_ERROR ("could not integrate function", GSL_EFAILED);
    }

}

max_stars_repo_name: ruslankuzmin/julia


text:
/*
    Copyright (C) 2008 Wei Dong <wdong@princeton.edu>. All Rights Reserved.

    This file is part of LSHKIT.

    LSHKIT 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.

    LSHKIT 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 LSHKIT.  If not, see <http://www.gnu.org/licenses/>.
*/

/**
 *  \file fitdata.cpp
 *  \brief Gather statistics from dataset for MPLSH tuning.
 *
 *  This program gahters statistical data from a small sample dataset
 *  for automatic MPLSH parameter tuning.  It carries out the following
 *  steps:
 *  -# Sample N points from the dataset. Only those N points will be used for future computation.
 *  -# Sample P pairs of points from the sample, calculate the distance for each pair.
 *  -# Sample Q points from the sample as queries points.
 *  -# Divide the sample into F folds.
 *  -# For i = 1 to F, take i folds and run K-NN search, so the query points
 *     will be searched against sample datasets of N/F, 2N/F, ..., N/F points.
 *
 *  The statistical data is printed to standard output after the progress display.
 *
\verbatim
Allowed options:
  -h [ --help ]          produce help message.
  -N [ -- ] arg (=0)     number of points to use
  -P [ -- ] arg (=50000) number of pairs to sample
  -Q [ -- ] arg (=1000)  number of queries to sample
  -K [ -- ] arg (=100)   search for K nearest neighbors
  -F [ -- ] arg (=10)    divide the sample to F folds
  -D [ --data ] arg      data file
\endverbatim
 */

#ifndef _MPLSH_FITDATA_H_
#define _MPLSH_FITDATA_H_

#include "logging.h"

#include <cstdlib>
#include <sstream>
#include <gsl/gsl_multifit.h>

namespace lshkit {

inline
bool is_good_value (double v) {
    return ((v > -std::numeric_limits<double>::max()) &&
            (v < std::numeric_limits<double>::max()));
}

inline
std::string FitData(const FloatMatrix& data,
                    unsigned N,            // number of points to use
                    unsigned P,            // number of pairs to sample
                    unsigned Q,            // number of queries to sample
                    unsigned K,            // search for K neighbors neighbors
                    unsigned F             // divide the sample to F folds
                   )
{
    LOG(LIB_INFO) << "started running FitData" << std::endl;

    std::vector<unsigned> idx(data.getSize());
    for (size_t i = 0; i < idx.size(); ++i) idx[i] = i;
    std::random_shuffle(idx.begin(), idx.end());

    if (N > 0 && N < static_cast<unsigned>(data.getSize())) idx.resize(N);

    metric::l2sqr<float> l2sqr(data.getDim());

    DefaultRng rng;
    boost::variate_generator<DefaultRng &, UniformUnsigned> gen(rng,
             UniformUnsigned(0, idx.size()-1));

    double gM = 0.0;
    double gG = 0.0;
    {
        // sample P pairs of points
        for (unsigned k = 0; k < P; ++k)
        {
            double dist, logdist;
            for (;;)
            {
                unsigned i = gen();
                unsigned j = gen();
                if (i == j) continue;
                dist = l2sqr(data[idx[i]], data[idx[j]]);
                logdist = log(dist);
                if (is_good_value(logdist)) break;
            }
            gM += dist;
            gG += logdist;
        }
        gM /= P;
        gG /= P;
        gG = exp(gG);
    }

    if (Q > idx.size()) Q = idx.size();
    if (K > idx.size() - Q) K = idx.size() - Q;
    /* sample query */
    std::vector<unsigned> qry(Q);

    SampleQueries(&qry, idx.size(), rng);

    /* do the queries */
    std::vector<Topk<unsigned> > topks(Q);
    for (unsigned i = 0; i < Q; ++i) topks[i].reset(K);

    /* ... */
    gsl_matrix *X = gsl_matrix_alloc(F * K, 3);
    gsl_vector *yM = gsl_vector_alloc(F * K);
    gsl_vector *yG = gsl_vector_alloc(F * K);
    gsl_vector *pM = gsl_vector_alloc(3);
    gsl_vector *pG = gsl_vector_alloc(3);
    gsl_matrix *cov = gsl_matrix_alloc(3,3);

    std::vector<double> M(K);
    std::vector<double> G(K);

    unsigned m = 0;
    for (unsigned l = 0; l < F; l++)
    {
        // Scan
        for (unsigned i = l; i< idx.size(); i += F)
        {
            for (unsigned j = 0; j < Q; j++)
            {
                unsigned id = qry[j];
                if (i != id)
                {
                    float d = l2sqr(data[idx[id]], data[idx[i]]);
                    if (is_good_value(log(double(d))))
                      topks[j] << Topk<unsigned>::Element(i, d);
                }
            }
        }

        fill(M.begin(), M.end(), 0.0);
        fill(G.begin(), G.end(), 0.0);

        for (unsigned i = 0; i < Q; i++)
        {
            for (unsigned k = 0; k < K; k++)
            {
                M[k] += topks[i][k].dist;
                G[k] += log(topks[i][k].dist);
            }
        }

        for (unsigned k = 0; k < K; k++)
        {
            M[k] = log(M[k]/Q);
            G[k] /= Q;
            gsl_matrix_set(X, m, 0, 1.0);
            gsl_matrix_set(X, m, 1, log(double(data.getSize() * (l + 1)) / double(F)));
            gsl_matrix_set(X, m, 2, log(double(k + 1)));
            gsl_vector_set(yM, m, M[k]);
            gsl_vector_set(yG, m, G[k]);
            ++m;
        }

        //++progress;
    }

    gsl_multifit_linear_workspace *work = gsl_multifit_linear_alloc(F * K, 3);

    double chisq;

    gsl_multifit_linear(X, yM, pM, cov, &chisq, work);
    gsl_multifit_linear(X, yG, pG, cov, &chisq, work);

    std::stringstream ss;
    ss << gM << " " << gG << std::endl;
    ss << gsl_vector_get(pM, 0) << " "
       << gsl_vector_get(pM, 1) << " "
       << gsl_vector_get(pM, 2) << std::endl;
    ss << gsl_vector_get(pG, 0) << " "
       << gsl_vector_get(pG, 1) << " "
       << gsl_vector_get(pG, 2) << std::endl;

    gsl_matrix_free(X);
    gsl_matrix_free(cov);
    gsl_vector_free(pM);
    gsl_vector_free(pG);
    gsl_vector_free(yM);
    gsl_vector_free(yG);

    LOG(LIB_INFO) << ss.str();
    LOG(LIB_INFO) << "finished FitData";

    return ss.str();
}

}     // namespace lshkit

#endif