text:
/*
** gsl routines for float data type
**
** G.Lohmann, July 2004
*/
#include <viaio/Vlib.h>

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

#include <gsl/gsl_cblas.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_sort_vector.h>
#include "gsl_utils.h"

/*
** y = Ax
*/
gsl_vector_float *fmat_x_vector(gsl_matrix_float *A, gsl_vector_float *x, gsl_vector_float *y)
{
  if (y == NULL) {
    y = gsl_vector_float_alloc (A->size1);
  }

  gsl_blas_sgemv(CblasNoTrans, 1.0, A, x, 0.0, y);

  return y;
}

/*
**    C = A x B^T
*/
gsl_matrix_float *fmat_x_matT(gsl_matrix_float *A,gsl_matrix_float *B,gsl_matrix_float *C)
{
  if (C == NULL) {
    C = gsl_matrix_float_alloc( A->size1, B->size1 );
  }

  gsl_blas_sgemm( CblasNoTrans, CblasTrans, 1.0, A, B, 0.0, C );

  return C;
}

/*
**    C = A^T x B
*/
gsl_matrix_float *fmatT_x_mat(gsl_matrix_float *A,gsl_matrix_float *B,gsl_matrix_float *C)
{
  if (C == NULL) {
    C = gsl_matrix_float_alloc( A->size2, B->size2 );
  }

  gsl_blas_sgemm( CblasTrans, CblasNoTrans, 1.0, A, B, 0.0, C );

  return C;
}

/*
**    C = A x B
*/
gsl_matrix_float *fmat_x_mat(gsl_matrix_float *A,gsl_matrix_float *B,gsl_matrix_float *C)
{
  if (C == NULL) {
    C = gsl_matrix_float_alloc( A->size1, B->size2 );
  }

  gsl_blas_sgemm( CblasNoTrans, CblasNoTrans, 1.0, A, B, 0.0, C );

  return C;
}

/*
**  y = Ax,  A double, x float
*/
gsl_vector_float *dmat_x_fvector(gsl_matrix *A, gsl_vector_float *x, gsl_vector_float *y)
{
  int i,j,nrows,ncols;
  float *ptr2,*ptr3,sum;
  double *ptr1;

  nrows = A->size1;
  ncols = A->size2;

  if (y == NULL) {
     y = gsl_vector_float_alloc (nrows);
  }

  if (x->size != ncols || y->size != nrows) {
    fprintf(stderr," fmat_x_vect: incongruent dimensions\n");
    exit(0);
  }

  ptr1 = A->data;
  ptr3 = y->data;

  for (i=0; i<nrows; i++) {
    sum = 0;
    ptr2 = x->data;
    for (j=0; j<ncols; j++) {
      sum += (*ptr1++) * (*ptr2++);
    }
    *ptr3++ = sum;
  }
  return y;
}

/*
**  z = x^T y
*/
float fskalarproduct(gsl_vector_float *x,gsl_vector_float *y)
{
  int i,n;
  float *ptr1,*ptr2,sum;

  n = x->size;
  if (y->size != n) {
    fprintf(stderr," fskalarproduct: incongruent vector sizes: %d %ld",n,y->size);
    exit(0);
  }

  ptr1 = x->data;
  ptr2 = y->data;
  sum = 0;
  for (i=0; i<n; i++) {
    sum += (*ptr1) * (*ptr2);
    ptr1++;
    ptr2++;
  }
  return sum;
}

/*
**    printout
*/
void fmatprint(FILE *fp,gsl_matrix_float *A,const char *format)
{
  int i,j;
  for (i=0; i<A->size1; i++) {
    for (j=0; j<A->size2; j++) {
      fprintf(fp,format,fmget(A,i,j));
    }
    fprintf(fp,"\n");
  }
  fprintf(fp,"\n");

}

/*
**    B = A^T
*/
gsl_matrix_float *ftranspose(gsl_matrix_float *A,gsl_matrix_float *B)
{
  int i,j,n,m;

  n = A->size1;
  m = A->size2;

  if (B == NULL) {
    B = gsl_matrix_float_alloc(m,n);
  }
  else if (B->size1 != m || B->size2 != n) {
    gsl_matrix_float_free(B);
    B = gsl_matrix_float_alloc(m,n);
  }

  for (i=0; i<n; i++) {
    for (j=0; j<m; j++) {
      fmset(B,j,i,fmget(A,i,j));
    }
  }
  return B;
}

/*
**    B = A^-1 = A^+
*/
gsl_matrix_float *fmat_PseudoInv(gsl_matrix_float *A,gsl_matrix_float *B)
{
  int i,j,k,l,n,m;
  double u,x;
  double *dbl_pp;
  float *flt_pp;

  m = A->size1;
  n = A->size2;

  if (B == NULL) {
    B = gsl_matrix_float_alloc (n,m);
  }
  else if (B->size1 != n || B->size2 != m) {
    gsl_matrix_float_free(B);
    B = gsl_matrix_float_alloc (n, m);
  }

  gsl_matrix *U = gsl_matrix_alloc (m, n);
  gsl_matrix *V = gsl_matrix_alloc (n, n);
  gsl_matrix *X = gsl_matrix_alloc (n, n);
  gsl_vector *w = gsl_vector_alloc (n);
  gsl_vector *work = gsl_vector_alloc (n);

  /* singular value decomposition */
  flt_pp = A->data;
  dbl_pp = U->data;
  for (i=0; i<A->size1 * A->size2; i++)
    *dbl_pp++ = *flt_pp++;

  /* gsl_linalg_SV_decomp_jacobi(U,V,w); */
  gsl_linalg_SV_decomp_mod (U,X,V,w,work);

  /* check if singular */
  k=0;
  for (j=0; j<n; j++) {
    if (fabs(w->data[j]) > 1.0e-6) k++;
  }
  if (k < 2) VError(" Design matrix is singular");

  /* exception from Gaby */
  int j0 = 0;
  double xmin, xmax, tiny=10.0e-6;
  xmax = gsl_vector_get(w,0);
  xmin = tiny;
  for (j=n-1; j >= 0; j--) {
    u = gsl_vector_get(w,j);
    if (u > 0 && u/xmax > tiny) {
      j0 = j;
      goto skip;
    }
  }
 skip: ;
  if (j0 < n-1) {
    fprintf(stderr," Warning: Matrix almost singular\n");
    xmin = gsl_vector_get(w,j0) - tiny;
    if (xmin < 0) xmin = 0;
  }

  /* Fill the result matrix */
  gsl_matrix_float_set_zero(B);
  for (k=0; k<n; k++) {
    for (l=0; l<m; l++) {
      for (j=0; j<n; j++) {
	u = gsl_vector_get(w,j);
        if (fabs(u) > xmin) {
          x = gsl_matrix_float_get(B,k,l);
	  x += gsl_matrix_get(V,k,j)*gsl_matrix_get(U,l,j)/u;
          gsl_matrix_float_set(B,k,l,x);
        }
      }
    }
  }
  gsl_matrix_free(U);
  gsl_matrix_free(V);
  gsl_matrix_free(X);
  gsl_vector_free(w);
  gsl_vector_free(work);
  return B;
}

/*
** returns the trace of a matrix M
*/
float trace(gsl_matrix_float *M)
{
  gsl_vector_float_view diag = gsl_matrix_float_diagonal(M);
  float sum = 0;
  int i = 0;
  for(i = 0; i<diag.vector.size; i++) {
    sum += gsl_vector_float_get(&diag.vector,i);
  }
  return sum;
}

/*
**  rank(A)
*/
int rank(gsl_matrix_float* mat)
{
  int m = mat->size1;
  int n = mat->size2;
  int i;

  /* Create matrix buffer from source matrix. This prevents the overwriting of
   * matrix mat by gsl_linalg_SV_decomp().*/
  gsl_matrix* U = gsl_matrix_alloc(m,n);
  double* dpU =  U->data;
  float* fpMat = mat->data;
  for(i=0;i<mat->size1* mat->size2;i++) {
    *dpU++ = *fpMat++;
  }

  /* create additional matrix buffers */
  gsl_matrix *V;
  gsl_matrix *X;
  gsl_vector *S;
  gsl_vector *work;
  S = gsl_vector_alloc(n);
  V = gsl_matrix_alloc(n,n);
  X = gsl_matrix_alloc (n, n);
  work = gsl_vector_alloc (n);

  /* SVD */
  /* gsl_linalg_SV_decomp_jacobi(U,V,S); */
  gsl_linalg_SV_decomp_mod (U,X,V,S,work);

  int rank = 0;
  /* counting the nonzero singular values with a tolerance of EPSILON*/
  for(i=0;i<S->size;i++){
    if(S->data[i] > EPSILON) {
      rank++;
    }
  }

  gsl_matrix_free(U);
  gsl_matrix_free(V);
  gsl_matrix_free(X);
  gsl_vector_free(S);
  gsl_vector_free(work);

  return rank;
}

/*
**  fsum
*/
gsl_vector_float *fsum(gsl_matrix_float *matrix, int dim,gsl_vector_float *v)
{

  int row,col;
  float sum=0;

  /* build sum over all columns */
  if (dim == 1) {

    if (v == NULL) {
      v = gsl_vector_float_alloc(matrix->size2);
    }
    if(v->size != matrix->size2) {
      fprintf(stderr, "Warning in fsum: vector size doesn't match related matrix dimension. Resizing ..");
      gsl_vector_float_free(v);
      v = gsl_vector_float_alloc(matrix->size2);
    }

    for (col=0;col<matrix->size2;col++){
      sum=0;
      for(row=0;row<matrix->size1;row++) {
	sum+=matrix->data[col+row*matrix->size2];
      }
      v->data[col] = sum;
    }
  }

  /* build sum over all rows */
  else {

    if (v == NULL) {
      v = gsl_vector_float_alloc(matrix->size1);
    }

    if (v->size != matrix->size1) {
      fprintf(stderr, "Warning in fsum: vector size doesn't match related matrix dimension. Resizing ..");
      gsl_vector_float_free(v);
      v = gsl_vector_float_alloc(matrix->size1);
    }

    for (row = 0;row<matrix->size1;row++) {
      sum = 0;
      for(col=0;col<matrix->size2;col++) {
	sum += matrix->data[col+row*matrix->size2];
      }
      v->data[row] = sum;
    }
  }
  return v;
}

/*
**  funique
*/
gsl_vector_float *funique(gsl_vector_float* V)
{

  /* working copy */
  gsl_vector_float* v = gsl_vector_float_alloc(V->size);
  /* pointer to result */
  gsl_vector_float* res;
  /* holds value of last saved element */
  float max=0;
  /* counter */
  int i;

  gsl_vector_float_memcpy(v,V);

  /* sort elements from V. */
  gsl_sort_vector_float(v);

  /* count number of different elements */
  int nelements = 0;

  /* first run: count different values */
  float* p = v->data;
  for(i=0;i<v->size;i++) {
    /* the first element in the list is special */
    if(i==0) {
      max = *p;
      nelements++;
    }
    else {
      if(*p > max) {
	max = *p;
	nelements++;
      }
    }
    p++;
  }

  /* second run: allocate mem & copy every value exactly once. */
  res = gsl_vector_float_alloc(nelements);
  p = v->data;
  float *s = res->data;
  for(i=0;i<v->size;i++) {
    if(i==0) {
      max = *p;
      *s++ = *p;
    }
    else {
      if(*p > max) {
	max = *p;
	*s++ = *p;
      }
    }
    p++;
  }

  gsl_vector_float_free(v);

  return res;
}

/*
**  subcols
*/
gsl_matrix_float *fmat_subcols(gsl_matrix_float* mat,gsl_vector_float* cols)
{

  /* some tests */

  if((cols->size < 0) || (cols->size > mat->size2)) {
    fprintf(stderr,"column vector: invalid dimensions");
    exit(-1);
  }

  float min, max;
  gsl_vector_float_minmax(cols, &min, &max);

  if((min < 0) || (max > mat->size2)) {
    fprintf(stderr,"column vector values exceed matrix dimensions!");
    exit(-1);
  }

  /* TODO reduce column vector to unique values with funique. We will omit
   * this due to performance reasons. */

  /* the return value */
  gsl_matrix_float* ret = gsl_matrix_float_alloc(mat->size1, cols->size);
  /* the copy buffer */
  gsl_vector_float* buff = gsl_vector_float_alloc(mat->size1);
  /* counter */
  int i;

  /* copy columns from matrix to return buffer */
  for (i = 0; i < cols->size; ++i) {
    gsl_matrix_float_get_col(buff, mat, (int)cols->data[i]);
    gsl_matrix_float_set_col(ret, i, buff);
  }

  gsl_vector_float_free(buff);

  return ret;
}

