text:
/* specfunc/test_bessel.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 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.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_sf.h>
#include "test_sf.h"

static double J[100];
static double Y[100];
static double I[100];
static double K[100];

int test_bessel(void)
{
  gsl_sf_result r;
  int i;
  int s = 0;
  int sa;

  TEST_SF(s, gsl_sf_bessel_J0_e, (0.1, &r),     0.99750156206604003230,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J0_e, (2.0, &r),     0.22389077914123566805,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J0_e, (100.0, &r),   0.019985850304223122424,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J0_e, (1.0e+10, &r), 2.1755917502468917269e-06, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_J1_e, (0.1, &r),      0.04993752603624199756,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J1_e, (2.0, &r),      0.57672480775687338720,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J1_e, (100.0, &r),   -0.07714535201411215803,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_J1_e, (1.0e+10, &r), -7.676508175684157103e-06, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Jn_e, (4, 0.1, &r),     2.6028648545684032338e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (5, 2.0, &r),     0.007039629755871685484,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (10, 20.0, &r),   0.18648255802394508321,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (100, 100.0, &r), 0.09636667329586155967,    TEST_TOL0, GSL_SUCCESS);

  /* exercise the BUG#3 problem */
  TEST_SF(s, gsl_sf_bessel_Jn_e, (2, 900.0, &r), -0.019974345269680646400, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (2, 15000.0, &r), -0.0020455820181216382666, TEST_TOL4, GSL_SUCCESS);

#ifdef TEST_LARGE
  TEST_SF(s, gsl_sf_bessel_Jn_e, (0, 1.0e+10, &r), 2.1755917502468917269e-06, TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (1, 1.0e+10, &r), -7.676508175684157103e-06, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jn_e, (0, 20000, &r), 0.00556597490495494615709982972, TEST_TOL4, GSL_SUCCESS);
#endif

  /*  Testcase demonstrating long calculation time:
      Time spent in gsl_sf_bessel_J_CF1 for large x<1000 and n<5
      grows in proportion to x
      Jonny Taylor <j.m.taylor@durham.ac.uk>  */
  TEST_SF(s, gsl_sf_bessel_Jn_e, (45, 900.0, &r),     0.02562434700634278108,    TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Y0_e, (0.1, &r),         -1.5342386513503668441,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y0_e, (2, &r),            0.5103756726497451196,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y0_e, (256.0, &r),       -0.03381290171792454909 ,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y0_e, (4294967296.0, &r), 3.657903190017678681e-06, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Y1_e, (0.1, &r),         -6.45895109470202698800,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y1_e, (2, &r),           -0.10703243154093754689,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y1_e, (100.0, &r),       -0.020372312002759793305,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Y1_e, (4294967296.0, &r), 0.000011612249378370766284, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Yn_e, (4, 0.1, &r),            -305832.29793353160319,    TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Yn_e, (5, 2, &r),              -9.935989128481974981,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Yn_e, (100, 100.0, &r),        -0.16692141141757650654,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Yn_e, (100, 4294967296.0, &r),  3.657889671577715808e-06, TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Yn_e, (1000, 4294967296.0, &r), 3.656551321485397501e-06, 2.0e-05, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Yn_e, (2, 15000.0, &r), -0.006185217273358617849, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_I0_scaled_e, (1e-10, &r),   0.99999999990000000001,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_scaled_e, (0.1, &r),     0.90710092578230109640,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_scaled_e, (2, &r),       0.30850832255367103953,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_scaled_e, (100.0, &r),   0.03994437929909668265,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_scaled_e, (65536.0, &r), 0.0015583712551952223537, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_I1_scaled_e, (0.1, &r),     0.04529844680880932501,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I1_scaled_e, (2, &r),       0.21526928924893765916,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I1_scaled_e, (100.0, &r),   0.03974415302513025267,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I1_scaled_e, (65536.0, &r), 0.0015583593657207350452, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_In_scaled_e, (  -4,    0.1, &r), 2.3575258620054605307e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_scaled_e, (   4,    0.1, &r), 2.3575258620054605307e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_scaled_e, (   5,    2.0, &r), 0.0013297610941881578142, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_scaled_e, ( 100,  100.0, &r), 1.7266862628167695785e-22, TEST_TOL0, GSL_SUCCESS);

  /* BJG: the "exact" values in the following two tests were originally computed from the
     taylor series for I_nu using "long double" and rescaling.  The last few digits
     were inaccurate due to cumulative roundoff.

     BJG: 2006/05 I have now replaced these with the term asymptotic
     expansion from A&S 9.7.1 which should be fully accurate. */

  TEST_SF(s, gsl_sf_bessel_In_scaled_e, (   2,    1e7, &r), 1.261566024466416433e-4, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_scaled_e, (   2,    1e8, &r), 3.989422729212649531e-5, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_I0_e, (0.1, &r), 1.0025015629340956014, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_e, (2.0, &r), 2.2795853023360672674, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I0_e, (100.0, &r), 1.0737517071310738235e+42, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_I1_e, (0.1, &r), 0.05006252604709269211,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I1_e, (2.0, &r), 1.59063685463732906340,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_I1_e, (100.0, &r), 1.0683693903381624812e+42, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_In_e, (   4,    0.1, &r), 2.6054690212996573677e-07, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_e, (   5,    2.0, &r), 0.009825679323131702321,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_In_e, ( 100,  100.0, &r), 4.641534941616199114e+21,  TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_K0_scaled_e, (0.1, &r), 2.6823261022628943831, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K0_scaled_e, (2.0, &r), 0.8415682150707714179, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K0_scaled_e, (100.0, &r), 0.1251756216591265789, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_K1_scaled_e, (0.1, &r), 10.890182683049696574, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K1_scaled_e, (2.0, &r), 1.0334768470686885732, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K1_scaled_e, (100.0, &r), 0.1257999504795785293, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Kn_scaled_e, (   4,    0.1, &r), 530040.2483725626207, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Kn_scaled_e, (   5,    2.0, &r), 69.68655087607675118, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Kn_scaled_e, ( 100,  100.0, &r), 2.0475736731166756813e+19, TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_K0_e, (0.1, &r), 2.4270690247020166125, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K0_e, (2.0, &r), 0.11389387274953343565, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K0_e, (100.0, &r), 4.656628229175902019e-45, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_K1_e, (0.1, &r), 9.853844780870606135,       TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K1_e, (2.0, &r), 0.13986588181652242728,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_K1_e, (100.0, &r), 4.679853735636909287e-45, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Kn_e, (   4,    0.1, &r), 479600.2497925682849,     TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Kn_e, (   5,    2.0, &r), 9.431049100596467443,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Kn_e, ( 100,  100.0, &r), 7.617129630494085416e-25, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_j0_e, (-10.0, &r), -0.05440211108893698134, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j0_e, (0.001, &r), 0.9999998333333416667, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j0_e, (  1.0, &r), 0.84147098480789650670, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j0_e, ( 10.0, &r), -0.05440211108893698134, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j0_e, (100.0, &r), -0.005063656411097587937, TEST_TOL1, GSL_SUCCESS);
#ifdef FIXME
  TEST_SF(s,  gsl_sf_bessel_j0_e, (1048576.0, &r), 3.1518281938718287624e-07, TEST_TOL2, GSL_SUCCESS);
#endif

  /* these values are from Mathematica */
#ifdef FIXME
  TEST_SF(s,  gsl_sf_bessel_j0_e, (1.0e18, &r), -9.9296932074040507620955e-19, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j0_e, (1.0e20, &r), -6.4525128526578084420581e-21, TEST_TOL0, GSL_SUCCESS);
