text:
#pragma once
#define NOMINMAX
#define WEBRTC_WIN
#define ABSL_USES_STD_STRING_VIEW

#include <unknwn.h>
#include <ppltasks.h>
#include <pplawait.h>
#include <gsl/gsl>
#include <winrt/Windows.Foundation.h>
#include <winrt/Windows.Foundation.Diagnostics.h>
#include <winrt/Windows.Foundation.Collections.h>
#include <winrt/Windows.Data.Json.h>
#include <winrt/Windows.Networking.Sockets.h>
#include <winrt/Windows.System.Threading.h>
#include <winrt/Windows.Storage.Streams.h>
#include <wrl/client.h>

ext: h

max_stars_repo_name: WamWooWam/Unicord


text:
/* specfunc/expint.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author: G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_expint.h"

#include "error.h"

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

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

/*
 Chebyshev expansions: based on SLATEC e1.f, W. Fullerton

 Series for AE11       on the interval -1.00000D-01 to  0.
					with weighted error   1.76E-17
					 log weighted error  16.75
			       significant figures required  15.70
				    decimal places required  17.55

 Series for AE12       on the interval -2.50000D-01 to -1.00000D-01
					with weighted error   5.83E-17
					 log weighted error  16.23
			       significant figures required  15.76
				    decimal places required  16.93

 Series for E11        on the interval -4.00000D+00 to -1.00000D+00
					with weighted error   1.08E-18
					 log weighted error  17.97
			       significant figures required  19.02
				    decimal places required  18.61

 Series for E12        on the interval -1.00000D+00 to  1.00000D+00
					with weighted error   3.15E-18
					 log weighted error  17.50
			approx significant figures required  15.8
				    decimal places required  18.10

 Series for AE13       on the interval  2.50000D-01 to  1.00000D+00
					with weighted error   2.34E-17
					 log weighted error  16.63
			       significant figures required  16.14
				    decimal places required  17.33

 Series for AE14       on the interval  0.	    to  2.50000D-01
					with weighted error   5.41E-17
					 log weighted error  16.27
			       significant figures required  15.38
				    decimal places required  16.97
*/

static double AE11_data[39] = {
   0.121503239716065790,
  -0.065088778513550150,
   0.004897651357459670,
  -0.000649237843027216,
   0.000093840434587471,
   0.000000420236380882,
  -0.000008113374735904,
   0.000002804247688663,
   0.000000056487164441,
  -0.000000344809174450,
   0.000000058209273578,
   0.000000038711426349,
  -0.000000012453235014,
  -0.000000005118504888,
   0.000000002148771527,
   0.000000000868459898,
  -0.000000000343650105,
  -0.000000000179796603,
   0.000000000047442060,
   0.000000000040423282,
  -0.000000000003543928,
  -0.000000000008853444,
  -0.000000000000960151,
   0.000000000001692921,
   0.000000000000607990,
  -0.000000000000224338,
  -0.000000000000200327,
  -0.000000000000006246,
   0.000000000000045571,
   0.000000000000016383,
  -0.000000000000005561,
  -0.000000000000006074,
  -0.000000000000000862,
   0.000000000000001223,
   0.000000000000000716,
  -0.000000000000000024,
  -0.000000000000000201,
  -0.000000000000000082,
   0.000000000000000017
};
static cheb_series AE11_cs = {
  AE11_data,
  38,
  -1, 1,
  20
};

static double AE12_data[25] = {
   0.582417495134726740,
  -0.158348850905782750,
  -0.006764275590323141,
   0.005125843950185725,
   0.000435232492169391,
  -0.000143613366305483,
  -0.000041801320556301,
  -0.000002713395758640,
   0.000001151381913647,
   0.000000420650022012,
   0.000000066581901391,
   0.000000000662143777,
  -0.000000002844104870,
  -0.000000000940724197,
  -0.000000000177476602,
  -0.000000000015830222,
   0.000000000002905732,
   0.000000000001769356,
   0.000000000000492735,
   0.000000000000093709,
   0.000000000000010707,
  -0.000000000000000537,
  -0.000000000000000716,
  -0.000000000000000244,
  -0.000000000000000058
};
static cheb_series AE12_cs = {
  AE12_data,
  24,
  -1, 1,
  15
};

static double E11_data[19] = {
  -16.11346165557149402600,
    7.79407277874268027690,
   -1.95540581886314195070,
    0.37337293866277945612,
   -0.05692503191092901938,
    0.00721107776966009185,
   -0.00078104901449841593,
    0.00007388093356262168,
   -0.00000620286187580820,
    0.00000046816002303176,
   -0.00000003209288853329,
    0.00000000201519974874,
   -0.00000000011673686816,
    0.00000000000627627066,
   -0.00000000000031481541,
    0.00000000000001479904,
   -0.00000000000000065457,
    0.00000000000000002733,
   -0.00000000000000000108
};
static cheb_series E11_cs = {
  E11_data,
  18,
  -1, 1,
  13
};

static double E12_data[16] = {
  -0.03739021479220279500,
   0.04272398606220957700,
  -0.13031820798497005440,
   0.01441912402469889073,
  -0.00134617078051068022,
   0.00010731029253063780,
  -0.00000742999951611943,
   0.00000045377325690753,
  -0.00000002476417211390,
   0.00000000122076581374,
  -0.00000000005485141480,
   0.00000000000226362142,
  -0.00000000000008635897,
   0.00000000000000306291,
  -0.00000000000000010148,
   0.00000000000000000315
};
static cheb_series E12_cs = {
  E12_data,
  15,
  -1, 1,
  10
};

static double AE13_data[25] = {
  -0.605773246640603460,
  -0.112535243483660900,
   0.013432266247902779,
  -0.001926845187381145,
   0.000309118337720603,
  -0.000053564132129618,
   0.000009827812880247,
  -0.000001885368984916,
   0.000000374943193568,
  -0.000000076823455870,
   0.000000016143270567,
  -0.000000003466802211,
   0.000000000758754209,
  -0.000000000168864333,
   0.000000000038145706,
  -0.000000000008733026,
   0.000000000002023672,
  -0.000000000000474132,
   0.000000000000112211,
  -0.000000000000026804,
   0.000000000000006457,
  -0.000000000000001568,
   0.000000000000000383,
  -0.000000000000000094,
   0.000000000000000023
};
static cheb_series AE13_cs = {
  AE13_data,
  24,
  -1, 1,
  15
};

static double AE14_data[26] = {
  -0.18929180007530170,
  -0.08648117855259871,
   0.00722410154374659,
  -0.00080975594575573,
   0.00010999134432661,
  -0.00001717332998937,
   0.00000298562751447,
  -0.00000056596491457,
   0.00000011526808397,
  -0.00000002495030440,
   0.00000000569232420,
  -0.00000000135995766,
   0.00000000033846628,
  -0.00000000008737853,
   0.00000000002331588,
  -0.00000000000641148,
   0.00000000000181224,
  -0.00000000000052538,
   0.00000000000015592,
  -0.00000000000004729,
   0.00000000000001463,
  -0.00000000000000461,
   0.00000000000000148,
  -0.00000000000000048,
   0.00000000000000016,
  -0.00000000000000005
};
static cheb_series AE14_cs = {
  AE14_data,
  25,
  -1, 1,
  13
};

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

