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

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_permute_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_permutation.h>

static int
test_QR_decomp_dim(const gsl_matrix * m, double eps)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  gsl_matrix * qr = gsl_matrix_alloc(M, N);
  gsl_matrix * a  = gsl_matrix_alloc(M, N);
  gsl_matrix * q  = gsl_matrix_alloc(M, M);
  gsl_matrix * r  = gsl_matrix_alloc(M, N);
  gsl_vector * d = gsl_vector_alloc(N);

  gsl_matrix_memcpy(qr, m);

  s += gsl_linalg_QR_decomp(qr, d);
  s += gsl_linalg_QR_unpack(qr, d, q, r);

  /* compute a = q r */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, q, r, 0.0, a);

  for(i=0; i<M; i++) {
    for(j=0; j<N; j++) {
      double aij = gsl_matrix_get(a, i, j);
      double mij = gsl_matrix_get(m, i, j);
      int foo = check(aij, mij, eps);
      if(foo) {
        printf("(%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n", M, N, i,j, aij, mij);
      }
      s += foo;
    }
  }

  gsl_vector_free(d);
  gsl_matrix_free(qr);
  gsl_matrix_free(a);
  gsl_matrix_free(q);
  gsl_matrix_free(r);

  return s;
}

static int
test_QR_decomp(void)
{
  int f;
  int s = 0;

  f = test_QR_decomp_dim(m35, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp m(3,5)");
  s += f;

  f = test_QR_decomp_dim(m53, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp m(5,3)");
  s += f;

  f = test_QR_decomp_dim(hilb2, 2 * 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(2)");
  s += f;

  f = test_QR_decomp_dim(hilb3, 2 * 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(3)");
  s += f;

  f = test_QR_decomp_dim(hilb4, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(4)");
  s += f;

  f = test_QR_decomp_dim(hilb12, 2 * 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp hilbert(12)");
  s += f;

  f = test_QR_decomp_dim(vander2, 8.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(2)");
  s += f;

  f = test_QR_decomp_dim(vander3, 64.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(3)");
  s += f;

  f = test_QR_decomp_dim(vander4, 1024.0 * GSL_DBL_EPSILON);
  gsl_test(f, "  QR_decomp vander(4)");
  s += f;

  f = test_QR_decomp_dim(vander12, 0.0005); /* FIXME: bad accuracy */
  gsl_test(f, "  QR_decomp vander(12)");
  s += f;

  return s;
}

static int
test_QR_decomp_r_eps(const gsl_matrix * m, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = m->size1;
  const size_t N = m->size2;
  size_t i, j;

  gsl_matrix * QR = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * A  = gsl_matrix_alloc(M, N);
  gsl_matrix * R  = gsl_matrix_alloc(N, N);
  gsl_matrix * Q  = gsl_matrix_alloc(M, M);
  gsl_matrix_view Q1 = gsl_matrix_submatrix(Q, 0, 0, M, N);

  gsl_matrix_memcpy(QR, m);

  s += gsl_linalg_QR_decomp_r(QR, T);
  s += gsl_linalg_QR_unpack_r(QR, T, Q, R);

  /* compute A = Q R */
  gsl_matrix_memcpy(A, &Q1.matrix);
  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, A);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(A, i, j);
          double mij = gsl_matrix_get(m, i, j);

          gsl_test_rel(aij, mij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  if (M > N)
    {
      gsl_matrix * R_alt  = gsl_matrix_alloc(M, N);
      gsl_matrix * Q_alt  = gsl_matrix_alloc(M, M);
      gsl_vector_view tau = gsl_matrix_diagonal(T);

      /* test that Q2 was computed correctly by comparing with Level 2 algorithm */
      gsl_linalg_QR_unpack(QR, &tau.vector, Q_alt, R_alt);

      for (i = 0; i < M; i++)
        {
          for (j = 0; j < M; j++)
            {
              double aij = gsl_matrix_get(Q, i, j);
              double bij = gsl_matrix_get(Q_alt, i, j);

              gsl_test_rel(aij, bij, eps, "%s Q (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                           desc, M, N, i, j, aij, bij);
            }
        }

      gsl_matrix_free(R_alt);
      gsl_matrix_free(Q_alt);
    }

  gsl_matrix_free(QR);
  gsl_matrix_free(T);
  gsl_matrix_free(A);
  gsl_matrix_free(Q);
  gsl_matrix_free(R);

  return s;
}

static int
test_QR_decomp_r(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

  for (M = 1; M <= 50; ++M)
    {
      for (N = 1; N <= M; ++N)
        {
          gsl_matrix * A = gsl_matrix_alloc(M, N);

          create_random_matrix(A, r);
          s += test_QR_decomp_r_eps(A, 1.0e6 * M * GSL_DBL_EPSILON, "QR_decomp_r random");

          gsl_matrix_free(A);
        }
    }

  s += test_QR_decomp_r_eps(m53,   1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r m(5,3)");
  s += test_QR_decomp_r_eps(hilb2, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(2)");
  s += test_QR_decomp_r_eps(hilb3, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(3)");
  s += test_QR_decomp_r_eps(hilb4, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(4)");
  s += test_QR_decomp_r_eps(hilb12, 1.0e2 * GSL_DBL_EPSILON, "QR_decomp_r hilbert(12)");
  s += test_QR_decomp_r_eps(vander2, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(2)");
  s += test_QR_decomp_r_eps(vander3, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(3)");
  s += test_QR_decomp_r_eps(vander4, 1.0e1 * GSL_DBL_EPSILON, "QR_decomp_r vander(4)");

  return s;
}

static int
test_QR_QTmat_r_eps(const gsl_matrix * A, const gsl_matrix * B, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  const size_t K = B->size2;
  size_t i, j;

  gsl_matrix * QR = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * work  = gsl_matrix_alloc(N, K);
  gsl_matrix * B1 = gsl_matrix_alloc(M, K);
  gsl_matrix * B2 = gsl_matrix_alloc(M, K);
  gsl_vector_view tau = gsl_matrix_diagonal(T);

  gsl_matrix_memcpy(QR, A);
  gsl_matrix_memcpy(B1, B);
  gsl_matrix_memcpy(B2, B);

  s += gsl_linalg_QR_decomp_r(QR, T);

  /* compute Q^T B with both recursive and non-recursive methods and compare */
  s += gsl_linalg_QR_QTmat_r(QR, T, B1, work);
  s += gsl_linalg_QR_QTmat(QR, &tau.vector, B2);

  for (i = 0; i < M; i++)
    {
      for (j = 0; j < K; j++)
        {
          double aij = gsl_matrix_get(B1, i, j);
          double bij = gsl_matrix_get(B2, i, j);

          gsl_test_rel(aij, bij, eps, "%s (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, K, i,j, aij, bij);
        }
    }

  gsl_matrix_free(QR);
  gsl_matrix_free(T);
  gsl_matrix_free(B1);
  gsl_matrix_free(B2);
  gsl_matrix_free(work);

  return s;
}

static int
test_QR_QTmat_r(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

  for (M = 1; M <= 50; ++M)
    {
      for (N = 1; N <= M; ++N)
        {
          const size_t K = GSL_MAX(N / 2, 1);
          gsl_matrix * A = gsl_matrix_alloc(M, N);
          gsl_matrix * B = gsl_matrix_alloc(M, K);

          create_random_matrix(A, r);
          create_random_matrix(B, r);
          s += test_QR_QTmat_r_eps(A, B, 1.0e6 * M * GSL_DBL_EPSILON, "QR_QTmat_r random");

          gsl_matrix_free(A);
          gsl_matrix_free(B);
        }
    }

  return s;
}

static int
test_QR_solve_r_eps(const gsl_matrix * A, const gsl_vector * rhs, const gsl_vector * sol,
                    const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix * QR = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_vector * x  = gsl_vector_alloc(N);

  gsl_matrix_memcpy(QR, A);

  s += gsl_linalg_QR_decomp_r(QR, T);
  s += gsl_linalg_QR_solve_r(QR, T, rhs, x);

  for (i = 0; i < M; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(sol, i);

      gsl_test_rel(xi, yi, eps, "%s (%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, N, i, xi, yi);
    }

  gsl_matrix_free(QR);
  gsl_matrix_free(T);
  gsl_vector_free(x);

  return s;
}

static int
test_QR_solve_r(gsl_rng * r)
{
  int s = 0;
  size_t N;

  for (N = 1; N <= 50; ++N)
    {
      gsl_matrix * A = gsl_matrix_alloc(N, N);
      gsl_vector * sol = gsl_vector_alloc(N);
      gsl_vector * rhs = gsl_vector_alloc(N);

      create_random_matrix(A, r);
      create_random_vector(sol, r);
      gsl_blas_dgemv(CblasNoTrans, 1.0, A, sol, 0.0, rhs);

      s += test_QR_solve_r_eps(A, rhs, sol, 1.0e5 * N * GSL_DBL_EPSILON, "QR_solve_r random");

      gsl_matrix_free(A);
      gsl_vector_free(sol);
      gsl_vector_free(rhs);
    }

  return s;
}

static int
test_QR_lssolve_r_eps(const gsl_matrix * A, const gsl_vector * b, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i;

  gsl_matrix * QR = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_alloc(M);
  gsl_vector * work  = gsl_vector_alloc(N);

  gsl_matrix * U = gsl_matrix_alloc(M, N);
  gsl_matrix * V = gsl_matrix_alloc(N, N);
  gsl_vector * S = gsl_vector_alloc(N);
  gsl_vector * x_svd = gsl_vector_alloc(N);

  gsl_vector * residual = gsl_vector_alloc(M);

  gsl_matrix_memcpy(QR, A);
  s += gsl_linalg_QR_decomp_r(QR, T);
  s += gsl_linalg_QR_lssolve_r(QR, T, b, x, work);

  gsl_matrix_memcpy(U, A);
  gsl_linalg_SV_decomp(U, V, S, work);
  gsl_linalg_SV_solve(U, V, S, b, x_svd);

  /* compare QR with SVD solution */
  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(x_svd, i);

      gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, yi);
    }

  if (M > N)
    {
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
      gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
      double norm = gsl_blas_dnrm2(&x2.vector);
      double norm_expected;

      /* compute residual and check || x(N+1:end) || = || b - A x || */
      gsl_vector_memcpy(residual, b);
      gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
      norm_expected = gsl_blas_dnrm2(residual);

      gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, M, N, norm, norm_expected);
    }

  gsl_matrix_free(QR);
  gsl_matrix_free(T);
  gsl_vector_free(x);
  gsl_matrix_free(U);
  gsl_matrix_free(V);
  gsl_vector_free(x_svd);
  gsl_vector_free(work);
  gsl_vector_free(S);
  gsl_vector_free(residual);

  return s;
}

static int
test_QR_lssolve_r(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

  for (M = 1; M <= 30; ++M)
    {
      for (N = 1; N <= M; ++N)
        {
          gsl_matrix * A = gsl_matrix_alloc(M, N);
          gsl_vector * b = gsl_vector_alloc(M);

          create_random_matrix(A, r);
          create_random_vector(b, r);

          s += test_QR_lssolve_r_eps(A, b, 1.0e5 * M * GSL_DBL_EPSILON, "QR_lssolve_r random");

          gsl_matrix_free(A);
          gsl_vector_free(b);
        }
    }

  return s;
}

static int
test_QR_UR_decomp_eps(const gsl_matrix * S, const gsl_matrix * A, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i, j;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * V = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * B1  = gsl_matrix_calloc(N, N);
  gsl_matrix * B2  = gsl_matrix_alloc(M, N);

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, S);
  gsl_matrix_memcpy(V, A);

  s += gsl_linalg_QR_UR_decomp(R, V, T);

  /*
   * compute B = Q R = [ R - T R ]
   *                   [ -V~ T R ]
   */

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);
  gsl_matrix_memcpy(B2, V);

  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1);   /* B1 = T R */
  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, B1, B2); /* B2 = -V~ T R */
  gsl_matrix_sub (B1, R);                                                            /* B1 = T R - R */
  gsl_matrix_scale (B1, -1.0);                                                       /* B1 = R - T R */

  /* test S = B1 */
  for (i = 0; i < N; i++)
    {
      for (j = i; j < N; j++)
        {
          double aij = gsl_matrix_get(B1, i, j);
          double mij = gsl_matrix_get(S, i, j);

          gsl_test_rel(aij, mij, eps, "%s B1 (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  /* test A = B2 */
  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(B2, i, j);
          double mij = gsl_matrix_get(A, i, j);

          gsl_test_rel(aij, mij, eps, "%s B2 (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  gsl_matrix_free(R);
  gsl_matrix_free(V);
  gsl_matrix_free(T);
  gsl_matrix_free(B1);
  gsl_matrix_free(B2);

  return s;
}

static int
test_QR_UR_decomp(gsl_rng * r)
{
  int s = 0;
  size_t M, N;

  for (M = 1; M <= 50; ++M)
    {
      for (N = 1; N <= M; ++N)
        {
          gsl_matrix * S = gsl_matrix_alloc(N, N);
          gsl_matrix * A = gsl_matrix_alloc(M, N);

          create_random_matrix(A, r);
          create_random_matrix(S, r);
          s += test_QR_UR_decomp_eps(S, A, 1.0e6 * M * GSL_DBL_EPSILON, "QR_UR_decomp random");

          gsl_matrix_free(S);
          gsl_matrix_free(A);
        }
    }

  return s;
}

static int
test_QR_UZ_decomp_eps(const gsl_matrix * S, const gsl_matrix * A, const double eps, const char * desc)
{
  int s = 0;
  const size_t M = A->size1;
  const size_t N = A->size2;
  size_t i, j;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * V = gsl_matrix_alloc(M, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * B1  = gsl_matrix_calloc(N, N);
  gsl_matrix * B2  = gsl_matrix_alloc(M, N);

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, S);
  gsl_matrix_memcpy(V, A);

  s += gsl_linalg_QR_UZ_decomp(R, V, T);

  /*
   * compute B = Q R = [ R - T R ]
   *                   [ -V~ T R ]
   */

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);

  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1);   /* B1 = T R */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, -1.0, V, B1, 0.0, B2);                 /* B2 = -V~ T R */
  gsl_matrix_sub (B1, R);                                                            /* B1 = T R - R */
  gsl_matrix_scale (B1, -1.0);                                                       /* B1 = R - T R */

  /* test S = B1 */
  for (i = 0; i < N; i++)
    {
      for (j = i; j < N; j++)
        {
          double aij = gsl_matrix_get(B1, i, j);
          double mij = gsl_matrix_get(S, i, j);

          gsl_test_rel(aij, mij, eps, "%s B1 (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  /* test A = B2 */
  for (i = 0; i < M; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(B2, i, j);
          double mij = gsl_matrix_get(A, i, j);

          gsl_test_rel(aij, mij, eps, "%s B2 (%3lu,%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, M, N, i,j, aij, mij);
        }
    }

  gsl_matrix_free(R);
  gsl_matrix_free(V);
  gsl_matrix_free(T);
  gsl_matrix_free(B1);
  gsl_matrix_free(B2);

  return s;
}

static int
test_QR_UZ_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 20;
  const size_t M_max = 2*N_max;
  gsl_matrix * U1 = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * U2 = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * A = gsl_matrix_alloc(M_max, N_max);
  size_t M, N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view S = gsl_matrix_submatrix(U1, 0, 0, N, N);
      gsl_matrix_view T = gsl_matrix_submatrix(U2, 0, 0, N, N);

      for (M = N; M <= 2*N_max; ++M)
        {
          gsl_matrix_view B = gsl_matrix_submatrix(A, 0, 0, M, N);
          gsl_matrix_view Bu = gsl_matrix_submatrix(&B.matrix, M - N, 0, N, N);

          gsl_matrix_set_zero(&B.matrix);

          if (M > N)
            {
              gsl_matrix_view Bd = gsl_matrix_submatrix(&B.matrix, 0, 0, M - N, N);
              create_random_matrix(&Bd.matrix, r);
            }

          create_random_matrix(&S.matrix, r);
          create_random_matrix(&T.matrix, r);
          gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &Bu.matrix, &T.matrix);

          s += test_QR_UZ_decomp_eps(&S.matrix, &B.matrix, 1.0e6 * M * GSL_DBL_EPSILON,
                                     "QR_UZ_decomp random");
        }
    }

  gsl_matrix_free(U1);
  gsl_matrix_free(U2);
  gsl_matrix_free(A);

  return s;
}

static int
test_QR_UU_decomp_eps(const gsl_matrix * U, const gsl_matrix * S, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = U->size2;
  size_t i, j;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * V = gsl_matrix_calloc(N, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * B1  = gsl_matrix_calloc(N, N);
  gsl_matrix * B2  = gsl_matrix_calloc(N, N);

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, V, S);

  s += gsl_linalg_QR_UU_decomp(R, V, T);

  /*
   * compute B = Q R = [ R - T R ]
   *                   [ -V~ T R ]
   */

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);

  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1);   /* B1 = T R */
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B2, B1);                               /* B2 := T R */
  gsl_blas_dtrmm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, V, B2);   /* B2 = -V~ T R */
  gsl_matrix_sub (B1, R);                                                            /* B1 = T R - R */
  gsl_matrix_scale (B1, -1.0);                                                       /* B1 = R - T R */

  /* test U = B1 */
  for (i = 0; i < N; i++)
    {
      for (j = i; j < N; j++)
        {
          double aij = gsl_matrix_get(B1, i, j);
          double mij = gsl_matrix_get(U, i, j);

          gsl_test_rel(aij, mij, eps, "%s B1 (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i, j, aij, mij);
        }
    }

  /* test S = B2 */
  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(B2, i, j);
          double mij = gsl_matrix_get(S, i, j);

          gsl_test_rel(aij, mij, eps, "%s B2 (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i,j, aij, mij);
        }
    }

  gsl_matrix_free(R);
  gsl_matrix_free(V);
  gsl_matrix_free(T);
  gsl_matrix_free(B1);
  gsl_matrix_free(B2);

  return s;
}

static int
test_QR_UU_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 50;
  gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * T = gsl_matrix_alloc(N_max, N_max);
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
      gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
      gsl_matrix_view c = gsl_matrix_submatrix(T, 0, 0, N, N);

      create_random_matrix(&c.matrix, r);
      gsl_matrix_set_zero(&a.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &c.matrix);

      create_random_matrix(&c.matrix, r);
      gsl_matrix_set_zero(&b.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &b.matrix, &c.matrix);

      s += test_QR_UU_decomp_eps(&a.matrix, &b.matrix, 1.0e6 * N * GSL_DBL_EPSILON,
                                 "QR_UU_decomp random");
    }

  gsl_matrix_free(U);
  gsl_matrix_free(S);
  gsl_matrix_free(T);

  return s;
}

static int
test_QR_UU_lssolve_eps(const gsl_matrix * U, const gsl_matrix * S, const gsl_vector * b,
                       const double eps, const char * desc)
{
  int s = 0;
  const size_t N = U->size1;
  const size_t M = 2 * N;
  size_t i;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * Y = gsl_matrix_calloc(N, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_alloc(M);
  gsl_vector * work  = gsl_vector_alloc(N);

  gsl_matrix * A = gsl_matrix_calloc(M, N);
  gsl_matrix * U_svd = gsl_matrix_alloc(M, N);
  gsl_matrix * V_svd = gsl_matrix_alloc(N, N);
  gsl_vector * S_svd = gsl_vector_alloc(N);
  gsl_vector * x_svd = gsl_vector_alloc(N);

  gsl_vector * residual = gsl_vector_alloc(M);

  gsl_matrix_view tmp;

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, Y, S);

  s += gsl_linalg_QR_UU_decomp(R, Y, T);
  s += gsl_linalg_QR_UU_lssolve(R, Y, T, b, x, work);

  tmp = gsl_matrix_submatrix(A, 0, 0, N, N);
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, U);

  tmp = gsl_matrix_submatrix(A, N, 0, N, N);
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, S);

  gsl_matrix_memcpy(U_svd, A);
  gsl_linalg_SV_decomp(U_svd, V_svd, S_svd, work);
  gsl_linalg_SV_solve(U_svd, V_svd, S_svd, b, x_svd);

  /* compare QR with SVD solution */
  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(x_svd, i);

      gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, yi);
    }

  if (M > N)
    {
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
      gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
      double norm = gsl_blas_dnrm2(&x2.vector);
      double norm_expected;

      /* compute residual and check || x(N+1:end) || = || b - A x || */
      gsl_vector_memcpy(residual, b);
      gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
      norm_expected = gsl_blas_dnrm2(residual);

      gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, M, N, norm, norm_expected);
    }

  gsl_matrix_free(A);
  gsl_matrix_free(R);
  gsl_matrix_free(Y);
  gsl_matrix_free(T);
  gsl_vector_free(x);
  gsl_matrix_free(U_svd);
  gsl_matrix_free(V_svd);
  gsl_vector_free(x_svd);
  gsl_vector_free(work);
  gsl_vector_free(S_svd);
  gsl_vector_free(residual);

  return s;
}