/*
**  toeplitz
*/
gsl_matrix_float *fmat_toeplitz(gsl_vector_float* v, gsl_matrix_float* A)
{
  int i,j;

  if(A == NULL)
    A = gsl_matrix_float_alloc(v->size, v->size);
  else {
    if((A->size1 != v->size)  || (A->size2 != v->size)) {
      fprintf(stderr, "Warning fmat_toeplitz: incongruent matrix dimensions. Trying to\
                    correct it.");
      gsl_matrix_float_free(A);
      A = gsl_matrix_float_alloc(v->size, v->size);
    }
  }

  for(i=0;i<A->size1;i++){
    for(j=0;j<A->size2;j++) {
      gsl_matrix_float_set(A,i,j,gsl_vector_float_get(v, fabs(i-j)));
    }
  }

  return A;
}

/*
**  A^-1
*/
gsl_matrix_float *fInv(gsl_matrix_float *M, gsl_matrix_float *result)
{
  int s; /*permutation signum for LU-Decomposition*/
  int i;

  int m = M->size1;
  int n = M->size2;

  if(m != n) {
    fprintf(stderr, "dInv: not a square matrix\n");
    exit(0);
  }

  if(result == NULL) {
    result = gsl_matrix_float_alloc(m,m);
  }
  if((result->size1 != m) || (result->size2 != n)) {
    fprintf(stderr,"dInv: incongruent matrix dimensions.\n");
    exit(0);
  }

  gsl_matrix *ludecomp = gsl_matrix_alloc(m,m);
  gsl_permutation *perm = gsl_permutation_alloc(m);
  gsl_matrix *res = gsl_matrix_alloc(m,m);

  /* convert from float to double */
  double *d_ptr = ludecomp->data;
  float *f_ptr = M->data;
  for(i=0;i<M->size1*M->size2;i++) {
    *d_ptr = (double)*f_ptr;
    d_ptr++;f_ptr++;
  }

  /* LU decomposition */
  gsl_linalg_LU_decomp(ludecomp,perm,&s);

  /* inverting matrix */
  gsl_linalg_LU_invert(ludecomp,perm,res);

  /* convert from double to float */
  d_ptr = res->data;
  f_ptr = result->data;
  for(i=0;i<res->size1*res->size2;i++) {
    *f_ptr = (float)*d_ptr;
    d_ptr++;f_ptr++;
  }

  gsl_matrix_free(ludecomp);
  gsl_matrix_free(res);
  gsl_permutation_free(perm);
  return result;
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
//        Copyright The Authors 2018.
//    Distributed under the 3-Clause BSD License.
//    (See accompanying file LICENSE or copy at
//   https://opensource.org/licenses/BSD-3-Clause)

#pragma once

#include <cstdint>   // for uint8_t
#include <gsl/gsl>  // for gsl::span
#include <string>    // for std::string

namespace blocxxi {
namespace codec {
namespace hex {

std::string Encode(gsl::span<const uint8_t> src, bool reverse = false,
                   bool lower_case = false);

void Decode(gsl::span<const char> src, gsl::span<uint8_t> dest,
            bool reverse = false);

}  // namespace hex
}  // namespace codec
}  // namespace blocxxi

ext: h


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_core_errors_h_
#define SQ_INCLUDE_GUARD_core_errors_h_

#include "core/Primitive.fwd.h"
#include "core/Token.fwd.h"
#include "core/typeutil.h"

#include <cstddef>
#include <filesystem>
#include <fmt/format.h>
#include <fmt/ostream.h>
#include <gsl/gsl>
#include <stdexcept>
#include <string_view>
#include <system_error>

namespace sq {

/**
 * Base class for errors thrown by the SQ code.
 */
class Exception : public std::runtime_error {
public:
  using std::runtime_error::runtime_error;
};

/**
 * Indicates that a required parameter of a field is missing.
 */
class ArgumentMissingError : public Exception {
public:
  /**
   * @param arg_name the name of the missing argument.
   * @param arg_type the type of the missing argument.
   */
  ArgumentMissingError(std::string_view arg_name, std::string_view arg_type);
};

/**
 * Indicates that a given parameter is of incorrect type.
 */
class ArgumentTypeError : public Exception {
public:
  /**
   * @param received the value of the given parameter.
   * @param type_expected the name of the type expected for the parameter.
   */
  ArgumentTypeError(const Primitive &received, std::string_view type_expected);
};

/**
 * Indicates a programming error in SQ.
 *
 * InternalError should never actually be thrown - they are used in places that
 * the programmer believes are dead code, but where the C++ language still
 * requires e.g. a return statement.
 */
class InternalError : public Exception {
public:
  using Exception::Exception;
};

class InvalidConversionError : public Exception {
public:
  using Exception::Exception;
  InvalidConversionError(std::string_view from, std::string_view to);
};

/**
 * Indicates attempted access of a non-existent field.
 */
class InvalidFieldError : public Exception {
public:
  /**
   * @param sq_type the SQ type of the parent of the missing field.
   * @param field the name of the field that was requested.
   *
   * E.g. if the query is "a.b" and "b" is not a field of "a" then sq_type
   * should be the type of "a" and field should be "b".
   */
  InvalidFieldError(std::string_view sq_type, std::string_view field);
};

/**
 * Indicates incorrect grammar in a query.
 */
class ParseError : public Exception {
public:
  using Exception::Exception;

  /**
   * Create a ParseError for when an unexpected token is found.
   *
   * @param token the unexpected token.
   * @param expecting the set of tokens that would have been valid in place
   * of the unexpected token.
   */
  ParseError(const Token &token, const TokenKindSet &expecting);
};

/**
 * Indicates a failure to interpret part of the input query as a token.
 */
class LexError : public ParseError {
public:
  /**
   * @param pos position in the input query (in characters) at which the lex
   *      error occured.
   * @param query the full input query.
   */
  LexError(gsl::index pos, std::string_view query);
};

/**
 * Indicates that an array operation has been requested on a non-array type.
 */
class NotAScalarError : public Exception {
public:
  using Exception::Exception;
};

/**
 * Indicates that an array operation has been requested on a non-array type.
 */
class NotAnArrayError : public Exception {
public:
  using Exception::Exception;
};
/**
 * Indicates that a requested feature has not been implemented.
 */
class NotImplementedError : public Exception {
public:
  using Exception::Exception;
};

/**
 * Indicates that a request to access an element outside of an allowed range.
 */
class OutOfRangeError : public Exception {
public:
  using Exception::Exception;
  /**
   * @param token the token in the query where the access was requested.
   * @param message details about the requested access.
   */
  OutOfRangeError(const Token &token, std::string_view message);
};

/**
 * Indicates that a pullup type field access has been requested for a field
 * access with siblings.
 */
class PullupWithSiblingsError : public Exception {
  using Exception::Exception;
};

/**
 * Indicates an error was received from a system library.
 */
class SystemError : public Exception {
public:
  using Exception::Exception;

  /**
   * Create a SystemError object.
   *
   * @param operation the operation that failed.
   * @param code the system error code associated with the error.
   */
  SystemError(std::string_view operation, std::error_code code);

  /**
   * Get the system error code associated with this error.
   */
  SQ_ND std::error_code code() const;

private:
  std::error_code code_;
};

/**
 * Indicates that the udev library returned an error.
 */
class UdevError : public SystemError {
  using SystemError::SystemError;
};

/**
 * Indicates that a filesystem error occurred.
 */
class FilesystemError : public SystemError {
public:
  using SystemError::SystemError;

  /**
   * Create a FilesystemError object.
   *
   * @param operation the operation that failed.
   * @param path the path for which the operation failed.
   * @param code the system error code associated with the error.
   */
  FilesystemError(std::string_view operation, const std::filesystem::path &path,
                  std::error_code code);
};

class NarrowingError : public OutOfRangeError {
public:
  using OutOfRangeError::OutOfRangeError;

  NarrowingError(auto &&target, auto &&value, auto &&kind_of_value_format,
                 auto &&...format_args);

  NarrowingError(auto &&target, auto &&value);
};

} // namespace sq

#include "core/errors.inl.h"

#endif // SQ_INCLUDE_GUARD_core_errors_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
#ifndef UBLAS_LIB_PLASMA_H
#define UBLAS_LIB_PLASMA_H

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <pthread.h>
#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

#include <ublas_types.h>

int ubf_init(void **ctx);
int ubf_gemm(void *ctx, ublas_matrix *a, ublas_matrix *b, ublas_matrix *c, double alpha, double beta);
int ubf_free(void *ctx);

#endif

max_stars_repo_name: brgmnn/uob-ublas

max_issues_repo_name: brgmnn/uob-ublas

max_forks_repo_name: brgmnn/uob-ublas


text:
/* code by: iperetta@ieee.org */
#ifndef MYMATRIX_H_
#define MYMATRIX_H_

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <stdexcept>
#include <limits>
#include <lapacke.h>
#include "myfun.h"

using namespace std;

#ifndef MyEPS_
#define MyEPS_
static double EPS = numeric_limits<double>::epsilon();
#endif
// =====================================================================
template<typename T>
struct brace {
	T first, second;
};

template<typename T> brace<T> make_brace(const T &m1, const T &m2) {
	brace<T> ans;
	ans.first = m1;
	ans.second = m2;
	return ans;
}

template<typename T>
struct triplet {
	T first, second, third;
};

template<typename T> triplet<T> make_triplet(const T &m1, const T &m2,
		const T &m3) {
	triplet<T> ans;
	ans.first = m1;
	ans.second = m2;
	ans.third = m3;
	return ans;
}
// =====================================================================
class Matrix;
Matrix column_vector(const vector<double> &vec);
Matrix row_vector(const vector<double> &vec);
Matrix identityMatrix(const size_t &N);
Matrix diagonalMatrix(const vector<double> &vec);

class Matrix { // Matrix A_mXn
	friend Matrix operator+(const double &lhs, const Matrix &rhs);
	friend Matrix operator-(const double &lhs, const Matrix &rhs);
	friend Matrix operator*(const double &lhs, const Matrix &rhs);
	friend Matrix operator/(const double &lhs, const Matrix &rhs);
	friend ostream& operator<<(ostream &os, const Matrix& M);
  protected:
	vector<double> elements;
	void badConvergence(const string &str) const;
	void invalidArgument(const int &errid, const string &str) const;
	void checkDim(const size_t &i_, const string &str = "") const;
	void checkDim(const size_t &i_, const size_t &j_,
		const string &str = "") const;
	void checkSize(const size_t &i_, const size_t &j_,
		const string &str = "") const;
	void ensureSquare(const string &str = "") const;
	void ensureSameSize(const Matrix &B, const string &str = "") const;
	void ensureNonSingularity(const string &str = "") const;
	void ensureVector(const string &str) const;
  public:
	size_t m; // number of rows of the matrix
	size_t n; // number of columns of the matrix
	Matrix(); //constructor
    Matrix(const size_t &m_, const size_t &n_, const double &cns = 0.);
    Matrix(const size_t &m_, const size_t &n_, const double *elem,
		const size_t &sz_elem);
	Matrix(const size_t &m_, const size_t &n_,
		const vector<double> &elem);
    void swap(Matrix &);
    ~Matrix() {};//destructor
    void print() const;
    string str() const;
    double& at(const size_t &i_);
    const double& at(const size_t &i_) const;
    double& at(const size_t &i_, const size_t &j_);
    const double& at(const size_t &i_, const size_t &j_) const;
    vector<double> as_vector() const;
    vector<double> row(const size_t &r_) const;
    void assign2row(const size_t &r_, const double &value);
    void assign2row(const size_t &r_, const vector<double> &values);
    void assign2row(const size_t &r_, const double *values,
		const size_t &sz);
    vector<double> column(const size_t &c_) const;
    void assign2column(const size_t &c_, const double &value);
    void assign2column(const size_t &c_, const vector<double> &values);
    void assign2column(const size_t &c_, const double *values,
		const size_t &sz);
    vector<double> diagonal() const;
    void assign2diagonal(const double &value);
    void assign2diagonal(const vector<double> &values);
    void assign2diagonal(const double *values, const size_t &sz);
	void swap_rows(const size_t &r1, const size_t &r2);
	void swap_columns(const size_t &c1, const size_t &c2);

    void self_concat_row(const double &cte);
    void self_concat_row(const vector<double> &values);
	void self_concat_row(const Matrix &B);
	void self_concat_column(const double &cte);
	void self_concat_column(const vector<double> &values);
	void self_concat_column(const Matrix &B);
	void self_delete_row(const size_t &r_);
	void self_delete_column(const size_t &c_);

	Matrix concat_row(const double &cte) const;
    Matrix concat_row(const vector<double> &values) const;
	Matrix concat_row(const Matrix &B) const;
	Matrix concat_column(const double &cte) const;
	Matrix concat_column(const vector<double> &values) const;
	Matrix concat_column(const Matrix &B) const;
	Matrix delete_row(const size_t &r_) const;
	Matrix delete_column(const size_t &c_) const;

	size_t rank() const;
	double determinant() const;
	double trace() const;
	double norm() const;
	double infinity_norm() const;
	double frobenius_norm() const;
	double condition_number() const;
	double cofactor(const size_t &r_, const size_t &c_) const;
	double highest_eigenvalue() const;
    Matrix transpose() const;
    Matrix adjoint() const;
    Matrix inverse() const;
    brace<Matrix> eigen() const;
    triplet<Matrix> lu_p_factorization() const;
    brace<Matrix> qr_factorization() const;
    brace<Matrix> lq_factorization() const;
    triplet<Matrix> svd() const;
    Matrix pseudoinverse() const;

    Matrix operator+(const double &value) const;
    void operator+=(const double &value);
    Matrix operator-(const double &value) const;
    void operator-=(const double &value);
    Matrix operator*(const double &value) const;
    void operator*=(const double &value);
    Matrix operator/(const double &value) const;
    void operator/=(const double &value);
    Matrix operator+(const Matrix &value) const;
    void operator+=(const Matrix &value);
    Matrix operator-(const Matrix &value) const;
    void operator-=(const Matrix &value);
    Matrix operator*(const Matrix &value) const;
    void operator*=(const Matrix &value);
    bool operator==(const Matrix &B) const;
    bool operator!=(const Matrix &B) const;

    double min_all(const bool &absolute) const;
    double max_all(const bool &absolute) const;

    bool isFullRank() const;
    bool isOrthogonal() const;
    bool isDiagonal() const;
    bool isPositiveDefinite() const;
    bool isSingular() const;
    bool isSymmetric() const;
    bool isSkewSymmetric() const;
    bool isValid() const;

    void adjust_elements(const double &tol);

    void saveAs(const string &namefile) const;
};

Matrix::Matrix() : m(0), n(0) {
}

Matrix::Matrix(const size_t &m_, const size_t &n_, const double &cns) :
		m(m_), n(n_) {
	this->elements.assign(m*n,cns);
}

Matrix::Matrix(const size_t &m_, const size_t &n_, const double *elem,
		const size_t &sz_elem) : m(m_), n(n_) {
	this->checkSize(this->m*this->n, sz_elem,"(CONSTRUCTOR)");
	this->elements.assign(elem, elem+sz_elem);
}

Matrix::Matrix(const size_t &m_, const size_t &n_, const vector<double>
		&elem) : m(m_), n(n_) {
	this->checkSize(this->m*this->n, elem.size(),"(CONSTRUCTOR)");
	this->elements = elem;
}

void Matrix::swap(Matrix &other) {
    // need to swap all data members
    std::swap(this->elements, other.elements);
    std::swap(this->m, other.m);
    std::swap(this->n, other.n);
}
// =====================================================================
void Matrix::badConvergence(const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"the algorithm failed to converge.";
	throw std::runtime_error(MSG.c_str());
}

void Matrix::invalidArgument(const int &errid, const string &str)
		const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"invalid argument";
	cerr << "[ " << errid << " ] ";
	throw std::invalid_argument(MSG.c_str());
}