#endif

  TEST_SF(s,  gsl_sf_bessel_j1_e, (-10.0, &r), -0.07846694179875154709, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j1_e, (0.01, &r), 0.003333300000119047399, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j1_e, (  1.0, &r), 0.30116867893975678925, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j1_e, ( 10.0, &r), 0.07846694179875154709, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j1_e, (100.0, &r), -0.008673825286987815220, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j1_e, (1048576.0, &r), -9.000855242905546158e-07, TEST_TOL0, GSL_SUCCESS);

  /*TEST_SF(s,  gsl_sf_bessel_j1_e, (1.0e18, &r), -1.183719902187107336049e-19, TEST_TOL0, GSL_SUCCESS);*/

  TEST_SF(s,  gsl_sf_bessel_j2_e, (-10.0, &r), 0.07794219362856244547, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j2_e, (0.01, &r), 6.666619047751322551e-06, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j2_e, (  1.0, &r), 0.06203505201137386110, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j2_e, ( 10.0, &r), 0.07794219362856244547, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j2_e, (100.0, &r), 0.004803441652487953480, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_j2_e, (1048576.0, &r), -3.1518539455252413111e-07, TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_jl_e, (0, 0.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_jl_e, (1,       10.0, &r),   0.07846694179875154709000, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (5,        1.0, &r),   0.00009256115861125816357, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (10,      10.0, &r),   0.06460515449256426427,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (100,    100.0, &r),   0.010880477011438336539,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (2000, 1048576.0, &r), 7.449384239168568534e-07,  TEST_SQRT_TOL0, GSL_SUCCESS);

  /* related to BUG#3 problem */
  TEST_SF(s, gsl_sf_bessel_jl_e, (2, 900.0, &r),   -0.0011089115568832940086,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_jl_e, (2, 15000.0, &r), -0.00005955592033075750554, TEST_TOL4, GSL_SUCCESS);

  /* Bug report by Mario Santos, value computed from AS 10.1.8 */
  TEST_SF(s,  gsl_sf_bessel_jl_e, (100, 1000.0, &r), -0.00025326311230945818285, TEST_TOL4, GSL_SUCCESS);

  /* Bug reported by Koichi Takahashi <ktakahashi@molsci.org>,
     computed from Pari besseljh(n,x) and AS 10.1.1 */

  TEST_SF(s,  gsl_sf_bessel_jl_e, (30, 3878.62, &r), -0.00023285567034330878410434732790, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (49, 9912.63, &r), 5.2043354544842669214485107019E-5 , TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (49, 9950.35, &r), 5.0077368819565969286578715503E-5 , TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_jl_e, (52, 9930.51, &r), -7.4838588266727718650124475651E-6 , TEST_TOL4, GSL_SUCCESS);

  /* bug report #37209 */
  TEST_SF(s,  gsl_sf_bessel_jl_e, (364, 36.62, &r), 1.118907148986954E-318, TEST_TOL0, GSL_SUCCESS);
  /*TEST_SF(s,  gsl_sf_bessel_jl_e, (149, 1.0, &r), 2.6599182755508469E-307, TEST_TOL0, GSL_SUCCESS);*/

  TEST_SF(s,  gsl_sf_bessel_y0_e, (0.001, &r), -999.99950000004166670, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y0_e, (  1.0, &r), -0.5403023058681397174, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y0_e, ( 10.0, &r), 0.08390715290764524523, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y0_e, (100.0, &r), -0.008623188722876839341, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y0_e, (65536.0, &r), 0.000011014324202158573930, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y0_e, (4294967296.0, &r), 2.0649445131370357007e-10, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_y1_e, ( 0.01, &r), -10000.499987500069444, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y1_e, (  1.0, &r), -1.3817732906760362241, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y1_e, ( 10.0, &r), 0.06279282637970150586, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y1_e, (100.0, &r), 0.004977424523868819543, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y1_e, (4294967296.0, &r), 1.0756463271573404688e-10, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_y2_e, ( 0.01, &r), -3.0000500012499791668e+06, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y2_e, (  1.0, &r), -3.605017566159968955, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y2_e, ( 10.0, &r), -0.06506930499373479347, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y2_e, (100.0, &r), 0.008772511458592903927, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_y2_e, (4294967296.0, &r), -2.0649445123857054207e-10, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s,  gsl_sf_bessel_yl_e, (0,        0.01, &r), -99.995000041666528,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (0,        1.0, &r),  -0.54030230586813972,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (1,       10.0, &r),   0.062792826379701506,   TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (5,        1.0, &r),  -999.44034339223641,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (10,       0.01, &r), -6.5473079797378378e+30, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (10,      10.0, &r),  -0.172453672088057849,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (100,      1.0, &r),  -6.6830794632586775e+186, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (100,    100.0, &r),  -0.0229838504915622811,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s,  gsl_sf_bessel_yl_e, (2000, 1048576.0, &r), 5.9545201447146155e-07,  TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_i0_scaled_e, (0.0, &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i0_scaled_e, (0.1, &r), 0.9063462346100907067, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i0_scaled_e, (2.0, &r), 0.24542109027781645493, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i0_scaled_e, (100.0, &r), 0.005000000000000000000, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_i1_scaled_e, (0.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i1_scaled_e, (0.1, &r), 0.030191419289002226846, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i1_scaled_e, (2.0, &r), 0.131868364583275317610, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i1_scaled_e, (100.0, &r), 0.004950000000000000000, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_i2_scaled_e, (0.0, &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i2_scaled_e, (0.1, &r), 0.0006036559400239012567, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i2_scaled_e, (2.0, &r), 0.0476185434029034785100, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_i2_scaled_e, (100.0, &r), 0.0048515000000000000000, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   0, 0.0,   &r), 1.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   1, 0.0,   &r), 0.0, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   4, 0.001, &r), 1.0571434341190365013e-15, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   4,   0.1, &r), 9.579352242057134927e-08,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   5,   2.0, &r), 0.0004851564602127540059,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, (   5, 100.0, &r), 0.004300446777500000000,   TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_il_scaled_e, ( 100, 100.0, &r), 1.3898161964299132789e-23, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_k0_scaled_e, (0.1, &r), 15.707963267948966192, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k0_scaled_e, (2.0, &r), 0.7853981633974483096, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k0_scaled_e, (100.0, &r), 0.015707963267948966192, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_k1_scaled_e, (0.1, &r), 172.78759594743862812, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k1_scaled_e, (2.0, &r), 1.1780972450961724644, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k1_scaled_e, (100.0, &r), 0.015865042900628455854, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_k2_scaled_e, (0.1, &r), 5199.335841691107810, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k2_scaled_e, (2.0, &r), 2.5525440310417070063, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_k2_scaled_e, (100.0, &r), 0.016183914554967819868, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_kl_scaled_e, (   4, 1.0/256.0, &r), 1.8205599816961954439e+14, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_kl_scaled_e, (   4, 1.0/8.0, &r),   6.1173217814406597530e+06, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_kl_scaled_e, (   5,   2.0, &r),     138.10735829492005119,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_kl_scaled_e, ( 100, 100.0, &r),     3.985930768060258219e+18,  TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Jnu_e, (0.0001, 1.0, &r),         0.7652115411876708497,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (0.0001, 10.0, &r),       -0.2459270166445205,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (0.0009765625, 10.0, &r), -0.2458500798634692,     TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (0.75, 1.0, &r),           0.5586524932048917478,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (0.75, 10.0, &r),         -0.04968928974751508135, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, ( 1.0, 0.001, &r), 0.0004999999375000026,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, ( 1.0,   1.0, &r), 0.4400505857449335160,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, ( 1.75,  1.0, &r), 0.168593922545763170103,     TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (30.0,   1.0, &r), 3.482869794251482902e-42,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (30.0, 100.0, &r), 0.08146012958117222297,    TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (10.0,   1.0, &r), 2.6306151236874532070e-10, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (10.0, 100.0, &r), -0.05473217693547201474,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (10.2, 100.0, &r), -0.03548919161046526864,   TEST_TOL2, GSL_SUCCESS);

  /* related to BUG#3 problem */
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (2.0, 900.0, &r),   -0.019974345269680646400,  TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Jnu_e, (2.0, 15000.0, &r), -0.0020455820181216382666, TEST_TOL4, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Ynu_e, (0.0001, 1.0, &r),  0.08813676933044478439,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (0.0001,10.0, &r),  0.05570979797521875261,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, ( 0.75,  1.0, &r), -0.6218694174429746383,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, ( 0.75, 10.0, &r),  0.24757063446760384953,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, ( 1.0, 0.001, &r), -636.6221672311394281,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, ( 1.0,   1.0, &r), -0.7812128213002887165,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (30.0,   1.0, &r), -3.0481287832256432162e+39, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (30.0, 100.0, &r),  0.006138839212010033452,   TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (10.0,   1.0, &r), -1.2161801427868918929e+08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (10.0, 100.0, &r),  0.05833157423641492875,    TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Ynu_e, (10.2, 100.0, &r),  0.07169383985546287091,    TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (0.0001,10.0, &r), 0.12783333709581669672,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, ( 1.0, 0.001, &r), 0.0004995003123542213370,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, ( 1.0,   1.0, &r), 0.20791041534970844887,    TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (30.0,   1.0, &r), 1.3021094983785914437e-42, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (30.0, 100.0, &r), 0.0004486987756920986146,  TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (10.0,   1.0, &r), 1.0127529864692066036e-10, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (10.0, 100.0, &r), 0.024176682718258828365,   TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_scaled_e, (10.2, 100.0, &r), 0.023691628843913810043,   TEST_TOL3, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Inu_e, (0.0001,10.0, &r), 2815.7166269770030352,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, ( 1.0, 0.001, &r), 0.0005000000625000026042,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, ( 1.0,   1.0, &r), 0.5651591039924850272,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, (30.0,   1.0, &r), 3.539500588106447747e-42,  TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, (30.0, 100.0, &r), 1.2061548704498434006e+40, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, (10.0,   1.0, &r), 2.7529480398368736252e-10, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, (10.0, 100.0, &r), 6.498975524720147799e+41,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Inu_e, (10.2, 100.0, &r), 6.368587361287030443e+41,  TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (0.0001,10.0, &r), 0.3916319346235421817, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, ( 1.0, 0.001, &r), 1000.9967345590684524, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, ( 1.0,   1.0, &r), 1.6361534862632582465, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (30.0,   1.0, &r), 1.2792629867539753925e+40, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (30.0, 100.0, &r), 10.673443449954850040, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (10.0,   1.0, &r), 4.912296520990198599e+08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (10.0, 100.0, &r), 0.20578687173955779807, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (10.0, 1000.0, &r), 0.04165905142800565788, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (10.0, 1.0e+8, &r), 0.00012533147624060789938, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_scaled_e, (10.2, 100.0, &r), 0.20995808355244385075, TEST_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_Knu_e, (0.0001,0.001, &r), 7.023689431812884141,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (0.0001,10.0, &r), 0.000017780062324654874306, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, ( 1.0, 0.001, &r), 999.9962381560855743,      TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, ( 1.0,   1.0, &r), 0.6019072301972345747,     TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (10.0, 0.001, &r), 1.8579455483904008064e+38, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (10.0,   1.0, &r), 1.8071328990102945469e+08, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (10.0, 100.0, &r), 7.655427977388100611e-45,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (10.2, 100.0, &r), 7.810600225948217841e-45,  TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (30.0,   1.0, &r), 4.706145526783626883e+39,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_Knu_e, (30.0, 100.0, &r), 3.970602055959398739e-43,  TEST_TOL2, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (0.0001,1.0e-100, &r), 5.439794449319847, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (0.0001,0.0001, &r), 2.232835507214331, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (0.0001,10.0, &r), -10.93743282256098, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, ( 1.0, 1.0e-100, &r), 230.2585092994045, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, ( 1.0, 1.0e-10, &r), 23.025850929940456840, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, ( 1.0, 0.001, &r), 6.907751517131146, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, ( 1.0,   1.0, &r), -0.5076519482107523309, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (30.0, 1.0e-100, &r), 6999.113586185543475, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (30.0,   1.0, &r), 91.34968784026325464, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (30.0, 100.0, &r), -97.63224126416760932, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (100.0, 1.0e-100, &r), 23453.606706185466825, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (100.0, 1.0, &r), 427.7532510250188083, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (100.0, 100.0, &r), -55.53422771502921431, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (1000.0, 1.0e-100, &r), 236856.183755993135, TEST_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (10000.0, 1.0e-100, &r), 2.39161558914890695e+06, TEST_TOL0, GSL_SUCCESS);

  /* [bug #31528] gsl_sf_bessel_lnKnu overflows for large nu */

  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (180.0, 2.2, &r), 735.1994170369583930752590258, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_lnKnu_e, (3500.5, 1500.0, &r), 1731.220077116482710070986699, TEST_TOL1, GSL_SUCCESS);

  sa = 0;
  gsl_sf_bessel_Jn_array(3, 38, 1.0, J);
  sa += ( test_sf_frac_diff(J[0],  0.0195633539826684059190  ) > TEST_TOL1);
  sa += ( test_sf_frac_diff(J[1],  0.0024766389641099550438  ) > TEST_TOL1);
  sa += ( test_sf_frac_diff(J[10], 1.9256167644801728904e-14 ) > TEST_TOL1);
  sa += ( test_sf_frac_diff(J[35], 6.911232970971626272e-57  ) > TEST_TOL1);
  gsl_test(sa, "  gsl_sf_bessel_Jn_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_Yn_array(3, 38, 1.0, Y);
  sa += ( test_sf_frac_diff(Y[0],  -5.821517605964728848      ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[1],  -33.27842302897211870      ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[10], -1.2753618701519837595e+12 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[35], -1.2124435663593357154e+54 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_Yn_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_In_scaled_array(3, 38, 1.0, I);
  sa += ( test_sf_frac_diff(I[0],  0.0081553077728142938170  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(I[1],  0.0010069302573377758637  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(I[10], 7.341518665628926244e-15  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(I[35], 2.5753065298357542893e-57 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_In_scaled_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_In_array(3, 38, 1.0, Y);
  sa += ( test_sf_frac_diff(Y[0],   0.0221684249243319024760  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[1],   0.0027371202210468663251  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[10],  1.9956316782072007564e-14 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[35],  7.000408942764452901e-57  ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_In_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_Kn_array(3, 38, 1.0, K);
  sa += ( test_sf_frac_diff(K[0],  7.101262824737944506 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[1],  44.23241584706284452 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[10], 1.9215763927929940846e+12 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[35], 1.8789385023806051223e+54 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_Kn_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_Kn_scaled_array(3, 38, 1.0, K);
  sa += ( test_sf_frac_diff(K[0],  19.303233695596904277 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[1],  120.23617222591483717 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[10], 5.223386190525076473e+12 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[35], 5.107484387813251411e+54 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_Kn_scaled_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_jl_array(50, 1.0, J);
  sa += ( test_sf_frac_diff(J[0],  0.84147098480789650670   ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(J[1],  0.30116867893975678925   ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(J[10], 7.116552640047313024e-11 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(J[50], 3.615274717489787311e-81 ) > TEST_TOL2 );
  gsl_test(sa, "  gsl_sf_bessel_jl_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_jl_steed_array(99, 1.0, J);
  sa += ( test_sf_frac_diff(J[0],  0.84147098480789650670   ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(J[1],  0.30116867893975678925   ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(J[10], 7.116552640047313024e-11 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(J[50], 3.615274717489787311e-81 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(J[80], 1.136352423414503264e-144 ) > TEST_TOL1 );
  gsl_test(sa, "  gsl_sf_bessel_jl_steed_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_yl_array(50, 1.0, Y);
  sa += ( test_sf_frac_diff(Y[0],  -0.5403023058681397174 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[1],  -1.3817732906760362241 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[10], -6.722150082562084436e+08  ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(Y[50], -2.7391922846297571576e+78 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_yl_array");
  s += sa;

  {
    double Y0[1];
    sa = 0;
    gsl_sf_bessel_yl_array(0, 1.0, Y0);
    sa += ( test_sf_frac_diff(Y0[0],  -0.5403023058681397174 ) > TEST_TOL0 );
    gsl_test(sa, "  gsl_sf_bessel_yl_array (lmax=0)");
    s += sa;
  }

  sa = 0;
  gsl_sf_bessel_il_scaled_array(50, 1.0, I);
  sa += ( test_sf_frac_diff(I[0],  0.43233235838169365410 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[1],  0.13533528323661269189 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[10], 2.7343719371837065460e-11 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[50], 1.3429606061892023653e-81 ) > TEST_TOL2 );
  gsl_test(sa, "  gsl_sf_bessel_il_scaled_array");
  s += sa;

  sa = 0;
  gsl_sf_bessel_il_scaled_array(50, 0.0, I);
  sa += ( test_sf_frac_diff(I[0],  1.0 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[1],  0.0 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[10], 0.0 ) > TEST_TOL2 );
  sa += ( test_sf_frac_diff(I[50], 0.0 ) > TEST_TOL2 );
  gsl_test(sa, "  gsl_sf_bessel_il_scaled_array (L=0)");
  s += sa;

  sa = 0;
  gsl_sf_bessel_kl_scaled_array(50, 1.0, K);
  sa += ( test_sf_frac_diff(K[0],  1.5707963267948966192     ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[1],  3.1415926535897932385     ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[10], 2.7231075458948147010e+09 ) > TEST_TOL0 );
  sa += ( test_sf_frac_diff(K[50], 1.1578440432804522544e+79 ) > TEST_TOL0 );
  gsl_test(sa, "  gsl_sf_bessel_kl_scaled_array");
  s += sa;

  {
    double K0[1];
    sa = 0;
    gsl_sf_bessel_kl_scaled_array(0, 1.0, K0);
    sa += ( test_sf_frac_diff(K[0],  1.5707963267948966192     ) > TEST_TOL0 );
    gsl_test(sa, "  gsl_sf_bessel_kl_scaled_array (lmax=0)");
    s += sa;
  }

  sa = 0;
  sa += ( gsl_sf_bessel_zero_J0_e(0, &r) != GSL_EINVAL );
  sa += ( r.val != 0.0 );
  s += sa;
  TEST_SF(s, gsl_sf_bessel_zero_J0_e, ( 1,  &r),  2.404825557695771, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J0_e, ( 2,  &r),  5.520078110286304, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J0_e, (20,  &r), 62.048469190227081, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J0_e, (25,  &r), 77.756025630388058, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J0_e, (100, &r), 313.37426607752784, TEST_TOL1, GSL_SUCCESS);

  sa = 0;
  sa += ( gsl_sf_bessel_zero_J1_e(0, &r) != GSL_SUCCESS );
  sa += ( r.val != 0.0 );
  s += sa;
  TEST_SF(s, gsl_sf_bessel_zero_J1_e, ( 1,  &r), 3.831705970207512, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J1_e, ( 2,  &r), 7.015586669815619, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J1_e, (20,  &r), 63.61135669848124, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J1_e, (25,  &r), 79.32048717547630, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_J1_e, (100, &r), 314.9434728377672, TEST_TOL2, GSL_SUCCESS);

  sa = 0;
  sa += ( gsl_sf_bessel_zero_Jnu_e(0.0, 0, &r) != GSL_EINVAL );
  sa += ( r.val != 0.0 );
  s += sa;
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (0.0,  1,  &r),  2.404825557695771, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (0.0,  2,  &r),  5.520078110286304, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (0.0, 20,  &r), 62.048469190227081, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (0.0, 25,  &r), 77.756025630388058, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (0.0, 100, &r), 313.37426607752784, TEST_TOL1, GSL_SUCCESS);

  sa = 0;
  sa += ( gsl_sf_bessel_zero_Jnu_e(1.0, 0, &r) != GSL_SUCCESS );
  sa += (r.val != 0.0);
  s += sa;
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 1, &r),  4.4934094579090641, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 1, &r),  8.7714838159599540, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 2, &r),  7.7252518369377072, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 2, &r),  12.338604197466944, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 3, &r),  10.904121659428900, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 3, &r),  15.700174079711671, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 4, &r),  14.066193912831473, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 4, &r),  18.980133875179921, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 5, &r),  17.220755271930768, TEST_TOL1, GSL_SUCCESS);

  /* Something wrong with the tolerances on these */
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 5, &r),  22.217799896561268, TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 8.0, 5, &r),  26.266814641176644, TEST_SQRT_TOL0, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (20.0, 5, &r),  41.413065513892636, TEST_SQRT_TOL0, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 6, &r),  20.371302959287563, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 6, &r),  25.430341154222704, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 8.0, 6, &r),  29.545659670998550, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 7, &r),  23.519452498689007, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 7, &r),  28.626618307291138, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 8.0, 7, &r),  32.795800037341462, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 8, &r),  26.666054258812674, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 8, &r),  31.811716724047763, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (10.0, 8, &r),  38.761807017881651, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 9, &r),  29.811598790892959, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 9, &r),  34.988781294559295, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (10.0, 9, &r),  42.004190236671805, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 1.5, 10, &r),  32.956389039822477, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 10, &r),  38.159868561967132, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 10, &r),  52.017241278881633, TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 11,  &r), 41.326383254047406, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 11,  &r), 55.289204146560061, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 12,  &r), 44.4893191232197314, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 12,  &r), 58.5458289043850856, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 13,  &r), 47.6493998066970948, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 13,  &r), 61.7897598959450550, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 14,  &r), 50.8071652030063595, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 14,  &r), 65.0230502510422545, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 15,  &r), 53.9630265583781707, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 15,  &r), 68.2473219964207837, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 16,  &r), 57.1173027815042647, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 16,  &r), 71.4638758850226630, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 17,  &r), 60.2702450729428077, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 17,  &r), 74.6737687121404241, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 18,  &r), 63.4220540458757799, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 18,  &r), 77.8778689734863729, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 19,  &r), 66.5728918871182703, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 19,  &r), 81.0768977206328326, TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 5.0, 20,  &r), 69.722891161716742,  TEST_TOL1, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (15.0, 20,  &r), 84.271459069716442,  TEST_TOL1, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 23.0, 11,  &r), 65.843393469524653, TEST_TOL6, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 30.0, 11,  &r), 74.797306585175426, TEST_TOL6, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 32.0, 15,  &r), 90.913637691861741, TEST_TOL6, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 50.0, 15,  &r), 113.69747988073942, TEST_TOL6, GSL_SUCCESS);

  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (  5.0, 22,  &r), 76.020793430591605, TEST_TOL2, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 10.0, 22,  &r), 83.439189796105756, TEST_TOL3, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, ( 12.0, 22,  &r), 86.345496520534055, TEST_TOL6, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (100.0, 22,  &r), 199.82150220122519, TEST_TOL4, GSL_SUCCESS);
  TEST_SF(s, gsl_sf_bessel_zero_Jnu_e, (500.0, 22,  &r), 649.34132440891735, TEST_TOL2, GSL_SUCCESS);

  sa = 0;
  for(i=0; i<100; i++) {
    J[i] = i/10.0;
  }
  gsl_sf_bessel_sequence_Jnu_e(2.0, GSL_MODE_DEFAULT, 100, J);
  sa += ( test_sf_frac_diff(J[1],  0.0012489586587999188454 ) > TEST_TOL1 );
  sa += ( test_sf_frac_diff(J[20], 0.3528340286156377192 ) > TEST_TOL4 );
  sa += ( test_sf_frac_diff(J[50], 0.0465651162777522155 ) > TEST_TOL4 );
  gsl_test(sa, "  gsl_sf_sequence_Jnu_e(2)");
  s += sa;

  sa = 0;
  for(i=0; i<100; i++) {
    J[i] = i;
  }
  gsl_sf_bessel_sequence_Jnu_e(12.0, GSL_MODE_DEFAULT, 100, J);
  sa += ( test_sf_frac_diff(J[1],   4.999718179448405289e-13 ) > TEST_TOL1 );
  sa += ( test_sf_frac_diff(J[5],   7.627813166084551355e-05 ) > TEST_TOL1 );
  sa += ( test_sf_frac_diff(J[7],   2.655620035894568062e-03 ) > TEST_TOL3 );
  sa += ( test_sf_frac_diff(J[10],  6.337025497015601509e-02 ) > TEST_TOL3 );
  sa += ( test_sf_frac_diff(J[15],  0.23666584405476805591 ) > TEST_TOL3 );
  sa += ( test_sf_frac_diff(J[30],  0.14825335109966010021 ) > TEST_TOL3 );
  sa += ( test_sf_frac_diff(J[70],  0.04109913716555364262 ) > TEST_TOL4 );
  sa += ( test_sf_frac_diff(J[99], -0.0015052760501176728  ) > TEST_TOL5 );
  gsl_test(sa, "  gsl_sf_sequence_Jnu_e(12)");
  s += sa;

  sa = 0;
  for(i=0; i<100; i++) {
    J[i] = i * 20;
  }
  gsl_sf_bessel_sequence_Jnu_e(1000.0, GSL_MODE_DEFAULT, 100, J);
  sa += ( test_sf_frac_diff(J[50],  0.04473067294796404088 ) > TEST_TOL5 );
  sa += ( test_sf_frac_diff(J[99],  0.01400619760349853902 ) > TEST_TOL6 );
  gsl_test(sa, "  gsl_sf_sequence_Jnu_e(1000)");
  s += sa;

  return s;
}

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include <utilities.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_stack_code.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_geo3d_vector.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"
#include "tz_neurotrace.h"
#include "tz_stack_threshold.h"
#include "tz_iarray.h"
#include "tz_darray.h"

