test.c

开放gsl矩阵运算

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

    double alpha = -0.9 ;
    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) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f1) singular abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f1) singular resabs") ;    
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f1) singular 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-7,"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") ;
  }

  /* Test the basic Gauss-Kronrod rules with a smooth oscillating
     function, over an unsymmetric range. This should find any
     discrepancies in the abscissae. */

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575483799046E-01;
    double exp_abserr = 8.760080200939757174E-06;
    double exp_resabs = 1.165564172429140788E+00;
    double exp_resasc = 9.334560307787327371E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_integration_qk15 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk15(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk15(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk15(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk15(f3) oscill resasc") ;

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

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575482959717E-01;
    double exp_abserr = 7.999213141433641888E-11;
    double exp_resabs = 1.150829032708484023E+00;
    double exp_resasc = 9.297591249133687619E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);
    
    gsl_integration_qk21 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk21(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) oscill resasc") ;

    gsl_integration_qk21 (&f, 2.71, 0.3,
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk21(f3) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qk21(f3) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk21(f3) reverse resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk21(f3) reverse resasc") ;
  }

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575482959717E-01;
    double exp_abserr = 1.285805464427459261E-14;
    double exp_resabs = 1.158150602093290571E+00;
    double exp_resasc = 9.277828092501518853E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_integration_qk31 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk31(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk31(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk31(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk31(f3) oscill resasc") ;

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

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575482959717E-01;
    double exp_abserr = 1.286535726271015626E-14;
    double exp_resabs = 1.158808363486595328E+00;
    double exp_resasc = 9.264382258645686985E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_integration_qk41 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk41(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk41(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk41(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk41(f3) oscill resasc") ;

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

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575482961938E-01;
    double exp_abserr = 1.285290995039385778E-14;
    double exp_resabs = 1.157687209264406381E+00;
    double exp_resasc = 9.264666884071264263E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_integration_qk51 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk51(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk51(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk51(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk51(f3) oscill resasc") ;

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

  {
    double result = 0, abserr = 0, resabs = 0, resasc = 0 ;
    double exp_result =-7.238969575482959717E-01;
    double exp_abserr = 1.286438572027470736E-14;
    double exp_resabs = 1.158720854723590099E+00;
    double exp_resasc = 9.270469641771273972E-01;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    gsl_integration_qk61 (&f, 0.3, 2.71, 
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,exp_result,1e-15,"qk61(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) oscill abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) oscill resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) oscill resasc") ;

    gsl_integration_qk61 (&f, 2.71, 0.3,
				  &result, &abserr, &resabs, &resasc) ;
    gsl_test_rel(result,-exp_result,1e-15,"qk61(f3) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qk61(f3) reverse abserr") ;
    gsl_test_rel(resabs,exp_resabs,1e-15,"qk61(f3) reverse resabs") ;
    gsl_test_rel(resasc,exp_resasc,1e-15,"qk61(f3) reverse resasc") ;
  }

  /* Test the non-adaptive gaussian integrator QNG */

  {
    int status = 0; size_t neval = 0 ;
    double result = 0, abserr = 0 ;
    double exp_result = 7.716049379303083211E-02;
    double exp_abserr = 9.424302199601294244E-08;
    int exp_neval  =  21;
    int exp_ier    =   0;

    double alpha = 2.6 ;
    gsl_function f = make_function(&f1, &alpha);
    
    status = gsl_integration_qng (&f, 0.0, 1.0, 1e-1, 0.0,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) smooth abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) smooth neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) smooth status") ;

    status = gsl_integration_qng (&f, 1.0, 0.0, 1e-1, 0.0,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) reverse abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) reverse neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) reverse status") ;
  }

  {
    int status = 0; size_t neval = 0 ;
    double result = 0, abserr = 0 ;

    double exp_result = 7.716049382706505200E-02;
    double exp_abserr = 2.666893044866214501E-12;
    int exp_neval  =  43;
    int exp_ier    =   0;

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

    status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-9,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qng(f1) smooth 43pt result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qng(f1) smooth 43pt abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) smooth 43pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) smooth 43pt status") ;

    status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-9,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse 43pt result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qng(f1) reverse 43pt abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) reverse 43pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) reverse 43pt status") ;
  }

  {
    int status; size_t neval = 0 ;
    double result = 0, abserr = 0 ;
    double exp_result =-7.238969575482961938E-01;
    double exp_abserr = 1.277676889520056369E-14;
    int exp_neval  =  43;
    int exp_ier    =   0;

    double alpha = 1.3 ;
    gsl_function f = make_function(&f3, &alpha);

    status = gsl_integration_qng (&f, 0.3, 2.71, 0.0, 1e-12,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qnq(f3) oscill result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) oscill abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f3) oscill neval") ;
    gsl_test_int(status,exp_ier,"qng(f3) oscill status") ;

    status = gsl_integration_qng (&f, 2.71, 0.3, 0.0, 1e-12,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qnq(f3) reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f3) reverse abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f3) reverse neval") ;
    gsl_test_int(status,exp_ier,"qng(f3) reverse status") ;
  }

  {
    int status = 0; size_t neval = 0 ;
    double result = 0, abserr = 0 ;

    double exp_result = 7.716049382716029525E-02;
    double exp_abserr = 8.566535680046930668E-16;
    int exp_neval  =  87;
    int exp_ier    =   0;

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

    status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-13,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qng(f1) 87pt smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt smooth abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt smooth neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) 87pt smooth status") ;

    status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-13,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qng(f1) 87pt reverse result") ;
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) 87pt reverse abserr") ;
    gsl_test_int((int)neval,exp_neval,"qng(f1) 87pt reverse neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) 87pt reverse status") ;
  }

  {
    int status = 0; size_t neval = 0 ;
    double result = 0, abserr = 0 ;

    double exp_result = 3.222948711817264211E+01;
    double exp_abserr = 2.782360287710622870E+01;
    int exp_neval  =  87;
    int exp_ier    =  GSL_ETOL;

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

    status = gsl_integration_qng (&f, 0.0, 1.0, 0.0, 1e-3,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,exp_result,1e-15,"qng(f1) sing beyond 87pt result");
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) sing beyond 87pt abserr");
    gsl_test_int((int)neval,exp_neval,"qng(f1) sing beyond 87pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) sing beyond 87pt status") ;

    status = gsl_integration_qng (&f, 1.0, 0.0, 0.0, 1e-3,
				  &result, &abserr, &neval) ;
    gsl_test_rel(result,-exp_result,1e-15,"qng(f1) reverse beyond 87pt result");
    gsl_test_rel(abserr,exp_abserr,1e-7,"qng(f1) rev beyond 87pt abserr");
    gsl_test_int((int)neval,exp_neval,"qng(f1) rev beyond 87pt neval") ;  
    gsl_test_int(status,exp_ier,"qng(f1) rev beyond 87pt status") ;
  }

  /* Test the adaptive integrator QAG */

  {
    int status = 0, i; struct counter_params p;
    double result = 0, abserr=0;

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    double exp_result = 7.716049382715854665E-02 ;
    double exp_abserr = 6.679384885865053037E-12 ;
    int exp_neval  =     165;
    int exp_ier    =       0;
    int exp_last   =       6;

    double a[6] = { 0, 0.5, 0.25, 0.125, 0.0625, 0.03125 } ;
    double b[6] = { 0.03125, 1, 0.5, 0.25, 0.125, 0.0625 } ;
    double r[6] = { 3.966769831709074375E-06, 5.491842501998222409E-02,
		    1.909827770934243926E-02, 2.776531175604360531E-03,
		    3.280661030752063693E-04, 3.522704932261797744E-05 } ;
    double e[6] = { 6.678528276336181873E-12, 6.097169993333454062E-16,
		    2.120334764359736934E-16, 3.082568839745514608E-17,
		    3.642265412331439511E-18, 3.910988124757650942E-19 } ;
    int order[6] = { 1, 2, 3, 4, 5, 6 } ;

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

    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qag (&fc, 0.0, 1.0, 0.0, 1e-10, w->limit,
				  GSL_INTEG_GAUSS15, w,
				  &result, &abserr) ;

    gsl_test_rel(result,exp_result,1e-15,"qag(f1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qag(f1) smooth abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qag(f1) smooth neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qag(f1) smooth last") ;  
    gsl_test_int(status,exp_ier,"qag(f1) smooth status") ;

    for (i = 0; i < 6 ; i++) 
	gsl_test_rel(w->alist[i],a[i],1e-15,"qag(f1) smooth alist") ;

    for (i = 0; i < 6 ; i++) 
	gsl_test_rel(w->blist[i],b[i],1e-15,"qag(f1) smooth blist") ;

    for (i = 0; i < 6 ; i++)