Segmentation fault Core dump

[rdekou]

I'm writing a code using least mean square to fit a sample of data. The code is running till an iteration where it crashes and returns a segmentation fault core dumped (the memory used is always closed to 2g mem=1805804kb,vmem=1926132kb, seems like the code can't go beyond) . But If i try to fit at this last iteration the code works perfectly, so i think that I'm not trying to access a variable at a wrong index but it might be dynamic memory allocation ( a pointer to a matrix that I didn't free or something like that). The code is complex and calls too many routines, I'm combining cimg , c++ and old c. How can I track this error?

[kbw]

Difficult to guess. But if you let it crash, then take a look at the call stack, you'll be able to identify where the crash happened and work backwards from there.

[rdekou]

Thanks for your reply,

Sorry but I don't understand, how can I have a look at the stack? (what is the stack?)

Here is a piece of my code my code:

//variables

{
int   max_iter=100, count; 

double somme1, somme2;

float corell, tol=1e-3;


float** covar= new float*; //the pointer should be initialised otherwise it doesn't work


CImg<double> s_m1, s_m2; 
CImg<float> a, ecarti ,ecarts, ZZ, YY, YY2, ZZ2, matrixcor, matrixparam; // a is the parameter matrix

CImg<float> lamb, sigma, lamb_ind, lamb_azi;



size_t nba=6, nbf=2; //nba nber of parameters, nbf nbre of directions




s_m1.assign(3); //2 iterations, s_m1(1 && 2)
s_m2.assign(3);
a.assign(nba);
ecarti.assign(nba);//low bound delimiter of trust region for each parameter
ecarts.assign(nba); //high bound -----/----------------------------------
matrixcor.assign(dim_piv_t, dim_piv_y, dim_piv_z);
matrixparam.assign(dim_piv_t, dim_piv_y, dim_piv_z, nba);



cimg_forXYZ(G1,tt, yy, zz){

//int tt=175, yy=38, zz=51;

if(abs(G1(tt,yy, zz))>0.73 && radius(tt,yy,zz) >=0.5){ //>=1 so that the mask is at least 1

int im=0, icy =0, icz=0;



mask = 2*radius(tt,yy,zz);

cout<<"radius"<<radius(tt,yy,zz)<<endl;


cimg_forX(dy,i){
dy(i) = 200;
dz(i) = 200;
}



if(yy+ mask >=dim_piv_y){


dy(icy) = dim_piv_y-1-yy;
icy++;
}




if(yy-mask<0){

dy(icy) = yy;
}


if(zz+ mask>=dim_piv_z){

dz(icz)=dim_piv_z-1-zz;
icz++;
}


if(zz-mask<0){


dz(icz) = zz;
}



if(!(dy(0)==200) || !(dy(1)==200)){

d1 = dy.min();
}

else d1=200; 


if(!(dz(0)==200) || !(dz(1)==200)){

d2 = dz.min();
}

else d2 =200;

if(!(dy(0)==200) || !(dy(1)==200) || !(dz(0)==200) || !(dz(1)==200)){

cout<<"a vortex close to the boundary"<<endl;

mask =  min(d1,d2);

}

else{
cout<<"a vortex not too close  the boundary"<<endl;
}





if(yy+ mask<dim_piv_y && yy-mask>=0 && zz+mask<dim_piv_z && zz-mask>=0){


do{




if(yy+ mask<dim_piv_y && yy-mask>=0 && zz+mask<dim_piv_z && zz-mask>=0){

  if(im==0){
  cout<<"****************start*********************"<<endl;
  }



cout<<"**mask="<<mask<<"**tt="<<tt<<"**yy="<<yy<<"**zz="<<zz<<"**flag="<<flag(tt,yy,zz)<<"**domy="<<domainy(tt,yy,zz)<<"**domz="<<domainz(tt,yy,zz)<<endl;



size_t  nbz= pow(2*mask+1,2);



lamb.assign(nbz, nbf); 
sigma.assign(nbz, nbf);

YY.assign(nbz);
ZZ.assign(nbz);


  cimg_forXY(sigma,i, j){
   sigma(i, j) =  1; 
  }
 
   
  //initialisation

  if(im==0){

  a(0) = vy(tt,yy,zz);  //u_conv is the velocity at the peak 
  a(1) = vz(tt,yy,zz);  //v_conv
  a(2) = yy; 
  a(3) = zz;

  
  a(4) = 0.5*sgn(G1(tt,yy,zz))*abs(  (vz(tt,yy+1,zz) - vz(tt,yy-1,zz)) -  (vy(tt,yy,zz+1) - vz(tt,yy,zz-1)))*(4/(2*g_Pi));    //circulation c'est wx.ds= (dvz/dy -dvy/dvz), wx compute with center finite difference and ds =4 , avec ou sans sgn(G1(t,y,z) 
      
  
  a(5) = pow(radius(tt,yy,zz),2);  
 }

  
  else{ //im==1;
 
  a(0) = matrixparam(tt,yy,zz,0); 
  a(1) = matrixparam(tt,yy,zz,1);
  a(2) = matrixparam(tt,yy,zz,2);
  a(3) = matrixparam(tt,yy,zz,3);
  a(4) = matrixparam(tt,yy,zz,4);
  a(5) = matrixparam(tt,yy,zz,5);
  }  

 
  cout<<"Im="<<im<<endl; 
  printf("Inputs\n"); 
  printf("a[0]=%e\n",a[0]);
  printf("a[1]=%e\n",a[1]);
  printf("a[2]=%e\n",a[2]);
  printf("a[3]=%e\n",a[3]);
  printf("a[4]=%e\n",a[4]);
  printf("a[5]=%e\n",a[5]);
  
  
  //ecarts
  
  ecarti(0) = -1; //to include 0 as a convection velocity
  ecarts(0) = 1;
  ecarti(1) = -1;
  ecarts(1) = 1; //in practice we will take the max velocity of the whole field
  ecarti(2) = a(2)-5;
  ecarts(2) = a(2)+5;
  ecarti(3) = a(3)-5;
  ecarts(3) = a(3)+5;
  ecarti(4) = -1;
  ecarts(4) = 1; //max circulation
  ecarti(5) = 0;
  ecarts(5) = pow(50,2); 

 

// cout<<"initial grid points"<<endl;
  
  count = 0;   
 
  for(int j = -mask; j <=mask; j++){
     for(int k = -mask; k <=mask; k++){

   YY(count) = a(2)+ j; //y
   ZZ(count) = a(3)+ k; //z
 
    
    count++;
    }
   }
 
  
// cout<<"Creating a Lamb oseen vortex with input parameters"<<endl;


   for(int i=0; i<=nbz-1; i++){

     lamb(i,0) =   vy(tt, YY(i),  ZZ(i));
     lamb(i,1) =   vz(tt, YY(i),  ZZ(i));
  }
  

  float* lamb1 = lamb.data(); 
  float* sigma1 = sigma.data();
  float** lamb2 = &lamb1; 
  float** sigma2 = &sigma1; 

    
  corell=     levenberg_marquardt_nf_nv( YY, ZZ, lamb2, sigma2, nbz,  nbf, a, ecarti, ecarts,  nba,  tol, covar,  max_iter, somme1, somme2); 

  
  s_m1(im+1) = somme1; 
  s_m2(im+1) = somme2;

  cout<<s_m1(im+1)<<"***"<<s_m2(im+1)<<endl;

  cout<<"Correlation="<<pow(corell,0.25)<<endl;

  
  
  matrixcor(tt,yy,zz) = corell;

  for(int k=0; k<nba; k++){
  matrixparam(tt,yy,zz,k) = a(k); 
  }
  
  if(im==1 ||  corell==-1|| pow(corell,0.25)<=0.75 || matrixparam(tt,yy,zz,2)+ mask>=dim_piv_y || matrixparam(tt,yy,zz,2)-mask<0 || matrixparam(tt,yy,zz,3)+mask>=dim_piv_z  || matrixparam(tt,yy,zz,3)-mask<0 ){
  cout<<"******************end*****************"<<endl; 
  }
 

  //destructor
 
lamb.assign(); //not needed 
sigma.assign();
ZZ.assign();
YY.assign();
  
   
  im++;
  

  }// end do
  
   
  while( im<=1 && pow(corell,0.25)>0.75 && matrixparam(tt,yy,zz,2)+ mask<dim_piv_y && matrixparam(tt,yy,zz,2)-mask>=0 && matrixparam(tt,yy,zz,3)+mask<dim_piv_z && matrixparam(tt,yy,zz,3)-mask>=0 ); //(s_m2(im)!=s_m2(im)) continue if nan

  }//end if in domain

  }//end if G1>0.73
 }//cimg_forXYZ  

In the function levenberg_marquardt_nf_nv this is how my variables are allocated and freed:

float levenberg_marquardt_nf_nv(float *x, float *y, float **z, float **sigma, size_t nbz, size_t nbf, float *a, float *ecarti, float *ecarts, size_t nba, float tol, float **covar, int max_iter, double &somme1, double &somme2)
{

float **alpha, **alpha_tmp, **alpha_shift;


//for matrix 
    alpha = new float*[nba*sizeof(float)];
  
           
              alpha_tmp = new float*[nba*sizeof(float)];

              
                for(i=0;i<nba;i++)
          {
            alpha[i] = new float[nba*sizeof(float)]; //so we can access the pointer alpha[i][j] (2d allocation)
           alpha_tmp[i] = new float[nba*sizeof(float)];

          }


 
        for (int i = 0; i < nba; i++) {
         delete [] alpha[i];
         delete [] alpha_tmp[i];
                  }

        delete[] alpha;   
        delete[] alpha_tmp;

 // for vectors

	
        float *beta=  new float[nba*sizeof(float)];

         delete[] beta;

}

The code is very long, I hope theses lines will make my question understandable.

Last: here is the output error file:

  /var/lib/torque/mom_priv/jobs/163208.server.SC: line 47: 26845 Segmentation fault      (core dumped) ./finalfit -fiv $fivort -fil $filab -fo2 $fo -blk $blk > outf_FP05.dat 2> outf_FP05.dat

These are the last lines of the ouptut file:

a vortex not too close  the boundary
****************start*********************
**mask=22**tt=786**yy=29**zz=51**flag=1**domy=0**domz=0
Im=0
Inputs
a[0]=-1.315598e-02
a[1]=1.944575e-02
a[2]=2.900000e+01
a[3]=5.100000e+01
a[4]=-6.695143e-03
a[5]=1.254704e+02
5.11594***-nan
Correlation=-nan
******************end*****************
radius21.9397
a vortex close to the boundary
****************start*********************
**mask=34**tt=423**yy=34**zz=51**flag=1**domy=0**domz=0
Im=0
Inputs
a[0]=-8.660687e-03
a[1]=5.812600e-03
a[2]=3.400000e+01
a[3]=5.100000e+01
a[4]=-4.657662e-03
a[5]=4.813512e+02
success
a[0]=-8.671571e-03
a[1]=5.801764e-03
a[2]=3.469041e+01
a[3]=5.239419e+01
a[4]=-5.366362e-03
a[5]=3.075504e+02
jter=7
covar=1.666882e-06
correl=1.133464e-01
deleted, low velocity mag
soo=1.321122e-04
szz=1.682346e+01
16.804***16.7918
Correlation=0
******************end*****************
radius10.9344
a vortex not too close  the boundary
****************start*********************
**mask=21**tt=175**yy=38**zz=51**flag=1**domy=0**domz=0

If I set the mask constant the code is running till the end but if I choose mask=2*radius (much more cases are fitted) as above the code is crashing after a number of iterations.

[ne555]

you could also run through valgrind to detect invalid read/write
and some logic errors as using uninitialized variables.


> the code is complex and calls too many routines
simplify it

[MikeyBoy]

1) Please use code tags when posting code to make it readable. It surely can't have escaped your notice that this is what other people do here?

http://www.cplusplus.com/articles/z13hAqkS/

2) I strongly recommend you learn how to use a debugger. This will enable you to examine the call stack at the point of the crash. It will also enable you to step through your code and investigate the contents of the memory at any given point, to help you analyse why problems are occurring.

