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

/*
 * Author:  G. Jungman
 */
#ifndef __GSL_BLAS_H__
#define __GSL_BLAS_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 <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

#include <gsl/gsl_blas_types.h>

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

__BEGIN_DECLS

/* ========================================================================
 * Level 1
 * ========================================================================
 */

GSL_FUN int gsl_blas_sdsdot (float alpha,
                     const gsl_vector_float * X,
                     const gsl_vector_float * Y,
                     float * result
                     );

GSL_FUN int gsl_blas_dsdot (const gsl_vector_float * X,
                    const gsl_vector_float * Y,
                    double * result
                    );

GSL_FUN int gsl_blas_sdot (const gsl_vector_float * X,
                   const gsl_vector_float * Y,
                   float * result
                   );

GSL_FUN int gsl_blas_ddot (const gsl_vector * X,
                   const gsl_vector * Y,
                   double * result
                   );

GSL_FUN int  gsl_blas_cdotu (const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_complex_float * dotu);

GSL_FUN int  gsl_blas_cdotc (const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_complex_float * dotc);

GSL_FUN int  gsl_blas_zdotu (const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_complex * dotu);

GSL_FUN int  gsl_blas_zdotc (const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_complex * dotc);

GSL_FUN float  gsl_blas_snrm2  (const gsl_vector_float * X);
GSL_FUN float  gsl_blas_sasum  (const gsl_vector_float * X);
GSL_FUN double gsl_blas_dnrm2  (const gsl_vector * X);
GSL_FUN double gsl_blas_dasum  (const gsl_vector * X);
GSL_FUN float  gsl_blas_scnrm2 (const gsl_vector_complex_float * X);
GSL_FUN float  gsl_blas_scasum (const gsl_vector_complex_float * X);
GSL_FUN double gsl_blas_dznrm2 (const gsl_vector_complex * X);
GSL_FUN double gsl_blas_dzasum (const gsl_vector_complex * X);

GSL_FUN CBLAS_INDEX_t gsl_blas_isamax (const gsl_vector_float * X);
GSL_FUN CBLAS_INDEX_t gsl_blas_idamax (const gsl_vector * X);
GSL_FUN CBLAS_INDEX_t gsl_blas_icamax (const gsl_vector_complex_float * X);
GSL_FUN CBLAS_INDEX_t gsl_blas_izamax (const gsl_vector_complex * X);

GSL_FUN int  gsl_blas_sswap (gsl_vector_float * X,
                     gsl_vector_float * Y);

GSL_FUN int  gsl_blas_scopy (const gsl_vector_float * X,
                     gsl_vector_float * Y);

GSL_FUN int  gsl_blas_saxpy (float alpha,
                     const gsl_vector_float * X,
                     gsl_vector_float * Y);

GSL_FUN int  gsl_blas_dswap (gsl_vector * X,
                     gsl_vector * Y);

GSL_FUN int  gsl_blas_dcopy (const gsl_vector * X,
                     gsl_vector * Y);

GSL_FUN int  gsl_blas_daxpy (double alpha,
                     const gsl_vector * X,
                     gsl_vector * Y);

GSL_FUN int  gsl_blas_cswap (gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

GSL_FUN int  gsl_blas_ccopy (const gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

GSL_FUN int  gsl_blas_caxpy (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     gsl_vector_complex_float * Y);

GSL_FUN int  gsl_blas_zswap (gsl_vector_complex * X,
                     gsl_vector_complex * Y);

GSL_FUN int  gsl_blas_zcopy (const gsl_vector_complex * X,
                     gsl_vector_complex * Y);

GSL_FUN int  gsl_blas_zaxpy (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     gsl_vector_complex * Y);

GSL_FUN int  gsl_blas_srotg (float a[], float b[], float c[], float s[]);

GSL_FUN int  gsl_blas_srotmg (float d1[], float d2[], float b1[], float b2, float P[]);

GSL_FUN int  gsl_blas_srot (gsl_vector_float * X,
                    gsl_vector_float * Y,
                    float c, float s);

GSL_FUN int  gsl_blas_srotm (gsl_vector_float * X,
                     gsl_vector_float * Y,
                     const float P[]);

GSL_FUN int  gsl_blas_drotg (double a[], double b[], double c[], double s[]);

GSL_FUN int  gsl_blas_drotmg (double d1[], double d2[], double b1[],
                      double b2, double P[]);

GSL_FUN int  gsl_blas_drot (gsl_vector * X,
                    gsl_vector * Y,
                    const double c, const double s);

GSL_FUN int  gsl_blas_drotm (gsl_vector * X,
                     gsl_vector * Y,
                     const double P[]);

GSL_FUN void gsl_blas_sscal  (float  alpha, gsl_vector_float * X);
GSL_FUN void gsl_blas_dscal  (double alpha, gsl_vector * X);
GSL_FUN void gsl_blas_cscal  (const gsl_complex_float alpha, gsl_vector_complex_float * X);
GSL_FUN void gsl_blas_zscal  (const gsl_complex alpha, gsl_vector_complex * X);
GSL_FUN void gsl_blas_csscal (float  alpha, gsl_vector_complex_float * X);
GSL_FUN void gsl_blas_zdscal (double alpha, gsl_vector_complex * X);

/* ===========================================================================
 * Level 2
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_FUN int  gsl_blas_sgemv (CBLAS_TRANSPOSE_t TransA,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_vector_float * X,
                     float beta,
                     gsl_vector_float * Y);

GSL_FUN int  gsl_blas_strmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_float * A,
                     gsl_vector_float * X);

GSL_FUN int  gsl_blas_strsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_float * A,
                     gsl_vector_float * X);

GSL_FUN int  gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_vector * X,
                     double beta,
                     gsl_vector * Y);

GSL_FUN int  gsl_blas_dtrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix * A,
                     gsl_vector * X);

GSL_FUN int  gsl_blas_dtrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix * A,
                     gsl_vector * X);

GSL_FUN int  gsl_blas_cgemv (CBLAS_TRANSPOSE_t TransA,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_vector_complex_float * X,
                     const gsl_complex_float beta,
                     gsl_vector_complex_float * Y);

GSL_FUN int  gsl_blas_ctrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex_float * A,
                     gsl_vector_complex_float * X);

GSL_FUN int  gsl_blas_ctrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex_float * A,
                     gsl_vector_complex_float * X);

GSL_FUN int  gsl_blas_zgemv (CBLAS_TRANSPOSE_t TransA,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_vector_complex * X,
                     const gsl_complex beta,
                     gsl_vector_complex * Y);

GSL_FUN int  gsl_blas_ztrmv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex * A,
                     gsl_vector_complex * X);

GSL_FUN int  gsl_blas_ztrsv (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t TransA, CBLAS_DIAG_t Diag,
                     const gsl_matrix_complex * A,
                     gsl_vector_complex *X);

/*
 * Routines with S and D prefixes only
 */
GSL_FUN int  gsl_blas_ssymv (CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_vector_float * X,
                     float beta,
                     gsl_vector_float * Y);

GSL_FUN int  gsl_blas_sger (float alpha,
                    const gsl_vector_float * X,
                    const gsl_vector_float * Y,
                    gsl_matrix_float * A);

GSL_FUN int  gsl_blas_ssyr (CBLAS_UPLO_t Uplo,
                    float alpha,
                    const gsl_vector_float * X,
                    gsl_matrix_float * A);

GSL_FUN int  gsl_blas_ssyr2 (CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_vector_float * X,
                     const gsl_vector_float * Y,
                     gsl_matrix_float * A);

GSL_FUN int  gsl_blas_dsymv (CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_vector * X,
                     double beta,
                     gsl_vector * Y);
GSL_FUN int  gsl_blas_dger (double alpha,
                    const gsl_vector * X,
                    const gsl_vector * Y,
                    gsl_matrix * A);

GSL_FUN int  gsl_blas_dsyr (CBLAS_UPLO_t Uplo,
                    double alpha,
                    const gsl_vector * X,
                    gsl_matrix * A);

GSL_FUN int  gsl_blas_dsyr2 (CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_vector * X,
                     const gsl_vector * Y,
                     gsl_matrix * A);

/*
 * Routines with C and Z prefixes only
 */

GSL_FUN int  gsl_blas_chemv (CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_vector_complex_float * X,
                     const gsl_complex_float beta,
                     gsl_vector_complex_float * Y);

GSL_FUN int  gsl_blas_cgeru (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

GSL_FUN int  gsl_blas_cgerc (const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

GSL_FUN int  gsl_blas_cher (CBLAS_UPLO_t Uplo,
                    float alpha,
                    const gsl_vector_complex_float * X,
                    gsl_matrix_complex_float * A);

GSL_FUN int  gsl_blas_cher2 (CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_vector_complex_float * X,
                     const gsl_vector_complex_float * Y,
                     gsl_matrix_complex_float * A);

GSL_FUN int  gsl_blas_zhemv (CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_vector_complex * X,
                     const gsl_complex beta,
                     gsl_vector_complex * Y);

GSL_FUN int  gsl_blas_zgeru (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

GSL_FUN int  gsl_blas_zgerc (const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

GSL_FUN int  gsl_blas_zher (CBLAS_UPLO_t Uplo,
                    double alpha,
                    const gsl_vector_complex * X,
                    gsl_matrix_complex * A);

GSL_FUN int  gsl_blas_zher2 (CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_vector_complex * X,
                     const gsl_vector_complex * Y,
                     gsl_matrix_complex * A);

/*
 * ===========================================================================
 * Prototypes for level 3 BLAS
 * ===========================================================================
 */

/*
 * Routines with standard 4 prefixes (S, D, C, Z)
 */
GSL_FUN int  gsl_blas_sgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_matrix_float * B,
                     float beta,
                     gsl_matrix_float * C);

GSL_FUN int  gsl_blas_ssymm (CBLAS_SIDE_t Side, CBLAS_UPLO_t Uplo,
                     float alpha,
                     const gsl_matrix_float * A,
                     const gsl_matrix_float * B,
                     float beta,
                     gsl_matrix_float * C);

GSL_FUN int  gsl_blas_ssyrk (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                     float alpha,
                     const gsl_matrix_float * A,
                     float beta,
                     gsl_matrix_float * C);

GSL_FUN int  gsl_blas_ssyr2k (CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t Trans,
                      float alpha,
                      const gsl_matrix_float * A,
                      const gsl_matrix_float * B,
                      float beta,
                      gsl_matrix_float * C);

GSL_FUN int  gsl_blas_strmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     float alpha,
                     const gsl_matrix_float * A,
                     gsl_matrix_float * B);

GSL_FUN int  gsl_blas_strsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     float alpha,
                     const gsl_matrix_float * A,
                     gsl_matrix_float * B);

GSL_FUN int  gsl_blas_dgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_matrix * B,
                     double beta,
                     gsl_matrix * C);

GSL_FUN int  gsl_blas_dsymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     double alpha,
                     const gsl_matrix * A,
                     const gsl_matrix * B,
                     double beta,
                     gsl_matrix * C);

GSL_FUN int  gsl_blas_dsyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     double alpha,
                     const gsl_matrix * A,
                     double beta,
                     gsl_matrix * C);

GSL_FUN int  gsl_blas_dsyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      double alpha,
                      const  gsl_matrix * A,
                      const  gsl_matrix * B,
                      double beta,
                      gsl_matrix * C);

GSL_FUN int  gsl_blas_dtrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     double alpha,
                     const gsl_matrix * A,
                     gsl_matrix * B);

GSL_FUN int  gsl_blas_dtrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     double alpha,
                     const gsl_matrix * A,
                     gsl_matrix * B);

GSL_FUN int  gsl_blas_cgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_csymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_csyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_csyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex_float alpha,
                      const gsl_matrix_complex_float * A,
                      const gsl_matrix_complex_float * B,
                      const gsl_complex_float beta,
                      gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_ctrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     gsl_matrix_complex_float * B);

GSL_FUN int  gsl_blas_ctrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     gsl_matrix_complex_float * B);