int gsl_sf_expint_E1_e(const double x, gsl_sf_result * result)
{
  const double xmaxt = -GSL_LOG_DBL_MIN;      /* XMAXT = -LOG (R1MACH(1)) */
  const double xmax  = xmaxt - log(xmaxt);    /* XMAX = XMAXT - LOG(XMAXT) */

  /* CHECK_POINTER(result) */

  if(x < -xmax) {
    OVERFLOW_ERROR(result);
  }
  else if(x <= -10.0) {
    const double s = exp(-x)/x;
    gsl_sf_result result_c;
    cheb_eval_e(&AE11_cs, 20.0/x+1.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * (fabs(x) + 1.0) * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= -4.0) {
    const double s = exp(-x)/x;
    gsl_sf_result result_c;
    cheb_eval_e(&AE12_cs, (40.0/x+7.0)/3.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= -1.0) {
    const double ln_term = -log(fabs(x));
    gsl_sf_result result_c;
    cheb_eval_e(&E11_cs, (2.0*x+5.0)/3.0, &result_c);
    result->val  = ln_term + result_c.val;
    result->err  = result_c.err + GSL_DBL_EPSILON * fabs(ln_term);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x == 0.0) {
    DOMAIN_ERROR(result);
  }
  else if(x <= 1.0) {
    const double ln_term = -log(fabs(x));
    gsl_sf_result result_c;
    cheb_eval_e(&E12_cs, x, &result_c);
    result->val  = ln_term - 0.6875 + x + result_c.val;
    result->err  = result_c.err + GSL_DBL_EPSILON * fabs(ln_term);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= 4.0) {
    const double s = exp(-x)/x;
    gsl_sf_result result_c;
    cheb_eval_e(&AE13_cs, (8.0/x-5.0)/3.0, &result_c);
    result->val  = s * (1.0 + result_c.val);
    result->err  = s * result_c.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x <= xmax) {
    const double s = exp(-x)/x;
    gsl_sf_result result_c;
    cheb_eval_e(&AE14_cs, 8.0/x-1.0, &result_c);
    result->val  = s * (1.0 +  result_c.val);
    result->err  = s * (GSL_DBL_EPSILON + result_c.err);
    result->err += 2.0 * (x + 1.0) * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

int gsl_sf_expint_E2_e(const double x, gsl_sf_result * result)
{
  const double xmaxt = -GSL_LOG_DBL_MIN;
  const double xmax  = xmaxt - log(xmaxt);

  /* CHECK_POINTER(result) */

  if(x < -xmax) {
    OVERFLOW_ERROR(result);
  }
  else if(x < 100.0) {
    const double ex = exp(-x);
    gsl_sf_result result_E1;
    int stat_E1 = gsl_sf_expint_E1_e(x, &result_E1);
    result->val  = ex - x*result_E1.val;
    result->err  = fabs(x) * (GSL_DBL_EPSILON*ex + result_E1.err);
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_E1;
  }
  else if(x < xmax) {
    const double c1  = -2.0;
    const double c2  =  6.0;
    const double c3  = -24.0;
    const double c4  =  120.0;
    const double c5  = -720.0;
    const double c6  =  5040.0;
    const double c7  = -40320.0;
    const double c8  =  362880.0;
    const double c9  = -3628800.0;
    const double c10 =  39916800.0;
    const double c11 = -479001600.0;
    const double c12 =  6227020800.0;
    const double c13 = -87178291200.0;
    const double y = 1.0/x;
    const double sum6 = c6+y*(c7+y*(c8+y*(c9+y*(c10+y*(c11+y*(c12+y*c13))))));
    const double sum  = y*(c1+y*(c2+y*(c3+y*(c4+y*(c5+y*sum6)))));
    result->val = exp(-x) * (1.0 + sum)/x;
    result->err = 2.0 * (x + 1.0) * GSL_DBL_EPSILON * result->val;
    return GSL_SUCCESS;
  }
  else {
    UNDERFLOW_ERROR(result);
  }
}

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

  {
    int status = gsl_sf_expint_E1_e(-x, result);
    result->val = -result->val;
    return status;
  }
}

#if 0
static double recurse_En(int n, double x, double E1)
{
  int i;
  double En = E1;
  double ex = exp(-x);
  for(i=2; i<=n; i++) {
    En = 1./(i-1) * (ex - x * En);
  }
  return En;
}
#endif

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

#include "eval.h"

double gsl_sf_expint_E1(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E1_e(x, &result));
}

double gsl_sf_expint_E2(const double x)
{
  EVAL_RESULT(gsl_sf_expint_E2_e(x, &result));
}

double gsl_sf_expint_Ei(const double x)
{
  EVAL_RESULT(gsl_sf_expint_Ei_e(x, &result));
}

ext: c


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

#include "tests.h"

void
test_gemm (void) {
const double flteps = 1e-4, dbleps = 1e-6;
  {
   int order = 101;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 1.0f;
   float beta = 0.0f;
   float A[] = { 0.199f, 0.237f, 0.456f, 0.377f };
   int lda = 4;
   float B[] = { 0.842f, -0.734f, 0.323f, -0.957f, -0.303f, -0.873f, -0.871f, -0.819f };
   int ldb = 2;
   float C[] = { 0.498f, -0.925f };
   int ldc = 2;
   float C_expected[] = { -0.222426f, -1.07973f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1466)");
     }
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 1.0f;
   float beta = 0.0f;
   float A[] = { -0.83f, 0.922f, -0.228f, -0.003f };
   int lda = 1;
   float B[] = { 0.072f, 0.345f, 0.944f, -0.39f, -0.577f, 0.656f, -0.693f, -0.453f };
   int ldb = 4;
   float C[] = { 0.583f, 0.522f };
   int ldc = 1;
   float C_expected[] = { 0.044268f, 1.24311f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1467)");
     }
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 0.1f;
   float beta = 0.1f;
   float A[] = { -0.838f, 0.622f, -0.494f, 0.304f };
   int lda = 4;
   float B[] = { 0.147f, 0.134f, 0.169f, 0.734f, -0.7f, 0.541f, -0.794f, -0.256f };
   int ldb = 4;
   float C[] = { -0.632f, -0.559f };
   int ldc = 2;
   float C_expected[] = { -0.0532188f, 0.0678514f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1468)");
     }
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 0.1f;
   float beta = 0.1f;
   float A[] = { -0.937f, 0.635f, 0.596f, -0.51f };
   int lda = 1;
   float B[] = { -0.688f, -0.265f, 0.049f, 0.133f, -0.918f, -0.147f, 0.977f, -0.21f };
   int ldb = 2;
   float C[] = { 0.844f, 0.999f };
   int ldc = 1;
   float C_expected[] = { 0.0474373f, 0.135125f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1469)");
     }
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { -0.165f, 0.638f, 0.346f, -0.697f };
   int lda = 1;
   float B[] = { 0.499f, -0.73f, 0.262f, 0.759f, 0.664f, 0.997f, -0.702f, -0.839f };
   int ldb = 2;
   float C[] = { 0.17f, 0.425f };
   int ldc = 2;
   float C_expected[] = { -0.224158f, -0.417831f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1470)");
     }
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = -0.3f;
   float beta = 0.1f;
   float A[] = { -0.603f, -0.714f, -0.893f, 0.046f };
   int lda = 4;
   float B[] = { 0.859f, -0.694f, -0.868f, -0.98f, -0.103f, 0.567f, -0.277f, -0.734f };
   int ldb = 4;
   float C[] = { 0.517f, -0.622f };
   int ldc = 1;
   float C_expected[] = { -0.160575f, -0.0234604f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1471)");
     }
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 0.1f;
   float beta = 1.0f;
   float A[] = { -0.087f, -0.047f, -0.051f, -0.615f };
   int lda = 1;
   float B[] = { -0.722f, -0.077f, 0.563f, 0.501f, 0.855f, 0.605f, 0.556f, -0.627f };
   int ldb = 4;
   float C[] = { -0.181f, -0.89f };
   int ldc = 2;
   float C_expected[] = { -0.208039f, -0.864557f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1472)");
     }
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha = 0.1f;
   float beta = 1.0f;
   float A[] = { -0.753f, -0.074f, -0.247f, -0.19f };
   int lda = 4;
   float B[] = { 0.061f, 0.743f, 0.22f, -0.682f, 0.733f, 0.417f, 0.772f, 0.665f };
   int ldb = 2;
   float C[] = { -0.253f, 0.972f };
   int ldc = 1;
   float C_expected[] = { -0.291994f, 0.898164f };
   cblas_sgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], flteps, "sgemm(case 1473)");
     }
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = 0;
   double beta = 0;
   double A[] = { 0.017, 0.191, 0.863, -0.97 };
   int lda = 4;
   double B[] = { -0.207, -0.916, -0.278, 0.403, 0.885, 0.409, -0.772, -0.27 };
   int ldb = 2;
   double C[] = { -0.274, -0.858 };
   int ldc = 2;
   double C_expected[] = { 0.0, 0.0 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1474)");
     }
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = 0;
   double beta = 0;
   double A[] = { 0.571, 0.081, 0.109, 0.988 };
   int lda = 1;
   double B[] = { -0.048, -0.753, -0.8, -0.89, -0.535, -0.017, -0.018, -0.544 };
   int ldb = 4;
   double C[] = { -0.876, -0.792 };
   int ldc = 1;
   double C_expected[] = { 0.0, 0.0 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1475)");
     }
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { 0.939, 0.705, 0.977, 0.4 };
   int lda = 4;
   double B[] = { -0.089, -0.822, 0.937, 0.159, 0.789, -0.413, -0.172, 0.88 };
   int ldb = 4;
   double C[] = { -0.619, 0.063 };
   int ldc = 2;
   double C_expected[] = { -0.7137904, -0.1270986 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1476)");
     }
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { -0.795, 0.81, 0.388, 0.09 };
   int lda = 1;
   double B[] = { -0.847, 0.031, -0.938, 0.09, -0.286, -0.478, -0.981, 0.881 };
   int ldb = 2;
   double C[] = { -0.242, -0.02 };
   int ldc = 1;
   double C_expected[] = { -0.1562981, -0.0026243 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1477)");
     }
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -1;
   double beta = 0;
   double A[] = { -0.556, 0.532, 0.746, 0.673 };
   int lda = 1;
   double B[] = { -0.525, 0.967, 0.687, -0.024, 0.527, 0.485, 0.109, -0.46 };
   int ldb = 2;
   double C[] = { -0.495, 0.859 };
   int ldc = 2;
   double C_expected[] = { -1.123883, 0.49819 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1478)");
     }
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -1;
   double beta = 0;
   double A[] = { -0.358, 0.224, -0.941, 0.513 };
   int lda = 4;
   double B[] = { -0.201, -0.159, -0.586, -0.016, -0.324, 0.411, 0.115, -0.229 };
   int ldb = 4;
   double C[] = { 0.558, 0.596 };
   int ldc = 1;
   double C_expected[] = { -0.57956, 0.017636 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1479)");
     }
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { -0.586, 0.809, 0.709, -0.524 };
   int lda = 1;
   double B[] = { 0.768, 0.7, 0.619, -0.478, -0.129, -0.778, -0.432, 0.454 };
   int ldb = 4;
   double C[] = { 0.042, 0.252 };
   int ldc = 2;
   double C_expected[] = { -0.1996785, 0.5813976 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1480)");
     }
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha = -0.3;
   double beta = 1;
   double A[] = { -0.164, 0.522, 0.948, -0.624 };
   int lda = 4;
   double B[] = { -0.142, 0.778, 0.359, 0.622, -0.637, -0.757, -0.282, -0.805 };
   int ldb = 2;
   double C[] = { -0.09, 0.183 };
   int ldc = 1;
   double C_expected[] = { -0.0248334, 0.1884672 };
   cblas_dgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[i], C_expected[i], dbleps, "dgemm(case 1481)");
     }
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { -0.082f, -0.281f, -0.096f, 0.913f, 0.974f, -0.706f, -0.773f, 0.522f };
   int lda = 4;
   float B[] = { 0.745f, -0.664f, 0.352f, -0.733f, 0.304f, -0.555f, -0.493f, -0.089f, 0.188f, 0.631f, 0.235f, 0.152f, -0.299f, -0.731f, -0.686f, -0.332f };
   int ldb = 2;
   float C[] = { -0.179f, -0.284f, -0.996f, -0.414f };
   int ldc = 2;
   float C_expected[] = { -1.06679f, 1.47116f, 0.599689f, 0.933532f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1482) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1482) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.0f};
   float A[] = { 0.044f, -0.33f, 0.279f, 0.712f, -0.363f, -0.788f, -0.768f, -0.551f };
   int lda = 1;
   float B[] = { 0.138f, 0.927f, -0.178f, -0.864f, 0.888f, 0.844f, -0.199f, 0.706f, -0.034f, 0.483f, 0.499f, 0.664f, 0.648f, 0.324f, 0.97f, 0.609f };
   int ldb = 4;
   float C[] = { -0.129f, 0.842f, 0.214f, -0.626f };
   int ldc = 1;
   float C_expected[] = { 1.81122f, 1.76205f, 1.0574f, -0.564966f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1483) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1483) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.812f, -0.471f, 0.241f, 0.795f, 0.439f, 0.131f, -0.636f, 0.531f };
   int lda = 4;
   float B[] = { 0.062f, 0.807f, 0.873f, 0.372f, 0.239f, 0.804f, 0.537f, -0.954f, -0.396f, 0.838f, 0.081f, 0.15f, 0.489f, -0.438f, 0.165f, 0.429f };
   int ldb = 4;
   float C[] = { 0.868f, 0.329f, -0.509f, 0.724f };
   int ldc = 2;
   float C_expected[] = { -0.868f, -0.329f, 0.509f, -0.724f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1484) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1484) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.832f, 0.198f, 0.794f, -0.522f, -0.319f, 0.578f, 0.332f, 0.746f };
   int lda = 1;
   float B[] = { -0.361f, 0.187f, -0.163f, -0.781f, 0.536f, 0.888f, -0.969f, 0.899f, 0.961f, -0.583f, 0.753f, 0.29f, -0.997f, 0.729f, -0.352f, -0.2f };
   int ldb = 2;
   float C[] = { 0.864f, 0.735f, -0.074f, -0.228f };
   int ldc = 1;
   float C_expected[] = { -0.864f, -0.735f, 0.074f, 0.228f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1485) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1485) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { 0.149f, 0.187f, 0.263f, -0.715f, -0.882f, -0.907f, 0.87f, -0.527f };
   int lda = 4;
   float B[] = { -0.915f, -0.249f, -0.986f, -0.799f, -0.136f, 0.712f, 0.964f, 0.799f, -0.569f, 0.686f, 0.603f, 0.758f, 0.161f, -0.698f, -0.263f, -0.256f };
   int ldb = 4;
   float C[] = { 0.622f, -0.824f, -0.482f, -0.161f };
   int ldc = 2;
   float C_expected[] = { -0.246901f, 0.083044f, 1.25556f, 0.009106f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1486) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1486) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 1.0f};
   float beta[2] = {0.0f, 0.1f};
   float A[] = { 0.963f, -0.943f, -0.734f, -0.253f, 0.832f, 0.545f, -0.815f, -0.434f };
   int lda = 1;
   float B[] = { 0.23f, -0.211f, 0.906f, 0.232f, -0.339f, 0.597f, -0.919f, 0.793f, 0.535f, 0.526f, 0.119f, 0.053f, 0.751f, 0.044f, 0.752f, -0.469f };
   int ldb = 2;
   float C[] = { 0.483f, -0.266f, -0.224f, -0.692f };
   int ldc = 1;
   float C_expected[] = { -0.047537f, 0.667177f, 1.02025f, 0.823778f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1487) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1487) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { -0.657f, -0.497f, -0.293f, -0.168f, -0.943f, -0.181f, 0.569f, 0.91f };
   int lda = 1;
   float B[] = { -0.047f, 0.796f, -0.913f, 0.998f, 0.365f, 0.467f, -0.627f, -0.523f, 0.885f, 0.234f, -0.494f, 0.071f, -0.361f, -0.154f, -0.055f, -0.32f };
   int ldb = 2;
   float C[] = { 0.956f, 0.268f, 0.152f, 0.717f };
   int ldc = 2;
   float C_expected[] = { -0.668685f, 0.134477f, -0.715786f, -0.478065f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1488) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1488) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {-0.3f, 0.1f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.394f, -0.482f, 0.631f, -0.833f, 0.221f, 0.672f, 0.2f, 0.967f };
   int lda = 4;
   float B[] = { 0.708f, 0.695f, 0.111f, -0.912f, 0.376f, 0.606f, -0.997f, -0.741f, 0.349f, 0.543f, 0.372f, -0.563f, 0.129f, -0.295f, -0.672f, -0.95f };
   int ldb = 4;
   float C[] = { 0.436f, 0.752f, 0.074f, 0.209f };
   int ldc = 1;
   float C_expected[] = { -0.325083f, -0.301952f, -0.283022f, 0.339919f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1489) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1489) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.827f, -0.862f, 0.373f, -0.265f, -0.9f, 0.892f, -0.319f, 0.151f };
   int lda = 1;
   float B[] = { 0.603f, 0.816f, -0.511f, 0.831f, -0.36f, -0.954f, -0.978f, 0.485f, 0.675f, 0.186f, 0.463f, 0.144f, 0.851f, -0.458f, 0.766f, -0.213f };
   int ldb = 4;
   float C[] = { -0.335f, 0.333f, -0.4f, 0.422f };
   int ldc = 2;
   float C_expected[] = { 2.7126f, 0.702111f, 0.437661f, 0.691294f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1490) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1490) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {1.0f, 0.0f};
   float beta[2] = {-0.3f, 0.1f};
   float A[] = { 0.966f, 0.476f, -0.013f, -0.655f, 0.773f, -0.543f, -0.231f, -0.353f };
   int lda = 4;
   float B[] = { -0.684f, 0.144f, 0.018f, -0.77f, -0.688f, 0.909f, -0.094f, -0.938f, -0.757f, 0.574f, -0.479f, 0.473f, 0.0f, 0.064f, -0.168f, 0.858f };
   int ldb = 2;
   float C[] = { -0.912f, 0.54f, 0.756f, 0.024f };
   int ldc = 1;
   float C_expected[] = { -0.156236f, 0.839112f, -0.230206f, -0.106256f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1491) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1491) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { 0.66f, -0.113f, -0.663f, -0.856f, 0.614f, -0.344f, -0.964f, -0.532f };
   int lda = 1;
   float B[] = { -0.606f, -0.965f, -0.279f, -0.312f, 0.63f, 0.967f, 0.041f, -0.557f, 0.663f, 0.619f, -0.134f, 0.261f, -0.388f, 0.525f, 0.222f, 0.538f };
   int ldb = 4;
   float C[] = { 0.114f, -0.376f, -0.851f, -0.682f };
   int ldc = 2;
   float C_expected[] = { 0.114f, -0.376f, -0.851f, -0.682f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1492) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1492) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { 0.212f, -0.752f, 0.679f, 0.49f, -0.029f, -0.488f, 0.567f, 0.374f };
   int lda = 4;
   float B[] = { -0.914f, 0.734f, -0.845f, 0.059f, -0.297f, 0.152f, -0.417f, -0.669f, 0.831f, -0.544f, 0.022f, 0.102f, -0.379f, -0.357f, -0.394f, -0.588f };
   int ldb = 2;
   float C[] = { -0.584f, 0.373f, 0.235f, 0.521f };
   int ldc = 1;
   float C_expected[] = { -0.584f, 0.373f, 0.235f, 0.521f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1493) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1493) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.135f, 0.128f, 0.909f, -0.963f, 0.299f, -0.944f, 0.944f, 0.942f };
   int lda = 1;
   float B[] = { 0.924f, -0.317f, -0.992f, -0.854f, -0.435f, 0.102f, 0.126f, 0.862f, 0.952f, 0.68f, 0.545f, 0.168f, 0.752f, 0.549f, 0.687f, -0.76f };
   int ldb = 2;
   float C[] = { -0.369f, -0.33f, 0.849f, -0.632f };
   int ldc = 2;
   float C_expected[] = { 0.326537f, 0.37603f, -0.86067f, 0.529817f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1494) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1494) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {-1.0f, 0.0f};
   float A[] = { 0.061f, -0.271f, -0.043f, -0.023f, 0.694f, 0.333f, 0.733f, -0.967f };
   int lda = 4;
   float B[] = { 0.088f, -0.607f, 0.589f, 0.375f, -0.897f, -0.954f, -0.216f, -0.195f, -0.865f, -0.511f, -0.219f, 0.535f, 0.976f, 0.582f, 0.464f, -0.041f };
   int ldb = 4;
   float C[] = { 0.533f, -0.63f, 0.405f, 0.667f };
   int ldc = 1;
   float C_expected[] = { -0.459906f, 0.552595f, -0.425391f, -0.533626f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1495) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1495) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { -0.676f, -0.116f, 0.707f, -0.256f, -0.893f, -0.966f, 0.159f, -0.246f };
   int lda = 1;
   float B[] = { 0.059f, 0.281f, -0.93f, -0.263f, 0.583f, -0.11f, 0.639f, -0.96f, -0.878f, 0.984f, 0.058f, 0.977f, -0.567f, 0.561f, -0.048f, -0.798f };
   int ldb = 4;
   float C[] = { 0.362f, -0.808f, 0.428f, -0.112f };
   int ldc = 2;
   float C_expected[] = { 0.362f, -0.808f, 0.428f, -0.112f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1496) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1496) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.0f};
   float beta[2] = {1.0f, 0.0f};
   float A[] = { -0.915f, 0.439f, 0.171f, -0.019f, 0.843f, 0.944f, -0.581f, 0.856f };
   int lda = 4;
   float B[] = { -0.284f, 0.207f, -0.27f, 0.832f, 0.894f, -0.626f, -0.305f, -0.006f, 0.562f, -0.744f, -0.533f, 0.126f, -0.375f, -0.333f, 0.275f, 0.748f };
   int ldb = 2;
   float C[] = { -0.763f, -0.829f, 0.708f, -0.613f };
   int ldc = 1;
   float C_expected[] = { -0.763f, -0.829f, 0.708f, -0.613f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1497) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1497) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { 0.496f, -0.9f, 0.825f, -0.678f, 0.41f, -0.585f, -0.264f, 0.308f };
   int lda = 1;
   float B[] = { 0.907f, 0.972f, -0.724f, 0.745f, -0.601f, 0.589f, 0.759f, -0.521f, -0.161f, -0.321f, 0.341f, -0.981f, -0.378f, -0.671f, -0.314f, -0.878f };
   int ldb = 4;
   float C[] = { -0.293f, 0.07f, 0.087f, -0.542f };
   int ldc = 2;
   float C_expected[] = { 0.10357f, -0.163927f, 0.444626f, -0.0076744f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1498) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1498) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   float alpha[2] = {0.0f, 0.1f};
   float beta[2] = {0.0f, 1.0f};
   float A[] = { -0.225f, -0.629f, -0.939f, -0.836f, -0.841f, -0.794f, 0.836f, -0.65f };
   int lda = 4;
   float B[] = { 0.869f, -0.453f, 0.8f, -0.947f, 0.545f, 0.716f, -0.507f, -0.228f, 0.722f, 0.372f, 0.77f, 0.317f, -0.153f, -0.524f, -0.465f, -0.684f };
   int ldb = 2;
   float C[] = { -0.896f, 0.91f, -0.973f, -0.269f };
   int ldc = 1;
   float C_expected[] = { -1.18974f, -1.0134f, 0.189027f, -1.14494f };
   cblas_cgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], flteps, "cgemm(case 1499) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], flteps, "cgemm(case 1499) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {-1, 0};
   double A[] = { -0.33, 0.457, 0.428, -0.19, 0.86, -0.53, 0.058, -0.942 };
   int lda = 4;
   double B[] = { 0.434, 0.653, -0.124, 0.191, -0.112, -0.84, -0.72, 0.075, -0.503, -0.109, 0.3, -0.898, 0.489, 0.384, 0.993, -0.804 };
   int ldb = 2;
   double C[] = { -0.792, -0.155, -0.608, -0.243 };
   int ldc = 2;
   double C_expected[] = { 0.042563, -0.465908, -0.649991, -1.621116 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1500) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1500) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {-1, 0};
   double A[] = { 0.726, -0.438, -0.23, -0.054, -0.019, 0.902, -0.883, -0.235 };
   int lda = 1;
   double B[] = { 0.159, -0.18, 0.386, -0.167, 0.971, -0.072, 0.87, -0.839, 0.474, 0.956, -0.235, 0.332, 0.826, -0.056, -0.941, 0.01 };
   int ldb = 4;
   double C[] = { -0.799, 0.973, -0.549, -0.177 };
   int ldc = 1;
   double C_expected[] = { -0.181084, 0.257841, 2.251901, 1.558195 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1501) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1501) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.109, 0.892, -0.723, 0.793, 0.109, -0.419, -0.534, 0.448 };
   int lda = 4;
   double B[] = { -0.875, -0.31, -0.027, 0.067, 0.274, -0.126, -0.548, 0.497, 0.681, 0.388, 0.909, 0.889, 0.982, -0.074, -0.788, 0.233 };
   int ldb = 4;
   double C[] = { 0.503, 0.067, 0.239, 0.876 };
   int ldc = 2;
   double C_expected[] = { 0.6553584, 0.0864583, 0.2559136, 0.7518389 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1502) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1502) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {1, 0};
   double A[] = { 0.334, 0.192, -0.992, -0.168, 0.154, -0.75, -0.797, -0.76 };
   int lda = 1;
   double B[] = { -0.82, 0.147, -0.237, 0.68, 0.317, 0.257, -0.406, -0.802, 0.058, 0.012, -0.832, 0.949, -0.263, -0.085, -0.064, 0.492 };
   int ldb = 2;
   double C[] = { 0.079, -0.602, -0.392, 0.316 };
   int ldc = 1;
   double C_expected[] = { 0.0980569, -0.6430449, -0.539207, 0.4226848 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1503) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1503) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 111;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0};
   double beta[2] = {-1, 0};
   double A[] = { -0.305, -0.698, -0.072, -0.383, 0.364, -0.656, 0.819, 0.194 };
   int lda = 4;
   double B[] = { 0.682, 0.498, -0.389, 0.923, -0.853, -0.558, -0.722, -0.085, -0.27, 0.026, -0.107, -0.036, 0.644, -0.327, -0.894, 0.34 };
   int ldb = 4;
   double C[] = { 0.981, -0.336, -0.377, -0.41 };
   int ldc = 2;
   double C_expected[] = { -0.981, 0.336, 0.377, 0.41 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1504) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1504) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 111;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0};
   double beta[2] = {-1, 0};
   double A[] = { -0.306, -0.709, -0.196, 0.285, 0.873, -0.802, 0.715, -0.179 };
   int lda = 1;
   double B[] = { 0.028, 0.109, 0.87, -0.446, 0.735, 0.731, 0.021, -0.186, 0.541, 0.97, -0.333, 0.002, -0.089, -0.01, 0.331, 0.851 };
   int ldb = 2;
   double C[] = { 0.902, -0.584, -0.695, -0.607 };
   int ldc = 1;
   double C_expected[] = { -0.902, 0.584, 0.695, 0.607 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1505) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1505) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.517, -0.136, 0.72, -0.237, 0.121, -0.66, 0.005, 0.759 };
   int lda = 1;
   double B[] = { -0.606, 0.049, 0.807, -0.236, -0.258, -0.412, 0.75, -0.659, 0.993, -0.029, -0.968, 0.707, -0.362, -0.005, 0.096, -0.241 };
   int ldb = 2;
   double C[] = { 0.63, 0.922, 0.025, -0.535 };
   int ldc = 2;
   double C_expected[] = { 1.117044, 1.983417, -1.276831, -0.447092 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1506) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1506) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-1, 0};
   double beta[2] = {1, 0};
   double A[] = { 0.064, 0.371, -0.01, -0.262, 0.143, -0.081, 0.1, -0.062 };
   int lda = 4;
   double B[] = { -0.749, 0.289, -0.239, -0.226, 0.284, 0.668, 0.305, 0.075, -0.36, 0.166, -0.416, 0.234, -0.267, 0.525, 0.116, -0.561 };
   int ldb = 4;
   double C[] = { 0.671, 0.763, 0.444, -0.246 };
   int ldc = 1;
   double C_expected[] = { 0.753107, 0.896395, 0.481996, -0.263126 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1507) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1507) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {-0.3, 0.1};
   double A[] = { -0.956, -0.751, 0.671, -0.633, 0.648, -0.042, 0.948, 0.826 };
   int lda = 1;
   double B[] = { 0.921, 0.506, -0.609, 0.817, -0.686, 0.991, 0.616, -0.482, -0.02, -0.34, 0.559, 0.976, 0.431, 0.385, -0.164, -0.778 };
   int ldb = 4;
   double C[] = { 0.074, -0.01, 0.165, 0.166 };
   int ldc = 2;
   double C_expected[] = { 0.166046, 0.491557, 1.473191, -0.033821 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1508) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1508) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {1, 0};
   double beta[2] = {-0.3, 0.1};
   double A[] = { -0.698, -0.062, 0.023, 0.704, 0.443, -0.46, 0.541, 0.296 };
   int lda = 4;
   double B[] = { 0.787, -0.199, 0.835, -0.276, -0.515, 0.467, -0.76, -0.483, 0.015, -0.394, -0.748, 0.02, 0.573, 0.3, -0.088, -0.238 };
   int ldb = 2;
   double C[] = { 0.935, -0.655, -0.797, 0.071 };
   int ldc = 1;
   double C_expected[] = { -1.070679, 0.178755, -0.344714, -0.308137 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1509) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1509) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 112;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { -0.202, -0.219, 0.741, 0.527, 0.054, 0.16, -0.359, 0.338 };
   int lda = 1;
   double B[] = { -0.872, 0.995, 0.722, 0.618, -0.27, 0.939, -0.743, 0.547, -0.864, 0.376, -0.997, -0.63, 0.887, -0.454, 0.436, -0.039 };
   int ldb = 4;
   double C[] = { -0.684, 0.463, -0.386, -0.524 };
   int ldc = 2;
   double C_expected[] = { 0.1423153, -0.066679, 0.1175618, 0.0012949 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1510) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1510) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 112;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0.1};
   double beta[2] = {-0.3, 0.1};
   double A[] = { -0.855, -0.173, -0.679, 0.824, 0.469, 0.786, 0.757, -0.109 };
   int lda = 4;
   double B[] = { 0.483, -0.888, -0.757, 0.551, -0.81, 0.23, -0.078, 0.725, -0.592, 0.394, 0.884, 0.802, -0.813, -0.016, -0.853, 0.783 };
   int ldb = 2;
   double C[] = { 0.181, -0.368, -0.864, -0.784 };
   int ldc = 1;
   double C_expected[] = { 0.1728438, 0.1183508, 0.2526999, 0.3004174 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1511) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1511) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-1, 0};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.446, -0.65, -0.724, 0.014, 0.792, -0.695, -0.81, -0.358 };
   int lda = 1;
   double B[] = { -0.08, 0.216, 0.689, 0.699, 0.073, -0.346, 0.821, -0.668, -0.798, 0.869, 0.451, -0.061, -0.41, 0.316, 0.104, -0.514 };
   int ldb = 2;
   double C[] = { -0.476, 0.211, -0.912, -0.243 };
   int ldc = 2;
   double C_expected[] = { 1.372475, -0.135616, 0.549353, -1.968747 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1512) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1512) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 111;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-1, 0};
   double beta[2] = {-0.3, 0.1};
   double A[] = { 0.669, 0.046, -0.094, 0.666, 0.23, 0.448, -0.795, -0.142 };
   int lda = 4;
   double B[] = { 0.037, -0.154, -0.739, 0.905, 0.793, -0.53, -0.34, 0.428, 0.072, -0.263, -0.603, -0.905, 0.681, -0.083, -0.511, -0.337 };
   int ldb = 4;
   double C[] = { 0.247, 0.575, -0.836, -0.883 };
   int ldc = 1;
   double C_expected[] = { -0.975939, 0.415528, 0.275533, 0.002716 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1513) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1513) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0};
   double beta[2] = {-1, 0};
   double A[] = { 0.369, 0.506, 0.217, -0.739, -0.395, 0.16, -0.329, -0.954 };
   int lda = 1;
   double B[] = { -0.622, -0.945, 0.416, -0.884, 0.797, -0.74, 0.519, -0.789, -0.348, 0.563, -0.398, -0.956, 0.227, 0.84, -0.079, 0.847 };
   int ldb = 4;
   double C[] = { 0.833, 0.761, 0.074, -0.448 };
   int ldc = 2;
   double C_expected[] = { -0.833, -0.761, -0.074, 0.448 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1514) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1514) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 112;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {0, 0};
   double beta[2] = {-1, 0};
   double A[] = { -0.141, 0.275, 0.717, 0.775, -0.701, -0.689, -0.883, -0.077 };
   int lda = 4;
   double B[] = { -0.526, -0.437, 0.133, -0.209, -0.83, 0.328, 0.916, -0.337, 0.762, -0.664, -0.566, 0.955, 0.168, 0.488, -0.172, -0.535 };
   int ldb = 2;
   double C[] = { -0.88, 0.945, 0.416, 0.99 };
   int ldc = 1;
   double C_expected[] = { 0.88, -0.945, -0.416, -0.99 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1515) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1515) imag");
     };
   };
  };

  {
   int order = 101;
   int transA = 113;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { -0.534, -0.013, -0.258, -0.31, -0.211, -0.883, -0.89, -0.499 };
   int lda = 1;
   double B[] = { -0.185, -0.798, -0.34, 0.716, 0.035, 0.968, -0.26, 0.784, -0.889, -0.344, -0.685, -0.647, -0.764, 0.03, 0.626, -0.989 };
   int ldb = 4;
   double C[] = { -0.793, -0.551, 0.182, 0.838 };
   int ldc = 2;
   double C_expected[] = { -0.5507177, -0.0286821, 0.2222276, 0.5197398 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1516) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1516) imag");
     };
   };
  };

  {
   int order = 102;
   int transA = 113;
   int transB = 113;
   int M = 1;
   int N = 2;
   int K = 4;
   double alpha[2] = {-0.3, 0.1};
   double beta[2] = {0, 0.1};
   double A[] = { 0.575, -0.128, -0.702, 0.758, 0.383, -0.914, 0.157, 0.368 };
   int lda = 4;
   double B[] = { 0.572, -0.841, 0.223, -0.334, -0.823, -0.84, 0.671, -0.871, 0.241, 0.927, -0.344, 0.281, -0.034, -0.104, 0.587, -0.329 };
   int ldb = 2;
   double C[] = { -0.612, 0.167, 0.647, 0.447 };
   int ldc = 1;
   double C_expected[] = { -0.7876717, 0.0341179, -0.0800018, 0.5717566 };
   cblas_zgemm(order, transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc);
   {
     int i;
     for (i = 0; i < 2; i++) {
       gsl_test_rel(C[2*i], C_expected[2*i], dbleps, "zgemm(case 1517) real");
       gsl_test_rel(C[2*i+1], C_expected[2*i+1], dbleps, "zgemm(case 1517) imag");
     };
   };
  };

}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#pragma once

