我正在使用外部库 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 条件失败,它必然会随之而来。