static int
test_QR_UU_lssolve(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 30;
  gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * T = gsl_matrix_alloc(N_max, N_max);
  gsl_vector * rhs = gsl_vector_alloc(2*N_max);
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
      gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
      gsl_matrix_view c = gsl_matrix_submatrix(T, 0, 0, N, N);
      gsl_vector_view rhsv = gsl_vector_subvector(rhs, 0, 2*N);

      create_random_vector(&rhsv.vector, r);

      create_random_matrix(&c.matrix, r);
      gsl_matrix_set_zero(&a.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &c.matrix);

      create_random_matrix(&c.matrix, r);
      gsl_matrix_set_zero(&b.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &b.matrix, &c.matrix);

      s += test_QR_UU_lssolve_eps(&a.matrix, &b.matrix, &rhsv.vector, 1.0e4 * N * GSL_DBL_EPSILON,
                                  "QR_UU_lssolve random");
    }

  gsl_matrix_free(U);
  gsl_matrix_free(S);
  gsl_matrix_free(T);
  gsl_vector_free(rhs);

  return s;
}

static int
test_QR_UD_decomp_eps(const gsl_matrix * U, const gsl_vector * D, const double eps, const char * desc)
{
  int s = 0;
  const size_t N = U->size2;
  size_t i, j;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * V = gsl_matrix_calloc(N, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_matrix * B1  = gsl_matrix_calloc(N, N);
  gsl_matrix * B2  = gsl_matrix_calloc(N, N);

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);

  s += gsl_linalg_QR_UD_decomp(R, D, V, T);

  /*
   * compute B = Q R = [ R - T R ]
   *                   [ -V~ T R ]
   */

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B1, T);

  gsl_blas_dtrmm (CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0, R, B1);   /* B1 = T R */
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, B2, B1);                               /* B2 := T R */
  gsl_blas_dtrmm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, -1.0, V, B2);   /* B2 = -V~ T R */
  gsl_matrix_sub (B1, R);                                                            /* B1 = T R - R */
  gsl_matrix_scale (B1, -1.0);                                                       /* B1 = R - T R */

  /* test U = B1 */
  for (i = 0; i < N; i++)
    {
      for (j = i; j < N; j++)
        {
          double aij = gsl_matrix_get(B1, i, j);
          double mij = gsl_matrix_get(U, i, j);

          gsl_test_rel(aij, mij, eps, "%s B1 (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i, j, aij, mij);
        }
    }

  /* test S = B2 */
  for (i = 0; i < N; i++)
    {
      for (j = 0; j < N; j++)
        {
          double aij = gsl_matrix_get(B2, i, j);
          double mij = (i == j) ? gsl_vector_get(D, i) : 0.0;

          gsl_test_rel(aij, mij, eps, "%s B2 (%3lu)[%lu,%lu]: %22.18g   %22.18g\n",
                       desc, N, i,j, aij, mij);
        }
    }

  gsl_matrix_free(R);
  gsl_matrix_free(V);
  gsl_matrix_free(T);
  gsl_matrix_free(B1);
  gsl_matrix_free(B2);

  return s;
}