[rdekou]

Thanks MikeyBoy,

I have changed my post and I add some tags. Do you have a link where I can learn how to use it faster?

[kbw]

What platform are you on?

[rdekou]

linux x86_64

[kbw]

install gdb and cgdb. cgdb is optional, but it's marginally nicer to use.

Build your program with debug (use the -g flag and remove any -O flags).

Let's say, your program is called test1. Run:

cgdb test

We just want to run the program and watch it crash; we'll run it, let it crash and show the call stack.

r

When it crashes, you'll want to see the call stack and the line number of the crash. Type:

bt

If you want to inspect different call frames, you specify the frame number. Say, we wanted to see frame 2:

frame 2

When you're done, type:

q

[ne555]

> float *beta= new float[nba*sizeof(float)];
¿why `sizeof'?
also, consider simply using http://www.cplusplus.com/reference/vector/vector/ or http://www.cplusplus.com/reference/valarray/valarray/


PS: I have no intention in manually checking that no-indented, no-commented code.

[rdekou]

kbw I am launching the code on the server and only gdb is install,

I write this small piece of code just to see how it works:

int main(int argc, char *argv[])
{


  CImg<float> lamb0,   lamb, lamb2; 

  lamb0.assign(100); 

  for(int i=100000; i<=100020; i++){

  lamb.assign(i);
  

  cimg_forX(lamb, xx){
  lamb(xx) = 1;  
  }
 
 
  float* i1 = lamb.data();
  float** i2 = &i1;
  

  lamb2.assign(i);   

  lamb.assign(); //free, useless? no?
  }     
   float *beta=  new float[100*sizeof(float)];

    delete[] beta; 

    
  
   float *beta_tmpv =  new float[100*sizeof(float)];//not free return error
   
    cout<<"end"<<endl;



	
	return 0;
}

I compile this code with:

g++ -g  example2.cpp   -lm -lpthread -lX11 -I/usr/include/netcdf-3/       -L/usr/lib64/ -lnetcdf -lnetcdf_c++ -o example2;

As it is a short code I ran it on the localhost with:

gdb example2

GNU gdb (GDB) Fedora (7.2-16.fc14)
Copyright (C) 2010 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-redhat-linux-gnu".
For bug reporting instructions, please see:
<http://www.gnu.org/software/gdb/bugs/>...
Reading symbols from /workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/example2...done.
(gdb) r
Starting program: /workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/example2 
[Thread debugging using libthread_db enabled]
end

Program exited normally.
Missing separate debuginfos, use: debuginfo-install cyrus-sasl-lib-2.1.23-12.fc14.x86_64 glibc-2.13-2.x86_64 hdf5-1.8.5.patch1-4.fc14.x86_64 keyutils-libs-1.2-6.fc12.x86_64 krb5-libs-1.8.2-6.fc14.x86_64 libX11-1.3.4-4.fc14.x86_64 libXau-1.0.6-1.fc14.x86_64 libcom_err-1.41.12-5.fc14.x86_64 libcurl-7.21.0-10.fc14.x86_64 libgcc-4.5.1-4.fc14.x86_64 libidn-1.18-1.fc14.x86_64 libselinux-2.0.96-5.fc14.x86_64 libssh2-1.2.4-1.fc14.x86_64 libstdc++-4.5.1-4.fc14.x86_64 libxcb-1.7-1.fc14.x86_64 netcdf-4.1.1-5.fc14.x86_64 nspr-4.8.6-1.fc14.x86_64 nss-3.12.7-6.fc14.x86_64 nss-softokn-freebl-3.12.10-1.fc14.x86_64 nss-util-3.12.7-2.fc14.x86_64 openldap-2.4.22-7.fc14.x86_64 openssl-1.0.0a-2.fc14.x86_64 zlib-1.2.5-2.fc14.x86_64
(gdb) bt
No stack.
(gdb) q

So is there a way to call gdb on a server? I launch my codes on the server (qsub ./submit) via a sh file where I set the required memory.

[rdekou]

ne555 valgrind seems to work for the test case

With the code:

int main(int argc, char *argv[])
{


  cout<<"program to test valgrind"<<endl;//http://www.cprogramming.com/debugging/valgrind.html

  CImg<float> lamb0,   lamb, lamb2; 

  lamb0.assign(100); 

  for(int i=100000; i<=100020; i++){

  lamb.assign(i);
  

  cimg_forX(lamb, xx){
  lamb(xx) = 1;  
  }
 
 
  float* i1 = lamb.data();
  float** i2 = &i1;
  

  lamb2.assign(i);   

  lamb.assign(); //free, useless? no?
  }     
   float *beta=  new float[100*sizeof(float)];

    delete[] beta; 

    
  
   float *beta_tmpv =  new float[100*sizeof(float)];//not free return error
   
    cout<<"end"<<endl;



	
	return 0;
}

launch with:

valgrind --tool=memcheck --leak-check=yes ./example1 > outex_FP05.dat 2>

a cat outex_FP05.dat returns:

program to test valgrind
end
rror detector
==9206== Copyright (C) 2002-2009, and GNU GPL'd, by Julian Seward et al.
==9206== Using Valgrind-3.5.0 and LibVEX; rerun with -h for copyright info
==9206== Command: ./example1
==9206== 
==9206== 
==9206== HEAP SUMMARY:
==9206==     in use at exit: 1,600 bytes in 1 blocks
==9206==   total heap usage: 47 allocs, 46 frees, 16,805,968 bytes allocated
==9206== 
==9206== 1,600 bytes in 1 blocks are definitely lost in loss record 1 of 1
==9206==    at 0x4A054F2: operator new[](unsigned long) (vg_replace_malloc.c:264)
==9206==    by 0x40C606: main (example1.cpp:63)
==9206== 
==9206== LEAK SUMMARY:
==9206==    definitely lost: 1,600 bytes in 1 blocks
==9206==    indirectly lost: 0 bytes in 0 blocks
==9206==      possibly lost: 0 bytes in 0 blocks
==9206==    still reachable: 0 bytes in 0 blocks
==9206==         suppressed: 0 bytes in 0 blocks
==9206== 
==9206== For counts of detected and suppressed errors, rerun with: -v
==9206== ERROR SUMMARY: 1 errors from 1 contexts (suppressed: 6 from 6)

But with the real code too much memory is asked and valgrind can't make check it I have:

9== Using Valgrind-3.5.0 and LibVEX; rerun with -h for copyright info
==20099== Command: ./finalfit -fiv /workdir/er2/dekou/dekou.copie/roux/data/FP05TST10/LSE/FP05TST10YHBLK0000-0533/multitime/raoul/rnine_lines/FP05TST10YHBLK0000-0533_rec_0.2khz_BLK0522_vortex.nc -fil /workdir/er2/dekou/dekou.copie/roux/data/FP05TST10/LSE/FP05TST10YHBLK0000-0533/multitime/raoul/rnine_lines/FP05TST10YHBLK0000-0533_rec_0.2khz_BLK0522_vortex_label.nc -fo2 /workdir/er2/foucaut/pourRaoul/FP05TST10/LSE/FP05TST10YHBLK0000-0533/multitime/raoul/rnine_lines/o2khz/ -blk 522
==20099== 
--20099:0:aspacem  Valgrind: FATAL: VG_N_SEGMENTS is too low.
--20099:0:aspacem    Increase it and rebuild.  Exiting now.

I have checked around and they recommand to rise VG_N_SEGMENTS in the source and rebuild valgrind (http://comments.gmane.org/gmane.comp.debugging.valgrind/11337)

I am waiting for my admin to do it but he said that it is tricky because the installation has too many erros.

[kbw]

So is there a way to call gdb on a server?

Build the server code for debug and run it in the debugger.

Your test didn't crash, so there's nothing to see when the process is not running. You can run it, then stop it with ctrl-c. You can you bt to see where it is, then continue execution with c.

[rdekou]

I just test my real code with gdb I have this in return:

radius19.9302
a vortex close to the boundary
****************start*********************
**mask=11**tt=128**yy=38**zz=11**flag=1**domy=0**domz=0
Im=0
Inputs
a[0]=3.954448e-03
a[1]=-1.161371e-02
a[2]=3.800000e+01
a[3]=1.100000e+01
a[4]=6.529239e-03
a[5]=3.972141e+02
success
a[0]=3.957173e-03
a[1]=-1.165525e-02
a[2]=4.021535e+01
a[3]=1.146497e+01
a[4]=1.238999e-02
a[5]=2.504108e+02
jter=9
covar=1.162150e-05
correl=6.809281e-01
deleted, low velocity mag
soo=7.406194e-05
szz=2.997775e+00
2.97564***2.97325
Correlation=0
******************end*****************
radius1.9794
a vortex not too close  the boundary
****************start*********************
**mask=3**tt=446**yy=40**zz=11**flag=1**domy=0**domz=0
Im=0
Inputs
a[0]=3.527341e-03
a[1]=-7.244069e-03
a[2]=4.000000e+01
a[3]=1.100000e+01
a[4]=7.067923e-03
a[5]=3.918014e+00

Program received signal SIGSEGV, Segmentation fault.
0x0000003094a776aa in _int_malloc () from /lib64/libc.so.6
Missing separate debuginfos, use: debuginfo-install cyrus-sasl-lib-2.1.23-12.fc14.x86_64 glibc-2.13-2.x86_64 hdf5-1.8.5.patch1-4.fc14.x86_64 keyutils-libs-1.2-6.fc12.x86_64 krb5-libs-1.8.2-6.fc14.x86_64 libX11-1.3.4-4.fc14.x86_64 libXau-1.0.6-1.fc14.x86_64 libcom_err-1.41.12-5.fc14.x86_64 libcurl-7.21.0-10.fc14.x86_64 libgcc-4.5.1-4.fc14.x86_64 libidn-1.18-1.fc14.x86_64 libselinux-2.0.96-5.fc14.x86_64 libssh2-1.2.4-1.fc14.x86_64 libstdc++-4.5.1-4.fc14.x86_64 libxcb-1.7-1.fc14.x86_64 netcdf-4.1.1-5.fc14.x86_64 nspr-4.8.6-1.fc14.x86_64 nss-3.12.7-6.fc14.x86_64 nss-softokn-freebl-3.12.10-1.fc14.x86_64 nss-util-3.12.7-2.fc14.x86_64 openldap-2.4.22-7.fc14.x86_64 openssl-1.0.0a-2.fc14.x86_64 zlib-1.2.5-2.fc14.x86_64
(gdb) 
(gdb) bt
#0  0x0000003094a776aa in _int_malloc () from /lib64/libc.so.6
#1  0x0000003094a79ffd in malloc () from /lib64/libc.so.6
#2  0x000000309eeba78d in operator new(unsigned long) ()
   from /usr/lib64/libstdc++.so.6
#3  0x000000309eeba8a9 in operator new[](unsigned long) ()
   from /usr/lib64/libstdc++.so.6
#4  0x000000000040682c in submatrix (a=0x685490, oldrl=0, oldrh=5, oldcl=0, oldch=5, 
    newrl=1, newcl=1) at ./utils.h:988
#5  0x000000000040716c in levenberg_marquardt_nf_nv (x=0x6852f0, y=0x6853c0, 
    z=0x7fffffffdb28, sigma=0x7fffffffdb20, nbz=49, nbf=2, a=0x681280, 
    ecarti=0x6812a0, ecarts=0x6813a0, nba=6, tol=0.00100000005, covar=0x681220, 
    max_iter=100, somme1=@0x7fffffffdb48, somme2=@0x7fffffffdb40) at ./utils.h:1163
#6  0x000000000040b2ae in main (argc=9, argv=0x7fffffffe028) at finalfit.cpp:584

I went in my fonction levenberg_marquardt_nf_nv

There is a line where I call

float levenberg_marquardt_nf_nv(float *x, float *y, float **z, float **s     igma, size_t nbz, size_t nbf, float *a, float *ecarti, float *ecarts, si     ze_t nba, float tol, float **covar, int max_iter, double &somme1, double      &somme2){

float **alpha, **alpha_tmp, **alpha_shift;

   alpha = new float*[nba*sizeof(float)];


              alpha_tmp = new float*[nba*sizeof(float)];


                for(i=0;i<nba;i++)
          {
            alpha[i] = new float[nba*sizeof(float)]; //so we can access the pointer alpha[i][j] (2d allocation)
           alpha_tmp[i] = new float[nba*sizeof(float)];

          }




alpha_shift=submatrix(alpha,0,nba-1,0, nba-1,1,1);

//and I delete it as: 

        for (int i = 0; i < nba; i++) {
         delete [] alpha[i];
         delete [] alpha_tmp[i];
                  }

        delete[] alpha;
        delete[] alpha_tmp;



 delete[] alpha_shift;

here is the routine submatrix:

float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
       long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */

{
        long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;


        /* allocate array of pointers to rows */

        float **m= new float*[(size_t) ((nrow+NR_END)*sizeof(float*))];


        if (!m) nrerror( "allocation failure in submatrix()");
        m += NR_END;
        m -= newrl;

        /* set pointers to rows */
        for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;

        /* return pointer to array of pointers to rows */
        return m;
}

[ne555]

I'll ask again, ¿why are you using `sizeof' there?
it seems that you are confusing `new' with `malloc()'