/*
 * trace_neuron - trace neuron from given seeds
 *
 * trace_neuron [!wtr] seed_file -Dsave_dir
 *   -r: write intermediate results
 *
 */
int main(int argc, char* argv[])
{
  static char *Spec[] = {
    "[!wtr] [-canvas <string>] [-mask <string>] [-res <string>] [-minr <int>]",
    "-minlen <double>",
    " <image:string> -S<string> -D<string>",
    NULL};

  Process_Arguments(argc, argv, Spec, 1);

  char *dir = Get_String_Arg("-D");

  char file_path[100];
  sprintf(file_path, "%s/%s", dir, Get_String_Arg("-S"));
  printf("%s\n", file_path);

  Geo3d_Scalar_Field *seed = Read_Geo3d_Scalar_Field(file_path);

  int idx;

  sprintf(file_path, "%s/%s.bn", dir, "max_r");
  double max_r;
  int tmp;
  if (fexist(file_path)) {
    darray_read2(file_path, &max_r, &tmp);
  } else {
    max_r = darray_max(seed->values, seed->size, &idx);
  }

  printf("%g\n", max_r);

  max_r *= 1.5;

  /*
  sprintf(file_path, "%s/%s", dir, "soma0.bn");
  if (!fexist(file_path)) {
    max_r *= 2.0;
  }
  */

  Set_Neuroseg_Max_Radius(max_r);

  Stack *signal = Read_Stack(Get_String_Arg("image"));

  dim_type dim[3];
  dim[0] = signal->width;
  dim[1] = signal->height;
  dim[2] = signal->depth;
  /*
  IMatrix *chord = Make_IMatrix(dim, 3);

  Stack *code = Make_Stack(GREY16,
			   signal->width, signal->height, signal->depth);
  */
  Rgb_Color color;
  Set_Color(&color, 255, 0, 0);

  Stack *canvas = NULL;

  char trace_file_path[100];
  sprintf(trace_file_path, "%s/%s", dir, Get_String_Arg("-canvas"));

  if (fexist(trace_file_path) == 1) {
    canvas = Read_Stack((char *) trace_file_path);
  } else {
    canvas = Copy_Stack(signal);
    Stretch_Stack_Value_Q(canvas, 0.999);
    Translate_Stack(canvas, COLOR, 1);
  }

  Stack *traced = NULL;

  char trace_mask_path[100];
  sprintf(trace_mask_path, "%s/%s", dir, Get_String_Arg("-mask"));

  if (fexist(trace_mask_path) == 1) {
    traced = Read_Stack((char *) trace_mask_path);
  } else {
    traced = Make_Stack(GREY, signal->width, signal->height, signal->depth);
    One_Stack(traced);
  }

  //Object_3d *obj = NULL;
  int seed_offset = -1;

  Neurochain *chain = NULL;

  double z_scale = 1.0;

  if (Is_Arg_Matched("-res")) {
    sprintf(file_path, "%s", Get_String_Arg("-res"));

    if (fexist(file_path)) {
      double res[3];
      int length;
      darray_read2(file_path, res, &length);
      if (res[0] != res[1]) {
	perror("Different X-Y resolutions.");
	TZ_ERROR(ERROR_DATA_VALUE);
      }
      z_scale = res[0] / res[2];
    }
  }

  //sprintf(file_path, "%s/%s", dir, Get_String_Arg("-M"));
  //Stack *stack = Read_Stack(file_path);

  tic();

  FILE *fp = NULL;
  char chain_file_path[100];
  char vrml_file_path[100];

  double min_chain_length = 25.0;

  if (Is_Arg_Matched("-minlen")) {
    min_chain_length = Get_Double_Arg("-minlen");
  }

  int *indices = iarray_malloc(seed->size);
  double *values = darray_malloc(seed->size);
  int i;

  Local_Neuroseg *locseg = (Local_Neuroseg *)
    malloc(seed->size * sizeof(Local_Neuroseg));

  int index = 0;
  for (i = 0; i < seed->size; i++) {
    printf("-----------------------------> seed: %d / %d\n", i, seed->size);
    indices[i] = i;
    index = i;
    int x = (int) seed->points[index][0];
    int y = (int) seed->points[index][1];
    int z = (int) seed->points[index][2];

    double width = seed->values[index];

    chain = New_Neurochain();

    seed_offset = Stack_Util_Offset(x, y, z, signal->width, signal->height,
				    signal->depth);

    if (width < 3.0) {
      width += 0.5;
    }
    Set_Neuroseg(&(locseg[i].seg), width, width, 12.0,
		 0.0, 0.0, 0.0);

    double cpos[3];
    cpos[0] = x;
    cpos[1] = y;
    cpos[2] = z;
    cpos[2] *= z_scale;

    Set_Neuroseg_Position(&(locseg[i]), cpos, NEUROSEG_CENTER);
    Stack_Fit_Score fs;
    fs.n = 1;
    fs.options[0] = 1;
    values[i] = Local_Neuroseg_Orientation_Search_C(&(locseg[i]), signal, z_scale, &fs);
  }

  darray_qsort(values, indices, seed->size);

  /*
  for (i = 0; i < seed->size; i++) {
    indices[i] = i;
  }
  darraycpy(values, seed->values, 0, seed->size);
  darray_qsort(values, indices, seed->size);
  */

  int counter = 0;

  //  for (i = seed->size - 1; i >= seed->size - 231; i--) {
  for (i = seed->size - 1; i >= 0; i--) {
    index = indices[i];

    printf("-----------------------------> seed: %d / %d\n", i, seed->size);

    sprintf(chain_file_path, "%s/chain%d.bn", dir, index);
    sprintf(vrml_file_path, "%s/chain%d.wrl", dir, index);

    if (fexist(chain_file_path) == 1) {
      chain = Read_Neurochain(chain_file_path);
      if (Neurochain_Geolen(chain) >= min_chain_length) {
	Write_Neurochain_Vrml(vrml_file_path, chain);
	Neurochain_Label(canvas, chain, z_scale);
	Neurochain_Erase_E(traced, chain, z_scale, 0,
			   Neurochain_Length(chain, FORWARD),
			   1.5, 0.0);
      }

      Free_Neurochain(chain);
      printf("chain exists\n");
      continue;
    }

    int x = (int) seed->points[index][0];
    int y = (int) seed->points[index][1];
    int z = (int) seed->points[index][2];

    if (*STACK_PIXEL_8(traced, x, y, z, 0) == 0) {
      printf("traced \n");
      continue;
    }

    double width = seed->values[index];

    if (width > max_r) {
      printf("too thick\n");
      continue;
    }

    if (Is_Arg_Matched("-minr")) {
      int max_level = (int) (width + 0.5);
      if (max_level <= Get_Int_Arg("-minr")) {
	printf("too thin\n");
	continue;
      }
    }
    /*
    seed_offset = Stack_Util_Offset(x, y, z, signal->width, signal->height,
				    signal->depth);
    */

    chain = New_Neurochain();
    /*
    Stack_Level_Code_Constraint(stack, code, chord->array, &seed_offset, 1,
				max_level + 1);

    Voxel_t v;
    v[0] = x;
    v[1] = y;
    v[2] = z;

    Stack *tmp_stack = Copy_Stack(stack);
    obj = Stack_Grow_Object_Constraint(tmp_stack, 1, v, chord, code,
				       max_level);
    Free_Stack(tmp_stack);

    Print_Object_3d_Info(obj);

    double vec[3];
    Object_3d_Orientation_Zscale(obj, vec, MAJOR_AXIS, z_scale);

    double theta, psi;
    Geo3d_Vector obj_vec;
    Set_Geo3d_Vector(&obj_vec, vec[0], vec[1], vec[2]);

    Geo3d_Vector_Orientation(&obj_vec, &theta, &psi);
    */

    /*
    if (width < 3.0) {
      width += 0.5;
    }
    Set_Neuroseg(&(chain->locseg.seg), width, width, 12.0,
		 0.0, 0.0, 0.0);

    double cpos[3];
    cpos[0] = x;
    cpos[1] = y;
    cpos[2] = z;
    cpos[2] *= z_scale;

    //Set_Neuroseg_Position(&(chain->locseg), cpos, NEUROSEG_BOTTOM);
    Set_Neuroseg_Position(&(chain->locseg), cpos, NEUROSEG_CENTER);
    Stack_Fit_Score fs;
    fs.n = 1;
    fs.options[0] = 1;
    Local_Neuroseg_Orientation_Search_C(&(chain->locseg), signal, z_scale,
					&fs);
    //fs.options[0] = 1;
    */

    Copy_Local_Neuroseg(&(chain->locseg), &(locseg[index]));
    Neurochain *chain_head = chain;

    if (Initialize_Tracing(signal, chain, NULL, z_scale) >= MIN_SCORE) {
      if ((Neuroseg_Hit_Traced(&(chain->locseg), traced, z_scale) == FALSE) &&
	  (chain->locseg.seg.r1 < max_r) &&
	  (chain->locseg.seg.r2 < max_r)) {
	//Initialize_Tracing(signal, chain, NULL, z_scale);
	chain = Trace_Neuron2(signal, chain, BOTH, traced, z_scale, 500);

	//Neurochain *chain_head = Neurochain_Head(chain);
	chain_head = Neurochain_Remove_Overlap_Segs(chain);
	chain_head = Neurochain_Remove_Turn_Ends(chain_head, 0.5);
	/*
	if (i == seed->size - 231) {
	  Print_Neurochain(chain_head);
	}
	*/

	fp = fopen(chain_file_path, "w");
	Neurochain_Fwrite(chain_head, fp);
	fclose(fp);
	if (Neurochain_Geolen(chain_head) >= min_chain_length) {
	  Write_Neurochain_Vrml(vrml_file_path, chain_head);

	  Neurochain_Erase_E(traced, chain_head, z_scale, 0,
			     Neurochain_Length(chain_head, FORWARD),
			     1.5, 0.0);
	  Neurochain_Label(canvas, chain_head, z_scale);

	  counter += Neurochain_Length(chain_head, FORWARD);
	  if (counter > 500) {
	    if (Is_Arg_Matched("-r")) {
	      Write_Stack((char *) trace_mask_path, traced);
	    }

	    if (Is_Arg_Matched("-r")) {
	    Write_Stack((char *) trace_file_path, canvas);
	    }

	    counter = 0;
	  }
	}
      }
    }

    Free_Neurochain(chain_head);

    //Kill_Object_3d(obj);
  }

  Write_Stack((char *) trace_file_path, canvas);
  if (Is_Arg_Matched("-r")) {
    Write_Stack((char *) trace_mask_path, traced);
  }

  Kill_Geo3d_Scalar_Field(seed);

  printf("Time passed: %lld\n", toc());

  return 0;
}

