Polylib_test.cpp File Reference

#include "Polylib.h"
#include <cmath>
#include <cstdio>

Go to the source code of this file.

Macros
#define	NPLOWER   5
#define	NPUPPER   15
#define	EPS   1e-12
#define	GAUSS_INT   1
#define	GAUSS_RADAUM_INT   1
#define	GAUSS_RADAUP_INT   1
#define	GAUSS_LOBATTO_INT   1
#define	GAUSS_DIFF   1
#define	GAUSS_RADAUM_DIFF   1
#define	GAUSS_RADAUP_DIFF   1
#define	GAUSS_LOBATTO_DIFF   1
#define	GAUSS_INTERP   1
#define	GAUSS_RADAUM_INTERP   1
#define	GAUSS_RADAUP_INTERP   1
#define	GAUSS_LOBATTO_INTERP   1

Functions
double	ddot (int, double *, int, double *, int)
double *	dvector (int, int)
int	test_gamma_fraction ()
int	main ()

Macro Definition Documentation

EPS

#define EPS   1e-12

Definition at line 96 of file Polylib_test.cpp.

GAUSS_DIFF

#define GAUSS_DIFF   1

Definition at line 102 of file Polylib_test.cpp.

GAUSS_INT

#define GAUSS_INT   1

Definition at line 98 of file Polylib_test.cpp.

GAUSS_INTERP

#define GAUSS_INTERP   1

Definition at line 106 of file Polylib_test.cpp.

GAUSS_LOBATTO_DIFF

#define GAUSS_LOBATTO_DIFF   1

Definition at line 105 of file Polylib_test.cpp.

GAUSS_LOBATTO_INT

#define GAUSS_LOBATTO_INT   1

Definition at line 101 of file Polylib_test.cpp.

GAUSS_LOBATTO_INTERP

#define GAUSS_LOBATTO_INTERP   1

Definition at line 109 of file Polylib_test.cpp.

GAUSS_RADAUM_DIFF

#define GAUSS_RADAUM_DIFF   1

Definition at line 103 of file Polylib_test.cpp.

GAUSS_RADAUM_INT

#define GAUSS_RADAUM_INT   1

Definition at line 99 of file Polylib_test.cpp.

GAUSS_RADAUM_INTERP

#define GAUSS_RADAUM_INTERP   1

Definition at line 107 of file Polylib_test.cpp.

GAUSS_RADAUP_DIFF

#define GAUSS_RADAUP_DIFF   1

Definition at line 104 of file Polylib_test.cpp.

GAUSS_RADAUP_INT

#define GAUSS_RADAUP_INT   1

Definition at line 100 of file Polylib_test.cpp.

GAUSS_RADAUP_INTERP

#define GAUSS_RADAUP_INTERP   1

Definition at line 108 of file Polylib_test.cpp.

NPLOWER

#define NPLOWER   5

Definition at line 94 of file Polylib_test.cpp.

NPUPPER

#define NPUPPER   15

Definition at line 95 of file Polylib_test.cpp.

Function Documentation

ddot()

double ddot	(	int	n,
		double *	x,
		int	incx,
		double *	y,
		int	incy
	)

Definition at line 564 of file Polylib_test.cpp.

{
    register double sum = 0.;
 
    while (n--)
    {
        sum += (*x) * (*y);
        x += incx;
        y += incy;
    }
    return sum;
}

Referenced by main().

dvector()

double * dvector	(	int	nl,
		int	nh
	)

Definition at line 577 of file Polylib_test.cpp.

{
    double *v;
 
    v = (double *)malloc((unsigned)(nh - nl + 1) * sizeof(double));
    return v - nl;
}

Referenced by main().

main()

int main

(

)

Definition at line 116 of file Polylib_test.cpp.