static int
test_QR_UD_decomp(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 30;
  gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
  gsl_vector * D = gsl_vector_alloc(N_max);
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
      gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
      gsl_vector_view diag = gsl_vector_subvector(D, 0, N);

      create_random_matrix(&b.matrix, r);
      gsl_matrix_set_zero(&a.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &b.matrix);

      create_random_vector(&diag.vector, r);

      s += test_QR_UD_decomp_eps(&a.matrix, &diag.vector, 1.0e4 * N * GSL_DBL_EPSILON,
                                 "QR_UD_decomp random");
    }

  gsl_matrix_free(U);
  gsl_matrix_free(S);
  gsl_vector_free(D);

  return s;
}

static int
test_QR_UD_lssolve_eps(const gsl_matrix * U, const gsl_vector * D, const gsl_vector * b,
                       const double eps, const char * desc)
{
  int s = 0;
  const size_t N = U->size1;
  const size_t M = 2 * N;
  size_t i;

  gsl_matrix * R = gsl_matrix_calloc(N, N);
  gsl_matrix * Y = gsl_matrix_calloc(N, N);
  gsl_matrix * T = gsl_matrix_alloc(N, N);
  gsl_vector * x = gsl_vector_alloc(M);
  gsl_vector * work  = gsl_vector_alloc(N);

  gsl_matrix * A = gsl_matrix_calloc(M, N);
  gsl_matrix * U_svd = gsl_matrix_alloc(M, N);
  gsl_matrix * V_svd = gsl_matrix_alloc(N, N);
  gsl_vector * S_svd = gsl_vector_alloc(N);
  gsl_vector * x_svd = gsl_vector_alloc(N);

  gsl_vector * residual = gsl_vector_alloc(M);

  gsl_matrix_view tmp;
  gsl_vector_view diag;

  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, R, U);

  s += gsl_linalg_QR_UD_decomp(R, D, Y, T);
  s += gsl_linalg_QR_UD_lssolve(R, Y, T, b, x, work);

  tmp = gsl_matrix_submatrix(A, 0, 0, N, N);
  gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &tmp.matrix, U);

  tmp = gsl_matrix_submatrix(A, N, 0, N, N);
  diag = gsl_matrix_diagonal(&tmp.matrix);
  gsl_vector_memcpy(&diag.vector, D);

  gsl_matrix_memcpy(U_svd, A);
  gsl_linalg_SV_decomp(U_svd, V_svd, S_svd, work);
  gsl_linalg_SV_solve(U_svd, V_svd, S_svd, b, x_svd);

  /* compare QR with SVD solution */
  for (i = 0; i < N; i++)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(x_svd, i);

      gsl_test_rel(xi, yi, eps, "%s (%3lu,%3lu)[%lu]: %22.18g   %22.18g\n",
                   desc, M, N, i, xi, yi);
    }

  if (M > N)
    {
      gsl_vector_view x1 = gsl_vector_subvector(x, 0, N);
      gsl_vector_view x2 = gsl_vector_subvector(x, N, M - N);
      double norm = gsl_blas_dnrm2(&x2.vector);
      double norm_expected;

      /* compute residual and check || x(N+1:end) || = || b - A x || */
      gsl_vector_memcpy(residual, b);
      gsl_blas_dgemv(CblasNoTrans, -1.0, A, &x1.vector, 1.0, residual);
      norm_expected = gsl_blas_dnrm2(residual);

      gsl_test_rel(norm, norm_expected, eps, "%s rnorm (%3lu,%3lu): %22.18g   %22.18g\n",
                   desc, M, N, norm, norm_expected);
    }

  gsl_matrix_free(A);
  gsl_matrix_free(R);
  gsl_matrix_free(Y);
  gsl_matrix_free(T);
  gsl_vector_free(x);
  gsl_matrix_free(U_svd);
  gsl_matrix_free(V_svd);
  gsl_vector_free(x_svd);
  gsl_vector_free(work);
  gsl_vector_free(S_svd);
  gsl_vector_free(residual);

  return s;
}

static int
test_QR_UD_lssolve(gsl_rng * r)
{
  int s = 0;
  const size_t N_max = 30;
  gsl_matrix * U = gsl_matrix_alloc(N_max, N_max);
  gsl_matrix * S = gsl_matrix_alloc(N_max, N_max);
  gsl_vector * D = gsl_vector_alloc(N_max);
  gsl_vector * rhs = gsl_vector_alloc(2*N_max);
  size_t N;

  for (N = 1; N <= N_max; ++N)
    {
      gsl_matrix_view a = gsl_matrix_submatrix(U, 0, 0, N, N);
      gsl_matrix_view b = gsl_matrix_submatrix(S, 0, 0, N, N);
      gsl_vector_view rhsv = gsl_vector_subvector(rhs, 0, 2*N);
      gsl_vector_view diag = gsl_vector_subvector(D, 0, N);

      create_random_vector(&rhsv.vector, r);
      create_random_vector(&diag.vector, r);

      create_random_matrix(&b.matrix, r);
      gsl_matrix_set_zero(&a.matrix);
      gsl_matrix_tricpy(CblasUpper, CblasNonUnit, &a.matrix, &b.matrix);

      s += test_QR_UD_lssolve_eps(&a.matrix, &diag.vector, &rhsv.vector, 1.0e4 * N * GSL_DBL_EPSILON,
                                  "QR_UD_lssolve random");
    }

  gsl_matrix_free(U);
  gsl_matrix_free(S);
  gsl_vector_free(D);
  gsl_vector_free(rhs);

  return s;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/**
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You 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.
 */
#ifndef LIBMINIFI_INCLUDE_UTILS_GSL_H_
#define LIBMINIFI_INCLUDE_UTILS_GSL_H_

#include <type_traits>

#include <gsl-lite/gsl-lite.hpp>

namespace org {
namespace apache {
namespace nifi {
namespace minifi {

namespace gsl = ::gsl_lite;

namespace utils {
namespace detail {
template<typename T>
using remove_cvref_t = typename std::remove_cv<typename std::remove_reference<T>::type>::type;
}  // namespace detail

template<typename Container, typename T>
Container span_to(gsl::span<T> span) {
  static_assert(std::is_constructible<Container, typename gsl::span<T>::iterator, typename gsl::span<T>::iterator>::value,
      "The destination container must have an iterator (pointer) range constructor");
  return Container(std::begin(span), std::end(span));
}
template<template<typename...> class Container, typename T>
Container<detail::remove_cvref_t<T>> span_to(gsl::span<T> span) {
  static_assert(std::is_constructible<Container<detail::remove_cvref_t<T>>, typename gsl::span<T>::iterator, typename gsl::span<T>::iterator>::value,
      "The destination container must have an iterator (pointer) range constructor");
  return span_to<Container<detail::remove_cvref_t<T>>>(span);
}
}  // namespace utils

}  // namespace minifi
}  // namespace nifi
}  // namespace apache
}  // namespace org

#endif  // LIBMINIFI_INCLUDE_UTILS_GSL_H_

ext: h

max_stars_repo_name: dtrodrigues/nifi-minifi-cpp

max_issues_repo_name: dtrodrigues/nifi-minifi-cpp

max_forks_repo_name: dtrodrigues/nifi-minifi-cpp


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

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_exp.h>
#include <gsl/gsl_sf_bessel.h>

#include "error.h"

#include "chebyshev.h"
#include "cheb_eval.c"

/*-*-*-*-*-*-*-*-*-*-*-* Private Section *-*-*-*-*-*-*-*-*-*-*-*/

/* based on SLATEC besk1(), besk1e() */

/* chebyshev expansions

 series for bk1        on the interval  0.          to  4.00000d+00
                                        with weighted error   7.02e-18
                                         log weighted error  17.15
                               significant figures required  16.73
                                    decimal places required  17.67

 series for ak1        on the interval  1.25000d-01 to  5.00000d-01
                                        with weighted error   6.06e-17
                                         log weighted error  16.22
                               significant figures required  15.41
                                    decimal places required  16.83

 series for ak12       on the interval  0.          to  1.25000d-01
                                        with weighted error   2.58e-17
                                         log weighted error  16.59
                               significant figures required  15.22
                                    decimal places required  17.16
*/

static double bk1_data[11] = {
   0.0253002273389477705,
  -0.3531559607765448760,
  -0.1226111808226571480,
  -0.0069757238596398643,
  -0.0001730288957513052,
  -0.0000024334061415659,
  -0.0000000221338763073,
  -0.0000000001411488392,
  -0.0000000000006666901,
  -0.0000000000000024274,
  -0.0000000000000000070
};

static cheb_series bk1_cs = {
  bk1_data,
  10,
  -1, 1,
  8
};

static double ak1_data[17] = {
   0.27443134069738830,
   0.07571989953199368,
  -0.00144105155647540,
   0.00006650116955125,
  -0.00000436998470952,
   0.00000035402774997,
  -0.00000003311163779,
   0.00000000344597758,
  -0.00000000038989323,
   0.00000000004720819,
  -0.00000000000604783,
   0.00000000000081284,
  -0.00000000000011386,
   0.00000000000001654,
  -0.00000000000000248,
   0.00000000000000038,
  -0.00000000000000006
};
static cheb_series ak1_cs = {
  ak1_data,
  16,
  -1, 1,
  9
};

static double ak12_data[14] = {
   0.06379308343739001,
   0.02832887813049721,
  -0.00024753706739052,
   0.00000577197245160,
  -0.00000020689392195,
   0.00000000973998344,
  -0.00000000055853361,
   0.00000000003732996,
  -0.00000000000282505,
   0.00000000000023720,
  -0.00000000000002176,
   0.00000000000000215,
  -0.00000000000000022,
   0.00000000000000002
};
static cheb_series ak12_cs = {
  ak12_data,
  13,
  -1, 1,
  7
};

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int gsl_sf_bessel_K1_scaled_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*GSL_DBL_MIN) {
    OVERFLOW_ERROR(result);
  }
  else if(x <= 2.0) {
    const double lx = log(x);
    const double ex = exp(x);
    int stat_I1;
    gsl_sf_result I1;
    gsl_sf_result c;
    cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
    stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
    result->val  = ex * ((lx-M_LN2)*I1.val + (0.75 + c.val)/x);
    result->err  = ex * (c.err/x + fabs(lx)*I1.err);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_I1;
  }
  else if(x <= 8.0) {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak1_cs, (16.0/x-5.0)/3.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    const double sx = sqrt(x);
    gsl_sf_result c;
    cheb_eval_e(&ak12_cs, 16.0/x-1.0, &c);
    result->val  = (1.25 + c.val) / sx;
    result->err  = c.err / sx;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
}

