text:
#ifndef BASIS_PURSUIT_E_H
#define BASIS_PURSUIT_E_H

#include "Matrix.h"
#include "LADMM.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem of the form f(x)+g(x) under the constraint Mx=c by IPALM_APPROX;
// where f(x)=0
//and g(x)=\frac{lambda2}{2}\|x\|_2+lambda1\|x\|_1.

template<typename L, typename D>
class Basis_pursuit_e: public LADMM<L, D>
{
private:

D lambda1;

D lambda2;

Matrix<L,D> my_M;

Matrix<L,D> my_A;

D val_lambda_f;

protected:

public:

  Basis_pursuit_e(const char* Matrix_file,D val_lambda1, D val_lambda2)
  :LADMM<L,D>(), my_M(Matrix_file)
  {
    lambda1=val_lambda1;
    lambda2=val_lambda2;
    this->lambda_f= 0;
  }

  L get_n(){return my_M.nfeatures;}
  L get_m(){return my_M.nsamples;}

  inline D value_of_g_j(D x, L i){
        return lambda2*x*x/2+lambda1*fabs(x);
      }

  inline D value_of_f_j(D x, L i){return 0;}

  inline D value_of_h_j(D x, L i){
	  return 0;
  }

  inline D gradient_of_f_j(D x, L i){return 0;}

  inline D prox_of_h_j(D x1,D x2, L i){
  	return my_M.b[i];
  }

  inline D prox_of_g_j(D x1,D x2, L i){
        D new_x;
        if(x1*x2> lambda1)
         new_x=(x1*x2- lambda1)/(x2+lambda2);
        else if(x1*x2< -lambda1)
         new_x=(x1*x2+ lambda1)/(x2+ lambda2);
        else
         new_x=0;
        return new_x;
      }

  inline void set_matrix_M(){
    this->data_M=my_M;
  }

  inline void set_matrix_A(){
    this->data_A.nsamples=0;
    this->data_A.nfeatures=this->data_M.nfeatures;
    this->data_A.nnz=0;
    this->data_A.ptr.resize(1,0);
    this->data_A.ptr_t.resize(this->data_M.nfeatures+1,0);
  }

  void LADMM_solver(D beta_0,  D rho,vector<D> & x0,vector<D> & y0, vector<D> & lambda0, L max_nb_outer, L p_N_1, string filename1, D time){
    this->ADMM_solve_with_Linear(beta_0, rho,x0,y0, lambda0, max_nb_outer,  p_N_1, filename1, time);
}

};

#endif


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

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

/* Of the two methods provided below, I think the Polar method is more
 * efficient, but only when you are actually producing two random
 * deviates.  We don't produce two, because then we'd have to save one
 * in a static variable for the next call, and that would screws up
 * re-entrant or threaded code, so we only produce one.  This makes
 * the Ratio method suddenly more appealing.
 *
 * [Added by Charles Karney] We use Leva's implementation of the Ratio
 * method which avoids calling log() nearly all the time and makes the
 * Ratio method faster than the Polar method (when it produces just one
 * result per call).  Timing per call (gcc -O2 on 866MHz Pentium,
 * average over 10^8 calls)
 *
 *   Polar method: 660 ns
 *   Ratio method: 368 ns
 *
 */

/* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */

double
gsl_ran_gaussian (const gsl_rng * r, const double sigma)
{
  double x, y, r2;

  do
    {
      /* choose x,y in uniform square (-1,-1) to (+1,+1) */
      x = -1 + 2 * gsl_rng_uniform_pos (r);
      y = -1 + 2 * gsl_rng_uniform_pos (r);

      /* see if it is in the unit circle */
      r2 = x * x + y * y;
    }
  while (r2 > 1.0 || r2 == 0);

  /* Box-Muller transform */
  return sigma * y * sqrt (-2.0 * log (r2) / r2);
}

/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130.
 * K+M, ACM Trans Math Software 3 (1977) 257-260.
 *
 * [Added by Charles Karney] This is an implementation of Leva's
 * modifications to the original K+M method; see:
 * J. L. Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */

double
gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma)
{
  double u, v, x, y, Q;
  const double s = 0.449871;    /* Constants from Leva */
  const double t = -0.386595;
  const double a = 0.19600;
  const double b = 0.25472;
  const double r1 = 0.27597;
  const double r2 = 0.27846;

  do                            /* This loop is executed 1.369 times on average  */
    {
      /* Generate a point P = (u, v) uniform in a rectangle enclosing
         the K+M region v^2 <= - 4 u^2 log(u). */

      /* u in (0, 1] to avoid singularity at u = 0 */
      u = 1 - gsl_rng_uniform (r);

      /* v is in the asymmetric interval [-0.5, 0.5).  However v = -0.5
         is rejected in the last part of the while clause.  The
         resulting normal deviate is strictly symmetric about 0
         (provided that v is symmetric once v = -0.5 is excluded). */
      v = gsl_rng_uniform (r) - 0.5;

      /* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
         much (for efficiency). */
      v *= 1.7156;

      /* Compute Leva's quadratic form Q */
      x = u - s;
      y = fabs (v) - t;
      Q = x * x + y * (a * y - b * x);

      /* Accept P if Q < r1 (Leva) */
      /* Reject P if Q > r2 (Leva) */
      /* Accept if v^2 <= -4 u^2 log(u) (K+M) */
      /* This final test is executed 0.012 times on average. */
    }
  while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));

  return sigma * (v / u);       /* Return slope */
}

double
gsl_ran_gaussian_pdf (const double x, const double sigma)
{
  double u = x / fabs (sigma);
  double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
  return p;
}

double
gsl_ran_ugaussian (const gsl_rng * r)
{
  return gsl_ran_gaussian (r, 1.0);
}

double
gsl_ran_ugaussian_ratio_method (const gsl_rng * r)
{
  return gsl_ran_gaussian_ratio_method (r, 1.0);
}

double
gsl_ran_ugaussian_pdf (const double x)
{
  return gsl_ran_gaussian_pdf (x, 1.0);
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#pragma once

#include <stdlib.h>
#include <stdio.h>

#include <MathLib/MathLibDll.h>

#include <MathLib/MathLib.h>
#include <MathLib/ThreeTuple.h>

#include <gsl/matrix/gsl_matrix.h>

#define MATRIX_AT(m, i, j) (*((m->data + ((i) * m->tda + (j)))))

class Vector3d;

/*====================================================================================================================================================================*
 | This class will be used to represent matrices of arbitrary sizes (m rows by n columns) that have elements of type double. The underlying data strucutre used by    |
 | this class is gsl's (Gnu Scientific Library) matrix class. This class also makes use of the ATLAS implementation of BLAS for some operations such as matrix-matrix |
 | multiplication. This class is meant to improve performance, not necessarily ease of use.                                                                           |
 *====================================================================================================================================================================*/

class MATHLIB_DECLSPEC Matrix {
friend class Vector;
friend class ParticleEngine;
protected:

	//this data structure holds all the matrix data
	gsl_matrix *matrix;

public:
	/**
		constructor - creates an m rows by n columns matrix that is not initialized to any particular values
	*/
	Matrix(int m, int n);

	/**
		default constructor
	*/
	Matrix();

	/**
		copy constructor - performs a deep copy of the matrix passed in as a parameter.
	*/
	Matrix(const Matrix& other);

	/**
		destructor.
	*/
	~Matrix();

	/**
		loads the matrix with all zero values.
	*/
	void loadZero();

	/**
		loads the matrix with 1's on the diagonal, 0's everywhere else - note: the matrix doesn't have to be square.
	*/
	void loadIdentity();

	/**
		this method sets the current matrix to a 3x3 matrix that is equal to the outer product of the vectors a and b
	*/
	void setToOuterproduct(const Vector3d& a, const Vector3d& b);

	/**
		this method resizes the current matrix to have m rows and n cols. If not enough space is allocated for it, then a new matrix of
		correct dimensions is allocated. There is no guarantee with regards to the data that is contained in the matrix after a resize operation.
	*/
	void resizeTo(int m, int n);

	/**
		copy operator - performs a deep copy of the matrix passed in as a parameter.
	*/
	Matrix& operator=(const Matrix &other);

	/**
		this method performs a shallow copy of the matrix that is passed in as a parameter.
	*/
	void shallowCopy(const Matrix& other, int startRow = 0, int startCol = 0, int endRow = -1, int endCol = -1);

	/**
		this method performs a deep copy of the matrix that is passed in as a paramerer.
	*/
	void deepCopy(const Matrix& other);

	/**
		Returns the number of columns
	*/
	int getColumnCount() const;

	/**
		Returns the number of rows
	*/
	int getRowCount()const;

	/**
		Multiplies each element in the current matrix by a constant
	*/
	void multiplyBy(const double val);

	/**
		This method sets the current matrix to be equal to one of the products: a * b, a'*b, a*b' or a'*b'.
		The values of transA and transB indicate which of the matrices are tranposed and which ones are not.
	*/
	void setToProductOf(const Matrix& a, const Matrix& b, bool transA = false, bool transB = false);

	/**
		This method computes the inverse of the matrix a and writes it over the current matrix. The implementation
		for the inverse of a matrix was obtained from Graphite. The parameter t is used as a threshold value for
		determinants, etc so that we still get a result for the inverse of our matrix even if it is very poorly conditioned.
	*/
	void setToInverseOf(const Matrix &a, double t = 0);

	/**
		This method prints the contents of the matrix - testing purpose only.
	*/
	void printMatrix() const;

	/**
		This method sets the current matrix to be a sub-matrix (starting at (i,j) and ending at (i+rows, j+cols)
		of the one that is passed in as a parameter - shallow copy only.
	*/
	void setToSubmatrix(const Matrix &a, int i, int j, int rows, int cols);

	/**
		This method returns a copy of the value of the matrix at (i,j)
	*/
	double get(int i, int j) const;

	/**
		This method sets the value of the matrix at (i,j) to newVal.
	*/
	void set(int i, int j, double newVal);

	/**
		This method is used to set the values in the matrix to the ones that are passed in the array of doubles.
		It is assumed that the array contains the right number of elements and that there is no space between consecutive rows (tda == nrCols).
	*/
	void setValues(double* vals);

	/**
		this method returns a pointer to its internal matrix data structure
	*/
	gsl_matrix* getMatrixPointer() const;

	/**
		Implement this operator to have a quick way of multiplying 3x3 matrices by vectors - used for dynamics for instance
	*/
	Vector3d operator * (const Vector3d &other);

	/**
		Implement a potentially faster function for 3x3 matrix - vector3d multiplication. Result and v need to be different!!
	*/
	void postMultiplyVector(const Vector3d& v, Vector3d& result);

	/**
		Still need to do (if need be): get row/col vector. Initialize a matrix from a vector so that we can create a 1xn matrix easily.
		make sure dgemv works properly, implement vector class, etc, swap, add, scale matrices, etc
	*/

	/**
		this method adds the matrix that is passed in as a parameter to the current matrix. In addition, it scales, both matrices by the two
		numbers that are passed in as parameters:
		*this = a * *this + b * other.
	*/
	void add(const Matrix& other, double scaleA = 1.0, double scaleB = 1.0);

	/**
		this method subtracts the matrix that is passed in as a parameter from the current matrix. In addition, it scales, both matrices by the two
		numbers that are passed in as parameters:
		*this = a * *this - b * other.
	*/
	void sub(const Matrix& other, double scaleA = 1.0, double scaleB = 1.0);

};

MATHLIB_DECLSPEC void testMatrixClass();

ext: h

max_stars_repo_name: MontyThibault/centre-of-mass-awareness

max_issues_repo_name: MontyThibault/centre-of-mass-awareness

max_forks_repo_name: MontyThibault/centre-of-mass-awareness


text:
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <stdlib.h>

/* Representation */
gsl_complex *alloc_gsl_complex(void)
{
  gsl_complex *c = malloc(sizeof(gsl_complex));
  return c;
}

void free_gsl_complex(gsl_complex *c)
{
  free(c);
}

gsl_complex_float *alloc_gsl_complex_float(void)
{
  gsl_complex_float *c = malloc(sizeof(gsl_complex_float));
  return c;
}

void free_gsl_complex_float(gsl_complex_float *c)
{
  free(c);
}

void mgsl_complex_rect(double x, double y, gsl_complex *res)
{
  gsl_complex ret = gsl_complex_rect(x, y);
  *res = ret;
}

void mgsl_complex_float_rect(float x, float y, gsl_complex_float *res)
{
  GSL_SET_COMPLEX(res, x, y);
}

void mgsl_complex_polar(double r, double theta, gsl_complex *res)
{
  gsl_complex ret = gsl_complex_polar(r, theta);
  *res = ret;
}

void mgsl_complex_float_polar(float r, float theta, gsl_complex_float *res)
{
  gsl_complex ret = gsl_complex_polar((double)r, (double)theta);
  GSL_SET_COMPLEX(res, ret.dat[0], ret.dat[1]);
}

/* Properties */
double mgsl_complex_arg(gsl_complex *c)
{
  return gsl_complex_arg(*c);
}

double mgsl_complex_abs(gsl_complex *c)
{
  return gsl_complex_abs(*c);
}

double mgsl_complex_abs2(gsl_complex *c)
{
  return gsl_complex_abs2(*c);
}

double mgsl_complex_logabs(gsl_complex *c)
{
  return gsl_complex_logabs(*c);
}

/* Arithmetic operators */

void mgsl_complex_add(gsl_complex *a, gsl_complex *b, gsl_complex *res)
{
  *res = gsl_complex_add(*a, *b);
}

void mgsl_complex_sub(gsl_complex *a, gsl_complex *b, gsl_complex *res)
{
  *res = gsl_complex_sub(*a, *b);
}

void mgsl_complex_mul(gsl_complex *a, gsl_complex *b, gsl_complex *res)
{
  *res = gsl_complex_mul(*a, *b);
}

void mgsl_complex_div(gsl_complex *a, gsl_complex *b, gsl_complex *res)
{
  *res = gsl_complex_div(*a, *b);
}

void mgsl_complex_add_real(gsl_complex *a, double x, gsl_complex *res)
{
  *res = gsl_complex_add_real(*a, x);
}

void mgsl_complex_sub_real(gsl_complex *a, double x, gsl_complex *res)
{
  *res = gsl_complex_sub_real(*a, x);
}

void mgsl_complex_mul_real(gsl_complex *a, double x, gsl_complex *res)
{
  *res = gsl_complex_mul_real(*a, x);
}

void mgsl_complex_div_real(gsl_complex *a, double x, gsl_complex *res)
{
  *res = gsl_complex_div_real(*a, x);
}

void mgsl_complex_add_imag(gsl_complex *a, double y, gsl_complex *res)
{
  *res = gsl_complex_add_imag(*a, y);
}

void mgsl_complex_sub_imag(gsl_complex *a, double y, gsl_complex *res)
{
  *res = gsl_complex_sub_imag(*a, y);
}

void mgsl_complex_mul_imag(gsl_complex *a, double y, gsl_complex *res)
{
  *res = gsl_complex_mul_imag(*a, y);
}

void mgsl_complex_div_imag(gsl_complex *a, double y, gsl_complex *res)
{
  *res = gsl_complex_div_imag(*a, y);
}

void mgsl_complex_conjugate(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_conjugate(*z);
}

void mgsl_complex_inverse(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_inverse(*z);
}

void mgsl_complex_negative(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_negative(*z);
}

/* Elementary Complex Functions */
void mgsl_complex_sqrt(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_sqrt(*z);
}

void mgsl_complex_sqrt_real(double x, gsl_complex *res)
{
  *res = gsl_complex_sqrt_real(x);
}

void mgsl_complex_pow(gsl_complex *z, gsl_complex *a, gsl_complex *res)
{
  *res = gsl_complex_pow(*z, *a);
}

void mgsl_complex_pow_real(gsl_complex *z, double x, gsl_complex *res)
{
  *res = gsl_complex_pow_real(*z, x);
}

void mgsl_complex_exp(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_exp(*z);
}

void mgsl_complex_log(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_log(*z);
}

void mgsl_complex_log10(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_log10(*z);
}

void mgsl_complex_log_b(gsl_complex *z, gsl_complex *b, gsl_complex *res)
{
  *res = gsl_complex_log_b(*z, *b);
}

/* Complex Trigonometric Functions */
void mgsl_complex_sin(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_sin(*z);
}

void mgsl_complex_cos(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_cos(*z);
}

void mgsl_complex_tan(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_tan(*z);
}

void mgsl_complex_sec(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_sec(*z);
}

void mgsl_complex_csc(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_csc(*z);
}

void mgsl_complex_cot(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_cot(*z);
}

/* Inverse Complex Trigonometric Functions */
void mgsl_complex_arcsin(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arcsin(*z);
}

void mgsl_complex_arcsin_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arcsin_real(z);
}