{
    int np, n, i;
    double *z, *w, *p, sum = 0, alpha, beta, *d, *dt;
 
    z = dvector(0, NPUPPER - 1);
    w = dvector(0, NPUPPER - 1);
    p = dvector(0, NPUPPER - 1);
 
    d  = dvector(0, NPUPPER * NPUPPER - 1);
    dt = dvector(0, NPUPPER * NPUPPER - 1);
 
#if GAUSS_INT
    /* Gauss Integration */
    printf("Begin checking Gauss integration\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgj(z, w, np, alpha, beta);
                for (n = 2; n < 2 * np - 1; ++n)
                {
                    jacobfd(np, z, p, NULL, n, alpha, beta);
                    sum = ddot(np, w, 1, p, 1);
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "integal was %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
 
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Integration\n");
#endif
 
#if GAUSS_RADAUM_INT
    /* Gauss Radau Integration */
    printf("Begin checking Gauss Radau Integration\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjm(z, w, np, alpha, beta);
                for (n = 2; n < 2 * np - 2; ++n)
                {
                    jacobfd(np, z, p, NULL, n, alpha, beta);
                    sum = ddot(np, w, 1, p, 1);
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "integal was %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
 
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radau (z=-1) Integration\n");
#endif
 
#if GAUSS_RADAUP_INT
    /* Gauss Radau Integration */
    printf("Begin checking Gauss Radau Integration\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjp(z, w, np, alpha, beta);
                for (n = 2; n < 2 * np - 2; ++n)
                {
                    jacobfd(np, z, p, NULL, n, alpha, beta);
                    sum = ddot(np, w, 1, p, 1);
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "integal was %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
 
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radau (z=1) Integration\n");
#endif
 
#if GAUSS_LOBATTO_INT
    /* Gauss Lobatto Integration */
    printf("Begin checking Gauss Lobatto integration\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwglj(z, w, np, alpha, beta);
                for (n = 2; n < 2 * np - 3; ++n)
                {
                    jacobfd(np, z, p, NULL, n, alpha, beta);
                    sum = ddot(np, w, 1, p, 1);
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "integal was %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
 
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Lobatto Integration\n");
#endif
 
#if GAUSS_DIFF
    printf("Begin checking differentiation through Gauss points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgj(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    Dgj(d, z, np, alpha, beta);
 
                    for (i = 0; i < np; ++i)
                        p[i] = pow(z[i], n);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum += fabs(ddot(np, d + i * np, 1, p, 1) -
                                    n * pow(z[i], n - 1));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Jacobi differentiation\n");
#endif
 
#if GAUSS_RADAUM_DIFF
    printf("Begin checking differentiation through Gauss Radau points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjm(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    Dgrjm(d, z, np, alpha, beta);
                    for (i = 0; i < np; ++i)
                        p[i] = pow(z[i], n);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum += fabs(ddot(np, d + i * np, 1, p, 1) -
                                    n * pow(z[i], n - 1));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radau (z=-1) differentiation\n");
#endif
 
#if GAUSS_RADAUP_DIFF
    printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjp(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    Dgrjp(d, z, np, alpha, beta);
                    for (i = 0; i < np; ++i)
                        p[i] = pow(z[i], n);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum += fabs(ddot(np, d + i * np, 1, p, 1) -
                                    n * pow(z[i], n - 1));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radau (z=1) differentiation\n");
#endif
 
#if GAUSS_LOBATTO_DIFF
    printf("Begin checking differentiation through Gauss Lobatto points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwglj(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    Dglj(d, z, np, alpha, beta);
                    for (i = 0; i < np; ++i)
                        p[i] = pow(z[i], n);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum += fabs(ddot(np, d + i * np, 1, p, 1) -
                                    n * pow(z[i], n - 1));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Lobatto differentiation\n");
#endif
 
    /* check interpolation routines */
#if GAUSS_INTERP
    printf("Begin checking interpolation through Gauss points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgj(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    for (i = 0; i < np; ++i)
                    {
                        w[i] = 2.0 * i / (double)(np - 1) - 1.0;
                        p[i] = pow(z[i], n);
                    }
                    Imgj(d, z, w, np, np, alpha, beta);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum +=
                            fabs(ddot(np, d + i * np, 1, p, 1) - pow(w[i], n));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Jacobi interpolation\n");
#endif
 
#if GAUSS_RADAUM_INTERP
    printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjm(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    for (i = 0; i < np; ++i)
                    {
                        w[i] = 2.0 * i / (double)(np - 1) - 1.0;
                        p[i] = pow(z[i], n);
                    }
                    Imgrjm(d, z, w, np, np, alpha, beta);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum +=
                            fabs(ddot(np, d + i * np, 1, p, 1) - pow(w[i], n));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
#endif
#if GAUSS_RADAUP_INTERP
    printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwgrjp(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    for (i = 0; i < np; ++i)
                    {
                        w[i] = 2.0 * i / (double)(np - 1) - 1.0;
                        p[i] = pow(z[i], n);
                    }
                    Imgrjp(d, z, w, np, np, alpha, beta);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum +=
                            fabs(ddot(np, d + i * np, 1, p, 1) - pow(w[i], n));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Radau (z=1) interpolation\n");
#endif
 
#if GAUSS_LOBATTO_INTERP
    printf("Begin checking Interpolation through Gauss Lobatto points\n");
    alpha = -0.5;
    while (alpha <= 5.0)
    {
        beta = -0.5;
        while (beta <= 5.0)
        {
 
            for (np = NPLOWER; np <= NPUPPER; ++np)
            {
                zwglj(z, w, np, alpha, beta);
                for (n = 2; n < np - 1; ++n)
                {
                    for (i = 0; i < np; ++i)
                    {
                        w[i] = 2.0 * i / (double)(np - 1) - 1.0;
                        p[i] = pow(z[i], n);
                    }
                    Imglj(d, z, w, np, np, alpha, beta);
                    sum = 0;
                    for (i = 0; i < np; ++i)
                        sum +=
                            fabs(ddot(np, d + i * np, 1, p, 1) - pow(w[i], n));
                    sum /= np;
                    if (fabs(sum) > EPS)
                        printf("alpha = %lf, beta = %lf, np = %d, n = %d "
                               "difference %lg\n",
                               alpha, beta, np, n, sum);
                }
            }
            beta += 0.5;
        }
        printf("finished checking all beta values for alpha = %lf\n", alpha);
        alpha += 0.5;
    }
    printf("Finished checking Gauss Lobatto interploation\n");
#endif
    test_gamma_fraction();
 
    free(z);
    free(w);
    free(p);
    free(d);
    free(dt);
    return 0;
}

References Nektar::LibUtilities::beta, ddot(), Polylib::Dgj(), Polylib::Dglj(), Polylib::Dgrjm(), Polylib::Dgrjp(), dvector(), EPS, Polylib::Imgj(), Polylib::Imglj(), Polylib::Imgrjm(), Polylib::Imgrjp(), Polylib::jacobfd(), NPLOWER, NPUPPER, CellMLToNektar.cellml_metadata::p, test_gamma_fraction(), Polylib::zwgj(), Polylib::zwglj(), Polylib::zwgrjm(), and Polylib::zwgrjp().

test_gamma_fraction()

int test_gamma_fraction

(

)

Definition at line 585 of file Polylib_test.cpp.

{
    double a = 362880.;
    double c = Polylib::gammaFracGammaF(10, 0., 0, 1.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10., 1.,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(1, 0., 12, -2.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 1., 10.,
           c * a - 1.);
 
    a = 30.;
    c = Polylib::gammaFracGammaF(4, 3., 5, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5., 7.,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(5, 0., 7, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 7., 5.,
           c * a - 1.);
 
    a = 21651.09375;
    c = Polylib::gammaFracGammaF(7, 3.5, 5, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 5.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(5, 0.5, 10, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5, 10.5,
           c * a - 1.);
 
    a = 97429.921875;
    c = Polylib::gammaFracGammaF(10, 0.5, 5, -0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 4.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(5, -0.5, 11, -0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 4.5, 10.5,
           c * a - 1.);
 
    a = 2279.0625;
    c = Polylib::gammaFracGammaF(10, -0.5, 5, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 9.5, 5.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(5, 0.5, 10, -0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 5.5, 9.5,
           c * a - 1.);
 
    a = 639383.8623046875;
    c = Polylib::gammaFracGammaF(10, 0.5, 0, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 10.5, 0.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(0, 0.5, 10, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 0.5, 10.5,
           c * a - 1.);
 
    a = 5967498288235982848.;
    c = Polylib::gammaFracGammaF(100, 0., 90, 0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 100., 90.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(90, 0.5, 100, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 90.5, 100.,
           c * a - 1.);
 
    a = 5618641603777298169856.;
    c = Polylib::gammaFracGammaF(200, 0., 191, -0.5);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200., 190.5,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(190, 0.5, 200, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190.5, 200.,
           c * a - 1.);
 
    a = 77396694214720029196288.;
    c = Polylib::gammaFracGammaF(200, 0., 190, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 200., 190.,
           c / a - 1.);
    c = Polylib::gammaFracGammaF(190, 0., 200, 0.);
    printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n", 190., 200.,
           c * a - 1.);
 
    int Q = 2;
    for (double alpha = -0.5; alpha <= 150.; alpha += 0.5)
    {
        for (double beta = -0.5; beta <= 150.; beta += 0.5)
        {
            a = Polylib::gammaF(Q + alpha) / Polylib::gammaF(Q + beta);
            c = Polylib::gammaFracGammaF(Q, alpha, Q, beta) / a - 1.;
            if (fabs(c) > 5.e-15)
            {
                printf("alpha = %7.2lf, beta = %7.2lf, difference %9.2e\n",
                       Q + alpha, Q + beta, c);
            }
        }
    }
 
    return 0;
}

References Nektar::LibUtilities::beta, Polylib::gammaF(), and Polylib::gammaFracGammaF().

Referenced by main().