GSL_FUN int  gsl_blas_zgemm (CBLAS_TRANSPOSE_t TransA,
                     CBLAS_TRANSPOSE_t TransB,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

GSL_FUN int  gsl_blas_zsymm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

GSL_FUN int  gsl_blas_zsyrk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

GSL_FUN int  gsl_blas_zsyr2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex alpha,
                      const gsl_matrix_complex * A,
                      const gsl_matrix_complex * B,
                      const gsl_complex beta,
                      gsl_matrix_complex *C);

GSL_FUN int  gsl_blas_ztrmm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     gsl_matrix_complex * B);

GSL_FUN int  gsl_blas_ztrsm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo, CBLAS_TRANSPOSE_t TransA,
                     CBLAS_DIAG_t Diag,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     gsl_matrix_complex * B);

/*
 * Routines with prefixes C and Z only
 */
GSL_FUN int  gsl_blas_chemm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex_float alpha,
                     const gsl_matrix_complex_float * A,
                     const gsl_matrix_complex_float * B,
                     const gsl_complex_float beta,
                     gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_cherk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     float alpha,
                     const gsl_matrix_complex_float * A,
                     float beta,
                     gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_cher2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex_float alpha,
                      const gsl_matrix_complex_float * A,
                      const gsl_matrix_complex_float * B,
                      float beta,
                      gsl_matrix_complex_float * C);

GSL_FUN int  gsl_blas_zhemm (CBLAS_SIDE_t Side,
                     CBLAS_UPLO_t Uplo,
                     const gsl_complex alpha,
                     const gsl_matrix_complex * A,
                     const gsl_matrix_complex * B,
                     const gsl_complex beta,
                     gsl_matrix_complex * C);

GSL_FUN int  gsl_blas_zherk (CBLAS_UPLO_t Uplo,
                     CBLAS_TRANSPOSE_t Trans,
                     double alpha,
                     const gsl_matrix_complex * A,
                     double beta,
                     gsl_matrix_complex * C);

GSL_FUN int  gsl_blas_zher2k (CBLAS_UPLO_t Uplo,
                      CBLAS_TRANSPOSE_t Trans,
                      const gsl_complex alpha,
                      const gsl_matrix_complex * A,
                      const gsl_matrix_complex * B,
                      double beta,
                      gsl_matrix_complex * C);

__END_DECLS

#endif /* __GSL_BLAS_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#pragma once
#include <vector>
#include <array>
#include <iostream>
#include <stdexcept>
#include <cassert>

#include <gsl/gsl.h>

#include"utils/math.h"

namespace util{
namespace math{

template<typename T, int64_t M>
struct Vector{
    Vector(gsl::span<T, M> const vals) : _val(M, T{}) {
        // std::cerr<<M<<" Vector{gsl::span<M>}\n";
        std::copy_n(vals.cbegin(), M, _val.begin());
        assert(span[0]==_val[0]);
    }
    Vector() : _val(M, T{}) {
        // std::cerr<<"Default construct Vector{gsl::span<M>}\n";
        assert(span[0]==_val[0]);
    }
    //CAUTION: careful to correctly implement non-default constructor if necessary.
    Vector(Vector const & vec) : _val{vec._val} {
        if(!(span[0]==_val[0])) {
            std::cerr<<_val.size()<<" Copy construct Vector{gsl::span<M>}\n";
            std::cerr<<"Assert fails : " << span[0] << " vs " << _val[0] <<std::endl;
        }
        assert(span[0]==_val[0]);
        // std::cerr<<_val.size()<<" Copy construct Vector{gsl::span<M>}\n";
    };
    Vector& operator=(const Vector& vec) {
        // std::cerr<<"Copy assignment Vector{gsl::span<M>}\n";
        //  std::cerr<<"size: " << vec._val.size() << " vs " <<_val.size() <<std::endl;
        // assert(_val.size()==vec._val.size());
        // assert(span[0]==_val[0]);
        std::copy_n(vec.span.cbegin(), M, _val.begin());
        // span=gsl::span<T, M>{_val};
        assert(span[0]==_val[0]);
        // std::cerr<<"size: " << _val.size() <<std::endl;
        return *this;
    };
    Vector& operator+=(const Vector& vec) {
        auto n = span.size();
        for(decltype(n)i=0; i!=n; ++i) span[i] += vec.span[i];
        return *this;
    }
    Vector& operator-=(const Vector& vec) {
        span-=vec.span;
        return *this;
    }
    Vector& operator*=(T x) {
        for(auto &v : _val) v *= x;
        return *this;
    }
    Vector(Vector&& vec) : _val{std::move(vec._val)} {
        span=gsl::span<T, M>{_val};
        //Fail if nan
        assert(span[0]==_val[0]);
        // std::cerr<<"Move Vector{gsl::span<M>}\n";
    };

    auto l1_norm() const {
        T norm{0.0};
        for(auto x : _val) norm += std::abs(x);
        return norm;
    }

    std::vector<T> _val;
    gsl::span<T, M> span=_val; //this default is critical!
};
template<typename T, int64_t M>
struct VectorView{
    VectorView(gsl::span<T, M> const vals) :span{vals} {}
    gsl::span<T, M> span;
};

template<typename T, int64_t M, int64_t N>
struct Matrix{
    Matrix(gsl::span<T, M,N> const vals) : _val(M*N, T{}) {
        std::cerr<<M<<" "<<N<<" Matrix{gsl::span<M,N>}\n";
        std::copy_n(vals.begin(), M*N, _val.begin());
    }
    Matrix() :_val(M*N, T{}) {}
    Matrix(Matrix const & mat) : _val{mat._val}{};
    Matrix& operator=(const Matrix& mat) {
        std::copy_n(mat.span.cbegin(), M*N, _val.begin());
        return *this;
    }
    Matrix(Matrix&& vec) : _val{std::move(vec._val)} {
        span=gsl::span<T, M,N>{_val};
        //Fail if nan
        assert(span[0][0]==_val[0]);
    };
    std::vector<T> _val;
    gsl::span<T, M, N> span=_val; //this default is critical!
};
template<typename T, int64_t M, int64_t N>
struct MatrixView{
    MatrixView(gsl::span<T, M,N> const vals) : span{vals} {}
//    MatrixView(gsl::span<T> const vals) : span{vals.as_array_view<T,M,N>()} {}
    gsl::span<T, M, N> span;
};

template<typename T, int64_t L,int64_t M,int64_t N>
struct ThirdTensorView{
    ThirdTensorView(gsl::span<T, L,M,N> const vals) : span{vals} {}
//    MatrixView(gsl::span<T> const vals) : span{vals.as_array_view<T,M,N>()} {}
    gsl::span<T, L, M, N> span;
};

namespace factory{

//TODO: val can be const &
template<typename T, std::size_t M>
auto Vector(std::array<T,M> &val){
    const auto span = gsl::as_span(val);
    return util::math::Vector<T,int64_t{M}>{span};
}

// template<typename T>
// auto Vector(std::vector<T> const &val){
//     constexpr auto span = gsl::span<const T>{val};
//     return util::math::Vector<T,span.size()>{span};
// }

// template<typename T>
// auto Vector(std::vector<T> const &val){
//     return gsl::span<const T>{val};
// }
//
// template<typename T>
// auto Vector(T vals){
//     auto a = util::math::Vector{to_cspan(vals)};
//     return a;
// }
//

}//namespace util::math::factory

//TODO: use bi-dimensional index, gsl::index<2>??
template<typename T, int64_t M, int64_t N>
auto transpose(gsl::span<T,M,N> mat){
    std::vector<T> tr_mat(M*N);
    auto tr = gsl::span<T,N,M>{tr_mat};
    for (int64_t i=0; i<mat.extent(0); ++i) {
       for (int64_t j=0; j<mat.extent(1); ++j){
          tr[j][i] = mat[i][j];
      }
    }
    return Matrix<T,N,M>{tr};
}

template<typename T, int64_t M>
T dot(gsl::span<T,M> const x, gsl::span<T,M> const y){
    // return std::inner_product(x.cbegin(), x.cend(), y.cbegin(), T{});
    T sum{};
    #pragma clang loop vectorize(enable)
    for(decltype(M) i=0; i!=M; ++i){
        sum+=x[i]*y[i];
    }
    return sum;
}
//Following also works, without template, but less safe:
// auto dot = [](auto const &x, auto const &y){
//     typedef typename std::remove_reference<decltype(x[0])>::type A;
//     return std::inner_product(x.cbegin(), x.cend(), y.cbegin(), A{0});
// };

template<typename T, int64_t M, int64_t N>
auto dotdot(Vector<T,M> &x, Matrix<T,M,N> &m, Vector<T,N> &y){
    T sum{};
    for(decltype(M) i=0; i!=M; ++i){
        for(decltype(N) j=0; j!=N; ++j){
            sum += x.span[i]*m.span[i][j]*y.span[j];
        }
    }
    return sum;
}

auto mul_sum_vec=[](int64_t i, auto &out,
                auto const &a, auto const &b) {
    out+=a[i]*b[i];
};
auto mul_sum_mat=[](int64_t i,int64_t j, auto &out,
                auto const &a, auto const &b) {
    out+=a[i][j]*b[i][j];
};

}//namespace util::math
}//namespace util

ext: h

max_stars_repo_name: uphere-co/nlp-prototype

max_issues_repo_name: uphere-co/nlp-prototype

max_forks_repo_name: uphere-co/nlp-prototype


text:
//====---- Sudoku/Size.h                                              ----====//
//
// The size definition that rules them all.
//====--------------------------------------------------------------------====//
#pragma once

#include <gsl/gsl>

#include <limits>

namespace Sudoku
{
template<int N>
struct Size
{
	using index                 = gsl::index;
	static constexpr index base = N;           // 3 for default 9*9 board
	static constexpr index elem = base * base; // 9 for default
	static constexpr index full = elem * elem; // 81 for default

	static_assert(
		1 < base && base < elem && base < full && elem < full,
		"struct Size: Board size out of bounds.");
	static_assert(
		N < std::numeric_limits<index>::max() / base &&
			N < std::numeric_limits<index>::max() / elem &&
			N < std::numeric_limits<index>::max() / elem / base &&
			N < std::numeric_limits<index>::max() / full &&
			N < std::numeric_limits<index>::max() / full / base,
		"struct Size: Board size too large for system integer size.");
	// fails at N > 35 when int less than 32-bits.
};

template<int N>
inline constexpr gsl::index base_size = Size<N>::base; // default 3
template<int N>
inline constexpr gsl::index elem_size = Size<N>::elem; // default 9
template<int N>
inline constexpr gsl::index full_size = Size<N>::full; // default 81

} // namespace Sudoku

max_stars_repo_name: FeodorFitsner/fwkSudoku

max_issues_repo_name: FeodorFitsner/fwkSudoku

max_forks_repo_name: FeodorFitsner/fwkSudoku


text:
/*
 * Copyright (c) 1997-1999 Massachusetts Institute of Technology
 *
 * 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 */

/* This file was automatically generated --- DO NOT EDIT */
/* Generated on Sun Nov  7 20:44:36 EST 1999 */

#include <fftw-int.h>
#include <fftw.h>

/* Generated by: ./genfft -magic-alignment-check -magic-twiddle-load-all -magic-variables 4 -magic-loopi -twiddle 16 */

/*
 * This function contains 174 FP additions, 84 FP multiplications,
 * (or, 136 additions, 46 multiplications, 38 fused multiply/add),
 * 50 stack variables, and 64 memory accesses
 */
static const fftw_real K382683432 = FFTW_KONST(+0.382683432365089771728459984030398866761344562);
static const fftw_real K923879532 = FFTW_KONST(+0.923879532511286756128183189396788286822416626);
static const fftw_real K707106781 = FFTW_KONST(+0.707106781186547524400844362104849039284835938);

/*
 * Generator Id's :
 * $Id: exprdag.ml,v 1.41 1999/05/26 15:44:14 fftw Exp $
 * $Id: fft.ml,v 1.43 1999/05/17 19:44:18 fftw Exp $
 * $Id: to_c.ml,v 1.25 1999/10/26 21:41:32 stevenj Exp $
 */