#include <complex>
#include <cstddef>
#include <cstdint>
#include <iosfwd>
#include <string>
#include <utility>
#include <vector>

#include <gsl/gsl>

#include "chainerx/error.h"
#include "chainerx/float16.h"

namespace chainerx {

// NOTE: The dtype list is not fixed yet!!!
enum class Dtype {
    kBool = 1,
    kInt8,
    kInt16,
    kInt32,
    kInt64,
    kUInt8,
    kFloat16,
    kFloat32,
    kFloat64,
};

inline bool IsValidDtype(Dtype dtype) {
    using Underlying = std::underlying_type_t<Dtype>;
    auto value = static_cast<Underlying>(dtype);
    return 1 <= value && value <= static_cast<Underlying>(Dtype::kFloat64);
}

std::ostream& operator<<(std::ostream& os, Dtype dtype);

}  // namespace chainerx

namespace std {

template <>
struct hash<::chainerx::Dtype> {
    size_t operator()(::chainerx::Dtype dtype) const {
        using T = std::underlying_type_t<::chainerx::Dtype>;
        return hash<T>()(static_cast<T>(dtype));
    }
};

}  // namespace std

namespace chainerx {

// Kind of dtypes.
enum class DtypeKind {
    kBool = 0,
    kInt,
    kUInt,
    kFloat,
};

// Gets the single character identifier compatible to NumPy's dtype kind
inline char GetDtypeKindChar(DtypeKind kind) {
    switch (kind) {
        case DtypeKind::kBool:
            return 'b';
        case DtypeKind::kInt:
            return 'i';
        case DtypeKind::kUInt:
            return 'u';
        case DtypeKind::kFloat:
            return 'f';
        default:
            throw DtypeError{"invalid dtype kind"};
    }
}

template <typename T>
constexpr bool IsFloatingPointV = std::is_floating_point<T>::value || std::is_same<std::remove_const_t<T>, Float16>::value;

// Tag type used for dynamic dispatching with dtype value.
//
// This class template is used to resolve mapping from runtime dtype values to compile-time primitive types.
template <typename T>
struct PrimitiveType;

#define CHAINERX_DEFINE_PRIMITIVE_TYPE(name, code, dtype, kind, t, dst) \
    template <>                                                         \
    struct PrimitiveType<t> {                                           \
        using type = t;                                                 \
        using device_storage_type = dst;                                \
        static constexpr char kCharCode = code;                         \
        static constexpr Dtype kDtype = dtype;                          \
        static constexpr int64_t kElementSize = sizeof(type);           \
        static constexpr DtypeKind kKind = kind;                        \
        static const char* GetName() { return name; }                   \
    }

// TODO(niboshi): Char codes are mapped according to current development environment. They should be remapped depending on the executing
// environment, as in NumPy.
CHAINERX_DEFINE_PRIMITIVE_TYPE("bool", '?', Dtype::kBool, DtypeKind::kBool, bool, bool);
CHAINERX_DEFINE_PRIMITIVE_TYPE("int8", 'b', Dtype::kInt8, DtypeKind::kInt, int8_t, int8_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("int16", 'h', Dtype::kInt16, DtypeKind::kInt, int16_t, int16_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("int32", 'i', Dtype::kInt32, DtypeKind::kInt, int32_t, int32_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("int64", 'l', Dtype::kInt64, DtypeKind::kInt, int64_t, int64_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("uint8", 'B', Dtype::kUInt8, DtypeKind::kUInt, uint8_t, uint8_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("float16", 'e', Dtype::kFloat16, DtypeKind::kFloat, chainerx::Float16, uint16_t);
CHAINERX_DEFINE_PRIMITIVE_TYPE("float32", 'f', Dtype::kFloat32, DtypeKind::kFloat, float, float);
CHAINERX_DEFINE_PRIMITIVE_TYPE("float64", 'd', Dtype::kFloat64, DtypeKind::kFloat, double, double);

#undef CHAINERX_DEFINE_PRIMITIVE_TYPE

namespace dtype_detail {
template <typename To, typename From>
using WithConstnessOf = std::conditional_t<std::is_const<From>::value, std::add_const_t<To>, std::remove_const_t<To>>;
}

template <typename T>
using TypeToDeviceStorageType = dtype_detail::WithConstnessOf<typename PrimitiveType<std::remove_const_t<T>>::device_storage_type, T>;

// Dtype mapped from primitive type.
template <typename T>
constexpr Dtype TypeToDtype = PrimitiveType<std::remove_const_t<T>>::kDtype;

// Invokes a function by passing PrimitiveType<T> corresponding to given dtype value.
//
// For example,
//     VisitDtype(Dtype::kInt32, f, args...);
// is equivalent to
//     f(PrimtiveType<int>, args...);
// Note that the dtype argument can be a runtime value. This function can be used for dynamic dispatching based on dtype values.
//
// Note (beam2d): This function should be constexpr, but GCC 5.x does not allow it because of the throw statement, so currently not marked
// as constexpr.
template <typename F, typename... Args>
auto VisitDtype(Dtype dtype, F&& f, Args&&... args) {
    switch (dtype) {
        case Dtype::kBool:
            return std::forward<F>(f)(PrimitiveType<bool>{}, std::forward<Args>(args)...);
        case Dtype::kInt8:
            return std::forward<F>(f)(PrimitiveType<int8_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt16:
            return std::forward<F>(f)(PrimitiveType<int16_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt32:
            return std::forward<F>(f)(PrimitiveType<int32_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt64:
            return std::forward<F>(f)(PrimitiveType<int64_t>{}, std::forward<Args>(args)...);
        case Dtype::kUInt8:
            return std::forward<F>(f)(PrimitiveType<uint8_t>{}, std::forward<Args>(args)...);
        case Dtype::kFloat16:
            return std::forward<F>(f)(PrimitiveType<chainerx::Float16>{}, std::forward<Args>(args)...);
        case Dtype::kFloat32:
            return std::forward<F>(f)(PrimitiveType<float>{}, std::forward<Args>(args)...);
        case Dtype::kFloat64:
            return std::forward<F>(f)(PrimitiveType<double>{}, std::forward<Args>(args)...);
        default:
            throw DtypeError{"invalid dtype: ", static_cast<std::underlying_type_t<Dtype>>(dtype)};
    }
}

// Invokes a function by passing PrimitiveType<T> corresponding to given numeric dtype value.
// See VisitDtype for more detail.
template <typename F, typename... Args>
auto VisitNumericDtype(Dtype dtype, F&& f, Args&&... args) {
    switch (dtype) {
        case Dtype::kInt8:
            return std::forward<F>(f)(PrimitiveType<int8_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt16:
            return std::forward<F>(f)(PrimitiveType<int16_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt32:
            return std::forward<F>(f)(PrimitiveType<int32_t>{}, std::forward<Args>(args)...);
        case Dtype::kInt64:
            return std::forward<F>(f)(PrimitiveType<int64_t>{}, std::forward<Args>(args)...);
        case Dtype::kUInt8:
            return std::forward<F>(f)(PrimitiveType<uint8_t>{}, std::forward<Args>(args)...);
        case Dtype::kFloat16:
            return std::forward<F>(f)(PrimitiveType<chainerx::Float16>{}, std::forward<Args>(args)...);
        case Dtype::kFloat32:
            return std::forward<F>(f)(PrimitiveType<float>{}, std::forward<Args>(args)...);
        case Dtype::kFloat64:
            return std::forward<F>(f)(PrimitiveType<double>{}, std::forward<Args>(args)...);
        default:
            throw DtypeError{"invalid dtype"};
    }
}

// Invokes a function by passing PrimitiveType<T> corresponding to given floating-point dtype value.
// See VisitDtype for more detail.
template <typename F, typename... Args>
auto VisitFloatingPointDtype(Dtype dtype, F&& f, Args&&... args) {
    switch (dtype) {
        case Dtype::kFloat16:
            return std::forward<F>(f)(PrimitiveType<chainerx::Float16>{}, std::forward<Args>(args)...);
        case Dtype::kFloat32:
            return std::forward<F>(f)(PrimitiveType<float>{}, std::forward<Args>(args)...);
        case Dtype::kFloat64:
            return std::forward<F>(f)(PrimitiveType<double>{}, std::forward<Args>(args)...);
        default:
            throw DtypeError{"invalid dtype"};
    }
}

// Gets the single character identifier compatible to NumPy's char code
inline char GetCharCode(Dtype dtype) {
    return VisitDtype(dtype, [](auto pt) { return decltype(pt)::kCharCode; });
}

// Gets the element size of the dtype in bytes.
inline int64_t GetItemSize(Dtype dtype) {
    return VisitDtype(dtype, [](auto pt) { return decltype(pt)::kElementSize; });
}

// Gets the kind of dtype.
inline DtypeKind GetKind(Dtype dtype) {
    return VisitDtype(dtype, [](auto pt) { return decltype(pt)::kKind; });
}

Dtype PromoteTypes(Dtype dt1, Dtype dt2);

// const char* representation of dtype compatible to NumPy's dtype name.
inline const char* GetDtypeName(Dtype dtype) {
    return VisitDtype(dtype, [](auto pt) { return decltype(pt)::GetName(); });
}

// Gets the dtype of given name.
Dtype GetDtype(const std::string& name);

// Returns a vector of all possible dtype values.
std::vector<Dtype> GetAllDtypes();

// Throws an exception if two dtypes mismatch.
void CheckEqual(Dtype lhs, Dtype rhs);

}  // namespace chainerx

max_stars_repo_name: hikjik/chainer

max_issues_repo_name: hitsgub/chainer

max_forks_repo_name: hitsgub/chainer


text:
#ifndef STDEVSTOOL_H
#define STDEVSTOOL_H

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/time.h>

#include "pdevslib.h"
#include <algorithm>    // std::binary_search, std::sort
#include <vector>

//#include "sink/ParameterReader.h"

#define STDEVS_REPRODUCIBLE_SEED 0  // Random reproducible generation. reproducibility of results
#define STDEVS_UNIQUE_SEED -1  // Random unique generation
#define STDEVS_DEFAULT_SEED STDEVS_UNIQUE_SEED  // default in case no configuration

#define SCILAB_VARNAME_REPRODUCIBLE_SIMULATION "ReproducibleSimu" // name to be used in the configuration

//* Forward declare the functions in ParameterReader.h. Including creates circular inclusion problems*/
//#include "sinks/ParameterReader.h"
void printLog(int level, const char *fmt,...);
template<typename T> T readDefaultParameterValue(std::string paramName);

/*
 * Class to generate random numbers.
 * It provides 3 singleton which can be used (it is advised to use always the first one):
 *     1) Random number generation based on the configuration (probably one of the instances below or a custom seeded)
 *     2) Reproducible random numbers (uses always the same seed)
 *     3) need unique random numbers (generates a new seed for every simulation)
 * It is advisable to use the IDistributionParameter instances, so that random number generation is based on the configuration.
 * If you don't want to read from the configuration, it is advisable to use the singleton instances available in this class. Only create new instances if you know what you are doing
 */
class STDEVS {
public:
 // Use this instance in models which need reproducible random numbers (uses always the same seed)
 static STDEVS* getReproducibleGenerationInstance();

 // Use this instance in models which need unique random numbers (generates a new seed)
 static STDEVS* getUniqueGenerationInstance();

 // Unless for specific needs, always use this instance which will point one of the previous based on the configuration or to a custom seeded.
 static STDEVS* getConfigurationInstance();

 /*
  * It is advisable to use the IDistributionParameter instances, so that random number generation is based on the configuration.
  *  If you don't want to read from the configuration, it is advisable to use the singleton instances available in this class. Only create new instances if you know what you are doing.
  *  If you are creating several instances of STDEVS in a very short time (for example when using VectorialDEVS), evaluate using the singletons available in this class to avoid generating the same seed in several instances.
  *
  */
 STDEVS(int randSeed);

  ~STDEVS() {printLog(10, "STDEVS Destroyed \n");};

	unsigned int poisson(double mu) { return gsl_ran_poisson (this->r, mu); }
	double bernoulli (double p) { return gsl_ran_bernoulli (this->r, p); }
	double exponential(double mu) { return gsl_ran_exponential (this->r, mu); }
	double uniform () { return gsl_rng_uniform (this->r); }
	double normal(double mu, double var) { return gsl_ran_gaussian(this -> r, var) + mu;}
	double pareto(double order, double minValue) {	return gsl_ran_pareto(this -> r, order, minValue); }

	double getWeightedElement(std::vector<double> elements, std::vector<double> weights);
	std::vector<double> calculateAccumulatedWeights(std::vector<double> weights);
	std::vector<int> getSortedByWeightSubset(std::vector<double> elements, std::vector<double> weights, const int length);
	std::vector<int> getWeightedSubset(std::vector<double> elements, std::vector<double> weights, const uint length);
	double modifyCumsumAvoidIndex(const int avoidIndex,
			std::vector<double>& cumsum,
			std::vector<double>& weights,
			std::vector<double>& elements);
private:
	 // Use this instance in models which need reproducible random numbers (uses always the same seed)
	 static STDEVS* ReproducibleGenerationInstance;

	 // Use this instance in models which need unique random numbers (generates a new seed)
	 static STDEVS* UniqueGenerationInstance;

	 // Unless for specific needs, always use this instance which will point one of the previous based on the configuration or to a custom seeded.
	 static STDEVS* ConfigurationInstance;
	 static STDEVS* getConfiguredSTDEVS();

	gsl_rng * r; // random seed for the instance

	struct WeightedComparator
	{
	    const std::vector<double> & value_vector;

	    WeightedComparator(const std::vector<double> & val_vec):
	        value_vector(val_vec) {}

	    bool operator()(int i1, int i2)
	    {
	        return value_vector[i1] < value_vector[i2];
	    }
	};
};

#endif

max_stars_repo_name: jorexe/haikunet

max_issues_repo_name: jorexe/haikunet

max_forks_repo_name: jorexe/haikunet


text:
/*
  This file is part of p4est.
  p4est is a C library to manage a collection (a forest) of multiple
  connected adaptive quadtrees or octrees in parallel.

  Copyright (C) 2010 The University of Texas System
  Additional copyright (C) 2011 individual authors
  Written by Carsten Burstedde, Lucas C. Wilcox, and Tobin Isaac

  p4est 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.

  p4est 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 p4est; if not, write to the Free Software Foundation, Inc.,
  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/

#ifndef P4_TO_P8
#include <p4est_plex.h>
#include <p4est_extended.h>
#include <p4est_bits.h>
#else
#include <p8est_plex.h>
#include <p8est_extended.h>
#include <p8est_bits.h>
#endif

#ifdef P4EST_WITH_PETSC
#include <petsc.h>
#include <petscdmplex.h>
#include <petscsf.h>

const char          help[] = "test creating DMPlex from " P4EST_STRING "\n";

static void
locidx_to_PetscInt (sc_array_t * array)
{
  sc_array_t         *newarray;
  size_t              zz, count = array->elem_count;
  P4EST_ASSERT (array->elem_size == sizeof (p4est_locidx_t));

  if (sizeof (p4est_locidx_t) == sizeof (PetscInt)) {
    return;
  }

  newarray = sc_array_new_size (sizeof (PetscInt), array->elem_count);
  for (zz = 0; zz < count; zz++) {
    p4est_locidx_t      il = *((p4est_locidx_t *) sc_array_index (array, zz));
    PetscInt           *ip = (PetscInt *) sc_array_index (newarray, zz);

    *ip = (PetscInt) il;
  }

  sc_array_reset (array);
  sc_array_init_size (array, sizeof (PetscInt), count);
  sc_array_copy (array, newarray);
  sc_array_destroy (newarray);
}

static void
coords_double_to_PetscScalar (sc_array_t * array)
{
  sc_array_t         *newarray;
  size_t              zz, count = array->elem_count;
  P4EST_ASSERT (array->elem_size == 3 * sizeof (double));

  if (sizeof (double) == sizeof (PetscScalar)) {
    return;
  }

  newarray = sc_array_new_size (3 * sizeof (PetscScalar), array->elem_count);
  for (zz = 0; zz < count; zz++) {
    double             *id = (double *) sc_array_index (array, zz);
    PetscScalar        *ip = (PetscScalar *) sc_array_index (newarray, zz);

    ip[0] = (PetscScalar) id[0];
    ip[1] = (PetscScalar) id[1];
    ip[2] = (PetscScalar) id[2];
  }

  sc_array_reset (array);
  sc_array_init_size (array, 3 * sizeof (PetscScalar), count);
  sc_array_copy (array, newarray);
  sc_array_destroy (newarray);
}

static void
locidx_pair_to_PetscSFNode (sc_array_t * array)
{
  sc_array_t         *newarray;
  size_t              zz, count = array->elem_count;
  P4EST_ASSERT (array->elem_size == 2 * sizeof (p4est_locidx_t));

  newarray = sc_array_new_size (sizeof (PetscSFNode), array->elem_count);
  for (zz = 0; zz < count; zz++) {
    p4est_locidx_t     *il = (p4est_locidx_t *) sc_array_index (array, zz);
    PetscSFNode        *ip = (PetscSFNode *) sc_array_index (newarray, zz);

    ip->rank = (PetscInt) il[0];
    ip->index = (PetscInt) il[1];
  }

  sc_array_reset (array);
  sc_array_init_size (array, sizeof (PetscSFNode), count);
  sc_array_copy (array, newarray);
  sc_array_destroy (newarray);
}
#endif

#ifndef P4_TO_P8
static int          refine_level = 5;
#else
static int          refine_level = 3;
#endif

static int
refine_fn (p4est_t * p4est, p4est_topidx_t which_tree,
           p4est_quadrant_t * quadrant)
{
  int                 cid;

  if (which_tree == 2 || which_tree == 3) {
    return 0;
  }

  cid = p4est_quadrant_child_id (quadrant);

  if (cid == P4EST_CHILDREN - 1 ||
      (quadrant->x >= P4EST_LAST_OFFSET (P4EST_MAXLEVEL - 2) &&
       quadrant->y >= P4EST_LAST_OFFSET (P4EST_MAXLEVEL - 2)
#ifdef P4_TO_P8
       && quadrant->z >= P4EST_LAST_OFFSET (P4EST_MAXLEVEL - 2)
#endif
      )) {
    return 1;
  }
  if ((int) quadrant->level >= (refine_level - (int) (which_tree % 3))) {
    return 0;
  }
  if (quadrant->level == 1 && cid == 2) {
    return 1;
  }
  if (quadrant->x == P4EST_QUADRANT_LEN (2) &&
      quadrant->y == P4EST_LAST_OFFSET (2)) {
    return 1;
  }
  if (quadrant->y >= P4EST_QUADRANT_LEN (2)) {
    return 0;
  }

  return 1;
}

static int
refine_tree_one_fn (p4est_t * p4est, p4est_topidx_t which_tree,
                    p4est_quadrant_t * quadrant)
{
  return ! !which_tree;
}

static int
test_forest (int argc, char **argv, p4est_t * p4est, int overlap)
{
  sc_MPI_Comm         mpicomm;
  int                 mpiret;
  int                 mpisize, mpirank;
  sc_array_t         *points_per_dim, *cone_sizes, *cones,
    *cone_orientations, *coords,
    *children, *parents, *childids, *leaves, *remotes;
  p4est_locidx_t      first_local_quad = -1;

  /* initialize MPI */
  mpicomm = p4est->mpicomm;
  mpiret = sc_MPI_Comm_size (mpicomm, &mpisize);
  SC_CHECK_MPI (mpiret);
  mpiret = sc_MPI_Comm_rank (mpicomm, &mpirank);
  SC_CHECK_MPI (mpiret);

  points_per_dim = sc_array_new (sizeof (p4est_locidx_t));
  cone_sizes = sc_array_new (sizeof (p4est_locidx_t));
  cones = sc_array_new (sizeof (p4est_locidx_t));
  cone_orientations = sc_array_new (sizeof (p4est_locidx_t));
  coords = sc_array_new (3 * sizeof (double));
  children = sc_array_new (sizeof (p4est_locidx_t));
  parents = sc_array_new (sizeof (p4est_locidx_t));
  childids = sc_array_new (sizeof (p4est_locidx_t));
  leaves = sc_array_new (sizeof (p4est_locidx_t));
  remotes = sc_array_new (2 * sizeof (p4est_locidx_t));

  p4est_get_plex_data (p4est, P4EST_CONNECT_FULL, (mpisize > 1) ? overlap : 0,
                       &first_local_quad, points_per_dim, cone_sizes, cones,
                       cone_orientations, coords, children, parents, childids,
                       leaves, remotes);

#ifdef P4EST_WITH_PETSC
  {
    PetscErrorCode      ierr;
    DM                  plex, refTree;
    PetscInt            pStart, pEnd;
    PetscSection        parentSection;
    PetscSF             pointSF;
    size_t              zz, count;

    locidx_to_PetscInt (points_per_dim);
    locidx_to_PetscInt (cone_sizes);
    locidx_to_PetscInt (cones);
    locidx_to_PetscInt (cone_orientations);
    coords_double_to_PetscScalar (coords);
    locidx_to_PetscInt (children);
    locidx_to_PetscInt (parents);
    locidx_to_PetscInt (childids);
    locidx_to_PetscInt (leaves);
    locidx_pair_to_PetscSFNode (remotes);

    P4EST_GLOBAL_PRODUCTION ("Begin PETSc routines\n");
    ierr = PetscInitialize (&argc, &argv, 0, help);
    CHKERRQ (ierr);

    ierr = DMPlexCreate (mpicomm, &plex);
    CHKERRQ (ierr);
    ierr = DMSetDimension (plex, P4EST_DIM);
    CHKERRQ (ierr);
    ierr = DMSetCoordinateDim (plex, 3);
    CHKERRQ (ierr);
    ierr = DMPlexCreateFromDAG (plex, P4EST_DIM,
                                (PetscInt *) points_per_dim->array,
                                (PetscInt *) cone_sizes->array,
                                (PetscInt *) cones->array,
                                (PetscInt *) cone_orientations->array,
                                (PetscScalar *) coords->array);
    CHKERRQ (ierr);
    ierr = PetscSFCreate (mpicomm, &pointSF);
    CHKERRQ (ierr);
    ierr =
      DMPlexCreateDefaultReferenceTree (mpicomm, P4EST_DIM, PETSC_FALSE,
                                        &refTree);
    CHKERRQ (ierr);
    ierr = DMPlexSetReferenceTree (plex, refTree);
    CHKERRQ (ierr);
    ierr = DMDestroy (&refTree);
    CHKERRQ (ierr);
    ierr = PetscSectionCreate (mpicomm, &parentSection);
    CHKERRQ (ierr);
    ierr = DMPlexGetChart (plex, &pStart, &pEnd);
    CHKERRQ (ierr);
    ierr = PetscSectionSetChart (parentSection, pStart, pEnd);
    CHKERRQ (ierr);
    count = children->elem_count;
    for (zz = 0; zz < count; zz++) {
      PetscInt            child =
        *((PetscInt *) sc_array_index (children, zz));

      ierr = PetscSectionSetDof (parentSection, child, 1);
      CHKERRQ (ierr);
    }
    ierr = PetscSectionSetUp (parentSection);
    CHKERRQ (ierr);
    ierr =
      DMPlexSetTree (plex, parentSection, (PetscInt *) parents->array,
                     (PetscInt *) childids->array);
    CHKERRQ (ierr);
    ierr = PetscSectionDestroy (&parentSection);
    CHKERRQ (ierr);
    ierr =
      PetscSFSetGraph (pointSF, pEnd - pStart, (PetscInt) leaves->elem_count,
                       (PetscInt *) leaves->array, PETSC_COPY_VALUES,
                       (PetscSFNode *) remotes->array, PETSC_COPY_VALUES);
    CHKERRQ (ierr);
    ierr = DMSetPointSF (plex, pointSF);
    CHKERRQ (ierr);
    ierr = PetscSFDestroy (&pointSF);
    CHKERRQ (ierr);
    ierr = DMViewFromOptions (plex, NULL, "-dm_view");
    CHKERRQ (ierr);
    /* TODO: test with rigid body modes as in plex ex3 */
    ierr = DMDestroy (&plex);
    CHKERRQ (ierr);

    ierr = PetscFinalize ();
    P4EST_GLOBAL_PRODUCTION ("End   PETSc routines\n");
  }
#endif

  sc_array_destroy (points_per_dim);
  sc_array_destroy (cone_sizes);
  sc_array_destroy (cones);
  sc_array_destroy (cone_orientations);
  sc_array_destroy (coords);
  sc_array_destroy (children);
  sc_array_destroy (parents);
  sc_array_destroy (childids);
  sc_array_destroy (leaves);
  sc_array_destroy (remotes);

  return 0;
}

static int
test_big (int argc, char **argv)
{
  sc_MPI_Comm         mpicomm;
  int                 mpiret;
  p4est_t            *p4est;
  p4est_connectivity_t *conn;

  /* initialize MPI */
  mpicomm = sc_MPI_COMM_WORLD;

#ifndef P4_TO_P8
  conn = p4est_connectivity_new_moebius ();
#else
  conn = p8est_connectivity_new_rotcubes ();
#endif
  p4est = p4est_new_ext (mpicomm, conn, 0, 1, 1, 0, NULL, NULL);
  p4est_refine (p4est, 1, refine_fn, NULL);
  p4est_balance (p4est, P4EST_CONNECT_FULL, NULL);
  p4est_partition (p4est, 0, NULL);

  mpiret = test_forest (argc, argv, p4est, 2);
  if (mpiret) {
    return mpiret;
  }

  p4est_destroy (p4est);
  p4est_connectivity_destroy (conn);

  return 0;
}

static int
test_small (int argc, char **argv)
{
  sc_MPI_Comm         mpicomm;
  int                 mpiret;
  p4est_t            *p4est;
  p4est_connectivity_t *conn;

  /* initialize MPI */
  mpicomm = sc_MPI_COMM_WORLD;

#ifndef P4_TO_P8
  conn = p4est_connectivity_new_brick (2, 1, 0, 0);
#else
  conn = p8est_connectivity_new_brick (2, 1, 1, 0, 0, 0);
#endif
  p4est = p4est_new (mpicomm, conn, 0, NULL, NULL);
  p4est_refine (p4est, 0, refine_tree_one_fn, NULL);

  mpiret = test_forest (argc, argv, p4est, 0);
  if (mpiret) {
    return mpiret;
  }

  p4est_partition (p4est, 0, NULL);

  mpiret = test_forest (argc, argv, p4est, 0);
  if (mpiret) {
    return mpiret;
  }

  p4est_destroy (p4est);
  p4est_connectivity_destroy (conn);

  return 0;
}

int
main (int argc, char **argv)
{
  sc_MPI_Comm         mpicomm;
  int                 mpiret;

  /* initialize MPI */
  mpiret = sc_MPI_Init (&argc, &argv);
  SC_CHECK_MPI (mpiret);
  mpicomm = sc_MPI_COMM_WORLD;

  sc_init (mpicomm, 1, 1, NULL, SC_LP_DEFAULT);
  p4est_init (NULL, SC_LP_DEFAULT);

  mpiret = test_small (argc, argv);
  if (mpiret) {
    return mpiret;
  }
  mpiret = test_big (argc, argv);
  if (mpiret) {
    return mpiret;
  }

  sc_finalize ();

  mpiret = sc_MPI_Finalize ();
  SC_CHECK_MPI (mpiret);

  return 0;
}


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

#ifndef __GSL_SYS_H__
#define __GSL_SYS_H__

#include <gsl/gsl_types.h>

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

__BEGIN_DECLS

GSL_EXPORT double gsl_log1p (const double x);
GSL_EXPORT double gsl_expm1 (const double x);
GSL_EXPORT double gsl_hypot (const double x, const double y);
GSL_EXPORT double gsl_acosh (const double x);
GSL_EXPORT double gsl_asinh (const double x);
GSL_EXPORT double gsl_atanh (const double x);

GSL_EXPORT int gsl_isnan (const double x);
GSL_EXPORT int gsl_isinf (const double x);
GSL_EXPORT int gsl_finite (const double x);

GSL_EXPORT double gsl_nan (void);
GSL_EXPORT double gsl_posinf (void);
GSL_EXPORT double gsl_neginf (void);
GSL_EXPORT double gsl_fdiv (const double x, const double y);

GSL_EXPORT double gsl_coerce_double (const double x);
GSL_EXPORT float gsl_coerce_float (const float x);
GSL_EXPORT long double gsl_coerce_long_double (const long double x);

GSL_EXPORT double gsl_ldexp(const double x, const int e);
GSL_EXPORT double gsl_frexp(const double x, int * e);

GSL_EXPORT int gsl_fcmp (const double x1, const double x2, const double epsilon);

__END_DECLS

#endif /* __GSL_SYS_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
#ifndef LOGGER_INCLUDE
#define LOGGER_INCLUDE

#include <gsl/string_span>

#include "log.h"
BOOST_LOG_GLOBAL_LOGGER(
    exec_helper_log_logger,
    execHelper::log::LoggerType); // NOLINT(modernize-use-using)

static const gsl::czstring<> LOG_CHANNEL = "log";
#define LOG(x)                                                                 \
    BOOST_LOG_STREAM_CHANNEL_SEV(exec_helper_log_logger::get(), LOG_CHANNEL,   \
                                 x)                                            \
        << boost::log::add_value(fileLog, __FILE__)                            \
        << boost::log::add_value(lineLog, __LINE__)

#endif /* LOGGER_INCLUDE */

max_stars_repo_name: exec-helper/source

max_issues_repo_name: exec-helper/source

max_forks_repo_name: exec-helper/source


text:
#ifdef __cplusplus
extern "C"
{
#endif // __cplusplus
#ifdef USE_BLAS
#include <cblas.h>

int _use_blas()
{
    return 1;
}
#else // USE_BLAS
#include <xmmintrin.h>

#if defined(_MSC_VER)
#define ALIGNAS(byte_alignment) __declspec(align(byte_alignment))
#elif defined(__GNUC__)
#define ALIGNAS(byte_alignment) __attribute__((aligned(byte_alignment)))
#endif

float cblas_snrm2(const int N, const float *m1, const int incX)
{
    if (N % 4 != 0) {
        fprintf(stderr, "cblas_snrm2() expects N to be a multiple of 4.\n");
        exit(EXIT_FAILURE);
    }

    float norm = 0;
    ALIGNAS(16) float z[4];
    __m128 X;
    __m128 Z = _mm_setzero_ps();

    for (int i=0; i<N; i+=4) {
        X = _mm_load_ps(&m1[i]);
        X = _mm_mul_ps(X, X);
        Z = _mm_add_ps(X, Z);
    }

    _mm_store_ps(z, Z);
    norm += z[0] + z[1] + z[2] + z[3];
    return sqrtf(norm);
}

float cblas_sdot(const int N, const float *m1, const int incX,
                 const float *m2, const int incY)
{
    if (N % 4 != 0) {
        fprintf(stderr, "cblas_sdot() expects N to be a multiple of 4.\n");
        exit(EXIT_FAILURE);
    }

    float dot = 0;
    ALIGNAS(16) float z[4];
    __m128 X, Y;
    __m128 Z = _mm_setzero_ps();

    for (int i=0; i<N; i+=4) {
        X = _mm_load_ps(&m1[i]);
        Y = _mm_load_ps(&m2[i]);
        X = _mm_mul_ps(X, Y);
        Z = _mm_add_ps(X, Z);
    }

    _mm_store_ps(z, Z);
    dot += z[0] + z[1] + z[2] + z[3];
    return dot;
}

int _use_blas()
{
    return 0;
}
#endif // USE_BLAS
#ifdef __cplusplus
}
#endif // __cplusplus


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

#ifndef __GSL_BLOCK_INT_H__
#define __GSL_BLOCK_INT_H__

#include <stdlib.h>
#include <gsl/gsl_errno.h>

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

__BEGIN_DECLS

struct gsl_block_int_struct
{
  size_t size;
  int *data;
};

typedef struct gsl_block_int_struct gsl_block_int;

gsl_block_int *gsl_block_int_alloc (const size_t n);
gsl_block_int *gsl_block_int_calloc (const size_t n);
void gsl_block_int_free (gsl_block_int * b);

int gsl_block_int_fread (FILE * stream, gsl_block_int * b);
int gsl_block_int_fwrite (FILE * stream, const gsl_block_int * b);
int gsl_block_int_fscanf (FILE * stream, gsl_block_int * b);
int gsl_block_int_fprintf (FILE * stream, const gsl_block_int * b, const char *format);

int gsl_block_int_raw_fread (FILE * stream, int * b, const size_t n, const size_t stride);
int gsl_block_int_raw_fwrite (FILE * stream, const int * b, const size_t n, const size_t stride);
int gsl_block_int_raw_fscanf (FILE * stream, int * b, const size_t n, const size_t stride);
int gsl_block_int_raw_fprintf (FILE * stream, const int * b, const size_t n, const size_t stride, const char *format);

size_t gsl_block_int_size (const gsl_block_int * b);
int * gsl_block_int_data (const gsl_block_int * b);

__END_DECLS

#endif /* __GSL_BLOCK_INT_H__ */

max_stars_repo_name: snipekill/FPGen

max_issues_repo_name: snipekill/FPGen

max_forks_repo_name: snipekill/FPGen


text:
/* bccgls.c
   $Revision: 281 $ $Date: 2006-12-13 09:03:53 -0800 (Wed, 13 Dec 2006) $

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

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

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

   BCLS is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
   Public License for more details.

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

*/
/*!
  \file
  This file implements conjugate gradients for the normal equations
  (CGLS).
*/

#include <cblas.h>
#include <string.h>
#include <assert.h>
#include <float.h>

#include "bcls.h"
#include "bclib.h"
#include "bccgls.h"

/*!

  \brief Mat-vec routine called by CGLS.

  CGLS call this routine, which in turn calls bcls_aprod:
  - If      mode == BCLS_PROD_A,   y      = A dxFree.
  - else if mode == BCLS_PROD_At,  dxFree = A'y.

  This routine is declared "static" so that it won't be confused with
  the user's own Aprod routine.

  \param[in,out] ls      BCLS solver context.
  \param[in]     mode    Specifies if the product with A or A' is required.
  \param[in]     m       Number of rows in A.
  \param[in]     nFree   Length of dxFree (and ix which is recovered via ls).
  \param[in,out] dxFree  RHS or LHS (depending on the mode).
  \param[in,out] y       RHS or LHS (depending on the mode).

*/
static void
aprod_free( BCLS *ls, const int mode, const int m, const int nFree,
	    double dxFree[], double y[] )
{
    int    *ix = ls->ix;
    double *dx = ls->dx;  // Used as workspace.

    assert( mode == BCLS_PROD_A || mode == BCLS_PROD_At );

    if (mode == BCLS_PROD_A) {

        bcls_scatter( nFree, ix, dxFree, dx );          // dx(ix) = dxFree.
	bcls_aprod( ls, BCLS_PROD_A, nFree, ix, dx, y );// y = A dx.

    } else {

	bcls_aprod( ls, BCLS_PROD_At, nFree, ix, dx, y );// dx = A' y.
        bcls_gather( nFree, ix, dx, dxFree );            // dxFree = dx(ix);
    }

    return;
}

/*!

  \brief CG for least squares.

  This routine implements CG for solving the symmetric linear system

     ( N'N + damp I ) x = N'r + c

  where N is an m-by-n submatrix of the general matrix A.  This
  implementation is adapted from Hestenes and Stiefel's CG for
  least-squares method to accomodate the damping term and the
  additional vector on the RHS.  (See, eg, Bjork, Numerical Methods
  for Least-Squares Problems, 1996, p. 289.)

  \param[in,out] ls     BCLS solver context.
  \param[in]     m      No. of rows in A.
  \param[in]     n      No. of cols in A.
  \param[in]     kmax   Maximum no. of CG iterations.
  \param[in]     tol    Required optimality tolerance:
                        - norm( A'A x - A'b ) / norm(A'b)
  \param[in]     damp   Regularization term on x.
  \param[in]     c      Additional vector (set NULL if not present).
  \param[out]    x      Solution vector.
  \param[in,out] r      On entry, r must be defined.\n
                        On exit,  r is the residual.
  \param[in,out] s      Work vector.
  \param[in,out] p      Work vector.
  \param[out]    itns   No. of CG iterations.
  \param[out]    opt    Relative residual (see tol, above).
  \param[in]     nout   File for log output (set NULL if not needed).

  \returns
  - 0:  Required accuracy was achieved.
  - 1:  The iteration limit (kmax) was reached.
  - 2:  N'N appears to be singular.

  \todo Update normx rather than recompute at each iteration.

*/
static int
bcls_cgls( BCLS *ls, int m, int n, int kmax, double tol, double damp,
	   double c[], double x[], double r[], double s[], double p[],
	   int *itns, double *opt, FILE *nout )
{
    int
        k = 0, info = 0, converged = 0, unstable = 0;
    double
        xmax = 0.0, resNE, Arnorm, Arnorm0,
        normp, norms, normx = 0.0, alpha, beta,
        gamma, gamma1, delta;
    const int
        damped = damp > 0.0;
    const double
        eps = DBL_EPSILON;

    // Initialize.
    bcls_dload( n, 0.0, x, 1 );                  // x  = 0
    aprod_free( ls, BCLS_PROD_At, m, n, s, r );  // s  = A'r
    if (c)
        cblas_daxpy( n, 1.0, c, 1, s, 1 );       // s += c
    cblas_dcopy( n, s, 1, p, 1 );                // p  = s
    Arnorm0    = cblas_dnrm2( n, s, 1 );
    gamma      = Arnorm0 * Arnorm0;

    if (nout) {
        fprintf(nout,"\n CGLS:\n");
        fprintf(nout,"  Rows............. %8d\t", m);
        fprintf(nout,"  Columns.......... %8d\n", n);
        fprintf(nout,"  Optimality tol... %8.1e\t", tol);
        fprintf(nout,"  Iteration limit.. %8d\n", kmax);
        fprintf(nout,"  Damp............. %8.1e\t", damp);
        fprintf(nout,"  Linear term...... %8s\n",
                (c) ? "yes" : "no" );
        fprintf(nout, "\n\n%5s %17s %17s %9s %10s\n",
                "k", "x(1)","x(n)","normx","resNE");
        fprintf(nout, "%5d %17.10e %17.10e %9.2e %10.2e\n",
                k, x[0], x[n-1], normx, 1.0 );
    }
    // -----------------------------------------------------------------
    // Main loop.
    // -----------------------------------------------------------------
    while ( k < kmax  &&  !converged  &&  !unstable ) {

        k++;

        aprod_free( ls, BCLS_PROD_A, m, n, p, s );   // s = A p

        norms = cblas_dnrm2( m, s, 1 );
        delta = norms * norms;
        if (damped) {
            normp  = cblas_dnrm2( n, p, 1 );
            delta += damp * normp * normp;
        }
        if (delta <= eps)
            delta =  eps;

        alpha = gamma / delta;

        cblas_daxpy( n,  alpha, p, 1, x, 1 );        // x  = x + alpha*p
        cblas_daxpy( m, -alpha, s, 1, r, 1 );        // r  = r - alpha*s
        aprod_free( ls, BCLS_PROD_At, m, n, s, r );  // s  = A'r
        if (damped)
            cblas_daxpy( n, -damp, x, 1, s, 1 );     // s += - damp*x
        if (c)
            cblas_daxpy( n,   1.0, c, 1, s, 1 );     // s += c

        Arnorm = cblas_dnrm2( n, s, 1 );
        gamma1 = gamma;
        gamma  = Arnorm * Arnorm;
        beta   = gamma / gamma1;

        cblas_dscal( n, beta, p, 1 );                // p  = beta*p
        cblas_daxpy( n, 1.0, s, 1, p, 1);            // p += s

        // Check convergence.
        normx     = cblas_dnrm2( n, x, 1 );
        xmax      = fmax( xmax, normx );
        converged = Arnorm <= Arnorm0 * tol;
        unstable  = normx * tol >= 1.0;

        // Output.
        resNE = Arnorm / Arnorm0;
        if (nout)
            fprintf(nout, "%5d %17.10e %17.10e %9.2e %10.2e\n",
                    k, x[0], x[n-1], normx, resNE );

    } // while

    double shrink = normx / xmax;
    if (converged)
        info = 0;
    else if (k >= kmax)
        info = 1;
    else if (shrink <= sqrt(tol))
        info = 2;

    // Output.
    if (nout) fprintf(nout,"\n");
    *itns = k;
    *opt  = resNE;
    return info;
}

/*!

  \brief Compute a Newton step using CGLS.

  See the description of #bcls_newton_step.

  \param[in,out] ls      BCLS solver context
  \param[in]     m       No. of rows in A.
  \param[in]     nFree   No. of free columns, i.e., no. of cols in N.
  \param[in]     ix      Index of vars for which a Newton step is computed.
  \param[in]     damp    Regularization term on dxFree.
  \param[in]     itnLim  Maximum no. of CG iterations.
  \param[in]     tol     Required optimality tolerance:
                         - norm( N'N x - N'b ) / norm(N'b).
  \param[in]     dxFree  Newton direction in the free vars (len = nFree).
  \param[in]     x       Current iterate.
  \param[in]     c       Linear term.  (Set to NULL if not present.)
  \param[in,out] r       On entry, r is current residual: r = b - Ax. \n
                         Will be used as workspace.
  \param[in]     itns    No. of CG iterations.
  \param[in]     opt     Relative residual (see tol, above).

  \return The termination flag returned by #bcls_cgls.

*/
int
bcls_newton_step_cgls( BCLS *ls, int m, int nFree, int ix[], double damp,
                       int itnLim, double tol, double dxFree[], double x[],
                       double c[], double r[], int *itns, double *opt )
{
    int
        j;
    double
        *u;
    const double
        damp2 = damp*damp;

    if (c) {
        u = ls->wrk_u;
        bcls_gather( nFree, ix, c, u );       // u  =   c(ix).
        cblas_dscal( nFree, -1.0, u, 1 );     // u  = - c(ix).
        if ( damp > 0.0 )                     // u += - damp2 x(ix).
            for (j = 0; j < nFree; j++)
                u[j] -= damp2 * x[ix[j]];
    }
    else
        u = NULL;

    return bcls_cgls( ls, m, nFree, itnLim, tol, damp2,
                      u, dxFree, r, ls->wrk_v, ls->wrk_w,
                      itns, opt, ls->minor_file );
}

max_stars_repo_name: echristakopoulou/glslim

max_issues_repo_name: echristakopoulou/glslim

max_forks_repo_name: echristakopoulou/glslim