void Matrix::checkDim(const size_t &i_, const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"index exceeds matrix dimension";
	if(i_ >= this->elements.size())
		throw std::out_of_range(MSG.c_str());
}

void Matrix::checkDim(const size_t &i_, const size_t &j_,
		const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"index exceeds matrix dimension";
	if(i_ >= this->m || j_ >= this->n)
		throw std::out_of_range(MSG.c_str());
}

void Matrix::checkSize(const size_t &i_, const size_t &j_,
		const string &str) 	const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"sizes of matrices don't match";
	if(i_ != j_)
		throw std::invalid_argument(MSG.c_str());
}

void Matrix::ensureSquare(const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"exclusive to square matrices";
	if(this->m != this->n)
		throw std::invalid_argument(MSG.c_str());
}

void Matrix::ensureSameSize(const Matrix &B, const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"matrices must be of same size";
	if(this->m != B.m || this->n != B.n)
		throw std::invalid_argument(MSG.c_str());
}

void Matrix::ensureNonSingularity(const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"matrix is singular";
	if(this->isSingular())
		throw std::invalid_argument(MSG.c_str());
}

void Matrix::ensureVector(const string &str) const {
	string MSG = (str.empty() ? "" : "in ") + str +
		(str.empty() ? "" : ", ") +
		"matrices must be row or column vectors";
	if(this->m != 1 && this->n != 1)
		throw std::invalid_argument(MSG.c_str());
}
// =====================================================================
void Matrix::print() const {
	cout << endl << ">> Matrix_" << this->m << "X" << this->n
		 << " =" << endl;
	for(size_t i = 0; i < this->m; i++) {
		for(size_t j = 0; j < this->n; j++)
			cout << this->elements.at(i*this->n + j) << '\t';
		cout << endl;
	}
	cout << endl;
}

string Matrix::str() const {
	string ST = "";
	for(size_t i = 0; i < this->m; i++) {
		for(size_t j = 0; j < this->n; j++)
			ST += my::num2str<double>(this->elements.at(i*this->n + j))
				+ ((j == this->n-1) ? '\0' : '\t');
		ST += '\n';
	}
	return ST;
}

ostream& operator<<(ostream &os, const Matrix& M) {
	return os << M.str();
}

double& Matrix::at(const size_t &i_) {
	this->checkDim(i_,"AT");
	return this->elements.at(i_);
}

const double& Matrix::at(const size_t &i_) const {
	this->checkDim(i_,"AT");
	return this->elements.at(i_);
}

double& Matrix::at(const size_t &i_, const size_t &j_) {
	this->checkDim(i_,j_,"AT");
	return this->elements.at(i_*(this->n)+j_);
}

const double& Matrix::at(const size_t &i_, const size_t &j_) const {
	this->checkDim(i_,j_,"AT");
	return this->elements.at(i_*(this->n)+j_);
}

vector<double> Matrix::as_vector() const {
	return this->elements;
}

vector<double> Matrix::row(const size_t &r_) const {
	this->checkDim(r_,0,"ROW");
	vector<double> row_(this->elements.begin()+(r_*this->n),
		this->elements.begin()+(r_*this->n+this->n));
	return row_;
}

void Matrix::assign2row(const size_t &r_, const double &value) {
	this->checkDim(r_,0,"ASSIGN2ROW");
	for(size_t i = 0; i < this->n; i++)
		this->elements.at(r_*this->n+i) = value;
}

void Matrix::assign2row(const size_t &r_,
		const vector<double> &values) {
	this->checkDim(r_,0,"ASSIGN2ROW");
	this->checkSize(this->n,values.size(),"ASSIGN2ROW");
	for(size_t i = 0; i < this->n; i++)
		this->elements.at(r_*this->n+i) = values.at(i);
}

void Matrix::assign2row(const size_t &r_, const double *values,
		const size_t &sz) {
	vector<double> vec(values,values+sz);
	this->assign2row(r_,vec);
}

vector<double> Matrix::column(const size_t &c_) const {
	this->checkDim(0,c_,"COLUMN");
	vector<double> column_(this->m,0.);
	for(size_t i = 0; i < this->m; i++)
		column_.at(i) = this->elements.at(i*this->n+c_);
	return column_;
}

void Matrix::assign2column(const size_t &c_, const double &value) {
	this->checkDim(0,c_,"ASSIGN2COLUMN");
	for(size_t i = 0; i < this->m; i++)
		this->elements.at(i*this->n+c_) = value;
}

void Matrix::assign2column(const size_t &c_,
		const vector<double> &values) {
	this->checkDim(0,c_,"ASSIGN2COLUMN");
	this->checkSize(this->m,values.size(),"ASSIGN2COLUMN");
	for(size_t i = 0; i < this->m; i++)
		this->elements.at(i*this->n+c_) = values.at(i);
}

void Matrix::assign2column(const size_t &c_, const double *values,
		const size_t &sz) {
	vector<double> vec(values,values+sz);
	this->assign2column(c_,vec);
}

vector<double> Matrix::diagonal() const {
	size_t dim = min(this->m,this->n);
	vector<double> diag_(dim,0.);
	for(size_t i = 0; i < dim; i++)
		diag_.at(i) = this->elements.at(i*this->n+i);
	return diag_;
}

void Matrix::assign2diagonal(const double &value) {
	size_t dim = min(this->m,this->n);
	for(size_t i = 0; i < dim; i++)
		this->elements.at(i*this->n+i) = value;
}

void Matrix::assign2diagonal(const vector<double> &values) {
	size_t dim = min(this->m,this->n);
	this->checkSize(dim,values.size(),"ASSIGN2DIAGONAL");
	for(size_t i = 0; i < dim; i++)
		this->elements.at(i*this->n+i) = values.at(i);
}

void Matrix::assign2diagonal(const double *values, const size_t &sz) {
	vector<double> vec(values,values+sz);
	this->assign2diagonal(vec);
}

void Matrix::swap_rows(const size_t &r1, const size_t &r2) {
	this->checkDim(r1,0,"SWAP_ROW");
	this->checkDim(r2,0,"SWAP_ROW");
	double aux;
	for(size_t j = 0; j < this->n; j++) {
		aux = this->elements.at(r1*this->n+j);
		this->elements.at(r1*this->n+j) =
			this->elements.at(r2*this->n+j);
		this->elements.at(r2*this->n+j) = aux;
	}
}

void Matrix::swap_columns(const size_t &c1, const size_t &c2) {
	this->checkDim(0,c1,"SWAP_COLUMN");
	this->checkDim(0,c2,"SWAP_COLUMN");
	double aux;
	for(size_t i = 0; i < this->m; i++) {
		aux = this->elements.at(i*this->n+c1);
		this->elements.at(i*this->n+c1) =
			this->elements.at(i*this->n+c2);
		this->elements.at(i*this->n+c2) = aux;
	}
}

void Matrix::self_concat_row(const double &cte) {
	Matrix M(this->concat_row(cte));
	this->swap(M);
}

void Matrix::self_concat_row(const vector<double> &values) {
	Matrix M(this->concat_row(values));
	this->swap(M);
}

void Matrix::self_concat_row(const Matrix &B) {
	Matrix M(this->concat_row(B));
	this->swap(M);
}

void Matrix::self_concat_column(const double &cte) {
	Matrix M(this->concat_column(cte));
	this->swap(M);
}

void Matrix::self_concat_column(const vector<double> &values) {
	Matrix M(this->concat_column(values));
	this->swap(M);
}

void Matrix::self_concat_column(const Matrix &B) {
	Matrix M(this->concat_column(B));
	this->swap(M);
}

void Matrix::self_delete_row(const size_t &r_) {
	Matrix M(this->delete_row(r_));
	this->swap(M);
}

void Matrix::self_delete_column(const size_t &c_) {
	Matrix M(this->delete_column(c_));
	this->swap(M);
}

Matrix Matrix::concat_row(const double &cte) const {
	vector<double> values(this->n,cte);
	return concat_row(values);
}

Matrix Matrix::concat_row(const vector<double> &values) const {
	this->checkSize(this->n,values.size(),"CONCAT_ROW");
	Matrix M(this->m+1,this->n);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < M.n; j++)
			M.elements.at(i*M.n+j) = this->elements.at(i*this->n+j);
	for(size_t j = 0; j < M.n; j++)
		M.elements.at(this->m*M.n+j) = values.at(j);
	return M;
}

Matrix Matrix::concat_row(const Matrix &B) const {
	this->checkSize(this->n,B.n,"CONCAT_ROW");
	Matrix M(this->m+B.m,this->n);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < M.n; j++)
			M.elements.at(i*M.n+j) = this->elements.at(i*this->n+j);
	for(size_t i = 0; i < B.m; i++)
		for(size_t j = 0; j < M.n; j++)
			M.elements.at((this->m+i)*M.n+j) =
				B.elements.at(i*B.n+j);
	return M;
}

Matrix Matrix::concat_column(const double &cte) const {
	vector<double> values(this->m,cte);
	return concat_column(values);
}

Matrix Matrix::concat_column(const vector<double> &values) const {
	this->checkSize(this->m,values.size(),"CONCAT_COLUMN");
	Matrix M(this->m,this->n+1);
	for(size_t i = 0; i < M.m; i++)
		for(size_t j = 0; j < this->n; j++)
			M.elements.at(i*M.n+j) = this->elements.at(i*this->n+j);
	for(size_t i = 0; i < M.m; i++)
		M.elements.at(i*M.n+this->n) = values.at(i);
	return M;
}

Matrix Matrix::concat_column(const Matrix &B) const {
	this->checkSize(this->m,B.m,"CONCAT_COLUMN");
	Matrix M(this->m,this->n+B.n);
	for(size_t i = 0; i < M.m; i++)
		for(size_t j = 0; j < this->n; j++)
			M.elements.at(i*M.n+j) = this->elements.at(i*this->n+j);
	for(size_t i = 0; i < M.m; i++)
		for(size_t j = 0; j < B.n; j++)
			M.elements.at(i*M.n+(this->n+j)) = B.elements.at(i*B.n+j);
	return M;
}

Matrix Matrix::delete_row(const size_t &r_) const {
	this->checkDim(r_,0,"DELETE_ROW");
	Matrix M(this->m-1,this->n);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			if(i != r_)
				M.elements.at(((i<r_)?i:i-1)*M.n+j) =
					this->elements.at(i*this->n+j);
	return M;
}

