c - 调用库函数后无法访问/修改所需的指针值



我正在使用外部库 https://people.sc.fsu.edu/~jburkardt/c_src/jacobi_eigenvalue/jacobi_eigenvalue.c 来计算矩阵的特征值。但是,在调用库函数后,当我尝试访问所需的指针值时,我收到了分段错误。

void asphericity(particle *p, int N, float aveLTail[MAX2][MAX2], double *Nominator, double *Denominator, double *Asp ,double *cmx, double *cmy, double *cmz){
int i, j, k,m,n,M;
double x[10000];
double y[10000];
double z[10000];
double *Q;
double gyration[3][3];
double v[M*M];
double d[M];
float ut, up;
float bs;
double hdistance;
double tdistance;
double rog; 
int hcount=0;
int tcount=0;
int it_max;
int it_num;
int rot_num;
double hx=0;
double hy=0;
double hz=0;
double tx=0;
double ty=0;
double tz=0;
m=0;

it_max=100;
for(i = 1; i<=N; i++) {
   if (p[i].mol == p[i-1].mol){
       if (p[i-1].type==1 ||p[i-1].type==2 || p[i-1].type==3 || p[i-1].type==6){
            hcount++;
            hx += p[i-1].x;
            hy += p[i-1].y;
            hz += p[i-1].z;
       }
       else if (p[i-1].type==4 && p[i-2].type==9){
            tcount++;
            tx += p[i-1].x;
            ty += p[i-1].y;
            tz += p[i-1].z;
            hcount += 2;
            hx += p[i-2].x+p[i-3].x;
            hy += p[i-2].y+p[i-3].y;
            hz += p[i-2].z+p[i-3].z;
        }  
       else if (p[i-1].type==4 && p[i-2].type!=9){
            tcount++;
            tx += p[i-1].x;
            ty += p[i-1].y;
            tz += p[i-1].z;
       }
   }
   else if (p[i].mol != p[i-1].mol && p[i-1].mol!=0) {
            tcount++;
            tx += p[i-1].x;
            ty += p[i-1].y;
            tz += p[i-1].z;
           (hx) /= hcount; /*geometric center of head beads*/
           (hy) /= hcount;
           (hz) /= hcount;
           (tx) /= tcount;
           (ty) /= tcount;
           (tz) /= tcount;
           hcount=0;
           tcount=0;
           hdistance=sqrt(pow((hx-*cmx),2)+pow((hy-*cmy),2)+pow((hz-*cmz),2));
           tdistance=sqrt(pow((tx-*cmx),2)+pow((ty-*cmy),2)+pow((tz-*cmz),2));
           if(hdistance>tdistance){
              x[m]=hx;
              y[m]=hy;
              z[m]=hz;
              m++;
             }
           hx=0;
           hy=0;
           hz=0;
           tx=0;
           ty=0;
           tz=0;
         }              
   } 

for(j = 0; j<m; j++) {
    gyration[0][0]+= pow((x[j]-*cmx),2);
    gyration[0][1]+=(x[j]-*cmx)*(y[j]-*cmy);
    gyration[0][2]+=(x[j]-*cmx)*(z[j]-*cmz); 
    gyration[1][0]+=(y[j]-*cmy)*(x[j]-*cmx);
    gyration[1][1]+=pow((y[j]-*cmy),2);
    gyration[1][2]+=(y[j]-*cmy)*(z[j]-*cmz);
    gyration[2][0]+=(z[j]-*cmz)*(x[j]-*cmx);
    gyration[2][1]+=(z[j]-*cmz)*(y[j]-*cmy);
    gyration[2][2]+=pow((z[j]-*cmz),2);
 }
for(i =0; i<3; i++){
   for(j=0; j<3; j++){
       gyration[i][j]/=m;
      }
}
 #define M 3
double Gyration[M*M] = { 
          gyration[0][0],gyration[0][1],gyration[0][2], 
          gyration[1][0],gyration[1][1],gyration[1][2],
          gyration[2][0],gyration[2][1],gyration[2][2] };
printf("%.6ft%.6ft%.6fn", *Denominator, *Nominator, *Asp);
jacobi_eigenvalue(M, Gyration, it_max, v, d, &it_num, &rot_num);
printf("%.6ft%.6ft%.6fn", *Denominator, *Nominator, *Asp);
rog = sqrt(d[0]+d[1]+d[2]);  
*Denominator= 2*pow(rog,4); 
*Nominator= pow((d[2]-d[1]),2)+pow((d[2]-d[0]),2)+pow((d[1]-d[0]),2);
* Asp = (*Nominator)/(*Denominator);
# undef M
 }

分母、提名人和 Asp 在 main 函数中被定义为双精度变量,并作为指针传递。但是,当我在库函数调用后尝试访问它们时,我收到了分段错误。我能够访问输出 v[] 和 d[]。

第一个 printf 显示 *分母、*提名符和 *Asp 零,而分割错误出现在第二个 printf 处。 因此,库函数中的某些内容不正确。 但是,我不确定问题到底在哪里以及如何解决它。

提前谢谢你。

代码的摘要视图揭示了几个潜在的"除以零"错误。关注*分母,以下是潜在的陷阱:

* Asp = (*Nominator)/(*Denominator);

如果*分母为零,那么您将崩溃(除以零)。此外,您在 for 循环中也有这个:

int hcount=0;
// some code
if (p[i].mol == p[i-1].mol){ 
  if (p[i-1].type==1 ||p[i-1].type==2 || p[i-1].type==3 || p[i-1].type==6){
     hcount++; 
  }
  else {
     // hcount incremented
  }
  // some code
}
else if (p[i].mol != p[i-1].mol && p[i-1].mol!=0) {
        tcount++;
        tx += p[i-1].x;
        ty += p[i-1].y;
        tz += p[i-1].z;
       (hx) /= hcount; /*geometric center of head beads*/
       (hy) /= hcount;
       (hz) /= hcount;
       (tx) /= tcount;
       (ty) /= tcount;
       (tz) /= tcount;
  }
请注意,在 else 中,你除以 hcount;

hcount 在声明时初始化为零,并且仅在 for 循环中的第一个 if 语句中递增;如果第一次迭代不满足第一个 if 条件并在 else 中找到它的方式,那么你将除以零。

另外,不需要

p[i].mol != p[i-1].mol

else if (p[i].mol != p[i-1].mol && p[i-1].mol!=0)

因为如果第一个 if 条件失败,它必然会随之而来。

最新更新