void fftw_twiddle_16(fftw_complex *A, const fftw_complex *W, int iostride, int m, int dist)
{
     int i;
     fftw_complex *inout;
     inout = A;
     for (i = m; i > 0; i = i - 1, inout = inout + dist, W = W + 15) {
	  fftw_real tmp7;
	  fftw_real tmp91;
	  fftw_real tmp180;
	  fftw_real tmp193;
	  fftw_real tmp18;
	  fftw_real tmp194;
	  fftw_real tmp94;
	  fftw_real tmp177;
	  fftw_real tmp77;
	  fftw_real tmp88;
	  fftw_real tmp161;
	  fftw_real tmp128;
	  fftw_real tmp144;
	  fftw_real tmp162;
	  fftw_real tmp163;
	  fftw_real tmp164;
	  fftw_real tmp123;
	  fftw_real tmp143;
	  fftw_real tmp30;
	  fftw_real tmp152;
	  fftw_real tmp100;
	  fftw_real tmp136;
	  fftw_real tmp41;
	  fftw_real tmp153;
	  fftw_real tmp105;
	  fftw_real tmp137;
	  fftw_real tmp54;
	  fftw_real tmp65;
	  fftw_real tmp156;
	  fftw_real tmp117;
	  fftw_real tmp141;
	  fftw_real tmp157;
	  fftw_real tmp158;
	  fftw_real tmp159;
	  fftw_real tmp112;
	  fftw_real tmp140;
	  ASSERT_ALIGNED_DOUBLE;
	  {
	       fftw_real tmp1;
	       fftw_real tmp179;
	       fftw_real tmp6;
	       fftw_real tmp178;
	       ASSERT_ALIGNED_DOUBLE;
	       tmp1 = c_re(inout[0]);
	       tmp179 = c_im(inout[0]);
	       {
		    fftw_real tmp3;
		    fftw_real tmp5;
		    fftw_real tmp2;
		    fftw_real tmp4;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp3 = c_re(inout[8 * iostride]);
		    tmp5 = c_im(inout[8 * iostride]);
		    tmp2 = c_re(W[7]);
		    tmp4 = c_im(W[7]);
		    tmp6 = (tmp2 * tmp3) - (tmp4 * tmp5);
		    tmp178 = (tmp4 * tmp3) + (tmp2 * tmp5);
	       }
	       tmp7 = tmp1 + tmp6;
	       tmp91 = tmp1 - tmp6;
	       tmp180 = tmp178 + tmp179;
	       tmp193 = tmp179 - tmp178;
	  }
	  {
	       fftw_real tmp12;
	       fftw_real tmp92;
	       fftw_real tmp17;
	       fftw_real tmp93;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp9;
		    fftw_real tmp11;
		    fftw_real tmp8;
		    fftw_real tmp10;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp9 = c_re(inout[4 * iostride]);
		    tmp11 = c_im(inout[4 * iostride]);
		    tmp8 = c_re(W[3]);
		    tmp10 = c_im(W[3]);
		    tmp12 = (tmp8 * tmp9) - (tmp10 * tmp11);
		    tmp92 = (tmp10 * tmp9) + (tmp8 * tmp11);
	       }
	       {
		    fftw_real tmp14;
		    fftw_real tmp16;
		    fftw_real tmp13;
		    fftw_real tmp15;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp14 = c_re(inout[12 * iostride]);
		    tmp16 = c_im(inout[12 * iostride]);
		    tmp13 = c_re(W[11]);
		    tmp15 = c_im(W[11]);
		    tmp17 = (tmp13 * tmp14) - (tmp15 * tmp16);
		    tmp93 = (tmp15 * tmp14) + (tmp13 * tmp16);
	       }
	       tmp18 = tmp12 + tmp17;
	       tmp194 = tmp12 - tmp17;
	       tmp94 = tmp92 - tmp93;
	       tmp177 = tmp92 + tmp93;
	  }
	  {
	       fftw_real tmp71;
	       fftw_real tmp124;
	       fftw_real tmp87;
	       fftw_real tmp121;
	       fftw_real tmp76;
	       fftw_real tmp125;
	       fftw_real tmp82;
	       fftw_real tmp120;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp68;
		    fftw_real tmp70;
		    fftw_real tmp67;
		    fftw_real tmp69;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp68 = c_re(inout[15 * iostride]);
		    tmp70 = c_im(inout[15 * iostride]);
		    tmp67 = c_re(W[14]);
		    tmp69 = c_im(W[14]);
		    tmp71 = (tmp67 * tmp68) - (tmp69 * tmp70);
		    tmp124 = (tmp69 * tmp68) + (tmp67 * tmp70);
	       }
	       {
		    fftw_real tmp84;
		    fftw_real tmp86;
		    fftw_real tmp83;
		    fftw_real tmp85;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp84 = c_re(inout[11 * iostride]);
		    tmp86 = c_im(inout[11 * iostride]);
		    tmp83 = c_re(W[10]);
		    tmp85 = c_im(W[10]);
		    tmp87 = (tmp83 * tmp84) - (tmp85 * tmp86);
		    tmp121 = (tmp85 * tmp84) + (tmp83 * tmp86);
	       }
	       {
		    fftw_real tmp73;
		    fftw_real tmp75;
		    fftw_real tmp72;
		    fftw_real tmp74;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp73 = c_re(inout[7 * iostride]);
		    tmp75 = c_im(inout[7 * iostride]);
		    tmp72 = c_re(W[6]);
		    tmp74 = c_im(W[6]);
		    tmp76 = (tmp72 * tmp73) - (tmp74 * tmp75);
		    tmp125 = (tmp74 * tmp73) + (tmp72 * tmp75);
	       }
	       {
		    fftw_real tmp79;
		    fftw_real tmp81;
		    fftw_real tmp78;
		    fftw_real tmp80;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp79 = c_re(inout[3 * iostride]);
		    tmp81 = c_im(inout[3 * iostride]);
		    tmp78 = c_re(W[2]);
		    tmp80 = c_im(W[2]);
		    tmp82 = (tmp78 * tmp79) - (tmp80 * tmp81);
		    tmp120 = (tmp80 * tmp79) + (tmp78 * tmp81);
	       }
	       {
		    fftw_real tmp126;
		    fftw_real tmp127;
		    fftw_real tmp119;
		    fftw_real tmp122;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp77 = tmp71 + tmp76;
		    tmp88 = tmp82 + tmp87;
		    tmp161 = tmp77 - tmp88;
		    tmp126 = tmp124 - tmp125;
		    tmp127 = tmp82 - tmp87;
		    tmp128 = tmp126 + tmp127;
		    tmp144 = tmp126 - tmp127;
		    tmp162 = tmp124 + tmp125;
		    tmp163 = tmp120 + tmp121;
		    tmp164 = tmp162 - tmp163;
		    tmp119 = tmp71 - tmp76;
		    tmp122 = tmp120 - tmp121;
		    tmp123 = tmp119 - tmp122;
		    tmp143 = tmp119 + tmp122;
	       }
	  }
	  {
	       fftw_real tmp24;
	       fftw_real tmp96;
	       fftw_real tmp29;
	       fftw_real tmp97;
	       fftw_real tmp98;
	       fftw_real tmp99;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp21;
		    fftw_real tmp23;
		    fftw_real tmp20;
		    fftw_real tmp22;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp21 = c_re(inout[2 * iostride]);
		    tmp23 = c_im(inout[2 * iostride]);
		    tmp20 = c_re(W[1]);
		    tmp22 = c_im(W[1]);
		    tmp24 = (tmp20 * tmp21) - (tmp22 * tmp23);
		    tmp96 = (tmp22 * tmp21) + (tmp20 * tmp23);
	       }
	       {
		    fftw_real tmp26;
		    fftw_real tmp28;
		    fftw_real tmp25;
		    fftw_real tmp27;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp26 = c_re(inout[10 * iostride]);
		    tmp28 = c_im(inout[10 * iostride]);
		    tmp25 = c_re(W[9]);
		    tmp27 = c_im(W[9]);
		    tmp29 = (tmp25 * tmp26) - (tmp27 * tmp28);
		    tmp97 = (tmp27 * tmp26) + (tmp25 * tmp28);
	       }
	       tmp30 = tmp24 + tmp29;
	       tmp152 = tmp96 + tmp97;
	       tmp98 = tmp96 - tmp97;
	       tmp99 = tmp24 - tmp29;
	       tmp100 = tmp98 - tmp99;
	       tmp136 = tmp99 + tmp98;
	  }
	  {
	       fftw_real tmp35;
	       fftw_real tmp102;
	       fftw_real tmp40;
	       fftw_real tmp103;
	       fftw_real tmp101;
	       fftw_real tmp104;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp32;
		    fftw_real tmp34;
		    fftw_real tmp31;
		    fftw_real tmp33;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp32 = c_re(inout[14 * iostride]);
		    tmp34 = c_im(inout[14 * iostride]);
		    tmp31 = c_re(W[13]);
		    tmp33 = c_im(W[13]);
		    tmp35 = (tmp31 * tmp32) - (tmp33 * tmp34);
		    tmp102 = (tmp33 * tmp32) + (tmp31 * tmp34);
	       }
	       {
		    fftw_real tmp37;
		    fftw_real tmp39;
		    fftw_real tmp36;
		    fftw_real tmp38;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp37 = c_re(inout[6 * iostride]);
		    tmp39 = c_im(inout[6 * iostride]);
		    tmp36 = c_re(W[5]);
		    tmp38 = c_im(W[5]);
		    tmp40 = (tmp36 * tmp37) - (tmp38 * tmp39);
		    tmp103 = (tmp38 * tmp37) + (tmp36 * tmp39);
	       }
	       tmp41 = tmp35 + tmp40;
	       tmp153 = tmp102 + tmp103;
	       tmp101 = tmp35 - tmp40;
	       tmp104 = tmp102 - tmp103;
	       tmp105 = tmp101 + tmp104;
	       tmp137 = tmp101 - tmp104;
	  }
	  {
	       fftw_real tmp48;
	       fftw_real tmp108;
	       fftw_real tmp64;
	       fftw_real tmp115;
	       fftw_real tmp53;
	       fftw_real tmp109;
	       fftw_real tmp59;
	       fftw_real tmp114;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp45;
		    fftw_real tmp47;
		    fftw_real tmp44;
		    fftw_real tmp46;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp45 = c_re(inout[iostride]);
		    tmp47 = c_im(inout[iostride]);
		    tmp44 = c_re(W[0]);
		    tmp46 = c_im(W[0]);
		    tmp48 = (tmp44 * tmp45) - (tmp46 * tmp47);
		    tmp108 = (tmp46 * tmp45) + (tmp44 * tmp47);
	       }
	       {
		    fftw_real tmp61;
		    fftw_real tmp63;
		    fftw_real tmp60;
		    fftw_real tmp62;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp61 = c_re(inout[13 * iostride]);
		    tmp63 = c_im(inout[13 * iostride]);
		    tmp60 = c_re(W[12]);
		    tmp62 = c_im(W[12]);
		    tmp64 = (tmp60 * tmp61) - (tmp62 * tmp63);
		    tmp115 = (tmp62 * tmp61) + (tmp60 * tmp63);
	       }
	       {
		    fftw_real tmp50;
		    fftw_real tmp52;
		    fftw_real tmp49;
		    fftw_real tmp51;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp50 = c_re(inout[9 * iostride]);
		    tmp52 = c_im(inout[9 * iostride]);
		    tmp49 = c_re(W[8]);
		    tmp51 = c_im(W[8]);
		    tmp53 = (tmp49 * tmp50) - (tmp51 * tmp52);
		    tmp109 = (tmp51 * tmp50) + (tmp49 * tmp52);
	       }
	       {
		    fftw_real tmp56;
		    fftw_real tmp58;
		    fftw_real tmp55;
		    fftw_real tmp57;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp56 = c_re(inout[5 * iostride]);
		    tmp58 = c_im(inout[5 * iostride]);
		    tmp55 = c_re(W[4]);
		    tmp57 = c_im(W[4]);
		    tmp59 = (tmp55 * tmp56) - (tmp57 * tmp58);
		    tmp114 = (tmp57 * tmp56) + (tmp55 * tmp58);
	       }
	       {
		    fftw_real tmp113;
		    fftw_real tmp116;
		    fftw_real tmp110;
		    fftw_real tmp111;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp54 = tmp48 + tmp53;
		    tmp65 = tmp59 + tmp64;
		    tmp156 = tmp54 - tmp65;
		    tmp113 = tmp48 - tmp53;
		    tmp116 = tmp114 - tmp115;
		    tmp117 = tmp113 - tmp116;
		    tmp141 = tmp113 + tmp116;
		    tmp157 = tmp108 + tmp109;
		    tmp158 = tmp114 + tmp115;
		    tmp159 = tmp157 - tmp158;
		    tmp110 = tmp108 - tmp109;
		    tmp111 = tmp59 - tmp64;
		    tmp112 = tmp110 + tmp111;
		    tmp140 = tmp110 - tmp111;
	       }
	  }
	  {
	       fftw_real tmp107;
	       fftw_real tmp131;
	       fftw_real tmp202;
	       fftw_real tmp204;
	       fftw_real tmp130;
	       fftw_real tmp203;
	       fftw_real tmp134;
	       fftw_real tmp199;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp95;
		    fftw_real tmp106;
		    fftw_real tmp200;
		    fftw_real tmp201;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp95 = tmp91 - tmp94;
		    tmp106 = K707106781 * (tmp100 - tmp105);
		    tmp107 = tmp95 + tmp106;
		    tmp131 = tmp95 - tmp106;
		    tmp200 = K707106781 * (tmp137 - tmp136);
		    tmp201 = tmp194 + tmp193;
		    tmp202 = tmp200 + tmp201;
		    tmp204 = tmp201 - tmp200;
	       }
	       {
		    fftw_real tmp118;
		    fftw_real tmp129;
		    fftw_real tmp132;
		    fftw_real tmp133;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp118 = (K923879532 * tmp112) + (K382683432 * tmp117);
		    tmp129 = (K382683432 * tmp123) - (K923879532 * tmp128);
		    tmp130 = tmp118 + tmp129;
		    tmp203 = tmp129 - tmp118;
		    tmp132 = (K382683432 * tmp112) - (K923879532 * tmp117);
		    tmp133 = (K382683432 * tmp128) + (K923879532 * tmp123);
		    tmp134 = tmp132 - tmp133;
		    tmp199 = tmp132 + tmp133;
	       }
	       c_re(inout[11 * iostride]) = tmp107 - tmp130;
	       c_re(inout[3 * iostride]) = tmp107 + tmp130;
	       c_re(inout[15 * iostride]) = tmp131 - tmp134;
	       c_re(inout[7 * iostride]) = tmp131 + tmp134;
	       c_im(inout[3 * iostride]) = tmp199 + tmp202;
	       c_im(inout[11 * iostride]) = tmp202 - tmp199;
	       c_im(inout[7 * iostride]) = tmp203 + tmp204;
	       c_im(inout[15 * iostride]) = tmp204 - tmp203;
	  }
	  {
	       fftw_real tmp139;
	       fftw_real tmp147;
	       fftw_real tmp196;
	       fftw_real tmp198;
	       fftw_real tmp146;
	       fftw_real tmp197;
	       fftw_real tmp150;
	       fftw_real tmp191;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp135;
		    fftw_real tmp138;
		    fftw_real tmp192;
		    fftw_real tmp195;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp135 = tmp91 + tmp94;
		    tmp138 = K707106781 * (tmp136 + tmp137);
		    tmp139 = tmp135 + tmp138;
		    tmp147 = tmp135 - tmp138;
		    tmp192 = K707106781 * (tmp100 + tmp105);
		    tmp195 = tmp193 - tmp194;
		    tmp196 = tmp192 + tmp195;
		    tmp198 = tmp195 - tmp192;
	       }
	       {
		    fftw_real tmp142;
		    fftw_real tmp145;
		    fftw_real tmp148;
		    fftw_real tmp149;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp142 = (K382683432 * tmp140) + (K923879532 * tmp141);
		    tmp145 = (K923879532 * tmp143) - (K382683432 * tmp144);
		    tmp146 = tmp142 + tmp145;
		    tmp197 = tmp145 - tmp142;
		    tmp148 = (K923879532 * tmp140) - (K382683432 * tmp141);
		    tmp149 = (K923879532 * tmp144) + (K382683432 * tmp143);
		    tmp150 = tmp148 - tmp149;
		    tmp191 = tmp148 + tmp149;
	       }
	       c_re(inout[9 * iostride]) = tmp139 - tmp146;
	       c_re(inout[iostride]) = tmp139 + tmp146;
	       c_re(inout[13 * iostride]) = tmp147 - tmp150;
	       c_re(inout[5 * iostride]) = tmp147 + tmp150;
	       c_im(inout[iostride]) = tmp191 + tmp196;
	       c_im(inout[9 * iostride]) = tmp196 - tmp191;
	       c_im(inout[5 * iostride]) = tmp197 + tmp198;
	       c_im(inout[13 * iostride]) = tmp198 - tmp197;
	  }
	  {
	       fftw_real tmp155;
	       fftw_real tmp167;
	       fftw_real tmp188;
	       fftw_real tmp190;
	       fftw_real tmp166;
	       fftw_real tmp189;
	       fftw_real tmp170;
	       fftw_real tmp185;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp151;
		    fftw_real tmp154;
		    fftw_real tmp186;
		    fftw_real tmp187;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp151 = tmp7 - tmp18;
		    tmp154 = tmp152 - tmp153;
		    tmp155 = tmp151 + tmp154;
		    tmp167 = tmp151 - tmp154;
		    tmp186 = tmp41 - tmp30;
		    tmp187 = tmp180 - tmp177;
		    tmp188 = tmp186 + tmp187;
		    tmp190 = tmp187 - tmp186;
	       }
	       {
		    fftw_real tmp160;
		    fftw_real tmp165;
		    fftw_real tmp168;
		    fftw_real tmp169;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp160 = tmp156 + tmp159;
		    tmp165 = tmp161 - tmp164;
		    tmp166 = K707106781 * (tmp160 + tmp165);
		    tmp189 = K707106781 * (tmp165 - tmp160);
		    tmp168 = tmp159 - tmp156;
		    tmp169 = tmp161 + tmp164;
		    tmp170 = K707106781 * (tmp168 - tmp169);
		    tmp185 = K707106781 * (tmp168 + tmp169);
	       }
	       c_re(inout[10 * iostride]) = tmp155 - tmp166;
	       c_re(inout[2 * iostride]) = tmp155 + tmp166;
	       c_re(inout[14 * iostride]) = tmp167 - tmp170;
	       c_re(inout[6 * iostride]) = tmp167 + tmp170;
	       c_im(inout[2 * iostride]) = tmp185 + tmp188;
	       c_im(inout[10 * iostride]) = tmp188 - tmp185;
	       c_im(inout[6 * iostride]) = tmp189 + tmp190;
	       c_im(inout[14 * iostride]) = tmp190 - tmp189;
	  }
	  {
	       fftw_real tmp43;
	       fftw_real tmp171;
	       fftw_real tmp182;
	       fftw_real tmp184;
	       fftw_real tmp90;
	       fftw_real tmp183;
	       fftw_real tmp174;
	       fftw_real tmp175;
	       ASSERT_ALIGNED_DOUBLE;
	       {
		    fftw_real tmp19;
		    fftw_real tmp42;
		    fftw_real tmp176;
		    fftw_real tmp181;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp19 = tmp7 + tmp18;
		    tmp42 = tmp30 + tmp41;
		    tmp43 = tmp19 + tmp42;
		    tmp171 = tmp19 - tmp42;
		    tmp176 = tmp152 + tmp153;
		    tmp181 = tmp177 + tmp180;
		    tmp182 = tmp176 + tmp181;
		    tmp184 = tmp181 - tmp176;
	       }
	       {
		    fftw_real tmp66;
		    fftw_real tmp89;
		    fftw_real tmp172;
		    fftw_real tmp173;
		    ASSERT_ALIGNED_DOUBLE;
		    tmp66 = tmp54 + tmp65;
		    tmp89 = tmp77 + tmp88;
		    tmp90 = tmp66 + tmp89;
		    tmp183 = tmp89 - tmp66;
		    tmp172 = tmp157 + tmp158;
		    tmp173 = tmp162 + tmp163;
		    tmp174 = tmp172 - tmp173;
		    tmp175 = tmp172 + tmp173;
	       }
	       c_re(inout[8 * iostride]) = tmp43 - tmp90;
	       c_re(inout[0]) = tmp43 + tmp90;
	       c_re(inout[12 * iostride]) = tmp171 - tmp174;
	       c_re(inout[4 * iostride]) = tmp171 + tmp174;
	       c_im(inout[0]) = tmp175 + tmp182;
	       c_im(inout[8 * iostride]) = tmp182 - tmp175;
	       c_im(inout[4 * iostride]) = tmp183 + tmp184;
	       c_im(inout[12 * iostride]) = tmp184 - tmp183;
	  }
     }
}