void mgsl_complex_arccos(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccos(*z);
}

void mgsl_complex_arccos_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arccos_real(z);
}

void mgsl_complex_arctan(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arctan(*z);
}

void mgsl_complex_arcsec(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arcsec(*z);
}

void mgsl_complex_arcsec_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arcsec_real(z);
}

void mgsl_complex_arccsc(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccsc(*z);
}

void mgsl_complex_arccsc_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arccsc_real(z);
}

void mgsl_complex_arccot(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccot(*z);
}

/* Complex Hyperbolic Functions */
void mgsl_complex_sinh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_sinh(*z);
}

void mgsl_complex_cosh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_cosh(*z);
}

void mgsl_complex_tanh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_tanh(*z);
}

void mgsl_complex_sech(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_sech(*z);
}

void mgsl_complex_csch(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_csch(*z);
}

void mgsl_complex_coth(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_coth(*z);
}

/* Inverse Complex Hyperbolic Functions */
void mgsl_complex_arcsinh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arcsinh(*z);
}

void mgsl_complex_arccosh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccosh(*z);
}

void mgsl_complex_arccosh_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arccosh_real(z);
}

void mgsl_complex_arctanh(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arctanh(*z);
}

void mgsl_complex_arctanh_real(double z, gsl_complex *res)
{
  *res = gsl_complex_arctanh_real(z);
}

void mgsl_complex_arcsech(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arcsech(*z);
}

void mgsl_complex_arccsch(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccsch(*z);
}

void mgsl_complex_arccoth(gsl_complex *z, gsl_complex *res)
{
  *res = gsl_complex_arccoth(*z);
}

max_stars_repo_name: frithnanth/raku-Math-Libgsl-Complex

max_issues_repo_name: frithnanth/raku-Math-Libgsl-Complex

max_forks_repo_name: frithnanth/raku-Math-Libgsl-Complex


text:
#include "linear.h"
#include <assert.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <lauxlib.h>
#include <cblas.h>
#include <lapacke.h>

/* matrix orders */
static const char * const ORDERS[] = { "row", "col", NULL };

/* checks an order */
static CBLAS_ORDER checkorder (lua_State *L, int index) {

	switch (luaL_checkoption(L, index, "row", ORDERS)) {
	case 0:
		return CblasRowMajor;

	case 1:
		return CblasColMajor;
	}

	/* not reached */
	assert(0);
	return (CBLAS_ORDER)0;
}

/* checks a transpose */
static CBLAS_TRANSPOSE checktranspose (lua_State *L, int index) {
	static const char * const TRANSPOSES[] = { "notrans", "trans", NULL };

	switch (luaL_checkoption(L, index, "notrans", TRANSPOSES)) {
	case 0:
		return CblasNoTrans;

	case 1:
		return CblasTrans;
	}

	/* not reached */
	assert(0);
	return (CBLAS_TRANSPOSE)0;
}

/* translates a transpose for LAPACK */
static char lapacktranspose (CBLAS_TRANSPOSE transpose) {
	switch (transpose) {
	case CblasNoTrans:
		return 'N';

	case CblasTrans:
		return 'T';

	default:
		/* not reached */
		assert(0);
		return '\0';
	}
}

/* returns an int value from a table */
static int intvalue (lua_State *L, const char *key, int dfl) {
	int result, isinteger;

	lua_getfield(L, -1, key);
	if (!lua_isnil(L, -1)) {
		result = lua_tointegerx(L, -1, &isinteger);
		if (!isinteger) {
			luaL_error(L, "bad field " LUA_QS, key);
		}
	} else {
		if (dfl < 0) {
			luaL_error(L, "missing field " LUA_QS, key);
		}
		result = dfl;
	}
	lua_pop(L, 1);
	return result;
}

/* returns an option value from a table */
static int optionvalue (lua_State *L, const char *key, const char *dfl,
		const char *options[]) {
	const char *str;
	int i;

	lua_getfield(L, -1, key);
	if (!lua_isnil(L, -1)) {
		str = lua_tostring(L, -1);
		if (str == NULL) {
			luaL_error(L, "bad field " LUA_QS, key);
		}
	} else {
		if (dfl == NULL) {
			luaL_error(L, "missing field " LUA_QS, key);
		}
		str = dfl;
	}
	lua_pop(L, 1);
	for (i = 0; options[i] != NULL; i++) {
		if (strcmp(options[i], str) == 0) {
			return i;
		}
	}
	luaL_error(L, "bad option " LUA_QS " in field " LUA_QS, str, key);
	return 0; /* not reached */
}

/* raises a linear argument error */
static int argerror (lua_State *L, int index) {
	return luaL_argerror(L, index, lua_pushfstring(L, "vector, or matrix "
			"expected, got %s", luaL_typename(L, index)));
}

/* pushes a new vector onto the stack */
static struct vector *newvector (lua_State *L, int size) {
	return lualinear_newvector(L, size);
}

/* pushes an existing vector onto the stack */
static struct vector *wrapvector (lua_State *L, int size, float *values) {
	return lualinear_wrapvector(L, size, values);
}

/* creates a new vector */
static int vector (lua_State *L) {
	int size;

	/* process arguments */
	size = luaL_checkinteger(L, 1);
	luaL_argcheck(L, size >= 1, 1, "bad dimension");

	/* create */
	newvector(L, size);
	return 1;
}

/* vector length implementation */
static int vector_len (lua_State *L) {
	struct vector *x;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	lua_pushinteger(L, x->size);
	return 1;
}

/* vector index implementation */
static int vector_index (lua_State *L) {
	struct vector *x;
	int index;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	index = luaL_checkinteger(L, 2);
	luaL_argcheck(L, index >= 1 && index <= x->size, 2, "bad index");
	lua_pushnumber(L, x->values[(size_t)(index - 1) * x->inc]);
	return 1;
}

/* matrix vector newindex implementation */
static int vector_newindex (lua_State *L) {
	struct vector *x;
	int index;
	float value;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	index = luaL_checkinteger(L, 2);
	luaL_argcheck(L, index >= 1 && index <= x->size, 2, "bad index");
	value = luaL_checknumber(L, 3);
	x->values[(size_t)(index - 1) * x->inc] = value;
	return 0;
}

/* vector next function */
static int vector_next (lua_State *L) {
	struct vector *x;
	int index;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	index = luaL_checkinteger(L, 2);
	if (index >= 0 && index < x->size) {
		lua_pushinteger(L, index + 1);
		lua_pushnumber(L, x->values[(size_t)index]);
		return 2;
	}
	lua_pushnil(L);
	return 1;
}