//#4  0x000000000040682c in submatrix (a=0x685490, oldrl=0, oldrh=5, oldcl=0, oldch=5, 
//newrl=1, newcl=1) at ./utils.h:988
        long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
//nrow = 6
//ncol = -1

        m += NR_END;
        m -= newrl;
//¿? ¿do you reset it after, so it can be properly deleted?


for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;
//In the first iteration
//m[1] = a[0] - 1
//which is probably non-dereferenceable 

you are invoking quite a lot undefined behaviour.


By the way, CImg probably has some way to operate with region of interest.
and again, consider using std::vector

[rdekou]

I took the routines from someone else and made some changes.

I think the lines

m += NR_END;
m -= newrl;

are useless as I have a #define NR_END 1 and newrl is always one:

alpha_shift=submatrix(alpha,0,nba-1,0, nba-1,1,1);

I also think that +ncol is useless as a[0]-1 can't be dereferenced, further in the code I have:

 ludcmp(alpha_shift, nba, indx-1, &d);
 lubksb(alpha_shift, nba, indx-1, beta-1);

with:

int ludcmp(float **a, int n, int *indx, float *d)
{
        int i,imax,j,k;
        float big,dum,sum,temp;
        float *vv;

        vv=vector2(1,n);
        *d=1.0;
        for (i=1;i<=n;i++) {
                big=0.0;
                for (j=1;j<=n;j++)
                        if ((temp=fabs(a[i][j])) > big) big=temp;
                if (big == 0.0) return 1; //nrerror("Singular matrix in routine ludcmp");
                vv[i]=1.0/big;
        }
        for (j=1;j<=n;j++) {
                for (i=1;i<j;i++) {
                        sum=a[i][j];
                        for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
                        a[i][j]=sum;
                }
                big=0.0;
                for (i=j;i<=n;i++) {
                        sum=a[i][j];
                        for (k=1;k<j;k++)
                                sum -= a[i][k]*a[k][j];
                        a[i][j]=sum;
                        if ( (dum=vv[i]*fabs(sum)) >= big) {
                                big=dum;
                                imax=i;
                        }
                }
                if (j != imax) {
                        for (k=1;k<=n;k++) {
                                dum=a[imax][k];
                                a[imax][k]=a[j][k];
                                a[j][k]=dum;
                        }
                        *d = -(*d);
                        vv[imax]=vv[j];
                }
                indx[j]=imax;
                if (a[j][j] == 0.0) a[j][j]=TINY;
                if (j != n) {
                        dum=1.0/(a[j][j]);
                        for (i=j+1;i<=n;i++) a[i][j] *= dum;
                }
        }
        free_vector(vv,1,n);
        return 0;
}