ext: c


text:
/*
 * Copyright (C) 2015  University of Oregon
 *
 * You may distribute under the terms of either the GNU General Public
 * License or the Apache License, as specified in the LICENSE file.
 *
 * For more information, see the LICENSE file.
 */
#ifndef VNMRFDIST_H
#define VNMRFDIST_H

#ifdef VNMR_GPL
#include <gsl/gsl_cdf.h>
#define f_distribution(f,n1,n2) gsl_cdf_fdist_P( (f), (n1), (n2) )

#else

extern double f_distribution(double f, double n1, double n2);
#endif

#endif

ext: h

max_stars_repo_name: DanIverson/OpenVnmrJ

max_issues_repo_name: timburrow/openvnmrj-source

max_forks_repo_name: timburrow/openvnmrj-source


text:
#include "lah.h"

#ifdef HAVE_LAPACK
#include <lapacke.h>

/* General Singular Value Decomposition*/
/* For now get all eigen values jobz ='A',
 * Can use this parameter later for truncated version */
lah_Return lah_SVD(lah_mat *A, lah_mat *U, lah_value *S, lah_mat *Vt)
{
    lapack_int res;

    if (A == NULL || U == NULL || S == NULL || Vt == NULL
        || A->nR != A->nC || A->nR != U->nR || A->nC != Vt->nC )
    {
        return lahReturnParameterError;
    }
    /*
    matrix_layout = LAH_LAPACK_LAYOUT(A);
    if (matrix_layout == LAPACK_COL_MAJOR)
    {
        lda = A->nR;
        ldu = U->nR;
        ldv = Vt->nR;
    }
    else
    {
        lda = A->nC;
        ldu = U->nC;
        ldv = Vt->nC;
    }
    */

    /* Do the actual SVD: */
    res = GESDD(LAH_LAPACK_LAYOUT, 'A', A->nC, A->nR, A->data, LAH_LEADING_DIM(A), S, U->data, LAH_LEADING_DIM(U), Vt->data, LAH_LEADING_DIM(Vt));
    return (res == 0) ? lahReturnOk : lahReturnExternError;
}
#endif


