如何修复我的程序中的“双重自由或腐败”错误？

Question

<>我用C++编写了一个程序，使用OpenBLAS和Booo:：Python和FFTW。当我执行程序时

*** Error in `python': double free or corruption (!prev): 0x0000561652580eb0 ***
======= Backtrace: =========
/usr/lib/libc.so.6(+0x72bdd)[0x7f4a039f0bdd]
/usr/lib/libc.so.6(+0x792ec)[0x7f4a039f72ec]
/usr/lib/libc.so.6(+0x7a6d1)[0x7f4a039f86d1]
...
/usr/lib/libpython2.7.so.1.0(PyRun_SimpleFileExFlags+0x19e)[0x7f4a0368bf8e]
/usr/lib/libpython2.7.so.1.0(Py_Main+0x5fe)[0x7f4a0368f92e]
/usr/lib/libc.so.6(__libc_start_main+0xea)[0x7f4a0399e4ca]
python(_start+0x2a)[0x56165079f7ba]

我不知道如何找出错误的来源，因为这是第一次造成错误。除此之外，我用Google搜索这个问题，但找不到任何解决方案。请教我一些注释和如何修复代码。在

我的代码是：

在快速傅里叶变换水电站在

^{pr2}$

快速傅里叶变换_图像.cpp在

#include "fft.hpp"


void static_fft_imaging(data uvdata, image imfits,
  double lambda_l1, double lambda_tv, double lambda_tsv)
{
  result imresult;

  /* input lambda_l1 and tv and tsv */
  imresult.lambda_l1 = lambda_l1;
  if (lambda_tv>0) {
    imresult.lambda_tv = lambda_tv;
  }
  if (lambda_tsv>0) {
    imresult.lambda_tsv = lambda_tsv;
  }

  std::cout<<"--- start a calculation of cost and gradcost ---"<<std::endl;
  for (int i=0; i<5; ++i) {
    std::cout<<"--- "<< i+1 <<" ---"<<std::endl;
    calc_cost(&imfits, &uvdata, &imresult);
  }
  std::cout<<"--- finished!! ---"<<std::endl;

  return;
}

void calc_cost(image *imfits, data *uvdata, result *imresult)
{
  int inc=1;
  double *in(NULL),*out2(NULL),*Vreal(NULL),*Vimag(NULL),*gradVamp(NULL),*gradVph(NULL);
  bool gradV_flag=false;
  dcomplex *out;
  fftw_plan fftwplan,ifftwrplan,ifftwiplan;

  /* check if we need to calculate gradient of amplitude and phase */
  if (uvdata->iscp || uvdata->isca) {
    gradV_flag = true;
  }

  /* initialize the cost and its gradient */
  imresult->gradcost = alloc_matrix<double>(imfits->NX, imfits->NY);

  imresult->cost = 0.;
  clear_matrix(imresult->gradcost, imfits->NX, imfits->NY);

  /* allocate arrays related with the FFT */
  Vreal = alloc_vector<double>(uvdata->Nsize);
  Vimag = alloc_vector<double>(uvdata->Nsize);

  clear_matrix(Vreal, 1, uvdata->Nsize);
  clear_matrix(Vimag, 1, uvdata->Nsize);
  if (gradV_flag) {
    gradVamp = alloc_matrix<double>(imfits->Npix, uvdata->Nuv);
    gradVph = alloc_matrix<double>(imfits->Npix, uvdata->Nuv);

    clear_matrix(gradVamp, imfits->Npix, uvdata->Nuv);
    clear_matrix(gradVph, imfits->Npix, uvdata->Nuv);
  }

  /* ---------------------------------------------- */
  /* Calculate chi-square and its gradient */
  /* ---------------------------------------------- */
  /* allocate arrays of the input and output for FFT */
  in = alloc_matrix<double>(imfits->NX, imfits->NY);
  out = alloc_matrix<dcomplex>(uvdata->NU, uvdata->NV);
  out2 = alloc_matrix<double>(uvdata->NU, uvdata->NV);

  fftwplan = fftw_plan_dft_r2c_2d(imfits->NX, imfits->NY,
                                   in, reinterpret_cast<fftw_complex*>(out), FFTW_MEASURE);
  ifftwrplan = fftw_plan_r2r_2d(uvdata->NU, uvdata->NV,
                                   out2, in, FFTW_REDFT01, FFTW_REDFT01, FFTW_MEASURE);
  ifftwiplan = fftw_plan_r2r_2d(uvdata->NU, uvdata->NV,
                                   out2, in, FFTW_RODFT01, FFTW_RODFT01, FFTW_MEASURE);

  /* FFTW */
  dcopy_(&(imfits->Npix), imfits->Iin, &inc, in, &inc);
  fftw_execute(fftwplan);

  /* copy output to some arrays */
  sort_data(out, uvdata->Nuv, Vreal, Vimag);

  /* calculate chi-square and its gradient */
  /* full-comp. visibility */
  if (uvdata->isfcv) {
    std::cout<<"--- full-comp. visibility ---"<<std::endl;
    chisq_fcv(uvdata, imresult, imfits->NX, imfits->NY, out2, in, &ifftwrplan, &ifftwiplan,
          Vreal, Vimag);
  }

  /* ------------------------------------------ */
  /* Calculate regularization functions */
  /* ------------------------------------------ */
  /* calculate cost function */
  /* l1 norm */
  if ((imresult->lambda_l1)>0) {
    imresult->cost += imresult->lambda_l1*dasum_(&(imfits->Npix), imfits->Iin, &inc);
  }
  /* TV term */
  if ((imresult->lambda_tv)>0) {
    imresult->cost += imresult->lambda_tv*tv(imfits->Iin, imfits->NX, imfits->NY);
  }
  /* TSV term */
  if ((imresult->lambda_tsv)>0) {
    imresult->cost += imresult->lambda_tsv*tsv(imfits->Iin, imfits->NX, imfits->NY);
  }

  /* calculate a gradient of cost function */
  for (int ipix=0; ipix<(imfits->Npix); ++ipix) {
    /* l1 norm */
    if ((imresult->lambda_l1)>0) {
      if (imfits->Iin[ipix]>tol) {
        imresult->gradcost[ipix] += imresult->lambda_l1;
      } else if (imfits->Iin[ipix]<-1*tol) {
        imresult->gradcost[ipix] -= imresult->lambda_l1;
      }
    }
    /* TV term */
    if ((imresult->lambda_tv)>0) {
      imresult->gradcost[ipix] += imresult->lambda_tv*gradtv(ipix, imfits->Iin, imfits->NX, imfits->NY);
    }
    /* TSV term */
    if ((imresult->lambda_tsv)>0) {
      imresult->gradcost[ipix] += imresult->lambda_tsv*gradtsv(ipix, imfits->Iin, imfits->NX, imfits->NY);
    }
  }

  //std::cout <<"cost = "<< imresult->cost << std::endl;
  //std::cout <<"gradcost = "<< (imresult->gradcost)[1500] << std::endl;

  /* free memories */  
  /* double */
  if (in!=NULL) {
    delete[] in;
    in = NULL;
  }
  if (out2!=NULL) {
    delete[] out2;
    out2 = NULL;
  }
  if (Vreal!=NULL) {
    delete[] Vreal;
    Vreal = NULL;
  }
  if (Vimag!=NULL) {
    delete[] Vimag;
    Vimag = NULL;
  }
  if (imresult->gradcost!=NULL) {
    delete[] imresult->gradcost;
    imresult->gradcost = NULL;
  }

  if (gradV_flag) {
    if (gradVamp!=NULL) {
      delete[] gradVamp;
      gradVamp = NULL;
    }
    if (gradVph!=NULL) {
      delete[] gradVph;
      gradVph = NULL;
    }
  }
  /* dcomplex */
  delete[] out;
  /* fftw_plan */
  fftw_destroy_plan(fftwplan);
  fftw_destroy_plan(ifftwrplan);
  fftw_destroy_plan(ifftwiplan);
}

快速傅里叶变换_工具.cpp在

#include "fft.hpp"


/* memmory allocation of matrix and vectors */
void clear_matrix(int* matrix, int width, int height) {
  for (int idx=0; idx<(width*height); ++idx) {
    matrix[idx] = int(0);
  }
};

void clear_matrix(double* matrix, int width, int height) {
  for (int idx=0; idx<(width*height); ++idx) {
    matrix[idx] = double(0);
  }
};

void clear_matrix(dcomplex* matrix, int width, int height) {
  dcomplex zero;

  for (int idx=0; idx<(width*height); ++idx) {
    matrix[idx] = zero;
  }
};

/* convert a position between one-dimention and two-dimention */
void ipix2ixy(int ipix, int *ix, int *iy, int NX) {
  /* calculate ix and iy from ipix */
  *ix = ipix%NX;
  *iy = ipix/NX;
}

void ixy2ipix(int ix, int iy, int *ipix, int NX) {
  /* calculate ipix from ix and iy */
  *ipix = ix+iy*NX;
}

/* sort and 0-padding (inverted-sort) data */
void sort_data(dcomplex *in, int Nuv, double *real, double *imag)
{
  /* sort data */
  for (int iuv=0; iuv<Nuv; ++iuv) {
    /* extract in[ipix] to some arrays */
    real[iuv] = std::real(in[iuv]);
    imag[iuv] = std::imag(in[iuv]);
  }
}

void isort_data(double *in, int *uidx, int *vidx, int NU, int Nuv, double *out)
{
  int ipix;

  /* 0-padding data */
  for (int iuv=0; iuv<Nuv; ++iuv) {
    ixy2ipix(uidx[iuv], vidx[iuv], &ipix, NU);

    /* substitute in to out */
    out[ipix] = in[iuv];
  }
}

/* calculate the A matrix */
void calc_A(double *x, double *y, double u, double v, int Npix, double *Ar, double *Ai)
{
  int inc=1;
  double uu,vv,*phase;

  /* initialize an array of phase */
  phase = alloc_vector<double>(Npix);
  clear_matrix(phase, 1, Npix);

  /* calculate phases */
  uu = 2*M_PI*u;
  daxpy_(&Npix, &uu, x, &inc, phase, &inc);
  vv = 2*M_PI*v;
  daxpy_(&Npix, &vv, y, &inc, phase, &inc);

  /* calculate the A matrix */
  for (int ipix=0; ipix<Npix; ++ipix) {
    Ar[ipix] = cos(phase[ipix]);
    Ai[ipix] = sin(phase[ipix]);
  }

  delete[] phase;
}

/* calculate a gradient of amplitude and phase */
void grad_ampph(image *fitsimage, data *uvdata, double *Vreal, double *Vimag,
  double *gradVamp, double *gradVph)
{
  int idx,inc=1;
  double fact,*Ar,*Ai,Vamp,Vampsq;

  /* initialize the A matrix */
  Ar = alloc_vector<double>(fitsimage->Npix);
  Ai = alloc_vector<double>(fitsimage->Npix);
  clear_matrix(Ar, 1, fitsimage->Npix);
  clear_matrix(Ai, 1, fitsimage->Npix);

  for (int iuv=0; iuv<(uvdata->Nuv); ++iuv) {
    /* calculate amplitude and its square */
    Vampsq = Vreal[iuv]*Vreal[iuv]+Vimag[iuv]*Vimag[iuv];
    Vamp = sqrt(Vampsq);

    /* calculate the A matrix */
    calc_A(fitsimage->x, fitsimage->y, (uvdata->u)[iuv], (uvdata->v)[iuv], fitsimage->Npix, Ar, Ai);

    /* calculate a gradient of amplitude and phase */
    idx = iuv*fitsimage->Npix;
    /* initialize gradVamp and gradVph */
    clear_matrix(&gradVamp[idx], 1, fitsimage->Npix);
    clear_matrix(&gradVph[idx], 1, fitsimage->Npix);
    /* visibility amplitude */
    fact = Vreal[iuv]/Vamp;
    daxpy_(&(fitsimage->Npix), &fact, Ar, &inc, &gradVamp[idx], &inc);
    fact = Vimag[iuv]/Vamp;
    daxpy_(&(fitsimage->Npix), &fact, Ar, &inc, &gradVamp[idx], &inc);
    /* phase */
    fact = Vreal[iuv]/Vampsq;
    daxpy_(&(fitsimage->Npix), &fact, Ai, &inc, &gradVph[idx], &inc);
    fact = -1*Vimag[iuv]/Vampsq;
    daxpy_(&(fitsimage->Npix), &fact, Ar, &inc, &gradVph[idx], &inc);
  }

  delete[] Ar;
  delete[] Ai;
}

/* calculate a cost and its gradient */
void chisq_fcv(data *uvdata, result *imageresult, int NX, int NY,
  double *in, double *out, fftw_plan *ifftwrplan, fftw_plan *ifftwiplan,
  double *Vreal, double *Vimag)
{
  int uvidx,inc=1,
     Npix=NX*NY;
  double fact,res1,tresr,tresi;
  double *res2r(NULL),*res2i(NULL),*res22r(NULL),*res22i(NULL);

  /* initialize matrix */
  res2r = alloc_vector<double>(uvdata->Nfcv);
  res2i = alloc_vector<double>(uvdata->Nfcv);
  res22r = alloc_matrix<double>(uvdata->NU, uvdata->NV);
  res22i = alloc_matrix<double>(uvdata->NU, uvdata->NV);
  clear_matrix(res2r, 1, uvdata->Nfcv);
  clear_matrix(res2i, 1, uvdata->Nfcv);
  clear_matrix(res22r, uvdata->NU, uvdata->NV);
  clear_matrix(res22i, uvdata->NU, uvdata->NV);

  for (int ifcv=0; ifcv<(uvdata->Nfcv); ++ifcv) {
    /* calculate residuals */
    uvidx = (uvdata->uvidxfcv)[ifcv];
    tresr = Vreal[uvidx]-(uvdata->Vfcvr)[ifcv];
    tresi = Vimag[uvidx]-(uvdata->Vfcvi)[ifcv];

    res1 = (tresr*tresr+tresi*tresi)/(uvdata->Varfcv)[ifcv];
    res2r[ifcv] = tresr/(uvdata->Varfcv)[ifcv];
    res2i[ifcv] = tresi/(uvdata->Varfcv)[ifcv];

    /* add a residual to cost function */
    imageresult->cost += res1;
  }

  /* 0-padding of res2r and res2i */
  isort_data(res2r, uvdata->uidx, uvdata->vidx, uvdata->NU, uvdata->Nfcv, res22r);
  isort_data(res2i, uvdata->uidx, uvdata->vidx, uvdata->NU, uvdata->Nfcv, res22i);

  /* cos-FFTW and sin-FFTW of real and imaginary */
  double *treal(NULL),*timag(NULL);
  /* initialize matrix */
  treal = alloc_matrix<double>(NX, NY);
  timag = alloc_matrix<double>(NX, NY);
  /* real part of residuals */
  dcopy_(&(uvdata->Nuv), res22r, &inc, in, &inc);
  fftw_execute(*ifftwrplan);
  dcopy_(&(uvdata->Nuv), out, &inc, treal, &inc);
  fftw_execute(*ifftwiplan);
  dcopy_(&(uvdata->Nuv), out, &inc, timag, &inc);
  /* imaginary part of residuals */
  dcopy_(&(uvdata->Nuv), res22i, &inc, in, &inc);
  fftw_execute(*ifftwiplan);
  fact = -1.;
  daxpy_(&(uvdata->Nuv), &fact, out, &inc, treal, &inc);
  fftw_execute(*ifftwrplan);
  fact = 1.;
  daxpy_(&(uvdata->Nuv), &fact, out, &inc, timag, &inc);

  /* copy result to gradcost */
  for (int ipix=0; ipix<Npix; ++ipix) {
    imageresult->gradcost[ipix] += -2*(treal[ipix]*treal[ipix]+timag[ipix]*timag[ipix]);
  }

  if (res2r!=NULL) {
    delete[] res2r;
    res2r = NULL;
  }
  if (res2i!=NULL) {
    delete[] res2i;
    res2i = NULL;
  }
  if (res22r!=NULL) {
    delete[] res22r;
    res22r = NULL;
  }
  if (res22i!=NULL) {
    delete[] res22i;
    res22i = NULL;
  }
  if (treal!=NULL) {
    delete[] treal;
    treal = NULL;
  }
  if (timag!=NULL) {
    delete[] timag;
    timag = NULL;
  }
}

double tv (double *Iin, int NX, int NY)
{
  int ipix2,ipix3,ix,iy,ix2,iy2,
    Npix=NX*NY;
  double dIx,dIy,tv;

  tv = 0.;
  for (int ipix=0; ipix<Npix; ++ipix) {
    /* calculate ix and iy with ipix */
    ipix2ixy(ipix, &ix, &iy, NX);
    ix2 = ix+1;
    iy2 = iy+1;
    /* calculate ipix2 (ipix3) with ix2 (ix) and iy (iy2) */
    ixy2ipix(ix2, iy, &ipix2, NX);
    ixy2ipix(ix, iy2, &ipix3, NX);

    /* calculate total variation */
    if (ix2>=NX) {
      dIx = 0.;
    } else {
      dIx = Iin[ipix2]-Iin[ipix];
    }
    if (iy2>=NY) {
      dIy = 0.;
    } else {
      dIy = Iin[ipix3]-Iin[ipix];
    }

    tv += sqrt(dIx*dIx+dIy*dIy);
  }

  return tv;
}

double tsv(double *Iin, int NX, int NY)
{
  int ipix2,ipix3,ix,iy,ix2,iy2,
    Npix=NX*NY;
  double dIx,dIy,tsv;

  tsv = 0.;
  for (int ipix=0; ipix<Npix; ++ipix) {
    /* calculate ix and iy with ipix */
    ipix2ixy(ipix, &ix, &iy, NX);
    ix2 = ix+1;
    iy2 = iy+1;
    /* calculate ipix2 (ipix3) with ix2 (ix) and iy (iy2) */
    ixy2ipix(ix2, iy, &ipix2, NX);
    ixy2ipix(ix, iy2, &ipix3, NX);

    /* calculate total variation */
    if (ix2>=NX) {
      dIx = 0.;
    } else {
      dIx = Iin[ipix2]-Iin[ipix];
    }
    if (iy2>=NY) {
      dIy = 0.;
    } else {
      dIy = Iin[ipix3]-Iin[ipix];
    }

    tsv += dIx*dIx+dIy*dIy;
  }

  return tsv;
}

double gradtv(int ipix, double *Iin, int NX, int NY)
{
  int ipix2,ipix3,ix,iy,ix2,iy2,ix3,iy3;
  double dIx,dIy,tv,gradtv;

  /* calculate ix and iy with ipix */
  ipix2ixy(ipix, &ix, &iy, NX);
  ix2 = ix+1;
  ix3 = ix-1;
  iy2 = iy+1;
  iy3 = iy-1;

  /* calculate gradtv */
  /* --------------------------- */
  /* 2*(ix,iy)-(ix2,iy)-(ix,iy2) */
  /* --------------------------- */
  gradtv = 0.;
  if (ix2>=NX) {
    dIx = 0.;
  } else {
    /* calculate ipix2 with ix2 and iy */
    ixy2ipix(ix2, iy, &ipix2, NX);

    dIx = Iin[ipix]-Iin[ipix2];
  }
  if (iy2>=NY) {
    dIy = 0.;
  } else {
    /* calculate ipix2 with ix and iy2 */
    ixy2ipix(ix, iy2, &ipix2, NX);

    dIy = Iin[ipix]-Iin[ipix2];
  }

  tv = sqrt(dIx*dIx+dIy*dIy);
  if (tv>tol) {
    gradtv += (dIx+dIy)/tv;
  }
  /* ------------------------------------ */
  /* (ix,iy)-(ix3,iy), (ix3,iy2)-(ix3,iy) */
  /* ------------------------------------ */
  if (ix3>=0) {
    /* calculate ipix2 with ix3 and iy */
    ixy2ipix(ix3, iy, &ipix2, NX);

    dIx = Iin[ipix]-Iin[ipix2];
    if (iy2>=NY) {
      dIy= 0.;
    } else {
      /* calculate ipix2 (ipix3) with ix3 and iy (iy2) */
      ixy2ipix(ix3, iy2, &ipix3, NX);

      dIy = Iin[ipix3]-Iin[ipix2];
    }

    tv = sqrt(dIx*dIx+dIy*dIy);
    if (tv>tol) {
      gradtv += dIx/tv;
    }
  }
  /* ------------------------------------ */
  /* (ix,iy)-(ix,iy3), (ix2,iy3)-(ix,iy3) */
  /* ------------------------------------ */
  if (iy3>=0) {
    /* calculate ipix2 with ix and iy3 */
    ixy2ipix(ix, iy3, &ipix2, NX);

    dIy = Iin[ipix]-Iin[ipix2];
    if (ix2>=NX) {
      dIx = 0.;
    } else {
      /* calculate ipix3 with ix2 and iy3 */
      ixy2ipix(ix2, iy3, &ipix3, NX);

      dIx = Iin[ipix3]-Iin[ipix2];
    }

    tv = sqrt(dIx*dIx+dIy*dIy);
    if (tv>tol) {
      gradtv += dIy/tv;
    }
  }

  return gradtv;
}

double gradtsv(int ipix, double *Iin, int NX, int NY)
{
  int ipix2,ipix3,ix,iy,ix2,iy2,ix3,iy3;
  double gradtsv;

  /* calculate ix and iy with ipix */
  ipix2ixy(ipix, &ix, &iy, NX);
  ix2 = ix+1;
  ix3 = ix-1;
  iy2 = iy+1;
  iy3 = iy-1;

  /* calculate gradtsv */
  gradtsv = 0.;
  if (ix2<NX && ix3>=0) {
    /* calculate ipix2 (ipix3) with ix2 (ix3) and iy */
    ixy2ipix(ix2, iy, &ipix2, NX);
    ixy2ipix(ix3, iy, &ipix3, NX);

    gradtsv += Iin[ipix2]-Iin[ipix3];
  }
  if (iy2<NY && iy3>=0) {
    /* calculate ipix2 (ipix3) with ix and iy2 (iy3) */
    ixy2ipix(ix, iy2, &ipix2, NX);
    ixy2ipix(ix, iy3, &ipix3, NX);

    gradtsv += Iin[ipix2]-Iin[ipix3];
  }

  return gradtsv;
}

在libfftw.cpp文件在

#include "fft.hpp"


BOOST_PYTHON_MODULE(libfft)
{
  Py_Initialize();
  npy::initialize();

  /* define classes */
  /* data */
  py::class_<data>("uvdata")
  .def("input_pos", &data::input_pos)
  .def("input_fcv", &data::input_fcv)
  .def("input_amp", &data::input_amp)
  .def("input_cp", &data::input_cp)
  .def("input_ca", &data::input_ca);
  /* image */
  py::class_<image>("imfits")
  .def("input_im", &image::input_im);

  py::def("static_fft_imaging", &static_fft_imaging);
}