void lubksb(float **a, int n, int *indx, float b[])
{
        int i,ii=0,ip,j;
        float sum;

        for (i=1;i<=n;i++) {
                ip=indx[i];
                sum=b[ip];
                b[ip]=b[i];
                if (ii)
                        for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
                else if (sum) ii=i;
                b[i]=sum;
        }
        for (i=n;i>=1;i--) {
                sum=b[i];
                for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
                b[i]=sum/a[i][i];
        }
}

What do you mean by reseting the pointer, so it can be proprely deleted?

do you mean doing float **m= new float*[(size_t) ((nrow+NR_END))]; again? no


I removed sizeof in new and then launch the code, it goes further and then crash:

**mask=11**tt=304**yy=11**zz=57**flag=1**domy=0**domz=0
Im=1
Inputs
a[0]=2.364184e-02
a[1]=-2.409885e-02
a[2]=1.208194e+01
a[3]=5.863020e+01
a[4]=-8.315042e-01
a[5]=4.338467e+01
success
a[0]=2.522490e-02
a[1]=-2.913467e-02
a[2]=1.221958e+01
a[3]=5.945376e+01
a[4]=-9.363349e-01
a[5]=5.342914e+01
jter=3
covar=2.700942e-03
correl=5.554838e-01
kept
soo=2.855380e+00
szz=5.148330e+00
2.29049***2.28852
Correlation=0.863312
******************end*****************
radius1.50393
a vortex not too close  the boundary
****************start*********************
**mask=3**tt=117**yy=14**zz=57**flag=1**domy=0**domz=0

