test.c

开放gsl矩阵运算

    int order[5] = { 2, 1, 3, 5, 4 } ;

    gsl_function f = make_function(&myfn1, 0);
    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qagi (&fc, 1.0e-7, 0.0, w->limit,
				   w, 
				   &result, &abserr) ;
    
    gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn1) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn1) smooth abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn1) smooth neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qagiu(myfn1) smooth last") ;  
    gsl_test_int(status,exp_ier,"qagiu(myfn1) smooth status") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn1) smooth alist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn1) smooth blist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn1) smooth rlist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn1) smooth elist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn1) smooth order");

    gsl_integration_workspace_free (w) ;

  }

  /* Test infinite range integral myfn1 using an absolute error bound */

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

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    /* All results are for GSL_IEEE_MODE=double-precision */

    double exp_result = 2.718281828459044647E+00;
    double exp_abserr = 1.588185109253204805E-10;
    int exp_neval  =      135;
    int exp_ier    =        0;
    int exp_last   =        5;

    double a[5] = { 0.000000000000000000E+00,
		    5.000000000000000000E-01,
		    2.500000000000000000E-01,
		    1.250000000000000000E-01,
		    6.250000000000000000E-02 } ;
    double b[5] = { 6.250000000000000000E-02,
		    1.000000000000000000E+00,
		    5.000000000000000000E-01,
		    2.500000000000000000E-01,
		    1.250000000000000000E-01 } ;
    double r[5] = { 8.315287189746029816E-07,
		    1.718281828459045091E+00,
		    8.646647167633871867E-01,
		    1.328565310599463256E-01,
		    2.477920647947255521E-03 } ;
    double e[5] = { 1.533437090413525935E-10,
		    4.117868247943567505E-12,
		    7.802455785301941044E-13,
		    5.395586026138397182E-13,
		    3.713312434866150125E-14 } ;
    int order[5] = { 1, 2, 3, 4, 5 } ;

    double alpha = 1.0 ;
    gsl_function f = make_function(&myfn2, &alpha);
    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qagil (&fc, 1.0, 1.0e-7, 0.0, w->limit,
				    w, 
				    &result, &abserr) ;
    
    gsl_test_rel(result,exp_result,1e-14,"qagiu(myfn2) smooth result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qagiu(myfn2) smooth abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qagiu(myfn2) smooth neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qagiu(myfn2) smooth last") ;  
    gsl_test_int(status,exp_ier,"qagiu(myfn2) smooth status") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->alist[i],a[i],1e-15,"qagiu(myfn2) smooth alist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->blist[i],b[i],1e-15,"qagiu(myfn2) smooth blist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->rlist[i],r[i],1e-14,"qagiu(myfn2) smooth rlist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_rel(w->elist[i],e[i],1e-4,"qagiu(myfn2) smooth elist") ;

    for (i = 0; i < 5 ; i++) 
	gsl_test_int((int)w->order[i],order[i]-1,"qagiu(myfn2) smooth order");

    gsl_integration_workspace_free (w) ;

  }

  /* Test integral f454 with integrable singular points */

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

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    /* All results are for GSL_IEEE_MODE=double-precision */

    double exp_result = 5.274080611672716401E+01;
    double exp_abserr = 1.755703848687062418E-04;
    int exp_neval  =        777;
    int exp_ier    =          0;
    int exp_last   =         20;

    double a[20] = { 9.687500000000000000E-01,
		    1.401269388548935790E+00,
		    1.414213562373095145E+00,
		    1.000000000000000000E+00,
		    0.000000000000000000E+00,
		    2.207106781186547462E+00,
		    1.810660171779821415E+00,
		    1.207106781186547462E+00,
		    5.000000000000000000E-01,
		    1.103553390593273731E+00,
		    1.612436867076458391E+00,
		    1.310660171779821415E+00,
		    7.500000000000000000E-01,
		    1.051776695296636976E+00,
		    1.513325214724776657E+00,
		    1.362436867076458391E+00,
		    8.750000000000000000E-01,
		    1.463769388548935790E+00,
		    1.388325214724776657E+00,
		    9.375000000000000000E-01} ;
    double b[20] = { 1.000000000000000000E+00,
		     1.414213562373095145E+00,
		     1.463769388548935790E+00,
		     1.051776695296636976E+00,
		     5.000000000000000000E-01,
		     3.000000000000000000E+00,
		     2.207106781186547462E+00,
		     1.310660171779821415E+00,
		     7.500000000000000000E-01,
		     1.207106781186547462E+00,
		     1.810660171779821415E+00,
		     1.362436867076458391E+00,
		     8.750000000000000000E-01,
		     1.103553390593273731E+00,
		     1.612436867076458391E+00,
		     1.388325214724776657E+00,
		     9.375000000000000000E-01,
		     1.513325214724776657E+00,
		     1.401269388548935790E+00,
		     9.687500000000000000E-01} ;
    double r[20] = { -1.125078814079027711E-01,
		     -1.565132123531515207E-01,
		     -4.225328513207429193E-01,
		     -1.830392049835374568E-01,
		     6.575875041899758092E-03,
		     4.873920540843067783E+01,
		     6.032891565603589079E+00,
		     -2.991531901645863023E-01,
		     -7.326282608704996063E-03,
		     -2.431894410706912923E-01,
		     5.911661670635662835E-01,
		     -2.236786562536174916E-01,
		     -5.647871991778510847E-02,
		     -1.305470403178642658E-01,
		     -1.721363984401322045E-01,
		     -1.589345454585119055E-01,
		     -7.406626263352669715E-02,
		     -2.208730668000830344E-01,
		     -1.048692749517999567E-01,
		     -6.302287584527696551E-02} ;
    double e[20] = { 2.506431410088378817E-02,
		     2.730454695485963826E-02,
		     1.017446081816190118E-01,
		     3.252808038935910834E-02,
		     7.300687878575027348E-17,
		     5.411138804637469780E-13,
		     6.697855121200013106E-14,
		     3.321267596107916554E-15,
		     1.417509685426979386E-16,
		     2.699945168224041491E-15,
		     6.573952690524728748E-15,
		     2.483331942899818875E-15,
		     6.270397525408045936E-16,
		     1.449363299575615261E-15,
		     1.911097929242846383E-15,
		     1.764527917763735212E-15,
		     8.223007012367522077E-16,
		     2.452183642810224359E-15,
		     1.164282836272345215E-15,
		     6.996944784151910810E-16} ;
    int order[20] = { 3, 4, 2, 1, 6, 7, 11, 8, 10, 12, 18,
		     15, 16, 14, 19, 17, 20, 13, 9, 5 } ;

    gsl_function f = make_function(&f454, 0);
    gsl_function fc = make_counter(&f, &p) ;

    double pts[4] ;

    pts[0] = 0.0;
    pts[1] = 1.0;
    pts[2] = sqrt(2.0);
    pts[3] = 3.0;

    status = gsl_integration_qagp (&fc, pts, 4,
				   0.0, 1.0e-3, w->limit,
				   w, 
				   &result, &abserr) ;
    
    gsl_test_rel(result,exp_result,1e-14,"qagp(f454) singular result") ;
    gsl_test_rel(abserr,exp_abserr,1e-5,"qagp(f454) singular abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qagp(f454) singular neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qagp(f454) singular last") ;  
    gsl_test_int(status,exp_ier,"qagp(f454) singular status") ;

    for (i = 0; i < 20 ; i++) 
	gsl_test_rel(w->alist[i],a[i],1e-15,"qagp(f454) singular alist") ;

    for (i = 0; i < 20 ; i++) 
	gsl_test_rel(w->blist[i],b[i],1e-15,"qagp(f454) singular blist") ;

    for (i = 0; i < 20 ; i++) 
	gsl_test_rel(w->rlist[i],r[i],1e-14,"qagp(f454) singular rlist") ;

    for (i = 0; i < 20 ; i++) 
	gsl_test_rel(w->elist[i],e[i],1e-4,"qagp(f454) singular elist") ;

    for (i = 0; i < 20 ; i++) 
	gsl_test_int((int)w->order[i],order[i]-1,"qagp(f454) singular order");

    gsl_integration_workspace_free (w) ;

  }


  /* Test cauchy integration using a relative error bound */

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

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    /* All results are for GSL_IEEE_MODE=double-precision */

    double exp_result = -8.994400695837000137E-02;
    double exp_abserr =  1.185290176227023727E-06;
    int exp_neval  =      215;
    int exp_ier    =        0;
    int exp_last   =        6;

    double a[6] = { -1.000000000000000000E+00,
		    2.500000000000000000E+00,
		    1.250000000000000000E+00,
		    6.250000000000000000E-01,
		    -5.000000000000000000E-01,
		    -7.500000000000000000E-01} ;
    double b[6] = { -7.500000000000000000E-01,
		    5.000000000000000000E+00,
		    2.500000000000000000E+00,
		    1.250000000000000000E+00,
		    6.250000000000000000E-01,
		    -5.000000000000000000E-01} ;
    double r[6] = { -1.234231128040012976E-01,
		    3.579970394639702888E-03,
		    2.249831615049339983E-02,
		    7.214232992127905808E-02,
		    2.079093855884046535E-02,
		    -8.553244917962132821E-02} ;
    double e[6] = { 1.172832717970022565E-06,
		    9.018232896137375412E-13,
		    1.815172652101790755E-12,
		    1.006998195150956048E-13,
		    1.245463873006391609E-08,
		    1.833082948207153514E-15 } ;
    int order[6] = { 1, 5, 3, 2, 4, 6 } ;

    double alpha = 1.0 ;
    gsl_function f = make_function(&f459, &alpha);
    gsl_function fc = make_counter(&f, &p) ;

    status = gsl_integration_qawc (&fc, -1.0, 5.0, 0.0, 0.0, 1.0e-3, w->limit,
				   w, 
				   &result, &abserr) ;
    
    gsl_test_rel(result,exp_result,1e-14,"qawc(f459) result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qawc(f459) last") ;  
    gsl_test_int(status,exp_ier,"qawc(f459) status") ;

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

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

    for (i = 0; i < 6 ; i++) 
	gsl_test_rel(w->rlist[i],r[i],1e-14,"qawc(f459) rlist") ;

    for (i = 0; i < 6 ; i++) 
	gsl_test_rel(w->elist[i],e[i],1e-5,"qawc(f459) elist") ;

    for (i = 0; i < 6 ; i++) 
	gsl_test_int((int)w->order[i],order[i]-1,"qawc(f459) order");

    p.neval = 0;
    status = gsl_integration_qawc (&fc, 5.0, -1.0, 0.0, 0.0, 1.0e-3, w->limit,
				   w, 
				   &result, &abserr) ;
    
    gsl_test_rel(result,-exp_result,1e-14,"qawc(f459) rev result") ;
    gsl_test_rel(abserr,exp_abserr,1e-6,"qawc(f459) rev abserr") ;
    gsl_test_int((int)(p.neval),exp_neval,"qawc(f459) rev neval") ;  
    gsl_test_int((int)(w->size),exp_last,"qawc(f459) rev last") ;  
    gsl_test_int(status,exp_ier,"qawc(f459) rev status") ;

    gsl_integration_workspace_free (w) ;

  }

  /* Test QAWS singular integration using a relative error bound */

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

    gsl_integration_qaws_table * t 
      = gsl_integration_qaws_table_alloc (0.0, 0.0, 1, 0);

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000) ;

    /* All results are for GSL_IEEE_MODE=double-precision */

    double exp_result = -1.892751853489401670E-01;
    double exp_abserr = 1.129133712015747658E-08;
    int exp_neval  =      280;
    int exp_ier    =        0;
    int exp_last   =        8;

    double a[8] = { 0.000000000000000000E+00,
		    5.000000000000000000E-01,
		    2.500000000000000000E-01,
		    1.250000000000000000E-01,
		    6.250000000000000000E-02,
		    3.125000000000000000E-02,
		    1.562500000000000000E-02,
		    7.812500000000000000E-03} ;
    double b[8] = { 7.812500000000000000E