erreur de compilation( cygwin)

Marsh Posté le 04-02-2011 à 11:09:20

Je suis un doctorant en finance et j’ai appris tout seul à programmer avec le langage C. Certes, je sais quelques notions dans la programmation mais je suis toujours débutant dans ce domaine. Je compile les programmes avec cygwin. De plus, j’utilise, lors de cette compilation, la librairie GSL (GNU scientific Librairy).
Bref, lors de l’exécution d’un programme (qui permet de calculer une intégration numérique et que j’ai nommé cdf_mdim), j’ai eu le message suivant (avec le débuggeur gdb) :

9 [main] cdf_Mdim 3816 _cygtls::handle_exceptions: Exception: STATUS_ACCESS_VIOLATION
1637 [main] cdf_Mdim 3816 open_stackdumpfile: Dumping stack trace to cdf_Mdim.exe.stackdump
521872 [main] cdf_Mdim 3816 _cygtls::handle_exceptions: Exception: STATUS_ACCESS_VIOLATION
578746 [main] cdf_Mdim 3816 _cygtls::handle_exceptions: Error while dumping state (probably corrupted stack)
Program received signal SIGSEGV, Segmentation fault.
0x61016525 in stack_info::walk () from /usr/bin/cygwin1.dll

Est que quelqu'un peut m’aider à résoudre ce problème (je peux vous envoyer mon code par mail dans ce cas). Est-ce que c’est un problème de mémoire de mon ordinateur (processeur centrino duo 1.73 ghz et une RAM de 2 Go) ?
Merci d’avance

Marsh Posté le 04-02-2011 à 11:09:20

Marsh Posté le 04-02-2011 à 14:04:07

C'est un problème de pile corrompu.

Cela peut venir, par exemple :

- un tableau qui est déclaré avec une taille trop petite par rapport à ce que l'on veut y mettre,
- une chaîne de caractères qui n'est pas terminée par un caractère nul,
- un indice qui est négatif (toto[-1]),
- une mauvaise libération des ressources (par exemple un oubli d'un free; par exemple un free() en trop (if (toto != NULL) free(toto);, mais toto n'a pas été initialisé, et contient autre chose que NULL par hasard); par exemple, un mélange entre delete du C++ et free du C; etc).

Marsh Posté le 04-02-2011 à 14:13:40

Voici mon programme principal:

Code :

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_sf_exp.h>
#include "cubature.h"
int count;
int n;
double bivar_cdf(unsigned dim, const double * x1, void *data)
{
int i,j;
int s;
double det,val, S;
double *R;
double P =dim;
R = (double *) malloc(1 * sizeof(double));
gsl_permutation* p = gsl_permutation_alloc(dim);
gsl_vector * x = gsl_vector_alloc (dim);
gsl_vector * xs = gsl_vector_alloc (dim);
gsl_matrix * A = gsl_matrix_alloc (dim, dim);
gsl_matrix * At = gsl_matrix_alloc (dim, dim);
gsl_matrix * SIG = gsl_matrix_alloc (dim, dim);
gsl_matrix *SIGi = gsl_matrix_alloc (dim, dim);
++count;
for (i=0; i < dim; i++){
for (j=0; j < dim; j++){
if (i < j) gsl_matrix_set (A,i,j,0);
else gsl_matrix_set (A,i,j,1);
}}
for (i=0; i < dim; i++) gsl_vector_set (x,i,x1[i]);
gsl_matrix_transpose_memcpy(At,A);
gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,1.0, A, At, 0.0, SIG);
gsl_linalg_LU_decomp( SIG, p, &s );
gsl_linalg_LU_invert(SIG,p,SIGi);
det = gsl_linalg_LU_det(SIG,s);
gsl_blas_dgemv(CblasNoTrans, 1.0, SIGi, x, 0.0, xs);
gsl_blas_ddot(x,xs,R);
S = exp ((-R[0]/2));
val =(1/(pow((2*M_PI),(P/2)))*(pow(abs(det),0.5)))*S;
return val;
}
int main()
{
int i;
double val, err;
double tol = 1.e-6;
double *xmin, *xmax;
unsigned dim = 6;
unsigned maxEval = 1000000;
printf("#==================================================================\n" );
printf("# Parameters:\n" );
printf("# dim = %d tol = %f maximum number of evals = %d\n",dim,tol,maxEval);
xmin = (double *) malloc(dim * sizeof(double));
xmax = (double *) malloc(dim * sizeof(double));
for (i = 0; i < dim; i++) {
xmin[i] = -99;
xmax[i] = 1.96;
}
count = 0;
adapt_integrate(bivar_cdf, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err);
printf("#------------------------------------------------------------------\n" );
printf("# Results:\n" );
printf("# Integration value = %g\n",val);
printf("# Est. err = %g\n", err);
printf("# Number of function evals required = %d\n", count);
printf("#==================================================================\n" );
free(xmax);
free(xmin);
return 0;
}