Program received signal SIGSEGV, Segmentation fault.
0x0000003094a776aa in _int_malloc () from /lib64/libc.so.6
Missing separate debuginfos, use: debuginfo-install cyrus-sasl-lib-2.1.23-12.fc14.x86_64 glibc-2.13-2.x86_64 hdf5-1.8.5.patch1-4.fc14.x86_64 keyutils-libs-1.2-6.fc12.x86_64 krb5-libs-1.8.2-6.fc14.x86_64 libX11-1.3.4-4.fc14.x86_64 libXau-1.0.6-1.fc14.x86_64 libcom_err-1.41.12-5.fc14.x86_64 libcurl-7.21.0-10.fc14.x86_64 libgcc-4.5.1-4.fc14.x86_64 libidn-1.18-1.fc14.x86_64 libselinux-2.0.96-5.fc14.x86_64 libssh2-1.2.4-1.fc14.x86_64 libstdc++-4.5.1-4.fc14.x86_64 libxcb-1.7-1.fc14.x86_64 netcdf-4.1.1-5.fc14.x86_64 nspr-4.8.6-1.fc14.x86_64 nss-3.12.7-6.fc14.x86_64 nss-softokn-freebl-3.12.10-1.fc14.x86_64 nss-util-3.12.7-2.fc14.x86_64 openldap-2.4.22-7.fc14.x86_64 openssl-1.0.0a-2.fc14.x86_64 zlib-1.2.5-2.fc14.x86_64
(gdb) bt
#0  0x0000003094a776aa in _int_malloc () from /lib64/libc.so.6
#1  0x0000003094a79ffd in malloc () from /lib64/libc.so.6
#2  0x000000309eeba78d in operator new(unsigned long) ()
   from /usr/lib64/libstdc++.so.6