text:
// lasp_math_raw.h
//
// Author: J.A. de Jong - ASCEE
//
// Description: Raw math routines working on raw arrays of floats and
// complex numbers.
//////////////////////////////////////////////////////////////////////
#pragma once
#ifndef LASP_MATH_RAW_H
#define LASP_MATH_RAW_H
#include "lasp_assert.h"
#include "lasp_tracer.h"
#include "lasp_types.h"
#include <float.h>
#include <math.h>

#if LASP_USE_BLAS == 1
#include <cblas.h>
#elif LASP_USE_BLAS == 0
#else
#error "LASP_USE_BLAS should be set to either 0 or 1"
#endif

#ifdef LASP_DOUBLE_PRECISION
#define c_real creal
#define c_imag cimag
#define d_abs fabs
#define c_abs cabs
#define c_conj conj
#define d_atan2 atan2
#define d_acos acos
#define d_sqrt sqrt
#define c_exp cexp
#define d_sin sin
#define d_cos cos
#define d_pow pow
#define d_log10 log10
#define d_ln log
#define d_epsilon (DBL_EPSILON)

#else  // LASP_DOUBLE_PRECISION not defined
#define c_conj conjf
#define c_real crealf
#define c_imag cimagf
#define d_abs fabsf
#define c_abs cabsf
#define d_atan2 atan2f
#define d_acos acosf
#define d_sqrt sqrtf
#define c_exp cexpf
#define d_sin sinf
#define d_cos cosf
#define d_pow powf
#define d_log10 log10f
#define d_ln logf
#define d_epsilon (FLT_EPSILON)
#endif // LASP_DOUBLE_PRECISION