RQ: si maxEval = 100000 le programme s'exécute normalement.

Message édité par gilou le 04-02-2011 à 15:57:28

Marsh Posté le 04-02-2011 à 14:37:48

unsigned maxEval = 1000000;

Avec quelle valeur de maxEval est-ce que cela ne va plus ?

Par ailleurs, il serait plus lisible d'écrire unsigned int ou unsigned long int ou unsigned char.
Si c'est un projet pour DOS, il me semble que unsigned int est sur 16 bit donc que sa valeur maximale est 65535.

double *R;
...
R = (double *) malloc(1 * sizeof(double));
...
gsl_blas_ddot(x,xs,R);

R est alloué avec malloc, mais n'est pas libéré.

D'ailleurs, il serait plus facile et plus performant de remplacer cette allocation, par une simple déclaration.

double R;
...
gsl_blas_ddot(x,xs,&R);

Message édité par billgatesanonym le 04-02-2011 à 14:39:04

Marsh Posté le 04-02-2011 à 14:46:35

Juste un point de détail:

Citation :

Par ailleurs, il serait plus lisible d'écrire unsigned int ou unsigned long int ou unsigned char.
Si c'est un projet pour DOS, il me semble que unsigned int est sur 16 bit donc que sa valeur maximale est 65535.

En C, si ça ne tient pas dans un int, le compilateur doit automatiquement le mettre dans un long. Donc à priori, ça devrait tenir.
EDIT: C'était en K&R C, mais apparemment, l'ANSI C impose unsigned = unsigned int en ce cas.

Si le pb est avec la taille de maxEval, comme on ne sait pas comment la valeur est utilisée dans adapt_integrate(bivar_cdf, 0, dim, xmin, xmax, maxEval, 0, tol, &val, &err); on ne peut pas dire grand chose.

A+,

Message édité par gilou le 04-02-2011 à 16:25:09

---------------
There's more than what can be linked! -- Iyashikei Anime Forever! -- AngularJS c'est un framework d'engulé! --

Marsh Posté le 04-02-2011 à 15:27:04

voici la fonction adapt-integrate:

Code :