#3  0x000000309eeba8a9 in operator new[](unsigned long) ()
   from /usr/lib64/libstdc++.so.6
#4  0x0000000000417e7d in cimg_library::CImg<float>::assign (this=0x7fffffffd920, 
    size_x=49, size_y=1, size_z=1, size_c=1)
    at /workdir/er2/dekou/dekou.copie/roux/code.LSE/CImg-1.5.7/CImg.h:10027
#5  0x000000000040a775 in main (argc=9, argv=0x7fffffffe028) at finalfit.cpp:487
(gdb) 

which correspond to ZZ.assign(nbz);

If I remove the +ncol the code crash at the beginning:

**mask=4**tt=589**yy=12**zz=4**flag=1**domy=0**domz=0
Im=0
Inputs
a[0]=8.656941e-03
a[1]=4.227287e-04
a[2]=1.200000e+01
a[3]=4.000000e+00
a[4]=-1.177092e-02
a[5]=4.383344e+01
*** glibc detected *** /workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit: free(): invalid next size (fast): 0x0000000000681c80 ***
======= Backtrace: =========
/lib64/libc.so.6[0x3094a76d63]
/workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit[0x407f63]
/workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit[0x40b277]
/lib64/libc.so.6(__libc_start_main+0xfd)[0x3094a1ee5d]
/workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit[0x403b19]
======= Memory map: ========
00400000-0047e000 r-xp 00000000 08:41 539390113                          /workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit
0067d000-0067e000 rw-p 0007d000 08:41 539390113                          /workdir/er2/dekou/dekou.copie/roux/code.LSE/fit-levenberg2/finalfit
0067e000-00710000 rw-p 00000000 00:00 0                                  [heap]
3094600000-309461f000 r-xp 00000000 08:03 294921                         /lib64/ld-2.13.so
309481e000-309481f000 r--p 0001e000 08:03 294921                         /lib64/ld-2.13.so
309481f000-3094820000 rw-p 0001f000 08:03 294921                         /lib64/ld-2.13.so
3094820000-3094821000 rw-p 00000000 00:00 0 
3094a00000-3094b91000 r-xp 00000000 08:03 294924                         /lib64/libc-2.13.so
3094b91000-3094d91000 ---p 00191000 08:03 294924                         /lib64/libc-2.13.so
3094d91000-3094d95000 r--p 00191000 08:03 294924                         /lib64/libc-2.13.so
3094d95000-3094d96000 rw-p 00195000 08:03 294924                         /lib64/libc-2.13.so
3094d96000-3094d9c000 rw-p 00000000 00:00 0 
3094e00000-3094e84000 r-xp 00000000 08:03 294948                         /lib64/libm-2.13.so
3094e84000-3095083000 ---p 00084000 08:03 294948                         /lib64/libm-2.13.so
3095083000-3095084000 r--p 00083000 08:03 294948                         /lib64/libm-2.13.so

