test.c

开放gsl矩阵运算

/* integration/test.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>

#include "tests.h"

gsl_function make_function (double (* f) (double, void *), double * p);

gsl_function make_function (double (* f) (double, void *), double * p)
{
  gsl_function f_new;

  f_new.function = f ;
  f_new.params = p ;

  return f_new;
}

struct counter_params {
  gsl_function * f;
  int neval;
} ;

double counter (double x, void * params);
gsl_function make_counter (gsl_function * f, struct counter_params * p);

double 
counter (double x, void * params)
{
  struct counter_params * p = (struct counter_params *) params;
  p->neval++ ; /* increment counter */
  return GSL_FN_EVAL(p->f, x);
}

gsl_function make_counter (gsl_function * f, struct counter_params * p)
{
  gsl_function f_new;

  p->f = f;
  p->neval = 0 ;
  
  f_new.function = &counter ;
  f_new.params = p ;

  return f_new;
}

void my_error_handler (const char *reason, const char *file,
		       int line, int err);

int main (void)
{
  gsl_ieee_env_setup ();
  gsl_set_error_handler (&my_error_handler); 

  /* Test the basic Gauss-Kronrod rules with a smooth positive function. */

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049357767090777E-02;
    double exp_abserr = 2.990224871000550874E-06;
    double exp_resabs = 7.716049357767090777E-02;
    double exp_resasc = 4.434273814139995384E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha) ;

    gsl_integration_qk15 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk15(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) smooth resasc") ;

    gsl_integration_qk15 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;

    gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049379303084599E-02;
    double exp_abserr = 9.424302194248481445E-08;
    double exp_resabs = 7.716049379303084599E-02;
    double exp_resasc = 4.434311425038358484E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk21 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk21(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) smooth resasc") ;

    gsl_integration_qk21 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049382494900855E-02;
    double exp_abserr = 1.713503193600029893E-09;
    double exp_resabs = 7.716049382494900855E-02;
    double exp_resasc = 4.427995051868838933E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk31 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk31(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) smooth resasc") ;

    gsl_integration_qk31 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049382681375302E-02;
    double exp_abserr = 9.576386660975511224E-11;
    double exp_resabs = 7.716049382681375302E-02;
    double exp_resasc = 4.421521169637691873E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk41 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk41(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) smooth resasc") ;

    gsl_integration_qk41 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049382708510540E-02;
    double exp_abserr = 1.002079980317363772E-11;
    double exp_resabs = 7.716049382708510540E-02;
    double exp_resasc = 4.416474291216854892E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk51 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk51(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) smooth resasc") ;

    gsl_integration_qk51 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk51(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 7.716049382713800753E-02;
    double exp_abserr = 1.566060362296155616E-12;
    double exp_resabs = 7.716049382713800753E-02;
    double exp_resasc = 4.419287685934316506E-02;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk61 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk61(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) smooth abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) smooth resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) smooth resasc") ;

    gsl_integration_qk61 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk61(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk61(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) reverse resasc") ;
  }

  /* Now test the basic rules with a positive function that has a
     singularity. This should give large values of abserr which would
     find discrepancies in the abserr calculation. */

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 1.555688196612745777E+01;
    double exp_abserr = 2.350164577239293706E+01;
    double exp_resabs = 1.555688196612745777E+01;
    double exp_resasc = 2.350164577239293706E+01;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk15 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk15(f1) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) singular resasc") ;

    gsl_integration_qk15 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk15(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 1.799045317938126232E+01;
    double exp_abserr = 2.782360287710622515E+01;
    double exp_resabs = 1.799045317938126232E+01;
    double exp_resasc = 2.782360287710622515E+01;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk21 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk21(f1) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) singular resasc") ;

    gsl_integration_qk21 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk21(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk21(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 2.081873305159121657E+01;
    double exp_abserr = 3.296500137482590276E+01;
    double exp_resabs = 2.081873305159121301E+01;
    double exp_resasc = 3.296500137482590276E+01;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk31 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk31(f1) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) singular resasc") ;

    gsl_integration_qk31 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk31(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 2.288677623903126701E+01;
    double exp_abserr = 3.671538820274916048E+01;
    double exp_resabs = 2.288677623903126701E+01;
    double exp_resasc = 3.671538820274916048E+01;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk41 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk41(f1) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) singular resasc") ;

    gsl_integration_qk41 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk41(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f1) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result = 2.449953612016972215E+01;
    double exp_abserr = 3.967771249391228849E+01;
    double exp_resabs = 2.449953612016972215E+01;
    double exp_resasc = 3.967771249391228849E+01;

    double alpha = -0.9 ;
    gsl_function f = make_function(&f1, &alpha);

    gsl_integration_qk51 (&f, 0.0, 1.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk51(f1) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) singular resasc") ;

    gsl_integration_qk51 (&f, 1.0, 0.0, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk51(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f1) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f1) reverse resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f1) reverse resasc") ;
  }