Matrix Matrix::delete_column(const size_t &c_) const {
	this->checkDim(0,c_,"DELETE_COLUMN");
	Matrix M(this->m,this->n-1);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			if(j != c_)
				M.elements.at(i*M.n+((j<c_)?j:j-1)) =
					this->elements.at(i*this->n+j);
	return M;
}
// ===================================
size_t Matrix::rank() const {
	size_t rk = 0;
	this->checkSize((this->m)*(this->n),this->elements.size(),"RANK");
	lapack_int m = this->m, n = this->n;
	lapack_int lda = n, ldu = m, ldvt = n, info;
	double *superb = new double[min(m,n)-1];
	double *a = new double[m*n],
		   *s = new double[min(m,n)],
		   *u = new double[m*m],
		   *vt = new double[n*n];
	for(size_t i = 0; i < this->elements.size(); i ++)
		a[i] = this->elements.at(i);
	info = LAPACKE_dgesvd( LAPACK_ROW_MAJOR, 'N', 'N', m, n, a, lda,
		s, u, ldu, vt, ldvt, superb );
    if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"RANK");
	if(info > 0) 	// Check for convergence
		this->badConvergence("RANK");
	for(int i = 0; i < min(m,n); i++)
		rk += (my::isLikelyZero(s[i],max(this->m,this->n))) ? 0 : 1;
	delete [] superb;
	delete [] a;
	delete [] s;
	delete [] u;
	delete [] vt;
	return rk;
}

double Matrix::determinant() const {
	this->ensureSquare("DETERMINANT");
	double det;
	switch(this->m) {
		case 1:
			det = this->elements.at(0);
			break;
		case 2:
			det = this->elements.at(0)*this->elements.at(3) -
				this->elements.at(1)*this->elements.at(2);
			break;
		case 3:
			det = this->elements.at(0)*this->elements.at(4)*
					this->elements.at(8)
				+ this->elements.at(1)*this->elements.at(5)*
					this->elements.at(6)
				+ this->elements.at(2)*this->elements.at(3)*
					this->elements.at(7)
				- this->elements.at(2)*this->elements.at(4)*
					this->elements.at(6)
				- this->elements.at(0)*this->elements.at(5)*
					this->elements.at(7)
				- this->elements.at(1)*this->elements.at(3)*
					this->elements.at(8);
			break;
		default:
			det = 0;
			for(size_t j = 0; j < this->n; j++)
				det += this->cofactor(0,j)*this->elements.at(j); // i=0
	}
	return det;
}

double Matrix::trace() const {
	this->ensureSquare("TRACE");
	vector<double> diag(this->diagonal());
	double total(0.);
	for(vector<double>::iterator it = diag.begin();
		it < diag.end(); it++) total += *it;
	return total;
}

double Matrix::norm() const {
	// for square matrices, this is the same as the "spectral norm",
	// the "induced matrix 2-norm" or the "Euclidean norm"
	if(this->m != 1 && this->n != 1) {
		this->checkSize((this->m)*(this->n),this->elements.size(),
			"NORM");
		lapack_int m = this->m, n = this->n;
		lapack_int lda = n, ldu = m, ldvt = n, info;
		double *superb = new double[min(m,n)-1];
		double *a = new double[m*n],
			   *s = new double[min(m,n)],
			   *u = new double[m*m],
			   *vt = new double[n*n];
		for(size_t i = 0; i < this->elements.size(); i ++)
			a[i] = this->elements.at(i);
		info = LAPACKE_dgesvd( LAPACK_ROW_MAJOR, 'N', 'N', m, n, a, lda,
			s, u, ldu, vt, ldvt, superb );
		if(info < 0) 	// Check for invalid (NaN)
			this->invalidArgument((int)info,"NORM");
		if(info > 0) 	// Check for convergence
			this->badConvergence("NORM");
		double hSingVal = s[0];
		for(int i = 1; i < min(m,n); i++)
			if(s[i] > hSingVal)
				hSingVal = s[i];
		delete [] superb;
		delete [] a;
		delete [] s;
		delete [] u;
		delete [] vt;
		return hSingVal;
	} else
		return sqrt(((this->m == 1) ? (*this)*(this->transpose()) :
			(this->transpose())*(*this)).elements.at(0));
}

double Matrix::infinity_norm() const {
	// Maximum absolute row sum norm
	double infnorm = 0., sum_;
	for(size_t i = 0; i < this->m; i++) {
		sum_ = 0.;
		for(size_t j = 0; j < this->n; j++)
			sum_ += abs(this->elements.at(i*this->n+j));
		if(sum_ > infnorm)
			infnorm = sum_;
	}
	return infnorm;
}

double Matrix::frobenius_norm() const {
	this->checkSize((this->m)*(this->n),this->elements.size(),
		"FRBNORM");
	double sum = 0.;
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			sum += pow(this->elements.at(i*this->n+j),2.);
	return sqrt(sum);
}

double Matrix::condition_number() const
{
	triplet<Matrix> SVD = this->svd();
	vector<double> singular_values = SVD.second.diagonal();
	return my::max(singular_values)/my::min(singular_values);
}
double Matrix::cofactor(const size_t &r_, const size_t &c_) const {
	this->ensureSquare("COFACTOR");
	Matrix CF(*this);
	CF.self_delete_row(r_);
	CF.self_delete_column(c_);
	return (((r_+c_)%2)?-1:1)*CF.determinant();
}

Matrix Matrix::transpose() const {
	Matrix T(this->n,this->m);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			T.elements.at(j*this->m+i) = this->elements.at(i*this->n+j);
	return T;
}

Matrix Matrix::adjoint() const {
	this->ensureSquare("ADJOINT");
	Matrix A(this->m,this->n);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			A.elements.at(i*this->n+j) = cofactor(i,j);
	return A.transpose();
}

Matrix Matrix::inverse() const {
	this->ensureSquare("INVERSE");
	this->ensureNonSingularity("INVERSE");
	return this->adjoint()/this->determinant();
}

brace<Matrix> Matrix::eigen() const {
	this->checkSize((this->n)*(this->n),this->elements.size(),"EIGEN");
	this->ensureSquare("EIGEN");
	lapack_int n = this->n, info;
	double *a = new double[n*n],
		   *b = new double[n*n],
		   *alphar = new double[n],
		   *alphai = new double[n],
		   *beta = new double[n],
		   *vl = new double,
		   *vr = new double[n*n];
	for(size_t i = 0; i < this->elements.size(); i ++) {
		a[i] = this->elements.at(i);
		b[i] = (i%(n+1) == 0) ? 1. : 0.;
	}
	info = LAPACKE_dggev( LAPACK_ROW_MAJOR, 'N', 'V', n, a, n, b, n,
		alphar, alphai, beta, vl, 1, vr, n );
	if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"EIGEN");
	if(info > 0) 	// Check for convergence
		this->badConvergence("EIGEN");
	Matrix EVAL(n,n),
		   EVEC(n,n,vr,n*n);
	vector<double> lambda(n,0.);
	for(int i = 0; i < n; i++) {
		lambda.at(i) = (abs(alphai[i]) <= EPS) ? ((abs(beta[i]) <= EPS)
			? NAN : alphar[i]/beta[i] ) : -NAN;
	}
	EVAL.assign2diagonal(lambda);
	delete [] a;
	delete [] b;
	delete [] alphar;
	delete [] alphai;
	delete [] beta;
	delete vl;
	delete [] vr;
	return make_brace<Matrix>(EVAL,EVEC);
}

triplet<Matrix> Matrix::lu_p_factorization() const {
	lapack_int m = this->m, n = this->n, info;
	double *a = new double[m*n];
	lapack_int *ipiv = new lapack_int[max(1,min(m,n))];
	for(size_t i = 0; i < this->elements.size(); i ++) {
		a[i] = this->elements.at(i);
	}
	info = LAPACKE_dgetrf( LAPACK_ROW_MAJOR, m, n, a, n, ipiv );
	if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"LU_P_FACTORIZATION");
	if(info > 0) 	// Check for convergence
		cerr << ">> Warning! Matrix U is singular." << endl;
	Matrix L(m,min(m,n)),
		   U(min(m,n),n),
		   P(m,m);
	vector<double> ones(min(m,n),1.);
	L.assign2diagonal(ones);
	for(int i = 0; i < m; i++)
		for(int j = 0; j < n; j++)
			if(j >= i)
				U.at(i,j) = a[i*n+j];
			else
				L.at(i,j) = a[i*n+j];
	vector<size_t> seq;
	size_t aux;
	for(int i = 0; i < m; i++)
		seq.push_back(i);
	for(int i = 0; i < min(m,n); i++) {
		aux = seq.at(i);
		seq.at(i) = seq.at(ipiv[i]-1);
		seq.at(ipiv[i]-1) = aux;
	}
	for(int i = 0; i < m; i++)
		P.at(i,seq.at(i)) = 1;
	delete [] a;
	delete [] ipiv;
	return make_triplet<Matrix>(L,U,P);
}

brace<Matrix> Matrix::qr_factorization() const {
	lapack_int m = this->m, n = this->n, info;
	double *a = new double[m*n];
	double *tau = new double[max(1,min(m,n))];
	for(size_t i = 0; i < this->elements.size(); i ++) {
		a[i] = this->elements.at(i);
	}
    info = LAPACKE_dgeqrf( LAPACK_ROW_MAJOR, m, n, a, n, tau );
	if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"QR_FACTORIZATION");
	/* The matrix Q is represented as a product of elementary reflectors
	*
	*     Q = H(1) H(2) . . . H(k), where k = min(m,n).
	*
	*  Each H(i) has the form
	*
	*     H(i) = I - tau * v * v'
	*
	*  where tau is a real scalar, and v is a real vector with
	*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
	*  A(i+1:m,i), and tau in TAU(i).*/
	Matrix Q(m,min(m,n)),
		   R(min(m,n),n),
		   H(identityMatrix(m)),
		   v(m,1);
	vector<double> ones(min(m,n),1.);
	Q.assign2diagonal(ones);
	for(int i = 0; i < m; i++)
		for(int j = 0; j < n; j++)
			if(j >= i)
				R.at(i,j) = a[i*n+j];
			else
				Q.at(i,j) = a[i*n+j];
	for(int i = 0; i < min(m,n); i++) {
		v = column_vector(Q.column(i));
		H *= identityMatrix(m)-(v*v.transpose())*tau[i];
	}
	Q = H;
	if(m > n)
		for(int i = n; i < m; i++)
			Q.self_delete_column(Q.n-1);
	delete [] a;
	delete [] tau;
	return make_brace<Matrix>(Q,R);
}

brace<Matrix> Matrix::lq_factorization() const {
	lapack_int m = this->m, n = this->n, info;
	double *a = new double[m*n];
	double *tau = new double[max(1,min(m,n))];
	for(size_t i = 0; i < this->elements.size(); i ++) {
		a[i] = this->elements.at(i);
	}
    info = LAPACKE_dgelqf( LAPACK_ROW_MAJOR, m, n, a, n, tau );
	if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"LQ_FACTORIZATION");
	Matrix L(m,min(m,n)),
		   Q(min(m,n),n),
		   H(identityMatrix(n)),
		   v(n,1);
	vector<double> ones(min(m,n),1.);
	Q.assign2diagonal(ones);
	for(int i = 0; i < m; i++)
		for(int j = 0; j < n; j++)
			if(j <= i)
				L.at(i,j) = a[i*n+j];
			else
				Q.at(i,j) = a[i*n+j];
	for(int i = 0; i < min(m,n); i++) {
		v = row_vector(Q.row(i));
		H *= identityMatrix(n)-(v.transpose()*v)*tau[i];
	}
	Q = H.transpose();
	if(m < n)
		for(int i = m; i < n; i++)
			Q.self_delete_row(Q.m-1);
	delete [] a;
	delete [] tau;
	return make_brace<Matrix>(L,Q);
}

triplet<Matrix> Matrix::svd() const {
	this->checkSize((this->m)*(this->n),this->elements.size(),"SVD");
	lapack_int m = this->m, n = this->n;
	lapack_int lda = n, ldu = m, ldvt = n, info;
	double *superb = new double[min(m,n)-1];
	double *a = new double[m*n],
		   *s = new double[min(m,n)],
		   *u = new double[m*m],
		   *vt = new double[n*n];
	for(size_t i = 0; i < this->elements.size(); i ++)
		a[i] = this->elements.at(i);
	info = LAPACKE_dgesvd( LAPACK_ROW_MAJOR, 'A', 'A', m, n, a, lda,
		s, u, ldu, vt, ldvt, superb );
    if(info < 0) 	// Check for invalid (NaN)
		this->invalidArgument((int)info,"SVD");
	if(info > 0) 	// Check for convergence
		this->badConvergence("SVD");
	// Adjust zeros
	for(size_t i = 0; i < min(this->m,this->n); i++)
		if(my::isLikelyZero(s[i],max(this->m,this->n)))
			s[i] = 0.;
	for(int i = 0; i < m*m; i++)
		if(my::isLikelyZero(u[i],max(this->m,this->n)))
			u[i] = 0.;
	for(int i = 0; i < n*n; i++)
		if(my::isLikelyZero(vt[i],max(this->m,this->n)))
			vt[i] = 0.;
	// Create matrices
	Matrix S(m,n),
		   U(m,m,u,m*m),
		   VT(n,n,vt,n*n);
	S.assign2diagonal(s,min(m,n));
	delete [] superb;
	delete [] a;
	delete [] s;
	delete [] u;
	delete [] vt;
	// Return matrices U, Sigma and V
	// To retrieve original matrix: U*Sigma*V.transpose()
	return make_triplet<Matrix>(U,S,VT.transpose());
}