int gsl_sf_bessel_K1_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x <= 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x < 2.0*GSL_DBL_MIN) {
    OVERFLOW_ERROR(result);
  }
  else if(x <= 2.0) {
    const double lx = log(x);
    int stat_I1;
    gsl_sf_result I1;
    gsl_sf_result c;
    cheb_eval_e(&bk1_cs, 0.5*x*x-1.0, &c);
    stat_I1 = gsl_sf_bessel_I1_e(x, &I1);
    result->val  = (lx-M_LN2)*I1.val + (0.75 + c.val)/x;
    result->err  = c.err/x + fabs(lx)*I1.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_I1;
  }
  else {
    gsl_sf_result K1_scaled;
    int stat_K1 = gsl_sf_bessel_K1_scaled_e(x, &K1_scaled);
    int stat_e  = gsl_sf_exp_mult_err_e(-x, 0.0,
                                           K1_scaled.val, K1_scaled.err,
                                           result);
    result->err = fabs(result->val) * (GSL_DBL_EPSILON*fabs(x) + K1_scaled.err/K1_scaled.val);
    return GSL_ERROR_SELECT_2(stat_e, stat_K1);
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_bessel_K1_scaled(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K1_scaled_e(x, &result));
}

double gsl_sf_bessel_K1(const double x)
{
  EVAL_RESULT(gsl_sf_bessel_K1_e(x, &result));
}

ext: c

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
/*
** CCM - connectivity consistency mapping
**
** G.Lohmann, July 2010
*/
#include <viaio/Vlib.h>
#include <viaio/file.h>
#include <viaio/mu.h>
#include <viaio/option.h>
#include <viaio/os.h>
#include <viaio/VImage.h>

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

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

#define SQR(x) ((x) * (x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

extern double VOCCC(gsl_matrix *,int,int);
extern double VICC(gsl_matrix *,int);
extern double VKendall_W(gsl_matrix *,int,int);
extern double SpearmanCorr(const float *data1,const float *data2, int n);

int metric = 1;
int verbose = 0;

VDictEntry TypeDict[] = {
  { "kendall", 0 },
  { "occc", 1 },
  { NULL }
};

double Correlation(const float *arr1,const float *arr2,int n)
{
  int i;
  double sx,sy,sxx,syy,sxy,rx;
  double tiny=1.0e-4;

  sxx = syy = sxy = sx = sy = 0;
  for (i=0; i<n; i++) {
    const double u = (double)arr1[i];
    const double v = (double)arr2[i];
    sx  += u;
    sy  += v;
    sxx += u*u;
    syy += v*v;
    sxy += u*v;
  }
  const double nx= n;
  const double u = nx*sxx - sx*sx;
  const double v = nx*syy - sy*sy;
  rx = 0;
  if (u*v > tiny)
    rx = (nx*sxy - sx*sy)/sqrt(u*v);
  return rx;
}

VImage VCCM(VImage **src, VImage mask, int nsubjects, int nslices,int first,int length,int type,int otype)
{
  VImage map=NULL,dest=NULL;
  int i,j,k,s,nvoxels,nt,len,last,b,r,c,bb,rr,cc,rad2,nrows,ncols;
  double u=0,v=0;
  gsl_matrix *data=NULL;
  gsl_matrix_float **mat=NULL;

  nrows = VImageNRows(mask);
  ncols = VImageNColumns(mask);

  /* count number of voxels, number of timesteps */
  nvoxels = nt = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (VGetPixel(mask,b,r,c) < 0.001) continue;
	if (VImageNBands(src[0][b]) > nt)
	  nt = VImageNBands(src[0][b]);
	nvoxels++;
      }
    }
  }

  /* get shortest length of time series across subjects */
  nt = 99999;
  for (s=0; s<nsubjects; s++) {
    for (b=0; b<nslices; b++) {
      for (r=0; r<nrows; r++) {
	for (c=0; c<ncols; c++) {
	  if (VGetPixel(mask,b,r,c) < 0.001) continue;
	  if (VImageNBands(src[s][b]) < nt)
	    nt = VImageNBands(src[s][b]);
	}
      }
    }
  }

  /* get time steps to include */
  if (length < 1) length = nt-2;
  last = first + length -1;
  if (last >= nt) last = nt-1;
  if (first < 0) first = 1;

  nt = last - first + 1;
  if (nt < 2) VError(" not enough timesteps, nt= %d",nt);

  fprintf(stderr," first= %d, last= %d, ntimesteps= %d\n",
          (int)first,(int)last,(int)nt);
  fprintf(stderr," nsubjects= %d,  nvoxels= %d\n",nsubjects,nvoxels);
  fprintf(stderr," type= %s\n",TypeDict[type].keyword);

  /*
  ** store voxel addresses
  */
  map = VCreateImage(1,3,nvoxels,VIntegerRepn);
  if (map == NULL) VError(" error allocating addr map");
  VFillImage(map,VAllBands,0);

  i = 0;
  for (b=0; b<nslices; b++) {
    for (r=0; r<nrows; r++) {
      for (c=0; c<ncols; c++) {
	if (VGetPixel(mask,b,r,c) < 0.001) continue;
	if (i >= nvoxels) VError(" nvox= %d",i);
	VPixel(map,0,0,i,VInteger) = b;
	VPixel(map,0,1,i,VInteger) = r;
	VPixel(map,0,2,i,VInteger) = c;
	i++;
      }
    }
  }

  /*
  ** avoid casting to float, copy data to matrix
  */
  mat = (gsl_matrix_float **) VCalloc(nsubjects,sizeof(gsl_matrix_float *));
  if (!mat) VError(" err allocating mat");

  for (s=0; s<nsubjects; s++) {
    mat[s] = gsl_matrix_float_calloc(nvoxels,nt);
    if (!mat[s]) VError(" err allocating mat[s]");

    for (i=0; i<nvoxels; i++) {
      b = VPixel(map,0,0,i,VInteger);
      r = VPixel(map,0,1,i,VInteger);
      c = VPixel(map,0,2,i,VInteger);

      float *ptr = gsl_matrix_float_ptr(mat[s],i,0);
      int k;
      for (k=first; k<=last; k++) {
	*ptr++ = (float) VGetPixel(src[s][b],k,r,c);
      }
    }
  }

  /*
  ** for each seed voxel assess intersubject consistency
  */
  dest = VCreateImage(nslices,nrows,ncols,VFloatRepn);
  VFillImage(dest,VAllBands,0);
  VCopyImageAttrs (src[0][0], dest);

  /*
  ** get correlation maps
  */
  data = gsl_matrix_calloc(nsubjects,nvoxels);

  rad2 = 3*9;
  for (i=0; i<nvoxels; i++) {
    if (i%10 == 0) fprintf(stderr," i: %7d  of  %d\r",i,nvoxels);
    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);

    len = 0;
    for (s=0; s<nsubjects; s++) {
      const float *arr1 = gsl_matrix_float_ptr(mat[s],i,0);

      k = 0;
      for (j=0; j<nvoxels; j++) {
	if (i == j) continue;
	bb = VPixel(map,0,0,j,VInteger);
	rr = VPixel(map,0,1,j,VInteger);
	cc = VPixel(map,0,2,j,VInteger);
	/* exclude because of spatial smoothness: */
	if ((SQR(b-bb) + SQR(r-rr) + SQR(c-cc)) < rad2) continue;

	const float *arr2 = gsl_matrix_float_ptr(mat[s],j,0);
	u = Correlation(arr1,arr2,nt);

	gsl_matrix_set(data,s,k,(double)u);
	k++;
      }
      len = k;
    }
    if (len >= data->size2) VError(" len= %d %d",len,data->size2);

    if (type == 0)
      v = VKendall_W(data,len,otype);
    else if (type == 1)
      v = VOCCC(data,len,verbose);
    else
      VError(" illegal type");

    b = VPixel(map,0,0,i,VInteger);
    r = VPixel(map,0,1,i,VInteger);
    c = VPixel(map,0,2,i,VInteger);
    VPixel(dest,b,r,c,VFloat) = v;
  }
  fprintf(stderr,"\n");
  return dest;
}

