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

#ifndef __GSL_MULTIFIT_H__
#define __GSL_MULTIFIT_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_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_types.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t nmax;         /* maximum number of observations */
  size_t pmax;         /* maximum number of parameters */
  size_t n;            /* number of observations in current SVD decomposition */
  size_t p;            /* number of parameters in current SVD decomposition */
  gsl_matrix * A;      /* least squares matrix for SVD, n-by-p */
  gsl_matrix * Q;
  gsl_matrix * QSI;
  gsl_vector * S;
  gsl_vector * t;
  gsl_vector * xt;
  gsl_vector * D;
  double rcond;        /* reciprocal condition number */
}
gsl_multifit_linear_workspace;

GSL_FUN gsl_multifit_linear_workspace *
gsl_multifit_linear_alloc (const size_t n, const size_t p);

GSL_FUN void
gsl_multifit_linear_free (gsl_multifit_linear_workspace * w);

GSL_FUN int
gsl_multifit_linear (const gsl_matrix * X,
                     const gsl_vector * y,
                     gsl_vector * c,
                     gsl_matrix * cov,
                     double * chisq,
                     gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_tsvd (const gsl_matrix * X,
                          const gsl_vector * y,
                          const double tol,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double * chisq,
                          size_t * rank,
                          gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_svd (const gsl_matrix * X,
                         gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_bsvd (const gsl_matrix * X,
                          gsl_multifit_linear_workspace * work);

GSL_FUN size_t
gsl_multifit_linear_rank(const double tol, const gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_solve (const double lambda,
                           const gsl_matrix * X,
                           const gsl_vector * y,
                           gsl_vector * c,
                           double *rnorm,
                           double *snorm,
                           gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_applyW(const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           gsl_matrix * WX,
                           gsl_vector * Wy);

GSL_FUN int
gsl_multifit_linear_stdform1 (const gsl_vector * L,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_wstdform1 (const gsl_vector * L,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_L_decomp (gsl_matrix * L, gsl_vector * tau);

GSL_FUN int
gsl_multifit_linear_stdform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              gsl_matrix * Xs,
                              gsl_vector * ys,
                              gsl_matrix * M,
                              gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_wstdform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               gsl_matrix * Xs,
                               gsl_vector * ys,
                               gsl_matrix * M,
                               gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_genform1 (const gsl_vector * L,
                              const gsl_vector * cs,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_genform2 (const gsl_matrix * LQR,
                              const gsl_vector * Ltau,
                              const gsl_matrix * X,
                              const gsl_vector * y,
                              const gsl_vector * cs,
                              const gsl_matrix * M,
                              gsl_vector * c,
                              gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_wgenform2 (const gsl_matrix * LQR,
                               const gsl_vector * Ltau,
                               const gsl_matrix * X,
                               const gsl_vector * w,
                               const gsl_vector * y,
                               const gsl_vector * cs,
                               const gsl_matrix * M,
                               gsl_vector * c,
                               gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_lreg (const double smin, const double smax,
                          gsl_vector * reg_param);

GSL_FUN int
gsl_multifit_linear_lcurve (const gsl_vector * y,
                            gsl_vector * reg_param,
                            gsl_vector * rho, gsl_vector * eta,
                            gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_lcurvature (const gsl_vector * y,
                                const gsl_vector * reg_param,
                                const gsl_vector * rho,
                                const gsl_vector * eta,
                                gsl_vector * kappa,
                                gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_lcorner(const gsl_vector *rho,
                            const gsl_vector *eta,
                            size_t *idx);

GSL_FUN int
gsl_multifit_linear_lcorner2(const gsl_vector *reg_param,
                             const gsl_vector *eta,
                             size_t *idx);

GSL_FUN int
gsl_multifit_linear_Lk(const size_t p, const size_t k, gsl_matrix *L);

GSL_FUN int
gsl_multifit_linear_Lsobolev(const size_t p, const size_t kmax,
                             const gsl_vector *alpha, gsl_matrix *L,
                             gsl_multifit_linear_workspace *work);

GSL_FUN int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double * chisq,
                      gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_wlinear_tsvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           const double tol,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double * chisq,
                           size_t * rank,
                           gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq,
                          gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_wlinear_usvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           double tol,
                           size_t * rank,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double *chisq,
                           gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_est (const gsl_vector * x,
                         const gsl_vector * c,
                         const gsl_matrix * cov, double *y, double *y_err);

GSL_FUN double
gsl_multifit_linear_rcond (const gsl_multifit_linear_workspace * w);

GSL_FUN int
gsl_multifit_linear_residuals (const gsl_matrix *X, const gsl_vector *y,
                               const gsl_vector *c, gsl_vector *r);

/* gcv.c */
GSL_FUN int
gsl_multifit_linear_gcv_init(const gsl_vector * y,
                             gsl_vector * reg_param,
                             gsl_vector * UTy,
                             double * delta0,
                             gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_gcv_curve(const gsl_vector * reg_param,
                              const gsl_vector * UTy,
                              const double delta0,
                              gsl_vector * G,
                              gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_gcv_min(const gsl_vector * reg_param,
                            const gsl_vector * UTy,
                            const gsl_vector * G,
                            const double delta0,
                            double * lambda,
                            gsl_multifit_linear_workspace * work);

GSL_FUN double
gsl_multifit_linear_gcv_calc(const double lambda,
                             const gsl_vector * UTy,
                             const double delta0,
                             gsl_multifit_linear_workspace * work);

GSL_FUN int
gsl_multifit_linear_gcv(const gsl_vector * y,
                        gsl_vector * reg_param,
                        gsl_vector * G,
                        double * lambda,
                        double * G_lambda,
                        gsl_multifit_linear_workspace * work);

typedef struct
{
  const char * name;     /* method name */
  int (*wfun)(const gsl_vector *r, gsl_vector *w);
  int (*psi_deriv)(const gsl_vector *r, gsl_vector *dpsi);
  double tuning_default; /* default tuning constant */
} gsl_multifit_robust_type;

typedef struct
{
  double sigma_ols;    /* OLS estimate of sigma */
  double sigma_mad;    /* MAD estimate of sigma */
  double sigma_rob;    /* robust estimate of sigma */
  double sigma;        /* final estimate of sigma */
  double Rsq;          /* R^2 coefficient of determination */
  double adj_Rsq;      /* degree of freedom adjusted R^2 */
  double rmse;         /* root mean squared error */
  double sse;          /* residual sum of squares */
  size_t dof;          /* degrees of freedom */
  size_t numit;        /* number of iterations */
  gsl_vector *weights; /* final weights */
  gsl_vector *r;       /* final residuals y - X c */
} gsl_multifit_robust_stats;

typedef struct
{
  size_t n;            /* number of observations */
  size_t p;            /* number of parameters */
  size_t numit;        /* number of iterations */
  size_t maxiter;      /* maximum iterations */
  const gsl_multifit_robust_type *type;
  double tune;         /* tuning parameter */

  gsl_vector *r;       /* residuals at current iteration */
  gsl_vector *weights; /* weights at current iteration */
  gsl_vector *c_prev;  /* coefficients from previous iteration */
  gsl_vector *resfac;  /* multiplicative factors for residuals */

  gsl_vector *psi;     /* psi(r) */
  gsl_vector *dpsi;    /* psi'(r) */

  gsl_matrix *QSI;     /* Q S^{-1} of original matrix X */
  gsl_vector *D;       /* balancing parameters of original matrix X */

  gsl_vector *workn;   /* workspace of length n */

  gsl_multifit_robust_stats stats; /* various statistics */

  gsl_multifit_linear_workspace *multifit_p;
} gsl_multifit_robust_workspace;

/* available types */
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_default;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_bisquare;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_cauchy;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_fair;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_huber;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_ols;
GSL_VAR const gsl_multifit_robust_type * gsl_multifit_robust_welsch;

GSL_FUN gsl_multifit_robust_workspace *gsl_multifit_robust_alloc(const gsl_multifit_robust_type *T,
                                                         const size_t n, const size_t p);
GSL_FUN void gsl_multifit_robust_free(gsl_multifit_robust_workspace *w);
GSL_FUN int gsl_multifit_robust_tune(const double tune,
                             gsl_multifit_robust_workspace *w);
GSL_FUN int gsl_multifit_robust_maxiter(const size_t maxiter,
                                gsl_multifit_robust_workspace *w);
GSL_FUN const char *gsl_multifit_robust_name(const gsl_multifit_robust_workspace *w);
GSL_FUN gsl_multifit_robust_stats gsl_multifit_robust_statistics(const gsl_multifit_robust_workspace *w);
GSL_FUN int gsl_multifit_robust_weights(const gsl_vector *r, gsl_vector *wts,
                                gsl_multifit_robust_workspace *w);
GSL_FUN int gsl_multifit_robust(const gsl_matrix * X, const gsl_vector * y,
                        gsl_vector * c, gsl_matrix *cov,
                        gsl_multifit_robust_workspace *w);
GSL_FUN int gsl_multifit_robust_est(const gsl_vector * x, const gsl_vector * c,
                            const gsl_matrix * cov, double *y, double *y_err);
GSL_FUN int gsl_multifit_robust_residuals(const gsl_matrix * X,
                                  const gsl_vector * y,
                                  const gsl_vector * c, gsl_vector * r,
                                  gsl_multifit_robust_workspace * w);

__END_DECLS

#endif /* __GSL_MULTIFIT_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA


text:
/**
 *
 * @file core_spamm.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Dulceneia Becker
 * @date 2011-06-14
 * @generated s Tue Jan  7 11:44:49 2014
 *
 **/
#include <cblas.h>
#include <lapacke.h>
#include "common.h"

static inline int CORE_spamm_a2(PLASMA_enum side, PLASMA_enum trans, PLASMA_enum uplo,
                                int M, int N, int K, int L,
                                int vi2, int vi3,
                                      float *A2, int LDA2,
                                const float *V,  int LDV,
                                      float *W,  int LDW);

static inline int CORE_spamm_w(PLASMA_enum side, PLASMA_enum trans, PLASMA_enum uplo,
                               int M, int N, int K, int L,
                               int vi2, int vi3,
                               const float *A1, int LDA1,
                                     float *A2, int LDA2,
                               const float *V,  int LDV,
                                     float *W,  int LDW);

/***************************************************************************//**
 *
 * @ingroup CORE_float
 *
 *  ZPAMM  performs one of the matrix-matrix operations
 *
 *                    LEFT                      RIGHT
 *     OP PlasmaW  :  W  = A1 + op(V) * A2  or  W  = A1 + A2 * op(V)
 *     OP PlasmaA2 :  A2 = A2 - op(V) * W   or  A2 = A2 - W * op(V)
 *
 *  where  op( V ) is one of
 *
 *     op( V ) = V   or   op( V ) = V**T   or   op( V ) = V**T,
 *
 *  A1, A2 and W are general matrices, and V is:
 *
 *        l = k: rectangle + triangle
 *        l < k: rectangle + trapezoid
 *        l = 0: rectangle
 *
 *  Size of V, both rowwise and columnwise, is:
 *
 *         ----------------------
 *          side   trans    size
 *         ----------------------
 *          left     N     M x K
 *                   T     K x M
 *          right    N     K x N
 *                   T     N x K
 *         ----------------------
 *
 *  LEFT (columnwise and rowwise):
 *
 *              |    K    |                 |         M         |
 *           _  __________   _              _______________        _
 *              |    |    |                 |             | \
 *     V:       |    |    |            V':  |_____________|___\    K
 *              |    |    | M-L             |                  |
 *           M  |    |    |                 |__________________|   _
 *              |____|    |  _
 *              \    |    |                 |    M - L    | L  |
 *                \  |    |  L
 *           _      \|____|  _
 *
 *
 *  RIGHT (columnwise and rowwise):
 *
 *          |         K         |                   |    N    |
 *          _______________        _             _  __________   _
 *          |             | \                       |    |    |
 *     V':  |_____________|___\    N        V:      |    |    |
 *          |                  |                    |    |    | K-L
 *          |__________________|   _             K  |    |    |
 *                                                  |____|    |  _
 *          |    K - L    | L  |                    \    |    |
 *                                                    \  |    |  L
 *                                               _      \|____|  _
 *
 *  Arguments
 *  ==========
 *
 * @param[in] op
 *
 *         OP specifies which operation to perform:
 *
 *         @arg PlasmaW  : W  = A1 + op(V) * A2  or  W  = A1 + A2 * op(V)
 *         @arg PlasmaA2 : A2 = A2 - op(V) * W   or  A2 = A2 - W * op(V)
 *
 * @param[in] side
 *
 *         SIDE specifies whether  op( V ) multiplies A2
 *         or W from the left or right as follows:
 *
 *         @arg PlasmaLeft  : multiply op( V ) from the left
 *                            OP PlasmaW  :  W  = A1 + op(V) * A2
 *                            OP PlasmaA2 :  A2 = A2 - op(V) * W
 *
 *         @arg PlasmaRight : multiply op( V ) from the right
 *                            OP PlasmaW  :  W  = A1 + A2 * op(V)
 *                            OP PlasmaA2 :  A2 = A2 - W * op(V)
 *
 * @param[in] storev
 *
 *         Indicates how the vectors which define the elementary
 *         reflectors are stored in V:
 *
 *         @arg PlasmaColumnwise
 *         @arg PlasmaRowwise
 *
 * @param[in] M
 *         The number of rows of the A1, A2 and W
 *         If SIDE is PlasmaLeft, the number of rows of op( V )
 *
 * @param[in] N
 *         The number of columns of the A1, A2 and W
 *         If SIDE is PlasmaRight, the number of columns of op( V )
 *
 * @param[in] K
 *         If SIDE is PlasmaLeft, the number of columns of op( V )
 *         If SIDE is PlasmaRight, the number of rows of op( V )
 *
 * @param[in] L
 *         The size of the triangular part of V
 *
 * @param[in] A1
 *         On entry, the M-by-N tile A1.
 *
 * @param[in] LDA1
 *         The leading dimension of the array A1. LDA1 >= max(1,M).
 *
 * @param[in,out] A2
 *         On entry, the M-by-N tile A2.
 *         On exit, if OP is PlasmaA2 A2 is overwritten
 *
 * @param[in] LDA2
 *         The leading dimension of the tile A2. LDA2 >= max(1,M).
 *
 * @param[in] V
 *         The matrix V as described above.
 *         If SIDE is PlasmaLeft : op( V ) is M-by-K
 *         If SIDE is PlasmaRight: op( V ) is K-by-N
 *
 * @param[in] LDV
 *         The leading dimension of the array V.
 *
 * @param[in,out] W
 *         On entry, the M-by-N matrix W.
 *         On exit, W is overwritten either if OP is PlasmaA2 or PlasmaW.
 *         If OP is PlasmaA2, W is an input and is used as a workspace.
 *
 * @param[in] LDW
 *         The leading dimension of array WORK.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 ******************************************************************************/
int
CORE_spamm(int op, PLASMA_enum side, PLASMA_enum storev,
           int M, int N, int K, int L,
           const float *A1, int LDA1,
                 float *A2, int LDA2,
           const float *V, int LDV,
                 float *W, int LDW)
{

    int vi2, vi3, uplo, trans, info;

    /* Check input arguments */
    if ((op != PlasmaW) && (op != PlasmaA2)) {
        coreblas_error(1, "Illegal value of op");
        return -1;
    }
    if ((side != PlasmaLeft) && (side != PlasmaRight)) {
        coreblas_error(2, "Illegal value of side");
        return -2;
    }
    if ((storev != PlasmaColumnwise) && (storev != PlasmaRowwise)) {
        coreblas_error(3, "Illegal value of storev");
        return -3;
    }
    if (M < 0) {
        coreblas_error(4, "Illegal value of M");
        return -4;
    }
    if (N < 0) {
        coreblas_error(5, "Illegal value of N");
        return -5;
    }
    if (K < 0) {
        coreblas_error(6, "Illegal value of K");
        return -6;
    }
    if (L < 0) {
        coreblas_error(7, "Illegal value of L");
        return -7;
    }
    if (LDA1 < 0) {
        coreblas_error(9, "Illegal value of LDA1");
        return -9;
    }
    if (LDA2 < 0) {
        coreblas_error(11, "Illegal value of LDA2");
        return -11;
    }
    if (LDV < 0) {
        coreblas_error(13, "Illegal value of LDV");
        return -13;
    }
    if (LDW < 0) {
        coreblas_error(15, "Illegal value of LDW");
        return -15;
    }

    /* Quick return */
    if ((M == 0) || (N == 0) || (K == 0))
        return PLASMA_SUCCESS;

    /*
     * TRANS is set as:
     *
     *        -------------------------------------
     *         side   direct     PlasmaW  PlasmaA2
     *        -------------------------------------
     *         left   colwise       T        N
     *                rowwise       N        T
     *         right  colwise       N        T
     *                rowwise       T        N
     *        -------------------------------------
     */

    /* Columnwise*/
    if (storev == PlasmaColumnwise) {
        uplo = CblasUpper;
        if (side == PlasmaLeft) {
            trans = op == PlasmaA2 ? PlasmaNoTrans : PlasmaTrans;
            vi2 = trans == PlasmaNoTrans ? M - L : K - L;
        }
        else {
            trans = op == PlasmaW ? PlasmaNoTrans : PlasmaTrans;
            vi2 = trans == PlasmaNoTrans ? K - L : N - L;
        }
        vi3 = LDV * L;
    }

    /* Rowwise */
    else {
        uplo = CblasLower;
        if (side == PlasmaLeft) {
            trans = op == PlasmaW ? PlasmaNoTrans : PlasmaTrans;
            vi2 = trans == PlasmaNoTrans ? K - L : M - L;
        }
        else {
            trans = op == PlasmaA2 ? PlasmaNoTrans : PlasmaTrans;
            vi2 = trans == PlasmaNoTrans ? N - L : K - L;
        }
        vi2 *= LDV;
        vi3  = L;
    }

    /**/
    if (op==PlasmaW) {
        info = CORE_spamm_w(
                side, trans, uplo, M, N, K, L, vi2, vi3,
                A1, LDA1, A2, LDA2, V, LDV, W, LDW);
        if (info != 0)
            return info;
    } else if (op==PlasmaA2) {
        info = CORE_spamm_a2(
                side, trans, uplo, M, N, K, L, vi2, vi3,
                A2, LDA2, V, LDV, W, LDW);
        if (info != 0)
            return info;
    }

    return PLASMA_SUCCESS;
}

/***************************************************************************/
static inline int
CORE_spamm_w(PLASMA_enum side, PLASMA_enum trans, PLASMA_enum uplo,
             int M, int N, int K, int L,
             int vi2, int vi3,
             const float *A1, int LDA1,
                   float *A2, int LDA2,
             const float *V, int LDV,
                   float *W, int LDW)
{

   /*
    * W = A1 + op(V) * A2  or  W = A1 + A2 * op(V)
    */

    int j;
    static float zone  =  1.0;
    static float zzero =  0.0;

    if (side == PlasmaLeft) {

        if (((trans == PlasmaTrans) && (uplo == CblasUpper)) ||
            ((trans == PlasmaNoTrans) && (uplo == CblasLower))) {

            /*
             * W = A1 + V' * A2
             */

            /* W = A2_2 */
            LAPACKE_slacpy_work(LAPACK_COL_MAJOR,
                lapack_const(PlasmaUpperLower),
                L, N,
                &A2[K-L], LDA2, W, LDW);

            /* W = V_2' * W + V_1' * A2_1 (ge+tr, top L rows of V') */
            if (L > 0) {
                /* W = V_2' * W */
                cblas_strmm(
                    CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, L, N,
                    (zone), &V[vi2], LDV,
                    W, LDW);

                /* W = W + V_1' * A2_1 */
                if (K > L) {
                    cblas_sgemm(
                        CblasColMajor, (CBLAS_TRANSPOSE)trans, CblasNoTrans,
                        L, N, K-L,
                        (zone), V, LDV,
                        A2, LDA2,
                        (zone), W, LDW);
                }
            }

            /* W_2 = V_3' * A2: (ge, bottom M-L rows of V') */
            if (M > L) {
                cblas_sgemm(
                    CblasColMajor, (CBLAS_TRANSPOSE)trans, CblasNoTrans,
                    (M-L), N, K,
                    (zone), &V[vi3], LDV,
                    A2, LDA2,
                    (zzero), &W[L], LDW);
            }

            /* W = A1 + W */
            for(j = 0; j < N; j++) {
                cblas_saxpy(
                        M, (zone),
                        &A1[LDA1*j], 1,
                        &W[LDW*j], 1);
            }
        }
        else {
            printf("Left Upper/NoTrans & Lower/ConjTrans not implemented yet\n");
            return PLASMA_ERR_NOT_SUPPORTED;

        }
    }
    else { //side right

        if (((trans == PlasmaTrans) && (uplo == CblasUpper)) ||
            ((trans == PlasmaNoTrans) && (uplo == CblasLower))) {
            printf("Right Upper/ConjTrans & Lower/NoTrans not implemented yet\n");
            return PLASMA_ERR_NOT_SUPPORTED;

        }
        else {

            /*
             * W = A1 + A2 * V
             */

            if (L > 0) {

                /* W = A2_2 */
                LAPACKE_slacpy_work(LAPACK_COL_MAJOR,
                    lapack_const(PlasmaUpperLower),
                    M, L,
                    &A2[LDA2*(K-L)], LDA2, W, LDW);

                /* W = W * V_2 --> W = A2_2 * V_2 */
                cblas_strmm(
                    CblasColMajor, CblasRight, (CBLAS_UPLO)uplo,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, M, L,
                    (zone), &V[vi2], LDV,
                    W, LDW);

                /* W = W + A2_1 * V_1 */
                if (K > L) {
                    cblas_sgemm(
                        CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)trans,
                        M, L, K-L,
                        (zone), A2, LDA2,
                        V, LDV,
                        (zone), W, LDW);
                }

            }

            /* W = W + A2 * V_3 */
            if (N > L) {
                cblas_sgemm(
                    CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)trans,
                    M, N-L, K,
                    (zone), A2, LDA2,
                    &V[vi3], LDV,
                    (zzero), &W[LDW*L], LDW);
            }

            /* W = A1 + W */
            for (j = 0; j < N; j++) {
                cblas_saxpy(
                        M, (zone),
                        &A1[LDA1*j], 1,
                        &W[LDW*j], 1);
            }
        }
    }

    return PLASMA_SUCCESS;
}

/***************************************************************************/
static inline int
CORE_spamm_a2(PLASMA_enum side, PLASMA_enum trans, PLASMA_enum uplo,
              int M, int N, int K, int L,
              int vi2, int vi3,
                    float *A2, int LDA2,
              const float *V, int LDV,
                    float *W, int LDW)
{

   /*
    *  A2 = A2 + op(V) * W  or  A2 = A2 + W * op(V)
    */

    int j;
    static float zone  =  1.0;
    static float mzone  =  -1.0;

    if (side == PlasmaLeft) {

        if (((trans == PlasmaTrans) && (uplo == CblasUpper)) ||
            ((trans == PlasmaNoTrans) && (uplo == CblasLower))) {

            printf("Left Upper/ConjTrans & Lower/NoTrans not implemented yet\n");
            return PLASMA_ERR_NOT_SUPPORTED;

        }
        else {  //trans

            /*
             * A2 = A2 - V * W
             */

            /* A2_1 = A2_1 - V_1  * W_1 */
            if (M > L) {
                cblas_sgemm(
                    CblasColMajor, (CBLAS_TRANSPOSE)trans, CblasNoTrans,
                    M-L, N, L,
                    (mzone), V, LDV,
                    W, LDW,
                    (zone), A2, LDA2);
            }

            /* W_1 = V_2 * W_1 */
            cblas_strmm(
                CblasColMajor, CblasLeft, (CBLAS_UPLO)uplo,
                (CBLAS_TRANSPOSE)trans, CblasNonUnit, L, N,
                (zone), &V[vi2], LDV,
                W, LDW);

            /* A2_2 = A2_2 - W_1 */
            for(j = 0; j < N; j++) {
                cblas_saxpy(
                    L, (mzone),
                    &W[LDW*j], 1,
                    &A2[LDA2*j+(M-L)], 1);
            }

            /* A2 = A2 - V_3  * W_2 */
            if (K > L) {
                cblas_sgemm(
                    CblasColMajor, (CBLAS_TRANSPOSE)trans, CblasNoTrans,
                    M, N, (K-L),
                    (mzone), &V[vi3], LDV,
                    &W[L], LDW,
                    (zone), A2, LDA2);
            }

        }
    }
    else { //side right

        if (((trans == PlasmaTrans) && (uplo == CblasUpper)) ||
            ((trans == PlasmaNoTrans) && (uplo == CblasLower))) {

            /*
             * A2 = A2 - W * V'
             */

            /* A2 = A2 - W_2 * V_3' */
            if (K > L) {
                cblas_sgemm(
                    CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)trans,
                    M, N, K-L,
                    (mzone), &W[LDW*L], LDW,
                    &V[vi3], LDV,
                    (zone), A2, LDA2);
            }

            /* A2_1 = A2_1 - W_1 * V_1' */
            if (N > L) {
                cblas_sgemm(
                    CblasColMajor, CblasNoTrans, (CBLAS_TRANSPOSE)trans,
                    M, N-L, L,
                    (mzone), W, LDW,
                    V, LDV,
                    (zone), A2, LDA2);
            }

            /* A2_2 =  A2_2 -  W_1 * V_2' */
            if (L > 0) {
                cblas_strmm(
                    CblasColMajor, CblasRight, (CBLAS_UPLO)uplo,
                    (CBLAS_TRANSPOSE)trans, CblasNonUnit, M, L,
                    (mzone), &V[vi2], LDV,
                    W, LDW);

                for (j = 0; j < L; j++) {
                    cblas_saxpy(
                        M, (zone),
                        &W[LDW*j], 1,
                        &A2[LDA2*(N-L+j)], 1);
                }
            }

        }
        else {
            printf("Right Upper/NoTrans & Lower/ConjTrans not implemented yet\n");
            return PLASMA_ERR_NOT_SUPPORTED;
        }
    }

    return PLASMA_SUCCESS;
}

ext: c


text:
/* expfit.c -- model functions for exponential + background */

#include <math.h>
#include <stddef.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

struct usertype {
  size_t n;
  double * y;
  double * sigma;
};

int expb_f (const gsl_vector * x, void *data,
        gsl_vector * f) {
  size_t n = ((struct usertype *)data)->n;
  double *y = ((struct usertype *)data)->y;
  double *sigma = ((struct usertype *) data)->sigma;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);
  double b = gsl_vector_get (x, 2);

  size_t i;

  for (i = 0; i < n; i++)
    {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp (-lambda * t) + b;
      gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
    }

  return GSL_SUCCESS;
}

int expb_df (const gsl_vector * x, void *params,
         gsl_matrix * J) {

  struct usertype *myparams = (struct usertype *) params;

  int n = myparams->n;
  double *sigma = myparams->sigma;

  double A = gsl_vector_get (x, 0);
  double lambda = gsl_vector_get (x, 1);

  size_t i;

  /*for (i = 0; i<10; i++){
    printf("sigma[%z] = %5.3f\n",i,sigma[i]);
    }*/

  for (i = 0; i < n; i++)
    {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double s = sigma[i];
      double e = exp(-lambda * t);
      gsl_matrix_set (J, i, 0, e/s);
      gsl_matrix_set (J, i, 1, -t * A * e/s);
      gsl_matrix_set (J, i, 2, 1/s);
    }
  return GSL_SUCCESS;
}

int expb_fdf (const gsl_vector * x, void *data,
          gsl_vector * f, gsl_matrix * J) {

  expb_f (x, data, f);
  expb_df (x, data, J);

  return GSL_SUCCESS;
}

void eval_F ( int fstatus, const int n, const int m,
	      const double *x, double *f, const void *params){

  // first, convert x from an array into a gsl_vector...
  gsl_vector * x_gsl = gsl_vector_alloc(n);
  int i;
  for (i = 0; i < n; i++) {
    gsl_vector_set(x_gsl,i,x[i]);
  }

  // then call the expb_f function
  gsl_vector * f_gsl =  gsl_vector_alloc(m);

  expb_f ( x_gsl,  params,  f_gsl);

  // then convert f back into an array...
  //  f = f_gsl->data;
  for (i = 0; i < m; i++){
    f[i] = gsl_vector_get( f_gsl,i);
    //    f[i] = f_gsl->data[i];
  }

  gsl_vector_free(x_gsl);
  gsl_vector_free(f_gsl);

}

void eval_J  ( int fstatus, const int n, const int m,
	       const double *x, double *J, const void *params){

  // first, convert x from an array into a gsl_vector...
  gsl_vector * x_gsl = gsl_vector_alloc(n);
  int i;
  int jj;
  for (i = 0; i < n; i++) {
    gsl_vector_set(x_gsl,i,x[i]);
  }

  // then call the expb_f function
  gsl_matrix * J_gsl = gsl_matrix_alloc(40,3);
  expb_df (x_gsl, params, J_gsl);

  // then convert J into an array...
  // (find better way...)
  gsl_matrix * J_gsl_t = gsl_matrix_alloc(3,40);
  gsl_matrix_transpose_memcpy(J_gsl_t, J_gsl);
  for ( i = 0; i < m*n; i++){
    J[i] = J_gsl_t->data[i];
  }
  // the following code loses the pointer to J (maybe...)
  // J = J_gsl_t->data;
  gsl_vector_free(x_gsl);
  gsl_matrix_free(J_gsl);
  gsl_matrix_free(J_gsl_t);

}

void eval_HF ( int fstatus, const int n, const int m,
	       const double *x, const double *f, double *hf, const void *params){

}

max_stars_repo_name: andpic/RALFit

max_issues_repo_name: andpic/RALFit

max_forks_repo_name: andpic/RALFit


text:
/*
 * Dice Tester - Counter tool to display statistics about the randomness of dice
 */

#include <gsl/gsl_cdf.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int
getSides(int argc, char **argv) {
  int sides = 6;
  if (argc >= 2 && atoi(argv[1]) > 0)
    sides = atoi(argv[1]);
  return sides;
}

double
chiSquared(double *observed, double *expected, int sides) {
  double sum = 0.0;
  for (int i = 0; i < sides; i++)
    sum += pow((observed[i] - expected[i]), 2) / expected[i];
  return sum;
}

int
main(int argc, char **argv)
{
  char input[8];
  int num, total, sides, memSize;
  double *observed, *expected;

  total = 0;

  sides = getSides(argc, argv);
  memSize = sizeof(double) * sides;
  observed = malloc(memSize);
  expected = malloc(memSize);
  printf("HERE %d\n", sides);
  for (int i = 0; i < sides; i++)
    observed[i] = 0.0;

  printf("HERE %d\n", sides);
  printf("Obs: %f\n", observed[0]);
  /* Load observed data */
  while (scanf("%s", input) != EOF) {
    num = atoi(input);

    if (num >= 1 && num <= sides) {
      printf("%d\n", num);
      observed[num-1] += 1;
      total += 1;
    }
    getchar();
  }

  /* Random = flat distribution */
  for (int i = 0; i < sides; i++)
    expected[i] = (double) total / sides;

  double chi2 = chiSquared(observed, expected, sides);
  double pval = 1-gsl_cdf_chisq_P(chi2, sides-1);

  fputs("\n########## RESULTS ##########\n", stderr);
  fputs("Num:\tExp:\tObs:\tPer:\t\n", stderr);
  for (int i = 0; i < sides; i++)
    fprintf(stderr, "%d:\t %2.0f\t %2.0f\t %3.0f%%\n", i+1, expected[i], observed[i], (observed[i] / total)*100);
  if (pval >= 0.1)
    fputs("\nDie is fair :)\n", stderr);
  else if (pval > 0.05)
    fputs("\nDie is probably fair\n", stderr);
  else
    fputs("\nDie not fair :(\n", stderr);

  free(observed);
  free(expected);

  return 0;
}


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_poly.h>
#include <gsl/gsl_interp.h>

typedef struct
{
  double *d;
  double *coeff;
  double *work;
}
polynomial_state_t;

static void *
polynomial_alloc (size_t size)
{
  polynomial_state_t *state =
    (polynomial_state_t *) malloc (sizeof (polynomial_state_t));

  if (state == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for polynomial state",
                      GSL_ENOMEM);
    }

  state->d = (double *) malloc (sizeof (double) * size);

  if (state->d == 0)
    {
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
    }

  state->coeff = (double *) malloc (sizeof (double) * size);

  if (state->coeff == 0)
    {
      free (state->d);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
    }

  state->work = (double *) malloc (sizeof (double) * size);

  if (state->work == 0)
    {
      free (state->coeff);
      free (state->d);
      free (state);
      GSL_ERROR_NULL ("failed to allocate space for d", GSL_ENOMEM);
    }

  return state;
}

static int
polynomial_init (void *vstate,
                 const double xa[], const double ya[], size_t size)
{
  polynomial_state_t *state = (polynomial_state_t *) vstate;

  int status = gsl_poly_dd_init (state->d, xa, ya, size);

  return status;
}

static int
polynomial_eval (const void *vstate,
                 const double xa[], const double ya[], size_t size, double x,
                 gsl_interp_accel * acc, double *y)
{
  const polynomial_state_t *state = (const polynomial_state_t *) vstate;

  *y = gsl_poly_dd_eval (state->d, xa, size, x);

  return GSL_SUCCESS;
}

static int
polynomial_deriv (const void *vstate,
                  const double xa[], const double ya[], size_t size, double x,
                  gsl_interp_accel * acc, double *y)
{
  const polynomial_state_t *state = (const polynomial_state_t *) vstate;

  gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);

  *y = state->coeff[1];

  return GSL_SUCCESS;
}

static int
polynomial_deriv2 (const void *vstate,
                   const double xa[], const double ya[], size_t size,
                   double x, gsl_interp_accel * acc, double *y)
{
  const polynomial_state_t *state = (const polynomial_state_t *) vstate;

  gsl_poly_dd_taylor (state->coeff, x, state->d, xa, size, state->work);

  *y = 2.0 * state->coeff[2];

  return GSL_SUCCESS;
}

static int
polynomial_integ (const void *vstate, const double xa[], const double ya[],
                  size_t size, gsl_interp_accel * acc, double a, double b,
                  double *result)
{
  const polynomial_state_t *state = (const polynomial_state_t *) vstate;
  size_t i;
  double sum;

  gsl_poly_dd_taylor (state->coeff, 0.0, state->d, xa, size, state->work);

  sum = state->coeff[0] * (b - a);

  for (i = 1; i < size; i++)
    {
      sum += state->coeff[i] * (pow (b, i + 1) - pow (a, i + 1)) / (i + 1.0);
    }

  *result = sum;

  return GSL_SUCCESS;
}

static void
polynomial_free (void *vstate)
{
  polynomial_state_t *state = (polynomial_state_t *) vstate;
  free (state->d);
  free (state->coeff);
  free (state->work);
  free (state);
}

static const gsl_interp_type polynomial_type = {
  "polynomial",
  3,
  &polynomial_alloc,
  &polynomial_init,
  &polynomial_eval,
  &polynomial_deriv,
  &polynomial_deriv2,
  &polynomial_integ,
  &polynomial_free,
};

const gsl_interp_type *gsl_interp_polynomial = &polynomial_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
// matrix.h

//  This matrix class is a C++ wrapper for the GNU Scientific Library
//  Copyright (C) 2001 Ramin Nakisa

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

#if !defined( _matrix_int_h )
#define _matrix_int_h

#ifdef __HP_aCC //for aCC B3910B A.01.27
#include <iostream.h>
#include <fstream.h>
#include <iomanip.h>
#else //for gcc3
#include <iostream>
#include <fstream>
#include <iomanip>
#endif

#include <math.h>
#include <stdlib.h>
#include <assert.h>
///
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix_int.h>
#include <gsl/gsl_linalg.h>
#include <gslwrap/permutation.h>
#include <gslwrap/vector_int.h>

#define type_is_int
#ifdef  type_is
#define type_is_double
#endif

namespace gsl
{

///
class matrix_int
{
#ifdef type_is_double
	friend class matrix_float;
	friend class matrix_int;
#endif
public:
	typedef int value_type;
	typedef vector_int vector_type;

	///
	matrix_int();
	///
	matrix_int( size_t new_rows, size_t new_cols, bool clear = true );

	template<class oclass>
	void copy(const oclass &other)
	{
		if ( static_cast<const void *>( this ) == static_cast<const void *>( &other ) )
			return;

		set_dimensions( other.get_rows(), other.get_cols() );
		for ( size_t i = 0; i < get_rows(); i++ )
		{
			for ( size_t j = 0; j < get_cols(); j++ )
			{
				gsl_matrix_int_set( m, i, j, (int)other(i,j));
			}
		}
	}
/*    	template<>  */
/*    	void copy<matrix_int>(const matrix_int &other)  */
/*    	{  */
/*    		set_dimensions(other.size1(),other.size2());  */
/*    		gsl_matrix_int_memcpy( m, other.m );  */
/*    	}  */
	// copy constructor for type matrix_int
	matrix_int( const matrix_int &other ):m(NULL) {copy(other);}
	///
	template<class oclass>
	matrix_int( const oclass &other ):m(NULL) {copy(other);}

	///
	~matrix_int();
	///
//	matrix_int( const char *Filename );
	///
	size_t get_rows() const {return m->size1;}
	///
	size_t get_cols() const {return m->size2;}
	///
	size_t size1() const {return m->size1;}
	///
	size_t size2() const {return m->size2;}

	///
	void dimensions( size_t *num_rows, size_t *num_cols ) const;
	///
int        get_element ( size_t row, size_t col ) const {return  gsl_matrix_int_get( m, row, col ) ;}
	const int &operator()( size_t row, size_t col ) const {return *gsl_matrix_int_ptr( m, row, col ) ;}
int       &operator()( size_t row, size_t col )       {return *gsl_matrix_int_ptr( m, row, col ) ;}
	///
	void set_element( size_t row, size_t col, const int &v ){ gsl_matrix_int_set( m, row, col, v );}
	///
	void set_elements( const int & new_value );
	void set_all ( const int & new_value ) {gsl_matrix_int_set_all ( m, new_value );}
	void set_zero() {gsl_matrix_int_set_zero( m );}
	///
	void set_dimensions( size_t new_rows, size_t new_cols );
	///
	void load( const char *filename );
	///
	void save( const char *filename ) const;
	///
	friend ostream& operator<< ( ostream& os, const matrix_int& m );
	//This function writes the elements of the matrix m to the stream stream in binary format. The return value is 0 for success and GSL_EFAILED if there was a problem writing to the file. Since the data is written in the native binary format it may not be portable between different architectures.
	int fwrite (FILE * stream) const {return gsl_matrix_int_fwrite (stream, m);}

//This function reads into the matrix m from the open stream stream in binary format. The matrix m must be preallocated with the correct dimensions since the function uses the size of m to determine how many bytes to read. The return value is 0 for success and GSL_EFAILED if there was a problem reading from the file. The data is assumed to have been written in the native binary format on the same architecture.
	int fread (FILE * stream) {return gsl_matrix_int_fread (stream, m);}

    ///
	void load_binary( const char *filename );
	///
	void save_binary( const char *filename ) const;
	///
	bool operator==( const matrix_int &other ) const;
	bool operator!=( const matrix_int &other ) const {return !((*this)==other);}

	matrix_int& operator=( const matrix_int &other ) {copy( other );return *this;}
	/// converts from any other matrix type
	template<class omatrix>
	matrix_int &operator=( const omatrix& other )
	{
			copy(other);
			return *this;
	}
   ///
	matrix_int operator+( const matrix_int &other ) const;
	///
	matrix_int operator+( const int &f ) const;
	///
	friend matrix_int operator+( const int &f, const matrix_int &other );
	///
	matrix_int &operator+=( const int &f );
	///
	matrix_int &operator+=( const matrix_int &other );
	///
	matrix_int operator-( const matrix_int &other ) const;
	///
	matrix_int operator-( const int &f ) const;
	///
	friend matrix_int operator-( const int &f, const matrix_int &other );
	///
	matrix_int &operator-=( const int &f );
	///
	matrix_int &operator-=( const matrix_int &other );
	///
	matrix_int operator*( const matrix_int &other ) const;
	///
	matrix_int operator*( const int &f ) const;
	///
	friend matrix_int operator*( const int &f, const matrix_int &other );
	///
	matrix_int &operator*=( const int &f );
	///
	matrix_int &operator*=( const matrix_int &other );
	///
	matrix_int operator/( const int &) const;
	///
	matrix_int &operator/=( const int &);
	///
	matrix_int transpose() const;
	///
	matrix_int LU_decomp(gsl::permutation *perm=NULL,int *psign=NULL) const;
	///
	matrix_int LU_invert() const;

	// return a submatrix of the this from row_min to row_max (not included!)
	matrix_int submatrix(size_t row_min, size_t row_max, size_t col_min, size_t col_max) const
		{
			matrix_int m(row_max - row_min, col_max - col_min);
			for (size_t i = row_min ; i < row_max ; i++)
			{
				for (size_t j = col_min ; j < col_max ; j++)
				{
					m(i - row_min,j - col_min) = (*this)(i,j);
				}
			}
			return m;
		}
private:
	///
	void LU_decomp( gsl_matrix_int **a,
					gsl_permutation **permutation,
					int *sign ) const;
public:
	/** returns sum of all the matrix elements. */
    int sum() const;
	/** returns logarithm of the determinant of the matrix. */
	double LU_lndet() const;

	/** returns a vector_int_view of a single row of the matrix. */
	vector_int_view       row( size_t rowindex );
	const vector_int_view row( size_t rowindex ) const ;
	/** returns a vector_int_view of a single column of the matrix. */
	vector_int_view       column( size_t colindex );
	const vector_int_view column( size_t colindex ) const;
	/** returns a vector_int_view of the diagonal elements of the matrix. */
	vector_int_view       diagonal();
	const vector_int_view diagonal() const;

	/** returns a column matrix containing a single row of the matrix. */
	matrix_int get_row( size_t rowindex ) const;
	/** returns a column matrix containing a single column of the matrix. */
	matrix_int get_col( size_t colindex ) const;
	/** calculates sum of rows returned as a column matrix. */
	matrix_int row_sum() const;
	/** calculates sum of columns returned as a row matrix. */
	matrix_int column_sum() const;
	/** returns trace (diagonal sum) of a square matrix. */
	double trace() const;
	/** calculates cholesky decomposition of the matrix, returning success if matrix is positive definite. */
	int cholesky_decomp( matrix_int &a ) const;
//  	/** returns index of nearest row in matrix to vector argument. */
//  	int nearest_row_index( const matrix_int &v ) const;
	/** calculates covariance of the matrix columns. */
	matrix_int covariance() const;
	/** returns 1 if matrix is square, 0 otherwise. */
	bool is_square() const;
	/** diag operator (sets the diagonal elements of the matrix to the elements of v */
	void diag(const vector_int& v);
	/** set diagonal elements of a square matrix to f. */
	void set_diagonal( int f );
	/** sets matrix to a k dimensional unit matrix. */
	void identity( size_t k );
	/** returns sum of nth power of all elements. */
	double norm( double n ) const;

/*  Function: double gsl_matrix_max (const gsl_matrix * m)  */
/*      This function returns the maximum value in the matrix m.  */
	double max() const {return gsl_matrix_int_max(m);}
/*  Function: double gsl_matrix_min (const gsl_matrix * m)  */
/*      This function returns the minimum value in the matrix m.  */
	double min()const{return gsl_matrix_int_min(m);}

	/** This function returns 1 if all the elements of the matrix m are zero, and 0 otherwise. */
	bool isnull() const { return gsl_matrix_int_isnull(m);}
/*  Function: void gsl_matrix_minmax (const gsl_matrix * m, double * min_out, double * max_out)  */
/*      This function returns the minimum and maximum values in the matrix m, storing them in min_out and max_out.  */

/*  Function: void gsl_matrix_max_index (const gsl_matrix * m, size_t * imax, size_t * jmax)  */
/*      This function returns the indices of the maximum value in the matrix m, storing them in imax and jmax. When there are several equal maximum elements then the first element found */
/*      is returned.  */

/*  Function: void gsl_matrix_min_index (const gsl_matrix * m, size_t * imax, size_t * jmax)  */
/*      This function returns the indices of the minimum value in the matrix m, storing them in imax and jmax. When there are several equal minimum elements then the first element found */
/*      is returned.  */

/*  Function: void gsl_matrix_minmax_index (const gsl_matrix * m, size_t * imin, size_t * imax)  */
/*      This function returns the indices of the minimum and maximum values in the matrix m, storing them in (imin,jmin) and (imax,jmax). When there are several equal minimum or */
/*      maximum elements then the first elements found are returned.  */

	/** for interfacing with gsl c */
/*  	gsl_matrix_int       *gslobj()       {if (!m){cout << "matrix_int::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;} */
/*  	const gsl_matrix_int *gslobj() const {if (!m){cout << "matrix_int::gslobj ERROR, data not initialized!! " << endl; exit(-1);}return m;} */
	gsl_matrix_int       *gslobj()       {assert(m);return m;}
	const gsl_matrix_int *gslobj() const {assert(m);return m;}
private:
	///
   gsl_matrix_int *m;

};
}
#undef type_is_int
#undef type_is_double

#endif // _matrix_int_h

ext: h

max_stars_repo_name: entn-at/GlottDNN

max_issues_repo_name: entn-at/GlottDNN

max_forks_repo_name: entn-at/GlottDNN


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

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

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

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

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

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

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t
 *
 *  PLASMA_zcgels - Solves overdetermined or underdetermined linear systems involving an M-by-N
 *  matrix A using the QR or the LQ factorization of A.  It is assumed that A has full rank.
 *  The following options are provided:
 *
 *  # trans = PlasmaNoTrans and M >= N: find the least squares solution of an overdetermined
 *    system, i.e., solve the least squares problem: minimize || B - A*X ||.
 *
 *  # trans = PlasmaNoTrans and M < N:  find the minimum norm solution of an underdetermined
 *    system A * X = B.
 *
 *  Several right hand side vectors B and solution vectors X can be handled in a single call;
 *  they are stored as the columns of the M-by-NRHS right hand side matrix B and the N-by-NRHS
 *  solution matrix X.
 *
 *  PLASMA_zcgels first attempts to factorize the matrix in COMPLEX and use this
 *  factorization within an iterative refinement procedure to produce a
 *  solution with COMPLEX*16 normwise backward error quality (see below).
 *  If the approach fails the method switches to a COMPLEX*16
 *  factorization and solve.
 *
 *  The iterative refinement is not going to be a winning strategy if
 *  the ratio COMPLEX performance over COMPLEX*16 performance is too
 *  small. A reasonable strategy should take the number of right-hand
 *  sides and the size of the matrix into account. This might be done
 *  with a call to ILAENV in the future. Up to now, we always try
 *  iterative refinement.
 *
 *  The iterative refinement process is stopped if ITER > ITERMAX or
 *  for all the RHS we have: RNRM < N*XNRM*ANRM*EPS*BWDMAX
 *  where:
 *
 *  - ITER is the number of the current iteration in the iterative refinement process
 *  - RNRM is the infinity-norm of the residual
 *  - XNRM is the infinity-norm of the solution
 *  - ANRM is the infinity-operator-norm of the matrix A
 *  - EPS is the machine epsilon returned by DLAMCH('Epsilon').
 *
 *  Actually, in its current state (PLASMA 2.1.0), the test is slightly relaxed.
 *
 *  The values ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.
 *
 *  We follow Bjorck's algorithm proposed in "Iterative Refinement of Linear
 *  Least Squares solutions I", BIT, 7:257-278, 1967.
 *
 *******************************************************************************
 *
 * @param[in] trans
 *          Intended usage:
 *          = PlasmaNoTrans:   the linear system involves A;
 *          = PlasmaConjTrans: the linear system involves A**H.
 *          Currently only PlasmaNoTrans is supported.
 *
 * @param[in] M
 *          The number of rows of the matrix A. M >= 0.
 *
 * @param[in] N
 *          The number of columns of the matrix A. N >= 0.
 *
 * @param[in] NRHS
 *          The number of right hand sides, i.e., the number of columns of the matrices B and X.
 *          NRHS >= 0.
 *
 * @param[in] A
 *          The M-by-N matrix A. This matrix is not modified.
 *
 * @param[in] LDA
 *          The leading dimension of the array A. LDA >= max(1,M).
 *
 * @param[in] B
 *          The M-by-NRHS matrix B of right hand side vectors, stored columnwise. Not modified.
 *
 * @param[in] LDB
 *          The leading dimension of the array B. LDB >= MAX(1,M,N).
 *
 * @param[out] X
 *          If return value = 0, the solution vectors, stored columnwise.
 *          if M >= N, rows 1 to N of X contain the least squares solution vectors; the residual
 *          sum of squares for the solution in each column is given by the sum of squares of the
 *          modulus of elements N+1 to M in that column;
 *          if M < N, rows 1 to N of X contain the minimum norm solution vectors;
 *
 * @param[in] LDX
 *          The leading dimension of the array X. LDX >= MAX(1,M,N).
 *
 * @param[out] ITER
 *          The number of the current iteration in the iterative refinement process
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *          \retval <0 if -i, the i-th argument had an illegal value
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcgels_Tile
 * @sa PLASMA_zcgels_Tile_Async
 * @sa PLASMA_dsgels
 * @sa PLASMA_zgels
 *
 ******************************************************************************/
int PLASMA_zcgels(PLASMA_enum trans, int M, int N, int NRHS,
                  PLASMA_Complex64_t *A, int LDA,
                  PLASMA_Complex64_t *B, int LDB,
                  PLASMA_Complex64_t *X, int LDX, int *ITER)
{
    int i, j;
    int NB, NBNB, MT, NT, NTRHS;
    int status;
    PLASMA_desc  descA;
    PLASMA_desc  descB;
    PLASMA_desc *descT;
    PLASMA_desc  descX;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;

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

    *ITER = 0;

    /* Check input arguments */
    if (trans != PlasmaNoTrans   &&
        trans != PlasmaConjTrans &&
        trans != PlasmaTrans )
    {
        plasma_error("PLASMA_zcgels", "illegal value of trans");
        return -1;
    }
    if (trans != PlasmaNoTrans) {
        plasma_error("PLASMA_zcgels", "only PlasmaNoTrans supported");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (M < 0) {
        plasma_error("PLASMA_zcgels", "illegal value of M");
        return -2;
    }
    if (N < 0) {
        plasma_error("PLASMA_zcgels", "illegal value of N");
        return -3;
    }
    if (NRHS < 0) {
        plasma_error("PLASMA_zcgels", "illegal value of NRHS");
        return -4;
    }
    if (LDA < max(1, M)) {
        plasma_error("PLASMA_zcgels", "illegal value of LDA");
        return -6;
    }
    if (LDB < max(1, max(M, N))) {
        plasma_error("PLASMA_zcgels", "illegal value of LDB");
        return -9;
    }
    if (LDX < max(1, max(M, N))) {
        plasma_error("PLASMA_zcgels", "illegal value of LDX");
        return -10;
    }
    /* Quick return */
    if (min(M, min(N, NRHS)) == 0) {
        for (i = 0; i < max(M, N); i++)
            for (j = 0; j < NRHS; j++)
                B[j*LDB+i] = 0.0;
        return PLASMA_SUCCESS;
    }

    /* Tune NB & IB depending on M, N & NRHS; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_ZCGELS, M, N, NRHS);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcgels", "plasma_tune() failed");
        return status;
    }

    /* Set MT, NT & NTRHS */
    NB    = PLASMA_NB;
    NBNB  = NB*NB;
    NT    = (N%NB==0) ? (N/NB) : (N/NB+1);
    MT    = (M%NB==0) ? (M/NB) : (M/NB+1);
    NTRHS = (NRHS%NB==0) ? (NRHS/NB) : (NRHS/NB+1);
    printf("M %d, N %d, NRHS %d, NB %d, MT %d, NT %d, NTRHS %d\n", M, N, NRHS, NB, MT, NT, NTRHS);

    plasma_sequence_create(plasma, &sequence);

    descA = plasma_desc_init(
                PlasmaComplexDouble,
                NB, NB, NBNB,
                M, N, 0, 0, M, N);

    if (M >= N) {
        descB = plasma_desc_init(
            PlasmaComplexDouble,
            NB, NB, NBNB,
            M, NRHS, 0, 0, M, NRHS);

        descX = plasma_desc_init(
            PlasmaComplexDouble,
            NB, NB, NBNB,
            M, NRHS, 0, 0, M, NRHS);

    }
    else {
        descB = plasma_desc_init(
            PlasmaComplexDouble,
            NB, NB, NBNB,
            N, NRHS, 0, 0, N, NRHS);

        descX = plasma_desc_init(
            PlasmaComplexDouble,
            NB, NB, NBNB,
            N, NRHS, 0, 0, N, NRHS);
    }

    /* DOUBLE PRECISION INITIALIZATION */
    /* Allocate memory for matrices in block layout */
    if (plasma_desc_mat_alloc(&descA) || plasma_desc_mat_alloc(&descB) || plasma_desc_mat_alloc(&descX)) {
        plasma_error("PLASMA_zcgels", "plasma_shared_alloc() failed");
        plasma_desc_mat_free(&descA);
        plasma_desc_mat_free(&descB);
        plasma_desc_mat_free(&descX);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }

    plasma_parallel_call_5(plasma_pzlapack_to_tile,
        PLASMA_Complex64_t*, A,
        int, LDA,
        PLASMA_desc, descA,
        PLASMA_sequence*, sequence,
        PLASMA_request*, &request);

    plasma_parallel_call_5(plasma_pzlapack_to_tile,
        PLASMA_Complex64_t*, B,
        int, LDB,
        PLASMA_desc, descB,
        PLASMA_sequence*, sequence,
        PLASMA_request*, &request);

    /* Allocate workspace */
    PLASMA_Alloc_Workspace_zgels_Tile(M, N, &descT);

    /* Call the native interface */
    status = PLASMA_zcgels_Tile_Async(PlasmaNoTrans, &descA, descT, &descB, &descX, ITER,
                                      sequence, &request);

    if (status == PLASMA_SUCCESS) {
        plasma_parallel_call_5(plasma_pztile_to_lapack,
            PLASMA_desc, descX,
            PLASMA_Complex64_t*, X,
            int, LDX,
            PLASMA_sequence*, sequence,
            PLASMA_request*, &request);
    }
    plasma_dynamic_sync();

    PLASMA_Dealloc_Handle_Tile(&descT);
    plasma_sequence_destroy(plasma, sequence);
    plasma_desc_mat_free(&descA);
    plasma_desc_mat_free(&descB);
    plasma_desc_mat_free(&descX);
    return status;
}

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile
 *
 *  PLASMA_zcgels_Tile - Solves overdetermined or underdetermined linear system of equations
 *  using the tile QR or the tile LQ factorization and mixed-precision iterative refinement.
 *  Tile equivalent of PLASMA_zcgesv().
 *  Operates on matrices stored by tiles.
 *  All matrices are passed through descriptors.
 *  All dimensions are taken from the descriptors.
 *
 *******************************************************************************
 *
 * @param[in] trans
 *          Intended usage:
 *          = PlasmaNoTrans:   the linear system involves A;
 *          = PlasmaConjTrans: the linear system involves A**H.
 *          Currently only PlasmaNoTrans is supported.
 *
 * @param[in,out] A
 *          - If the iterative refinement converged, A is not modified;
 *          - otherwise, it fell back to double precision solution, and
 *          on exit the M-by-N matrix A contains:
 *          if M >= N, A is overwritten by details of its QR factorization as returned by
 *                     PLASMA_zgeqrf;
 *          if M < N, A is overwritten by details of its LQ factorization as returned by
 *                      PLASMA_zgelqf.
 *
 * @param[out] T
 *          On exit:
 *          - if the iterative refinement converged, T is not modified;
 *          - otherwise, it fell back to double precision solution,
 *          and then T is an auxiliary factorization data.
 *
 * @param[in,out] B
 *          On entry, the M-by-NRHS matrix B of right hand side vectors, stored columnwise;
 *          On exit, if return value = 0, B is overwritten by the solution vectors, stored
 *          columnwise:
 *          if M >= N, rows 1 to N of B contain the least squares solution vectors; the residual
 *          sum of squares for the solution in each column is given by the sum of squares of the
 *          modulus of elements N+1 to M in that column;
 *          if M < N, rows 1 to N of B contain the minimum norm solution vectors;
 *
 * @param[in] B
 *          The descriptor of the M-by-NRHS matrix B of right hand side vectors, stored columnwise. Not modified.
 *
 * @param[in,out] X
 *          On entry, it's only the descriptor where to store the result.
 *          On exit, if return value = 0, X is the solution vectors, stored columnwise:
 *          if M >= N, rows 1 to N of X contain the least squares solution vectors; the residual
 *          sum of squares for the solution in each column is given by the sum of squares of the
 *          modulus of elements N+1 to M in that column;
 *          if M < N, rows 1 to N of X contain the minimum norm solution vectors;
 *
 * @param[out] ITER
 *          If > 0, ITER is the number of the current iteration in the iterative refinement process.
 *          -ITERMAX-1, if the refinment step failed and the double precision factorization has been used.
 *
 *******************************************************************************
 *
 * @return
 *          \retval PLASMA_SUCCESS successful exit
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcgels
 * @sa PLASMA_zcgels_Tile_Async
 * @sa PLASMA_dsgels_Tile
 * @sa PLASMA_zgels_Tile
 *
 ******************************************************************************/
int PLASMA_zcgels_Tile(PLASMA_enum trans, PLASMA_desc *A, PLASMA_desc *T,
                       PLASMA_desc *B, PLASMA_desc *X, int *ITER)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

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

/***************************************************************************//**
 *
 * @ingroup PLASMA_Complex64_t_Tile_Async
 *
 *  PLASMA_zcgels_Tile_Async - Solves overdetermined or underdetermined linear
 *  system of equations using the tile QR or the tile LQ factorization and
 *  mixed-precision iterative refinement.
 *  Non-blocking equivalent of PLASMA_zcgels_Tile().
 *  May return before the computation is finished.
 *  Allows for pipelining of operations at runtime.
 *
 *******************************************************************************
 *
 * @param[in] sequence
 *          Identifies the sequence of function calls that this call belongs to
 *          (for completion checks and exception handling purposes).
 *
 * @param[out] request
 *          Identifies this function call (for exception handling purposes).
 *
 *******************************************************************************
 *
 * @sa PLASMA_zcgels
 * @sa PLASMA_zcgels_Tile
 * @sa PLASMA_dsgels_Tile_Async
 * @sa PLASMA_zgels_Tile_Async
 *
 ******************************************************************************/
int PLASMA_zcgels_Tile_Async(PLASMA_enum trans, PLASMA_desc *A, PLASMA_desc *T,
                             PLASMA_desc *B, PLASMA_desc *X, int *ITER,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    int M, N, NRHS, NB, NBNB, MT, NT, NTRHS;
    PLASMA_desc descA;
    PLASMA_desc descT;
    PLASMA_desc descB;
    PLASMA_desc descX;
    plasma_context_t *plasma;
    double *work;

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

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

    /* Check descriptors for correctness */
    if (plasma_desc_check(A) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcgels_Tile", "invalid first descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descA = *A;
    }
    if (plasma_desc_check(T) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcgels_Tile", "invalid second descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descT = *T;
    }
    if (plasma_desc_check(B) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcgels_Tile", "invalid third descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descB = *B;
    }
    if (plasma_desc_check(X) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zcgels_Tile", "invalid fourth descriptor");
        return PLASMA_ERR_ILLEGAL_VALUE;
    } else {
        descX = *X;
    }
    /* Check input arguments */
    if (descA.nb != descA.mb || descB.nb != descB.mb || descX.nb != descX.mb) {
        plasma_error("PLASMA_zcgels_Tile", "only square tiles supported");
        return PLASMA_ERR_ILLEGAL_VALUE;
    }
    if (trans != PlasmaNoTrans) {
        plasma_error("PLASMA_zcgels_Tile", "only PlasmaNoTrans supported");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    /* Quick return  - currently NOT equivalent to LAPACK's:
    if (min(M, min(N, NRHS)) == 0) {
        for (i = 0; i < max(M, N); i++)
            for (j = 0; j < NRHS; j++)
                B[j*LDB+i] = 0.0;
        return PLASMA_SUCCESS;
    }
    */

    if (0 == 0) {
    // START SPECIFIC

    /* Set M, M, NRHS, NB, MT, NT & NTRHS */
    M    = descA.lm;
    N    = descA.ln;
    NRHS = descB.ln;
    NB   = descA.nb;
    NBNB = NB*NB;

    MT = (M%NB==0) ? (M/NB) : (M/NB+1);
    NT = (N%NB==0) ? (N/NB) : (N/NB+1);
    NTRHS = (NRHS%NB==0) ? (NRHS/NB) : (NRHS/NB+1);
    printf("M %d, N %d, NRHS %d, NB %d, MT %d, NT %d, NTRHS %d\n", M, N, NRHS, NB, MT, NT, NTRHS);

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

    PLASMA_desc descR = plasma_desc_init(
        PlasmaComplexDouble,
        NB, NB, NBNB,
        M, NRHS, 0, 0, M, NRHS);

    if (plasma_desc_mat_alloc(&descR)) {
        plasma_error("PLASMA_zcgesv", "plasma_shared_alloc() failed");
        plasma_desc_mat_free(&descR);
        plasma_shared_free(plasma, work);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }

    PLASMA_desc descSA = plasma_desc_init(
            PlasmaComplexFloat,
            NB, NB, NBNB,
            M, N, 0, 0, M, N);

    PLASMA_desc descST = plasma_desc_init(
            PlasmaComplexFloat,
            IB, NB, IBNB,
            M, N, 0, 0, M, N);

    PLASMA_desc descSX = plasma_desc_init(
            PlasmaComplexFloat,
            NB, NB, NBNB,
            M, NRHS, 0, 0, M, NRHS);

    /* Allocate memory for single precision matrices in block layout */
    if (plasma_desc_mat_alloc(&descSA) || plasma_desc_mat_alloc(&descST) || plasma_desc_mat_alloc(&descSX)) {
        plasma_error("PLASMA_zcgesv", "plasma_shared_alloc() failed");
        plasma_desc_mat_free(&descSA);
        plasma_desc_mat_free(&descST);
        plasma_desc_mat_free(&descSX);
        plasma_desc_mat_free(&descR);
        plasma_shared_free(plasma, work);
        return PLASMA_ERR_OUT_OF_RESOURCES;
    }

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

    printf("Anorm=%e, cte=%e\n", Anorm, cte);

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

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

    if (descSA.m >= descSA.n) {

        /* Compute the QR factorization of SA */
        printf("Facto\n"); fflush(stdout);
        plasma_parallel_call_4(plasma_pcgeqrf,
            PLASMA_desc, descSA,
            PLASMA_desc, descST,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        printf("Solve\n"); fflush(stdout);
        plasma_parallel_call_5(plasma_pcunmqr,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_desc, descST,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_9(plasma_pctrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaUpper,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex32_t, 1.0,
            PLASMA_desc, plasma_desc_submatrix(descSA, 0, 0, descSA.n, descSA.n),
            PLASMA_desc, plasma_desc_submatrix(descSX, 0, 0, descSA.n, descSX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    }
    else {
        plasma_parallel_call_3(plasma_pztile_zero,
            PLASMA_desc, plasma_desc_submatrix(descSX, descSA.m, 0, descSA.n-descSA.m, descSX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_4(plasma_pcgelqf,
            PLASMA_desc, descSA,
            PLASMA_desc, descST,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_9(plasma_pctrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaLower,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex32_t, 1.0,
            PLASMA_desc, plasma_desc_submatrix(descSA, 0, 0, descSA.m, descSA.m),
            PLASMA_desc, plasma_desc_submatrix(descSX, 0, 0, descSA.m, descSX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_5(plasma_pcunmlq,
            PLASMA_desc, descSA,
            PLASMA_desc, descSX,
            PLASMA_desc, descST,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    }

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

    /* Compute R = B - AX. */
    printf("R = B - Ax\n"); fflush(stdout);
    printf("R = B - Ax ... cpy\n"); fflush(stdout);
    PLASMA_zlacpy(descB,descR);
    printf("R = B - Ax ... gemm\n"); fflush(stdout);

    plasma_parallel_call_9(plasma_pzgemm,
        PLASMA_enum, PlasmaNoTrans,
        PLASMA_enum, PlasmaNoTrans,
        PLASMA_Complex64_t, negone,
        PLASMA_desc, descA,
        PLASMA_desc, descX,
        PLASMA_Complex64_t, one,
        PLASMA_desc, descR,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* Check whether the NRHS normwise backward error satisfies the
       stopping criterion. If yes return. Note that ITER=0 (already set). */
    printf("Norm of X and R\n"); fflush(stdout);
    PLASMA_zlange(PlasmaInfNorm, descX, Xnorm, work);
    PLASMA_zlange(PlasmaInfNorm, descR, Rnorm, work);

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

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

    printf("Rnorm=%e, Xnorm * cte=%e, Rnorm=%e, cte=%e\n", Rnorm, Xnorm * cte, Rnorm, cte);

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

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

        /* Solve the system SA*SX = SR */
        if (descSA.m >= descSA.n) {

            plasma_parallel_call_5(plasma_pcunmqr,
                PLASMA_desc, descSA,
                PLASMA_desc, descSX,
                PLASMA_desc, descST,
                PLASMA_sequence*, sequence,
                PLASMA_request*, request);

            plasma_parallel_call_9(plasma_pctrsm,
                PLASMA_enum, PlasmaLeft,
                PLASMA_enum, PlasmaUpper,
                PLASMA_enum, PlasmaNoTrans,
                PLASMA_enum, PlasmaNonUnit,
                PLASMA_Complex32_t, 1.0,
                PLASMA_desc, plasma_desc_submatrix(descSA, 0, 0, descSA.n, descSA.n),
                PLASMA_desc, plasma_desc_submatrix(descSX, 0, 0, descSA.n, descSX.n),
                PLASMA_sequence*, sequence,
                PLASMA_request*, request);
        } else {
            plasma_parallel_call_9(plasma_pctrsm,
                PLASMA_enum, PlasmaLeft,
                PLASMA_enum, PlasmaLower,
                PLASMA_enum, PlasmaNoTrans,
                PLASMA_enum, PlasmaNonUnit,
                PLASMA_Complex32_t, 1.0,
                PLASMA_desc, plasma_desc_submatrix(descSA, 0, 0, descSA.m, descSA.m),
                PLASMA_desc, plasma_desc_submatrix(descSX, 0, 0, descSA.m, descSX.n),
                PLASMA_sequence*, sequence,
                PLASMA_request*, request);

            plasma_parallel_call_5(plasma_pcunmlq,
                PLASMA_desc, descSA,
                PLASMA_desc, descSX,
                PLASMA_desc, descST,
                PLASMA_sequence*, sequence,
                PLASMA_request*, request);
        }

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

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

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

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

        printf("Rnorm=%e, Xnorm * cte=%e, Rnorm=%e, cte=%e\n", Rnorm, Xnorm * cte, Rnorm, cte);

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

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

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

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

    printf("Go back DOUBLE\n");
    // END SPECIFIC
    }

    /* Single-precision iterative refinement failed to converge to a
       satisfactory solution, so we resort to double precision. */
    PLASMA_zlacpy(descB, descX);

    if (descA.m >= descA.n) {
        plasma_parallel_call_4(plasma_pzgeqrf,
            PLASMA_desc, descA,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_5(plasma_pzunmqr,
            PLASMA_desc, descA,
            PLASMA_desc, descX,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_9(plasma_pztrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaUpper,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex64_t, 1.0,
            PLASMA_desc, plasma_desc_submatrix(descA, 0, 0, descA.n, descA.n),
            PLASMA_desc, plasma_desc_submatrix(descX, 0, 0, descA.n, descX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    }
    else {
        plasma_parallel_call_3(plasma_pztile_zero,
            PLASMA_desc, plasma_desc_submatrix(descX, descA.m, 0, descA.n-descA.m, descX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_4(plasma_pzgelqf,
            PLASMA_desc, descA,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_9(plasma_pztrsm,
            PLASMA_enum, PlasmaLeft,
            PLASMA_enum, PlasmaLower,
            PLASMA_enum, PlasmaNoTrans,
            PLASMA_enum, PlasmaNonUnit,
            PLASMA_Complex64_t, 1.0,
            PLASMA_desc, plasma_desc_submatrix(descA, 0, 0, descA.m, descA.m),
            PLASMA_desc, plasma_desc_submatrix(descX, 0, 0, descA.m, descX.n),
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);

        plasma_parallel_call_5(plasma_pzunmlq,
            PLASMA_desc, descA,
            PLASMA_desc, descX,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    }
    return PLASMA_SUCCESS;
}


text:
/*! \file Function.h
 *  \brief header file that defines several function structures
 */
// Function.h

#ifndef NBODYFUNCTION_H
#define NBODYFUNCTION_H

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

namespace Math
{

/// \struct Math::gsl_math_function
/// \brief gsl function wrapper
struct gsl_math_function
{
  int (* function) (const gsl_vector *params, void * data, gsl_vector *output);
};

/// \struct Math::gsl_math_function
/// \brief gsl function wrapper differentiating wrt params, storing info in Jacobian matrix
struct gsl_math_function_df
{
  int (* function) (const gsl_vector * params, void * data, gsl_matrix * outputJacobian);
};

/// \struct Math::gsl_fitting_data
struct gsl_fitting_data {
  size_t n;
  double * x;
  double * y;
  size_t npar;
  int * iparindex;
  size_t nallparams;
  double *allparams;
};

/// \struct Math::math_function
/// \brief single dependent variable function
struct math_function
{
  Double_t (* function) (Double_t x, void * params);
  int (* gsl_function) (const gsl_vector *params, void * data, gsl_vector *output);
  int (* gsl_function_df) (const gsl_vector *params, void * data, gsl_matrix * outputJacobian);
  void * params;
};

/// \struct Math::math_multidim_function
/// \brief multi-dimensional function
struct math_multidim_function
{
  Double_t (* function) (Double_t *x, int ndim, void * params);
  int ndim;
  void * params;
};
}
#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#include <petsc.h>

typedef struct _msc_t {
  PetscInt nmasks;
  PetscInt* masks;
  PetscInt* mask_offsets;
  PetscInt* signs;
  PetscScalar* coeffs;
} msc_t;

typedef struct _shell_context {
  PetscInt nmasks;
  PetscInt* masks;
  PetscInt* mask_offsets;
  PetscInt* signs;
  PetscReal* real_coeffs;     // we store only the real or complex part -- whichever is nonzero
  void *left_subspace_data;
  void *right_subspace_data;
  PetscReal nrm;
} shell_context;

max_stars_repo_name: GregDMeyer/dynamite

max_issues_repo_name: GregDMeyer/dynamite

max_forks_repo_name: GregDMeyer/dynamite


text:
#include <stdio.h>
#include <pthread.h>
#include <time.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>

#if defined(PARSEC_HAVE_SCHED_SETAFFINITY) && defined (__USE_GNU) || 1
#include <linux/unistd.h>
#define gettid() syscall(__NR_gettid)
#else
#warning Not using sched_setaffinity
#endif  /* PARSEC_HAVE_SCHED_SETAFFINITY */

#include <sched.h>

#include <cblas.h>

#define MAX_THREADS 128

static int NB = 120;
static volatile int running = 1;

#define TYPE      float
#define GEMM_FUNC cblas_sgemm
/*
#define TYPE      double
#define GEMM_FUNC cblas_dgemm
*/
static void run_gemm(const TYPE *A, const TYPE *B, TYPE *C)
{
    GEMM_FUNC( (const enum CBLAS_ORDER)CblasColMajor,
               (const enum CBLAS_TRANSPOSE)CblasNoTrans, (const enum CBLAS_TRANSPOSE)CblasNoTrans,
               NB /* A.rows */,
               NB /* A.cols */,
               NB /* B.rows */,
               1.0, (TYPE*)A, NB /* A.data_stride */,
                    (TYPE*)B, NB /* B.data_stride */,
               1.0, (TYPE*)C, NB /* C.data_stride */);
}

static TYPE *init_matrix(void)
{
    TYPE *res;
    int i, j;
    res = (TYPE*)calloc(NB*NB, sizeof(TYPE));
    for(i = 0; i < NB; i++)
        for(j = 0; j < NB; j++)
            res[i*NB+j] = (TYPE)rand() / (TYPE)RAND_MAX;
    return res;
}

static void *thread_loop(void *_proc)
{
    TYPE *A, *B, *C;
    unsigned long int proc = (unsigned long int)_proc;
    int i;
    unsigned long long int time;
    struct timespec start, end;

    A = init_matrix();
    B = init_matrix();
    C = init_matrix();

#if defined(PARSEC_HAVE_SCHED_SETAFFINITY) && defined (__USE_GNU)
    {
        cpu_set_t cpuset;

        CPU_ZERO(&cpuset);
        CPU_SET(proc, &cpuset);

        if( -1 == sched_setaffinity(gettid(), sizeof(cpu_set_t), &cpuset) ) {
            printf( "Unable to set the thread affinity (%s)\n", strerror(errno) );
        }
    }
#endif  /* PARSEC_HAVE_SCHED_SETAFFINITY */

    if( proc == 0 ) {
        for(i = 0; i < 100; i++) {
            /* Ensures that all threads run cblas_dgemm on another core */
            /* while warming up */
            run_gemm(A, B, C);
        }

        for(i = 0; i < 1000; i++) {
            /* Take the time */
            clock_gettime(CLOCK_REALTIME, &start);
            run_gemm(A, B, C);
            clock_gettime(CLOCK_REALTIME, &end);
            time = end.tv_nsec - start.tv_nsec + ( (end.tv_sec - start.tv_sec) * 1000000000ULL );
            printf("NB = %d TIME = %llu ns  %f GFlops\n", NB, time,
                   (2*(NB/1e3)*(NB/1e3)*(NB/1e3)) / ((double)time  / 1e9));
        }

        running = 0;
    } else {
        while( running ) {
            run_gemm(A, B, C);
        }
    }
    return NULL;
}

int main(int argc, char *argv[])
{
    unsigned long int i, nbcores;
    pthread_t threads[MAX_THREADS];

    if( argc != 3 ) {
        fprintf(stderr, "usage: %s <Matrix Size> <Nb Cores>\n", argv[0]);
        return 1;
    }

    NB = atoi(argv[1]);
    nbcores = atoi(argv[2]);

    if( NB <= 1 || nbcores < 1 ) {
        fprintf(stderr, "usage: %s <Matrix Size> <Nb Cores>\n", argv[0]);
        return 1;
    }

    for(i = 0; i < nbcores; i++) {
        pthread_create(&threads[i], NULL, thread_loop, (void*)(i+1));
    }

    thread_loop((void*)0);

    for(i = 0; i < nbcores; i++) {
        pthread_join(threads[i], NULL);
    }
    return 0;
}

max_stars_repo_name: NLAFET/ABFT

max_issues_repo_name: NLAFET/ABFT

max_forks_repo_name: NLAFET/ABFT


text:
// -*- c++ -*-

#if 1
#include <string>
#include <vector>
#include <utility> // for std::pair
#include <stdint.h>

#include <iostream>
class Foo
{
public:
  static std::string GetName() { return "Foo"; }

  Foo();
  ~Foo();

  void setDouble(double d) { m_d = d;}
  void setDouble(double d, float x) { m_d = d+x;}
  //void setFoo(const std::string &s) {std::cout << s << std::endl;}
  void setFoo(const std::string &s="f", double x=1) {std::cout << s << x << std::endl;}
  void setInt(int i) { std::cout << i << std::endl;}
  void setInt(long i, float x=2) {std::cout << i+x << std::endl;}
  double getDouble() { return m_d; }

  void setUint64(uint64_t v) { m_uint64 = v; }
  uint64_t getUint64() { return m_uint64; }

  std::vector<Foo> getCollection() { return std::vector<Foo>(1, *this); }
  std::vector<std::vector<std::pair<int, int> > > getMatrix() { return std::vector<std::vector<std::pair<int,int> > >(); }

  std::vector< int>::iterator getIterator() { return std::vector<int>().begin(); }

  static void frobnicate() { std::cout << "frobnicate"; }
  //void frobnicate() { std::cout << "frobnicate -- non static"; }

  Foo operator+(const Foo& a) { return Foo(); }
  Foo& operator+=(const Foo& a) {
    this->m_d += a.m_d;
    this->m_uint64 += a.m_uint64;
    return *this;
  }

  Foo& operator=(const Foo& rhs) {
    this->m_d = rhs.m_d;
    this->m_uint64 = rhs.m_uint64;
    return *this;
  }

  void save(std::ostream &out) { out << m_d << m_uint64; }

  double getme(double d) { return d; }
  int    getme(int i) { return i; }
  void   getme() const { return ; }
  void   getme()       { return ; }

private:
  double m_d;
  uint64_t m_uint64;
};

Foo::Foo()
{}

Foo::~Foo()
{}

class IFoo
{
public:

  virtual
  void virtual_meth() const {}

  virtual
  void pure_virtual_meth() const = 0;

};

//template class std::vector<Foo>;
namespace Math {
  void do_hello() {
    std::cout << "helllo" << std::endl;
  }
  double do_add(double i, double j=2) { return i+j; }
  double do_add(int i, int j=2, int k=3, int l=4) { return i+j+k+l; }
  double do_add() { return 42; }

  const std::string& say_hello() { static std::string hi="hi"; return hi;}

  double adder(double i, double j) { return i+j; }
}

namespace Math2 {
  // std::string do_hello() {
  //   return "hi";
  // }
  std::string do_hello(const std::string& name = "you") {
    return std::string("hello ") + name;
  }
  std::string do_hello(const char* name) {
    return std::string("hello -- ") + std::string(name);
  }

  std::string do_hello_c(const char* name) {
    return std::string("hello -- ") + std::string(name);
  }
}

namespace NS {
  class Bar {
  public:
    Bar() {}
    void syHello() {
      std::cout << "hello" << std::endl;
    }
  };

  template<class F> F tmpl_fct() { return F(); }

  template<> int tmpl_fct<int>() { return 42; }

  template<int N> int tmpl_fct_n() { return N; }
  template<> int tmpl_fct_n<42>() { return 42; }
  template<> int tmpl_fct_n<-1>() { return -1; }
}

namespace TT {
  typedef int    foo_t;
  typedef foo_t* bar_t;
  typedef const bar_t baz_t;
  //'int *'    'foo *'    'bar'
}

#endif

#define BIT(n)       (1ULL << (n))
enum MyEnum {
  kValue = BIT(2)
};

typedef int            Ssiz_t;      //String size (int)
struct LongStr_t
{
  Ssiz_t    fCap;    // Max string length (including null)
  Ssiz_t    fSize;   // String length (excluding null)
  char     *fData;   // Long string data
};

enum Enum0 {
  kMinCap = (sizeof(LongStr_t) - 1)/sizeof(char) <= 2 ?
            2 : (sizeof(LongStr_t) - 1)/sizeof(char),
  kMin1 = 1 >= 2 ? 1 : 2
};

typedef void (*Func_t)();

class TVersionCheck {
public:
   TVersionCheck(int versionCode);  // implemented in TSystem.cxx
};

//static TVersionCheck gVersionCheck(335360);
//#define ROOT_VERSION_CODE 335360
//static TVersionCheck gVersionCheck = ROOT_VERSION_CODE;

class EnumNs
{
public:
   enum ESTLtype { kSTL       = 300 /* TVirtualStreamerInfo::kSTL */,
                   kSTLbitset =  8
   };
   // TStreamerElement status bits
   enum {
      kHasRange     = BIT(6),
      kCache        = BIT(9),
      kRepeat       = BIT(10),
      kRead         = BIT(11),
      kWrite        = BIT(12),
      kDoNotDelete  = BIT(13)
   };

};

class Base
{
public:
  virtual void initialize() = 0;
};

class Base2
{
public:
  virtual void execute() = 0;
};

class Alg : public Base, virtual public Base2
{
public:
  virtual void initialize() { std::cout << "::initialize();\n"; }
  virtual void execute() { std::cout << "::execute();\n"; }
};

class WithPrivateBase: public Base, private Base2
{
public:
  virtual void initialize() { std::cout << "::initialize();\n"; }
  virtual void execute() { std::cout << "::execute();\n"; }

  WithPrivateBase() {}
  WithPrivateBase(int /*i*/) {}

  Enum0 myenum;

  operator double() { return myenum; }

  enum Enum1 {
    kMinCap = (sizeof(LongStr_t) - 1)/sizeof(char) <= 2 ?
              2 : (sizeof(LongStr_t) - 1)/sizeof(char),
    kMin1 = 1 >= 2 ? 1 : 2
  };

private:
  void some_private_method() { std::cout << "::private()\n"; }
};

// #include <cblas.h>

max_stars_repo_name: sbinet/go-cxxdict

max_issues_repo_name: sbinet/go-cxxdict

max_forks_repo_name: sbinet/go-cxxdict