text:
/*
GENETIC - A simple genetic algorithm.

Copyright 2014, Javier Burguete Tolosa.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

1. Redistributions of source code must retain the above copyright notice, this
	list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above copyright notice,
	this list of conditions and the following disclaimer in the documentation
	and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY Javier Burguete Tolosa ``AS IS'' AND ANY EXPRESS
OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
EVENT SHALL Javier Burguete Tolosa OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/

/**
 * \file reproduction.c
 * \brief Source file to define the reproduction functions.
 * \author Javier Burguete Tolosa.
 * \copyright Copyright 2014 Javier Burguete Tolosa. All rights reserved.
 */
#define _GNU_SOURCE
#include <string.h>
#include <glib.h>
#include <gsl/gsl_rng.h>
//#include "config.h"
#include "bits.h"
#include "entity.h"
#include "reproduction.h"

void (*reproduction) (Entity *, Entity *, Entity *, unsigned int, gsl_rng *);
  ///< Pointer to the function to apply the reproduction operation.

/**
 * Function to mate two genotypes by single-point reproduction of each genome.
 */
static void
reproduction_singlepoints (Entity * father,     ///< Father.
                           Entity * mother,     ///< Mother.
                           Entity * son,        ///< Son.
                           unsigned int nbits,  ///< Genome nbits.
                           gsl_rng * rng)
                           ///< GSL random numbers generator.
{
  int i;
  i = gsl_rng_uniform_int (rng, nbits);
  bit_copy (son->genome, mother->genome, 0, 0, i);
  bit_copy (son->genome, father->genome, i, i, nbits - i);
}

/**
 * Function to mate two genotypes by double-point reproduction of each genome.
 */
static void
reproduction_doublepoints (Entity * father,     ///< Father.
                           Entity * mother,     ///< Mother.
                           Entity * son,        ///< Son.
                           unsigned int nbits,  ///< Genome nbits.
                           gsl_rng * rng)
                           ///< GSL random numbers generator.
{
  unsigned int i, j, k;
  i = gsl_rng_uniform_int (rng, nbits);
  j = gsl_rng_uniform_int (rng, nbits);
  if (i > j)
    k = i, i = j, j = k;
  bit_copy (son->genome, mother->genome, 0, 0, i);
  bit_copy (son->genome, father->genome, i, i, j - i);
  bit_copy (son->genome, mother->genome, j, j, nbits - j);
}

/**
 * Function to mate two genotypes by random mixing both genomes.
 */
static void
reproduction_mixing (Entity * father,   ///< Father.
                     Entity * mother,   ///< Mother.
                     Entity * son,      ///< Son.
                     unsigned int nbits,        ///< Genome nbits.
                     gsl_rng * rng)     ///< GSL random numbers generator.
{
  unsigned int i;
  for (i = 0; i < nbits; i++)
    {
      if (gsl_rng_uniform_int (rng, 2))
        {
          bit_get (father->genome, i) ?
            bit_set (son->genome, i) : bit_clear (son->genome, i);
        }
      else
        {
          bit_get (mother->genome, i) ?
            bit_set (son->genome, i) : bit_clear (son->genome, i);
        }
    }
}

/**
 * Function to select the reproduction operations.
 */
void
reproduction_init (unsigned int type)   ///< Type of reproduction operations.
{
  switch (type)
    {
    case REPRODUCTION_TYPE_SINGLEPOINTS:
      reproduction = &reproduction_singlepoints;
      break;
    case REPRODUCTION_TYPE_DOUBLEPOINTS:
      reproduction = &reproduction_doublepoints;
      break;
    default:
      reproduction = &reproduction_mixing;
    }
}

max_stars_repo_name: jburguete/genetic

max_issues_repo_name: jburguete/genetic

max_forks_repo_name: jburguete/genetic


text:
#ifndef CWANNIER_ECACHE_H
#define CWANNIER_ECACHE_H

#include <stdlib.h>
#include <stdbool.h>
#include <gsl/gsl_vector.h>
#include "input.h"

typedef struct {
    int na;
    int nb;
    int nc;
    int num_bands;
    int G_order[3];
    int G_neg[3];
    InputFn Efn;
    bool use_cache;
    gsl_vector **energies;
} EnergyCache;

#include "submesh.h"

EnergyCache* init_EnergyCache(int na, int nb, int nc, int num_bands, int G_order[3], int G_neg[3], InputFn Efn, bool use_cache);

void free_EnergyCache(EnergyCache *Ecache);

void energy_from_cache(EnergyCache *Ecache, int i, int j, int k, gsl_vector *energies);

void verify_sort(gsl_vector *energies);

#endif // CWANNIER_ECACHE_H

ext: h

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
#pragma once

#include <vector>
#include <string>
#include <memory>
#include <functional>
#include <gsl>

struct Vec3 { float x,y,z; };
struct Vec2 { float x,y; };

#include "types.h"
#include "automaticId.h"

max_stars_repo_name: racerxdl/bedhock

max_issues_repo_name: racerxdl/bedhock

max_forks_repo_name: racerxdl/bedhock


text:
#if !defined(REGEXMOVE_H_INCLUDED)
#define REGEXMOVE_H_INCLUDED

#include "Utils.h"

#include <filesystem>
#include <gsl/gsl>
#include <iosfwd>
#include <regex>
#include <string>

class RegExMove
{
public:
	static int usage(::std::ostream& strm, const ::std::string& progName,
		const char* pMsg);

	RegExMove(::gsl::span<const char*const> args);
	int run() const;

	RegExMove(const RegExMove&) = delete;
	RegExMove& operator=(const RegExMove&) = delete;
	RegExMove(RegExMove&&) = delete;
	RegExMove& operator=(RegExMove&&) = delete;

PRIVATE_EXCEPT_IN_TEST:
	using Path = ::std::filesystem::path;
	using DirEntry = ::std::filesystem::directory_entry;

	template <typename DirIter> void processDirectoryEntries() const;
	void processDirectoryEntry(DirEntry const& dirEntry) const;
	void renamePath(Path const& p) const;

	bool				m_caseSensitive;
	bool				m_renameDirectories;
	bool				m_renameFiles;
	bool				m_recursiveSearch;
	bool				m_verboseOutput;
	bool				m_allowOverwriteOnNameCollision;
	Path				m_rootDir;
	::std::string	m_patternStr;
	::std::regex	m_pattern;
	::std::string	m_replacement;
};

#endif // REGEXMOVE_H_INCLUDED

max_stars_repo_name: IanEmmons/CmdLineUtil

max_issues_repo_name: IanEmmons/CmdLineUtil

max_forks_repo_name: IanEmmons/CmdLineUtil


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

#ifndef __GSL_MULTISET_H__
#define __GSL_MULTISET_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.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

struct gsl_multiset_struct
{
  size_t n;
  size_t k;
  size_t *data;
};

typedef struct gsl_multiset_struct gsl_multiset;

gsl_multiset *gsl_multiset_alloc (const size_t n, const size_t k);
gsl_multiset *gsl_multiset_calloc (const size_t n, const size_t k);
void gsl_multiset_init_first (gsl_multiset * c);
void gsl_multiset_init_last (gsl_multiset * c);
void gsl_multiset_free (gsl_multiset * c);
int gsl_multiset_memcpy (gsl_multiset * dest, const gsl_multiset * src);

int gsl_multiset_fread (FILE * stream, gsl_multiset * c);
int gsl_multiset_fwrite (FILE * stream, const gsl_multiset * c);
int gsl_multiset_fscanf (FILE * stream, gsl_multiset * c);
int gsl_multiset_fprintf (FILE * stream, const gsl_multiset * c, const char *format);

size_t gsl_multiset_n (const gsl_multiset * c);
size_t gsl_multiset_k (const gsl_multiset * c);
size_t * gsl_multiset_data (const gsl_multiset * c);

int gsl_multiset_valid (gsl_multiset * c);
int gsl_multiset_next (gsl_multiset * c);
int gsl_multiset_prev (gsl_multiset * c);

INLINE_DECL size_t gsl_multiset_get (const gsl_multiset * c, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
size_t
gsl_multiset_get (const gsl_multiset * c, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= c->k)) /* size_t is unsigned, can't be negative */
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return c->data[i];
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_MULTISET_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
static char help[] = "A structured-grid Poisson solver using DMDA+KSP.\n\n";

#include <petsc.h>

extern PetscErrorCode formMatrix(DM, Mat);
extern PetscErrorCode formExact(DM, Vec);
extern PetscErrorCode formRHS(DM, Vec);

//STARTMAIN
int main(int argc,char **args) {
    PetscErrorCode ierr;
    DM            da;
    Mat           A;
    Vec           b,u,uexact;
    KSP           ksp;
    PetscReal     errnorm;
    DMDALocalInfo info;

    ierr = PetscInitialize(&argc,&args,NULL,help); if (ierr) return ierr;

    // change default 9x9 size using -da_grid_x M -da_grid_y N
    ierr = DMDACreate2d(PETSC_COMM_WORLD,
                 DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
                 9,9,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);

    // create linear system matrix A
    ierr = DMSetFromOptions(da); CHKERRQ(ierr);
    ierr = DMSetUp(da); CHKERRQ(ierr);
    ierr = DMCreateMatrix(da,&A); CHKERRQ(ierr);
    ierr = MatSetFromOptions(A); CHKERRQ(ierr);

    // create RHS b, approx solution u, exact solution uexact
    ierr = DMCreateGlobalVector(da,&b); CHKERRQ(ierr);
    ierr = VecDuplicate(b,&u); CHKERRQ(ierr);
    ierr = VecDuplicate(b,&uexact); CHKERRQ(ierr);

    // fill vectors and assemble linear system
    ierr = formExact(da,uexact); CHKERRQ(ierr);
    ierr = formRHS(da,b); CHKERRQ(ierr);
    ierr = formMatrix(da,A); CHKERRQ(ierr);

    // create and solve the linear system
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,u); CHKERRQ(ierr);

    // report on grid and numerical error
    ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uxact
    ierr = VecNorm(u,NORM_INFINITY,&errnorm); CHKERRQ(ierr);
    ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,
                "on %d x %d grid:  error |u-uexact|_inf = %g\n",
                info.mx,info.my,errnorm); CHKERRQ(ierr);

    VecDestroy(&u);  VecDestroy(&uexact);  VecDestroy(&b);
    MatDestroy(&A);  KSPDestroy(&ksp);  DMDestroy(&da);
    return PetscFinalize();
}
//ENDMAIN