static const int twiddle_order[] =
{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15};
fftw_codelet_desc fftw_twiddle_16_desc =
{
     "fftw_twiddle_16",
     (void (*)()) fftw_twiddle_16,
     16,
     FFTW_FORWARD,
     FFTW_TWIDDLE,
     352,
     15,
     twiddle_order,
};

ext: c

max_stars_repo_name: albertsgrc/ftdock-opt

max_issues_repo_name: albertsgrc/ftdock-opt

max_forks_repo_name: albertsgrc/ftdock-opt


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

#ifndef __GSL_STATISTICS_FLOAT_H__
#define __GSL_STATISTICS_FLOAT_H__

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

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

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_float_mean (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_variance (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_sd (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_variance_with_fixed_mean (const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_sd_with_fixed_mean (const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_absdev (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_skew (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_kurtosis (const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_lag1_autocorrelation (const float data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_float_covariance (const float data1[], const size_t stride1,const float data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_float_variance_m (const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_sd_m (const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_absdev_m (const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_skew_m_sd (const float data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_float_kurtosis_m_sd (const float data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_float_lag1_autocorrelation_m (const float data[], const size_t stride, const size_t n, const double mean);

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

/* DEFINED FOR FLOATING POINT TYPES ONLY */

GSL_EXPORT double gsl_stats_float_wmean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_wvariance (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_wsd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_wvariance_with_fixed_mean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_wsd_with_fixed_mean (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_float_wabsdev (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_wskew (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_float_wkurtosis (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_float_wvariance_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_float_wsd_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_float_wabsdev_m (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean);
GSL_EXPORT double gsl_stats_float_wskew_m_sd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean, const double wsd);
GSL_EXPORT double gsl_stats_float_wkurtosis_m_sd (const float w[], const size_t wstride, const float data[], const size_t stride, const size_t n, const double wmean, const double wsd);

/* END OF FLOATING POINT TYPES */

GSL_EXPORT double gsl_stats_float_pvariance (const float data1[], const size_t stride1, const size_t n1, const float data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_float_ttest (const float data1[], const size_t stride1, const size_t n1, const float data2[], const size_t stride2, const size_t n2);

GSL_EXPORT float gsl_stats_float_max (const float data[], const size_t stride, const size_t n);
GSL_EXPORT float gsl_stats_float_min (const float data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_float_minmax (float * min, float * max, const float data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_float_max_index (const float data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_float_min_index (const float data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_float_minmax_index (size_t * min_index, size_t * max_index, const float data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_float_median_from_sorted_data (const float sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_float_quantile_from_sorted_data (const float sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_FLOAT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/*
 * author: Achim Gaedke
 * created: May 2001
 * file: pygsl/src/ieeemodule.c
 * $Id: ieeemodule.c,v 1.7 2008/10/25 18:56:35 schnizer Exp $
 *
 */
#include <pygsl/error_helpers.h>
#include <pygsl/general_helpers.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_errno.h>
/*
 * constants definitions
 */

typedef struct {
  char *name;
  int flag;
} const_int_names;

static const_int_names ieeeConsts[] = {
  /* used as precision argument in set_mode*/
  {"single_precision",GSL_IEEE_SINGLE_PRECISION},
  {"double_precision",GSL_IEEE_DOUBLE_PRECISION},
  {"extended_precision",GSL_IEEE_EXTENDED_PRECISION},
  /* used as round argument in set_mode*/
  {"round_to_nearest",GSL_IEEE_ROUND_TO_NEAREST},
  {"round_down",GSL_IEEE_ROUND_DOWN},
  {"round_up",GSL_IEEE_ROUND_UP},
  {"round_to_zero",GSL_IEEE_ROUND_TO_ZERO},
  /* used as exception argument in set_mode*/
  {"mask_invalid",GSL_IEEE_MASK_INVALID},
  {"mask_denormalized",GSL_IEEE_MASK_DENORMALIZED},
  {"mask_division_by_zero",GSL_IEEE_MASK_DIVISION_BY_ZERO},
  {"mask_overflow",GSL_IEEE_MASK_OVERFLOW},
  {"mask_underflow",GSL_IEEE_MASK_UNDERFLOW},
  {"mask_all",GSL_IEEE_MASK_ALL},
  {"trap_inexact",GSL_IEEE_TRAP_INEXACT},
  /* used as type in ieee_*_rep */
  {"type_nan",GSL_IEEE_TYPE_NAN},
  {"type_inf",GSL_IEEE_TYPE_INF},
  {"type_normal",GSL_IEEE_TYPE_NORMAL},
  {"type_denormal",GSL_IEEE_TYPE_DENORMAL},
  {"type_zero",GSL_IEEE_TYPE_ZERO},
  /* sentinel */
  {NULL,0}
};

static void define_const_ints(PyObject* module)
{
  int i=0;

  while (ieeeConsts[i].name!=NULL)
    {
      PyModule_AddIntConstant(module,ieeeConsts[i].name,ieeeConsts[i].flag);
      i++;
    }
  return;
}

/*
 * method definitions
 */

static PyObject*
set_mode(PyObject *self, PyObject *args)
{
  int precision=2;
  int rounding=1;
  int exception_mask=GSL_IEEE_MASK_ALL;
  int flag;
  if (!PyArg_ParseTuple(args, "|iii", &precision, &rounding,  &exception_mask))
    return NULL;

  flag = gsl_ieee_set_mode(precision, rounding, exception_mask);
  if (GSL_SUCCESS != PyGSL_error_flag(flag))
    return NULL;

  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject*
env_setup(PyObject *self, PyObject *args)
{
  gsl_ieee_env_setup();
  Py_INCREF(Py_None);
  return Py_None;
}

static PyObject*
bin_repr(PyObject *self, PyObject *args)
{
  double x;
  gsl_ieee_double_rep r;
  if (!PyArg_ParseTuple(args,"d", &x))
      return NULL;
  gsl_ieee_double_to_rep(&x,&r);
  return Py_BuildValue("(isii)",r.sign,r.mantissa,r.exponent,r.type);
}

static PyObject*
ieee_isnan(PyObject *self, PyObject *arg)
{
     double tmp;
     if(PyGSL_PYFLOAT_TO_DOUBLE(arg, &tmp, NULL) != GSL_SUCCESS)
	  return NULL;
     return PyInt_FromLong(gsl_isnan(tmp));
}

static PyObject*
ieee_isinf(PyObject *self, PyObject *arg)
{
     double tmp;
     if(PyGSL_PYFLOAT_TO_DOUBLE(arg, &tmp, NULL) != GSL_SUCCESS)
	  return NULL;
     return PyInt_FromLong(gsl_isinf(tmp));

}

static PyObject*
ieee_finite(PyObject *self, PyObject *arg)
{
     double tmp;
     if(PyGSL_PYFLOAT_TO_DOUBLE(arg, &tmp, NULL) != GSL_SUCCESS)
	  return NULL;
     return PyInt_FromLong(gsl_finite(tmp));

}

static PyObject* ieee_nan(PyObject *self)
{
  return PyFloat_FromDouble(GSL_NAN);
}

static PyObject* ieee_neginf(PyObject *self)
{
  return PyFloat_FromDouble(GSL_NEGINF);
}

static PyObject* ieee_posinf(PyObject *self)
{
  return PyFloat_FromDouble(GSL_POSINF);
}

/* setup gsl mode and ieee modes from environment variable GSL_IEEE_MODE */

static const char env_setup_doc[] =
"Set the IEEE mode, the GSL library should use, using the GSL_IEEE_MODE\n\
environement variable\n\
You can also use os.environ['GSL_IEEE_MODE']='<desired_behaviour>' to\n\
set the variable\n";

static const char bin_repr_doc[] =
"takes a double and returns its  sign, mantissa, exponent and type\n";

static const char set_mode_doc[] =
"Sets the ieee float handling mode\n\
\tUsage:\n\
\tset_mode(precision, rounding, exception_mask)\n\
Each flag can be a combination of the constants defined in this module\n\
Raises  gsl_NoHardwareSupportError if support is not available.\n\
\n\
WARNING: Use with care! It can abort your programm abnormaly!";

static const char  ieee_nan_doc[] = "Returns not a number\n";
static const char  ieee_neginf_doc[] = "Returns a negative infinite\n";
static const char  ieee_posinf_doc[] = "Returns a positive infinite\n";
static PyMethodDef ieeeMethods[] = {
  {"set_mode", set_mode, METH_VARARGS, (char*)set_mode_doc},
  {"env_setup", env_setup, METH_NOARGS,  (char *)env_setup_doc},
  {"bin_repr", bin_repr, METH_VARARGS, NULL},
  /* tests on special IEEE-FP meanings */
  {"isnan",ieee_isnan,METH_O, NULL},
  {"isinf",ieee_isinf,METH_O, NULL},
  {"finite",ieee_finite,METH_O, NULL},
  /* some special ieee-numbers */
  {"nan",(PyCFunction)ieee_nan, METH_NOARGS, (char *)ieee_nan_doc},
  {"neginf",(PyCFunction)ieee_neginf, METH_NOARGS, (char *)ieee_neginf_doc},
  {"posinf",(PyCFunction)ieee_posinf, METH_NOARGS, (char *)ieee_posinf_doc},
  {NULL,     NULL}        /* Sentinel */
};

DL_EXPORT(void) initieee(void)
{
  PyObject* m;

  m=Py_InitModule("ieee", ieeeMethods);
  init_pygsl();
  define_const_ints(m);

  return;
}

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
#include "defines.h"
#include <string.h>
#include "lib.h"
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <sys/types.h>
#include <sys/stat.h>

extern REAL * opa;
extern int NTERMS;

/**
 * @param float * w
 * @param float * sig
 * @param float * spectro
 * @param int nspectro
 * @param float * spectra
 * @param float * d_spectra
 * @param float * beta
 * @param float * alpha
 *
 * */
int covarm(REAL *w,REAL *sig,float *spectro,int nspectro,REAL *spectra,REAL  *d_spectra,REAL *beta,REAL *alpha){

	int j,i,bt_nf,bt_nc,aux_nf,aux_nc;

	REAL AP[NTERMS*NTERMS*NPARMS],BT[NPARMS*NTERMS];

	REAL *BTaux,*APaux;

	//printf("\nVALORES DEL SIGMA SQUARE\n");

	for(j=0;j<NPARMS;j++){
		for(i=0;i<nspectro;i++){
			opa[i]= w[j]*(spectra[i+nspectro*j]-spectro[i+nspectro*j]);
		}

		BTaux=BT+(j*NTERMS);
		APaux=AP+(j*NTERMS*NTERMS);
		multmatrixIDLValueSigma(opa,nspectro,1,d_spectra+j*nspectro*NTERMS,NTERMS,nspectro,BTaux,&bt_nf,&bt_nc,sig+(nspectro*j)); //bt de tam NTERMS x 1
		multmatrix_transpose_sigma(d_spectra+j*nspectro*NTERMS,NTERMS,nspectro,d_spectra+j*nspectro*NTERMS,NTERMS,nspectro,APaux,&aux_nf,&aux_nc,w[j], sig+(nspectro*j));//ap de tam NTERMS x NTERMS

	}

	totalParcialf(BT,NPARMS,NTERMS,beta); //beta de tam 1 x NTERMS
	totalParcialMatrixf(AP,NTERMS,NTERMS,NPARMS,alpha); //alpha de tam NTERMS x NTERMS

	return 1;
}

/**
 * @param REAL * W
 * @param REAL * sig
 * @param float * spectro
 * @param int nspectro
 * @param REAL * spectra
 * @param REAL * d_spectra
 * @param REAL * beta
 * @param REAL * alpha
 *
 * */
int covarm2(REAL *w,REAL *sig,float *spectro,int nspectro,REAL *spectra,REAL  *d_spectra,REAL *beta,REAL *alpha){

	int j,i,h,k,aux_nf,aux_nc;

	REAL AP[NTERMS*NTERMS*NPARMS],BT[NPARMS*NTERMS];

	REAL *BTaux,*APaux;
	REAL sum,sum2;

	for(j=0;j<NPARMS;j++){
		for(i=0;i<nspectro;i++){
			opa[i]= w[j]*(spectra[i+nspectro*j]-spectro[i+nspectro*j]);
		}

		BTaux=BT+(j*NTERMS);
		APaux=AP+(j*NTERMS*NTERMS);

		for ( i = 0; i < NTERMS; i++){
		    for ( h = 0; h < NTERMS; h++){
				sum=0;
				if(i==0){
					sum2=0;
				}
				for ( k = 0;  k < nspectro; k++){
					REAL dAux = (*(d_spectra+j*nspectro*NTERMS+h*nspectro+k));
					sum += (*(d_spectra+j*nspectro*NTERMS+i*nspectro+k) * dAux ) * (w[j]/sig[nspectro*j+k]);

					if(i==0){
						sum2+= ((opa[k] * dAux  ))/sig[nspectro*j+k];
					}
				}

				APaux[NTERMS*i+h] = sum;
				if(i==0){
					BTaux[h] = sum2;
				}
     		}

		}
	}

	totalParcialf(BT,NPARMS,NTERMS,beta); //beta de tam 1 x NTERMS
	totalParcialMatrixf(AP,NTERMS,NTERMS,NPARMS,alpha); //alpha de tam NTERMS x NTERMS

	return 1;
}

/**
 * @param spectra: array with synthetic spectro
 * @param nspectro: size of spectro
 * @param spectro: original spectro
 * @param w: array of weight for I,Q,U,V
 * @param sig: array with sigma for I,Q,U,V
 * @param nfree: (nspectro * NPARMS) - NTERMS, NPARAMs is 4 and NTERMS 11.
 * */
REAL fchisqr(REAL * spectra,int nspectro,float *spectro,REAL *w,REAL *sig,REAL nfree){

	REAL TOT,dif1,dif2,dif3,dif4;
	REAL opa1,opa2,opa3,opa4;
	int i,j;

	TOT=0;
	opa1=0;
	opa2=0;
	opa3=0;
	opa4=0;
	for(i=0;i<nspectro;i++){
		dif1=spectra[i]-spectro[i];
		dif2=spectra[i+nspectro]-spectro[i+nspectro];
		dif3=spectra[i+nspectro*2]-spectro[i+nspectro*2];
		dif4=spectra[i+nspectro*3]-spectro[i+nspectro*3];

		opa1+= (((dif1*dif1)*w[0])/(sig[i]));
		opa2+= (((dif2*dif2)*w[1])/(sig[i+nspectro]));
		opa3+= (((dif3*dif3)*w[2])/(sig[i+nspectro*2]));
		opa4+= (((dif4*dif4)*w[3])/(sig[i+nspectro*3]));
	}
	TOT+= opa1+opa2+opa3+opa4;

	return TOT/nfree;

}

/*

	Multiplica la matriz a (tamaño naf,nac)
	por la matriz b (de tamaño nbf,nbc)
	al estilo IDL, es decir, filas de a por columnas de b,
	el resultado se almacena en resultOut (de tamaño fil,col)

	El tamaño de salida (fil,col) corresponde con (nbf,nac).

	El tamaño de columnas de b, nbc, debe de ser igual al de filas de a, naf.

*/

/**
 * @param REAL * a
 * @param int naf
 * @param int nac
 * @param REAL * b
 * @param int nbf
 * @param int nbc
 * @param REAL * result
 * @param int * fil
 * @param int * col
 * @param REAL value
 * */
int multmatrixIDLValue(REAL *a,int naf,int nac,REAL *b,int nbf,int nbc,REAL *result,int *fil,int *col,REAL value){

   int i,j,k;
   REAL sum;

	if(naf==nbc){
		(*fil)=nbf;
		(*col)=nac;

		for ( i = 0; i < nbf; i++){
		    for ( j = 0; j < nac; j++){
				sum=0;
				for ( k = 0;  k < naf; k++){
					sum += a[k*nac+j] * b[i*nbc+k];
				}
				result[((nac)*i)+j] = sum/value;
      		}
		}
		return 1;
	}
	else
		printf("\n \n Error en multmatrixIDLValue no coinciden nac y nbf!!!! ..\n\n");
	return 0;
}

/**
 *
 * */
int multmatrixIDLValueSigma(REAL *a,int naf,int nac,REAL *b,int nbf,int nbc,REAL *result,int *fil,int *col, REAL * sigma){

   int i,j,k;
   REAL sum;

	if(naf==nbc){
		(*fil)=nbf;
		(*col)=nac;

		for ( i = 0; i < nbf; i++){
				sum=0;
				for ( k = 0;  k < naf; k++){
						sum += (((a[k] * b[i*nbc+k])))/sigma[k];
				}
				result[i] = sum;
		}
		return 1;
	}
	else
		printf("\n \n Error en multmatrixIDLValue no coinciden nac y nbf!!!! ..\n\n");
	return 0;
}

/**
 * @param REAL * A
 * @param int f
 * @param int c
 * @param REAL * result
 * */
void totalParcialf(REAL * A, int f,int c,REAL * result){

	int i,j;
	REAL sum;
	for(i=0;i<c;i++){
		sum = 0;
		for(j=0;j<f;j++){
			sum+=A[j*c+i];
		}
		result[i] = sum;
	}
}

/**
 * @param REAL * A
 * @param int f
 * @param int c
 * @param int p
 * @param REAL * result
 * */
void totalParcialMatrixf(REAL * A, int f,int c,int p,REAL *result){

	int i,j,k;
	REAL sum;
	for(i=0;i<f;i++)
		for(j=0;j<c;j++){
			sum=0;
			for(k=0;k<p;k++)
				sum+=A[i*c+j+f*c*k];
			result[i*c+j] = sum;
		}
}

/**
 * @param PRECISION * a
 * @param int naf
 * @param int nac
 * @param PRECISION * b
 * @param int nbf
 * @param int nbc
 * @param PRECISION * result
 * @param int * fil
 * @param int * col
 *
 * Multiply matrix "a"(naf,nac) with matrix "b" (nbf,nbc). Algebraic matrix multiplication style, that is,
 * columns of "a" by rows of "b". The result is stored in "result"(fil,col).
 * The size of columns of "a", nac, must be the same of rows of "b", nbf.
 * */
int multmatrix(PRECISION *a,int naf,int nac, PRECISION *b,int nbf,int nbc,PRECISION *result,int *fil,int *col){

    int i,j,k;
    PRECISION sum;

	if(nac==nbf){
		(*fil)=naf;
		(*col)=nbc;

		for ( i = 0; i < naf; i++)
		    for ( j = 0; j < nbc; j++){
				sum=0;
				for ( k = 0;  k < nbf; k++){
					sum += a[i*nac+k] * b[k*nbc+j];
				}
				result[(*col)*i+j] = sum;
      		}
		return 1;
	}
	return 0;

}

/**
 * @param REAL * a
 * @param int naf
 * @param int nac
 * @param REAL *b
 * @param int nbf
 * @param int nbc
 * @param REAL * result
 * @param int * fil
 * @param int * col
 * @param REAL value
 * */
int multmatrix_transpose(REAL *a,int naf,int nac, REAL *b,int nbf,int nbc,REAL *result,int *fil,int *col,REAL value){

    int i,j,k;
    REAL sum;

	if(nac==nbc){
		(*fil)=naf;
		(*col)=nbf;

		for ( i = 0; i < naf; i++){
		    for ( j = 0; j < nbf; j++){
				sum=0;
				for ( k = 0;  k < nbc; k++){
					sum += a[i*nac+k] * b[j*nbc+k];
				}

				result[(*col)*i+j] = (sum)*value;
     		}

		}
		return 1;
	}else{
		printf("\n \n Error en multmatrix_transpose no coinciden nac y nbc!!!! ..\n\n");
	}

	return 0;
}

/**
 * @param REAL * a
 * @param int naf
 * @param int nac
 * @param REAL *b
 * @param int nbf
 * @param int nbc
 * @param REAL * result
 * @param int * fil
 * @param int * col
 * @param REAL weigth
 * @param REAL * sigma
 *
 * */
int multmatrix_transpose_sigma(REAL *a,int naf,int nac, REAL *b,int nbf,int nbc,REAL *result,int *fil,int *col,REAL weigth, REAL * sigma){

    int i,j,k;
    REAL sum;

	if(nac==nbc){
		(*fil)=naf;
		(*col)=nbf;

		for ( i = 0; i < naf; i++){
		    for ( j = 0; j < nbf; j++){
				sum=0;
				for ( k = 0;  k < nbc; k++){
						sum += (a[i*nac+k] * b[j*nbc+k]) * (weigth/sigma[k]);
				}

				result[(*col)*i+j] = sum;
     		}

		}
		return 1;
	}else{
		printf("\n \n Error en multmatrix_transpose no coinciden nac y nbc!!!! ..\n\n");
	}

	return 0;
}
//Media de un vector de longitud numl

/**
 * @param int nspectro
 */
int CalculaNfree(int nspectro)
{
	int nfree;
	nfree = 0;

	nfree = (nspectro * NPARMS) - NTERMS;

	return nfree;
}

/**
 * @param const char * path
 * */
int isDirectory(const char *path) {
   struct stat statbuf;
   if (stat(path, &statbuf) != 0)
       return 0;
   return S_ISDIR(statbuf.st_mode);
}

/**
 * @param tpuntero * cabeza
 * @param char * fileName
 * */
void insert_in_linked_list (tpuntero *cabeza, char * fileName){
    tpuntero nuevo;
    nuevo = malloc(sizeof(tnodo));
    strcpy(nuevo->d_name,fileName);
    nuevo->next = *cabeza;
    *cabeza = nuevo;
}

/**
 * @param tpuntero  cabeza
 * @param char * fileName
 * */
int checkNameInLista(tpuntero cabeza,char * fileName){
	int found = 0;
    while(cabeza != NULL && !found){ //Mientras cabeza no sea NULL
		if(strcmp(cabeza->d_name,fileName)==0)
			found = 1;
		else
        	cabeza = cabeza->next; //Pasamos al siguiente nodo
    }
	return found;
}

/**
 * @param tpuntero * cabeza
 * */

void deleteList(tpuntero *cabeza){
    tpuntero actual; //Puntero auxiliar para eliminar correctamente la lista

    while(*cabeza != NULL){ //Mientras cabeza no sea NULL
        actual = *cabeza; //Actual toma el valor de cabeza
        *cabeza = (*cabeza)->next; //Cabeza avanza 1 posicion en la lista
        free(actual); //Se libera la memoria de la posicion de Actual (el primer nodo), y cabeza queda apuntando al que ahora es el primero
    }
}


text:
/* BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*!
 * \file sim.c
 * \brief generate mocked 2d data.
 */
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <stddef.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

/////////////////////////////////////////////////////////////////////////////////
#ifdef SpecAstro
/* baselines for 3C 273 dataset observed by the GRAVITY (Nature, 2020, 563, 657),
 * unit is meter */
int n_base_sa_3c273 = 24;
double base_sa_3c273[]={
  -39.847287,   18.757261,
  -72.901618,  -13.864206,
  -105.767452, -59.897942,
  -32.453674,  -32.669331,
  -64.714063,  -78.769286,
  -32.064961,  -46.108491,
  -53.774627,   21.045077,
  -71.062525,  -10.421538,
  -79.670158,  -55.520144,
  -16.867270,  -31.447448,
  -25.151512,  -76.555774,
   -8.846202,  -45.116789,
  -54.064283,   19.430318,
  -86.349699,  -12.827258,
  -115.711019, -58.525511,
  -32.484928,  -32.180585,
  -61.579109,  -77.881063,
  -29.231494,  -45.704607,
  -52.779207,   21.321411,
  -83.468788,  -12.802497,
  -107.827509, -59.908816,
  -28.601448,  -33.449429,
  -50.961674,  -81.626333,
  -23.699820,  -46.979455
};
#endif
//////////////////////////////////////////////////////////////////////////////////////

void *model;

void sim()
{
  if(thistask != roottask)
    return;

  FILE *fp;
  char fname[200];
  int i, j, incr;

  sim_init();

  double *pm = (double *)model, error, fcon;

  smooth_init(parset.n_vel_recon, TransV);

  if(parset.flag_dim == -1)
  {
    /* note that here use sigma_hat = sigma/sqrt(tau) */
    printf("sim with ln(sigma) = %f and  ln(taud) = %f.\n", var_param[1], var_param[2]);
    reconstruct_con_from_varmodel(exp(var_param[1]), exp(var_param[2]), 1.0, 0.0);
  }
  else
  {
    con_scale = 1.0;
    line_scale = 1.0;
    line_error_mean = con_error_mean = 0.01;
    /* arguments: sigma_hat, tau, alapha, and syserr */
    create_con_from_random(0.03, 45.0, 1.0, 0.0);
  }
  calculate_con_rm(model);

  sprintf(fname, "%s/%s", parset.file_dir, "/data/sim_con_full.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  for(i=0; i<parset.n_con_recon; i++)
  {
    fprintf(fp, "%e %e %e\n", Tcon[i]*(1.0+parset.redshift), Fcon[i]/con_scale, Fcerrs[i]/con_scale);
  }
  fclose(fp);

  sprintf(fname, "%s/%s", parset.file_dir, "/data/sim_con.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  if(parset.flag_dim != -2)
  {
    gsl_interp_init(gsl_linear, Tcon, Fcon, parset.n_con_recon);
    for(i=0; i<n_con_data; i++)
    {
      //fprintf(fp, "%f %f %f\n", Tcon[i], Fcon[i]/con_scale, Fcerrs[i]/con_scale);
      fcon = gsl_interp_eval(gsl_linear, Tcon, Fcon, Tcon_data[i], gsl_acc);
      fprintf(fp, "%e %e %e\n", Tcon_data[i]*(1.0+parset.redshift),
            (fcon+gsl_ran_ugaussian(gsl_r)*con_error_mean)/con_scale, con_error_mean/con_scale);
    }
  }
  else
  {
    incr = fmax(0.5/(Tcon[1]-Tcon[0]), 1.0); //cadence to be 0.5day or increasement to be 1
    for(i=0; i<parset.n_con_recon; i+=incr)
    {
      if(Tcon[i] >= 0.0)
      {
        fprintf(fp, "%e %e %e\n", Tcon[i]*(1.0+parset.redshift), Fcon[i]/con_scale, Fcerrs[i]/con_scale);
      }
    }
  }
  fclose(fp);

  gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);
  transfun_1d_cal(model, 0);
  calculate_line_from_blrmodel(model, Tline, Fline, parset.n_line_recon);

  sprintf(fname, "%s/%s", parset.file_dir, "/data/sim_line.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  if(parset.flag_dim == -1)
  {
    error = line_error_mean * sqrt(n_line_data) * (Vline_data[1] - Vline_data[0]);
  }
  else
  {
    error = line_error_mean;
  }
  for(i=0; i<parset.n_line_recon; i++)
  {
    fprintf(fp, "%e %e %e\n", Tline[i]*(1.0+parset.redshift), Fline[i]/line_scale + gsl_ran_ugaussian(gsl_r)*error/line_scale,
      error/line_scale);
  }

  // output transfer function.
  sprintf(fname, "%s/%s", parset.file_dir, parset.tran_out_file);
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  for(i=0; i<parset.n_tau; i++)
  {
    fprintf(fp, "%e %e\n", TransTau[i], Trans1D[i]);
  }
  fclose(fp);

  transfun_2d_cal(model, TransV, Trans2D, parset.n_vel_recon, 1);
  calculate_line2d_from_blrmodel(model, Tline, TransV,
          Trans2D, Fline2d, parset.n_line_recon, parset.n_vel_recon);

  sprintf(fname, "%s/%s", parset.file_dir, "/data/sim_line2d.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  fprintf(fp, "# %d %d\n", parset.n_line_recon, parset.n_vel_recon);
  for(i=0; i<parset.n_line_recon; i++)
  {
    fprintf(fp, "# %f\n", Tline[i]*(1.0+parset.redshift));
    for(j=0; j<parset.n_vel_recon; j++)
    {
      fprintf(fp, "%e %e %e\n", TransW[j],
        (Fline2d[i*parset.n_vel_recon + j] + gsl_ran_ugaussian(gsl_r)*line_error_mean*0.3)/line_scale, line_error_mean/line_scale);
    }

    fprintf(fp, "\n");
  }
  fclose(fp);

  sprintf(fname, "%s/%s", parset.file_dir, "/data/sim_broadening.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }

  for(i=0; i<parset.n_line_recon; i++)
  {
    fprintf(fp, "%f %f\n", parset.InstRes * VelUnit, parset.InstRes_err * VelUnit);
  }
  fclose(fp);

  // output 2d transfer function
  sprintf(fname, "%s/%s", parset.file_dir, parset.tran2d_out_file);
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s\n", fname);
    exit(-1);
  }
  fprintf(fp, "# %d %d\n", parset.n_tau, parset.n_vel_recon);
  for(i=0; i<parset.n_tau; i++)
  {
    for(j=0; j<parset.n_vel_recon; j++)
    {
      fprintf(fp, "%e %e %e\n", TransV[j]*VelUnit, TransTau[i], Trans2D[i*parset.n_vel_recon + j]);
    }

    fprintf(fp, "\n");
  }
  fclose(fp);

#ifdef SpecAstro
  double *sa_pm;
  sa_pm = (double *)pm + num_params_blr;

  sa_smooth_init(parset.n_sa_vel_recon, vel_sa, parset.sa_InstRes);

  gen_sa_cloud_sample((void *)sa_pm, 3, 0);
  calculate_sa_sim_with_sample(pm, vel_sa, parset.n_sa_vel_recon, base_sa, parset.n_sa_base_recon,
                                   phase_sa, Fline_sa);

  sprintf(fname, "%s/%s", parset.file_dir, "data/sim_sa.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
    exit(0);
  }
  // output sa line
  fprintf(fp, "# %d %d %d\n", 1, parset.n_sa_vel_recon, parset.n_sa_base_recon);
  for(j=0; j<parset.n_sa_vel_recon; j++)
  {
    fprintf(fp, "%e %e %e\n", wave_sa[j], Fline_sa[j] + gsl_ran_ugaussian(gsl_r)*sa_line_error_mean, sa_line_error_mean);
  }
  fprintf(fp, "\n");
  for(i=0; i<parset.n_sa_base_recon; i++)
  {
    fprintf(fp, "# %f %f\n", base_sa[i*2], base_sa[i*2+1]);
    for(j=0; j<parset.n_sa_vel_recon; j++)
    {
      fprintf(fp, "%e %e %e\n", wave_sa[j],
       phase_sa[i*parset.n_sa_vel_recon + j]/(PhaseFactor * wave_sa[j]) + gsl_ran_ugaussian(gsl_r)*sa_phase_error_mean,
       sa_phase_error_mean);
    }
    fprintf(fp, "\n");
  }
  fclose(fp);
  sa_smooth_end();

#endif

  smooth_end();
  sim_end();
}

void sim_init()
{
  int i, j, idx;
  double dT, Tspan;
  double *pm, Rblr, mbh;

  switch(parset.flag_blrmodel)
  {
    case -1:
      num_params_blr_model = num_params_MyTransfun2d;
      transfun_1d_cal = transfun_1d_cal_mytransfun;
      transfun_2d_cal = transfun_2d_cal_mytransfun;
      break;

    case 0:
      num_params_blr_model = num_params_MyBLRmodel2d;
      gen_cloud_sample = gen_cloud_sample_mymodel;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 1:
      num_params_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
    case 2:
      num_params_blr_model = sizeof(BLRmodel2)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model2;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
    case 3:
      num_params_blr_model = sizeof(BLRmodel3)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model3;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
    case 4:
      num_params_blr_model = sizeof(BLRmodel4)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model4;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
    case 5:
      num_params_blr_model = sizeof(BLRmodel5)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model5;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 6:
      num_params_blr_model = sizeof(BLRmodel6)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model6;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 7:
      num_params_blr_model = sizeof(BLRmodel7)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model7;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 8:
      num_params_blr_model = sizeof(BLRmodel8)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model8;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 9:
      num_params_blr_model = sizeof(BLRmodel9)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model9;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    default:
      num_params_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_1d_cal = transfun_1d_cal_cloud;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
  }

#ifdef SpecAstro
  set_sa_blr_model();
  /* SA */
  num_params_sa_blr = num_params_sa_blr_model + num_params_sa_extpar;
  num_params_sa = num_params_sa_blr;

#endif

  /* use epoch-independent broadening */
  if(parset.flag_InstRes > 1)
  {
    num_params_res = 1;
    parset.InstRes = 220.0/VelUnit;
    parset.InstRes_err = 0.0;
  }
  else
  {
    parset.InstRes /= VelUnit;
    parset.InstRes_err = 0.0;
  }

  if(parset.flag_narrowline > 1)
  {
    printf("# set flag_narrowline to 1.\n");
    parset.flag_narrowline = 1;
  }
  else if(parset.flag_narrowline == 0)
  {
    parset.width_narrowline = 0.0;
  }

  parset.flag_linecenter = 0;
  num_params_linecenter = 0;

  parset.num_particles = 1;
  which_particle_update = 0;
  force_update = 1;
  which_parameter_update = -1;

  num_params_blr = num_params_blr_model + num_params_nlr
                 + num_params_res + num_params_linecenter + 1; /* include line sys err */
  num_params_var = num_params_drw + num_params_trend + num_params_difftrend + num_params_resp;

  num_params_blr_tot = num_params_blr;
#ifdef SpecAstro
  num_params_blr_tot  += num_params_sa_blr;
#endif
  num_params = num_params_blr_tot + num_params_var + parset.n_con_recon;

  /* index of A and Ag */
  idx_resp = num_params_blr_tot + num_params_drw + num_params_trend;
  /* index of different trend */
  idx_difftrend = idx_resp + num_params_resp;

  model = malloc(num_params * sizeof(double));
  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));

  /* setup parameters */
  pm = (double *)model;
  if(parset.flag_blrmodel == -1)
  {
    set_par_value_mytransfun_sim(pm);
  }
  else
  {
    set_par_value_sim(pm, parset.flag_blrmodel);
    set_par_fix_blrmodel();
    for(i=0; i<num_params_blr_model; i++)
    {
      if(par_fix[i] == 1)
      {
        pm[i] = par_fix_val[i];
      }
    }
  }

  /* spectral broadening, note this is a deviation from the input value */
  pm[num_params_blr_model + num_params_nlr ] = 0.0;

  pm[idx_resp + 0] = log(1.0); //A
  pm[idx_resp + 1] = 0.0;      //Ag

#ifdef SpecAstro
  double *sa_model = pm + num_params_blr;
  set_idx_par_mutual();
  set_par_value_sim(sa_model, parset.flag_sa_blrmodel);
  set_par_fix_sa_blrmodel();
  for(i=0; i<num_params_sa_blr_model; i++)
  {
    if(par_fix[num_params_blr + i] == 1)
    {
      pm[num_params_blr + i] = par_fix_val[num_params_blr + i];
    }
  }

  /* set the same mbh and inc */
  sa_model[idx_sa_par_mutual[0]] = pm[idx_rm_par_mutual[0]]; //mbh
  sa_model[idx_sa_par_mutual[1]] = pm[idx_rm_par_mutual[1]]; //inc

  sa_model[num_params_sa_blr_model]     = log(550.0); //DA
  sa_model[num_params_sa_blr_model + 1] = 0.0;        //PA
  sa_model[num_params_sa_blr_model + 2] = 0.0;        //FA
  sa_model[num_params_sa_blr_model + 3] = 0.0;        //CO
#endif

  Fcon = malloc(parset.n_con_recon * sizeof(double));
  Fcon_rm = malloc(parset.n_con_recon * sizeof(double));

  idx = get_idx_mbh_from_blrmodel();
  mbh = exp(pm[idx]);

  if(parset.flag_blrmodel != 8) /* model 8 is particular */
  {
    Rblr = exp(pm[0]);
  }
  else
  {
    Rblr = exp(pm[9]);
  }

  if(parset.flag_dim == -1)
  {
    Tspan = Tcon_data[n_con_data -1] - Tcon_data[0];

  /* set time array for continuum */
    Tcon_min = Tcon_data[0] - time_back_set - 10.0;

    Tcon_max = fmax(Tcon_data[n_con_data-1], Tline_data[n_line_data-1]) + fmax(0.05*Tspan, 10.0);
    dT = (Tcon_max - Tcon_min)/(parset.n_con_recon -1);

    for(i=0; i<parset.n_con_recon; i++)
    {
      Tcon[i] = Tcon_min + i*dT;
    }
  }
  else
  {
    Tspan = Rblr*10.0;

    rcloud_min_set = 0.0;
    rcloud_max_set = Tspan/2.0;
    if(parset.rcloud_max > 0.0)
      rcloud_max_set = fmin(rcloud_max_set, parset.rcloud_max);
    printf("RM rcloud_min_max_set: %f %f\n", rcloud_min_set, rcloud_max_set);

    time_back_set = 2.0*rcloud_max_set;
    Tcon_min = 0.0 - time_back_set - 10.0;
    Tcon_max = Tspan + 10.0;

    dT = (Tcon_max - Tcon_min)/(parset.n_con_recon -1);

    for(i=0; i<parset.n_con_recon; i++)
    {
      Tcon[i] = Tcon_min + i*dT;
    }
  }

  /* set Larr_rec */
  for(i=0;i<parset.n_con_recon;i++)
  {
    Larr_rec[i*nq + 0]=1.0;
    for(j=1; j<nq; j++)
      Larr_rec[i*nq + j] = pow(Tcon[i], j);
  }

  if(parset.flag_dim == -1)
  {
    parset.n_line_recon = n_line_data;
    parset.n_vel_recon = n_vel_data;
  }

  TransTau = malloc(parset.n_tau * sizeof(double));
  TransV = malloc(parset.n_vel_recon * sizeof(double));
  TransW = malloc(parset.n_vel_recon * sizeof(double));
  Trans1D = malloc(parset.n_tau * sizeof(double));
  Trans2D = malloc(parset.n_tau * parset.n_vel_recon * sizeof(double));
  Tline = malloc(parset.n_line_recon * sizeof(double));
  Fline = malloc(parset.n_line_recon * sizeof(double));
  Fline2d = malloc(parset.n_line_recon * parset.n_vel_recon * sizeof(double));

  if(parset.flag_dim == -1)
  {
    memcpy(Tline, Tline_data, n_line_data * sizeof(double));
    memcpy(TransV, Vline_data, n_vel_data * sizeof(double));
    memcpy(TransW, Wline_data, n_vel_data * sizeof(double));
  }
  else
  {
    Tline_min = 0.0;
    Tline_max = Tcon_max - 1.0;

    dT = (Tline_max - Tline_min)/(parset.n_line_recon - 1);

    for(i=0; i<parset.n_line_recon; i++)
    {
      Tline[i] = Tline_min + i*dT;
    }

    double vel_max_set, vel_min_set;
    vel_max_set = sqrt(pow(2.0*sqrt(mbh/Rblr), 2.0) + pow(2.0*parset.InstRes, 2.0));
    vel_min_set = - vel_max_set;
    double dVel = (vel_max_set- vel_min_set)/(parset.n_vel_recon -1.0);

    for(i=0; i<parset.n_vel_recon; i++)
    {
      TransV[i] = vel_min_set + dVel*i;
      TransW[i] = (1.0 + TransV[i]/C_Unit) * parset.linecenter * (1.0+parset.redshift);
    }
  }

  clouds_tau = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_weight = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_vel = malloc(parset.n_cloud_per_task * parset.n_vel_per_cloud * sizeof(double));

  if(parset.flag_save_clouds && thistask == roottask)
  {
    if(parset.n_cloud_per_task <= 1000)
      icr_cloud_save = 1;
    else
      icr_cloud_save = parset.n_cloud_per_task/1000;

    char fname[200];
    sprintf(fname, "%s/%s", parset.file_dir, parset.cloud_out_file);
    fcloud_out = fopen(fname, "w");
    if(fcloud_out == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s\n", fname);
      exit(-1);
    }
  }

#ifdef SpecAstro
  double saRblr;

  if(parset.flag_dim == -1)
  {
    parset.n_sa_vel_recon = n_vel_sa_data;
    parset.n_sa_base_recon = n_base_sa_data;
  }
  else
  {
    sa_flux_norm = 1.0;
    parset.n_sa_vel_recon = 40;
    if(parset.flag_gravity == 1)
      parset.n_sa_base_recon = n_base_sa_3c273;
    else
      parset.n_sa_base_recon = 20;

    sa_phase_error_mean = 0.01;
    sa_line_error_mean = 0.01;
  }

  vel_sa = malloc(parset.n_sa_vel_recon * sizeof(double));
  wave_sa = malloc(parset.n_sa_vel_recon * sizeof(double));
  Fline_sa = malloc(parset.n_sa_vel_recon * sizeof(double));
  base_sa = malloc(parset.n_sa_base_recon * 2 * sizeof(double));
  phase_sa = malloc(parset.n_sa_vel_recon * parset.n_sa_base_recon * sizeof(double));

  clouds_alpha = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_beta = malloc(parset.n_cloud_per_task * sizeof(double));

  workspace_phase = malloc(parset.n_sa_vel_recon * 3 * sizeof(double));

  if(parset.flag_sa_blrmodel != 8)
  {
    saRblr = exp(pm[num_params_blr]);
  }
  else
  {
    saRblr = exp(pm[num_params_blr + 9]);
  }
  rcloud_max_set = fmax(rcloud_max_set, saRblr * 5.0);
  printf("SA rcloud_min_max_set: %f %f\n", rcloud_min_set, rcloud_max_set);

  if(parset.flag_dim == -1)
  {
    memcpy(vel_sa, vel_sa_data, parset.n_sa_vel_recon * sizeof(double));
    memcpy(wave_sa, wave_sa_data, parset.n_sa_vel_recon * sizeof(double));

    memcpy(base_sa, base_sa_data, parset.n_sa_base_recon * 2 * sizeof(double));
  }
  else
  {
    double vel_max_set, vel_min_set;
    vel_max_set = sqrt(pow(2.0*sqrt(mbh/saRblr), 2.0) + pow(2.0*parset.sa_InstRes, 2.0));
    vel_min_set = - vel_max_set;
    double dVel = (vel_max_set- vel_min_set)/(parset.n_sa_vel_recon -1.0);
    for(i=0; i<parset.n_sa_vel_recon; i++)
    {
      vel_sa[i] = vel_min_set + dVel*i;
      wave_sa[i] = (1.0 + vel_sa[i]/C_Unit) * parset.sa_linecenter * (1.0+parset.redshift);
    }

    if(parset.flag_gravity == 1)
      memcpy(base_sa, base_sa_3c273, n_base_sa_3c273*2*sizeof(double));
    else
    {
      double phi;
      for(i=0; i<parset.n_sa_base_recon; i++)
      {
        phi = -PI/2.0 + PI/parset.n_sa_base_recon * i;
        base_sa[i*2+0] = 100.0*cos(phi);
        base_sa[i*2+1] = 100.0*sin(phi);
      }
    }

  }
#endif

  return;
}

void sim_end()
{
  free(model);
  free(par_fix);
  free(par_fix_val);
  free(Fcon);
  free(Fcon_rm);

  free(TransTau);
  free(TransV);
  free(TransW);
  free(Tline);
  free(Fline);
  free(Fline2d);
  free(Trans2D);
  free(Trans1D);

  free(clouds_tau);
  free(clouds_weight);
  free(clouds_vel);

  if(parset.flag_save_clouds && thistask == roottask)
  {
    fclose(fcloud_out);
  }

#ifdef SpecAstro

  free(vel_sa);
  free(wave_sa);
  free(Fline_sa);
  free(base_sa);
  free(phase_sa);

  free(clouds_alpha);
  free(clouds_beta);

  free(workspace_phase);
#endif
}

/*
 * set parameter values for either RM or SA
 */
void set_par_value_sim(double *pm, int flag_model)
{
  int i;
  switch(flag_model)
  {
    case 0:
      set_par_value_mymodel_sim(pm);
      break;

    case 1:
      i=0;
      pm[i++] = log(4.0);
      pm[i++] = 0.9;
      pm[i++] = 0.2;
      pm[i++] = cos(20.0/180.0*PI);
      pm[i++] = 40.0;
      pm[i++] = 0.0;
      pm[i++] = log(3.0);
      pm[i++] = 0.1;
      pm[i++] = 0.5;
      break;

    case 2:
      i=0;
      pm[i++] = log(4.0);
      pm[i++] = 0.9;
      pm[i++] = 0.2;
      pm[i++] = cos(20.0/180.0*PI);
      pm[i++] = 40.0;
      pm[i++] = 0.0;
      pm[i++] = log(3.0);
      pm[i++] = log(0.01);
      pm[i++] = log(0.1);
      break;

    case 3:
      i=0;
      pm[i++] = log(3.0);
      pm[i++] = log(5.0);
      pm[i++] = -1.0;
      pm[i++] = cos(20.0/180.0*PI);
      pm[i++] = 40.0;
      pm[i++] = 0.0;
      pm[i++] = log(3.0);
      pm[i++] = 0.5;
      pm[i++] = 0.5;
      break;

    case 4:
      i=0;
      pm[i++] = log(3.0);
      pm[i++] = log(5.0);
      pm[i++] = -1.0;
      pm[i++] = cos(20.0/180.0*PI);
      pm[i++] = 40.0;
      pm[i++] = 0.0;
      pm[i++] = log(3.0);
      pm[i++] = 0.5;
      pm[i++] = 0.5;
      break;

    case 5:
      i=0;
      pm[i++] = log(4.0);  //mu
      pm[i++] = 0.5;       //Fin
      pm[i++] = log(2.0);  //Fout
      pm[i++] = 1.5;       //alpha
      pm[i++] = cos(20.0/180.0*PI);      //inc
      pm[i++] = 40.0;      //opn
      pm[i++] = 0.5;       //k
      pm[i++] = 2.0;       //gam
      pm[i++] = 0.5;      //xi
      pm[i++] = log(2.0); //mbh
      pm[i++] = 0.5;      //fellip
      pm[i++] = 0.5;      //fflow
      pm[i++] = log(0.01);
      pm[i++] = log(0.1);
      pm[i++] = log(0.01);
      pm[i++] = log(0.1);
      pm[i++] = 0.0;      //theta_rot
      pm[i++] = log(0.001); //sig_turb
      break;

    case 6:
      i=0;
      pm[i++] = log(4.0);   // mu
      pm[i++] = 1.0;        // beta
      pm[i++] = 0.25;        // F
      pm[i++] = cos(20.0/180.0*PI); // inc
      pm[i++] = 40.0;       // opn
      pm[i++] = 0.0;        // kappa
      pm[i++] = 1.0;        // gamma
      pm[i++] = 1.0;        // obscuration
      pm[i++] = log(2.0);  //mbh
      pm[i++] = 0.5;       //fellip
      pm[i++] = 0.4;       //fflow
      pm[i++] = log(0.01); //
      pm[i++] = log(0.1);  //
      pm[i++] = log(0.01); //
      pm[i++] = log(0.1);  //
      pm[i++] = 0.0;       // theta_rot
      pm[i++] = log(0.001);  // sig_turb
      break;

    case 7:
      i=0;
      pm[i++] = log(4.0);  //mu
      pm[i++] = 0.8;       //beta
      pm[i++] = 0.1;       //F
      pm[i++] = cos(20.0/180.0*PI);//inc
      pm[i++] = 40.0;      //opn
      pm[i++] = 0.0;       //kappa
      pm[i++] = 1.0;       //gamma
      pm[i++] = 0.0;       //xi

      pm[i++] = 0.5;      //fsh
      pm[i++] = log(8.0); //mu_un
      pm[i++] = 1.2;      //beta_un
      pm[i++] = 0.1;      //F_un
      pm[i++] = 20.0;     //opn_un

      pm[i++] = log(2.0); //mbh
      pm[i++] = 0.5;      //fellip
      pm[i++] = 0.4;      //fflow
      pm[i++] = log(0.01);
      pm[i++] = log(0.1);
      pm[i++] = log(0.01);
      pm[i++] = log(0.1);
      pm[i++] = 0.0;      //theta_rot
      pm[i++] = 0.5;      //fellip_un
      pm[i++] = 0.4;      //fflow_un
      pm[i++] = log(0.001); // sig_turb
      break;

    case 8:
      i=0;
      pm[i++] = 50.0;   //theta_min
      pm[i++] = 20.0;   //dtheta_max
      pm[i++] = log(1.0); // r_min
      pm[i++] = 4.0;      // fr_max
      pm[i++] = 1.0;      // gamma
      pm[i++] = 1.0;      // alpha
      pm[i++] = -2.0;      // lambda
      pm[i++] = -0.5;      // k
      pm[i++] = 0.0;      // xi
      pm[i++] = log(30.0);     // Rv
      pm[i++] = log(20.0);     // Rblr, should larger than r_max=r_max * fr_max
      pm[i++] = cos(30.0/180.0*PI);     // inc
      pm[i++] = log(4.0); // mbh
      break;

    case 9:
      i=0;
      pm[i++] = log(4.0);   //mu
      pm[i++] = 1.0;        //beta
      pm[i++] = 0.2;        //F
      pm[i++] = cos(30.0/180.0*PI);     // inc
      pm[i++] = 30.0;       //opn
      pm[i++] = log(4.0);   // mbh
      break;
  }
  return;
}

/*
 * get index of mbh from a BLR model.
 *
 */
int get_idx_mbh_from_blrmodel()
{
  int idx = -1;
  switch(parset.flag_blrmodel)
  {
    case 0:
      idx = offsetof(MyBLRmodel, mbh);
      break;

    case 1:
      idx = offsetof(BLRmodel1, mbh);
      break;

    case 2:
      idx = offsetof(BLRmodel2, mbh);
      break;

    case 3:
      idx = offsetof(BLRmodel3, mbh);
      break;

    case 4:
      idx = offsetof(BLRmodel4, mbh);
      break;

    case 5:
      idx = offsetof(BLRmodel5, mbh);
      break;

    case 6:
      idx = offsetof(BLRmodel6, mbh);
      break;

    case 7:
      idx = offsetof(BLRmodel7, mbh);
      break;

    case 8:
      idx = offsetof(BLRmodel8, mbh);
      break;

    case 9:
      idx = offsetof(BLRmodel9, mbh);
      break;
  }
  return idx / sizeof(double);
}

max_stars_repo_name: LiyrAstroph/BRAINS

max_issues_repo_name: LiyrAstroph/BRAINS

max_forks_repo_name: LiyrAstroph/BRAINS