关于 segmentation fault
Hi 各位大牛，我第一次来论坛，我在写一个Householder_arnoldi算法，用来找到krylov space的基，但我对C语言知道的不多，尤其是指针和内存操作，所以特把代码贴上，希望有C语言大牛帮我一下，感激不尽！
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
int main(
int argc,
char *argv[])
{
FILE *infileptr;
int dim1, dim2,i, j,k;
double **matA;
infileptr = fopen(argv[1],"r");
if(infileptr == NULL)
{
printf("FILE %s does not exist, abort\n", argv[1]);
printf("Usage: exe input_matrix_file_name\n");
exit(-1);
}
fread(&dim1,sizeof(int), 1, infileptr);
fread(&dim2,sizeof(int), 1, infileptr);
printf("\n\nm = %d, n = %d\n", dim1, dim2);
int M=dim1,N=dim2;
matA = (double **)malloc(sizeof(double *)*M);
matA[0] = (double *)malloc(sizeof(double)*M*N);
for(i = 1; i < M; i++)
{
matA[i] = matA[i-1] + N;
}
///// The code below is to read in the matrix data from the data file
// Matrix data is saved in matA
for(i = 0; i < M; i++)
{
fread(matA[i],sizeof(double), N,infileptr);
printf("ROW %d: ", i);
for(j =0; j < N; j++)
{
printf("%g ", matA[i][j]);
}
printf("\n");
}
fclose(infileptr);
// Implement your alg. here
double householder(double *, int , double **, int);
void reflector(double *, double **, int);
double **z;
double **w;
double **h;
double **v;
double **I;
double m=5; // krylov space dimension
z = (double **)malloc(sizeof(double *)*(m+1));
z[0] = (double *)malloc(sizeof(double)*(m+1)*M);
w = (double **)malloc(sizeof(double *)*(m+1));
w[0] = (double *)malloc(sizeof(double)*(m+1)*M);
h = (double **)malloc(sizeof(double *)*(m+1));
h[0] = (double *)malloc(sizeof(double)*(m+1)*M);
v = (double **)malloc(sizeof(double *)*(m+1));
v[0] = (double *)malloc(sizeof(double)*(m+1)*M);
I = (double **)malloc(sizeof(double *)*M);
I[0] = (double *)malloc(sizeof(double)*M*M);
for(i = 1; i < (m+1); i++)
{
z[i] = z[i-1] + M;
w[i] = w[i-1] + M;
h[i] = h[i-1] + M;
v[i] = v[i-1] + M;
}
for(i=0; i<(m+1);i++)
{
for (j=0;j<M;j++)
{ z[i][j]=0;
w[i][j]=0;
h[i][j]=0;
v[i][j]=0;
}
}
z[0][0]=1;//z1=e1
// generate identity matrix
for (i=1; i<M; i++)
{h[i]=h[i-1]+M;}
for (i=0; i<M; i++)
{
for (j=0;j<M;j++)
{ I[i][j]=0;}
I[i][i]=1;
}
for (j=0;j<(m+1);j++)
{
h[j][j]= householder(z[j], M-j, &w[j],M); // try to find wj for zj
for (i=0;i<j;i++)
{h[j][i]=z[j][i];}
v[j]=I[j]; // compute vj=P1*P2*...*Pj*ej
for (k=j-1;k>=0;k--)
{reflector(w[k],&v[j],M);}
if (j<=m)
{for(i=0;i<M;i++) //compute zj+1=Pj*Pj-1*...*P1*A*vj
{
for(k=0;k<M;k++)
{ z[j+1][i] += matA[i][k]*v[j][k];}
}
for (k=0;k<j;k++)
{reflector(w[k],&z[j+1],M);}
}
}
for(j=0; j < (m+1); j++)
{
for(i =0; i < M; i++)
{
printf("%g ", v[j][i]);
}
printf("\n");
}
free(z);
free(w);
free(h);
free(v);
free(I);
return 1;
}
double householder(double *x, int k, double **u, int n)
{
double *v;
double norm;
double alpha =0;
double xmax;
int i,j;
// compute xmax
xmax=x[0];
for (i=n-k;i<n;i++)
{ if (x[i]<x[i+1])
xmax=x[i+1];
}
// compute vj and norm of vector v
for (j=n-k;j<n;j++)
{ v[j]=x[j]/xmax;
norm=norm+v[j]*v[j];
}
norm=sqrt(norm);
// compute alpha and v[1]
if (v[n-k]>0)
{v[n-k]=v[n-k]+norm;
alpha=-norm*xmax;}
else
{v[n-k]=v[n-k]-norm;
alpha=norm*xmax;}
// compute u(w)
for (i=n-k;i<n;i++)
{*u[i]=v[i]/norm;}
return alpha;
}
void reflector(double *u, double **x, int n)
{
int i,j;
double r=0;
for (i=0;i<n;i++)
{r=r+2*u[i]*(*x[i]);}
for (j=0;j<n;j++)
{*x[j]=*x[j]-u[j]*r;}
}