int main (int argc, char *argv[])
{
  static VArgVector in_files;
  static VString out_filename;
  static VString mask_filename;
  static VShort type   = 0;
  static VShort otype  = 0;
  static VShort first  = 2;
  static VShort length = 0;
  static VOptionDescRec options[] = {
    {"in", VStringRepn, 0, & in_files, VRequiredOpt, NULL,"Input files" },
    {"out", VStringRepn, 1, & out_filename, VRequiredOpt, NULL,"Output file" },
    {"mask",VStringRepn,1,(VPointer) &mask_filename,VRequiredOpt,NULL,"Mask file"},
    {"first",VShortRepn,1,(VPointer) &first,VOptionalOpt,NULL,"First timestep to use"},
    {"length",VShortRepn,1,(VPointer) &length,VOptionalOpt,NULL,
     "Length of time series to use, '0' to use full length"},
    {"type",VShortRepn,1,(VPointer) &type,VOptionalOpt,TypeDict,"Type of metric"},
    {"gauss",VShortRepn,1,(VPointer) &otype,VOptionalOpt,NULL,"Use gaussian version of kendall"},
  };
  VString in_filename;
  VAttrList list=NULL,geolist=NULL;
  int b,r,c,i=0;
  float u;
  char *prg = GetLipsiaName("vccm");

  /*
  ** parse command line
  */
  if (! VParseCommand (VNumber (options), options, & argc, argv)) {
    VReportUsage (argv[0], VNumber (options), options, NULL);
    exit (EXIT_FAILURE);
  }
  if (argc > 1) {
    VReportBadArgs (argc, argv);
    exit (EXIT_FAILURE);
  }

  /*
  ** read mask
  */
  VAttrList mask_list = VReadAttrList(mask_filename,0L,TRUE,FALSE);
  if (mask_list == NULL) VError(" error reading %s",mask_filename);
  VImage mask = VReadImage(mask_list);
  if (mask == NULL) VError(" no mask found");

  /*
  ** read input images
  */
  int nslices=0,nt=0,nrows=0,ncols=0;
  int xslices=0,xnt=0,xrows=0,xcols=0;
  int nsubjects = (int)in_files.number;
  fprintf(stderr," nsubjects= %d\n",nsubjects);
  if (nsubjects < 2) VError(" not enough input files (%d), CCM should be used on >= 2 data sets",nsubjects);
  VImage **src = (VImage **) VCalloc(nsubjects,sizeof(VImage *));

  for (i=0; i<nsubjects; i++) {
    in_filename = ((VString *) in_files.vector)[i];
    list    = VReadAttrList(in_filename,0L,TRUE,FALSE);
    xslices = VAttrListNumImages(list);
    if (i==0) nslices = xslices;
    else if (nslices != xslices) VError(" inconsistent dimensions in input files");
    src[i] = VAttrListGetImages(list,nslices);
    if (i==0) VImageDimensions(src[i],nslices,&nt,&nrows,&ncols);
    else {
      VImageDimensions(src[i],xslices,&xnt,&xrows,&xcols);
      if (nrows != xrows) VError(" inconsistent dimensions in input files");
      if (ncols != xcols) VError(" inconsistent dimensions in input files");
      if (nt != xnt) VError(" inconsistent dimensions in input files");
    }

    /* use geometry info from 1st file */
    if (geolist == NULL) geolist = VGetGeoInfo(list);
  }

  /*
  **  CCM
  */
  VImage dest = VCCM(src,mask,nsubjects,nslices,(int)first,(int)length,(int)type,(int)otype);

  /* invert to get discordant map */
  VImage disc = VCreateImageLike(dest);
  VFillImage(disc,VAllBands,0);
  for (b=0; b<VImageNBands(dest); b++) {
    for (r=0; r<VImageNRows(dest); r++) {
      for (c=0; c<VImageNColumns(dest); c++) {
        if (VGetPixel(mask,b,r,c) < 1) continue;
        u = VPixel(dest,b,r,c,VFloat);
        VPixel(disc,b,r,c,VFloat) = 1-u;
      }
    }
  }

  /*
  ** output
  */
  VAttrList out_list = VCreateAttrList();
  if (geolist != NULL) {
    double *D = VGetGeoDim(geolist,NULL);
    D[0] = 3;  /* 3D */
    D[4] = 1;  /* just one timestep */
    VSetGeoDim(geolist,D);
    VSetGeoInfo(geolist,out_list);
  }
  VAppendAttr (out_list,"concordant",NULL,VImageRepn,dest);
  VAppendAttr (out_list,"discordant",NULL,VImageRepn,disc);

  VHistory(VNumber(options),options,prg,&list,&out_list);
  FILE *fp = VOpenOutputFile (out_filename, TRUE);
  if (! VWriteFile (fp, out_list)) exit (1);
  fclose(fp);
  fprintf (stderr, "%s: done.\n", argv[0]);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


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

/***************************************************************************//**
 *
 **/
void QUARK_CORE_slaset2(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int M, int N,
                       float alpha, float *A, int LDA)
{
    DAG_CORE_LASET;
    QUARK_Insert_Task(quark, CORE_slaset2_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,  VALUE,
        sizeof(int),                        &M,     VALUE,
        sizeof(int),                        &N,     VALUE,
        sizeof(float),         &alpha, VALUE,
        sizeof(float)*M*N,     A,      OUTPUT,
        sizeof(int),                        &LDA,   VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_slaset2_quark = PCORE_slaset2_quark
#define CORE_slaset2_quark PCORE_slaset2_quark
#endif
void CORE_slaset2_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int M;
    int N;
    float alpha;
    float *A;
    int LDA;

    quark_unpack_args_6(quark, uplo, M, N, alpha, A, LDA);
    CORE_slaset2(uplo, M, N, alpha, A, LDA);
}


text:
#ifndef H_GPC_UTILS
#define H_GPC_UTILS

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_poly.h>

/* The GSL is a pain to work with. The library DOES NOT HAVE a determinant() function
  or an inv() function: you have to write your own routines. */

#define gmg gsl_matrix_get
#define gms gsl_matrix_set

void m_trans(const gsl_matrix*A, gsl_matrix*A_t);
void m_mult(const gsl_matrix*A, const gsl_matrix*B, gsl_matrix*AB);
void m_add(const gsl_matrix*A, const gsl_matrix*B, gsl_matrix*ApB);
void m_add_to(const gsl_matrix*A, gsl_matrix*B);
void m_scale(double m, gsl_matrix*A);
double m_dot(const gsl_matrix*A,const gsl_matrix*B);
void m_inv(const gsl_matrix*A, gsl_matrix*invA);
double m_det(const gsl_matrix*A);

/* Returns the real part of the roots in roots */
int poly_real_roots(unsigned int n, const double*a, double *roots);

int poly_greatest_real_root(unsigned int n, const double*a, double *root);

void m_display(const char*str, gsl_matrix*m);

#endif


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

#include <config.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>

#include "schur.h"

/*
 * This module contains some routines related to manipulating the
 * Schur form of a matrix which are needed by the eigenvalue solvers
 *
 * This file contains routines based on original code from LAPACK
 * which is distributed under the modified BSD license. The LAPACK
 * routine used is DLANV2.
 */

static inline void schur_standard_form(gsl_matrix *A, gsl_complex *eval1,
                                       gsl_complex *eval2, double *cs,
                                       double *sn);

/*
gsl_schur_standardize()
  Wrapper function for schur_standard_form - convert a 2-by-2 eigenvalue
block to standard form and then update the Schur form and
Schur vectors.

Inputs: T        - Schur form
        row      - row of T of 2-by-2 block to be updated
        eval1    - where to store eigenvalue 1
        eval2    - where to store eigenvalue 2
        update_t - 1 = update the entire matrix T with the transformation
                   0 = do not update rest of T
        Z        - (optional) if non-null, accumulate transformation
*/

void
gsl_schur_standardize(gsl_matrix *T, size_t row, gsl_complex *eval1,
                      gsl_complex *eval2, int update_t, gsl_matrix *Z)
{
  const size_t N = T->size1;
  gsl_matrix_view m;
  double cs, sn;

  m = gsl_matrix_submatrix(T, row, row, 2, 2);
  schur_standard_form(&m.matrix, eval1, eval2, &cs, &sn);

  if (update_t)
    {
      gsl_vector_view xv, yv, v;

      /*
       * The above call to schur_standard_form transformed a 2-by-2 block
       * of T into upper triangular form via the transformation
       *
       * U = [ CS -SN ]
       *     [ SN  CS ]
       *
       * The original matrix T was
       *
       * T = [ T_{11} | T_{12} | T_{13} ]
       *     [   0*   |   A    | T_{23} ]
       *     [   0    |   0*   | T_{33} ]
       *
       * where 0* indicates all zeros except for possibly
       * one subdiagonal element next to A.
       *
       * After schur_standard_form, T looks like this:
       *
       * T = [ T_{11} | T_{12}  | T_{13} ]
       *     [   0*   | U^t A U | T_{23} ]
       *     [   0    |    0*   | T_{33} ]
       *
       * since only the 2-by-2 block of A was changed. However,
       * in order to be able to back transform T at the end,
       * we need to apply the U transformation to the rest
       * of the matrix T since there is no way to apply a
       * similarity transformation to T and change only the
       * middle 2-by-2 block. In other words, let
       *
       * M = [ I 0 0 ]
       *     [ 0 U 0 ]
       *     [ 0 0 I ]
       *
       * and compute
       *
       * M^t T M = [ T_{11} | T_{12} U |   T_{13}   ]
       *           [ U^t 0* | U^t A U  | U^t T_{23} ]
       *           [   0    |   0* U   |   T_{33}   ]
       *
       * So basically we need to apply the transformation U
       * to the i x 2 matrix T_{12} and the 2 x (n - i + 2)
       * matrix T_{23}, where i is the index of the top of A
       * in T.
       *
       * The BLAS routine drot() is suited for this.
       */

      if (row < (N - 2))
        {
          /* transform the 2 rows of T_{23} */

          v = gsl_matrix_row(T, row);
          xv = gsl_vector_subvector(&v.vector,
                                    row + 2,
                                    N - row - 2);

          v = gsl_matrix_row(T, row + 1);
          yv = gsl_vector_subvector(&v.vector,
                                    row + 2,
                                    N - row - 2);

          gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
        }

      if (row > 0)
        {
          /* transform the 2 columns of T_{12} */

          v = gsl_matrix_column(T, row);
          xv = gsl_vector_subvector(&v.vector,
                                    0,
                                    row);

          v = gsl_matrix_column(T, row + 1);
          yv = gsl_vector_subvector(&v.vector,
                                    0,
                                    row);

          gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
        }
    } /* if (update_t) */

  if (Z)
    {
      gsl_vector_view xv, yv;

      /*
       * Accumulate the transformation in Z. Here, Z -> Z * M
       *
       * So:
       *
       * Z -> [ Z_{11} | Z_{12} U | Z_{13} ]
       *      [ Z_{21} | Z_{22} U | Z_{23} ]
       *      [ Z_{31} | Z_{32} U | Z_{33} ]
       *
       * So we just need to apply drot() to the 2 columns
       * starting at index 'row'
       */

      xv = gsl_matrix_column(Z, row);
      yv = gsl_matrix_column(Z, row + 1);

      gsl_blas_drot(&xv.vector, &yv.vector, cs, sn);
    } /* if (Z) */
} /* gsl_schur_standardize() */

/*******************************************************
 *            INTERNAL ROUTINES                        *
 *******************************************************/

/*
schur_standard_form()
  Compute the Schur factorization of a real 2-by-2 matrix in
standard form:

[ A B ] = [ CS -SN ] [ T11 T12 ] [ CS SN ]
[ C D ]   [ SN  CS ] [ T21 T22 ] [-SN CS ]

where either:
1) T21 = 0 so that T11 and T22 are real eigenvalues of the matrix, or
2) T11 = T22 and T21*T12 < 0, so that T11 +/- sqrt(|T21*T12|) are
   complex conjugate eigenvalues

Inputs: A     - 2-by-2 matrix
        eval1 - where to store eigenvalue 1
        eval2 - where to store eigenvalue 2
        cs    - where to store cosine parameter of rotation matrix
        sn    - where to store sine parameter of rotation matrix

Notes: based on LAPACK routine DLANV2
*/

static inline void
schur_standard_form(gsl_matrix *A, gsl_complex *eval1, gsl_complex *eval2,
                    double *cs, double *sn)
{
  double a, b, c, d; /* input matrix values */
  double tmp;
  double p, z;
  double bcmax, bcmis, scale;
  double tau, sigma;
  double cs1, sn1;
  double aa, bb, cc, dd;
  double sab, sac;

  a = gsl_matrix_get(A, 0, 0);
  b = gsl_matrix_get(A, 0, 1);
  c = gsl_matrix_get(A, 1, 0);
  d = gsl_matrix_get(A, 1, 1);

  if (c == 0.0)
    {
      /*
       * matrix is already upper triangular - set rotation matrix
       * to the identity
       */
      *cs = 1.0;
      *sn = 0.0;
    }
  else if (b == 0.0)
    {
      /* swap rows and columns to make it upper triangular */

      *cs = 0.0;
      *sn = 1.0;

      tmp = d;
      d = a;
      a = tmp;
      b = -c;
      c = 0.0;
    }
  else if (((a - d) == 0.0) && (GSL_SIGN(b) != GSL_SIGN(c)))
    {
      /* the matrix has complex eigenvalues with a == d */
      *cs = 1.0;
      *sn = 0.0;
    }
  else
    {
      tmp = a - d;
      p = 0.5 * tmp;
      bcmax = GSL_MAX(fabs(b), fabs(c));
      bcmis = GSL_MIN(fabs(b), fabs(c)) * GSL_SIGN(b) * GSL_SIGN(c);
      scale = GSL_MAX(fabs(p), bcmax);
      z = (p / scale) * p + (bcmax / scale) * bcmis;

      if (z >= 4.0 * GSL_DBL_EPSILON)
        {
          /* real eigenvalues, compute a and d */

          z = p + GSL_SIGN(p) * fabs(sqrt(scale) * sqrt(z));
          a = d + z;
          d -= (bcmax / z) * bcmis;

          /* compute b and the rotation matrix */

          tau = gsl_hypot(c, z);
          *cs = z / tau;
          *sn = c / tau;
          b -= c;
          c = 0.0;
        }
      else
        {
          /*
           * complex eigenvalues, or real (almost) equal eigenvalues -
           * make diagonal elements equal
           */

          sigma = b + c;
          tau = gsl_hypot(sigma, tmp);
          *cs = sqrt(0.5 * (1.0 + fabs(sigma) / tau));
          *sn = -(p / (tau * (*cs))) * GSL_SIGN(sigma);

          /*
           * Compute [ AA BB ] = [ A B ] [ CS -SN ]
           *         [ CC DD ]   [ C D ] [ SN  CS ]
           */
          aa = a * (*cs) + b * (*sn);
          bb = -a * (*sn) + b * (*cs);
          cc = c * (*cs) + d * (*sn);
          dd = -c * (*sn) + d * (*cs);

          /*
           * Compute [ A B ] = [ CS SN ] [ AA BB ]
           *         [ C D ]   [-SN CS ] [ CC DD ]
           */
          a = aa * (*cs) + cc * (*sn);
          b = bb * (*cs) + dd * (*sn);
          c = -aa * (*sn) + cc * (*cs);
          d = -bb * (*sn) + dd * (*cs);

          tmp = 0.5 * (a + d);
          a = d = tmp;

          if (c != 0.0)
            {
              if (b != 0.0)
                {
                  if (GSL_SIGN(b) == GSL_SIGN(c))
                    {
                      /*
                       * real eigenvalues: reduce to upper triangular
                       * form
                       */
                      sab = sqrt(fabs(b));
                      sac = sqrt(fabs(c));
                      p = GSL_SIGN(c) * fabs(sab * sac);
                      tau = 1.0 / sqrt(fabs(b + c));
                      a = tmp + p;
                      d = tmp - p;
                      b -= c;
                      c = 0.0;

                      cs1 = sab * tau;
                      sn1 = sac * tau;
                      tmp = (*cs) * cs1 - (*sn) * sn1;
                      *sn = (*cs) * sn1 + (*sn) * cs1;
                      *cs = tmp;
                    }
                }
              else
                {
                  b = -c;
                  c = 0.0;
                  tmp = *cs;
                  *cs = -(*sn);
                  *sn = tmp;
                }
            }
        }
    }

  /* set eigenvalues */

  GSL_SET_REAL(eval1, a);
  GSL_SET_REAL(eval2, d);
  if (c == 0.0)
    {
      GSL_SET_IMAG(eval1, 0.0);
      GSL_SET_IMAG(eval2, 0.0);
    }
  else
    {
      tmp = sqrt(fabs(b) * fabs(c));
      GSL_SET_IMAG(eval1, tmp);
      GSL_SET_IMAG(eval2, -tmp);
    }

  /* set new matrix elements */

  gsl_matrix_set(A, 0, 0, a);
  gsl_matrix_set(A, 0, 1, b);
  gsl_matrix_set(A, 1, 0, c);
  gsl_matrix_set(A, 1, 1, d);
} /* schur_standard_form() */


text:
#ifndef ALLVARS_H
#define ALLVARS_H

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "core_simulation.h"

#ifdef HDF5
#include <hdf5.h>
#define MODELNAME        "SAGE"
#endif

#define NDIM 3

#define ABORT(sigterm)                                                  \
do {                                                                \
  printf("Error in file: %s\tfunc: %s\tline: %i\n", __FILE__, __FUNCTION__, __LINE__); \
  myexit(sigterm);                                                \
} while(0)

#define  STEPS 10         /* Number of integration intervals between two snapshots */
#define  MAXGALFAC 1
#define  ALLOCPARAMETER 10.0
#define  MAX_NODE_NAME_LEN 50
#define  ABSOLUTEMAXSNAPS 1000  /* The largest number of snapshots for any simulation */
#define  MAXTAGS          300  /* Max number of parameters */

#define  GRAVITY     6.672e-8
#define  SOLAR_MASS  1.989e33
#define  SOLAR_LUM   3.826e33
#define  RAD_CONST   7.565e-15
#define  AVOGADRO    6.0222e23
#define  BOLTZMANN   1.3806e-16
#define  GAS_CONST   8.31425e7
#define  C           2.9979e10
#define  PLANCK      6.6262e-27
#define  CM_PER_MPC  3.085678e24
#define  PROTONMASS  1.6726e-24
#define  HUBBLE      3.2407789e-18   /* in h/sec */

#define  SEC_PER_MEGAYEAR   3.155e13
#define  SEC_PER_YEAR       3.155e7

#define  MAX_STRING_LEN     1024 /* Max length of a string containing a name */

// This structure contains the properties that are output
struct GALAXY_OUTPUT
{
  int   SnapNum;
  int   Type;

  long long   GalaxyIndex;
  long long   CentralGalaxyIndex;
  int   SAGEHaloIndex;
  int   SAGETreeIndex;
  long long   SimulationHaloIndex;

  int   mergeType;  /* 0=none; 1=minor merger; 2=major merger; 3=disk instability; 4=disrupt to ICS */
  int   mergeIntoID;
  int   mergeIntoSnapNum;
  float dT;

  /* (sub)halo properties */
  float Pos[3];
  float Vel[3];
  float Spin[3];
  int   Len;
  float Mvir;
  float CentralMvir;
  float Rvir;
  float Vvir;
  float Vmax;
  float VelDisp;

  /* baryonic reservoirs */
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float BlackHoleMass;
  float ICS;

  /* metals */
  float MetalsColdGas;
  float MetalsStellarMass;
  float MetalsBulgeMass;
  float MetalsHotGas;
  float MetalsEjectedMass;
  float MetalsICS;

  /* to calculate magnitudes */
  float SfrDisk;
  float SfrBulge;
  float SfrDiskZ;
  float SfrBulgeZ;

  /* misc */
  float DiskScaleRadius;
  float Cooling;
  float Heating;
  float QuasarModeBHaccretionMass;
  float TimeOfLastMajorMerger;
  float TimeOfLastMinorMerger;
  float OutflowRate;

  /* infall properties */
  float infallMvir;
  float infallVvir;
  float infallVmax;
};

/* This structure contains the properties used within the code */
struct GALAXY
{
  int   SnapNum;
  int   Type;

  int   GalaxyNr;
  int   CentralGal;
  int   HaloNr;
  long long MostBoundID;

  int   mergeType;  /* 0=none; 1=minor merger; 2=major merger; 3=disk instability; 4=disrupt to ICS */
  int   mergeIntoID;
  int   mergeIntoSnapNum;
  float dT;

  /* (sub)halo properties */
  float Pos[3];
  float Vel[3];
  int   Len;
  float Mvir;
  float deltaMvir;
  float CentralMvir;
  float Rvir;
  float Vvir;
  float Vmax;

  /* baryonic reservoirs */
  float ColdGas;
  float StellarMass;
  float BulgeMass;
  float HotGas;
  float EjectedMass;
  float BlackHoleMass;
  float ICS;

  /* metals */
  float MetalsColdGas;
  float MetalsStellarMass;
  float MetalsBulgeMass;
  float MetalsHotGas;
  float MetalsEjectedMass;
  float MetalsICS;

  /* to calculate magnitudes */
  float SfrDisk[STEPS];
  float SfrBulge[STEPS];
  float SfrDiskColdGas[STEPS];
  float SfrDiskColdGasMetals[STEPS];
  float SfrBulgeColdGas[STEPS];
  float SfrBulgeColdGasMetals[STEPS];

  /* misc */
  float DiskScaleRadius;
  float MergTime;
  double Cooling;
  double Heating;
  float r_heat;
  float QuasarModeBHaccretionMass;
  float TimeOfLastMajorMerger;
  float TimeOfLastMinorMerger;
  float OutflowRate;
  float TotalSatelliteBaryons;

  /* infall properties */
  float infallMvir;
  float infallVvir;
  float infallVmax;
}
*Gal, *HaloGal;

/* auxiliary halo data */
struct halo_aux_data
{
  int DoneFlag;
  int HaloFlag;
  int NGalaxies;
  int FirstGalaxy;
}
*HaloAux;

extern int    FirstFile;    /* first and last file for processing */
extern int    LastFile;

extern int    Ntrees;      /* number of trees in current file  */
extern int    NumGals;     /* Total number of galaxies stored for current tree */
extern int    MaxGals;     /* Maximum number of galaxies allowed for current tree */
extern int    FoF_MaxGals;

extern int    GalaxyCounter;     /* unique galaxy ID for main progenitor line in tree */

extern int    LastSnapShotNr;

extern char   OutputDir[MAX_STRING_LEN];
extern char   FileNameGalaxies[MAX_STRING_LEN];
extern char   TreeName[MAX_STRING_LEN];
extern char   TreeExtension[MAX_STRING_LEN]; // If the trees are in HDF5, they will have a .hdf5 extension. Otherwise they have no extension.
extern char   SimulationDir[MAX_STRING_LEN];
extern char   FileWithSnapList[MAX_STRING_LEN];

extern int    TotHalos;
extern int    TotGalaxies[ABSOLUTEMAXSNAPS];
extern int    *TreeNgals[ABSOLUTEMAXSNAPS];

extern int    *FirstHaloInSnap;

extern int    *TreeNHalos;
extern int    *TreeFirstHalo;

#ifdef MPI
extern int ThisTask, NTask, nodeNameLen;
extern char *ThisNode;
#endif

extern double Omega;
extern double OmegaLambda;
extern double PartMass;
extern double Hubble_h;
extern double EnergySNcode, EnergySN;
extern double EtaSNcode, EtaSN;

/* recipe flags */
extern int    ReionizationOn;
extern int    SupernovaRecipeOn;
extern int    DiskInstabilityOn;
extern int    AGNrecipeOn;
extern int    SFprescription;

/* recipe parameters */
extern int    NParam;
extern char   ParamTag[MAXTAGS][50];
extern int    ParamID[MAXTAGS];
extern void   *ParamAddr[MAXTAGS];

extern double RecycleFraction;
extern double Yield;
extern double FracZleaveDisk;
extern double ReIncorporationFactor;
extern double ThreshMajorMerger;
extern double BaryonFrac;
extern double SfrEfficiency;
extern double FeedbackReheatingEpsilon;
extern double FeedbackEjectionEfficiency;
extern double RadioModeEfficiency;
extern double QuasarModeEfficiency;
extern double BlackHoleGrowthRate;
extern double Reionization_z0;
extern double Reionization_zr;
extern double ThresholdSatDisruption;

extern double UnitLength_in_cm,
  UnitTime_in_s,
  UnitVelocity_in_cm_per_s,
  UnitMass_in_g,
  RhoCrit,
  UnitPressure_in_cgs,
  UnitDensity_in_cgs,
  UnitCoolingRate_in_cgs,
  UnitEnergy_in_cgs,
  UnitTime_in_Megayears,
  G,
  Hubble,
  a0, ar;

extern int    ListOutputSnaps[ABSOLUTEMAXSNAPS];

extern double ZZ[ABSOLUTEMAXSNAPS];
extern double AA[ABSOLUTEMAXSNAPS];
extern double *Age;

extern int    MAXSNAPS;
extern int    NOUT;
extern int    Snaplistlen;

extern gsl_rng *random_generator;

extern int TreeID;
extern int FileNum;

#ifdef HDF5
extern char          *core_output_file;
extern size_t         HDF5_dst_size;
extern size_t        *HDF5_dst_offsets;
extern size_t        *HDF5_dst_sizes;
extern const char   **HDF5_field_names;
extern hid_t         *HDF5_field_types;
extern int            HDF5_n_props;
#endif

#define DOUBLE 1
#define STRING 2
#define INT 3

enum Valid_TreeTypes
{
  genesis_lhalo_hdf5 = 0,
  lhalo_binary = 1,
  num_tree_types
};
enum Valid_TreeTypes TreeType;

#endif  /* #ifndef ALLVARS_H */

max_stars_repo_name: darrencroton/sage

max_issues_repo_name: darrencroton/SAGE

max_forks_repo_name: darrencroton/SAGE


text:
#include <bindings.cmacros.h>
#include <gsl/gsl_vector.h>

BC_INLINE2(GSL_COMPLEX_AT,gsl_vector_complex*,size_t,gsl_complex*)
BC_INLINE2(GSL_COMPLEX_FLOAT_AT,gsl_vector_complex_float*,size_t,gsl_complex_float*)
/* BC_INLINE2(GSL_COMPLEX_LONG_DOUBLE_AT,gsl_vector_complex_long_double*,size_t,gsl_complex_long_double*) */


text:
#ifndef CTETRA_PARTIAL_H
#define CTETRA_PARTIAL_H

#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include "input.h"
#include "evcache.h"
#include "submesh.h"

double Gauss_PartialDos_Eig(double E, double sigma, int orig_index, int eig_index, EvecCache *evCache);

double Gauss_PartialDos(double E, double sigma, int orig_index, EvecCache *evCache);

double** Gauss_PartialDosList(UEInputFn UEfn, int na, int nb, int nc, double sigma, int num_bands, gsl_matrix *R, double **Es, int num_dos);

#endif //CTETRA_PARTIAL_H

max_stars_repo_name: tflovorn/ctetra

max_issues_repo_name: tflovorn/ctetra

max_forks_repo_name: tflovorn/ctetra


text:
/*! \file Integrate.h
 *  \brief header file for integrating functions or data.
 */

#ifndef INTEGRATE_H
#define INTEGRATE_H

#include <cmath>
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>
#include <Precision.h>
#include <Function.h>
#include <Interpolate.h>
#include <Random.h>

using namespace std;
namespace Math
{
    /// Integrate function f(x) from x[lower] to x[upper].  Uses some
    /// algorithm I found somewhere, sometime.
    Double_t IntegrateData(const Double_t x[], const Double_t f[], const int lower,
                         const int upper);

    /// Integrate function f(x) from x[lower] to x[upper].  Uses trapeziod
    /// algorithm.
    Double_t IntegrateTrap(const Double_t x[], const Double_t f[], const int a,
                         const int b);
    ///uses adaptive trapezoidal rule based on numerical recipes
    Double_t IntegrateTrapezoidal(const math_function * f, const Double_t a,
                         const Double_t b, const int n);
    Double_t IntegrateQTrap(const math_function * f, const Double_t a,
                         const Double_t b, const Double_t epsrel, int IMAX=16);
    ///uses simple trapezoidal
    Double_t IntegrateSimpleTrapezoidal(const math_function * f, const Double_t a,
                         const Double_t b, const int n);

    ///extended closed integration
    Double_t IntegrateClosed(const math_function * f, const Double_t a,
                         const Double_t b, const int n);

    ///simpsons 3/8 rule
    Double_t IntegrateSimpson(const math_function * f, const Double_t a,
                         const Double_t b);

    ///Romberg integration based on numerical recipes
    Double_t IntegrateRomberg(const math_function * f, const Double_t a,
                         const Double_t b, const Double_t epsrel, const int n, const int K);

    ///integration based vegas monte carlo integrator. can also return error and chisq
    ///vegas monte is either gsl or NR based vegas integrator.
    Double_t IntegrateVegasMonte(math_multidim_function * f, Double_t *a, Double_t *b, const int n, Double_t *error=0,Double_t*chisq=0,int interations=6);
    double IntegrateVegasMonte(gsl_monte_function * f, double *a, double *b, const int n, double *error=0,double *chisq=0);

    ///Romberg integration based on numerical recipes but uses vegas monte carlo integrator.
    Double_t IntegrateRombergMonte(math_multidim_function * f, Double_t *a,
                         Double_t *b, const Double_t epsrel, const int sampling, const int n, const int K, int interations=6);

    ///subroutine used in my nr vegas integrator.
    void rebin(Double_t rc, int nd, Double_t r[], Double_t xin[], Double_t xi[]);
    ///numerical recipes vegas integrator. for simplicity use IntegrateVegasMonte as interface
    void vegas(math_multidim_function *fxn, Double_t regn[], int ndim, int init,
    unsigned long ncall, int itmx, int nprn, long int idum, Double_t *tgral, Double_t *sd,
    Double_t *chi2a);
}

#endif

max_stars_repo_name: ICRAR/NBodylib

max_issues_repo_name: ICRAR/NBodylib

max_forks_repo_name: ICRAR/NBodylib


text:
#ifndef __bvd_rng_h_
#define __bvd_rng_h_

#include "Cow.h"
#include <cmath>
#include <gsl/gsl_rng.h>
#include <time.h>
#include "Model_Constants.h"
#include <cmath>
#include <gsl/gsl_rng.h>
#include <time.h>

class Random_Number_Generator
{
 public:
	Random_Number_Generator();
	Random_Number_Generator( unsigned int seed );
	void init( unsigned int seed);
	/// draw a random integer/double between lo (inclusive) and hi (exclusive)
	/// from a uniform distribution.
	int          ran_unif_int(    int hi   , int lo=0);
	double       ran_unif_double( double hi, double lo=0 );

	int ran_triangular_int(double lo,double hi,double mod);
	double ran_triangular_double(double lo,double hi,double mod);

	int ran_poisson_int(int r, double mu);

	// true or false, depending on whether a calf ( age < bvd_const::age_threshold_calf ) just infected with BVD will die or not
	bool         will_TI_calf_die();
	double       time_of_death_infected_calf();
	double       time_of_next_infection(double rate);
	Calf_Status calf_outcome_from_infection      ( double  time_of_pregnancy );
	double       time_of_abortion_due_to_infection( double  time_of_pregnancy );
	double       time_of_rest_after_calving(int calving_number);

	// true or false if a birth is a deadbirth, assuming that the infection status of neither mother nor calf are known.
	// This depends on the age of the mother, as indicated in table 4.1 of the handout ("Anteil nicht lebend ausgetragener Kälber, Färsen/Kühle")
	// The wording of the handout is not clear whether in the figures of the handout also abortions are included.
	// In this function, abortions are *not* included.
	bool         is_this_a_deadbirth( bool first_pregnancy );

	//How long is a calf protected by the maternal antibodies in the colostrum.
	double       duration_of_MA();

	//Duration of an infection with BVD ( time that the infection state 'transiently infected' lasts )
	double       duration_of_infection();
	double         duration_of_pregnancy();
	double staggering_first_inseminations();
	int          number_of_calvings();
	virtual bool         is_calf_female();
	//  double         lifetime_male_calves();
	double       first_insemination_age();
	double       lifetime_PI();
	double       time_of_death_as_calf();
	double       insemination_result( bool first_pregnancy , bool* conception );
	double       life_expectancy_male_cow();

	//It has to be decided whether a beginning pregnancy will lead to an abortion or a birth.
	// Also here, pointer *birth is meant for output. if false, abortion is the case. If true birth is the case.
	// The return value is the time until the event (in case of birth -> duration of pregnancy, in case of abortion -> time until abortion ).
	// afaik, in the statistics we have, only the total abortion rates are given.
	// We don't know yet (as of 10.12.2015) whether there is a correlation between abortion happening and the infection status of the mother or the age of the mother.
	// So for now, treat it as uncorrelated. And use the total abortion rates given in the handout.
	double       conception_result( double age_of_mother, Infection_Status is_of_mother, bool* birth );

	void getNRandomNumbersInRange(int n, int rangeLeft, int rangeRight, int* resultArray);
	//double cowWellTimeOfBirth(double time);
	bool bloodTestRightResult();

	double removeTimeAfterFirstTest();
	double removeTimeAfterSecondTest();
	double retestTime();
	unsigned int getSeed();
	bool 		 cowGetsASecondChance();
	double 		timeOfFirstTest();
	virtual bool 		vaccinationWorks();
 private:
  gsl_rng * generator;
unsigned int seed ;
};

#endif


text:
/*! \file vhartree.h
    \brief Hartreeポテンシャルを求めるクラスの宣言

    Copyright ©  2015 @dc1394 All Rights Reserved.
    This software is released under the BSD 2-Clause License.
*/

#ifndef _VHARTREE_H_
#define _VHARTREE_H_

#pragma once

#include "diffdata.h"
#include "property.h"
#include <memory>           // for std::unique_ptr
#include <gsl/gsl_spline.h> // for gsl_interp_accel, gsl_interp_accel_free, gsl_spline, gsl_spline_free

namespace schrac {
    //! A class.
    /*!
        Hartreeポテンシャルを求めるクラス
    */
    class Vhartree final {
        // #region コンストラクタ・デストラクタ

    public:
        //! A constructor.
        /*!
            唯一のコンストラクタ
            \param r_mesh_ rのメッシュ
        */
        Vhartree(std::vector<double> const & r_mesh);

        //! A private copy constructor.
        /*!
            コピーコンストラクタ
        */
        Vhartree(Vhartree const & rhs);

        //! A destructor.
        /*!
            デフォルトデストラクタ
        */
        ~Vhartree() = default;

        // #region メンバ関数

        //!  A public member function (const).
        /*!
            Hartreeポテンシャルの微分値を返す
            \param r 極座標のr
            \return Hartreeポテンシャルの微分値
        */
        double dvhartree_dr(double r) const;

        //!  A public member function.
        /*!
            Hartreeポテンシャルが境界条件を満たすようにセットする
            \param Z 原子核の電荷
        */
        void set_vhartree_boundary_condition(double Z);

        //!  A public member function.
        /*!
            Hartreeポテンシャルを初期化する
        */
        void vhart_init();

        //!  A public member function (const).
        /*!
            Hartreeポテンシャルの値を返す
            \param r 極座標のr
            \return Hartreeポテンシャルの値
        */
        double vhartree(double r) const;

        // #endregion メンバ関数

        // #region プロパティ

    public:
        //! A property.
        /*!
            Hartreeポテンシャルが格納された可変長配列へのプロパティ
        */
        Property<std::vector<double>> Vhart;

        // #endregion プロパティ

        // #region メンバ変数

    private:
        //! A private member variable.
        /*!
            gsl_interp_accelへのスマートポインタ
        */
        std::unique_ptr<gsl_interp_accel, decltype(&gsl_interp_accel_free)> const acc_;

        //! A private member variable.
        /*!
            rのメッシュが格納された可変長配列
        */
        std::vector<double> const r_mesh_;

        //! A private member variable.
        /*!
            gsl_interp_typeへのスマートポインタ
        */
        std::unique_ptr<gsl_spline, decltype(&gsl_spline_free)> const spline_;

        //! A private member variable.
        /*!
            Hartreeポテンシャルが格納された可変長配列
        */
        std::vector<double> vhart_;

        // #endregion メンバ変数

        // #region 禁止されたコンストラクタ・メンバ関数

    private:
        //! A private constructor (deleted).
        /*!
            デフォルトコンストラクタ（禁止）
        */
        Vhartree() = delete;

        //! A private member function (deleted).
        /*!
            operator=()の宣言（禁止）
            \param コピー元のオブジェクト（未使用）
            \return コピー元のオブジェクト
        */
        Vhartree & operator=(Vhartree const &) = delete;

        // #endregion 禁止されたコンストラクタ・メンバ関数
    };
}

#endif  // _VHARTREE_H_


text:
#pragma once
#include "Hypo/Graphics/Exports.h"
#include "Hypo/System/DataTypes/ObjPtr.h"
#include "Hypo/Window/Context/GraphicsContext.h"
#include <gsl/span>

namespace Hypo
{

	enum class TextureType
	{
		None,
		Rgb,
		Rgba,
		sRgb,
		sRgba,
		Stencil,
		Depth
	};

	HYPO_GRAPHICS_API uInt32 GetTexturePixelSize(TextureType type);
	HYPO_GRAPHICS_API uInt32 GetTextureSize(uInt32 width, uInt32 height, TextureType type);

	class TextureData
	{
	public:
		TextureData(uInt32 width, uInt32 height, TextureType type);

		TextureData(uInt32 width, uInt32 height, TextureType type, std::vector<Byte> pixels);

		TextureData() = default;
		TextureData(const TextureData&) = default;

		uInt32 GetWidth() const { return m_Width; }
		uInt32 GetHeight() const { return m_Height; }
		TextureType GetType()const { return m_Type; }
		uInt32 GetPixelsSize() const { return m_Pixels.size(); }
		const std::vector<Byte>& GetPixels() const { return m_Pixels; }
	private:
		uInt32 m_Width = 0;
		uInt32 m_Height = 0;
		TextureType m_Type = TextureType::None;
		std::vector<Byte> m_Pixels;
	};

	HYPO_GRAPHICS_API TextureData TextureFromFile(std::string path);
}

max_stars_repo_name: TheodorLindberg/Hypo

max_issues_repo_name: TheodorLindberg/Hypo

max_forks_repo_name: TheodorLindberg/Hypo


text:
//Causal FIR filtering of each vector in X along dim.
//FIR impulse response is given in vector B with length Q+1.
//(I use P for IIR filter order, since same as polynomial order.)

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

#ifdef __cplusplus
namespace codee {
extern "C" {
#endif

int fir_s (float *Y, const float *X, const float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim);
int fir_d (double *Y, const double *X, const double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim);
int fir_c (float *Y, const float *X, const float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim);
int fir_z (double *Y, const double *X, const double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim);

int fir_s (float *Y, const float *X, const float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in fir_s: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    //Initialize Y
    float b = *B++;
    for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b * *X; }
    X -= N; Y -= N - 1u;

    if (N==0u || L==1u || Q==0u) {}
    else if (L==N)
    {
        if (L<30000u)
        {
            for (size_t q=1u; q<=Q; ++q, ++B, X-=L-q+1u, Y-=L-q)
            {
                b = *B;
                for (size_t l=q; l<L; ++l, ++X, ++Y) { *Y = fmaf(b,*X,*Y); }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, ++B, ++Y) { cblas_saxpy((int)(L-q),*B,X,1,Y,1); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            if (L<30000u)
            {
                for (size_t v=V; v>0u; --v, B-=Q, X+=Q, ++Y)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B)
                    {
                        b = *B;
                        for (size_t l=q; l<L; ++l, ++X, ++Y) { *Y = fmaf(b,*X,*Y); }
                        if (q<Q) { X -= L-q; Y -= L-q-1u; }
                    }
                }
            }
            else
            {
                for (size_t v=V; v>0u; --v, B-=Q, X+=L, Y+=L-Q)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B, ++Y) { cblas_saxpy((int)(L-q),*B,X,1,Y,1); }
                }
            }
        }
        else
        {
            Y += K - 1u;
            for (size_t g=G; g>0u; --g, X+=BS*(L-1u), Y+=BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, ++X, B-=Q, Y-=K*Q-1u)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B, Y+=K) { cblas_saxpy((int)(L-q),*B,X,(int)K,Y,(int)K); }
                }
            }
        }
    }

    return 0;
}