#include "cubature_copyright.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include <float.h>
#include <string.h>
void Usage();
typedef double (*integrand) (unsigned ndim, const double *x, void *);
int adapt_integrate(integrand f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
unsigned maxEval,
double reqAbsError, double reqRelError,
double *val, double *estimated_error);
typedef struct {
double val, err;
} esterr;
/*
* Function prototypes
*/
static double relError(esterr ee);
typedef struct {
unsigned dim;
double *center; /* length 2*dim = center followed by half-width */
double *width; /* length 2*dim = center followed by half-width */
double vol; /* cache volume = product of widths */
} hypercube;
typedef struct {

Message édité par gilou le 04-02-2011 à 15:57:45

Marsh Posté le 04-02-2011 à 15:32:39

C'est son proto, mais ça ne nous dit toujours pas ce qu'elle fait de maxEval
Par contre effectivement, quand tu passes ça en paramètre de fonction, le unsigned est probablement casté en unsigned int, et la ta valeur 100000 pourrait ne pas passer.
A+,

Message édité par gilou le 04-02-2011 à 15:34:15

---------------
There's more than what can be linked! -- Iyashikei Anime Forever! -- AngularJS c'est un framework d'engulé! --

Marsh Posté le 04-02-2011 à 15:56:26

OK gilou, je vais essayer
voici le programme cubature.c

Code :

#include "cubature.h"
#include "cubature_debug.h"
/*
Adaptive multidimensional integration on hypercubes (or, really,
hyper-rectangles) using cubature rules.
A cubature rule takes a function and a hypercube and evaluates
the function at a small number of points, returning an estimate
of the integral as well as an estimate of the error, and also
a suggested dimension of the hypercube to subdivide.
Given such a rule, the adaptive integration is simple:
1) Evaluate the cubature rule on the hypercube(s).
Stop if converged.
2) Pick the hypercube with the largest estimated error,
and divide it in two along the suggested dimension.
3) Goto (1).
*/
/*
* Integrate the function f from xmin[dim] to xmax[dim], with at most
* maxEval function evaluations (0 for no limit), until the given
* absolute is achieved relative error. val returns the integral, and
* estimated_error returns the estimate for the absolute error in val.
* The return value of the function is 0 on success and non-zero if there
* was an error.
*/
/*
*=======================================================
* cubature subroutines themselves follow.
*=======================================================
*/
static double relError(esterr ee)
{
return (ee.val == 0.0 ? HUGE_VAL : fabs(ee.err / ee.val));
}
/*
* Compute volume of hypercube from the product of its widths
*/
static double compute_vol(const hypercube * h)
{
unsigned i;
double vol = 1;
for (i = 0; i < h->dim; i++)
vol *= 2.0 * h->width[i];
return vol;
}
/*
* This should be more efficient than looping to copy contents and still
* requires a single malloc. Saves on indirect addressing arithmetic
* later, easier to understand the code.
*/
static hypercube make_hypercube(unsigned dim, const double *center,
const double *width)
{
unsigned long int i,dlen;
hypercube h;
h.dim = dim;
dlen = sizeof(double)*dim;
h.center = (double *) malloc(2*dlen);
h.width = &h.center[dim];
memcpy(h.center,center,dlen);
memcpy(h.width,width,dlen);
h.vol = compute_vol(&h);
if(debug(D_CUBATURE)){
printf("# make_hypercube(): h.center = (" );
for(i = 0;i<dim;i++){
printf("%4.2f",h.center[i]);
if(i != dim-1) printf("," );
}
printf(" ) h.width = (" );
for(i = 0;i<dim;i++){
printf("%4.2f",h.width[i]);
if(i != dim-1) printf("," );
}
printf(" )\n" );
}
return h;
}
static hypercube make_hypercube_range(unsigned dim, const double *xmin,
const double *xmax)
{
unsigned long int i,dlen;
hypercube h;
h.dim = dim;
dlen = sizeof(double)*dim;
h.center = (double *) malloc(2*dlen);
h.width = &h.center[dim];
for (i = 0; i < dim; ++i) {
h.center[i] = 0.5 * (xmin[i] + xmax[i]);
h.width[i] = 0.5 * (xmax[i] - xmin[i]);
}
h.vol = compute_vol(&h);
/* printf("Made ranged hypercube at %0x, %0x\n",h.center,h.hwidth); */
return h;
}
/*
* Destroy a hypercube. Note that the data was malloc'd only
* to h->center, so we don't have to do anything with h->hwidth.
*/
static void destroy_hypercube(hypercube *h)
{
/* printf("About to free hypercube %0x\n",h->center); */
free(h->center);
h->dim = 0;
}
static region make_region(const hypercube * h)
{
region R;
if(debug(D_CUBATURE)){
printf("# make_region(): Entry.\n" );
}
R.h = make_hypercube(h->dim, h->center, h->width);
R.splitDim = 0;
return R;
}
static void destroy_region(region * R)
{
destroy_hypercube(&R->h);
}
/*
* This takes a region with a splitDim defined and creates two regions
* (hypercubes) with half the width in the selected dimension. The original
* region R contains the left-half region, the region R2 contains the
* right hand region on return.
*/
static void cut_region(region * R, region * R2)
{
unsigned d = R->splitDim, dim = R->h.dim;
/* Halve the width and volume */
R->h.width[d] *= 0.5;
R->h.vol *= 0.5;
/* Make a new hypercube from the center to the new width. */
R2->h = make_hypercube(dim, R->h.center, R->h.width);
/* Shift the center of the original back to the left */
R->h.center[d] -= R->h.width[d];
/* Shift the center of the new one over to the right */
R2->h.center[d] += R->h.width[d];
}
static void destroy_rule(rule * r)
{
if (r->destroy)
r->destroy(r);
free(r);
}
/*
* This routine returns a region with its "worst direction" set.
* At this point the integral for the region is already set as well.
*/
static region eval_region(region R, integrand f, void *fdata, rule * r)
{
if(debug(D_CUBATURE)){
printf("# eval_region(): Entry.\n" );
}
R.splitDim = r->evalError(r, f, fdata, &R.h, &R.ee);
return R;
}
/*
*============================================================
* Functions to loop over points in a hypercube.
*============================================================
*/
/*
* Based on orbitrule.cpp in HIntLib-0.0.10
*/
/*
* ls0 returns the least-significant 0 bit of n (e.g. it returns 0 if the LSB
* is 0, it returns 1 if the 2 LSBs are 01, etcetera).
*/
#if (defined(__GNUC__) || defined(__ICC)) && (defined(__i386__) || defined (__x86_64__))
/*
* use x86 bit-scan instruction, based on count_trailing_zeros() macro in
* GNU GMP's longlong.h.
*/
static unsigned ls0(unsigned n)
{
unsigned count;
n = ~n;
__asm__("bsfl %1,%0": "=r"(count):"rm"(n));
return count;
}
#else
static unsigned ls0(unsigned n)
{
const unsigned bits[] = {
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 7,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 6,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 5,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 4,
0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 8,
};
unsigned bit;
bit = 0;
while ((n & 0xff) == 0xff) {
n >>= 8;
bit += 8;
}
return bit + bits[n & 0xff];
}
#endif
/*
* This rule is used for sum5/lambda5 displacements. Note that "r" is
* actually plus or minus the (half) width scaled by lambda(s).
*
* Evaluate the integral on all 2^n points (+/-r,...+/-r) * A Gray-code
* ordering is used to minimize the number of coordinate updates in p.
*/
static double evalR_Rfs(integrand f, void *fdata, unsigned dim, double *p,
const double *c, const double *r)
{
double sum = 0;
unsigned i, ip;
unsigned signs = 0; /* 0/1 bit = +/- for corresponding element of r[] */
if(debug(D_CUBATURE)){
printf("# evalR_Rfs(): Entry.\n" );
}
/*
* We start with the point where r is ADDed in every coordinate (This implies
* signs=0). The points[i] are center[i] + width[i].
*/
for (i = 0; i < dim; ++i)
p[i] = c[i] + r[i];
/* Loop through the points in gray-code ordering */
for (i = 0;; ++i) {
unsigned mask;
if(debug(D_CUBATURE)){
printf("# evalR_Rfs(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum += f(dim, p, fdata);
unsigned d = ls0(i); /* which coordinate to flip */
if (d >= dim)
break;
/* flip the d-th bit and add/subtract r[d] */
mask = 1U << d;
signs ^= mask;
p[d] = (signs & mask) ? c[d] - r[d] : c[d] + r[d];
}
if(debug(D_CUBATURE)){
printf("# evalR_Rfs(): Done.\n" );
}
return sum;
}
/*
* This one is used to evaluate sum4 at lambda4, which is really
* lambda. Go figure.
*/
static double evalRR0_0fs(integrand f, void *fdata, unsigned dim,
double *p, const double *c, const double *r)
{
unsigned i, j, ip;
double sum = 0;
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Entry.\n" );
}
for (i = 0; i < dim - 1; i++) {
p[i] = c[i] - r[i];
for (j = i + 1; j < dim; j++) {
p[j] = c[j] - r[j];
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum += f(dim, p, fdata);
p[i] = c[i] + r[i];
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum += f(dim, p, fdata);
p[j] = c[j] + r[j];
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum += f(dim, p, fdata);
p[i] = c[i] - r[i];
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum += f(dim, p, fdata);
p[j] = c[j]; /* Done with j -> Restore p[j] */
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Restoring p to (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
}
p[i] = c[i]; /* Done with i -> Restore p[i] */
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Restoring p to (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
}
if(debug(D_CUBATURE)){
printf("# evalRR0_0fs(): Done.\n" );
}
return sum;
}
/*
* This evaluates sum1, sum2 and sum3, mislabelled for no good
* reason (relative to their call) by one. They correspond as:
* center(p)->sum1(0), lambda2(Lambda2)->sum2(1), lambda4(Lambda)->sum3(2)
* which is just incredibly Evil. I'll really think hard about fixing
* this ONCE I've figured out why the routine is calling points that
* are WAY past the integration limits, which is almost certainly a simple
* bug.
*/
static double evalR0_0fs4d(integrand f, void *fdata, unsigned dim,
double *p, const double *c, double *sum0_,
const double *r1, double *sum1_,
const double *r2, double *sum2_)
{
double maxdiff = 0;
unsigned i, ip, dimDiffMax = 0;
double sum0, sum1 = 0, sum2 = 0; /* copies for aliasing, performance */
double ratio = r1[0] / r2[0];
ratio *= ratio;
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Entry.\n" );
}
/*
* Evaluate at the center of the current hypercube
*/
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",c[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum0 = f(dim, c, fdata);
for (i = 0; i < dim; i++) {
double f1a, f1b, f2a, f2b, diff;
p[i] = c[i] - r1[i];
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum1 += (f1a = f(dim, p, fdata));
p[i] = c[i] + r1[i];
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum1 += (f1b = f(dim, p, fdata));
p[i] = c[i] - r2[i];
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum2 += (f2a = f(dim, p, fdata));
p[i] = c[i] + r2[i];
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Evaluating f at (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
sum2 += (f2b = f(dim, p, fdata));
p[i] = c[i];
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Setting p back to (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",p[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
diff = fabs(f1a + f1b - 2 * sum0 - ratio * (f2a + f2b - 2 * sum0));
if (diff > maxdiff) {
maxdiff = diff;
dimDiffMax = i;
}
}
*sum0_ += sum0;
*sum1_ += sum1;
*sum2_ += sum2;
if(debug(D_CUBATURE)){
printf("# evalR0_0fs4d(): Done.\n" );
}
return dimDiffMax;
}
#define num0_0(dim) (1U)
#define numR0_0fs(dim) (2 * (dim))
#define numRR0_0fs(dim) (2 * (dim) * (dim-1))
#define numR_Rfs(dim) (1U << (dim))
/*
*=======================================================================
* Based on rule75genzmalik.cpp in HIntLib-0.0.10: An embedded cubature
* rule of degree 7 (embedded rule degree 5) due to A. C. Genz and A. A.
* Malik. See: A. C. Genz and A. A. Malik, "An imbedded [sic] family of
* fully symmetric numerical integration rules," SIAM J. Numer. Anal. 20
* (3), 580-588 (1983).
*=======================================================================
*/
#define real(x) ((double)(x))
#define to_int(n) ((int)(n))
static int isqr(int x)
{
return x * x;
}
static void destroy_rule75genzmalik(rule * r_)
{
rule75genzmalik *r = (rule75genzmalik *) r_;
free(r->p);
}
static unsigned rule75genzmalik_evalError(rule * r_, integrand f,
void *fdata, const hypercube * h,
esterr * ee)
{
/* lambda2 = sqrt(9/70), lambda4 = sqrt(9/10), lambda5 = sqrt(9/19) */
const double lambda2 = 0.3585685828003180919906451539079374954541;
const double lambda4 = 0.9486832980505137995996680633298155601160;
const double lambda5 = 0.6882472016116852977216287342936235251269;
const double weight2 = 980. / 6561.;
const double weight4 = 200. / 19683.;
const double weightE2 = 245. / 486.;
const double weightE4 = 25. / 729.;
rule75genzmalik *r = (rule75genzmalik *) r_;
unsigned i, dimDiffMax, dim = r_->dim;
double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4, sum5, result, res5th;
if(debug(D_CUBATURE)){
printf("# rule75genzmalik_evalError(): Entry.\n" );
}
for (i = 0; i < dim; i++){
r->p[i] = h->center[i];
r->widthLambda[i] = h->width[i] * lambda4;
r->widthLambda2[i] = h->width[i] * lambda2;
r->widthLambda5[i] = h->width[i] * lambda5;
if(debug(D_CUBATURE)){
printf("# i=%d: center = %f, wlam = %f, wlam2 = %f, wlam5 = %f\n",i,
r->p[i],r->widthLambda[i],r->widthLambda2[i],r->widthLambda5[i]);
}
}
/*
* Evaluate function at the center, at f(lambda2,0,...,0) and at
* f(lambda3=lambda4, 0,...,0). Estimate dimension with largest error
*/
dimDiffMax =
evalR0_0fs4d(f, fdata, dim, r->p, h->center, &sum1, r->widthLambda2,
&sum2, r->widthLambda, &sum3);
/*
* Calculate sum4 for f(lambda4, lambda4, 0, ...,0)
*/
sum4 = evalRR0_0fs(f, fdata, dim, r->p, h->center, r->widthLambda);
/*
* Calculate sum5 for f(lambda5, lambda5, ..., lambda5)
*/
sum5 = evalR_Rfs(f, fdata, dim, r->p, h->center, r->widthLambda5);
/* Calculate fifth and seventh order results */
result =
h->vol * (r->weight1 * sum1 + weight2 * sum2 + r->weight3 * sum3 +
weight4 * sum4 + r->weight5 * sum5);
res5th =
h->vol * (r->weightE1 * sum1 + weightE2 * sum2 +
r->weightE3 * sum3 + weightE4 * sum4);
ee->val = result;
ee->err = fabs(res5th - result);
return dimDiffMax;
}
static rule *make_rule75genzmalik(unsigned dim)
{
rule75genzmalik *r;
if (dim < 2)
return 0; /* this rule does not support 1d integrals */
r = (rule75genzmalik *) malloc(sizeof(rule75genzmalik));
r->parent.dim = dim;
r->weight1 =
(real(12824 - 9120 * to_int(dim) + 400 * isqr(to_int(dim)))
/ real(19683));
r->weight3 = real(1820 - 400 * to_int(dim)) / real(19683);
r->weight5 = real(6859) / real(19683) / real(1U << dim);
r->weightE1 = (real(729 - 950 * to_int(dim) + 50 * isqr(to_int(dim)))
/ real(729));
r->weightE3 = real(265 - 100 * to_int(dim)) / real(1458);
r->p = (double *) malloc(sizeof(double) * dim * 4);
r->widthLambda = r->p + dim;
r->widthLambda2 = r->p + 2 * dim;
r->widthLambda5 = r->p + 3 * dim;
r->parent.num_points = num0_0(dim) + 2 * numR0_0fs(dim)
+ numRR0_0fs(dim) + numR_Rfs(dim);
r->parent.evalError = rule75genzmalik_evalError;
r->parent.destroy = destroy_rule75genzmalik;
return (rule *) r;
}
/***************************************************************************/
/*
* 1d 15-point Gaussian quadrature rule, based on qk15.c and qk.c in GNU
* GSL (which in turn is based on QUADPACK).
*/
static unsigned rule15gauss_evalError(rule * r, integrand f, void *fdata,
const hypercube * h, esterr * ee)
{
/*
* Gauss quadrature weights and kronrod quadrature abscissae and weights as
* evaluated with 80 decimal digit arithmetic by L. W. Fullerton, Bell Labs,
* Nov. 1981.
*/
const unsigned n = 8;
const double xgk[8] = { /* abscissae of the 15-point kronrod rule */
0.991455371120812639206854697526329,
0.949107912342758524526189684047851,
0.864864423359769072789712788640926,
0.741531185599394439863864773280788,
0.586087235467691130294144838258730,
0.405845151377397166906606412076961,
0.207784955007898467600689403773245,
0.000000000000000000000000000000000
/* xgk[1], xgk[3], ... abscissae of the 7-point gauss rule.
xgk[0], xgk[2], ... to optimally extend the 7-point gauss rule */
};
static const double wg[4] = { /* weights of the 7-point gauss rule */
0.129484966168869693270611432679082,
0.279705391489276667901467771423780,
0.381830050505118944950369775488975,
0.417959183673469387755102040816327
};
static const double wgk[8] = { /* weights of the 15-point kronrod rule */
0.022935322010529224963732008058970,
0.063092092629978553290700663189204,
0.104790010322250183839876322541518,
0.140653259715525918745189590510238,
0.169004726639267902826583426598550,
0.190350578064785409913256402421014,
0.204432940075298892414161999234649,
0.209482141084727828012999174891714
};
const double center = h->center[0];
const double width = h->width[0];
double fv1[n - 1], fv2[n - 1];
const double f_center = f(1, &center, fdata);
double result_gauss = f_center * wg[n / 2 - 1];
double result_kronrod = f_center * wgk[n - 1];
double result_abs = fabs(result_kronrod);
double result_asc, mean, err;
unsigned j;
for (j = 0; j < (n - 1) / 2; ++j) {
int j2 = 2 * j + 1;
double x, f1, f2, fsum, w = width * xgk[j2];
x = center - w;
fv1[j2] = f1 = f(1, &x, fdata);
x = center + w;
fv2[j2] = f2 = f(1, &x, fdata);
fsum = f1 + f2;
result_gauss += wg[j] * fsum;
result_kronrod += wgk[j2] * fsum;
result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
}
for (j = 0; j < n / 2; ++j) {
int j2 = 2 * j;
double x, f1, f2, w = width * xgk[j2];
x = center - w;
fv1[j2] = f1 = f(1, &x, fdata);
x = center + w;
fv2[j2] = f2 = f(1, &x, fdata);
result_kronrod += wgk[j2] * (f1 + f2);
result_abs += wgk[j2] * (fabs(f1) + fabs(f2));
}
ee->val = result_kronrod * width;
/* compute error estimate: */
mean = result_kronrod * 0.5;
result_asc = wgk[n - 1] * fabs(f_center - mean);
for (j = 0; j < n - 1; ++j)
result_asc += wgk[j] * (fabs(fv1[j] - mean) + fabs(fv2[j] - mean));
err = (result_kronrod - result_gauss) * width;
result_abs *= width;
result_asc *= width;
if (result_asc != 0 && err != 0) {
double scale = pow((200 * err / result_asc), 1.5);
if (scale < 1)
err = result_asc * scale;
else
err = result_asc;
}
if (result_abs > DBL_MIN / (50 * DBL_EPSILON)) {
double min_err = 50 * DBL_EPSILON * result_abs;
if (min_err > err)
err = min_err;
}
ee->err = err;
return 0; /* no choice but to divide 0th dimension */
}
static rule *make_rule15gauss(unsigned dim)
{
rule *r;
if (dim != 1)
return 0; /* this rule is only for 1d integrals */
r = (rule *) malloc(sizeof(rule));
r->dim = dim;
r->num_points = 15;
r->evalError = rule15gauss_evalError;
r->destroy = 0;
return r;
}
/*
*========================================================================
* binary heap implementation (ala _Introduction to Algorithms_ by Cormen,
* Leiserson, and Rivest), for use as a priority queue of regions to
* integrate.
*========================================================================
*/
static void heap_resize(heap * h, unsigned nalloc)
{
h->nalloc = nalloc;
h->items = (heap_item *) realloc(h->items, sizeof(heap_item) * nalloc);
}
static heap heap_alloc(unsigned nalloc)
{
heap h;
h.n = 0;
h.nalloc = 0;
h.items = 0;
h.ee.val = h.ee.err = 0;
heap_resize(&h, nalloc);
return h;
}
/*
* note that heap_free does not deallocate anything referenced by the items
*/
static void heap_free(heap * h)
{
h->n = 0;
heap_resize(h, 0);
}
static void heap_push(heap * h, heap_item hi)
{
int insert;
h->ee.val += hi.ee.val;
h->ee.err += hi.ee.err;
insert = h->n;
if (++(h->n) > h->nalloc)
heap_resize(h, h->n * 2);
while (insert) {
int parent = (insert - 1) / 2;
if (KEY(hi) <= KEY(h->items[parent]))
break;
h->items[insert] = h->items[parent];
insert = parent;
}
h->items[insert] = hi;
}
static heap_item heap_pop(heap * h)
{
heap_item ret;
int i, n, child;
if (!(h->n)) {
fprintf(stderr, "attempted to pop an empty heap\n" );
exit(EXIT_FAILURE);
}
ret = h->items[0];
h->items[i = 0] = h->items[n = --(h->n)];
while ((child = i * 2 + 1) < n) {
int largest;
heap_item swap;
if (KEY(h->items[child]) <= KEY(h->items[i]))
largest = i;
else
largest = child;
if (++child < n && KEY(h->items[largest]) < KEY(h->items[child]))
largest = child;
if (largest == i)
break;
swap = h->items[i];
h->items[i] = h->items[largest];
h->items[i = largest] = swap;
}
h->ee.val -= ret.ee.val;
h->ee.err -= ret.ee.err;
return ret;
}
/*
*========================================================================
* adaptive integration, analogous to adaptintegrator.cpp in HIntLib
*========================================================================
*/
static int ruleadapt_integrate(rule * r, integrand f, void *fdata,
unsigned dim, const hypercube * h, unsigned maxEval,
double reqAbsError, double reqRelError, esterr * ee)
{
unsigned initialRegions; /* number of initial regions (non-adaptive) */
unsigned minIter; /* minimum number of adaptive subdivisions */
unsigned maxIter; /* maximum number of adaptive subdivisions */
unsigned initialPoints;
heap regions;
unsigned i,j,ip;
int status = -1; /* = ERROR */
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Entry.\n" );
}
initialRegions = 1; /* or: use some percentage of maxEval/r->num_points */
initialPoints = initialRegions * r->num_points;
minIter = 0;
if (maxEval) {
if (initialPoints > maxEval) {
initialRegions = maxEval / r->num_points;
initialPoints = initialRegions * r->num_points;
}
maxEval -= initialPoints;
maxIter = maxEval / (2 * r->num_points);
} else
maxIter = UINT_MAX;
if (initialRegions == 0)
return status; /* ERROR */
/*
* Allocate a heap with room for a set of region evaluations
*/
regions = heap_alloc(initialRegions);
/*
* Evaluate a region and push it onto the heap.
*/
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Make hypercube into region and put it on heap\n" );
}
heap_push(&regions, eval_region(make_region(h), f, fdata, r));
/* or:
if (initialRegions != 1)
partition(h, initialRegions, EQUIDISTANT, &regions, f,fdata, r);
*/
/*
* Now we work to convergence or maxIter, whichever comes first
*/
for (i = 0; i < maxIter; ++i) {
region R, R2;
if (i >= minIter && (regions.ee.err <= reqAbsError
|| relError(regions.ee) <= reqRelError)) {
status = 0; /* converged! */
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Converged!\n" );
}
break;
}
/*
* Get the worst region (which might be the only region)
*/
R = heap_pop(&regions); /* get worst region */
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Worst hypercube is at h.center = (" );
for(ip = 0;ip<dim;ip++){
printf("%4.2f",R.h.center[ip]);
if(ip != dim-1) printf("," );
}
printf(" )\n" );
}
/* Split it along the working (worst) dimension */
cut_region(&R, &R2);
for(j = 0;j<R.h.dim;j++){
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Split worst dimension, R max = %f R2 max = %f\n",R.h.center[j]+R.h.width[j],
R2.h.center[j]+R2.h.width[j]);
}
}
/* printf("Split it into %0x and %0x\n",R.h.center,R2.h.center); */
/* Evaluate its left half */
if(debug(D_CUBATURE)){
printf("# # ruleadapt_integrate(): Evaluate left half onto heap.\n" );
}
heap_push(&regions, eval_region(R, f, fdata, r));
/* And its right half */
if(debug(D_CUBATURE)){
printf("# # ruleadapt_integrate(): Evaluate right half onto heap.\n" );
}
heap_push(&regions, eval_region(R2, f, fdata, r));
}
/*
* Re-sum the integral and errors for all the regions on the heap
* and then remove them (clearing the heap).
*/
ee->val = ee->err = 0;
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): About to destroy %d regions.\n",regions.n);
}
for (i = 0; i < regions.n; ++i) {
/* printf("i=%d: adding value = %f, error = %f\n",i,regions.items[i].ee.val,
regions.items[i].ee.err); */
ee->val += regions.items[i].ee.val;
ee->err += regions.items[i].ee.err;
destroy_region(&regions.items[i]);
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): Destroyed region %d.\n",i);
}
}
if(debug(D_CUBATURE)){
printf("# ruleadapt_integrate(): freeing regions and done.\n" );
}
heap_free(&regions);
return status;
}
int adapt_integrate(integrand f, void *fdata,
unsigned dim, const double *xmin, const double *xmax,
unsigned maxEval, double reqAbsError,
double reqRelError, double *val,
double *estimated_error)
{
rule *r;
hypercube h;
esterr ee;
if (dim == 0) { /* trivial integration */
ee.val = f(0, xmin, fdata);
ee.err = 0;
return 0;
}
if(debug(D_CUBATURE)){
printf("# Debugging routines in cubature.c\n" );
}
/*
* Use rule15gauss for 1d integrals, rule75genzmalik for d>1
*/
r = dim == 1 ? make_rule15gauss(dim) : make_rule75genzmalik(dim);
if(debug(D_CUBATURE)){
printf("# adapt_integrate(): Make a ranged hypercube over integration volume.\n" );
}
h = make_hypercube_range(dim, xmin, xmax);
if(debug(D_CUBATURE)){
printf("# adapt_integrate(): calling ruleadapt_integrate with this hypercube.\n" );
}
/*
* We need to pass dim for debugging only
*/
int status = ruleadapt_integrate(r, f, fdata, dim, &h,maxEval, reqAbsError,
reqRelError, &ee);
if(debug(D_CUBATURE)){
printf("# adapt_integrate(): Done. Cleaning up.\n" );
}
*val = ee.val;
*estimated_error = ee.err;
destroy_hypercube(&h);
destroy_rule(r);
return status;
}