/// Positive infinite
#define d_inf (INFINITY)

#ifdef M_PI
static const d number_pi = M_PI;
#else
static const d number_pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
#endif

/**
 * Set all elements in an array equal to val
 *
 * @param to
 * @param val
 * @param size
 */
static inline void d_set(d to[],d val,us size) {
    for(us i=0;i<size;i++) {
        to[i]=val;
    }
}
/**
 * Set all elements in an array equal to val
 *
 * @param to
 * @param val
 * @param size
 */
static inline void c_set(c to[],c val,us size) {
    for(us i=0;i<size;i++) {
        to[i]=val;
    }
}

/**
 * Return the maximum of two doubles
 *
 * @param a value 1
 * @param b value 2
 *
 * @returns the maximum of value 1 and 2
 */
static inline d d_max(const d a,const d b) {
    return a>b?a:b;
}

/**
 * Return the dot product of two arrays, one of them complex-valued,
 * the other real-valued
 *
 * @param a the complex-valued array
 * @param b the real-valued array
 * @param size the size of the arrays. *Should be equal-sized!*
 *
 * @return the dot product
 */
static inline c cd_dot(const c a[],const d b[],us size){
    c result = 0;
    us i;
    for(i=0;i<size;i++){
        result+=a[i]*b[i];
    }
    return result;
}
/**
 * Return the dot product of two complex-valued arrays. Wraps BLAS
 * when LASP_USE_BLAS == 1.
 *
 * @param a complex-valued array
 * @param b complex-valued array
 * @param size the size of the arrays. *Should be equal-sized!*
 *
 * @return the dot product
 */
static inline c cc_dot(const c a[],const c b[],us size){
    #if LASP_USE_BLAS == 1
    WARN("CBlas zdotu not yet tested");
    #if LASP_DOUBLE_PRECISION
    // assert(0);
    return cblas_zdotu(size,(d*) a,1,(d*) b,1);
    #else
    return cblas_cdotu(size,(d*) a,1,(d*) b,1);
    #endif
    #else
    c result = 0;
    us i;
    for(i=0;i<size;i++){
        result+=a[i]*b[i];
    }
    return result;
    #endif
}

/**
 * Compute the dot product of two real arrays.
 *
 * @param a First array.
 * @param b Second array.
 * @param size Size of the arrays.
 * @return The result.
 */
static inline d d_dot(const d a[],const d b[],const us size){
    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    return cblas_ddot(size,a,1,b,1);
    #else  // Single precision function
    return cblas_sdot(size,a,1,b,1);
    #endif
    #else  // No BLAS, do it manually
    d result = 0;
    us i;
    for(i=0;i<size;i++){
        result+=a[i]*b[i];
    }
    return result;
    #endif
}

/**
 * Copy array of floats.
 *
 * @param to : Array to write to
 * @param from : Array to read from
 * @param size : Size of arrays
 * @param to_inc : Mostly equal to 1, the stride of the array to copy to
 * @param from_inc : Mostly equal to 1, the stride of the array to copy from
 */
static inline void d_copy(d to[],
                          const d from[],
                          const us size,
                          const us to_inc,
                          const us from_inc){
    #if LASP_USE_BLAS == 1
    cblas_dcopy(size,from,from_inc,to,to_inc);
    #else
    us i;
    for(i=0;i<size;i++)
        to[i*to_inc] = from[i*from_inc];
    #endif
}

/**
 * Copy array of floats to array of complex floats. Imaginary part set
 * to zero.
 *
 * @param to : Array to write to
 * @param from : Array to read from
 * @param size : Size of arrays
 */
static inline void cd_copy(c to[],const d from[],const us size) {
    us i;
    for(i=0;i<size;i++) {
        to[i] = (c) (from[i]);
        dbgassert(cimag(to[i]) == 0,"Imaginary part not equal to zero");
    }
}