int fir_d (double *Y, const double *X, const double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in fir_d: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;

    //Initialize Y
    double b = *B++;
    for (size_t n=N; n>0u; --n, ++X, ++Y) { *Y = b * *X; }
    X -= N; Y -= N - 1u;

    if (N==0u || L==1u || Q==0u) {}
    else if (L==N)
    {
        if (L<30000u)
        {
            for (size_t q=1u; q<=Q; ++q, ++B, X-=L-q+1u, Y-=L-q)
            {
                b = *B;
                for (size_t l=q; l<L; ++l, ++X, ++Y) { *Y = fma(b,*X,*Y); }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, ++B, ++Y) { cblas_daxpy((int)(L-q),*B,X,1,Y,1); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            if (L<30000u)
            {
                for (size_t v=V; v>0u; --v, B-=Q, X+=Q, ++Y)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B)
                    {
                        b = *B;
                        for (size_t l=q; l<L; ++l, ++X, ++Y) { *Y = fma(b,*X,*Y); }
                        if (q<Q) { X -= L-q; Y -= L-q-1u; }
                    }
                }
            }
            else
            {
                for (size_t v=V; v>0u; --v, B-=Q, X+=L, Y+=L-Q)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B, ++Y) { cblas_daxpy((int)(L-q),*B,X,1,Y,1); }
                }
            }
        }
        else
        {
            Y += K - 1u;
            for (size_t g=G; g>0u; --g, X+=BS*(L-1u), Y+=BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, ++X, B-=Q, Y-=K*Q-1u)
                {
                    for (size_t q=1u; q<=Q; ++q, ++B, Y+=K) { cblas_daxpy((int)(L-q),*B,X,(int)K,Y,(int)K); }
                }
            }
        }
    }

    return 0;
}

int fir_c (float *Y, const float *X, const float *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in fir_c: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    float xr, xi, br, bi;

    //Initialize Y
    br = *B++; bi = *B++;
    for (size_t n=N; n>0u; --n)
    {
        xr = *X++; xi = *X++;
        *Y++ = br*xr - bi*xi;
        *Y++ = br*xi + bi*xr;
    }
    X -= 2u*N; Y -= 2u*N - 2u;

    if (N==0u || L==1u || Q==0u) {}
    else if (L==N)
    {
        if (L<30000u)
        {
            for (size_t q=1u; q<=Q; ++q, X-=2u*(L-q+1u), Y-=2u*(L-q))
            {
                br = *B++; bi = *B++;
                for (size_t l=q; l<L; ++l)
                {
                    xr = *X++; xi = *X++;
                    *Y++ += br*xr - bi*xi;
                    *Y++ += br*xi + bi*xr;
                }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, B+=2u, Y+=2u) { cblas_caxpy((int)(L-q),B,X,1,Y,1); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            if (L<30000u)
            {
                for (size_t v=V; v>0u; --v, B-=2u*Q, X+=2u*Q, Y+=2u)
                {
                    for (size_t q=1u; q<=Q; ++q)
                    {
                        br = *B++; bi = *B++;
                        for (size_t l=q; l<L; ++l)
                        {
                            xr = *X++; xi = *X++;
                            *Y++ += br*xr - bi*xi;
                            *Y++ += br*xi + bi*xr;
                        }
                        if (q<Q) { X -= 2u*(L-q); Y -= 2u*(L-q-1u); }
                    }
                }
            }
            else
            {
                for (size_t v=V; v>0u; --v, X+=2u*L, B-=2u*Q, Y+=2u*(L-Q-1u))
                {
                    for (size_t q=0u; q<=Q; ++q, B+=2u, Y+=2u) { cblas_caxpy((int)(L-q),B,X,1,Y,1); }
                }
            }
        }
        else
        {
            Y += 2u*K - 2u;
            for (size_t g=G; g>0u; --g, X+=2u*BS*(L-1u), Y+=2u*BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, X+=2u, B-=2u*Q, Y-=2u*K*Q-2u)
                {
                    for (size_t q=1u; q<=Q; ++q, B+=2u, Y+=2u*K) { cblas_caxpy((int)(L-q),B,X,(int)K,Y,(int)K); }
                }
            }
        }
    }

    return 0;
}