Matrix Matrix::pseudoinverse() const {
	triplet<Matrix> SVD = this->svd();
	vector<double> sv = SVD.second.diagonal();
	double maxsv = 0;
	for(size_t i = 0; i < sv.size(); i++)
		if(maxsv < sv.at(i))
			maxsv = sv.at(i);
	for(size_t i = 0; i < sv.size(); i++) {
		if(sv.at(i) <= this->m*this->n*EPS)
			sv.at(i) = 0.;
		else
			sv.at(i) = 1./sv.at(i);
	}
	SVD.second.assign2diagonal(sv);
	return (SVD.first*SVD.second*SVD.third.transpose()).transpose();
}
// =========
double Matrix::min_all(const bool &absolute) const
{
	double minelem = (absolute)? abs(this->elements.front()) :
		this->elements.front();
	for(size_t i = 1; i < this->elements.size(); i++)
		if((absolute)? abs(this->elements.at(i)) < minelem :
			this->elements.at(i) < minelem)
				minelem = (absolute)? abs(this->elements.at(i)) :
					this->elements.at(i);
	return minelem;
}
double Matrix::max_all(const bool &absolute) const
{
	double maxelem = (absolute)? abs(this->elements.front()) :
		this->elements.front();
	for(size_t i = 1; i < this->elements.size(); i++)
		if((absolute)? abs(this->elements.at(i)) > maxelem :
			this->elements.at(i) > maxelem)
				maxelem = (absolute)? abs(this->elements.at(i)) :
					this->elements.at(i);
	return maxelem;
}
// =========
Matrix Matrix::operator+(const double &value) const {
	Matrix A(*this);
	A += value;
	return A;
}

void Matrix::operator+=(const double &value) {
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) += value;
}

Matrix Matrix::operator-(const double &value) const {
	Matrix A(*this);
	A -= value;
	return A;
}

void Matrix::operator-=(const double &value) {
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) -= value;
}

Matrix Matrix::operator*(const double &value) const {
	Matrix A(*this);
	A *= value;
	return A;
}

void Matrix::operator*=(const double &value) {
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) *= value;
}

Matrix Matrix::operator/(const double &value) const {
	Matrix A(*this);
	A /= value;
	return A;
}

void Matrix::operator/=(const double &value) {
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) /= value;
}

Matrix Matrix::operator+(const Matrix &B) const {
	Matrix A(*this);
	A += B;
	return A;
}

void Matrix::operator+=(const Matrix &B) {
	this->ensureSameSize(B,"OPERATOR+");
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) += B.elements.
				at(i*this->n+j);
}

Matrix Matrix::operator-(const Matrix &B) const {
	Matrix A(*this);
	A -= B;
	return A;
}

void Matrix::operator-=(const Matrix &B) {
	this->ensureSameSize(B,"OPERATOR-");
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			this->elements.at(i*this->n+j) -= B.elements.
				at(i*this->n+j);
}

Matrix Matrix::operator*(const Matrix &B) const {
	this->checkSize(this->n,B.m,"OPERATOR*");
	Matrix R(this->m,B.n);
	for(size_t i = 0; i < this->m; i++)
		for(size_t j = 0; j < this->n; j++)
			for(size_t k = 0; k < B.n; k++) {
				R.elements.at(i*B.n+k) += this->elements.
					at(i*this->n+j)*B.elements.at(j*B.n+k);
				if(abs(R.elements.at(i*B.n+k)) <= this->m*EPS)
					R.elements.at(i*B.n+k) = 0.;
			}
	return R;
}

void Matrix::operator*=(const Matrix &B) {
	Matrix A = (*this)*B;
	this->swap(A);
}

bool Matrix::operator==(const Matrix &B) const {
	if(this->m != B.m || this->n != B.n)
		return false;
	else {
		bool chk = true;
		for(size_t i = 0; i < this->m; i++) {
			for(size_t j = 0; j < this->n; j++) {
				if(abs(this->elements.at(i*this->n+j)-
						B.elements.at(i*B.n+j)) > 0.5*EPS) {
					chk = false;
					break;
				}
			}
			if(!chk) break;
		}
		return chk;
	}
}

bool Matrix::operator!=(const Matrix &B) const {
	return !((*this) == B);
}
// ======
Matrix operator+(const double &lhs, const Matrix &rhs) {
	return rhs+lhs;
}
Matrix operator-(const double &lhs, const Matrix &rhs) {
	return rhs*(-1.)+lhs;
}
Matrix operator*(const double &lhs, const Matrix &rhs) {
	return rhs*lhs;
}
Matrix operator/(const double &lhs, const Matrix &rhs) {
	// This is not an inverse! It's an element-wise division
	Matrix result(rhs);
	for(size_t i = 0; i < result.m; i++)
		for(size_t j = 0; j < result.n; j++)
			result.at(i,j) = lhs/result.at(i,j);
	return result;
}
// ======
bool Matrix::isFullRank() const {
	return this->rank() == min(this->m,this->n);
}

bool Matrix::isOrthogonal() const {
	Matrix O = (*this)*(this->transpose());
	return (O == identityMatrix(O.m));
}

bool Matrix::isDiagonal() const {
	Matrix D = diagonalMatrix(this->diagonal());
	if(this->m != this->n) {
		if(this->m > this->n) {
			vector<double> zeros(this->n,0.);
			for(size_t i = this->n; i < this->m; i++)
				D.self_concat_row(zeros);
		} else {
			vector<double> zeros(this->m,0.);
			for(size_t i = this->m; i < this->n; i++)
				D.self_concat_column(zeros);
		}
	}
	return (*this == D);
}

bool Matrix::isPositiveDefinite() const {
	//real symmetric matrix
	if(!(this->isSymmetric()))
		return false;
	else {
		// if all its eigenvalues are positive.
		cout << "not implemented yet!" << endl;
		return false;
	}
}

bool Matrix::isSingular() const {
	return (my::isLikelyZero(this->determinant()));
}

bool Matrix::isSymmetric() const {
	return (*this == this->transpose());
}

bool Matrix::isSkewSymmetric() const {
	return (*this == this->transpose()*(-1.));
}

bool Matrix::isValid() const {
	bool chk = true;
	for(vector<double>::const_iterator it = this->elements.begin();
		it != this->elements.end(); it++)
	{
		if(std::isnan(*it) || my::isinf(*it))
		{
			chk = false;
			break;
		}
	}
	return chk;
}

void Matrix::adjust_elements(const double &tol)
{
	for(size_t i = 0; i < this->elements.size(); i++)
		if(my::isLikelyZero(this->elements.at(i),this->m*this->n) ||
				abs(this->elements.at(i)) < tol)
			this->elements.at(i) = 0.;
}

void Matrix::saveAs(const string &namefile) const
{
	ofstream matrixfile;
	matrixfile.open(namefile);
	if (matrixfile.is_open())
	{
		for(size_t i = 0; i < this->m; i++)
		{
			for(size_t j = 0; j < this->n; j++)
				matrixfile << this->elements.at(i*(this->n)+j)
						   << ((j+1 == this->n) ? "" : "\t");
			matrixfile << endl;
		}
	}
	else
	{
		cerr << "Unable to open " << namefile << " file!" << endl;
	}
	matrixfile.close();
}
// =====================================================================
Matrix column_vector(const vector<double> &vec) {
	Matrix A(vec.size(),1,vec);
	return A;
}

Matrix column_vector(const size_t &qty, const double &val) {
	vector<double> vec(qty,val);
	Matrix A(qty,1,vec);
	return A;
}

Matrix row_vector(const vector<double> &vec) {
	Matrix A(1,vec.size(),vec);
	return A;
}

Matrix row_vector(const size_t &qty, const double &val) {
	vector<double> vec(qty,val);
	Matrix A(1,qty,vec);
	return A;
}

Matrix subMatrix(const Matrix &M, const size_t &fstr,
		const size_t &fstc, const size_t &lstr, const size_t &lstc) {
	cout << "row: " << fstr << " " << lstr << endl;
	cout << "col: " << fstc << " " << lstc << endl;
	cout << "M: " << M.m << " " << M.n << endl;
	if(fstr > lstr || lstr >= M.m || fstc > lstc || lstc >= M.n) {
		string MSG = "in SUBMATRIX, index exceeds matrix dimension";
		throw std::out_of_range(MSG.c_str());
	}
	Matrix A(lstr-fstr+1,lstc-fstc+1);
	cout << A.m << " " << A.n << endl;
	for(size_t i = 0; i < A.m; i++)
		for(size_t j = 0; j < A.n; j++)
			A.at(i,j) = M.at(i+fstr,j+fstc);
	return A;
}

Matrix identityMatrix(const size_t &N) {
	Matrix A(N,N);
	A.assign2diagonal(1.);
	return A;
}

Matrix diagonalMatrix(const vector<double> &vec) {
	Matrix A(vec.size(),vec.size());
	A.assign2diagonal(vec);
	return A;
}

Matrix linear_AASTtoRowVector_VAR(const AAST &tree, const size_t &n_var)
{
	Matrix A, b(n_var, 1);
	vector<double> aux(n_var, 0.);
	double cnst = tree.eval_with_variables(aux);
	aux.at(0) = 1.;
	A = row_vector(aux);
	b.at(0) = tree.eval_with_variables(aux) - cnst;
	for(size_t i = 1; i < n_var; i++)
	{
		aux.assign(n_var, 0.);
		aux.at(i) = 1.;
		A.self_concat_row(row_vector(aux));
		b.at(i) = tree.eval_with_variables(aux) - cnst;
	}
	return ((A.pseudoinverse()*b).transpose()).concat_column(cnst);
}

Matrix linear_AASTtoRowVector_UNK(const AAST &tree, const size_t &n_unk)
{
	Matrix A, b(n_unk, 1);
	vector<double> aux(n_unk, 0.);
	double cnst = tree.eval_with_unknowns(aux);
	aux.at(0) = 1.;
	A = row_vector(aux);
	b.at(0) = tree.eval_with_unknowns(aux) - cnst;
	for(size_t i = 1; i < n_unk; i++)
	{
		aux.assign(n_unk, 0.);
		aux.at(i) = 1.;
		A.self_concat_row(row_vector(aux));
		b.at(i) = tree.eval_with_unknowns(aux) - cnst;
	}
	return ((A.pseudoinverse()*b).transpose()).concat_column(cnst);
}
// =====================================================================
#endif

max_stars_repo_name: iperetta/system-modelling

max_issues_repo_name: iperetta/system-modelling

max_forks_repo_name: iperetta/system-modelling


text:
//
//  mclib_pluto.c
//  to read in pluto AMR chombo files
//
//  Created by Tyler Parsotan on 11/26/19.
//

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdbool.h>
#include <unistd.h>
#include <dirent.h>
#include <limits.h>
#include "hdf5.h"
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_bessel.h>
#include "mclib.h"

#include <omp.h>

void readPlutoChombo( char pluto_file[200], double r_inj, double fps, double **x, double **y, double **szx, double **szy, double **r,\
double **theta, double **velx, double **vely, double **dens, double **pres, double **gamma, double **dens_lab, double **temp, int *number, int ph_inj_switch, double min_r, double max_r, double min_theta, double max_theta, FILE *fPtr)
{
    hid_t  file, dset, space, group, attr;
    herr_t status;
    hsize_t dims[1]={0}; //hold number of processes in each level
    int i=0, j=0, k=0, l=0, m=0, r_count=0, num_dims=0, num_levels=0, num_vars=4, logr=0;
    int nbx=0, nby=0, nbz=0, total_size=0, total_box_size=0, offset=0, elem_factor=0;
    box2d prob_domain[1]={0,0,0,0};
    char level[200]="";
    double ph_rmin=0, ph_rmax=0, ph_thetamin=0, ph_thetamax=0, r_grid_innercorner=0, r_grid_outercorner=0, theta_grid_innercorner=0, theta_grid_outercorner=0;
    int *level_dims=NULL, *box_offset=NULL;
    double *radii=NULL, *angles=NULL, *dradii=NULL, *dangles=NULL;
    double *dombeg1=NULL, *dombeg2=NULL, *dombeg3=NULL, *dx=NULL, *g_x2stretch=NULL, *g_x3stretch=NULL;
    double *all_data=NULL, *x1_buffer=NULL, *x2_buffer=NULL, *dx1_buffer=NULL, *dx2_buffer=NULL;
    double *dens_buffer=NULL, *pres_buffer=NULL, *vel_x2_buffer=NULL, *vel_x1_buffer=NULL;

    if (ph_inj_switch==0)
    {
        ph_rmin=min_r;
        ph_rmax=max_r;
        ph_thetamin=min_theta-2*0.017453292519943295; //min_theta - 2*Pi/180 (2 degrees)
        ph_thetamax=max_theta+2*0.017453292519943295; //max_theta + 2*Pi/180 (2 degrees)
    }

    //define dataset for boxes of each level
    hid_t box_dtype = H5Tcreate (H5T_COMPOUND, sizeof(box2d));
    H5Tinsert(box_dtype, "lo_i", HOFFSET(box2d, lo_i), H5T_NATIVE_INT);
    H5Tinsert(box_dtype, "lo_j", HOFFSET(box2d, lo_j), H5T_NATIVE_INT);
    H5Tinsert(box_dtype, "hi_i", HOFFSET(box2d, hi_i), H5T_NATIVE_INT);
    H5Tinsert(box_dtype, "hi_j", HOFFSET(box2d, hi_j), H5T_NATIVE_INT);

    //open the pluto file
    file = H5Fopen (pluto_file, H5F_ACC_RDONLY, H5P_DEFAULT);

    fprintf(fPtr, ">> MCRaT: Reading positional, density, pressure, and velocity information...\n");
    fflush(fPtr);

    //1. read in the number of dims
    group = H5Gopen (file, "/Chombo_global", H5P_DEFAULT);

    attr = H5Aopen (group, "SpaceDim", H5P_DEFAULT);
    status = H5Aread (attr, H5T_NATIVE_INT, &num_dims);

    status = H5Aclose (attr);
    status = H5Gclose (group);

    //printf("readPlutoChombo spacedim: %d\n", num_dims);

    //2. get the number of levels
    attr = H5Aopen (file, "num_levels", H5P_DEFAULT);
    status = H5Aread (attr, H5T_NATIVE_INT, &num_levels);

    status = H5Aclose (attr);
    //printf("readPlutoChombo num_levels: %d\n", num_levels);

    dombeg1=malloc(num_levels*sizeof(double));
    dombeg2=malloc(num_levels*sizeof(double));
    dombeg3=malloc(num_levels*sizeof(double));
    dx=malloc(num_levels*sizeof(double));
    g_x2stretch=malloc(num_levels*sizeof(double));
    g_x3stretch=malloc(num_levels*sizeof(double));
    level_dims=malloc(num_levels*sizeof(int));

    //3. get number of variables to read in (should be 4 in non-MHD case)
    attr = H5Aopen (file, "num_components", H5P_DEFAULT);
    status = H5Aread (attr, H5T_NATIVE_INT, &num_vars);

    status = H5Aclose (attr);
    //printf("readPlutoChombo num_vars: %d\n", num_vars);

    //get the total number of values that I need to allocate memory for
    for (i=0;i<num_levels;i++)
    {
        snprintf(level, sizeof(level), "level_%d", i);
        //printf("Opening level %d Boxes\n", i);

        group = H5Gopen(file, level, H5P_DEFAULT);

        dset= H5Dopen(group, "data:datatype=0", H5P_DEFAULT);

        //get dimensions of array and save it
        space = H5Dget_space (dset);
        H5Sget_simple_extent_dims(space, dims, NULL); //save dimesnions in dims

        total_size+=dims[0];
        *(level_dims+i)=dims[0];

        status = H5Sclose (space);
        H5Dclose(dset);
        H5Gclose(group);
    }
    //printf("The total number of elements is %d\n", total_size);

    //now allocate space to save data
    all_data=malloc(total_size*sizeof (double));

    //and allocate arrays to hold r and theta values to calculate them on the fly
    x1_buffer=malloc((total_size/num_vars)*sizeof (double));
    x2_buffer=malloc((total_size/num_vars)*sizeof (double));
    dx1_buffer=malloc((total_size/num_vars)*sizeof (double));
    dx2_buffer=malloc((total_size/num_vars)*sizeof (double));

    vel_x1_buffer= malloc ((total_size/num_vars) * sizeof (double));
    vel_x2_buffer=malloc ((total_size/num_vars) * sizeof (double));
    dens_buffer= malloc ((total_size/num_vars) * sizeof (double));
    pres_buffer=malloc ((total_size/num_vars) * sizeof (double));

    offset=0;
    //read in the data
    for (i=0;i<num_levels;i++)
    {
        snprintf(level, sizeof(level), "level_%d", i);
        //printf("Opening level %d Boxes\n", i);

        group = H5Gopen(file, level, H5P_DEFAULT);

        //read in the data
        dset= H5Dopen(group, "data:datatype=0", H5P_DEFAULT);
        status = H5Dread (dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,(all_data+offset));
        H5Dclose(dset);
        //printf("First few data %e %e %e Last few data %e %e\n", *(all_data+offset), *(all_data+offset+1), *(all_data+offset+2), *(all_data+offset+(*(level_dims+i))-2), *(all_data+offset+(*(level_dims+i))-1));

        //read in the box offsets in all_data
        dset= H5Dopen(group, "data:offsets=0", H5P_DEFAULT);

        space = H5Dget_space (dset);
        H5Sget_simple_extent_dims(space, dims, NULL); //save dimesnions in dims
        box_offset=malloc(dims[0]*sizeof(int));
        status = H5Sclose (space);

        status = H5Dread (dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, box_offset);
        H5Dclose(dset);

        //open various properties about that refinement level
        attr = H5Aopen (group, "prob_domain", H5P_DEFAULT);
        status = H5Aread (attr, box_dtype, &prob_domain);
        status = H5Aclose (attr);
        //printf("Prob_domain %d %d %d %d\n", prob_domain->lo_i, prob_domain->lo_j, prob_domain->hi_i, prob_domain->hi_j);

        attr = H5Aopen (group, "dx", H5P_DEFAULT);
        status = H5Aread (attr, H5T_NATIVE_DOUBLE, (dx+i));
        status = H5Aclose (attr);
        //printf("dx %e\n", *(dx+i));

        attr = H5Aopen (group, "logr", H5P_DEFAULT);
        status = H5Aread (attr, H5T_NATIVE_INT, &logr);
        status = H5Aclose (attr);
        //printf("logr %d\n", logr);

        attr = H5Aopen (group, "domBeg1", H5P_DEFAULT);
        status = H5Aread (attr, H5T_NATIVE_DOUBLE, (dombeg1+i));
        status = H5Aclose (attr);
        //printf("dombeg1 %e\n", *(dombeg1+i));

        //set default just in case
        *(dombeg2+i)=0;
        *(dombeg3+i)=0;

        if (num_dims==2)
        {
            attr = H5Aopen (group, "g_x2stretch", H5P_DEFAULT);
            status = H5Aread (attr, H5T_NATIVE_DOUBLE, (g_x2stretch+i));
            status = H5Aclose (attr);
            //printf("g_x2stretch %e\n", *(g_x2stretch+i));

            attr = H5Aopen (group, "domBeg2", H5P_DEFAULT);
            status = H5Aread (attr, H5T_NATIVE_DOUBLE, (dombeg2+i));
            status = H5Aclose (attr);
            //printf("dombeg2 %e\n", *(dombeg2+i));

        }
        else if (num_dims==3)
        {
            attr = H5Aopen (group, "g_x3stretch", H5P_DEFAULT);
            status = H5Aread (attr, H5T_NATIVE_DOUBLE, (g_x3stretch+i));
            status = H5Aclose (attr);

            attr = H5Aopen (group, "domBeg3", H5P_DEFAULT);
            status = H5Aread (attr, H5T_NATIVE_DOUBLE, (dombeg3+i));
            status = H5Aclose (attr);
        }

        //read in the boxes
        dset= H5Dopen(group, "boxes", H5P_DEFAULT);

        //get dimensions of array and save it
        space = H5Dget_space (dset);
        H5Sget_simple_extent_dims(space, dims, NULL); //save dimesnions in dims

        box2d box_data[dims[0]];

        status = H5Dread (dset, box_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT,box_data);
        status = H5Sclose (space);
        H5Dclose(dset);

        radii=malloc( (prob_domain->hi_i - prob_domain->lo_i +1) * sizeof (double ));
        angles=malloc( (prob_domain->hi_j - prob_domain->lo_j +1) * sizeof (double ));
        dradii=malloc( (prob_domain->hi_i - prob_domain->lo_i +1) * sizeof (double ));
        dangles=malloc( (prob_domain->hi_j - prob_domain->lo_j +1) * sizeof (double ));

        //create the arrays that hold the refinement level radii and angles
        for (j=0;j<(prob_domain->hi_i - prob_domain->lo_i +1);j++)
        {
            if (logr==0)
            {
                *(radii+j)=(*(dombeg1+i)) + (*(dx+i)) * (prob_domain->lo_i + j + 0.5);
                *(dradii+j)=(*(dx+i));
            }
            else
            {
                *(radii+j)=(*(dombeg1+i)) * 0.5 * (exp((*(dx+i)) * (prob_domain->lo_i + j + 1)) + exp((*(dx+i)) * (prob_domain->lo_i + j))   );
                *(dradii+j)=(*(dombeg1+i)) * (exp((*(dx+i)) * (prob_domain->lo_i + j + 1)) - exp((*(dx+i)) * (prob_domain->lo_i + j))   );
            }
            //if (i==2) printf("radii: %0.8e dr: %0.8e\n", *(radii+j), *(dradii+j));
        }

        for (j=0;j<(prob_domain->hi_j - prob_domain->lo_j +1);j++)
        {
            *(angles+j)=(*(dombeg2+i)) + (*(dx+i)) * (*(g_x2stretch+i)) * (prob_domain->lo_j + j + 0.5);
            *(dangles+j)=(*(dx+i))*(*(g_x2stretch+i));
        }

        //go through the boxes to create the r and theta arrays
        total_box_size=0;
        for (j=0; j<dims[0]; j++)
        {
            //printf("i %d %d %d %d %d \n", j, box_data[j].lo_i, box_data[j].lo_j, box_data[j].hi_i, box_data[j].hi_j);
            nbx=box_data[j].hi_i-box_data[j].lo_i+1;
            nby=1;
            nbz=1;

            if (num_dims >1)
            {
                nby=box_data[j].hi_j-box_data[j].lo_j+1;
            }
            /* for future 3D support
            else if (dim>2)
            {
                nbz=box_data[j].hi_k-box_data[j].lo_k+1;
            }
             */
            //loop over each variable values of box
            for (k=0;k<num_vars;k++)
            {
                //loop over the radii
                for (l=0; l<nby; l++)
                {
                    //loop over the angles
                    for (m=0 ;m<nbx ;m++)
                    {
                        switch (k)
                        {
                            case 0:
                                //when k is 0 save the radii, caus eonly need to do once and also save the densities
                                #if GEOMETRY == SPHERICAL
                                    *(x1_buffer+offset/num_vars+(*(box_offset+j))/num_vars+ l*nbx +m  )= (*(radii+box_data[j].lo_i+m))*HYDRO_L_SCALE;
                                    *(x2_buffer+(offset+(*(box_offset+j)))/num_vars + l*nbx +m )=(*(angles+box_data[j].lo_j+l));
                                    *(dx1_buffer+(offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(dradii+box_data[j].lo_i+m))*HYDRO_L_SCALE;
                                    *(dx2_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m )=(*(dangles+box_data[j].lo_j+l));
                                #elif GEOMETRY == CARTESIAN
                                    *(x1_buffer+offset/num_vars+(*(box_offset+j))/num_vars+ l*nbx +m  )= (*(radii+box_data[j].lo_i+m))*HYDRO_L_SCALE;
                                    *(x2_buffer+(offset+(*(box_offset+j)))/num_vars + l*nbx +m )=(*(angles+box_data[j].lo_j+l))*HYDRO_L_SCALE;
                                    *(dx1_buffer+(offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(dradii+box_data[j].lo_i+m))*HYDRO_L_SCALE;
                                    *(dx2_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m )=(*(dangles+box_data[j].lo_j+l))*HYDRO_L_SCALE;
                                #endif

                                *(dens_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(all_data+offset+(*(box_offset+j))+ k*nbx*nby + l*nbx +m  ));
                                break;
                            case 1:
                                *(vel_x1_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(all_data+offset+(*(box_offset+j))+ k*nbx*nby + l*nbx +m  ));
                                break;
                            case 2:
                                *(vel_x2_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(all_data+offset+(*(box_offset+j))+ k*nbx*nby + l*nbx +m  ));
                                break;
                            case 3:
                                *(pres_buffer+ (offset+(*(box_offset+j)))/num_vars + l*nbx +m)=(*(all_data+offset+(*(box_offset+j))+ k*nbx*nby + l*nbx +m  ));
                                break;
                            default:
                                break;
                        }

                    }
                }
            }
            /*
              // this was for testing, the actual nested loop is below
              if (i==2 && j==dims[0]-1)
              {

                  r_count=0;
                  for (k=0;k<1;k++)
                  {
                      //loop over the radii
                      for (l=0; l<nby; l++)
                      {
                          //loop over the angles
                          for (m=0 ;m<nbx ;m++)
                          {
                              printf("count: %d idx: %d all_data val: %0.8e r: %e theta %e\n", (*(box_offset+j))+r_count, (*(box_offset+j))+k*nbx*nby + l*nby +m, *(all_data+offset+(*(box_offset+j))+ k*nbx*nby + l*nbx +m  ), *(radii+box_data[j].lo_i+m), *(angles+box_data[j].lo_j+l));
                              r_count++;
                          }
                          printf("\n");
                      }
                      printf("\n");
                  }

                  for (k=0;k<64;k++)
                  {
                      printf("r: %e, theta %e dens data: %e\n",*(x1_buffer+ (offset+(*(box_offset+j)))/num_vars+k), *(x2_buffer+ (offset+(*(box_offset+j)))/num_vars +k), *(dens_buffer+ (offset+(*(box_offset+j)))/num_vars +k) );
                  }
                  printf("\n");

              }
        */

        }

        offset+=(*(level_dims+i));

        H5Gclose(group);
        free(radii); free(angles); free(dradii); free(dangles); free(box_offset);

    }
    status = H5Fclose (file);
    free(dombeg1); free(dombeg2); free(dombeg3); free(dx); free(g_x2stretch); free(g_x3stretch); free(all_data);

    //have all the data so need to go through and decide which values we will be keeping based on phtoon ranges we want to keep
    //fill in radius array and find in how many places r > injection radius
    #if SYNCHROTRON_SWITCH == ON
        elem_factor=2;
    #else
        elem_factor=0;
    #endif
    r_count=0;
    while (r_count==0)
    {
        r_count=0;
        elem_factor++;
        for (i=0;i<(total_size/num_vars);i++)
        {
            if (ph_inj_switch==0)
            {
                #if GEOMETRY == SPHERICAL
                    r_grid_innercorner = (*(x1_buffer+i)) - 0.5 * (*(dx1_buffer+i));
                    r_grid_outercorner = (*(x1_buffer+i)) + 0.5 * (*(dx1_buffer+i));

                    theta_grid_innercorner = (*(x2_buffer+i)) - 0.5 * (*(dx2_buffer+i));
                    theta_grid_outercorner = (*(x2_buffer+i)) + 0.5 * (*(dx2_buffer+i));
                #elif GEOMETRY == CARTESIAN
                    r_grid_innercorner = pow((*(x1_buffer+i) - *(dx1_buffer+i)/2.0) * ((*(x1_buffer+i) - *(dx1_buffer+i)/2.0))+(*(x2_buffer+i) - *(dx2_buffer+i)/2.0) * (*(x2_buffer+i) - *(dx2_buffer+i)/2.0),0.5);
                    r_grid_outercorner = pow((*(x1_buffer+i) + *(dx1_buffer+i)/2.0) * ((*(x1_buffer+i) + *(dx1_buffer+i)/2.0))+(*(x2_buffer+i) + *(dx2_buffer+i)/2.0) * (*(x2_buffer+i) + *(dx2_buffer+i)/2.0),0.5);

                    theta_grid_innercorner = acos( (*(x2_buffer+i) - *(dx2_buffer+i)/2.0) /r_grid_innercorner); //arccos of y/r for the bottom left corner
                    theta_grid_outercorner = acos( (*(x2_buffer+i) + *(dx2_buffer+i)/2.0) /r_grid_outercorner);
                #endif

                if (((ph_rmin - elem_factor*C_LIGHT/fps) <= r_grid_outercorner) && (r_grid_innercorner  <= (ph_rmax + elem_factor*C_LIGHT/fps) ) && (theta_grid_outercorner >= ph_thetamin) && (theta_grid_innercorner <= ph_thetamax) )
                {
                    r_count++;
                }

            }
            else
            {
                if ((*(x1_buffer+i)) > (0.95*r_inj) )
                {
                    r_count++;
                }
            }
        }
        //printf( "r_count: %d\n", r_count);
    }
    fprintf(fPtr, "Elem factor: %d Ph_rmin: %e rmax: %e Chosen FLASH min_r: %e max_r: %e min_theta: %e degrees max_theta: %e degrees\n", elem_factor, ph_rmin, ph_rmax, ph_rmin - (elem_factor*C_LIGHT/fps), ph_rmax + (elem_factor*C_LIGHT/fps), ph_thetamin*180/M_PI, ph_thetamax*180/M_PI);
    fflush(fPtr);

    //allocate memory to hold processed data
    (*pres)=malloc (r_count * sizeof (double ));
    (*velx)=malloc (r_count * sizeof (double ));
    (*vely)=malloc (r_count * sizeof (double ));
    (*dens)=malloc (r_count * sizeof (double ));
    (*x)=malloc (r_count * sizeof (double ));
    (*y)=malloc (r_count * sizeof (double ));
    (*r)=malloc (r_count * sizeof (double ));
    (*theta)=malloc (r_count * sizeof (double ));
    (*gamma)=malloc (r_count * sizeof (double ));
    (*dens_lab)=malloc (r_count * sizeof (double ));
    (*szx)=malloc (r_count * sizeof (double ));
    (*szy)=malloc (r_count * sizeof (double ));
    (*temp)=malloc (r_count * sizeof (double ));

    //assign values based on r> 0.95*r_inj
    j=0;
    for (i=0;i<(total_size/num_vars);i++)
    {
        if (ph_inj_switch==0)
        {
           #if GEOMETRY == SPHERICAL
               r_grid_innercorner = (*(x1_buffer+i)) - 0.5 * (*(dx1_buffer+i));
               r_grid_outercorner = (*(x1_buffer+i)) + 0.5 * (*(dx1_buffer+i));

               theta_grid_innercorner = (*(x2_buffer+i)) - 0.5 * (*(dx2_buffer+i));
               theta_grid_outercorner = (*(x2_buffer+i)) + 0.5 * (*(dx2_buffer+i));
           #elif GEOMETRY == CARTESIAN
               r_grid_innercorner = pow((*(x1_buffer+i) - *(dx1_buffer+i)/2.0) * ((*(x1_buffer+i) - *(dx1_buffer+i)/2.0))+(*(x2_buffer+i) - *(dx2_buffer+i)/2.0) * (*(x2_buffer+i) - *(dx2_buffer+i)/2.0),0.5);
               r_grid_outercorner = pow((*(x1_buffer+i) + *(dx1_buffer+i)/2.0) * ((*(x1_buffer+i) + *(dx1_buffer+i)/2.0))+(*(x2_buffer+i) + *(dx2_buffer+i)/2.0) * (*(x2_buffer+i) + *(dx2_buffer+i)/2.0),0.5);

               theta_grid_innercorner = acos( (*(x2_buffer+i) - *(dx2_buffer+i)/2.0) /r_grid_innercorner); //arccos of y/r for the bottom left corner
               theta_grid_outercorner = acos( (*(x2_buffer+i) + *(dx2_buffer+i)/2.0) /r_grid_outercorner);
           #endif

            if (((ph_rmin - elem_factor*C_LIGHT/fps) <= r_grid_outercorner) && (r_grid_innercorner  <= (ph_rmax + elem_factor*C_LIGHT/fps) ) && (theta_grid_outercorner >= ph_thetamin) && (theta_grid_innercorner <= ph_thetamax))
            {
                (*pres)[j]=(*(pres_buffer+i));
                (*velx)[j]= (*(vel_x1_buffer+i))*sin((*(x2_buffer+i))) + (*(vel_x2_buffer+i))*cos((*(x2_buffer+i))) ; //convert from spherical to cartesian
                (*vely)[j]= (*(vel_x1_buffer+i))*cos((*(x2_buffer+i))) - (*(vel_x2_buffer+i))*sin((*(x2_buffer+i))); //z axis is actually y axis in 2D

                (*dens)[j]=(*(dens_buffer+i));
                (*x)[j]=(*(x1_buffer+i))*sin((*(x2_buffer+i)));
                (*y)[j]=(*(x1_buffer+i))*cos((*(x2_buffer+i)));
                (*r)[j]=(*(x1_buffer+i));
                (*szx)[j]=(*(dx1_buffer+i));
                (*szy)[j]=(*(dx2_buffer+i));
                (*theta)[j]=(*(x2_buffer+i));//theta in radians in relation to jet axis
                (*gamma)[j]=1/pow(1.0-( (*(vel_x1_buffer+i))*(*(vel_x1_buffer+i)) + (*(vel_x2_buffer+i))*(*(vel_x2_buffer+i)) ),0.5); //v is in units of c
                (*dens_lab)[j]= (*(dens_buffer+i)) /pow(1.0-( (*(vel_x1_buffer+i))*(*(vel_x1_buffer+i)) + (*(vel_x2_buffer+i))*(*(vel_x2_buffer+i)) ),0.5);
                (*temp)[j]=pow(3*(*(pres_buffer+i))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
                j++;
            }
        }
        else
        {
            if ((*(x1_buffer+i)) > (0.95*r_inj) )
            {
                (*pres)[j]=(*(pres_buffer+i));
                (*velx)[j]= (*(vel_x1_buffer+i))*sin((*(x2_buffer+i))) + (*(vel_x2_buffer+i))*cos((*(x2_buffer+i))) ; //convert from spherical to cartesian
                (*vely)[j]= (*(vel_x1_buffer+i))*cos((*(x2_buffer+i))) - (*(vel_x2_buffer+i))*sin((*(x2_buffer+i))); //z axis is actually y axis in 2D

                (*dens)[j]=(*(dens_buffer+i));
                (*x)[j]=(*(x1_buffer+i))*sin((*(x2_buffer+i)));
                (*y)[j]=(*(x1_buffer+i))*cos((*(x2_buffer+i)));
                (*r)[j]=(*(x1_buffer+i));
                (*szx)[j]=(*(dx1_buffer+i));
                (*szy)[j]=(*(dx2_buffer+i));
                (*theta)[j]=(*(x2_buffer+i));//theta in radians in relation to jet axis
                (*gamma)[j]=1/pow(1.0-( (*(vel_x1_buffer+i))*(*(vel_x1_buffer+i)) + (*(vel_x2_buffer+i))*(*(vel_x2_buffer+i)) ),0.5); //v is in units of c
                (*dens_lab)[j]= (*(dens_buffer+i)) /pow(1.0-( (*(vel_x1_buffer+i))*(*(vel_x1_buffer+i)) + (*(vel_x2_buffer+i))*(*(vel_x2_buffer+i)) ),0.5);
                (*temp)[j]=pow(3*(*(pres_buffer+i))*pow(C_LIGHT,2.0)/(A_RAD) ,1.0/4.0);
                j++;
            }
        }
    }

    //fprintf(fPtr, "number: %d\n", r_count);

    *number=r_count;

    free(x1_buffer); free(x2_buffer); free(dx1_buffer); free(dx2_buffer); free(dens_buffer); free(pres_buffer); free(vel_x1_buffer); free(vel_x2_buffer);

}

void modifyPlutoName(char file[200], char prefix[200], int frame)
{
    int lim1=0, lim2=0, lim3=0;
    char test[200]="" ;

    //if (strcmp(DIM_SWITCH, dim_2d_str)==0)
    #if DIMENSIONS == 2
    {
        //2D case
        lim1=10;
        lim2=100;
        lim3=1000;
    }
    #else
    {
        //3d case
        lim1=100;
        lim2=1000;
        lim3=10000;
    }
    #endif

    if (frame<lim1)
    {
        snprintf(test,sizeof(test), "%s%.3d%d.hdf5",prefix,000,frame);
    }
    else if (frame<lim2)
    {
        snprintf(test,sizeof(test), "%s%.2d%d.hdf5",prefix,00,frame);
    }
    else if (frame<lim3)
    {
        snprintf(test,sizeof(test), "%s%d%d.hdf5",prefix,0,frame);
    }
    else
    {
        snprintf(test,sizeof(test), "%s%d.hdf5",prefix,frame);
    }
    strncpy(file, test, sizeof(test));//had to do this workaround for some weird reason
}

ext: c

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT


text:
#pragma once
#include "utility/types.h"
#include <gsl/gsl>
#include <memory>

struct nk_context;
struct nk_user_font;
struct nk_buffer;

///
namespace cws80 {

class GraphicsDevice;
struct FontRequest;

//
class NkScreen {
public:
    NkScreen();
    ~NkScreen();

    void init(GraphicsDevice &gdev, uint w, uint h, gsl::span<const FontRequest> fontreqs);
    void clear();
    bool should_render() const;
    void render(void *draw_context);
    void update();
    nk_context *context() const;
    nk_user_font *font(uint id) const;

    uint width() const;
    uint height() const;

private:
    struct Impl;
    std::unique_ptr<Impl> P;
};

}  // namespace cws80


text:
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_gamma.h>

double
f(double x, void * params)
{
  int m = *(int *) params;
  double f = gsl_pow_int(x, m) + 1.0;
  return f;
}

int
main (int argc, char *argv[])
{
  gsl_integration_fixed_workspace * w;
  const gsl_integration_fixed_type * T = gsl_integration_fixed_hermite;
  int m = 10;
  int n = 6;
  double expected, result;
  gsl_function F;

  if (argc > 1)
    m = atoi(argv[1]);

  if (argc > 2)
    n = atoi(argv[2]);

  w = gsl_integration_fixed_alloc(T, n, 0.0, 1.0, 0.0, 0.0);

  F.function = &f;
  F.params = &m;

  gsl_integration_fixed(&F, &result, w);

  if (m % 2 == 0)
    expected = M_SQRTPI + gsl_sf_gamma(0.5*(1.0 + m));
  else
    expected = M_SQRTPI;

  printf ("m               = %d\n", m);
  printf ("intervals       = %zu\n", gsl_integration_fixed_n(w));
  printf ("result          = % .18f\n", result);
  printf ("exact result    = % .18f\n", expected);
  printf ("actual error    = % .18f\n", result - expected);

  gsl_integration_fixed_free (w);

  return 0;
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/* Copyright 2019-2020 Canaan Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#pragma once
#include "compiler_defs.h"
#include <cstring>
#include <gsl/gsl-lite.hpp>

BEGIN_NS_NNCASE_RUNTIME

class span_reader
{
public:
    span_reader(gsl::span<const gsl::byte> span)
        : span_(span)
    {
    }

    bool empty() const noexcept { return span_.empty(); }
    size_t avail() const noexcept { return span_.size_bytes(); }

    template <class T>
    T read()
    {
        auto value = *reinterpret_cast<const T *>(span_.data());
        advance(sizeof(T));
        return value;
    }

    template <class T>
    T read_unaligned()
    {
        alignas(T) uint8_t storage[sizeof(T)];
        std::memcpy(storage, span_.data(), sizeof(T));
        advance(sizeof(T));
        return *reinterpret_cast<const T *>(storage);
    }

    template <class T>
    void read(T &value)
    {
        value = *reinterpret_cast<const T *>(span_.data());
        advance(sizeof(T));
    }

    template <class T>
    void read_span(gsl::span<const T> &span, size_t size)
    {
        span = { reinterpret_cast<const T *>(span_.data()), size };
        advance(sizeof(T) * size);
    }

    template <class T = gsl::byte>
    gsl::span<const T> read_span(size_t size)
    {
        gsl::span<const T> span(reinterpret_cast<const T *>(span_.data()), size);
        advance(sizeof(T) * size);
        return span;
    }

    void read_avail(gsl::span<const gsl::byte> &span)
    {
        span = span_;
        span_ = {};
    }

    gsl::span<const gsl::byte> read_avail()
    {
        auto span = span_;
        span_ = {};
        return span;
    }

    gsl::span<const gsl::byte> peek_avail()
    {
        return span_;
    }

    template <class T>
    T peek()
    {
        auto value = *reinterpret_cast<const T *>(span_.data());
        return value;
    }

    template <class T>
    T peek_unaligned()
    {
        T value;
        std::memcpy(&value, span_.data(), sizeof(T));
        return value;
    }

    template <class T>
    T peek_unaligned_with_offset(size_t offset)
    {
        T value;
        std::memcpy(&value, span_.data() + offset, sizeof(T));
        return value;
    }

    template <class T>
    const T *get_ref()
    {
        auto ptr = reinterpret_cast<const T *>(span_.data());
        advance(sizeof(T));
        return ptr;
    }

    template <class T>
    void get_ref(const T *&ptr)
    {
        ptr = get_ref<T>();
    }

    void skip(size_t count)
    {
        advance(count);
    }

private:
    void advance(size_t count)
    {
        span_ = span_.subspan(count);
    }

private:
    gsl::span<const gsl::byte> span_;
};

END_NS_NNCASE_RUNTIME


text:
/*
 * Copyright 2009-2019 The VOTCA-MPIP Development Team (http://www.votca.org)
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 * author: Denis Andrienko
 */

#ifndef __VOTCA_CTP_EIGEN_H
#define __VOTCA_CTP_EIGEN_H

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <boost/numeric/conversion/cast.hpp>

#include <votca/tools/eigen.h>
#include <Eigen/Eigenvalues>

/*
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_eigen.h>
*/

#ifdef DEBUG_LINALG
#include <chrono>
#endif

namespace boost { namespace numeric { namespace ublas {

// specialization for double precision
template<class F, class A>
inline matrix<double,F,A>    // prod( m1, m2 )
prod(const matrix<double,F,A> &m1, const matrix<double,F,A> &m2)
{

    #ifdef DEBUG_LINALG
    auto start = std::chrono::system_clock::now();
    #endif

    boost::numeric::ublas::matrix<double,F,A> AxB( m1.size1(), m2.size2() );

    // Note that a Map<Type> m(&data) gives a view into the data (can change it)
    // while assignments Matrix<Type> m = Map< Matrix<Type> > (&data)
    // lead to a data copy and lost connection between the reference and matrix
    Eigen::Map<const Eigen::MatrixXd> mA_( &(m1.data()[0]), m1.size2(), m1.size1() );
    Eigen::Map<const Eigen::MatrixXd> mB_( &(m2.data()[0]), m2.size2(), m2.size1() );
    Eigen::Map<Eigen::MatrixXd> AxB_( &(AxB.data()[0]), AxB.size2(), AxB.size1() );

    AxB_ = mB_*mA_;

    #ifdef DEBUG_LINALG
    auto end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end - start;
    std::cout << "\n... ... ... \x1b[0;34mEigen "
              << "[" << m1.size1() << "x" << m1.size2() << "]"
              << "*[" << m2.size1() << "x" << m2.size2() << "]"
              << " time: " << elapsed_seconds.count()
              << "s\x1b[0;39m";
    #endif

    return AxB;

            }

// specialization for single precision
template<class F, class A>
inline matrix<float, F, A> // prod( m1, m2 )
prod(const matrix<float, F, A> &m1, const matrix<float, F, A> &m2) {

    #ifdef DEBUG_LINALG
    auto start = std::chrono::system_clock::now();
    #endif

    boost::numeric::ublas::matrix<float, F, A> AxB(m1.size1(), m2.size2());

    Eigen::Map<const Eigen::MatrixXf> mA_(&(m1.data()[0]), m1.size2(), m1.size1());
    Eigen::Map<const Eigen::MatrixXf> mB_(&(m2.data()[0]), m2.size2(), m2.size1());
    Eigen::Map<Eigen::MatrixXf> AxB_(&(AxB.data()[0]), AxB.size2(), AxB.size1());

    AxB_ = mB_*mA_;

    #ifdef DEBUG_LINALG
    auto end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end - start;
    std::cout << "\n... ... ... \x1b[0;34mEigen "
            << "[" << m1.size1() << "x" << m1.size2() << "]"
            << "*[" << m2.size1() << "x" << m2.size2() << "]"
            << " time: " << elapsed_seconds.count()
            << "s\x1b[0;39m";
    #endif

    return AxB;
}

}}}

namespace votca { namespace ctp { namespace linalg {

inline void invert_symm( const ub::matrix<double> &A, ub::matrix<double> &V){
    throw "Inversion of a matrix, linalg_invert, is not implemented in eigen.h";
}

/**
 * \brief eigenvalues of a symmetric matrix A*x=E*x
 * @param E vector of eigenvalues
 * @param V input: matrix to diagonalize
 * @param V output: eigenvectors
 *
 * This function wraps Eigen::EigenSolver
 *
 */
inline void eigenvalues_symm(const ub::matrix<double> &A,
                                   ub::vector<double> &E,
                                   ub::matrix<double> &V)
{

    #ifdef DEBUG_LINALG
        auto start = std::chrono::system_clock::now();
    #endif

/*
        gsl_error_handler_t *handler = gsl_set_error_handler_off();
	const size_t N = A.size1();

        V.resize(N, N, false);

        // gsl does not handle conversion of a symmetric_matrix
        ub::matrix<double> _A( N,N );
        //make copy of A as A is destroyed by GSL
        _A = A;

	E.resize(N, false);
	V.resize(N, N, false);
	gsl_matrix_view A_view = gsl_matrix_view_array(&_A(0,0), N, N);
	gsl_vector_view E_view = gsl_vector_view_array(&E(0), N);
	gsl_matrix_view V_view = gsl_matrix_view_array(&V(0,0), N, N);
	gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(N);

	int status = gsl_eigen_symmv(&A_view.matrix, &E_view.vector, &V_view.matrix, w);
        assert(status == 0);
	gsl_eigen_symmv_sort(&E_view.vector, &V_view.matrix, GSL_EIGEN_SORT_VAL_ASC);
	gsl_eigen_symmv_free(w);
	gsl_set_error_handler(handler);

        std::cout << "\nREFERENCE";
        std::cout << "\n A" << A;
        std::cout << "\n E" << E;
        std::cout << "\n V" << V;
        std::cout << "\nREFERENCE\n";

*/
    Eigen::Map<const Eigen::MatrixXd> A_( &(A.data()[0]), A.size2(), A.size1() );

    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(A_.transpose());
    /*
      std::cout << "\nEIGEN";
      std::cout << "\n E\n" << es.eigenvalues();
      std::cout << "\n V\n" << es.eigenvectors();
      std::cout << "\nEIGEN";
    */

    // copy back to UBLAS objects
    Eigen::Map<Eigen::VectorXd>( &E(0), A.size2() ) = es.eigenvalues();
    Eigen::Map<Eigen::MatrixXd>( &V(0,0), A.size2(), A.size1() ) = es.eigenvectors().transpose();

    #ifdef DEBUG_LINALG
        auto end = std::chrono::system_clock::now();

        std::chrono::duration<double> elapsed_seconds = end - start;

        std::cout << "\n... ... ... \x1b[0;34mEigen Eigenvalues"
                  << "[" << A.size1() << "x" << A.size2() << "]"
                  << " time: " << elapsed_seconds.count()
                  << "s\x1b[0;39m";
    #endif
}

}}}

#endif  // __VOTCA_CTP_EIGEN_H

max_stars_repo_name: jimbach/ctp

max_issues_repo_name: jimbach/ctp

max_forks_repo_name: jimbach/ctp


text:
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_math.h>
#include "gsl_cblas.h"

#include "tests.h"

void
test_her2 (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float A[] = { -0.821f, 0.954f };
   float X[] = { 0.532f, 0.802f };
   int incX = -1;
   float Y[] = { 0.016f, -0.334f };
   int incY = -1;
   float A_expected[] = { -0.302288f, 0.0f };
   cblas_cher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cher2(case 1450) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cher2(case 1450) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float A[] = { -0.821f, 0.954f };
   float X[] = { 0.532f, 0.802f };
   int incX = -1;
   float Y[] = { 0.016f, -0.334f };
   int incY = -1;
   float A_expected[] = { -0.302288f, 0.0f };
   cblas_cher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cher2(case 1451) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cher2(case 1451) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float A[] = { -0.821f, 0.954f };
   float X[] = { 0.532f, 0.802f };
   int incX = -1;
   float Y[] = { 0.016f, -0.334f };
   int incY = -1;
   float A_expected[] = { -0.302288f, 0.0f };
   cblas_cher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cher2(case 1452) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cher2(case 1452) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   float alpha[2] = {-1.0f, 0.0f};
   float A[] = { -0.821f, 0.954f };
   float X[] = { 0.532f, 0.802f };
   int incX = -1;
   float Y[] = { 0.016f, -0.334f };
   int incY = -1;
   float A_expected[] = { -0.302288f, 0.0f };
   cblas_cher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], flteps, "cher2(case 1453) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], flteps, "cher2(case 1453) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-0.3, 0.1};
   double A[] = { -0.334, 0.286 };
   double X[] = { -0.14, -0.135 };
   int incX = -1;
   double Y[] = { 0.455, 0.358 };
   int incY = -1;
   double A_expected[] = { -0.264521, 0.0 };
   cblas_zher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zher2(case 1454) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zher2(case 1454) imag");
     };
   };
  };

  {
   int order = 101;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-0.3, 0.1};
   double A[] = { -0.334, 0.286 };
   double X[] = { -0.14, -0.135 };
   int incX = -1;
   double Y[] = { 0.455, 0.358 };
   int incY = -1;
   double A_expected[] = { -0.264521, 0.0 };
   cblas_zher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zher2(case 1455) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zher2(case 1455) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 121;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-0.3, 0.1};
   double A[] = { -0.334, 0.286 };
   double X[] = { -0.14, -0.135 };
   int incX = -1;
   double Y[] = { 0.455, 0.358 };
   int incY = -1;
   double A_expected[] = { -0.264521, 0.0 };
   cblas_zher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zher2(case 1456) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zher2(case 1456) imag");
     };
   };
  };

  {
   int order = 102;
   int uplo = 122;
   int N = 1;
   int lda = 1;
   double alpha[2] = {-0.3, 0.1};
   double A[] = { -0.334, 0.286 };
   double X[] = { -0.14, -0.135 };
   int incX = -1;
   double Y[] = { 0.455, 0.358 };
   int incY = -1;
   double A_expected[] = { -0.264521, 0.0 };
   cblas_zher2(order, uplo, N, alpha, X, incX, Y, incY, A, lda);
   {
     int i;
     for (i = 0; i < 1; i++) {
       gsl_test_rel(A[2*i], A_expected[2*i], dbleps, "zher2(case 1457) real");
       gsl_test_rel(A[2*i+1], A_expected[2*i+1], dbleps, "zher2(case 1457) imag");
     };
   };
  };

}