Message édité par gilou le 04-02-2011 à 15:58:07

Marsh Posté le 04-02-2011 à 16:29:48

Bon au fait, ça plante pour quelle type de valeur.
Parce que comme dans le code on a du static heap heap_alloc(unsigned nalloc)
Peut être qu'à un moment la pile va bouffer trop de mémoire, et je vois pas trop de tests de plus au niveau de l'alloc pour verifier que tout s'est bien passé :heink:
A+,

---------------
There's more than what can be linked! -- Iyashikei Anime Forever! -- AngularJS c'est un framework d'engulé! --

Marsh Posté le 04-02-2011 à 16:52:43

En fait, comme j'ai mentionné je suis pas un connaisseur en informatique et en programmation. Déjà, les deux derniers codes (cubature.c et cubature.h) c'est pas mon oeuvre. j'ai trouvé c'est deux codes sur le net.
Commet je fais des tests au niveau de l'allocation ?

Marsh Posté le 04-02-2011 à 16:52:43

Marsh Posté le 04-02-2011 à 19:20:48

Heap_alloc fait en fait appel à heap_realloc qui fait appel à realloc:

Code :

static void heap_resize(heap * h, unsigned nalloc)
{
h->nalloc = nalloc;
h->items = (heap_item *) realloc(h->items, sizeof(heap_item) * nalloc);
}

Si realloc échoue, ça va renvoyer NULL, et ça, ça devrait être testé, ce qui n'est pas le cas ici.
A+,

---------------
There's more than what can be linked! -- Iyashikei Anime Forever! -- AngularJS c'est un framework d'engulé! --

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed

Name *

Email *

URL

Message *