My program works fine in debug mode but when I try it in release, it always crash when arriving when I use a calloc. The weird thing is that I use exactly the same way to write my calloc several times in a row and the crash only occurs the third time...
This is a bit of my code for you to get what I mean :
#include <stdio.h>
#include <math.h>
#include <dir.h>
#include <stdlib.h>
#include "Normes.h"
#define EPS 1.0e-14
void atimes(int n, float x[], float r[], int itrnsp, int *ija, float *sa);
void dsprsax(float *sa, int *ija, float *x, float *b, int n);
void dsprstx(float* sa, int *ija, float *x, float *b, int n);
float produit_scalaire(float *x, float *y, int n);
void asolve(int n, float b[], float x[], int itrnsp, int *ija, float *sa);
int linbcg(int n, float b[], float x[], int itol, float tol, int itmax, int *iter, float *err, int *ija, float *sa)
{
*iter = 0;
int flag;
float bnrm2;
float omega;
float rho;
float beta;
float rho_1;
float alpha;
float resid;
int i;
float *r=NULL;
float *vecteur_temporaire=NULL;
float *r_tilde=NULL;
float *p=NULL;
float *p_hat=NULL;
float *v=NULL;
float *s=NULL;
float *s_hat=NULL;
float *t=NULL;
vecteur_temporaire=calloc((n+1), sizeof(float));
if (vecteur_temporaire==NULL)
{
exit(0);
}
printf("\nok1\n");
r=calloc((n+1), sizeof(float));
if (r==NULL)
{
exit(0);
}
printf("\nok2\n");
r_tilde=calloc((n+1), sizeof(float));
if (r_tilde==NULL)
{
exit(0);
}
printf("\nok3\n");
p=calloc((n+1), sizeof(float));
if (p==NULL)
{
exit(0);
}
p_hat=calloc((n+1), sizeof(float));
if (p_hat==NULL)
{
exit(0);
}
v=calloc((n+1), sizeof(float));
if (v==NULL)
{
exit(0);
}
s=calloc((n+1), sizeof(float));
if (s==NULL)
{
exit(0);
}
s_hat=calloc((n+1), sizeof(float));
if (s_hat==NULL)
{
exit(0);
}
t=calloc((n+1), sizeof(float));
if (t==NULL)
{
exit(0);
}
bnrm2 = norme2(b, n);
if (bnrm2==0)
{
bnrm2=1;
}
atimes(n, x, vecteur_temporaire, 0, ija, sa);
for (i=1; i<=n; i++)
{
r[i] = b[i] - vecteur_temporaire[i];
}
*err = norme2(r, n)/bnrm2;
if (*err < tol) return 0;
omega = 1.0;
for (i=1; i<=n; i++)
{
r_tilde[i] = r[i];
}
while (*iter<=itmax)
{
(*iter)++;
rho = produit_scalaire(r, r_tilde, n);
if (rho==0)
{
return 0;
}
if (*iter>1)
{
beta = (rho/rho_1)*(alpha/omega);
for (i=1; i<=n; i++)
{
p[i] = r[i]+beta*(p[i] - omega*v[i]);
}
}
else
{
for (i=1; i<=n; i++)
{
p[i] = r[i];
}
}
//printf("\nrho=%.2f iter=%d p[1]=%.2f", rho, *iter, p[1]);
asolve(n, p, p_hat, 0, ija, sa);
atimes(n, p_hat, v, 0, ija, sa);
alpha = rho/(produit_scalaire(r_tilde, v, n));
for (i=1; i<=n; i++)
{
s[i] = r[i] - alpha*v[i];
}
if (norme2(s, n)<tol)
{
for (i=1; i<=n; i++)
{
x[i] = x[i] + alpha*p_hat[i];
}
resid = norme2(s, n)/bnrm2;
}
asolve(n, s, s_hat, 0, ija, sa);
atimes(n, s_hat, t, 0, ija, sa);
//printf("\nnorme(s) = %.2f , norme(s_hat) = %.2f\n\n", norme2(s, n), norme2(s_hat, n));
//printf("\nt = (%.2f,%.2f,%.2f,%.2f) omega = %.2f", t[1], t[2], t[3], t[4], omega);
omega = produit_scalaire(t, s, n)/produit_scalaire(t, t, n);
//printf("\t... %.2f", omega);
for (i=1; i<=n; i++)
{
x[i] = x[i] + alpha*p_hat[i]+omega*s_hat[i];
r[i] = s[i] - omega*t[i];
}
*err = norme2(r, n)/bnrm2;
//printf("\niter = %d , err = %.2f , numerr = %.2f , omega = %.2f, t = (%.2f, %.2f, %.2f, %.2f)", *iter, *err, norme2(r, n), omega, t[1], t[2], t[3], t[4]);
if (*err<=tol)
{
return 0;
}
if (omega == 0)
{
return 0;
}
rho_1=rho;
}
if (*err<=tol) //manque une condition mais on essaie comme ça
{
flag =0;
}
else if (omega == 0.0)
{
flag = -2;
}
else if (rho == 0.0)
{
flag = -1;
}
else
{
flag =1;
}
free(r);
free(r_tilde);
free(p);
free(p_hat);
free(v);
free(s);
free(s_hat);
free(t);
free(vecteur_temporaire);
return flag;
}
/*flag indique la réussite de l'algorithme :
flag = 0 : solution trouvee au taux d'errer impose
flag = 1 : pas de convergence apres le nombre maximum de convergences autorisees
flag = -1 : breakdown : rho = 0
flag = -2 : breakdown : omega = 0
*/
void atimes(int n, float x[], float r[], int itrnsp, int *ija, float *sa)
{
void dsprsax(float sa[], int ija[], float x[], float b[], int n);
void dsprstx(float sa[], int ija[], float x[], float b[], int n);
//These are float versions of sprsax and sprstx. if (itrnsp) dsprstx(sa,ija,x,r,n);
if (itrnsp) dsprstx(sa, ija, x, r, n);
else dsprsax(sa, ija, x, r, n);
}
void dsprsax(float *sa, int *ija, float *x, float *b, int n)
{
void nrerror();
int i,k;
if (ija[1] != n+2) printf("dsprsax: mismatched vector and matrix");
for (i=1;i<=n;i++) {
b[i]=sa[i]*x[i];
for (k=ija[i];k<=ija[i+1]-1;k++) b[i] += sa[k]*x[ija[k]];
}
}
void dsprstx(float* sa, int *ija, float *x, float *b, int n)
{
void nrerror();
int i,j,k;
if (ija[1] != n+2) printf("mismatched vector and matrix in dsprstx");
for (i=1;i<=n;i++) b[i]=sa[i]*x[i];
for (i=1;i<=n;i++) {
for (k=ija[i];k<=ija[i+1]-1;k++) {
j=ija[k];
b[j] += sa[k]*x[i];
}
}
}
float produit_scalaire(float *x, float *y, int n)
{
float ans = 0;
int i;
for (i=1; i<=n; i++)
{
ans = ans + x[i]*y[i];
}
return ans;
}
void asolve(int n, float b[], float x[], int itrnsp, int *ija, float *sa)
{
int i;
for(i=1;i<=n;i++) x[i]=b[i]; /*(sa[i]!=0 ? b[i]/sa[i] : b[i]);*/
//The matrix A is the diagonal part of A, stored in the first n elements of sa. Since the transpose matrix has the same diagonal, the flag itrnsp is not used.
return;
}
So here I initialize several pointers then try to allocate them, it works fine for the first two (vecteur_temporaire and r) but fail when trying to allocate the third, r_tilde (ie. on my console, ok1 and ok2 are printed but the program crash before printing ok3).
I don't think I lack memory (here n = 21 so that's not so bad, and I tested printing something if r==NULL but nothing is printed).
Do you know what can be the cause of such a problem? I'm not really surprised it works in debug mode (apparently it's usual), but I can't get why it works twice but not the third time. Thanks for your help !