/* vector ipairs function */
static int vector_ipairs (lua_State *L) {
	luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	lua_pushcfunction(L, vector_next);
	lua_pushvalue(L, 1);
	lua_pushinteger(L, 0);
	return 3;
}

/* returns the string representation of a vector */
static int vector_tostring (lua_State *L) {
	struct vector *x;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	lua_pushfstring(L, "vector: %p", x);
	return 1;
}

/* frees a vector */
static int vector_free (lua_State *L) {
	struct vector *x;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x->ref == LUA_NOREF) {
		free(x->values);
	} else {
		luaL_unref(L, LUA_REGISTRYINDEX, x->ref);
	}
	return 0;
}

/* pushes a new matrix onto the stack */
static struct matrix *newmatrix (lua_State *L, int rows, int cols,
		CBLAS_ORDER order) {
	return lualinear_newmatrix(L, rows, cols, order);
}

/* pushes an existing matrix onto the stack */
static struct matrix *wrapmatrix (lua_State *L, int rows, int cols,
		CBLAS_ORDER order, float *values) {
	return lualinear_wrapmatrix(L, rows, cols, order, values);
}

/* creates a new matrix */
static int matrix (lua_State *L) {
	int rows, cols;
	CBLAS_ORDER order;

	/* process arguments */
	rows = luaL_checkinteger(L, 1);
	luaL_argcheck(L, rows >= 1, 1, "bad dimension");
	cols = luaL_checkinteger(L, 2);
	luaL_argcheck(L, cols >= 1, 2, "bad dimension");
	order = checkorder(L, 3);

	/* create */
	newmatrix(L, rows, cols, order);
	return 1;
}

/* returns the length of a matrix */
static int matrix_len (lua_State *L) {
	struct matrix *X;

	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	switch (X->order) {
	case CblasRowMajor:
		lua_pushinteger(L, X->rows);
		break;

	case CblasColMajor:
		lua_pushinteger(L, X->cols);
		break;
	}
	return 1;
}

/* matrix index implementation */
static int matrix_index (lua_State *L) {
	struct matrix *X;
	int index, size;
	struct vector *x;

	/* process arguments */
	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	index = luaL_checkinteger(L, 2);
	luaL_argcheck(L, index >= 1, 2, "bad index");
	switch (X->order) {
	case CblasRowMajor:
		luaL_argcheck(L, index <= X->rows, 2, "bad index");
		size = X->cols;
		break;

	case CblasColMajor:
		luaL_argcheck(L, index <= X->cols, 2, "bad index");
		size = X->rows;
		break;

	default:
		/* not reached */
		size = -1;
		assert(0);
	}

	/* create vector */
	x = wrapvector(L, size, &X->values[(size_t)(index - 1) * X->ld]);
	lua_pushvalue(L, 1);
	x->ref = luaL_ref(L, LUA_REGISTRYINDEX);
	return 1;
}

/* matrix next function */
static int matrix_next (lua_State *L) {
	struct matrix *X;
	int index, majorsize, minorsize;
	struct vector *x;

	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	index = luaL_checkinteger(L, 2);
	switch (X->order) {
	case CblasRowMajor:
		majorsize = X->rows;
		minorsize = X->cols;
		break;

	case CblasColMajor:
		majorsize = X->cols;
		minorsize = X->rows;
		break;

	default:
		/* not reached */
		assert(0);
		return 0;
	}
	if (index >= 0 && index < majorsize) {
		lua_pushinteger(L, index + 1);
		x = wrapvector(L, minorsize, &X->values[(size_t)index * X->ld]);
		lua_pushvalue(L, 1);
		x->ref = luaL_ref(L, LUA_REGISTRYINDEX);
		return 2;
	}
	lua_pushnil(L);
	return 1;
}

/* matrix ipairs function */
static int matrix_ipairs (lua_State *L) {
	luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	lua_pushcfunction(L, matrix_next);
	lua_pushvalue(L, 1);
	lua_pushinteger(L, 0);
	return 3;
}

/* returns the string representation of a matrix */
static int matrix_tostring (lua_State *L) {
	struct matrix *X;

	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	lua_pushfstring(L, "matrix: %p", X);
	return 1;
}

/* frees a matrix */
static int matrix_free (lua_State *L) {
	struct matrix *X;

	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X->ref == LUA_NOREF) {
		free(X->values);
	} else {
		luaL_unref(L, LUA_REGISTRYINDEX, X->ref);
	}
	return 0;
}

/* returns the type of a linear object */
static int type (lua_State *L) {
	if (luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE) != NULL) {
		lua_pushliteral(L, "vector");
		return 1;
	}
	if (luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE) != NULL) {
		lua_pushliteral(L, "matrix");
		return 1;
	}
	lua_pushnil(L);
	return 1;
}

/* returns the size of a linear object */
static int size (lua_State *L) {
	struct vector *x;
	struct matrix *X;

	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		lua_pushinteger(L, x->size);
		return 1;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		lua_pushinteger(L, X->rows);
		lua_pushinteger(L, X->cols);
		lua_pushstring(L, ORDERS[X->order == CblasRowMajor ? 0 : 1]);
		return 3;
	}
	return argerror(L, 1);
}

/* transposed vector */
static int tvector (lua_State *L) {
	struct matrix *X;
	int index, size;
	struct vector *x;

	/* process arguments */
	X = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	index = luaL_checkinteger(L, 2);
	luaL_argcheck(L, index >= 1, 2, "bad index");
	switch (X->order) {
	case CblasRowMajor:
		luaL_argcheck(L, index <= X->cols, 2, "bad index");
		size = X->rows;
		break;

	case CblasColMajor:
		luaL_argcheck(L, index <= X->rows, 2, "bad index");
		size = X->cols;
		break;

	default:
		/* not reached */
		size = -1;
		assert(0);
	}

	/* create vector */
	x = wrapvector(L, size, &X->values[index - 1]);
	x->inc = X->ld;
	lua_pushvalue(L, 1);
	x->ref = luaL_ref(L, LUA_REGISTRYINDEX);
	return 1;
}

/* subvector or submatrix */
static int sub (lua_State *L) {
	struct vector *x, *s;
	struct matrix *X, *S;

	/* process arguments */
	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		int start, end;

		start = luaL_optinteger(L, 2, 1);
		luaL_argcheck(L, start >= 1 && start <= x->size, 2,
				"bad index");
		end = luaL_optinteger(L, 3, x->size);
		luaL_argcheck(L, end >= start && end <= x->size, 3,
				"bad index");
		s = wrapvector(L, end - start + 1, &x->values[
				(size_t)(start - 1) * x->inc]);
		s->inc = x->inc;
		lua_pushvalue(L, 1);
		s->ref = luaL_ref(L, LUA_REGISTRYINDEX);
		return 1;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		int rowstart, rowend, colstart, colend;

		switch (X->order){
		case CblasRowMajor:
			rowstart = luaL_optinteger(L, 2, 1);
			luaL_argcheck(L, rowstart >= 1 && rowstart <= X->rows,
					2, "bad index");
			colstart = luaL_optinteger(L, 3, 1);
			luaL_argcheck(L, colstart >= 1 && colstart <= X->cols,
					3, "bad index");
			rowend = luaL_optinteger(L, 4, X->rows);
			luaL_argcheck(L, rowend >= rowstart && rowend
					<= X->rows, 4, "bad index");
			colend = luaL_optinteger(L, 5, X->cols);
			luaL_argcheck(L, colend >= colstart && colend
					<= X->cols, 5, "bad index");
			S = wrapmatrix(L, rowend - rowstart + 1, colend
					- colstart + 1, X->order, &X->values[
					(size_t)(rowstart - 1) * X->ld
					+ colstart - 1]);
			break;

		case CblasColMajor:
			colstart = luaL_optinteger(L, 2, 1);
			luaL_argcheck(L, colstart >= 1 && colstart <= X->cols,
					2, "bad index");
			rowstart = luaL_optinteger(L, 3, 1);
			luaL_argcheck(L, rowstart >= 1 && rowstart <= X->rows,
					3, "bad index");
			colend = luaL_optinteger(L, 4, X->cols);
			luaL_argcheck(L, colend >= colstart && colend
					<= X->cols, 4, "bad index");
			rowend = luaL_optinteger(L, 5, X->rows);
			luaL_argcheck(L, rowend >= rowstart && rowend
					<= X->rows, 5, "bad index");
			S = wrapmatrix(L, rowend - rowstart + 1, colend
					- colstart + 1, X->order, &X->values[
					(size_t)(colstart - 1) * X->ld
					+ rowstart - 1]);
			break;

		default:
			/* not reached */
			assert(0);
			return 0;
		}
		S->ld = X->ld;
		lua_pushvalue(L, 1);
		S->ref = luaL_ref(L, LUA_REGISTRYINDEX);
		return 1;
	}
	return argerror(L, 1);
}

/* unwinds matrices into a vector */
static int unwind (lua_State *L) {
	struct vector *x;
	int index, i, j, k;
	size_t base;
	struct matrix *X;

	if (lua_gettop(L) == 0) {
		return luaL_error(L, "wrong number of arguments");
	}
	x = luaL_checkudata(L, lua_gettop(L), LUALINEAR_VECTOR_METATABLE);
	index = 1;
	i = 0;
	while (i < x->size) {
		X = luaL_checkudata(L, index, LUALINEAR_MATRIX_METATABLE);
		luaL_argcheck(L, X->rows * X->cols <= x->size - i, index,
				"matrix too large");
		switch (X->order) {
		case CblasRowMajor:
			for (j = 0; j < X->rows; j++) {
				base = (size_t)j * X->ld;
				for (k = 0; k < X->cols; k++) {
					x->values[(size_t)i * x->inc]
							= X->values[base + k];
					i++;
				}
			}
			break;

		case CblasColMajor:
			for (j = 0; j < X->cols; j++) {
				base = (size_t)j * X->ld;
				for (k = 0; k < X->rows; k++) {
					x->values[(size_t)i * x->inc]
							= X->values[base + k];
					i++;
				}
			}
			break;
		}
		index++;
	}
	return 0;
}

/* reshapes a vector into matrices */
static int reshape (lua_State *L) {
	struct vector *x;
	int index, i, j, k;
	size_t base;
	struct matrix *X;

	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	index = 2;
	i = 0;
	while (i < x->size) {
		X = luaL_checkudata(L, index, LUALINEAR_MATRIX_METATABLE);
		luaL_argcheck(L, X->rows * X->cols <= x->size - i, index,
				"matrix too large");
		switch (X->order) {
		case CblasRowMajor:
			for (j = 0; j < X->rows; j++) {
				base = (size_t)j * X->ld;
				for (k = 0; k < X->cols; k++) {
					X->values[base + k] = x->values[
							(size_t)i * x->inc];
					i++;
				}
			}
			break;

		case CblasColMajor:
			for (j = 0; j < X->cols; j++) {
				base = (size_t)j * X->ld;
				for (k = 0; k < X->rows; k++) {
					X->values[base + k] = x->values[
							(size_t)i * x->inc];
					i++;
				}
			}
			break;
		}
		index++;
	}
	return 0;
}

/* converts a vector or matrix to a table */
static int totable (lua_State *L) {
	struct vector *x;
	struct matrix *X;
	int i, j;
	const float *value;

	/* check and process arguments */
	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		lua_createtable(L, 0, 3);
		lua_pushliteral(L, "vector");
		lua_setfield(L, -2, "type");
		lua_pushinteger(L, x->size);
		lua_setfield(L, -2, "length");
		lua_createtable(L, x->size, 0);
		value = x->values;
		for (i = 0; i < x->size; i++) {
			lua_pushnumber(L, *value);
			lua_rawseti(L, -2, i + 1);
			value += x->inc;
		}
		lua_setfield(L, -2, "values");
		return 1;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		lua_createtable(L, 0, 5);
		lua_pushliteral(L, "matrix");
		lua_setfield(L, -2, "type");
		lua_pushinteger(L, X->rows);
		lua_setfield(L, -2, "rows");
		lua_pushinteger(L, X->cols);
		lua_setfield(L, -2, "cols");
		switch (X->order) {
		case CblasRowMajor:
			lua_pushliteral(L, "rowmajor");
			lua_setfield(L, -2, "order");
			lua_createtable(L, X->rows, 0);
			for (i = 0; i < X->rows; i++) {
				lua_createtable(L, X->cols, 0);
				value = &X->values[(size_t)i * X->ld];
				for (j = 0; j < X->cols; j++) {
					lua_pushnumber(L, *value++);
					lua_rawseti(L, -2, j + 1);
				}
				lua_rawseti(L, -2, i + 1);
			}
			lua_setfield(L, -2, "values");
			break;

		case CblasColMajor:
			lua_pushliteral(L, "colmajor");
			lua_setfield(L, -2, "order");
			lua_createtable(L, X->cols, 0);
			for (i = 0; i < X->cols; i++) {
				lua_createtable(L, X->rows, 0);
				value = &X->values[(size_t)i * X->ld];
				for (j = 0; j < X->rows; j++) {
					lua_pushnumber(L, *value++);
					lua_rawseti(L, -2, j + 1);
				}
				lua_rawseti(L, -2, i + 1);
			}
			lua_setfield(L, -2, "values");
			break;
		}
		return 1;
	}
	return argerror(L, 1);
}

/* converts a table to a vector or matrix */
static int tolinear (lua_State *L) {
	static const char *types[] = { "vector", "matrix", NULL };
	static const char *orders[] = { "rowmajor", "colmajor", NULL };
	struct vector *x;
	struct matrix *X;
	int size, rows, cols, major, minor;
	CBLAS_ORDER order;
	int i, j;
	int isnum;
	float *value;

	/* check arguments */
	luaL_checktype(L, 1, LUA_TTABLE);
	lua_settop(L, 1);

	/* handle types */
	switch (optionvalue(L, "type", NULL, types)) {
	case 0: /* vector */
		size = intvalue(L, "length", -1);
		if (size < 1) {
			return luaL_error(L, "bad field " LUA_QS, "length");
		}
		x = newvector(L, size);
		lua_getfield(L, 1, "values");
		if (lua_type(L, -1) != LUA_TTABLE) {
			return luaL_error(L, "bad field " LUA_QS, "values");
		}
		value = x->values;
		for (i = 0; i < size; i++) {
			lua_rawgeti(L, -1, i + 1);
			*value++ = lua_tonumberx(L, -1, &isnum);
			if (!isnum) {
				return luaL_error(L, "bad value at index %d",
						i + 1);
			}
			lua_pop(L, 1);
		}
		lua_pop(L, 1);
		return 1;

	case 1: /* matrix */
		rows = intvalue(L, "rows", -1);
		if (rows < 1) {
			return luaL_error(L, "bad field " LUA_QS, "rows");
		}
		cols = intvalue(L, "cols", -1);
		if (cols < 1) {
			return luaL_error(L, "bad field " LUA_QS, "cols");
		}
		switch (optionvalue(L, "order", NULL, orders)) {
		case 0:
			order = CblasRowMajor;
			major = rows;
			minor = cols;
			break;

		case 1:
			order = CblasColMajor;
			major = cols;
			minor = rows;
			break;

		default:
			/* not reched */
			assert(0);
			return 0;
		}
		X = newmatrix(L, rows, cols, order);
		lua_getfield(L, 1, "values");
		if (lua_type(L, -1) != LUA_TTABLE) {
			return luaL_error(L, "bad field " LUA_QS, "values");
		}
		for (i = 0; i < major; i++) {
			value = &X->values[i * X->ld];
			lua_rawgeti(L, -1, i + 1);
			if (lua_type(L, -1) != LUA_TTABLE) {
				return luaL_error(L, "bad value at index %d",
						i + 1);
			}
			for (j = 0; j < minor; j++) {
				lua_rawgeti(L, -1, j + 1);
				*value++ = lua_tonumberx(L, -1, &isnum);
				if (!isnum) {
					return luaL_error(L, "bad value at "
							"index (%d,%d)", i + 1,
							j + 1);
				}
				lua_pop(L, 1);
			}
			lua_pop(L, 1);
		}
		lua_pop(L, 1);
		return 1;
	}

	/* not reached */
	assert(0);
	return 0;
}

/* invokes the DOT subprogram (x' y) */
static int dot (lua_State *L) {
	struct vector *x, *y;
	float dot;

	/* check and process arguments */
	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	y = luaL_checkudata(L, 2, LUALINEAR_VECTOR_METATABLE);
	luaL_argcheck(L, y->size == x->size, 2, "dimension mismatch");

	/* invoke subprogram */
	dot = cblas_sdot(x->size, x->values, x->inc, y->values, y->inc);
	lua_pushnumber(L, dot);
	return 1;
}

/* invokes the NRM2 subprogram (||x||_2) */
static int nrm2 (lua_State *L) {
	struct vector *x;
	float nrm2;

	/* check and process arguments */
	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);

	/* invoke subprogram */
	nrm2 = cblas_snrm2(x->size, x->values, x->inc);
	lua_pushnumber(L, nrm2);
	return 1;
}

/* invokes the ASUM subprogram (sigma |x|) */
static int asum (lua_State *L) {
	struct vector *x;
	float asum;

	/* check and process arguments */
	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);

	/* invoke subprogram */
	asum = cblas_sasum(x->size, x->values, x->inc);
	lua_pushnumber(L, asum);
	return 1;
}

/* invokes the IAMAX subprogram (argmax |x|) */
static int iamax (lua_State *L) {
	struct vector *x;
	int iamax;

	/* check and process arguments */
	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);

	/* invoke subprogram */
	iamax = cblas_isamax(x->size, x->values, x->inc);
	lua_pushinteger(L, iamax + 1);
	return 1;
}

/* sum implementation */
static float _sum (const float *values, int size, int inc) {
	float sum;
	int i;

	sum = 0.0;
	#pragma omp parallel for private(i) schedule(auto) \
			 if(size >= LUALINEAR_OMP_MINSIZE) reduction(+:sum)
	for (i = 0; i < size; i++) {
		sum += values[(size_t)i * inc];
	}
	return sum;
}

/* sum implementation (sigma x_i) */
static int sum (lua_State *L) {
	struct vector *x, *y;
	struct matrix *X;
	int i;

	/* check and process arguments */
	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		lua_pushnumber(L, _sum(x->values, x->size, x->inc));
		return 1;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		y = luaL_checkudata(L, 2, LUALINEAR_VECTOR_METATABLE);
		switch (checktranspose(L, 3)) {
		case CblasNoTrans:
			switch (X->order) {
			case CblasRowMajor:
				luaL_argcheck(L, y->size == X->rows, 2,
						"dimension mismatch");
				for (i = 0; i < X->rows; i++) {
					y->values[(size_t)i * y->inc] = _sum(
							&X->values[(size_t)i
							* X->ld], X->cols, 1);
				}
				break;

			case CblasColMajor:
				luaL_argcheck(L, y->size == X->cols, 2,
						"dimension mismatch");
				for (i = 0; i < X->cols; i++) {
					y->values[(size_t)i * y->inc] = _sum(
							&X->values[(size_t)i
							* X->ld], X->rows, 1);
				}
				break;
			}
			break;

		case CblasTrans:
			switch (X->order) {
			case CblasRowMajor:
				luaL_argcheck(L, y->size == X->cols, 2,
						"dimension mismatch");
				for (i = 0; i < X->cols; i++) {
					y->values[(size_t)i * y->inc] = _sum(
							&X->values[(size_t)i],
							X->rows, X->ld);
				}
				break;

			case CblasColMajor:
				luaL_argcheck(L, y->size == X->rows, 2,
						"dimension mismatch");
				for (i = 0; i < X->rows; i++) {
					y->values[(size_t)i * y->inc] = _sum(
							&X->values[(size_t)i],
							X->cols, X->ld);
				}
				break;
			}
			break;

		default:
			/* not reached */
			assert(0);
			break;
		}
		return 0;
	}
	return argerror(L, 1);
}

/* xy function */
typedef void(*xyfunction)(int, float *, int, float *, int, float);

/* invokes an (x,y) subproram */
static int xy (lua_State *L, xyfunction s, int hasy, int hasalpha) {
	int index, i;
	float alpha;
	struct vector *x, *y;
	struct matrix *X, *Y;

	/* check and process arguments */
	index = 2;
	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		if (hasy) {
			y = luaL_testudata(L, 2, LUALINEAR_VECTOR_METATABLE);
			Y = luaL_testudata(L, 2, LUALINEAR_MATRIX_METATABLE);
			if (y == NULL && Y == NULL) {
				return argerror(L, 2);
			}
			index++;
		} else {
			y = x;
			Y = NULL;
		}
		if (hasalpha) {
			alpha = luaL_optnumber(L, index, 1.0);
			index++;
		} else {
			alpha = 0.0;
		}
		if (y != NULL) {
			/* invoke subprogram on vector-vector */
			luaL_argcheck(L, y->size == x->size, 2,
					"dimension mismatch");
			s(x->size, x->values, x->inc, y->values, y->inc, alpha);
			return 0;
		}

		/* invoke subprogram on vector-matrix */
		switch (checktranspose(L, index)) {
		case CblasNoTrans:
			switch (Y->order) {
			case CblasRowMajor:
				luaL_argcheck(L, 1, x->size == Y->cols,
						"dimension mismatch");
				for (i = 0; i < Y->rows; i++) {
					s(x->size, x->values, x->inc,
							&Y->values[(size_t)i
							* Y->ld], 1, alpha);
				}
				break;

			case CblasColMajor:
				luaL_argcheck(L, 1, x->size == Y->rows,
						"dimension mismatch");
				for (i = 0; i < Y->cols; i++) {
					s(x->size, x->values, x->inc,
							&Y->values[(size_t)i
							* Y->ld], 1, alpha);
				}
				break;
			}
			break;

		case CblasTrans:
			switch (Y->order) {
			case CblasRowMajor:
				luaL_argcheck(L, 1, x->size == Y->rows,
						"dimension mismatch");
				for (i = 0; i < Y->rows; i++) {
					s(x->size, x->values, x->inc,
							&Y->values[(size_t)i],
							Y->ld, alpha);
				}
				break;

			case CblasColMajor:
				luaL_argcheck(L, 1, x->size == Y->cols,
						"dimension mismatch");
				for (i = 0; i < Y->cols; i++) {
					s(x->size, x->values, x->inc,
							&Y->values[(size_t)i],
							Y->ld, alpha);
				}
				break;
			}
			break;

		default:
			/* not reached */
			assert(0);
		}
		return 0;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		if (hasy) {
			Y = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
			luaL_argcheck(L, X->order == Y->order, 2,
					"order mismatch");
			luaL_argcheck(L, X->rows == Y->rows && X->cols
					== Y->cols, 2, "dimension mismatch");
			index++;
		} else {
			Y = X;
		}
		if (hasalpha) {
			alpha = luaL_optnumber(L, index, 1.0);
			index++;
		} else {
			alpha = 0.0;
		}

		/* invoke subprogram on matrix-matrix */
		switch (X->order) {
		case CblasRowMajor:
			for (i = 0; i < X->rows; i++) {
				s(X->cols, &X->values[(size_t)i * X->ld], 1,
						&Y->values[(size_t)i * Y->ld],
						1, alpha);
			}
			break;

		case CblasColMajor:
			for (i = 0; i < X->cols; i++) {
				s(X->rows, &X->values[(size_t)i * X->ld], 1,
						&Y->values[(size_t)i * Y->ld],
						1, alpha);
			}
			break;
		}
		return 0;
	}
	return argerror(L, 1);
}

/* wraps the SWAP subprogram */
static void _swap (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	(void)alpha;
	cblas_sswap(size, x, incx, y, incy);
}

/* invokes the SWAP subprogram (y <-> x) */
static int swap (lua_State *L) {
	return xy(L, _swap, 1, 0);
}

/* wraps the COPY subprogram */
static void _copy (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	(void)alpha;
	cblas_scopy(size, x, incx, y, incy);
}

/* invokes the COPY subprogram (y <- x) */
static int copy (lua_State *L) {
	return xy(L, _copy, 1, 0);
}

/* wraps the AXPY subprogram */
static void _axpy (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	cblas_saxpy(size, alpha, x, incx, y, incy);
}

/* invokes the AXPY subprogram (y <- alpha x + y) */
static int axpy (lua_State *L) {
	return xy(L, _axpy, 1, 1);
}

/* wraps the SCAL subprogram */
static void _scal (int size, float *x, int incx, float *y, int incy,
			float alpha) {
	(void)y;
	(void)incy;
	cblas_sscal(size, alpha, x, incx);
}

/* invokes the SCAL subprogram (x <- alpha x) */
static int scal (lua_State *L) {
	return xy(L, _scal, 0, 1);
}

/* set operation implementation */
static void _set (int size, float *x, int incx, float *y, int incy,
			float alpha) {
	int i;

	(void)y;
	(void)incy;
	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		x[(size_t)i * incx] = alpha;
	}
}

/* performs a set operation (x <- alpha) */
static int set (lua_State *L) {
	return xy(L, _set, 0, 1);
}

/* uniform RNG implementation */
static void _uniform (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	int i;

	(void)y;
	(void)incy;
	(void)alpha;
	for (i = 0; i < size; i++) {
		*x = (float)random() * (1.0 / ((float)RAND_MAX + 1.0));
		x += incx;
	}
}

/* performs a uniform operation (x <- uniform) */
static int uniform (lua_State *L) {
	return xy(L, _uniform, 0, 0);
}

/* normal RNG implementation */
static void _normal (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	int i;
	float u1, u2, r, s, c;

	(void)y;
	(void)incy;
	(void)alpha;
	for (i = 0; i < size - 1; i += 2) {
		do {
			u1 = (float)random() * (1.0 / (float)RAND_MAX);
			u2 = (float)random() * (1.0 / (float)RAND_MAX);
		} while (u1 <= -DBL_MAX);
		r = sqrt(-2.0 * logf(u1));
		sincosf(2 * M_PI * u2, &s, &c);
		*x = r * c;
		x += incx;
		*x = r * s;
		x += incx;
	}
	if (i < size) {
		do {
			u1 = (float)random() * (1.0 / (float)RAND_MAX);
			u2 = (float)random() * (1.0 / (float)RAND_MAX);
		} while (u1 <= -DBL_MAX);
		*x = sqrtf(-2.0 * logf(u1)) * cosf(2 * M_PI * u2);
		x += incx;
	}
}

/* performs a normal operation (x <- normal) */
static int normal (lua_State *L) {
	return xy(L, _normal, 0, 0);
}

/* inc operation implementation */
static void _inc (int size, float *x, int incx, float *y, int incy,
			float alpha) {
	int i;

	(void)y;
	(void)incy;
	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		x[(size_t)i * incx] += alpha;
	}
}

/* performs a inc operation (x <- x + alpha) */
static int inc (lua_State *L) {
	return xy(L, _inc, 0, 1);
}

/* element-wise multiplication implementation, alpha = 1 */
static void _mul1 (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	int i;

	(void)alpha;
	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		y[(size_t)i * incy] *= x[(size_t)i * incx];
	}
}

/* element-wise multiplication implementation, alpha = -1 */
static void _mulm1 (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	int i;

	(void)alpha;
	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		y[(size_t)i * incy] /= x[(size_t)i * incx];
	}
}

/* element-wise multiplication implementation, alpha = any */
static void _mul (int size, float *x, int incx, float *y, int incy,
		float alpha) {
	int i;

	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		y[(size_t)i * incy] *= pow(x[(size_t)i * incx], alpha);
	}
}

/* performs element-wise multiplication (y <- x^alpha .* y) */
static int mul (lua_State *L) {
	float alpha;

	alpha = luaL_optnumber(L, 3, 1.0);
	if (alpha == 1.0) {
		return xy(L, _mul1, 1, 1);
	}
	if (alpha == -1.0) {
		return xy(L, _mulm1, 1, 1);
	}
	return xy(L, _mul, 1, 1);
}

/* power raising operation implementation */
static void _pow (int size, float *x, int incx, float *y, int incy,
			float alpha) {
	int i;

	(void)y;
	(void)incy;
	#pragma omp parallel for private(i) schedule(auto) \
			if(size >= LUALINEAR_OMP_MINSIZE)
	for (i = 0; i < size; i++) {
		x[(size_t)i * incx] = pow(x[(size_t)i * incx], alpha);
	}
}

/* performs element-wise power raising (x <- x^alpha) */
static int powx (lua_State *L) {
	return xy(L, _pow, 0, 1);
}

/* apply function */
typedef float(*applyfunction)(float);

/* applies a function to a value */
static int apply (lua_State *L, applyfunction apply, int parallel) {
	struct vector *x;
	struct matrix *X;
	int i, j;
	size_t base;

	if (lua_type(L, 1) == LUA_TNUMBER) {
		lua_pushnumber(L, apply(lua_tonumber(L, 1)));
		return 1;
	}
	x = luaL_testudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	if (x != NULL) {
		#pragma omp parallel for private(i) schedule(auto) \
				if(parallel && x->size >= LUALINEAR_OMP_MINSIZE)
		for (i = 0; i < x->size; i++) {
			x->values[(size_t)i * x->inc] =
					apply(x->values[(size_t)i
					* x->inc]);
		}
		return 0;
	}
	X = luaL_testudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	if (X != NULL) {
		switch (X->order) {
		case CblasRowMajor:
			for (i = 0; i < X->rows; i++) {
				base = (size_t)i * X->ld;
				#pragma omp parallel for private(j) \
						schedule(auto) \
						if(parallel && X->cols \
						>= LUALINEAR_OMP_MINSIZE)
				for (j = 0; j < X->cols; j++) {
					X->values[base + j] = apply(
							X->values[base
							+ j]);
				}
			}
			break;

		case CblasColMajor:
			for (i = 0; i < X->cols; i++) {
				base = (size_t)i * X->ld;
				#pragma omp parallel for private(j) \
						schedule(auto) \
						if(parallel && X->rows \
						>= LUALINEAR_OMP_MINSIZE)
				for (j = 0; j < X->rows; j++) {
					X->values[base + j] = apply(
							X->values[base
							+ j]);
				}
			}
			break;
		}
		return 0;
	}
	return luaL_argerror(L, 1, lua_pushfstring(L, "number, vector, or "
			"matrix expected, got %s", luaL_typename(L, 1)));

}

/* sign function implementation */
static float _sign (float x) {
	if (x > 0) {
		return 1;
	}
	if (x < 0) {
		return -1;
	}
	return x;
}

/* sign function */
static int sign (lua_State *L) {
	return apply(L, _sign, 1);
}

/* abs function implementation */
static float _abs (float x) {
	return fabs(x);
}

/* abs function */
static int absx (lua_State *L) {
	return apply(L, _abs, 1);
}

/* exp function */
static int expx (lua_State *L) {
	return apply(L, expf, 1);
}

/* log function */
static int logx (lua_State *L) {
	return apply(L, logf, 1);
}

/* logistic function implementation */
static float _logistic (float z) {
	return 1.0 / (1.0 + expf(-z));
}

/* logistic function */
static int logistic (lua_State *L) {
	return apply(L, _logistic, 1);
}

/* tanh function */
static int tanhx (lua_State *L) {
	return apply(L, tanhf, 1);
}

/* softplus function implementation */
static float _softplus (float x) {
	return logf(1 + expf(x));
}

/* softplus function */
static int softplus (lua_State *L) {
	return apply(L, _softplus, 1);
}

/* rectifier function implementation */
static float _rectifier (float x) {
	return x > 0.0 ? x : 0.0;
}

/* rectifier function */
static int rectifier (lua_State *L) {
	return apply(L, _rectifier, 1);
}

/* current Lua state */
static __thread lua_State *TL;

/* apply function implementation */
static float _apply (float x) {
	float result;

	lua_pushvalue(TL, -1);
	lua_pushnumber(TL, x);
	lua_call(TL, 1, 1);
	result = lua_tonumber(TL, -1);
	lua_pop(TL, 1);
	return result;
}

/* apply function */
static int applyx (lua_State *L) {
	luaL_checktype(L, 2, LUA_TFUNCTION);
	lua_settop(L, 2);
	TL = L;
	return apply(L, _apply, 0);
}

/* invokes the GEMV subprogram (y <- alpha A x + b y) */
static int gemv (lua_State *L) {
	struct matrix *A;
	struct vector *x, *y;
	float alpha, beta;
	CBLAS_TRANSPOSE ta;
	int m, n;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	x = luaL_checkudata(L, 2, LUALINEAR_VECTOR_METATABLE);
	y = luaL_checkudata(L, 3, LUALINEAR_VECTOR_METATABLE);
	alpha = luaL_optnumber(L, 4, 1.0);
	beta = luaL_optnumber(L, 5, 0.0);
	ta = checktranspose(L, 6);
	m = ta == CblasNoTrans ? A->rows : A->cols;
	n = ta == CblasNoTrans ? A->cols : A->rows;
	luaL_argcheck(L, x->size == n, 2, "dimension mismatch");
	luaL_argcheck(L, y->size == m, 3, "dimension mismatch");

	/* invoke subprogram */
	cblas_sgemv(A->order, ta, A->rows, A->cols, alpha, A->values, A->ld,
			x->values, x->inc, beta, y->values, y->inc);
	return 0;
}

/* invokes the GER subprogram (A <- alpha x y' + A) */
static int ger (lua_State *L) {
	struct vector *x, *y;
	struct matrix *A;
	float alpha;

	/* check and process arguments */
	x = luaL_checkudata(L, 1, LUALINEAR_VECTOR_METATABLE);
	y = luaL_checkudata(L, 2, LUALINEAR_VECTOR_METATABLE);
	A = luaL_checkudata(L, 3, LUALINEAR_MATRIX_METATABLE);
	alpha = luaL_optnumber(L, 4, 1.0);
	luaL_argcheck(L, x->size == A->rows, 1, "dimension mismatch");
	luaL_argcheck(L, y->size == A->cols, 2, "dimension mismatch");

	/* invoke subprogram */
	cblas_sger(A->order, A->rows, A->cols, alpha, x->values, x->inc,
			y->values, y->inc, A->values, A->ld);
	return 0;
}

/* invokes the GEMM subprogram (C <- alpha A B + beta C) */
static int gemm (lua_State *L) {
	struct matrix *A, *B, *C;
	float alpha, beta;
	CBLAS_TRANSPOSE ta, tb;
	int m, n, ka, kb;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	B = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, B->order == A->order, 2, "order mismatch");
	C = luaL_checkudata(L, 3, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, C->order == A->order, 3, "order mismatch");
	alpha = luaL_optnumber(L, 4, 1.0);
	beta = luaL_optnumber(L, 5, 0.0);
	ta = checktranspose(L, 6);
	tb = checktranspose(L, 7);
	m = ta == CblasNoTrans ? A->rows : A->cols;
	n = tb == CblasNoTrans ? B->cols : B->rows;
	ka = ta == CblasNoTrans ? A->cols : A->rows;
	kb = tb == CblasNoTrans ? B->rows : B->cols;
	luaL_argcheck(L, ka == kb, 2, "dimension mismatch");

	/* invoke subprogramm */
	cblas_sgemm(A->order, ta, tb, m, n, ka, alpha, A->values, A->ld,
			B->values, B->ld, beta, C->values, C->ld);
	return 0;
}

/* invokes the GESV subprogram */
static int gesv (lua_State *L) {
	struct matrix *A, *B;
	int *ipiv, result;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, A->rows == A->cols, 1, "not square");
	B = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, B->order == A->order, 2, "order mismatch");
	luaL_argcheck(L, B->rows == A->rows, 2, "dimension mismatch");

	/* invoke subprogramm */
	ipiv = calloc(A->rows, sizeof(lapack_int));
	if (ipiv == NULL) {
		return luaL_error(L, "cannot allocate indexes");
	}
	result = LAPACKE_sgesv(A->order, A->rows, B->cols, A->values, A->ld,
			ipiv, B->values, B->ld);
	free(ipiv);
	lua_pushinteger(L, result);
	return 1;
}

/* invokes the GELS subprogram */
static int gels (lua_State *L) {
	struct matrix *A, *B;
	char ta;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	B = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, B->order == A->order, 2, "order mismatch");
	ta = lapacktranspose(checktranspose(L, 3));
	luaL_argcheck(L, B->rows == (A->rows >= A->cols ? A->rows : A->cols), 2,
			"dimension mismatch");

	/* invoke subprogramm */
	lua_pushinteger(L, LAPACKE_sgels(A->order, ta, A->rows, A->cols,
			B->cols, A->values, A->ld, B->values, B->ld));
	return 1;
}

/* calculates the inverse of a matrix */
static int inv (lua_State *L) {
	struct matrix *A;
	int *ipiv, result;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, A->rows == A->cols, 1, "not square");

	/* invoke subprograms */
	ipiv = calloc(A->rows, sizeof(lapack_int));
	if (ipiv == NULL) {
		return luaL_error(L, "cannot allocate indexes");
	}
	result = LAPACKE_sgetrf(A->order, A->rows, A->cols, A->values, A->ld,
			ipiv);
	if (result != 0) {
		free(ipiv);
		lua_pushinteger(L, result);
		return 1;
	}
	result = LAPACKE_sgetri(A->order, A->rows, A->values, A->ld, ipiv);
	free(ipiv);
	lua_pushinteger(L, result);
	return 1;
}

/* calculates the determinant of a matrix */
static int det (lua_State *L) {
	struct matrix *A;
	float *copy, *d, *s, det;
	int n, *ipiv, result, neg, i;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, A->rows == A->cols, 1, "not square");
	n = A->rows;

	/* copy matrix */
	copy = calloc((size_t)n * n, sizeof(float));
	if (copy == NULL) {
		return luaL_error(L, "cannot allocate values");
	}
	d = copy;
	s = A->values;
	for (i = 0; i < n; i++) {
		memcpy(d, s, (size_t)n * sizeof(float));
		d += n;
		s += A->ld;
	}

	/* invoke subprograms */
	ipiv = calloc(n, sizeof(lapack_int));
	if (ipiv == NULL) {
		free(copy);
		return luaL_error(L, "cannot allocate indexes");
	}
	result = LAPACKE_sgetrf(A->order, n, n, copy, n, ipiv);
	if (result != 0) {
		free(copy);
		free(ipiv);
		lua_pushnumber(L, 0.0);
		return 1;
	}

	/* calculate determinant */
	det = 1.0;
	neg = 0;
	for (i = 0; i < n; i++) {
		det *= copy[(size_t)i * n + i];
		if (ipiv[i] != i + 1) {
			neg = !neg;
		}
	}
	free(copy);
	free(ipiv);
	lua_pushnumber(L, neg ? -det : det);
	return 1;
}

/* calculates the covariance of a matrix */
static int cov (lua_State *L) {
	struct matrix *A, *B;
	int ddof, i, j, k;
	float *means, *v, *vi, *vj, sum;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	B = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, A->cols == B->rows, 2, "dimension mismatch");
	luaL_argcheck(L, B->rows == B->cols, 2, "not square");
	ddof = luaL_optinteger(L, 3, 0);
	luaL_argcheck(L, ddof >= 0 && ddof < A->rows, 3, "bad ddof");

	/* calculate means */
	means = calloc((size_t)A->cols, sizeof(float));
	if (means == NULL) {
		return luaL_error(L, "cannot allocate values");
	}
	switch (A->order) {
	case CblasRowMajor:
		#pragma omp parallel for private(i, j, sum, v) schedule(auto) \
				if(A->rows * A->cols >= LUALINEAR_OMP_MINSIZE)
		for (i = 0; i < A->cols; i++) {
			sum = 0.0;
			v = &A->values[i];
			for (j = 0; j < A->rows; j++) {
				sum += *v;
				v += A->ld;
			}
			means[i] = sum / A->rows;
		}
		break;

	case CblasColMajor:
		#pragma omp parallel for private(i, j, sum, v) schedule(auto) \
				if(A->rows * A->cols >= LUALINEAR_OMP_MINSIZE)
		for (i = 0; i < A->cols; i++) {
			sum = 0.0;
			v = &A->values[(size_t)i * A->ld];
			for (j = 0; j < A->rows; j++) {
				sum += *v;
				v++;
			}
			means[i] = sum / A->rows;
		}
		break;
	}

	/* calculate covariance */
	switch (A->order) {
	case CblasRowMajor:
		for (i = 0; i < A->cols; i++) {
			#pragma omp parallel for private(j, k, sum, vi, vj) \
					schedule(auto) if(A->rows * (A->cols \
					- i) >= LUALINEAR_OMP_MINSIZE)
			for (j = i; j < A->cols; j++) {
				sum = 0.0;
				vi = &A->values[i];
				vj = &A->values[j];
				for (k = 0; k < A->rows; k++) {
					sum += (*vi - means[i])
							* (*vj - means[j]);
					vi += A->ld;
					vj += A->ld;
				}
				B->values[(size_t)i * B->ld + j] = B->values[
						(size_t)j * B->ld + i]
						= sum / (A->rows - ddof);
			}
		}
		break;

	case CblasColMajor:
		for (i = 0; i < A->cols; i++) {
			#pragma omp parallel for private(j, k, sum, vi, vj) \
					schedule(auto) if(A->rows * (A->cols \
					- i) >= LUALINEAR_OMP_MINSIZE)
			for (j = i; j < A->cols; j++) {
				sum = 0.0;
				vi = &A->values[(size_t)i * A->ld];
				vj = &A->values[(size_t)j * A->ld];
				for (k = 0; k < A->rows; k++) {
					sum += (*vi - means[i])
							* (*vj - means[j]);
					vi++;
					vj++;
				}
				B->values[(size_t)i * B->ld + j] = B->values[
						(size_t)j * B->ld + i]
						= sum / (A->rows - ddof);
			}
		}
		break;
	}
	free(means);
	return 0;
}

/* calculates the correlation of a matrix */
static int corr (lua_State *L) {
	struct matrix *A, *B;
	int i, j, k;
	float *means, *stds, *v, *vi, *vj, sum;

	/* check and process arguments */
	A = luaL_checkudata(L, 1, LUALINEAR_MATRIX_METATABLE);
	B = luaL_checkudata(L, 2, LUALINEAR_MATRIX_METATABLE);
	luaL_argcheck(L, A->cols == B->rows, 2, "dimension mismatch");
	luaL_argcheck(L, B->rows == B->cols, 2, "not square");

	/* calculate means and stds */
	means = calloc((size_t)A->cols, sizeof(float));
	if (means == NULL) {
		return luaL_error(L, "cannot allocate values");
	}
	stds = calloc((size_t)A->cols, sizeof(float));
	if (stds == NULL) {
		free(means);
		return luaL_error(L, "cannot allocate values");
	}
	switch (A->order) {
	case CblasRowMajor:
		#pragma omp parallel for private(i, j, sum, v) schedule(auto) \
				if(A->rows * A->cols >= LUALINEAR_OMP_MINSIZE)
		for (i = 0; i < A->cols; i++) {
			sum = 0.0;
			v = &A->values[i];
			for (j = 0; j < A->rows; j++) {
				sum += *v;
				v += A->ld;
			}
			means[i] = sum / A->rows;
			sum = 0.0;
			v = &A->values[i];
			for (j = 0; j < A->rows; j++) {
				sum += (*v - means[i]) * (*v - means[i]);
				v += A->ld;
			}
			stds[i] = sqrt(sum);
		}
		break;

	case CblasColMajor:
		#pragma omp parallel for private(i, j, sum, v) schedule(auto) \
				if(A->rows * A->cols >= LUALINEAR_OMP_MINSIZE)
		for (i = 0; i < A->cols; i++) {
			sum = 0.0;
			v = &A->values[(size_t)i * A->ld];
			for (j = 0; j < A->rows; j++) {
				sum += *v;
				v++;
			}
			means[i] = sum / A->rows;
			sum = 0.0;
			v = &A->values[(size_t)i * A->ld];
			for (j = 0; j < A->rows; j++) {
				sum += (*v - means[i]) * (*v - means[i]);
				v++;
			}
			stds[i] = sqrt(sum);
		}
		break;
	}

	/* calculate correlation */
	switch (A->order) {
	case CblasRowMajor:
		for (i = 0; i < A->cols; i++) {
			#pragma omp parallel for private(j, k, sum, vi, vj) \
					schedule(auto) if(A->rows * (A->cols \
					- i) >= LUALINEAR_OMP_MINSIZE)
			for (j = i; j < A->cols; j++) {
				sum = 0.0;
				vi = &A->values[i];
				vj = &A->values[j];
				for (k = 0; k < A->rows; k++) {
					sum += (*vi - means[i])
							* (*vj - means[j]);
					vi += A->ld;
					vj += A->ld;
				}
				B->values[(size_t)i * B->ld + j] = B->values[
						(size_t)j * B->ld + i]
						= sum / (stds[i] * stds[j]);
			}
		}
		break;

	case CblasColMajor:
		for (i = 0; i < A->cols; i++) {
			#pragma omp parallel for private(j, k, sum, vi, vj) \
					schedule(auto) if(A->rows * (A->cols \
					- i) >= LUALINEAR_OMP_MINSIZE)
			for (j = i; j < A->cols; j++) {
				sum = 0.0;
				vi = &A->values[(size_t)i * A->ld];
				vj = &A->values[(size_t)j * A->ld];
				for (k = 0; k < A->rows; k++) {
					sum += (*vi - means[i])
							* (*vj - means[j]);
					vi++;
					vj++;
				}
				B->values[(size_t)i * B->ld + j] = B->values[
						(size_t)j * B->ld + i]
						= sum / (stds[i] * stds[j]);
			}
		}
		break;
	}
	free(means);
	free(stds);
	return 0;
}

/*
 * Exported functions.
 */

int luaopen_linear (lua_State *L) {
	static const luaL_Reg FUNCTIONS[] = {
		{ "vector", vector },
		{ "matrix", matrix },
		{ "type", type },
		{ "size", size },
		{ "tvector", tvector },
		{ "sub", sub },
		{ "unwind", unwind },
		{ "reshape", reshape },
		{ "totable", totable },
		{ "tolinear", tolinear },
		{ "dot", dot },
		{ "nrm2", nrm2 },
		{ "asum", asum },
		{ "iamax", iamax },
		{ "sum", sum },
		{ "swap", swap },
		{ "copy", copy },
		{ "axpy", axpy },
		{ "scal", scal },
		{ "set", set },
		{ "uniform", uniform },
		{ "normal", normal },
		{ "inc", inc },
		{ "mul", mul },
		{ "pow", powx },
		{ "sign", sign },
		{ "abs", absx },
		{ "exp", expx },
		{ "log", logx },
		{ "logistic", logistic },
		{ "tanh", tanhx },
		{ "softplus", softplus },
		{ "rectifier", rectifier },
		{ "apply", applyx },
		{ "gemv", gemv },
		{ "ger", ger },
		{ "gemm", gemm },
		{ "gesv", gesv },
		{ "gels", gels },
		{ "inv", inv },
		{ "det", det },
		{ "cov", cov },
		{ "corr", corr },
		{ NULL, NULL }
	};

	/* register functions */
	#if LUA_VERSION_NUM >= 502
	luaL_newlib(L, FUNCTIONS);
	#else
	luaL_register(L, luaL_checkstring(L, 1), FUNCTIONS);
	#endif

	/* vector metatable */
	luaL_newmetatable(L, LUALINEAR_VECTOR_METATABLE);
	lua_pushcfunction(L, vector_len);
	lua_setfield(L, -2, "__len");
	lua_pushcfunction(L, vector_index);
	lua_setfield(L, -2, "__index");
	lua_pushcfunction(L, vector_newindex);
	lua_setfield(L, -2, "__newindex");
	lua_pushcfunction(L, vector_ipairs);
	lua_setfield(L, -2, "__ipairs");
	lua_pushcfunction(L, vector_tostring);
	lua_setfield(L, -2, "__tostring");
	lua_pushcfunction(L, vector_free);
	lua_setfield(L, -2, "__gc");
	lua_pop(L, 1);

	/* matrix metatable */
	luaL_newmetatable(L, LUALINEAR_MATRIX_METATABLE);
	lua_pushcfunction(L, matrix_len);
	lua_setfield(L, -2, "__len");
	lua_pushcfunction(L, matrix_index);
	lua_setfield(L, -2, "__index");
	lua_pushcfunction(L, matrix_ipairs);
	lua_setfield(L, -2, "__ipairs");
	lua_pushcfunction(L, matrix_tostring);
	lua_setfield(L, -2, "__tostring");
	lua_pushcfunction(L, matrix_free);
	lua_setfield(L, -2, "__gc");
	lua_pop(L, 1);

	return 1;
}

max_stars_repo_name: anaef/lua-linear

max_issues_repo_name: anaef/lua-linear

max_forks_repo_name: anaef/lua-linear


text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "parmt_utils.h"
#ifdef PARMT_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "iscl/memory/memory.h"

/*!
 * @brief Converts the NED based fundamental faults to Green's functions which
 *        scale the NED based moment tensor terms s.t. \f$ u = G m \f$.
 *        where \f$ m = \{ m_{xx}, m_{yy}, m_{zz}, m_{xy}, m_{xz}, m_{yz} \} \f$
 *        are the NED moment tensor terms and \f$ u \f$ is the estimate for
 *        the icomp'th component with source receiver azimuth, az, and given
 *        channel orientation encapsulated by cmpaz and cmpinc.
 *
 * @param[in] npgrns     number of points in greens functions
 * @param[in] ldg        leading dimension of G (>= 6)
 * @param[in] icomp      =1 for vertical Green's functions.
 *                       =2 for north (channel 2) Green's functions.
 *                       =3 for east (channel 3) Green's functions.
 * @param[in] azin       source to receiver azimuth (degrees)
 * @param[in] bazin      receiver to source azimuth (degrees).  if you
 *                       do not know then set this to:
 *                         fmod(az + 180.0, 360.0)
 * @param[in] cmpaz      component azimuth (0 north, +90 east)
 * @param[in] cmpinc     component inclinantion (-90 up, 0 east/north, +90 down)
 * @param[in] ZDS        vertical greens fn for 90 degree dip slip [npgrns]
 * @param[in] ZSS        vertical greens fn for vertical strike slip [npgrns]
 * @param[in] ZDD        vertical greens fn for 45 degree dip slip [npgrns]
 * @param[in] ZEX        vertical greens fn for an explosion [npgrns]
 * @param[in] RDS        radial greens fn for 90 degree dip slip [npgrns]
 * @param[in] RSS        radial greens fn for vertical strike slip [npgrns]
 * @param[in] RDD        radial greens fn for 45 degree dip slip [npgrns]
 * @param[in] REX        radial greens fn for an explosion [npgrns]
 * @param[in] TDS        transverse greens fn for 90 degree dip slip [npgrns]
 * @param[in] TSS        transverse greens fn for vertical strike slip [npgrns]
 *
 * @param[out] G         row major Green's functions matrix [npgrns x ldg].
 *                       the columns are ordered:
 *                        \f$ \{ G_{xx}, G_{yy}, G_{zz},
 *                               G_{xy}, G_{xz}, G_{yz} \}  \f$
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_utils_ff2mtGreensMatrix64f(const int npgrns, const int ldg,
                                     const int icomp,
                                     const double azin, const double bazin,
                                     const double cmpaz, const double cmpinc,
                                     const double *__restrict__ ZDS,
                                     const double *__restrict__ ZSS,
                                     const double *__restrict__ ZDD,
                                     const double *__restrict__ ZEX,
                                     const double *__restrict__ RDS,
                                     const double *__restrict__ RSS,
                                     const double *__restrict__ RDD,
                                     const double *__restrict__ REX,
                                     const double *__restrict__ TDS,
                                     const double *__restrict__ TSS,
                                     double *__restrict__ G)
{
    double *Gxx, *Gyy, *Gzz, *Gxy, *Gxz, *Gyz;
    int ierr;
    ierr = 0;
    if (npgrns < 1 || ldg < 6 || G == NULL)
    {
        if (npgrns < 1)
        {
            fprintf(stderr, "%s: No points in greens fns\n", __func__);
        }
        if (ldg < 6){fprintf(stderr, "%s: ldg must be >= 6\n", __func__);}
        if (G == NULL){fprintf(stderr, "%s: Error G is NULL\n", __func__);}
        return -1;
    }
    Gxx = memory_calloc64f(npgrns);
    Gyy = memory_calloc64f(npgrns);
    Gzz = memory_calloc64f(npgrns);
    Gxy = memory_calloc64f(npgrns);
    Gxz = memory_calloc64f(npgrns);
    Gyz = memory_calloc64f(npgrns);
    ierr = parmt_utils_ff2mtGreens64f(npgrns, icomp, azin, bazin,
                                      cmpaz, cmpinc,
                                      ZDS, ZSS, ZDD, ZEX,
                                      RDS, RSS, RDD, REX,
                                      TDS, TSS,
                                      Gxx, Gyy, Gzz,
                                      Gxy, Gxz, Gyz);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to compute greens functions columns\n",
                __func__);
    }
    else
    {
        cblas_dcopy(npgrns, Gxx, 1, &G[0], ldg);
        cblas_dcopy(npgrns, Gyy, 1, &G[1], ldg);
        cblas_dcopy(npgrns, Gzz, 1, &G[2], ldg);
        cblas_dcopy(npgrns, Gxy, 1, &G[3], ldg);
        cblas_dcopy(npgrns, Gxz, 1, &G[4], ldg);
        cblas_dcopy(npgrns, Gyz, 1, &G[5], ldg);
    }
    memory_free64f(&Gxx);
    memory_free64f(&Gyy);
    memory_free64f(&Gzz);
    memory_free64f(&Gxy);
    memory_free64f(&Gxz);
    memory_free64f(&Gyz);
    return ierr;
}
//============================================================================//
/*!
 * @brief Converts the NED based fundamental faults to Green's functions which
 *        scale the NED based moment tensor terms s.t. \f$ u = G m \f$.
 *        where \f$ m = \{ m_{xx}, m_{yy}, m_{zz}, m_{xy}, m_{xz}, m_{yz} \} \f$
 *        are the NED moment tensor terms and \f$ u \f$ is the estimate for
 *        the icomp'th component with source receiver azimuth, az, and given
 *        channel orientation encapsulated by cmpaz and cmpinc.
 *
 * @param[in] npgrns     number of points in greens functions
 * @param[in] icomp      =1 for vertical Green's functions.
 *                       =2 for north (channel 2) Green's functions.
 *                       =3 for east (channel 3) Green's functions.
 * @param[in] azin       source to receiver azimuth (degrees)
 * @param[in] bazin      receiver to source azimuth (degrees).  if you
 *                       do not know then set this to:
 *                         fmod(az + 180.0, 360.0)
 * @param[in] cmpaz      component azimuth (0 north, +90 east)
 * @param[in] cmpinc     component inclinantion (-90 up, 0 east/north, +90 down)
 * @param[in] ZDS        vertical greens fn for 90 degree dip slip [npgrns]
 * @param[in] ZSS        vertical greens fn for vertical strike slip [npgrns]
 * @param[in] ZDD        vertical greens fn for 45 degree dip slip [npgrns]
 * @param[in] ZEX        vertical greens fn for an explosion [npgrns]
 * @param[in] RDS        radial greens fn for 90 degree dip slip [npgrns]
 * @param[in] RSS        radial greens fn for vertical strike slip [npgrns]
 * @param[in] RDD        radial greens fn for 45 degree dip slip [npgrns]
 * @param[in] REX        radial greens fn for an explosion [npgrns]
 * @param[in] TDS        transverse greens fn for 90 degree dip slip [npgrns]
 * @param[in] TSS        transverse greens fn for vertical strike slip [npgrns]
 *
 * @param[out] Gxx       greens function scaling mxx moment tensor term [npgrns]
 * @param[out] Gyy       greens function scaling myy moment tensor term [npgrns]
 * @param[out] Gzz       greens function scaling mzz moment tensor term [npgrns]
 * @param[out] Gxy       greens function scaling mxy moment tensor term [npgrns]
 * @param[out] Gxz       greens function scaling mxz moment tensor term [npgrns]
 * @param[out] Gyz       greens function scaling myz moment tensor term [npgrns]
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_utils_ff2mtGreens64f(const int npgrns, const int icomp,
                               const double azin, const double bazin,
                               const double cmpaz, const double cmpinc,
                               const double *__restrict__ ZDS,
                               const double *__restrict__ ZSS,
                               const double *__restrict__ ZDD,
                               const double *__restrict__ ZEX,
                               const double *__restrict__ RDS,
                               const double *__restrict__ RSS,
                               const double *__restrict__ RDD,
                               const double *__restrict__ REX,
                               const double *__restrict__ TDS,
                               const double *__restrict__ TSS,
                               double *__restrict__ Gxx,
                               double *__restrict__ Gyy,
                               double *__restrict__ Gzz,
                               double *__restrict__ Gxy,
                               double *__restrict__ Gxz,
                               double *__restrict__ Gyz)
{
    double az, baz, c2a, ca, caz, cos_caz, cost,
           gxx_e, gxy_e, gxz_e, gyy_e, gyz_e, gzz_e,
           gxx_n, gxy_n, gxz_n, gyy_n, gyz_n, gzz_n,
           gxx_r, gxy_r, gxz_r, gyy_r, gyz_r, gzz_r,
           gxx_t, gxy_t, gxz_t, gyy_t, gyz_t, gzz_t,
           gxx_z, gxy_z, gxz_z, gyy_z, gyz_z, gzz_z,
           rex, rdd, rds, rss, s2a, sa, sin_caz, sint,
           tds, tss, xscal,
           zex, zdd, zds, zss;
    int i;
    const double pi180 = M_PI/180.0;
    const double half = 0.5;
    const double sixth = 1.0/6.0;
    const double third = 1.0/3.0;
    //------------------------------------------------------------------------//
    if (icomp < 1 || icomp > 3)
    {
        fprintf(stderr, "%s: Invalid component: %d\n", __func__, icomp);
        return -1;
    }
    if (icomp == 1)
    {
        if (ZDS == NULL || ZSS == NULL || ZDD == NULL || ZEX == NULL)
        {
            if (ZDS == NULL){fprintf(stderr, "%s: zds is NULL\n", __func__);}
            if (ZSS == NULL){fprintf(stderr, "%s: zss is NULL\n", __func__);}
            if (ZDD == NULL){fprintf(stderr, "%s: zdd is NULL\n", __func__);}
            if (ZEX == NULL){fprintf(stderr, "%s: zex is NULL\n", __func__);}
            return -1;
        }
    }
    else
    {
        if (RDS == NULL || RSS == NULL || RDD == NULL || REX == NULL ||
            TDS == NULL || TSS == NULL)
        {
            if (RDS == NULL){fprintf(stderr, "%s: rds is NULL\n", __func__);}
            if (RSS == NULL){fprintf(stderr, "%s: rss is NULL\n", __func__);}
            if (RDD == NULL){fprintf(stderr, "%s: rdd is NULL\n", __func__);}
            if (REX == NULL){fprintf(stderr, "%s: rex is NULL\n", __func__);}
            if (TDS == NULL){fprintf(stderr, "%s: tds is NULL\n", __func__);}
            if (TSS == NULL){fprintf(stderr, "%s: tss is NULL\n", __func__);}
            return -1;
        }
    }
    // Get the backazimuth
    az = azin;
    baz = bazin; //fmod(az + 180.0, 360.0);
    // Convert to radians
    az = az*pi180;
    baz = baz*pi180;
    caz = cmpaz*pi180;
    // Compute the geometric factors
    sa = sin(az);
    ca = cos(az);
    s2a = sin(2.*az);
    c2a = cos(2.*az);
    sint = sin(baz);
    cost = cos(baz);
    cos_caz = cos(caz);
    sin_caz = sin(caz);
    // the cmpaz will come from the metadata so fix that
    if (icomp == 3)
    {
       cos_caz = cos(caz - M_PI/2.0);
       sin_caz = sin(caz - M_PI/2.0);
    }
    // Set the columns (Minson et. al. 2008)
    xscal = 1.0;
    if (fabs(cmpinc - 90.0) < 1.e-4){xscal =-1.0;}
    if (icomp == 1)
    {
        #pragma omp simd
        for (i=0; i<npgrns; i++)
        {
            zss = ZSS[i];
            zdd = ZDD[i];
            zds = ZDS[i];
            zex = ZEX[i];
            // compute the greens function corresponding to the mt
            gxx_z = half*(zss*c2a) - sixth*zdd + third*zex;
            gyy_z =-half*(zss*c2a) - sixth*zdd + third*zex;
            gzz_z = third*(zdd + zex);
            gxy_z = zss*s2a;
            gxz_z = zds*ca;
            gyz_z = zds*sa;
            // copy and include polarity so it matches the observation
            Gxx[i] = xscal*gxx_z;
            Gyy[i] = xscal*gyy_z;
            Gzz[i] = xscal*gzz_z;
            Gxy[i] = xscal*gxy_z;
            Gxz[i] = xscal*gxz_z;
            Gyz[i] = xscal*gyz_z;
        } // Loop on points
    }
    else
    {
        // North or 2 component
        if (icomp == 2)
        {
            // Loop on rows (data points)
            //#pragma omp simd
            for (i=0; i<npgrns; i++)
            {
                rss = RSS[i];
                rdd = RDD[i];
                rds = RDS[i];
                rex = REX[i];
                tds = TDS[i];
                tss = TSS[i];

                gxx_r = half*(rss*c2a) - sixth*rdd + third*rex;
                gyy_r =-half*(rss*c2a) - sixth*rdd + third*rex;
                gzz_r = third*(rdd + rex);
                gxy_r = rss*s2a;
                gxz_r = rds*ca;
                gyz_r = rds*sa;

                gxx_t = half*tss*s2a;
                gyy_t =-half*tss*s2a;
                gzz_t = 0.0;
                gxy_t =-tss*c2a;
                gxz_t = tds*sa;
                gyz_t =-tds*ca;

                // Rotate from (r, t) -> (n, e)
                //   n = cos(baz-180)*r - sin(baz-180)*t
                //     =-cos(baz)*r + sin(baz)*t
                //   e = sin(baz-180)*r + cos(baz-180)*t
                //     =-sin(baz)*r - cos(baz)*t
                gxx_n =-cost*gxx_r + sint*gxx_t;
                gxx_e =-sint*gxx_r - cost*gxx_t;
                gyy_n =-cost*gyy_r + sint*gyy_t;
                gyy_e =-sint*gyy_r - cost*gyy_t;
                gzz_n =-cost*gzz_r + sint*gzz_t;
                gzz_e =-sint*gzz_r - cost*gzz_t;
                gxy_n =-cost*gxy_r + sint*gxy_t;
                gxy_e =-sint*gxy_r - cost*gxy_t;
                gxz_n =-cost*gxz_r + sint*gxz_t;
                gxz_e =-sint*gxz_r - cost*gxz_t;
                gyz_n =-cost*gyz_r + sint*gyz_t;
                gyz_e =-sint*gyz_r - cost*gyz_t;
                // Rotate from (n, e) -> (1, 2) around az
                //   1 = cos(comp_az)*n + sin(comp_az)*e
                Gxx[i] = cos_caz*gxx_n + sin_caz*gxx_e;
                Gyy[i] = cos_caz*gyy_n + sin_caz*gyy_e;
                Gzz[i] = cos_caz*gzz_n + sin_caz*gzz_e;
                Gxy[i] = cos_caz*gxy_n + sin_caz*gxy_e;
                Gxz[i] = cos_caz*gxz_n + sin_caz*gxz_e;
                Gyz[i] = cos_caz*gyz_n + sin_caz*gyz_e;

            } // Loop on points
        }
        // East or 3 component
        else
        {
            // Loop on rows (data points)
            #pragma omp simd
            for (i=0; i<npgrns; i++)
            {
                rss = RSS[i];
                rdd = RDD[i];
                rds = RDS[i];
                rex = REX[i];
                tds = TDS[i];
                tss = TSS[i];

                gxx_r = half*(rss*c2a) - sixth*rdd + third*rex;
                gyy_r =-half*(rss*c2a) - sixth*rdd + third*rex;
                gzz_r = third*(rdd + rex);
                gxy_r = rss*s2a;
                gxz_r = rds*ca;
                gyz_r = rds*sa;

                gxx_t = half*tss*s2a;
                gyy_t =-half*tss*s2a;
                gzz_t = 0.0;
                gxy_t =-tss*c2a;
                gxz_t = tds*sa;
                gyz_t =-tds*ca;

                // Rotate from (r, t) -> (n, e)
                //   n = cos(baz-180)*r - sin(baz-180)*t
                //     =-cos(baz)*r + sin(baz)*t
                //   e = sin(baz-180)*r + cos(baz-180)*t
                //     =-sin(baz)*r - cos(baz)*t
                gxx_n =-cost*gxx_r + sint*gxx_t;
                gxx_e =-sint*gxx_r - cost*gxx_t;
                gyy_n =-cost*gyy_r + sint*gyy_t;
                gyy_e =-sint*gyy_r - cost*gyy_t;
                gzz_n =-cost*gzz_r + sint*gzz_t;
                gzz_e =-sint*gzz_r - cost*gzz_t;
                gxy_n =-cost*gxy_r + sint*gxy_t;
                gxy_e =-sint*gxy_r - cost*gxy_t;
                gxz_n =-cost*gxz_r + sint*gxz_t;
                gxz_e =-sint*gxz_r - cost*gxz_t;
                gyz_n =-cost*gyz_r + sint*gyz_t;
                gyz_e =-sint*gyz_r - cost*gyz_t;
                // Rotate from (n, e) -> (1, 2) around az
                //   2 =-sin(comp_az)*n + cos(comp_az)*e
                Gxx[i] =-sin_caz*gxx_n + cos_caz*gxx_e;
                Gyy[i] =-sin_caz*gyy_n + cos_caz*gyy_e;
                Gzz[i] =-sin_caz*gzz_n + cos_caz*gzz_e;
                Gxy[i] =-sin_caz*gxy_n + cos_caz*gxy_e;
                Gxz[i] =-sin_caz*gxz_n + cos_caz*gxz_e;
                Gyz[i] =-sin_caz*gyz_n + cos_caz*gyz_e;
            } // Loop on points
        } // End check on 1 or 2 component
    } // End check on component
    return 0;
}


text:
#pragma once

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

#include <tiny_dnn/tiny_dnn.h>

#include "defines.h"
#include "matrix_vector_ops.h"
#include "genann.h"

#define NUM_STATES 12
#define NUM_INPUTS 4

// Dynamics model from G. Gremillion, S. Humbert paper "System Identification of
// a Quadrotor Micro Air Vehicle" (equation 3 in paper)
// Augmented to include position dynamics in hover
// Dynamics and controls matrix:
// cos(ps)*u  -sin(ps)*v    0     0     0     0     0     0     0     0     0     0     0     0     0     0
// sin(ps)*u  cos(ps)*v     0     0     0     0     0     0     0     0     0     0     0     0     0     0
//  0     0     1     0     0     0     0     0     0     0     0     0     0     0     0     0
//  0     0     0   X_u*u   0     0     0     0     0     0  X_th*th  0     0     0     0     0
//  0     0     0     0   Y_v*v   0     0     0     0  Y_ph*ph  0     0     0     0     0     0
//  0     0     0     0     0   Z_w*w   0     0     0     0     0     0     0     0     0 Z_thr*d_thr
//  0     0     0     0     0     0   L_p*p   0     0  L_ph*ph  0     0 L_la*d_la 0     0     0
//  0     0     0     0     0     0     0   M_q*q   0     0   M_th*th 0     0 M_lo*d_lo 0     0
//  0     0     0     0     0     0     0     0   N_r*r   0     0     0     0     0 N_ya*d_ya 0
//  0     0     0     0     0     0   ph_p*p  0     0     0     0     0 ph_la*d_la 0    0     0
//  0     0     0     0     0     0     0   th_q*q  0     0     0     0     0 th_lo*d_lo 0    0
//  0     0     0     0     0     0     0     0   ps_r*r  0     0     0     0     0 ps_ya*d_ya 0
//
//  Parameters:
//  X_u   =  -0.27996
//  Y_v   =  -0.22566
//  Z_w   =  -1.2991
//  L_p   =  -2.5110
//  M_q   =  -2.4467
//  N_r   =  -0.4948
//  X_th  =  -10.067
//  Y_ph  =   9.8648
//  L_ph  =  -21.358
//  M_th  =  -18.664
//  ph_p  =   0.9655
//  th_q  =   0.9634
//  ps_r  =   0.6748
//  Z_thr =  -39.282
//  L_la  =   11.468
//  M_lo  =   9.5711
//  N_ya  =   3.5647
//  ph_la =   0.0744
//  th_lo =   0.0594
//  ps_ya =   0.0397
//
// State vector:
//    x y z u v w p q r phi theta psi
//
// Control vector:
//    del_lat del_lon del_yaw del_thrust
int dynamics(gsl_vector *dy, double t, const gsl_vector *y, const gsl_vector *u);
int setupDynamics();
int teardownDynamics();

gsl_vector *feedback(gsl_vector *yd, gsl_vector *y);
int setupFeedback();
int teardownFeedback();

gsl_vector *nnFeedback(gsl_vector *yd, gsl_vector *y);
int setupNnFeedback(int dataShape[2]);
int setupNnFeedback(FILE *in);
int teardownNnFeedback();
genann *getFeedbackNn();

max_stars_repo_name: umd-agrc/SimpleControlSim

max_issues_repo_name: umd-agrc/SimpleControlSim

max_forks_repo_name: umd-agrc/SimpleControlSim