/**
 * Copy float vector to complex vector. Imaginary part set
 * to zero.
 *
 * @param to : Vector to write to
 * @param from : Vector to read from
 */
static inline void c_copy(c to[],const c from[],const us size){

    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    cblas_zcopy(size,(d*) from,1,(d*) to,1);
    #else
    cblas_ccopy(size,(d*) from,1,(d*) to,1);
    #endif
    #else
    us i;
    for(i=0;i<size;i++)
        to[i] = from[i];
    #endif
}
/**
 * Multiply y with fac, and add result to x
 *
 * @param x[in,out] Array to add to
 * @param y[in] Array to add to x
 * @param[in] fac Factor with which to multiply y
 * @param[in] size Size of the arrays
 */
static inline void d_add_to(d x[],const d y[],
                            const d fac,const us size){
    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    cblas_daxpy(size,fac,y,1,x,1);
    #else
    cblas_saxpy(size,fac,y,1,x,1);
    #endif
    #else
    us i;
    for(i=0;i<size;i++)
        x[i]+=fac*y[i];
    #endif
}
/**
 * x = x + fac*y
 *
 * @param[in,out] x Array to add to
 * @param[in] y Array to add to x
 * @param[in] fac Factor with which to multiply y
 * @param[in] size Size of the arrays
 */
static inline void c_add_to(c x[],const c y[],
                            const c fac,const us size){
    fsTRACE(15);
    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    cblas_zaxpy(size,(d*) &fac,(d*) y,1,(d*) x,1);
    #else
    cblas_caxpy(size,(d*) &fac,(d*) y,1,(d*) x,1);
    #endif
    #else
    us i;
    for(i=0;i<size;i++)
        x[i]+=fac*y[i];
    #endif
    feTRACE(15);
}

/**
 * Scale an array of doubles
 *
 * @param a array
 * @param scale_fac scale factor
 * @param size size of the array
 */
static inline void d_scale(d a[],const d scale_fac,us size){
    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    cblas_dscal(size,scale_fac,a,1);
    #else
    cblas_sscal(size,scale_fac,a,1);
    #endif
    #else
    us i;
    for(i=0;i<size;i++)
        a[i]*=scale_fac;
    #endif
}

/**
 * Scale an array of complex floats
 *
 * @param a array
 * @param scale_fac scale factor
 * @param size size of the array
 */
static inline void c_scale(c a[],const c scale_fac,us size){
    #if LASP_USE_BLAS == 1
    // Complex argument should be given in as array of two double
    // values. The first the real part, the second the imaginary
    // part. Fortunately the (c) type stores the two values in this
    // order. To be portable and absolutely sure anything goes well,
    // we convert it explicitly here.

    #if LASP_DOUBLE_PRECISION
    cblas_zscal(size,(d*) &scale_fac,(d*) a,1);
    #else
    cblas_cscal(size,(d*) &scale_fac,(d*) a,1);
    #endif
    #else
    us i;
    for(i=0;i<size;i++)
        a[i]*=scale_fac;
    #endif
}

/**
 * Compute the maximum value of an array
 *
 * @param a array
 * @param size size of the array
 * @return maximum
 */
static inline d darray_max(const d a[],us size){
    us i;
    d max = a[0];
    for(i=1;i<size;i++){
        if(a[i] > max) max=a[i];
    }
    return max;
}
/**
 * Compute the minimum of an array
 *
 * @param a array
 * @param size size of the array
 * @return minimum
 */
static inline d d_min(const d a[],us size){
    us i;
    d min = a[0];
    for(i=1;i<size;i++){
        if(a[i] > min) min=a[i];
    }
    return min;
}

/**
 * Compute the \f$ L_2 \f$ norm of an array of doubles
 *
 * @param a Array
 * @param size Size of array
 */
static inline d d_norm(const d a[],us size){
    #if LASP_USE_BLAS == 1
    return cblas_dnrm2(size,a,1);
    #else
    d norm = 0;
    us i;
    for(i=0;i<size;i++){
        norm+=a[i]*a[i];
    }
    norm = d_sqrt(norm);
    return norm;
    #endif

}

/**
 * Compute the \f$ L_2 \f$ norm of an array of complex floats
 *
 * @param a Array
 * @param size Size of array
 */
static inline d c_norm(const c a[],us size){
    #if LASP_USE_BLAS == 1
    return cblas_dznrm2(size,(d*) a,1);
    #else
    d norm = 0;
    us i;
    for(i=0;i<size;i++){
        d absa = c_abs(a[i]);
        norm+=absa*absa;
    }
    norm = d_sqrt(norm);
    return norm;
    #endif

}

/**
 * Computes the element-wise vector product, or Hadamard product of
 * arr1 and arr2
 *
 * @param res Where the result will be stored
 * @param arr1 Array 1
 * @param vec2 Array 2
 * @param size: Size of the arrays
 */
void d_elem_prod_d(d res[],
                   const d arr1[],
                   const d arr2[],
                   const us size);

/**
 * Computes the element-wise vector product, or Hadamard product of
 * arr1 and arr2 for complex floats
 *
 * @param res Where the result will be stored
 * @param arr1 Array 1
 * @param vec2 Array 2
 * @param size: Size of the arrays
 */
void c_hadamard(c res[],
                const c arr1[],
                const c arr2[],
                const us size);

/**
 * Compute the complex conjugate of a complex vector and store the
 * result.
 *
 * @param res Result vector
 * @param in Input vector
 * @param size Size of the vector
 */
static inline void carray_conj(c res[],const c in[],const us size) {
    // First set the result vector to zero
    fsTRACE(15);
    c_set(res,0,size);
    #if LASP_USE_BLAS == 1
    #if LASP_DOUBLE_PRECISION
    // Cast as a float, scale all odd elements with minus one to find
    // the complex conjugate.
    cblas_daxpy(size ,1.0,(d*) in,2,(d*) res,2);
    cblas_daxpy(size,-1.0,&((d*) in)[1],2,&((d*) res)[1],2);
    #else
    cblas_faxpy(size ,1,(d*) in,2,(d*) res,2);
    cblas_faxpy(size,-1,&((d*) in)[1],2,&((d*) res)[1],2);
    #endif  // LASP_DOUBLE_PRECISION
    #else
    for(us i=0;i<size;i++) {
        res[i] = c_conj(in[i]);
    }
    #endif  // LASP_USE_BLAS
    feTRACE(15);
}
/**
 * In place complex conjugation
 *
 * @param res Result vector
 * @param size Size of the vector
 */
static inline void c_conj_inplace(c res[],us size) {
    #if LASP_USE_BLAS
    #if LASP_DOUBLE_PRECISION
    // Cast as a float, scale all odd elements with minus one to find
    // the complex conjugate.
    cblas_dscal(size,-1,&((d*) res)[1],2);
    #else
    cblas_sscal(size,-1,&((d*) res)[1],2);
    #endif  // LASP_DOUBLE_PRECISION
    #else
    for(us i=0;i<size;i++) {
        res[i] = c_conj(res[i]);
    }
    #endif  // LASP_USE_BLAS
}

#endif // LASP_MATH_RAW_H
//////////////////////////////////////////////////////////////////////

max_stars_repo_name: asceenl/lasp

max_issues_repo_name: asceenl/lasp

max_forks_repo_name: asceenl/lasp


text:
/*****************************************************************************
*
* Rokko: Integrated Interface for libraries of eigenvalue decomposition
*
* Copyright (C) 2012-2020 Rokko Developers https://github.com/t-sakashita/rokko
*
* Distributed under the Boost Software License, Version 1.0. (See accompanying
* file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
*
*****************************************************************************/

#ifndef ROKKO_SCALAPACK_SCALAPACK_INTERFACE_H
#define ROKKO_SCALAPACK_SCALAPACK_INTERFACE_H

#include <rokko/config.h>
#include <rokko/mangling.h>
#include <lapacke.h>