int fir_z (double *Y, const double *X, const double *B, const size_t R, const size_t C, const size_t S, const size_t H, const int iscolmajor, const size_t Q, const size_t dim)
{
    if (dim>3u) { fprintf(stderr,"error in fir_z: dim must be in [0 3]\n"); return 1; }

    const size_t N = R*C*S*H;
    const size_t L = (dim==0u) ? R : (dim==1u) ? C : (dim==2u) ? S : H;
    double xr, xi, br, bi;

    //Initialize Y
    br = *B++; bi = *B++;
    for (size_t n=N; n>0u; --n)
    {
        xr = *X++; xi = *X++;
        *Y++ = br*xr - bi*xi;
        *Y++ = br*xi + bi*xr;
    }
    X -= 2u*N; Y -= 2u*N - 2u;

    if (N==0u || L==1u || Q==0u) {}
    else if (L==N)
    {
        if (L<30000u)
        {
            for (size_t q=1u; q<=Q; ++q, X-=2u*(L-q+1u), Y-=2u*(L-q))
            {
                br = *B++; bi = *B++;
                for (size_t l=q; l<L; ++l)
                {
                    xr = *X++; xi = *X++;
                    *Y++ += br*xr - bi*xi;
                    *Y++ += br*xi + bi*xr;
                }
            }
        }
        else
        {
            for (size_t q=1u; q<=Q; ++q, B+=2u, Y+=2u) { cblas_zaxpy((int)(L-q),B,X,1,Y,1); }
        }
    }
    else
    {
        const size_t K = (iscolmajor) ? ((dim==0u) ? 1u : (dim==1u) ? R : (dim==2u) ? R*C : R*C*S) : ((dim==0u) ? C*S*H : (dim==1u) ? S*H : (dim==2u) ? H : 1u);
        const size_t BS = (iscolmajor && dim==0u) ? C*S*H : K;
        const size_t V = N/L, G = V/BS;

        if (K==1u && (G==1u || BS==1u))
        {
            if (L<30000u)
            {
                for (size_t v=V; v>0u; --v, B-=2u*Q, X+=2u*Q, Y+=2u)
                {
                    for (size_t q=1u; q<=Q; ++q)
                    {
                        br = *B++; bi = *B++;
                        for (size_t l=q; l<L; ++l)
                        {
                            xr = *X++; xi = *X++;
                            *Y++ += br*xr - bi*xi;
                            *Y++ += br*xi + bi*xr;
                        }
                        if (q<Q) { X -= 2u*(L-q); Y -= 2u*(L-q-1u); }
                    }
                }
            }
            else
            {
                for (size_t v=V; v>0u; --v, X+=2u*L, B-=2u*Q, Y+=2u*(L-Q-1u))
                {
                    for (size_t q=0u; q<=Q; ++q, B+=2u, Y+=2u) { cblas_zaxpy((int)(L-q),B,X,1,Y,1); }
                }
            }
        }
        else
        {
            Y += 2u*K - 2u;
            for (size_t g=G; g>0u; --g, X+=2u*BS*(L-1u), Y+=2u*BS*(L-1u))
            {
                for (size_t bs=BS; bs>0u; --bs, X+=2u, B-=2u*Q, Y-=2u*K*Q-2u)
                {
                    for (size_t q=1u; q<=Q; ++q, B+=2u, Y+=2u*K) { cblas_zaxpy((int)(L-q),B,X,(int)K,Y,(int)K); }
                }
            }
        }
    }

    return 0;
}

#ifdef __cplusplus
}
}
#endif