Thanks.

[rdekou]

This is where the routines come from https://www.cfa.harvard.edu/~sasselov/rec/code/nrutil.c:

The original version of submatrix with malloc is:

float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch,
	long newrl, long newcl)
/* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */
{
	long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl;
	float **m;

	/* allocate array of pointers to rows */
	m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*)));
	if (!m) nrerror("allocation failure in submatrix()");
	m += NR_END;
	m -= newrl;

	/* set pointers to rows */
	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol;

	/* return pointer to array of pointers to rows */
	return m;
}

[kbw]

ne555 has been looking at your code.

You still haven't answered why you're using constructs like:

alpha_tmp = new float*[nba*sizeof(float)];

How many elements are in that array? Don't you really mean:

alpha_tmp = new float*[nba];

Any why are you messing around with allocated arrays when you have a rich set of standard collections and supporting infrastructure?

If you're using notions like matrix and submatrix, why aren't they encalsulated in objects?

The larger C programs get, the more careful and regimented in your methods you have be. C++ lifts that burden from you but you have to use the relevant language constructs to achieve those benefits.

[rdekou]

Check my previous post, I removed sizeof everywhere the code goes further and then crash.



nba is the number of parameters i'm fitting I set it to 6.

I'm extending someone else code, writing it from the beginning or changing it completely will need more time.

why aren't they encalsulated in objects? sorry I didn't get it.

[kbw]

why aren't they encalsulated in objects? sorry I didn't get it.

You're processing matrix and submatrix things. If they were wrapped in classes, you could use them with confidence; knowing that they manage they're own memory correctly amonst other things.

But as all that code is effictively inlined, there's a co-mingling of matix memory management, submatrix memory management, their indexing methods, copy semantics. It's difficult to unpick what is what, and where the problem is.