//STARTMATRIX
PetscErrorCode formMatrix(DM da, Mat A) {
    PetscErrorCode ierr;
    DMDALocalInfo  info;
    MatStencil     row, col[5];
    PetscReal      hx, hy, v[5];
    PetscInt       i, j, ncols;

    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
    for (j = info.ys; j < info.ys+info.ym; j++) {
        for (i = info.xs; i < info.xs+info.xm; i++) {
            row.j = j;           // row of A corresponding to (x_i,y_j)
            row.i = i;
            col[0].j = j;        // diagonal entry
            col[0].i = i;
            ncols = 1;
            if (i==0 || i==info.mx-1 || j==0 || j==info.my-1) {
                v[0] = 1.0;      // on boundary: trivial equation
            } else {
                v[0] = 2*(hy/hx + hx/hy); // interior: build a row
                if (i-1 > 0) {
                    col[ncols].j = j;    col[ncols].i = i-1;
                    v[ncols++] = -hy/hx;
                }
                if (i+1 < info.mx-1) {
                    col[ncols].j = j;    col[ncols].i = i+1;
                    v[ncols++] = -hy/hx;
                }
                if (j-1 > 0) {
                    col[ncols].j = j-1;  col[ncols].i = i;
                    v[ncols++] = -hx/hy;
                }
                if (j+1 < info.my-1) {
                    col[ncols].j = j+1;  col[ncols].i = i;
                    v[ncols++] = -hx/hy;
                }
            }
            ierr = MatSetValuesStencil(A,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
        }
    }
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    return 0;
}
//ENDMATRIX

//STARTEXACT
PetscErrorCode formExact(DM da, Vec uexact) {
    PetscErrorCode ierr;
    PetscInt       i, j;
    PetscReal      hx, hy, x, y, **auexact;
    DMDALocalInfo  info;

    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
    ierr = DMDAVecGetArray(da, uexact, &auexact);CHKERRQ(ierr);
    for (j = info.ys; j < info.ys+info.ym; j++) {
        y = j * hy;
        for (i = info.xs; i < info.xs+info.xm; i++) {
            x = i * hx;
            auexact[j][i] = x*x * (1.0 - x*x) * y*y * (y*y - 1.0);
        }
    }
    ierr = DMDAVecRestoreArray(da, uexact, &auexact);CHKERRQ(ierr);
    return 0;
}

PetscErrorCode formRHS(DM da, Vec b) {
    PetscErrorCode ierr;
    PetscInt       i, j;
    PetscReal      hx, hy, x, y, f, **ab;
    DMDALocalInfo  info;

    ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
    hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
    ierr = DMDAVecGetArray(da, b, &ab);CHKERRQ(ierr);
    for (j=info.ys; j<info.ys+info.ym; j++) {
        y = j * hy;
        for (i=info.xs; i<info.xs+info.xm; i++) {
            x = i * hx;
            if (i==0 || i==info.mx-1 || j==0 || j==info.my-1) {
                ab[j][i] = 0.0;  // on boundary: 1*u = 0
            } else {
                f = 2.0 * ( (1.0 - 6.0*x*x) * y*y * (1.0 - y*y)
                    + (1.0 - 6.0*y*y) * x*x * (1.0 - x*x) );
                ab[j][i] = hx * hy * f;
            }
        }
    }
    ierr = DMDAVecRestoreArray(da, b, &ab); CHKERRQ(ierr);
    return 0;
}
//ENDEXACT

ext: c


text:
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>

int
main (void)
{
  int N = 4;
  double x[4] = {0.00, 0.10,  0.27,  0.30};
  double y[4] = {0.15, 0.70, -0.10,  0.15};
             /* Note: y[0] == y[3] for periodic data */

  gsl_interp_accel *acc = gsl_interp_accel_alloc ();
  const gsl_interp_type *t = gsl_interp_cspline_periodic;
  gsl_spline *spline = gsl_spline_alloc (t, N);

  int i; double xi, yi;

  printf ("#m=0,S=5\n");
  for (i = 0; i < N; i++)
    {
      printf ("%g %g\n", x[i], y[i]);
    }

  printf ("#m=1,S=0\n");
  gsl_spline_init (spline, x, y, N);

  for (i = 0; i <= 100; i++)
    {
      xi = (1 - i / 100.0) * x[0] + (i / 100.0) * x[N-1];
      yi = gsl_spline_eval (spline, xi, acc);
      printf ("%g %g\n", xi, yi);
    }

  gsl_spline_free (spline);
  gsl_interp_accel_free (acc);
  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <stdio.h>
#include <string.h>

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_matrix_double.h>

int main() {
	gsl_matrix* A;
	gsl_matrix* U;
    gsl_matrix* V;
    gsl_vector* S;
    gsl_vector* work;
    int i, j;
    int M, N;

    M = N = 2;

    work = gsl_vector_alloc(N);
    S = gsl_vector_alloc(N);
    A = gsl_matrix_alloc(M, N);
    V = gsl_matrix_alloc(N, N);

    gsl_matrix_set(A, 0, 0, -0.93);
    gsl_matrix_set(A, 0, 1,  1.80);
    gsl_matrix_set(A, 1, 0,  0);
    gsl_matrix_set(A, 1, 1,  0);

    gsl_linalg_SV_decomp(A, V, S, work);
    U = A;

    printf("S = [");
    for (i=0; i<N; i++)
        printf(" %g", gsl_vector_get(S, i));
    printf(" ]\n");

    printf("U = [");
    for (j=0; j<M; j++) {
        if (j)
            printf("    [");
        for (i=0; i<N; i++)
            printf(" %g", gsl_matrix_get(U, j, i));
        printf(" ]\n");
    }

    printf("V = [");
    for (j=0; j<N; j++) {
        if (j)
            printf("    [");
        for (i=0; i<N; i++)
            printf(" %g", gsl_matrix_get(V, j, i));
        printf(" ]\n");
    }

    gsl_matrix_set(A, 0, 0, -0.93);
    gsl_matrix_set(A, 0, 1,  1.80);
    gsl_matrix_set(A, 1, 0,  0);
    gsl_matrix_set(A, 1, 1,  0);

    gsl_linalg_SV_decomp_jacobi(A, V, S);

    printf("S = [");
    for (i=0; i<N; i++)
        printf(" %g", gsl_vector_get(S, i));
    printf(" ]\n");

    printf("U = [");
    for (j=0; j<M; j++) {
        if (j)
            printf("    [");
        for (i=0; i<N; i++)
            printf(" %g", gsl_matrix_get(U, j, i));
        printf(" ]\n");
    }

    printf("V = [");
    for (j=0; j<N; j++) {
        if (j)
            printf("    [");
        for (i=0; i<N; i++)
            printf(" %g", gsl_matrix_get(V, j, i));
        printf(" ]\n");
    }

	return 0;
}

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet


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

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_long_double.h>

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

__BEGIN_DECLS

typedef struct
{
  size_t size1;
  size_t size2;
  size_t tda;
  long double * data;
  gsl_block_long_double * block;
  int owner;
} gsl_matrix_long_double;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_view;

typedef _gsl_matrix_long_double_view gsl_matrix_long_double_view;

typedef struct
{
  gsl_matrix_long_double matrix;
} _gsl_matrix_long_double_const_view;

typedef const _gsl_matrix_long_double_const_view gsl_matrix_long_double_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_long_double *
gsl_matrix_long_double_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_long_double *
gsl_matrix_long_double_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_block (gsl_block_long_double * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_long_double *
gsl_matrix_long_double_alloc_from_matrix (gsl_matrix_long_double * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_long_double *
gsl_vector_long_double_alloc_row_from_matrix (gsl_matrix_long_double * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_long_double *
gsl_vector_long_double_alloc_col_from_matrix (gsl_matrix_long_double * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_long_double_free (gsl_matrix_long_double * m);

/* Views */

GSL_EXPORT
_gsl_matrix_long_double_view
gsl_matrix_long_double_submatrix (gsl_matrix_long_double * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_matrix_long_double_row (gsl_matrix_long_double * m, const size_t i);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_matrix_long_double_column (gsl_matrix_long_double * m, const size_t j);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_matrix_long_double_diagonal (gsl_matrix_long_double * m);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_matrix_long_double_subdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_EXPORT
_gsl_vector_long_double_view
gsl_matrix_long_double_superdiagonal (gsl_matrix_long_double * m, const size_t k);

GSL_EXPORT
_gsl_matrix_long_double_view
gsl_matrix_long_double_view_array (long double * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_long_double_view
gsl_matrix_long_double_view_array_with_tda (long double * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector (gsl_vector_long_double * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_long_double_view
gsl_matrix_long_double_view_vector_with_tda (gsl_vector_long_double * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_submatrix (const gsl_matrix_long_double * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_row (const gsl_matrix_long_double * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_column (const gsl_matrix_long_double * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_diagonal (const gsl_matrix_long_double * m);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_subdiagonal (const gsl_matrix_long_double * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_long_double_const_view
gsl_matrix_long_double_const_superdiagonal (const gsl_matrix_long_double * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array (const long double * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_array_with_tda (const long double * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector (const gsl_vector_long_double * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_long_double_const_view
gsl_matrix_long_double_const_view_vector_with_tda (const gsl_vector_long_double * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT long double   gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x);

GSL_EXPORT long double * gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_EXPORT const long double * gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_long_double_set_zero (gsl_matrix_long_double * m);
GSL_EXPORT void gsl_matrix_long_double_set_identity (gsl_matrix_long_double * m);
GSL_EXPORT void gsl_matrix_long_double_set_all (gsl_matrix_long_double * m, long double x);

GSL_EXPORT int gsl_matrix_long_double_fread (FILE * stream, gsl_matrix_long_double * m) ;
GSL_EXPORT int gsl_matrix_long_double_fwrite (FILE * stream, const gsl_matrix_long_double * m) ;
GSL_EXPORT int gsl_matrix_long_double_fscanf (FILE * stream, gsl_matrix_long_double * m);
GSL_EXPORT int gsl_matrix_long_double_fprintf (FILE * stream, const gsl_matrix_long_double * m, const char * format);

GSL_EXPORT int gsl_matrix_long_double_memcpy(gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);
GSL_EXPORT int gsl_matrix_long_double_swap(gsl_matrix_long_double * m1, gsl_matrix_long_double * m2);

GSL_EXPORT int gsl_matrix_long_double_swap_rows(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_long_double_swap_columns(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_long_double_swap_rowcol(gsl_matrix_long_double * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_long_double_transpose (gsl_matrix_long_double * m);
GSL_EXPORT int gsl_matrix_long_double_transpose_memcpy (gsl_matrix_long_double * dest, const gsl_matrix_long_double * src);

GSL_EXPORT long double gsl_matrix_long_double_max (const gsl_matrix_long_double * m);
GSL_EXPORT long double gsl_matrix_long_double_min (const gsl_matrix_long_double * m);
GSL_EXPORT void gsl_matrix_long_double_minmax (const gsl_matrix_long_double * m, long double * min_out, long double * max_out);

GSL_EXPORT void gsl_matrix_long_double_max_index (const gsl_matrix_long_double * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_long_double_min_index (const gsl_matrix_long_double * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_long_double_minmax_index (const gsl_matrix_long_double * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_long_double_isnull (const gsl_matrix_long_double * m);

GSL_EXPORT int gsl_matrix_long_double_add (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_EXPORT int gsl_matrix_long_double_sub (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_EXPORT int gsl_matrix_long_double_mul_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_EXPORT int gsl_matrix_long_double_div_elements (gsl_matrix_long_double * a, const gsl_matrix_long_double * b);
GSL_EXPORT int gsl_matrix_long_double_scale (gsl_matrix_long_double * a, const double x);
GSL_EXPORT int gsl_matrix_long_double_add_constant (gsl_matrix_long_double * a, const double x);
GSL_EXPORT int gsl_matrix_long_double_add_diagonal (gsl_matrix_long_double * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_long_double_get_row(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t i);
GSL_EXPORT int gsl_matrix_long_double_get_col(gsl_vector_long_double * v, const gsl_matrix_long_double * m, const size_t j);
GSL_EXPORT int gsl_matrix_long_double_set_row(gsl_matrix_long_double * m, const size_t i, const gsl_vector_long_double * v);
GSL_EXPORT int gsl_matrix_long_double_set_col(gsl_matrix_long_double * m, const size_t j, const gsl_vector_long_double * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
long double
gsl_matrix_long_double_get(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_long_double_set(gsl_matrix_long_double * m, const size_t i, const size_t j, const long double x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
long double *
gsl_matrix_long_double_ptr(gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (long double *) (m->data + (i * m->tda + j)) ;
}

extern inline
const long double *
gsl_matrix_long_double_const_ptr(const gsl_matrix_long_double * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const long double *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_LONG_DOUBLE_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


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

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

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

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

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


text:
#ifndef DERIVATIVE_H
#define DERIVATIVE_H

#include <memory>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_deriv.h>

#include "constants.h"

void deriv_discrete(gsl_vector *df, gsl_vector *f, double step);

#endif // DERIVATIVE_H

max_stars_repo_name: ichi-rika/glottal-inverse

max_issues_repo_name: ichi-rika/glottal-inverse

max_forks_repo_name: ichi-rika/glottal-inverse


text:
/* Siconos-Numerics, Copyright INRIA 2005-2011.
 * Siconos is a program dedicated to modeling, simulation and control
 * of non smooth dynamical systems.
 * Siconos is a free software; you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 * Siconos 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with Siconos; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 * Contact: Vincent ACARY, siconos-team@lists.gforge.inria.fr
 */

#ifndef SiconosBlas_H
#define SiconosBlas_H

#include "SiconosConfig.h"

#if defined(__cplusplus)
#include <algorithm>
using std::min;
using std::max;
#else
#define min(a,b) ((a)>(b)?(b):(a))
#define max(a,b) ((a) >= (b) ? (a) : (b))
#endif

#if defined(__cplusplus)
extern "C"
{
#endif

#if defined(HAS_MKL_CBLAS)
#include <mkl_cblas.h>
#elif defined(HAS_ATLAS_CBLAS)
#include <cblas.h>
#elif defined(HAS_ACCELERATE)
#include <Accelerate.h>
#else
#include <cblas.h>
#endif

#ifdef __cplusplus
}
#endif

#include "matrix_la.h"

#endif // SiconosBlas_H

max_stars_repo_name: siconos/siconos-deb

max_issues_repo_name: siconos/siconos-deb

max_forks_repo_name: siconos/siconos-deb


text:
/* ————————————————————————————————————————————————————————————————————————————————————————————————
Backwards elimination algorithm using Bidiag2.
———————————————————————————————————————————————————————————————————————————————————————————————————
This is a MEX function which uses the Bidiag2 PLS algorithm to estimate the cross-validated root-
mean-squared-error, RMSEcv, of an initial design matrix X, regressed upon y, together with every
possible single inactivation of a column. The function returns a matrix with RMSEcv as a function
of variable selection and PLS component where the first row corresponds to the initial variable
selection and every following row corresponds to the inactivation of variable 1 to n.

The function takes 5 inputs:
* Input 1: a [m x n] design matrix with observations as rows and variables as columns.
* Input 2: a [m x 1] vector with response values (multiple responses is not supported).
* Input 3: a [1 x 1] scalar 'A' specifying the maximum number of PLS components to calculate.
* Input 4: a [1 x 1] scalar 'MaxIter' specifying how many MonteCarlo CV iterations should be done.
* Input 5: a [1 x 1] scalar 'NoValObs' specifying how many unique responses should be in one validation fold.

The function outputs 1 variable:
* Output 1: a [(n+1) x A] matrix with RMSEcv. Variable selection along rows and components
along columns.

Example on how to compile and run from Matlab:
% Compile .C to .mexw64
>> mex -largeArrayDims -lmwblas BackwardsSelectionBidag2.c

% Run from Matlab when compiled:
>> X = rand(10000, 256);
>> y = rand(10000, 1);
>> A = 10;
>> MaxIters = 10;
>> Valobs = round(length(unique(y))/4);

>> [ RMSEcv ] = BackwardsSelectionBidag2( X , y, A, MaxIters, Valobs );

Example of compatible C compilers:
* Microsoft Visual C++ 2013 Professional (C)
* Microsoft Visual C++ 2015 Professional (C)
* Intel Parallel Studio XE 2017

Written 2017-08-14 by
petter.stefansson@nmbu.no
———————————————————————————————————————————————————————————————————————————————————————————————— */

#include <mex.h>	// Needed to communicate with matlab
#include <blas.h>	// Needed for blas functions
#include <string.h> // Needed to avoid compiler warning due to memcpy when using old compilers
#include <math.h>   // Needed to take the square-root (sqrt)
#include <time.h>   // Needed for counting CPU clock cycle which is used to set seed for rand()

/* ——————————————————————————————————— Function declarations ——————————————————————————————————— */
void PLS(const double *X, const double *y, int A, size_t m, size_t n, size_t p, double *beta,
	double *B, double *w, double *wn, double *W, double *rho, double *rhoi, double *d, double *tty,
	double *Xt, double *Ww, double *WWw, double *theta, double *thetai, double *Tt);

int randr(unsigned int min, unsigned int max);

void MarkAsVal(bool *IsVal, int NoValObs, size_t m, size_t *valm, size_t *trainm, const double *y);

void ExtractXandY(const double *X, const double *y, double *Xtrain, double *Xval, double *ytrain,
	double *yval, size_t m, size_t n, bool *IsVal);

void Pred(double *Xval, double *beta, double *yval, size_t n, size_t valm, int A, double *yhat,
	double *RMSEcv, int cvIter, int MaxIters, int ShavingIndex);

void SwapCols(double *Xtrain, double *Xval, size_t trainm, size_t valm, size_t n, int ColToShave);

/* ——————————————————————————————————— Matlab gateway start ———————————————————————————————————— */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	/* —————————————————————————— Variable type and name declaration ——————————————————————————— */
	const double *X, *y;
	double *Xtrain, *Xval, *ytrain, *yval, *beta, *yhat, *B, *w, *wn, *W, *rho, *rhoi, *d, *tty,
		   *Xt, *Ww, *WWw, *theta, *thetai, *Tt;

	int cvIter, NoValObs, MaxIters, A, ColToShave;
	bool *IsVal;
	size_t m, n, p, trainm, valm;

	/* ————————————————————  Get pointers to the input variables from Matlab ——————————————————— */
	X = mxGetPr(prhs[0]);			        // First input (X matrix).
	y = mxGetPr(prhs[1]);			        // Second input (Y vector).

	/* ————————————————————  Get input scalar values containing PLS settings ——————————————————— */
	A        = (int)mxGetScalar(prhs[2]);   // Third input (max number of components to calculate)
	MaxIters = (int)mxGetScalar(prhs[3]);   // Fourth input (Max cv-iters)
	NoValObs = (int)mxGetScalar(prhs[4]);   // Fifth input (Number of validation observations)

	/* ——————————————————————— Get the dimensions of the input variables ——————————————————————— */
	m = mxGetM(prhs[0]);                    // Number of rows in    X.
	n = mxGetN(prhs[0]);                    // Number of columns in X.
	p = mxGetN(prhs[1]);                    // Number of columns in y.

	/* ——————————————————————————— If input A is larger than n, let A = n —————————————————————— */
	if ((A > n) || (A < 1)) { A = n; }

	/* ——————————————————————————————— Specify Matlab outputs —————————————————————————————————— */
	double *RMSEcv;
	plhs[0] = mxCreateDoubleMatrix(n + 1, A, mxREAL);
	RMSEcv = mxGetPr(plhs[0]);

	/* ———— Allocate space for variables that are used multiple times and can be overwritten ——— */
	IsVal     = (bool*)malloc(sizeof(bool)     * m                  ); // [m-by-1]
	Xtrain    = (double*)malloc(sizeof(double) * m    *   n         ); // [m-by-n] (overallocated)
	Xval      = (double*)malloc(sizeof(double) * m    *   n         ); // [m-by-n] (overallocated)
	ytrain    = (double*)malloc(sizeof(double) * m                  ); // [m-by-1] (overallocated)
	yval      = (double*)malloc(sizeof(double) * m                  ); // [m-by-1] (overallocated)
	beta      = (double*)malloc(sizeof(double) * n    *   A         ); // [n-by-A]
	yhat      = (double*)malloc(sizeof(double) * m    *   A         ); // [m-by-A] (overallocated)
	B         = (double*)malloc(sizeof(double) * A    *   2         ); // [A-by-2]
	w         = (double*)malloc(sizeof(double) * n                  ); // [n-by-1]
	wn        = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	W         = (double*)malloc(sizeof(double) * n    * A           ); // [n-by-A]
    rho       = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	rhoi      = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	d         = (double*)malloc(sizeof(double) * n                  ); // [n-by-1]
	tty       = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	Xt        = (double*)malloc(sizeof(double) * n                  ); // [n-by-1]
	Ww        = (double*)malloc(sizeof(double) * A                  ); // [A-by-1] (overallocated)
	WWw       = (double*)malloc(sizeof(double) * n                  ); // [n-by-1]
	theta     = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	thetai    = (double*)malloc(sizeof(double)                      ); // [1-by-1]
	Tt        = (double*)malloc(sizeof(double) * A                  ); // [A-by-1] (overallocated)

	/* Before starting set the seed of the RNG to the number of clock cycles since start.		 */
	srand(clock());

	/* CV-loop starts here																	     */
	for (cvIter = 0; cvIter < MaxIters; cvIter++) {

		/* ————————————— Mark rows of X as either validation rows or training rows ————————————— */
		MarkAsVal(IsVal, NoValObs, m, &valm, &trainm, y);

		/* —————————— Extract data from X and y and place in Xtrain/Xval/ytrain/yval ——————————— */
		ExtractXandY(X, y, Xtrain, Xval, ytrain, yval, m, n, IsVal);

		/* ————————————————————————————————————————————————————————————————————————————————————— */
		/* Shaving step 1. Evaluate matrix with all variables included.							 */
		/* ————————— Call PLS function to estimate Beta using Xtrain & ytrain —————————————————— */
		PLS( Xtrain, ytrain, A, trainm, n, p, beta, B, w, wn, W, rho, rhoi, d, tty, Xt, Ww, WWw, theta, thetai, Tt);
		/* ——————— Use beta to estimate error when predicting yval using Xval —————————————————— */
		Pred( Xval, beta, yval, n, valm, A, yhat, RMSEcv, cvIter, MaxIters, 0);

		/* ————————————————————————————————————————————————————————————————————————————————————— */
		/* Shaving step 2. Evaluate matrix with all except the last variable included.			 */
		/* ————————— Call PLS function to estimate Beta using Xtrain & ytrain —————————————————— */
		PLS(Xtrain, ytrain, A, trainm, n - 1, p, beta, B, w, wn, W, rho, rhoi, d, tty, Xt, Ww, WWw, theta, thetai, Tt);
		/* ——————— Use beta to estimate error when predicting yval using Xval —————————————————— */
		Pred(Xval, beta, yval, n - 1, valm, A, yhat, RMSEcv, cvIter, MaxIters, n);
		/* ————————————————————————————————————————————————————————————————————————————————————— */

		/* ————————————————————————————————————————————————————————————————————————————————————— */
		/* Shaving step 3. Keep last column inactive and swap it one by one with cols 0 : n-1    */
		for (ColToShave = 0; ColToShave < (n - 1); ColToShave++) {

			/* Place ColToShave at the end, effectecly inactivating it.							 */
			SwapCols(Xtrain, Xval, trainm, valm, n, ColToShave + 1);
			/* ———————— Call PLS function to estimate Beta using Xtrain & ytrain ——————————————— */
			PLS(Xtrain, ytrain, A, trainm, n - 1, p, beta, B, w, wn, W, rho, rhoi, d, tty, Xt, Ww, WWw, theta, thetai, Tt);
			/* —————— Use beta to estimate error when predicting yval using Xval ——————————————— */
			Pred(Xval, beta, yval, n - 1, valm, A, yhat, RMSEcv, cvIter, MaxIters, ColToShave + 1);
		}
	}

	/* —————————————————————————————————— Free allocated memory ———————————————————————————————— */
	free(IsVal);
	free(Xtrain);
	free(Xval);
	free(ytrain);
	free(yval);
	free(beta);
	free(yhat);
	free(B);
	free(w);
	free(wn);
	free(W);
	free(rho);
	free(rhoi);
	free(d);
	free(tty);
	free(Xt);
	free(Ww);
	free(WWw);
	free(theta);
	free(thetai);
	free(Tt);
}

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for estimating Beta using the bidag2 algorithm										 */
void PLS(const double *X, const double *y, int A, size_t m, size_t n, size_t p, double *beta,
	double *B, double *w, double *wn, double *W, double *rho, double *rhoi, double *d,
	double *tty, double *Xt, double *Ww, double *WWw, double *theta, double *thetai, double *Tt){

	/* ——————————————————————————  Variable type and name declaration —————————————————————————— */
	double *t;
	double *T;
	double *Xw;
	double *TTt;

	int a, row;
	size_t st_a;

	mwSignedIndex IntConstOne = 1;
	double one = 1.0, zero = 0.0, none = -1.0;

	/* ————————————————— Allocate memory for variables used in the calculation ————————————————— */
	t      = (double*)malloc(sizeof(double)  * m                  ); // [m-by-1]
	T      = (double*)malloc(sizeof(double)  * m        * A       ); // [m-by-A]
	Xw     = (double*)malloc(sizeof(double)  * m                  ); // [m-by-1]
	TTt    = (double*)malloc(sizeof(double)  * m                  ); // [m-by-1]

	/* ————————————————————————————————————————————————————————————————————————————————————————— */
	/* ——————————————————————————————— Start of calculation part ——————————————————————————————— */
	/* ——————————————————————————————————————————————————————————————————————————————————————————
	Step 1: w = X'*y; 																             */
	/// X         = [m-by-n]
	/// y         = [m-by-1]
	/// w         = [n-by-1]

	// Matrix-Vector multiplication -> dgemv.
	dgemv("T", &m, &n, &one, X, &m, y, &IntConstOne, &zero, w,&IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 2: w = w / norm(w); 																     */
	/// w         = [n-by-1]
	/// wn        = [1-by-1]

	// Vector norm -> dnrm2.
	*wn = dnrm2(&n, w, &IntConstOne);
	// Inverse resulting scalar to enable multiplication instead of division.
	*wn = 1 / (*wn);
	// Vector-scalar multiplication -> dscal.
	dscal(&n, wn, w, &IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 3: W(:,CurrentComp) = w;																 */
	/// W         = [n-by-A]
	/// w         = [n-by-1]

	// Copy n doubles from w to W.
	memcpy(W, w, sizeof(double) * n);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 4: t = X*w;																             */
	/// X         = [m-by-n]
	/// w         = [n-by-1]
	/// t         = [m-by-1]

	// Matrix-Vector multiplication -> dgemv.
	dgemv("N", &m, &n, &one, X, &m, w, &IntConstOne, &zero, t, &IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 5: rho = norm(t);																         */
	/// t         = [m-by-1]
	/// rho       = [1-by-1]

	// Vector norm->dnrm2.
	*rho = dnrm2(&m, t, &IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 6: t = t / rho;																         */
	/// t         = [m-by-1]
	/// rho       = [1-by-1]
	/// rhoi      = [1-by-1]

	// Create inverse rho (rhoi) scalar to enable multiplication instead of division.
	*rhoi = 1 / (*rho);
	// Vector-scalar multiplication ->dscal.
	dscal(&m, rhoi, t, &IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 7: T(:,CurrentComp) = t;														         */
	/// T         = [m-by-A]
	/// t         = [m-by-p]

	// Copy m doubles from t to T.
	memcpy(T, t, sizeof(double) * m);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 8: B(1, 1) = rho;														                 */
	/// B         = [A-by-2]
	/// rho       = [1-by-1]

	B[0] = *rho;

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 9: d = w / rho;													                     */
	/// w         = [n-by-1]
	/// rho       = [1-by-1]
	/// rhoi      = [1-by-1]
	/// d         = [n-by-1]

	// Make d = w
	memcpy(d, w, sizeof(double) * n);

	// Then scale it with the inversed version of rho to make it w/rho.
	dscal(&n, rhoi, d, &IntConstOne);

	/*—————————————————————————————————————————————————————————————————————————————————————————————
	Step 9: beta(:, 1) = (t'*y)*d;													             */
	/// t         = [m-by-1]
	/// d         = [n-by-1]
	/// y         = [m-by-1]
	/// tty       = [1-by-1]
	/// beta(:,1) = [n-by-1]

	// (t'*y) vector-vector product -> ddot.
	*tty = ddot(&m, t, &IntConstOne, y, &IntConstOne);

	// (d*tty) vector-scalar product. (Here treated as matrix-matrix).
	dgemm("N", "N", &n, &p, &p, &one, d, &n, tty, &p, &zero, beta, &n);

	//* Start of component loop.															     */
	for (a = 1; a < A; a++) {
		/* size_t version of a (preferred by BLAS instead of int). 								 */
		st_a = (size_t)a;

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 10: w = X'*t - rho*w;													             */
		/// w         = [n-by-1]
		/// rho       = [1-by-1]
		/// X         = [m-by-n]
		/// t         = [m-by-1]
		/// Xt        = [n-by-1]

		// Xt = X'*t matrix-vector multiplication -> dgemv.
		dgemv("T", &m, &n, &one, X, &m, t, &IntConstOne, &zero, Xt, &IntConstOne);

		// The value of rho wont be used again, overwrite with -1*rho.
		*rho = -1 * (*rho);

		// w = rho*w scalar-vector multiplication -> dscal.
		dscal(&n, rho, w, &IntConstOne);

		// w = 1*Xt + w vector times constant plus vector -> daxpy.
		daxpy(&n, &one, Xt, &IntConstOne, w, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 11: w = w - W(:,1:a)*(W(:,1:a)'*w);											     */
		/// W(:,1:a)   = [n-by-a]
		/// w          = [n-by-1]
		/// Ww         = [a-by-1]
		/// WWw        = [n-by-1]

		// Ww = (W(:,1:a)'*w) matrix-vector multiplication -> dgemm.
		dgemv("T", &n, &st_a, &one, W, &n, w, &IntConstOne, &zero, Ww, &IntConstOne);

		// WWw = W(:,1:a)*(W(:,1:a)'*w)  matrix-vector multiplication -> dgemm.
		dgemv("N", &n, &st_a, &one, W, &n, Ww, &IntConstOne, &zero, WWw, &IntConstOne);

		//  w ←  w + -1 * WWw constant times a vector plus a vector -> daxpy.
		daxpy(&n, &none, WWw, &IntConstOne, w, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 12: theta = norm(w);													             */
		/// w         = [n-by-1]
		/// theta     = [1-by-1]

		// Vector norm -> dnrm2.
		*theta = dnrm2(&n, w, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 13: w = w / theta;													                 */
		/// w         = [n-by-1]
		/// theta     = [1-by-1]

		// w = w * (1/theta) scalar-vector multiplication -> dscal.
		*thetai = 1 / (*theta);
		dscal(&n, thetai, w, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 14: W(:,a) = w;													                 */
		/// w         = [n-by-1]
		/// W         = [n-by-A]

		// Copy n elemments from w into W at offset n * a.
		memcpy(W + n * a, w, sizeof(double) * n);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 15: t = X*w - theta*t;													             */
		/// X         = [m-by-n]
		/// w         = [n-by-1]
		/// theta     = [1-by-1]
		/// t         = [m-by-1]
		/// Xw        = [m-by-1]

		// Xw = X*w = matrix vector multiplication -> dgemv.
		dgemv("N", &m, &n, &one, X, &m, w, &IntConstOne, &zero, Xw, &IntConstOne);

		// t = Xw - theta*t Better to solve with daxpy?
		//for (row = 0; row < m; row++){
	    //	t[row] = Xw[row] - (*theta) * t[row];
		//}
		dscal(&m, theta, t, &IntConstOne);
		dscal(&m, &none, t, &IntConstOne);
		daxpy(&m, &one, Xw, &IntConstOne, t, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 16: t = t - T(:,1:a)*(T(:,1:a)'*t);												 */
		/// t         = [m-by-1]
		/// T(:,1:a)  = [m-by-a]
		/// Tt        = [a-by-1]
		/// TTt       = [m-by-1]

		// Tt = (T(:,1:a)'*t) matrix-vector multiplication -> dgemv.
		dgemv("T", &m, &st_a, &one, T, &m, t, &IntConstOne, &zero, Tt, &IntConstOne);

		// TTt = T(:,1:a)*(T(:,1:a)'*t) matrix-vector multiplication -> dgemv.
		dgemv("N", &m, &st_a, &one, T, &m, Tt, &IntConstOne, &zero, TTt, &IntConstOne);

		//  t ←  t + -1 * TTt constant times a vector plus a vector -> daxpy.
		daxpy(&m, &none, TTt, &IntConstOne, t, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 17: rho = norm(t);													                 */
		/// t         = [m-by-1]
		/// rho       = [1-by-1]

		// Vector norm -> dnrm2.
		*rho = dnrm2(&m, t, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 18: t = t/rho;													                     */
		/// t         = [m-by-1]
		/// rho       = [1-by-1]

		// Update inverse version of rho.
		*rhoi = 1 / (*rho);

		// t = t * (1/rho) scalar-vector multiplication -> dscal
		dscal(&m, rhoi, t, &IntConstOne);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 19: T(:, a) = t;													                 */
		/// T         = [m-by-A]
		/// t         = [m-by-1]

		// Copy m elemments from t into T at offset m * a
		memcpy(T + m * a, t, sizeof(double) * m);

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 20: B(a, 1) = rho;													                 */
		/// B         = [A-by-2]
		/// rho       = [1-by-1]

		B[a] = *rho;

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 21: B(a - 1, 2) = theta;													         */
		/// B         = [A-by-2]
		/// theta     = [1-by-1]

		B[(a-1) + A] = *theta;

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 22: d = (w - theta*d) / rho;													     */
		/// d         = [n-by-1]
		/// rho       = [1-by-1]
		/// theta     = [1-by-1]
		/// w         = [n-by-1]

		for (row = 0; row < n; row++){
			 d[row] = ( w[row] - (*theta) * d[row] ) / (*rho);
		}

		/*—————————————————————————————————————————————————————————————————————————————————————————
		Step 23: beta(:, a) = beta(:, a - 1) + (t'*y)*d;									     */
		/// beta      = [n-by-A]
		/// t         = [m-by-1]
		/// y         = [m-by-1]
		/// d         = [n-by-1]

		// (t'*y) vector-vector product -> ddot
		*tty = ddot(&m, t, &IntConstOne, y, &IntConstOne);

		// (tty*d) = scalar-vector product. (Here treated as matrix-matrix).
		dgemm("N", "N", &n, &p, &p, &one, d, &n, tty, &p, &zero, beta + n * a , &n);

		//  beta(:, a) ←  beta(:, a) + 1 * beta(:, a - 1) -> daxpy
		daxpy(&n, &one, beta + n * (a-1), &IntConstOne, beta + n * a, &IntConstOne);

		//for (row = 0; row < n; row++) {
			//beta[row + n * a] = beta[row + n * a] + beta[row + n * (a - 1)];
			//beta[row + n * a] = beta[row + n * (a - 1)];
		//}

    /* End of component loop.																	 */
	}

	free(t);
	free(T);
	free(Xw);
	free(TTt);
}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for drawing a random integer that lies within range.								 */
int randr(unsigned int min, unsigned int max) {
	return min + rand() / (RAND_MAX / (max - min + 1) + 1);
}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for marking rows in the data as validation or not during CV							 */
void MarkAsVal(bool *IsVal, int NoValObs, size_t m, size_t *valm, size_t *trainm, const double *y) {
	int activeobs, row, index;
	bool alreadyactive;

	/* Initialize all IsVal elements to false */
	for (row = 0; row < m; row++) {
		IsVal[row] = false;
	}

	activeobs = 0;
	while (activeobs < NoValObs) { // NOTE! this needs a check aginst input NoValObs larger than unique(y)!
								   /* Draw random rownumber between 0 and m */
		index = randr(0, m);
		alreadyactive = false;
		/* Mark all the rows with the same response as validation to keep replicates together */
		for (row = 0; row < m; row++) {
			if (y[row] == y[index]) {
				/* If IsVal was 0 it means this is a newly found observation. */
				if (IsVal[row] == false & alreadyactive == false) {
					activeobs += 1;
					alreadyactive = true;
				}
				/* If IsVal is already 1 stop looping and pick a new observation to save time */
				else if (IsVal[row] == true) {
					break;
				}
				IsVal[row] = true;
			}
		}
	}

	/* Calculate the number of rows in Xval/yval and Xtrain/ytrain so they can be allocated */
	*valm = 0;
	for (row = 0; row < m; row++) {
		if (IsVal[row] == true) {
			*valm += 1;
		}
	}
	*trainm = m - (*valm);
}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for populating Xtrain Xval ytrain yval given known IsVal vector					     */
void ExtractXandY(const double *X, const double *y, double *Xtrain, double *Xval, double *ytrain,
	double *yval, size_t m, size_t n, bool *IsVal) {

	int col, row, xte, xve, yve, yte;

	xte = xve = 0;
	yte = yve = 0;

	/* Outer loop over all columns of X */
	for (col = 0; col < n; col++) {
		/* Inner loop over all rows of X */
		for (row = 0; row < m; row++) {

			if (IsVal[row] == true) {
				Xval[xve] = X[row + col * m];
				xve += 1;
				if (col == 0) {
					yval[yve] = y[row];
					yve += 1;
				}
			}
			else {
				Xtrain[xte] = X[row + col * m];
				xte += 1;
				if (col == 0) {
					ytrain[yte] = y[row];
					yte += 1;
				}
			}
		}
	}
}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for calculating validation error using Xval and beta								 */
void Pred(double *Xval, double *beta, double *yval, size_t n, size_t valm, int A, double *yhat, double *RMSEcv, int cvIter, int MaxIters, int ShavingIndex) {

	size_t st_A;
	int col;
	double one = 1.0, zero = 0.0, none = -1.0;
	mwSignedIndex IntConstOne = 1;
	st_A = (size_t)A;

	/* yhat = Xval*beta														                     */
	dgemm("N", "N", &valm, &st_A, &n, &one, Xval, &valm, beta, &n, &zero, yhat, &valm);

	/* Loop the columns of yhat and calulate the error for each component.						 */
	for (col = 0; col < A; col++) {

		/* yhat is here overwritten and becomes epsilon	where epsilon = yhat + (-Yval).	         */
		daxpy(&valm, &none, yval, &IntConstOne, yhat + col * valm, &IntConstOne);

		/* Calculate the Mean Squared Error by taking the dot product of each error column with
		itself and then dividing by the number of observations (valm). At the same time, add it
		to the	existing value of the RMSE vector, thereby performing a summation of the errors
		of each cross-validation iteration at the same time (prevents the need of a 2D matrix).	 */
		if (ShavingIndex == 0) {
			RMSEcv[ShavingIndex + col*(n + 1)] += ddot(&valm, yhat + col * valm, &IntConstOne, yhat + col * valm, &IntConstOne) / valm;
		}
		else {
			RMSEcv[ShavingIndex + col*(n + 2)] += ddot(&valm, yhat + col * valm, &IntConstOne, yhat + col * valm, &IntConstOne) / valm;
		}
		/* If its the last cross-validation iteration, finalize RMSEcv (which right now
		represents sum(MSE,1) [folds-by-components] ) by dividing by the number of cv iterations
		and taking the square-root of the currently squared errors.								 */
		if (cvIter == MaxIters - 1) {
			if (ShavingIndex == 0) {
				RMSEcv[ShavingIndex + col*(n + 1)] = sqrt(RMSEcv[ShavingIndex + col*(n + 1)] / MaxIters);
			}
			else {
				RMSEcv[ShavingIndex + col*(n + 2)] = sqrt(RMSEcv[ShavingIndex + col*(n + 2)] / MaxIters);
			}
		}

	}

}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

/* ————————————————————————————————————————————————————————————————————————————————————————————— */
/* Function for placing a column at the and of a matrix (swapping it with the last column)       */
void SwapCols(double *Xtrain, double *Xval, size_t trainm, size_t valm, size_t n, int ColToShave) {

	mwSignedIndex IntConstOne = 1;

	/* Use dswap to swap location of two column vectors in Xtrain.								 */
	dswap(&trainm, Xtrain + (n - 1) * trainm, &IntConstOne,	Xtrain + ColToShave * trainm, &IntConstOne);

	/* Use dswap to swap location of two column vectors in Xval.								 */
	dswap(&valm, Xval + (n - 1) * valm, &IntConstOne, Xval + ColToShave * valm, &IntConstOne);
}
/* ————————————————————————————————————————————————————————————————————————————————————————————— */

max_stars_repo_name: Petter-s/RegressionFunctions

max_issues_repo_name: Petter-s/RegressionFunctions

max_forks_repo_name: Petter-s/RegressionFunctions


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include "sum.h"

int main(int argc, char *argv[]) {
    int num_bands = 2;
    double t = 1.0;
    double E0 = 6.0;
    double deltaE = 14.0;
    void Efn(double k[3], gsl_vector *energies) {
        double cx = cos(2.0 * M_PI * k[0]);
        double cy = cos(2.0 * M_PI * k[1]);
        double cz = cos(2.0 * M_PI * k[2]);
        double tk = -2.0 * t * (cx + cy + cz);
        int i;
        double E0_band;
        for (i = 0; i < num_bands; i++) {
            E0_band = E0 + ((double)i) * deltaE;
            gsl_vector_set(energies, i, E0_band + tk);
        }
    }
    int na = 8;
    int nb = 8;
    int nc = 8;
    bool use_cache = true;

    gsl_matrix *R = gsl_matrix_calloc(3, 3); // R -> all zeros
    // Overall scale for R doesn't matter.
    gsl_matrix_set(R, 0, 0, 1.0);
    gsl_matrix_set(R, 1, 1, 1.0);
    gsl_matrix_set(R, 2, 2, 1.0);

    double num_electrons = 1.0;
    double E_Fermi = 0.0;
    double energy = SumEnergy(&E_Fermi, Efn, na, nb, nc, num_bands, num_electrons, R, use_cache);
    printf("Got E_Fermi = %f\n", E_Fermi);
    printf("Got energy = %f\n", energy);

    double tol = 1e-8;
    double expected = 6.0;
    if (fabs(energy - expected) > tol) {
        printf("Incorrect energy; got %f, expected %f\n", energy, expected);
        return 1;
    }

    printf("Sum test passed.\n");
    return 0;
}


text:
/**
 *
 * @generated d Tue Jan  7 11:45:26 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#ifdef USE_MKL
#include <mkl_cblas.h>
#include <mkl_lapacke.h>
#else
#include <cblas.h>
#include <lapacke.h>
#endif
#include "plasma.h"
#include <core_blas.h>
#include "auxiliary.h"

/*--------------------------------------------------------------
 * Check the solution
 */

double d_check_solution(int M, int N, int NRHS, double *A, int LDA,
                      double *B,  double *X, int LDB)
{
/*     int info_solution; */
    double Rnorm = -1.00;
    double zone  =  1.0;
    double mzone = -1.0;
    double *work = (double *)malloc(max(M, N)* sizeof(double));

    LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, N,    A, LDA, work);
    LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work);
    LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, (zone), A, LDA, X, LDB, (mzone), B, LDB);

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work);

    free(work);

    return Rnorm;
}


text:
//
// Copyright (c) Microsoft. All rights reserved.
// This code is licensed under the MIT License (MIT).
// THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF
// ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
// IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
// PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
//
// Developed by Minigraph
//
// Author:  James Stanard
//

#pragma once

#pragma warning(disable:4201) // nonstandard extension used : nameless struct/union
#pragma warning(disable:4238) // nonstandard extension used : class rvalue used as lvalue
#pragma warning(disable:4324) // structure was padded due to __declspec(align())

#ifndef WIN32_LEAN_AND_MEAN
    #define WIN32_LEAN_AND_MEAN
#endif
#ifndef NOMINMAX
    #define NOMINMAX
#endif
#include <windows.h>

#include <d3d12.h>

#pragma comment(lib, "d3d12.lib")
#pragma comment(lib, "dxgi.lib")

#define MY_IID_PPV_ARGS IID_PPV_ARGS
#define D3D12_GPU_VIRTUAL_ADDRESS_NULL      ((D3D12_GPU_VIRTUAL_ADDRESS)0)
#define D3D12_GPU_VIRTUAL_ADDRESS_UNKNOWN   ((D3D12_GPU_VIRTUAL_ADDRESS)-1)

#pragma region
#include <d3dx12.h>
#include <d3d12shader.h>
#include <d3dcompiler.h>
#pragma endregion

#pragma region

#pragma comment (lib, "d3dcompiler.lib")

#pragma endregion

#include <cstdint>
#include <cstdio>
#include <cstdarg>
#include <vector>
#include <memory>
#include <string>
#include <exception>
#include <numeric>
#include <wrl.h>
#include <ppltasks.h>
#include <gsl\gsl>

#include "Utility.h"
#include "VectorMath.h"
#include "EngineTuning.h"
#include "EngineProfiling.h"

template< typename T >
using NotNull = gsl::not_null< T >;

#define SHADER_ARGS(shader) (shader), \
								 sizeof(shader)

template<typename T, std::size_t N>
struct Array :public std::array<T, N>
{
    static inline std::function<T(T, T)> multiply = [](T priview, T current) {
        return priview * current;
    };
public:
    T product()
    {
        return std::accumulate(begin(), end(), 1, multiply);
    }

};

using U32 = uint32_t;
using U32x1 = U32;
using U32x2 = Array<U32, 2>;
using U32x3 = Array<U32, 3>;
using F32 = float;
using F32x2 = Array<F32, 2>;
using F32x3 = Array<F32, 3>;

max_stars_repo_name: sienaiwun/DirectX-Graphics-Samples

max_issues_repo_name: sienaiwun/DirectX-Graphics-Samples

max_forks_repo_name: sienaiwun/DirectX-Graphics-Samples


text:
/* bcls.c
   $Revision: 282 $ $Date: 2006-12-17 17:38:00 -0800 (Sun, 17 Dec 2006) $

   ---------------------------------------------------------------------
   This file is part of BCLS (Bound-Constrained Least Squares).

   Copyright (C) 2006 Michael P. Friedlander, Department of Computer
   Science, University of British Columbia, Canada. All rights
   reserved. E-mail: <mpf@cs.ubc.ca>.

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

   BCLS 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 Lesser General
   Public License for more details.

   You should have received a copy of the GNU Lesser General Public
   License along with BCLS; if not, write to the Free Software
   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
   USA
   ---------------------------------------------------------------------
*/
/*!
  \file
  BCLS user-callable library routines.
*/

#include <cblas.h>
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
#include <assert.h>
#include <setjmp.h>

#include "bcls.h"
#include "bclib.h"
#include "bcsolver.h"
#include "bcversion.h"

/*!

  \brief Malloc wrapper.  Private to this file.

  Pointer to the newly allocated memory,
  or NULL if malloc returned an eror.

  \param[in]  len   Number of elements needed
  \param[in]  size  Memory needed for each element.
  \param[in]  who   Short character description of the memory needed.

  \return Pointer to the newly allocated memory.

*/
static void *
xmalloc (int len, size_t size, char * who)
{
    register void *value = malloc(len * size);
    if (value == NULL)
	fprintf( stderr, "Not enough memory to allocate %s.\n", who );
    return value;
}

/*!

  \brief Create a BCLS problem instance.

  This routine create a BCLS problem.  It *must* be called before any
  other BCLS routine.  It will create and initialize a new BCLS
  problem with enough space to accomodate the specified problem.  The
  problem is then initialized using bcls_init_prob.

  In order to re-initialize the BCLS problem instance, but not
  deallocate memory that's already sufficient for an equal or smaller
  size problem, use bcls_init_prob.

  \param[in] mmax  Maximum number of rows    in any A.
  \param[in] nmax  Maximum number of columns in any A.

  \return Return a pointer to a new BCLS problem.  If the return is
  NULL, then the initialization failed.

*/
BCLS *
bcls_create_prob( int mmax, int nmax )
{
    int mnmax = imax( nmax, mmax );

    assert( mmax > 0 && nmax > 0 );

    // Allocate the problem.
    BCLS *ls = (BCLS *)xmalloc(1, sizeof(BCLS), "ls" );

    // Check if the problem has been successfully allocated.
    if (ls == NULL) {
	fprintf( stderr, "XXX Could not allocate a BCLS problem.\n");
	return ls;
    }

    // -----------------------------------------------------------------
    // Allocate workspace vectors large enough to accomdate the problem.
    // -----------------------------------------------------------------

    // Record the maximum alloted workspace.
    ls->mmax = mmax;
    ls->nmax = nmax;

    // Residual: r(m+n).
    ls->r = (double *)xmalloc( mmax+nmax, sizeof(double), "r" );
    if (ls->r == NULL) goto error;

    // Gradient: g(n).
    ls->g = (double *)xmalloc( nmax, sizeof(double), "g" );
    if (ls->g == NULL) goto error;

    // Search direction, full space: dx(n).
    ls->dx = (double *)xmalloc( nmax, sizeof(double), "dx" );
    if (ls->dx == NULL) goto error;

    // Search direction, subspace: dxFree(n).
    ls->dxFree = (double *)xmalloc( nmax, sizeof(double), "dxFree" );
    if (ls->dxFree == NULL) goto error;

    // Step to each breakpoint: aBreak(n).
    ls->aBreak = (double *)xmalloc( nmax, sizeof(double), "aBreak" );
    if (ls->aBreak == NULL) goto error;

    // Indices of each breakpoint: iBreak(n).
    ls->iBreak = (int *)xmalloc( nmax, sizeof(int), "iBreak" );
    if (ls->iBreak == NULL) goto error;

    // Variable indices: ix(n).
    ls->ix = (int *)xmalloc( nmax, sizeof(int), "ix" );
    if (ls->ix == NULL) goto error;

    // Workspace: wrk_v( max(n, m) ).
    ls->wrk_u = (double *)xmalloc( mnmax, sizeof(double), "wrk_u" );
    if (ls->wrk_u == NULL) goto error;

    // Workspace: wrk_v( max(n, m) ).
    ls->wrk_v = (double *)xmalloc( mnmax, sizeof(double), "wrk_v" );
    if (ls->wrk_v == NULL) goto error;

    // Workspace: wrk_w( max(n, m) ).
    ls->wrk_w = (double *)xmalloc( mnmax, sizeof(double), "wrk_w" );
    if (ls->wrk_w == NULL) goto error;

    // -----------------------------------------------------------------
    // Initialize this new problem instance.
    // -----------------------------------------------------------------
    bcls_init_prob( ls );

    // -----------------------------------------------------------------
    // Exits.
    // -----------------------------------------------------------------

    // Successfull exit.
    return ls;

 error:
    // Unsuccessful exit.
    bcls_free_prob( ls );
    return ls;
}

/*!

  \brief Initialize a BCLS problem.

  Initialize a BCLS problem.  You can call this routine to "reset" a
  BCLS problem, i.e., reset all parameter values.  Note that it is
  automatically called by bcls_create_prob.

  \param[in,out] ls  BCLS problem context.

*/
void
bcls_init_prob( BCLS *ls )
{
    // Some quick error checking.
    assert( ls->mmax > 0 && ls->nmax > 0 );

    // Check if the problem has been successfully allocated.
    if (ls == NULL) {
	fprintf( stderr, "XXX The BCLS problem is NULL.\n");
	return;
    }

    // Initialize the problem structure.
    ls->print_info    =  NULL;
    ls->print_hook    =  NULL;
    ls->fault_info    =  NULL;
    ls->fault_hook    =  NULL;
    ls->Aprod         =  NULL;
    ls->Usolve        =  NULL;
    ls->CallBack      =  NULL;
    ls->UsrWrk        =  NULL;
    ls->anorm         =  NULL;
    ls->print_level   =  1;
    ls->proj_search   =  BCLS_PROJ_SEARCH_EXACT;
    ls->newton_step   =  BCLS_NEWTON_STEP_LSQR;
    ls->minor_file    =  NULL;
    ls->itnMaj        =  0;
    ls->itnMajLim     =  5  * (ls->nmax);
    ls->itnMin        =  0;
    ls->itnMinLim     =  10 * (ls->nmax);
    ls->nAprodT       =  0;
    ls->nAprodF       =  0;
    ls->nAprod1       =  0;
    ls->nUsolve       =  0;
    ls->m             =  0;
    ls->n             =  0;
    ls->unconstrained =  0;
    ls->damp          =  0.0;
    ls->damp_min      =  1.0e-4;
    ls->exit          =  BCLS_EXIT_UNDEF;
    ls->soln_rNorm    = -1.0;
    ls->soln_dInf     =  0.0;
    ls->soln_jInf     = -1;
    ls->soln_stat     =  BCLS_SOLN_UNDEF;
    ls->optTol        =  1.0e-6;
    ls->conlim        =  1.0 / ( 10.0 * sqrt( DBL_EPSILON ) );
    ls->mu            =  1.0e-2;
    ls->backtrack     =  1.0e-1;
    ls->backtrack_limit = 10;

    // Initialize timers.
    bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL] ),  BCLS_TIMER_INIT );
    ls->stopwatch[BCLS_TIMER_TOTAL].name = "Total time";

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_APROD] ),  BCLS_TIMER_INIT );
    ls->stopwatch[BCLS_TIMER_APROD].name = "Total time for Aprod";

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_USOLVE] ),  BCLS_TIMER_INIT );
    ls->stopwatch[BCLS_TIMER_USOLVE].name = "Total time for Usolve";

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_LSQR] ),  BCLS_TIMER_INIT );
    ls->stopwatch[BCLS_TIMER_LSQR].name = "Total time for LSQR";

    // Initialize constants.
    ls->eps         = DBL_EPSILON;
    ls->eps2        = pow( DBL_EPSILON, 0.50 );
    ls->eps3        = pow( DBL_EPSILON, 0.75 );
    ls->epsx        = ls->eps3;
    ls->epsfixed    = DBL_EPSILON;
    ls->BigNum      = BCLS_INFINITY;

    return;
}

/*!

  \brief Free and reinitialize all resources in an existing BCLS problem.

  \return
  - 0: no errors
  - 1: the BCLS problem is NULL (ie, does not exist).

*/
int
bcls_free_prob( BCLS *ls )
{
    // Report an error if it's already NULL.
    if (ls == NULL) return 1;

    // Deallocate the workspace.
    assert( ls->r       != NULL ); free( ls->r       );
    assert( ls->g       != NULL ); free( ls->g       );
    assert( ls->dx      != NULL ); free( ls->dx      );
    assert( ls->dxFree  != NULL ); free( ls->dxFree  );
    assert( ls->aBreak  != NULL ); free( ls->aBreak  );
    assert( ls->iBreak  != NULL ); free( ls->iBreak  );
    assert( ls->ix      != NULL ); free( ls->ix      );
    assert( ls->wrk_u   != NULL ); free( ls->wrk_u   );
    assert( ls->wrk_v   != NULL ); free( ls->wrk_v   );
    assert( ls->wrk_w   != NULL ); free( ls->wrk_w   );

    // Deallocate the BCLS problem.
    free( ls );

    return 0;
}

/*!

  \brief Install a print-hook routine.

  This routine installs a user-defined print-hook routine.

  The parameter info is a transit pointer passed to the hook routine.

  The parameter hook is an entry point to the user-defined print-hook
  routine. This routine is called by the routine "print" every time an
  informative message should be output. The routine "print" passes to
  the hook routine the transit pointer info and the character string
  msg, which contains the message. If the hook routine returns zero,
  the routine print prints the message in an usual way. Otherwise, if
  the hook routine returns non-zero, the message is not printed.

  In order to uninstall the hook routine the parameter hook should be
  specified as NULL (in this case the parameter info is ignored).

  \param[in,out] ls    BCLS problem context.
  \param[in]     info  Transit pointer passed to the hook routine.
  \param[in]     hook  Pointer to the print-hook routine.

*/
void
bcls_set_print_hook( BCLS *ls,
		     void *info,
		     int (*hook)(void *info, char *msg))
{
    ls->print_info = info;
    ls->print_hook = hook;
    return;
}

/*!

  \brief Install a fault-hook routine.

  This routine installs a user-defined fault-hook routine.

  The parameter info is a transit pointer passed to the hook routine.

  The parameter "hook" is an entry point to the user-defined
  fault-hook routine. This routine is called by the routine
  "bcls_fault" every time an error message should be output. The
  routine "bcls_fault" passes to the hook routine the transit pointer
  info and the character string msg, which contains the message. If
  the hook routine returns zero, the routine print prints the message
  in an usual way. Otherwise, if the hook routine returns non-zero,
  the message is not printed.

  In order to uninstall the hook routine the parameter hook should be
  specified as NULL (in this case the parameter info is ignored).

  \param[in,out] ls    BCLS problem context.
  \param[in]     info  Transit pointer passed to the hook routine.
  \param[in]     hook  Pointer to the fault-hook routine.

*/
void
bcls_set_fault_hook( BCLS *ls,
		     void *info,
		     int (*hook)(void *info, char *msg))
{
    ls->fault_info = info;
    ls->fault_hook = hook;
    return;
}

/*!

  \brief Give BCLS access to set of column weights.

  Define a set of column weights.  BCLS will use these to scale the
  steepest-descent steps.

  \param[in,out] ls     BCLS problem context.
  \param[in]     anorm  Array of columns norms of A.

*/
void
bcls_set_anorm( BCLS *ls,
                double anorm[] )
{
    ls->anorm = anorm;
    return;
}

/*!

  \brief Install a user-defined preconditioning routine.

  Replace each subproblem
  \f[
  \def\minimize{\displaystyle\mathop{\hbox{minimize}}}
  \minimize_x \| Ax - b \|
  \f]
  with
  \f[
  \def\minimize{\displaystyle\mathop{\hbox{minimize}}}
  \minimize_y \| AU^{-1}y - b \|
  \quad\mbox{with}\quad
  Ux = y.
  \f]

  \param[in,out]  ls     BCLS problem context.
  \param[in]      Usolve Pointer to the user's preconditioning routine.

*/
void
bcls_set_usolve( BCLS *ls,
                 int (*Usolve)( int mode, int m, int n, int nix,
                                int ix[], double v[], double w[],
                                void *UsrWrk ) )
{
    assert( Usolve != NULL );
    ls->Usolve = Usolve;
    return;
}

/*!

  \brief Compute the column norms of A.

  Compute the columns norms of A.  Each column of A is generated by
  multiplying A by a unit vector.  The two-norm squared of each
  column is stored in aprod.

  This routine must be called *after* bcls_create_prob.  (As with any
  other BCLS routine.)

  \param[in]      ls      BCLS problem context.
  \param[in]      m       Number of rows in A.
  \param[in]      n       Number of columns in A.
  \param[in]      Aprod   User's matrix-product routine.
  \param[in,out]  UsrWrk  Pointer to user's workspace (could be NULL).
  \param[out]     anorm   Vector of column norms of A.

  \return Returns any error codes returned by the user's aprod
  routine.

*/
int
bcls_compute_anorm( BCLS *ls, int n, int m,
                    int (*Aprod)
                    ( int mode, int m, int n, int nix,
                      int ix[], double x[], double y[], void *UsrWrk ),
                    void *UsrWrk,
                    double anorm[] )
{
    // Make sure that the problem instance has been created.
    assert( ls    != NULL );
    assert( anorm != NULL );

    int j;
    int err;
    const int mode = 1;         // Always:  aj = A*ej.
    const int nix  = 1;         // Only one column is ever needed.
    int    *ix = ls->ix;
    double *e  = ls->wrk_u;
    double *aj = ls->wrk_v;

    bcls_dload( n, 1.0, e, 1 );

    for (j = 0; j < n; j++) {

        // Multiply A times the jth unit vector.
        ix[0] = j;

        // aj <- A * ej.
        err = Aprod( mode, m, n, nix, ix, e, aj, UsrWrk );

        // Exit if Aprod returned an error.
        if (err) break;

        // Compute the norm of aj and store it in anorm.
        anorm[j] = cblas_dnrm2( m, aj, 1 );

        // Make sure that the column norm is too small.
        anorm[j] = fmax( BCLS_MIN_COLUMN_NORM, anorm[j] );
    }

    return err;
}

/*!

  \brief Load the problem into the BCLS data structure.

  The routiens loads the complete problem description into a BCLS
  problem instance.  No work is really being done: only the various
  structure elements are being set to point to the right places.

  \param[in,out]  ls      BCLS problem context.
  \param[in]      m       Number of rows in A.  Note: m <= mmax.
  \param[in]      n       Number of columns in A.  Note: n <= nmax.
  \param[in]      Aprod   Hook to user's matrix-vector product routine.
  \param[in,out]  UsrWrk  Transit pointer passed directly to Aprod/Usolve.
  \param[in]      damp    Regularization parameter.
  \param[in,out]  x       The starting point.
  \param[in]      b       The RHS vector.
  \param[in]      c       Defines a linear term (set to NULL if one doesn't exist).
  \param[in]      bl      Lower bounds on x.
  \param[in]      bu      Upper bounds on x.

*/
void
bcls_set_problem_data( BCLS *ls, int m, int n,
		       int (*Aprod)( int mode, int m, int n, int nix,
				     int ix[], double x[], double y[],
				     void *UsrWrk ),
		       void *UsrWrk, double damp,
		       double x[], double b[], double c[],
		       double bl[], double bu[] )
{
    assert( m <= ls->mmax );
    assert( n <= ls->nmax );

    assert( Aprod  != NULL ); ls->Aprod  = Aprod;
    assert( m      >= 0    ); ls->m      = m;
    assert( n      >= 0    ); ls->n      = n;
                              ls->UsrWrk = UsrWrk;
    assert( damp   >= 0.0  ); ls->damp   = damp;
    assert( x      != NULL ); ls->x      = x;
    assert( b      != NULL ); ls->b      = b;
                              ls->c      = c;
    assert( bl     != NULL ); ls->bl     = bl;
    assert( bu     != NULL ); ls->bu     = bu;

    return;
}

/*!

  \brief Return an exit message.

  \param[in] flag  The code of the error.

  \return Error messsage.

*/
char *
bcls_exit_msg( int flag )
{
    char *msg;

    if      (flag==BCLS_EXIT_CNVGD)  msg="Optimal solution found";
    else if (flag==BCLS_EXIT_MAJOR)  msg="Too many major iterations";
    else if (flag==BCLS_EXIT_MINOR)  msg="Too many minor iterations";
    else if (flag==BCLS_EXIT_UNBND)  msg="Found direction of infinite descent";
    else if (flag==BCLS_EXIT_INFEA)  msg="Bounds are inconsistent";
    else if (flag==BCLS_EXIT_APROD)  msg="Aprod requested immediate exit";
    else if (flag==BCLS_EXIT_USOLVE) msg="Usolve requested immediate exit";
    else if (flag==BCLS_EXIT_CALLBK) msg="CallBack requested immediate exit";
    else                             msg="Undefined exit";

    return msg;
}

/*!

  \brief Solve the current problem instance.

  \return
   -  #BCLS_EXIT_CNVGD (0) - Successful exit.  Found an optimal solution.
   -  #BCLS_EXIT_MAJOR     - Too many major iterations.
   -  #BCLS_EXIT_MINOR     - Too many inner iterations.
   -  #BCLS_EXIT_UNBND     - Found direction of infinite descent.
   -  #BCLS_EXIT_INFEA     - Bounds are inconsistent.
   -  #BCLS_EXIT_APROD     - Aprod  requested immediate exit.
   -  #BCLS_EXIT_USOLVE    - Usolve requested immediate exit.

*/
int
bcls_solve_prob( BCLS *ls )
{
    int err;
    int jpInf;
    const int preconditioning = ls->Usolve != NULL;
    double pInf, timeTot;
    double bNorm;  // Norm of the RHS (or 1.0, which ever is bigger).
    BCLS_timer watch;

    // Print the banner.
    //PRINT1(" ----------------------------------------------------\n");
    //PRINT1(" BCLS -- Bound-Constrained Least Squares, Version %s\n",
    //      bcls_version_info() );
    //PRINT1(" Compiled on %s\n",  bcls_compilation_info() );
    //PRINT1(" ----------------------------------------------------\n");

    // -----------------------------------------------------------------
    // Print parameters and diagnostic information.
    // -----------------------------------------------------------------
    bcls_print_params( ls );

    // Examine the bounds and print a summary about them.
    // Also check if the problem is feasible!
    err = bcls_examine_bnds( ls, ls->n, ls->bl, ls->bu );
    if (err) {
	ls->exit = err;
        goto direct_exit;
    }

    // Set the return for longjmp.  First call to setjmp returns 0.
    err = setjmp(ls->jmp_env);
    if (err) {
        ls->exit = err;
        goto direct_exit;
    }

    // Compute and print some statistics about the column scales.
    bcls_examine_column_scales( ls, ls->anorm );

    // Print a log header.
    //    PRINT1("\n %5s  %5s    %9s  %9s  %11s  %9s %7s %7s %7s\n",
    //	   "Major","Minor","Residual","xNorm","Optimal",
    //	   ls->newton_step == BCLS_NEWTON_STEP_LSQR
    //	   ? "LSQR" : "CGLS",
    //	   "nSteps", "Free", "OptFac");

    // -----------------------------------------------------------------
    // Call the BCLS solver.
    // -----------------------------------------------------------------
    bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL]), BCLS_TIMER_START );

    bcls_solver( ls, ls->m, ls->n, &bNorm,
		 ls->x, ls->b, ls->c, ls->bl, ls->bu, ls->r, ls->g,
		 ls->dx, ls->dxFree, ls->ix,
		 ls->aBreak, ls->iBreak, ls->anorm );

    bcls_timer( &(ls->stopwatch[BCLS_TIMER_TOTAL]), BCLS_TIMER_STOP );

    // =================================================================
    // Exits.
    // =================================================================
 direct_exit:
    // -----------------------------------------------------------------
    // Check feasibility and print a solution summary (iterations, etc.)
    // -----------------------------------------------------------------
    bcls_primal_inf( ls->n, ls->x, ls->bl, ls->bu, &pInf, &jpInf );

    //PRINT1("\n");
    //    PRINT1(" BCLS Exit %4d -- %s\n",ls->exit, bcls_exit_msg(ls->exit));
    //PRINT1("\n");
    //PRINT1(" No. of iterations      %8d %7s", ls->itnMin, "");
    //PRINT1(" Objective value       %17.9e\n", ls->soln_obj  );
    //PRINT1(" No. of major iterations%8d %7s", ls->itnMaj, "");
    //PRINT1(" Optimality           (%7d) %8.1e\n",
    //	   ls->soln_jInf,ls->soln_dInf);
    //PRINT1(" No. of calls to Aprod  %8d", ls->nAprodT);

    //if (ls->proj_search == BCLS_PROJ_SEARCH_EXACT)
    //   PRINT1(" (%5d)", ls->nAprod1 );
    //else
    //  PRINT1("  %5s ", "");
    // PRINT1(" Norm of RHS            %16.1e\n", bNorm);
    // if (pInf > ls->eps)
    //  PRINT1(" Feasibility          (%7d) %7.1e\n", jpInf, pInf);
    //if (preconditioning)
    //  PRINT1(" No. of calls to Usolve %8d\n", ls->nUsolve);
    //PRINT1("\n");

    // -----------------------------------------------------------------
    // Print timing statistics.
    // -----------------------------------------------------------------
    //timeTot = fmax(ls->eps,ls->stopwatch[BCLS_TIMER_TOTAL].total);

    //watch   = ls->stopwatch[BCLS_TIMER_TOTAL];
    //timeTot = fmax(1e-6, watch.total); // Safeguard against timeTot = 0.

    //PRINT1(" %-25s %5.1f (%4.2f) secs\n",
    //	   watch.name, watch.total, 1.0 );

    //watch = ls->stopwatch[BCLS_TIMER_LSQR];
    //PRINT1(" %-25s %5.1f (%4.2f) secs\n",
    //	   watch.name, watch.total, watch.total / timeTot );

    //    watch = ls->stopwatch[BCLS_TIMER_APROD];
    //PRINT1(" %-25s %5.1f (%4.2f) secs\n",
    //	   watch.name, watch.total, watch.total / timeTot );

    // Only print time for Usolve is user provided this routine.
    //if (preconditioning) {
      //  watch = ls->stopwatch[BCLS_TIMER_USOLVE];
      //PRINT1(" %-25s %5.1f (%4.2f) secs\n",
      //       watch.name, watch.total, watch.total / timeTot );
    //}

    return (ls->exit + err);
}

max_stars_repo_name: echristakopoulou/glslim

max_issues_repo_name: echristakopoulou/glslim

max_forks_repo_name: echristakopoulou/glslim


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT License.

#pragma once

#include "activation_info.h"
#include "attention_wrapper.h"
#include "core/providers/cpu/rnn/rnn_helpers.h"

#include "core/common/common.h"
#include "core/common/logging/logging.h"
#include "core/framework/allocator.h"

#include <gsl/span>

namespace onnxruntime {
namespace contrib {

namespace deepcpu = ::onnxruntime::rnn::detail::deepcpu;
namespace logging = ::onnxruntime::logging;

using ::onnxruntime::AllocatorPtr;
using ::onnxruntime::IAllocatorUniquePtr;
using ::onnxruntime::contrib::detail::ActivationInfo;
using ::onnxruntime::rnn::detail::ActivationFuncs;
using ::onnxruntime::rnn::detail::Direction;

namespace rnn {
namespace detail {

// copying the peephole values into UniDirectionalAttnLstm seems unnecessary. don't do that until proven necessary
#define LSTM_NO_PEEPHOLE_COPY

template <typename T>
class UniDirectionalAttnLstm {
 public:
  UniDirectionalAttnLstm(AllocatorPtr allocator,
                         const logging::Logger& logger,
                         const int seq_length,
                         const int batch_size,
                         const int input_size,
                         const int hidden_size,
                         Direction direction,
                         const bool input_forget,
                         AttentionWrapper<T>& attention_wrapper,
                         const gsl::span<const T>& bias,
                         const gsl::span<const T>& peephole_weights,
                         const gsl::span<const T>& initial_hidden_state,
                         const gsl::span<const T>& initial_cell_state,
                         const ActivationFuncs::Entry& activation_func_f,
                         const ActivationFuncs::Entry& activation_func_g,
                         const ActivationFuncs::Entry& activation_func_h,
                         const float clip,
                         onnxruntime::concurrency::ThreadPool* ttp);

  void Compute(const gsl::span<const T>& inputs,
               const gsl::span<const int>& sequence_lengths,
               const int num_directions,
               const gsl::span<const T>& input_weights,
               const gsl::span<const T>& recurrent_weights,
               gsl::span<T>& outputs,
               gsl::span<T>& final_hidden_state,
               gsl::span<T>& final_cell_state);

  ~UniDirectionalAttnLstm() = default;

  AttentionWrapper<T>& GetAttentionWrapper() {
    return attention_wrapper_;
  }

 private:
  using span_T_const_iter = typename gsl::span<T>::const_iterator;
  using span_T_iter = typename gsl::span<T>::iterator;

  void SetNumThreads();

  void GateComputations(span_T_iter& out, span_T_iter& out_end,
                        span_T_iter& C_prev, span_T_iter& C_prev_end,  // Ct-1 value not 'ct'. using 'C' for clarity
                        span_T_iter& C_prev_clipped, span_T_iter& C_prev_clipped_end,
                        span_T_iter& batched_output, span_T_iter& batched_output_end,
                        const gsl::span<const int>& seq_lengths,
                        const int min_sequence_length,
                        const int step,
                        const int row,
                        const int local_fused_hidden_rows,
                        bool output_sequence);

  void AllocateBuffers();

  void InitializeBuffers(const gsl::span<const T>& initial_hidden_state,
                         const gsl::span<const T>& initial_cell_state);

  void LoadPeepholeWeights(const gsl::span<const T>& peephole_weights);
  void LoadBias(const gsl::span<const T>& WbRb_values);

  AllocatorPtr allocator_;
  const logging::Logger& logger_;

  int seq_length_;
  int batch_size_;
  int input_size_;
  int hidden_size_;
  int attention_size_;
  int attention_context_size_;

  Direction direction_;
  bool input_forget_;
  float clip_;

  bool batch_parallel_;

  bool use_bias_;
  bool use_peepholes_;

  int hidden_num_threads_ = -1;

  IAllocatorUniquePtr<T> output_iofc_ptr_;
  IAllocatorUniquePtr<T> hidden0_ptr_, batched_hidden0_ptr_;
  gsl::span<T> output_iofc_;
  gsl::span<T> hidden0_, batched_hidden0_;

  IAllocatorUniquePtr<T> internal_memory_prev_ptr_, batched_internal_memory_prev_ptr_;
  IAllocatorUniquePtr<T> internal_memory_cur_ptr_, batched_internal_memory_cur_ptr_;
  IAllocatorUniquePtr<T> batched_internal_memory_clipped_ptr_;
  gsl::span<T> internal_memory_prev_, batched_internal_memory_prev_;
  gsl::span<T> internal_memory_cur_, batched_internal_memory_cur_;
  gsl::span<T> batched_internal_memory_clipped_;

  IAllocatorUniquePtr<T> bias_WRi_ptr_, bias_WRf_ptr_, bias_WRo_ptr_, bias_WRc_ptr_;
  IAllocatorUniquePtr<T> batched_bias_WRi_ptr_, batched_bias_WRf_ptr_, batched_bias_WRo_ptr_, batched_bias_WRc_ptr_;
  IAllocatorUniquePtr<T> peephole_i_ptr_, peephole_f_ptr_, peephole_o_ptr_;
  IAllocatorUniquePtr<T> inputs_reverse_ptr_, outputs_reverse_ptr_;
  gsl::span<T> bias_WRi_, bias_WRf_, bias_WRo_, bias_WRc_;
  gsl::span<T> batched_bias_WRi_, batched_bias_WRf_, batched_bias_WRo_, *batched_bias_WRc_;
  gsl::span<T> inputs_reverse_, outputs_reverse_;

#if defined(LSTM_NO_PEEPHOLE_COPY)
  gsl::span<const T> peephole_i_, peephole_f_, peephole_o_;
#else
  gsl::span<T> peephole_i_, peephole_f_, peephole_o_;
#endif

  IAllocatorUniquePtr<int> sequence_lengths_ptr_;
  gsl::span<int> sequence_lengths_;

  deepcpu::ClipWithBiasFuncPtr clip_with_bias_ptr_;

  ActivationInfo<deepcpu::ActivationFuncPtr> activation_f_;
  ActivationInfo<deepcpu::ActivationFuncPtr> activation_g_;
  ActivationInfo<deepcpu::LstmMergeGatesFuncPtr> activation_h_;

  AttentionWrapper<T>& attention_wrapper_;

  onnxruntime::concurrency::ThreadPool* ttp_;
};

}  // namespace detail
}  // namespace rnn
}  // namespace contrib
}  // namespace onnxruntime

max_stars_repo_name: hqucms/onnxruntime

max_issues_repo_name: hqucms/onnxruntime

max_forks_repo_name: hqucms/onnxruntime