#ifdef __cplusplus
extern "C" {
#endif

#define SCALAPACK_descinit ROKKO_GLOBAL(descinit,DESCINIT)
void SCALAPACK_descinit(int* desc, const int* m, const int* n, const int* mb, const int* nb,
                        const int* irsrc, const int* icsrc, const int* ictxt, const int* lld,
                        int* info);

#define SCALAPACK_indxg2p ROKKO_GLOBAL(indxg2p,INDXG2P)
int SCALAPACK_indxg2p(int* indxglob, int* nb, int* iproc, int* isrcproc, int* nprocs);

#define SCALAPACK_numroc ROKKO_GLOBAL(numroc,NUMROC)
int SCALAPACK_numroc(int* n, int* nb, int* iproc, int* isrcproc, int* nprocs);

#define SCALAPACK_pdelget ROKKO_GLOBAL(pdelget,PDELGET)
void SCALAPACK_pdelget(const char* scope, const char* top, double* alpha,
                       const double* A, const int* ia, const int* ja, const int* descA);

#define SCALAPACK_pdelset ROKKO_GLOBAL(pdelset,PDELSET)
void SCALAPACK_pdelset(double* A, const int* ia, const int* ja, const int* descA,
                       const double* alpha);

#define SCALAPACK_pdlamch ROKKO_GLOBAL(pdlamch,PDLAMCH)
double SCALAPACK_pdlamch(const int* icnt, const char* cmch);

#define SCALAPACK_pdlange ROKKO_GLOBAL(pdlange,PDLANGE)
double SCALAPACK_pdlange(const char* norm, int* m, int* n, const double* A,
                         const int* ia, const int* ja, const int* descA, double* work);

#define SCALAPACK_pdlaprnt ROKKO_GLOBAL(pdlaprnt,PDLAPRNT)
void SCALAPACK_pdlaprnt(const int* m, const int* n, const double* A, const int* ia, const int* ja,
                        const int* descA, const int* irprnt, const int* icprnt, const char* cmatnm,
                        const int* nout, double* work);

#define SCALAPACK_PSYEV_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* jobz, const char* uplo, const int* n, \
                                 TYPE* A, const int* ia, const int* ja, const int* descA, \
                                 TYPE_REAL* w, TYPE* Z, const int* iz, const int* jz, const int* descZ, \
                                 TYPE* work, const int* lwork, int* info);

SCALAPACK_PSYEV_DECL(pssyev, PSSYEV, float, float);
SCALAPACK_PSYEV_DECL(pdsyev, PDSYEV, double, double);
SCALAPACK_PSYEV_DECL(pcheev, PCHEEV, lapack_complex_float, float);
SCALAPACK_PSYEV_DECL(pzheev, PZHEEV, lapack_complex_double, double);

#undef SCALAPACK_PSYEV_DECL

#define SCALAPACK_PSYEVD_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* jobz, const char* uplo, const int* n, \
  TYPE* A, const int* ia, const int* ja, const int* descA, \
  TYPE_REAL* w, \
  TYPE* Z, const int* iz, const int* jz, const int* descZ, \
  TYPE* work, const int* lwork, int* iwork, const int* liwork, int* info);

SCALAPACK_PSYEVD_DECL(pssyevd, PSSYEVD, float, float);
SCALAPACK_PSYEVD_DECL(pdsyevd, PDSYEVD, double, double);
SCALAPACK_PSYEVD_DECL(pcheevd, PCHEEVD, lapack_complex_float, float);
SCALAPACK_PSYEVD_DECL(pzheevd, PZHEEVD, lapack_complex_double, double);

#undef SCALAPACK_PSYEVD_DECL

#define SCALAPACK_PSYEVX_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* jobz, const char* range, const char* uplo, const int* n, \
  TYPE* A, const int* iA, const int* jA, const int* descA, \
  const TYPE_REAL* vl, const TYPE_REAL* vu, const int* il, const int* iu, \
  const TYPE_REAL* abstol, int* m, int* nZ, TYPE_REAL* w, const TYPE_REAL* orfac, \
  TYPE* Z, const int* iZ, const int* jZ, const int* descZ, \
  TYPE* work, const int* lwork, int* iwork, const int* liwork, \
  int* ifail, int* iclustr, TYPE_REAL* gap, int* info);

SCALAPACK_PSYEVX_DECL(pssyevx, PSSYEVX, float, float);
SCALAPACK_PSYEVX_DECL(pdsyevx, PDSYEVX, double, double);
SCALAPACK_PSYEVX_DECL(pcheevx, PCHEEVX, lapack_complex_float, float);
SCALAPACK_PSYEVX_DECL(pzheevx, PZHEEVX, lapack_complex_double, double);

#undef SCALAPACK_PSYEVX_DECL

#ifdef ROKKO_HAVE_SCALAPACK_PDSYEVR
#define SCALAPACK_PSYEVR_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* jobz, const char* range, const char* uplo, const int* n, \
  TYPE* A, const int* iA, const int* jA, const int* descA, \
  const TYPE_REAL* vl, const TYPE_REAL* vu, const int* il, const int* iu, \
  int* m, int* nz, \
  TYPE_REAL* w, TYPE* Z, const int* iz, const int* jz, const int* descZ, \
  TYPE* work, const int* lwork, int* iwork, const int* liwork, int* info);

SCALAPACK_PSYEVR_DECL(pssyevr, PSSYEVR, float, float);
SCALAPACK_PSYEVR_DECL(pdsyevr, PDSYEVR, double, double);
SCALAPACK_PSYEVR_DECL(pcheevr, PCHEEVR, lapack_complex_float, float);
SCALAPACK_PSYEVR_DECL(pzheevr, PZHEEVR, lapack_complex_double, double);

#undef SCALAPACK_PSYEVR_DECL
#endif

#define SCALAPACK_pselget ROKKO_GLOBAL(pselget,PSELGET)
void SCALAPACK_pselget(const char* scope, const char* top, float* alpha,
                       const float* A, const int* ia, const int* ja, const int* descA);

#define SCALAPACK_pselset ROKKO_GLOBAL(pselset,PSELSET)
void SCALAPACK_pselset(float* A, const int* ia, const int* ja, const int* descA,
                       const float* alpha);

#define SCALAPACK_pslamch ROKKO_GLOBAL(pslamch,PSLAMCH)
float SCALAPACK_pslamch(const int* icnt, const char* cmch);

#define SCALAPACK_pslange ROKKO_GLOBAL(pslange,PSLANGE)
float SCALAPACK_pslange(const char* norm, int* m, int* n, const float* A,
                        const int* ia, const int* ja, const int* descA, float* work);

#define SCALAPACK_pslaprnt ROKKO_GLOBAL(pslaprnt,PSLAPRNT)
void SCALAPACK_pslaprnt(const int* m, const int* n, const float* A, const int* ia, const int* ja,
                        const int* descA, const int* irprnt, const int* icprnt, const char* cmatnm,
                        const int* nout, float* work);

#define SCALAPACK_PSTEBZ_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const int* ictxt, const char* range, const char* order, const int* n, \
                                 const TYPE* vl, const TYPE* vu, const int* il, const int* iu, \
                                 const TYPE* abstol, const TYPE* d, const TYPE* e, int* m, int* nsplit, \
                                 TYPE* w, int* iblock, int* isplit, \
                                 TYPE* work, const int* lwork, int* iwork, const int* liwork, int* info);

SCALAPACK_PSTEBZ_DECL(psstebz, PSSTEBZ, float);
SCALAPACK_PSTEBZ_DECL(pdstebz, PDSTEBZ, double);

#undef SCALAPACK_PSTEBZ_DECL

#define SCALAPACK_PSTEIN_DECL(NAMES, NAMEL, TYPE, TYPE_REAL) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const int* n, const TYPE_REAL* d, const TYPE_REAL* e, const int* m, \
                                 TYPE_REAL* w, const int* iblock, const int* isplit, const TYPE_REAL* orfac, \
                                 TYPE* Z, const int* iZ, const int* jZ, const int* descZ, \
                                 TYPE* work, const int* lwork, int* iwork, const int* liwork, \
                                 int* ifail, int* iclustr, TYPE_REAL* gap, int* info);

SCALAPACK_PSTEIN_DECL(psstein, PSSTEIN, float, float);
SCALAPACK_PSTEIN_DECL(pdstein, PDSTEIN, double, double);
SCALAPACK_PSTEIN_DECL(pcstein, PCSTEIN, lapack_complex_float, float);
SCALAPACK_PSTEIN_DECL(pzstein, PZSTEIN, lapack_complex_double, double);

#undef SCALAPACK_PSTEIN_DECL

#define SCALAPACK_PSTEDC_DECL(NAMES, NAMEL, TYPE) \
void ROKKO_GLOBAL(NAMES, NAMEL) (const char* compz, const int* n, TYPE* d, TYPE* e, \
                                 TYPE* Q, const int* iq, const int* jq, const int* descQ, \
                                 TYPE* work, const int* lwork, int* iwork, const int* liwork, int* info);

SCALAPACK_PSTEDC_DECL(psstedc, PSSTEDC, float);
SCALAPACK_PSTEDC_DECL(pdstedc, PDSTEDC, double);

#undef SCALAPACK_PSTEDC_DECL

#ifdef __cplusplus
}
#endif

#endif // ROKKO_SCALAPACK_SCALAPACK_INTERFACE_H

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko