Jacobi法のプログラムをCで作成したのですが、、、

[yu]

Segmentaion faultがでて上手く動きません。Input,A,Vは出力されるのですが、Iterarionの部分が上手くいってないようです。
しかし、どこを修正すれば良いのかが分かりません。どなたか教えていただけないでしょうか。
Cは初めて3ヶ月くらいです。まだまだわからないことがたくさんあります。よろしくお願いします。

OS:Mac

エラー文
(Input)
1.0000E+00 1.0000E+00
1.0000E+00 1.0000E+00
(A)
1.0000E+00 1.0000E+00
1.0000E+00 1.0000E+00
(V)
1.0000E+00 0.0000E+00
0.0000E+00 1.0000E+00

----- Iterarion 1 -----
Segmentation fault




プログラムコード
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define MAXROW 20 /* 行数 */
#define MAXCOL 20 /* 列数 */
#define NUMITEMS 10 /* 画面に表示する際の1行の要素 */
#define NEARLY_ZERO 1.E-16 /* ゼロ */
#define ZERO_TOLERANCE 1.E-6 /* ゼロ */
#define vect_squared(n,x) vect_prod((n),(x),(x))
#define L2_norm(n,x) sqrt(vect_squared((n),(x)))
typedef struct
{
double R;
double I;
} complex_t;
complex_t determinant(int n, complex_t A[][MAXCOL]);
double vect_prod(int n, double *x, double *y);
void mat_vect(int m, int n, double A[][MAXCOL], double *x, double *b);
void vect_mat(int m, int n, double *x, double A[][MAXCOL], double *b);
void mat_prod(int m, int n, int l,
double A[][MAXCOL], double B[][MAXCOL], double C[][MAXCOL]);
void vect_copy(double *b, double *a, int n);
void mat_copy(int m, int n, double A[][MAXCOL], double B[][MAXCOL]);
void mat_read(int m, int n, double A[][MAXCOL]);
void vect_print(char *s, int n, double *x);
void mat_print(char *s, int m, int n, double A[][MAXCOL]);
double cleanup(int n, double A[][MAXCOL]);
#include "matrix.h"
/* ベクトルの内積 */
double vect_prod(int n, double *x, double *y)
{
int i;
double prod = 0.0;
for(i = 0; i < n; i++)
prod += x * y;
return prod;
}

/* 行列とベクトルの積 */
void mat_vect(int m, int n, double A[][MAXCOL], double *x, double *b)
{
int i, j;
for(i = 0; i < m; i++)
{
b = 0.0;
for(j = 0; j < n; j++)
b += A[j] * x[j];
}
}
/* ベクトルと行列の積 */
void vect_mat(int m, int n, double *x, double A[][MAXCOL], double *b)
{
int i, j;
for(j = 0; j < n; j++)
{
b[j] = 0.0;
for(i = 0; i < m; i++)
b[j] += x * A[j];
}
}

/* 行列と行列の積 */
void mat_prod(int m, int n, int l,
double A[][MAXCOL], double B[][MAXCOL], double C[][MAXCOL])
{
int i;
for(i = 0; i < m; i++)
vect_mat(n, n, &A[0], B, &C[0]) ;
}

/* ベクトルのコピー */
void vect_copy(double *b, double *a, int n)
{
int i;
for(i = 0; i < n; i++)
b = a[i];
}

/* 行列のコピー */
void mat_copy(int m, int n, double A[][MAXCOL], double B[][MAXCOL])
{
int i;
for(i = 0; i < m; i++)
vect_copy(&B[i][0], &A[i][0], n);
}

/* 行列の読み込み */
void mat_read(int m, int n, double A[][MAXCOL])
{
int i, j;
for(i = 0; i < m; i++)
for(j = 0; j < n; j++)
scanf("%lf", &A[i][j]);
}

/* ベクトルの表示 */
void vect_print(char *s, int n, double *x)
{
int j;
if(s) printf("(%s)\n", s);
for(j = 0; j < n; j++)
printf("%11.4E%s", x[j], (j+1)% NUMITEMS ? " " : "\n");
if(j % NUMITEMS)
printf("\n");
}

/* 行列の表示 */
void mat_print(char *s, int m, int n, double A[][MAXCOL])
{
int i;
printf("(%s)\n", s);
for(i = 0; i < m; i++)
vect_print(NULL, n, &A[i][0]);
}

/* 微小要素の消去 */
double cleanup(int n, double A[][MAXCOL])
{
int i, j;
double max_norm = 0.0;
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
if(fabs(A[i][j]) > max_norm) max_norm = fabs(A[i][j]);
if(max_norm > NEARLY_ZERO)
{
for(i = 0; i < n; i++)
for(j = 0; j < n; j++)
if(fabs(A[i][j]) / max_norm < ZERO_TOLERANCE)
A[i][j] = 0.0;
}
return max_norm;
}

//行列の固有値(古典的Jacob法)
#include "matrix.h"
#define MAXITR 200 /* 繰り返し回数の上限　*/

/* 回転行列 */
typedef struct
{
int p, q; /* ゼロにする非対角要素の位置　*/
double c; /* cos θ */
double s; /* sin θ */
double t; /* tan θ */
double t2; /* s/(c+1) */
}JROTATE_t;

/* Aの更新 */
void update_A(int n, double A[][MAXCOL], JROTATE_t *J)
{
int i;
double aip, aiq;
for(i = 0; i < n; i++)
{
aip = A[i][J->p];
aiq = A[i][J->q];
if(i != J->p && i != J->q)
{
A[i][J->p] = aip - J->s * (J->t2 * aip + aiq);
A[i][J->q] = aiq + J->s * (aip - J->t2 * aiq);
A[J->p][i] = A[i][J->p];
A[J->q][i] = A[i][J->q];
}
else if(i == J->p)
{
A[i][i] = aip - J->t * aiq;
A[i][J->q] = 0.0;
}
else
{
A[i][i] = aiq + J->t * aip;
A[i][J->p] = 0.0;
}
}
}

/* Vの更新 */
void update_V(int n, double V[][MAXCOL], JROTATE_t *J)
{
int i;
double vip, viq;
for(i = 0; i < n; i++)
{
vip = V[i][J->p];
viq = V[i][J->q];
V[i][J->p] = vip - J->s * (J->t2 * vip + viq);
V[i][J->q] = viq + J->s * (vip - J->t2 * viq);
}
}

/* 次にゼロにする非対角要素の選択(古典的ヤコビ法)*/
void select_nondiag(int n, double A[][MAXCOL], JROTATE_t *J)
{
if(J->p == n - 2) /* 最終行*/
{
J->p = 0;
J->q = 1;
}
else if(J->q == n - 1) /* 最終列*/
{
J->p++;
J->q = J->p + 1;
}
else
J->q++;
}

/* 回転行列を求める*/
void calc_Jmat(double A[][MAXCOL], JROTATE_t *J)
{
if(fabs(A[J->p][J->q]) < NEARLY_ZERO)
{
J->c = 1.0;
J->s = 0.0;
J->t2 = 0.0;
}
else
{
double phi;
phi = 0.5 * (A[J->q][J->q] - A[J->p][J->p]) / A[J->p][J->q];
if(fabs(phi) < 1.E10)
J->t = 1.0 / (fabs(phi) + sqrt(phi*phi + 1.0));
else
J->t = 0.5 / fabs(phi);
J->t = phi < 0.0 ? -J->t : J->t;
J->c = 1.0 / sqrt(1.0 + J->t * J->t);
J->s = J->t * J->c;
J->t2 = J->s / (1.0 + J->c);
}
}

/* 非対角要素の２乗和*/
double nondiag_norm(int n, double A[][MAXCOL])
{
int i, j;
double norm = 0.0;
for(i = 0; i < n; i++)
{
for(j = 0; j < n; j++)
{
if(j == i)
continue;
norm += A[i][j] * A[i][j];
}
}
return sqrt(norm);
}

/*
* ヤコビ法
* 実対称行列の固有値、固有ベクトルの計算
*　　(入力)
*　　　　n : 行列サイズ
*　　　　A[][] : nxn実対称行列
*　　(出力)
*　　　　lamda[] : Aの固有値
*　　　　A[][] : 固有ベクトルを並べたもの
*　　(戻り値)
*　　　　0 : 正常に終了
*　　　　1 : 繰り返し回数の上限に達した
*/

int Jacobi(int n, double A[][MAXCOL], double lambda[])
{
int i, j, iteration = 0, ret_code;
double V[MAXROW][MAXCOL];
JROTATE_t J;
for(i = 0; i < n; i++)
{
for(j = 0; j < n; j++)
V[i][j] = 0.0;
V[i][i] = 1.0;
}
mat_print("A", n, n, A);
mat_print("V", n, n, V);
while(1)
{
if(nondiag_norm(n, A) < ZERO_TOLERANCE)
{
ret_code = 0; break; /* 正常終了 */
}
if(iteration++ >= MAXITR)
{
ret_code = 1; break; /* 繰り返し回数の上限 */
}
printf("\n----- Iterarion %d -----\n", iteration);
select_nondiag(n, A, &J);
calc_Jmat(A, &J);
printf("(p,q) = (%d, %d) c = %f s = %f\n", J.p, J.q, J.c, J.s);
update_A(n, A, &J);
update_V(n, V, &J);
cleanup(n, A); cleanup(n, V);
mat_print("A", n, n, A);
mat_print("V", n, n, V);
printf("norm of nondiagonals = %f\n", nondiag_norm(n, A));
}
for(i = 0; i < n; i++)
lambda[i] = A[i][i];
mat_copy(n, n, V, A);
return ret_code;
}

/* メイン関数 */
int main(void)
{
int n, ret_status;
double A[MAXROW][MAXCOL], Aorg[MAXROW][MAXCOL];
double lambda[MAXCOL];
scanf("%d", &n);
mat_read(n, n, A);
mat_print("Input", n, n, A);
mat_copy(n, n, A, Aorg);
ret_status = Jacobi(n, A, lambda);
vect_print("eigen values", n, lambda);
printf("status = %d\n", ret_status);
return 0;
}

[さかまき]

printf("\n----- Iterarion %d -----\n", iteration);
が表示されて
printf("(p,q) = (%d, %d) c = %f s = %f\n", J.p, J.q, J.c, J.s);
が表示されていないから
select_nondiag(n, A, &J);　か　calc_Jmat(A, &J);　があやしいことになります。
ところどころに printf をいれて、どこで発生しているかを絞り込んでいきましょう。

#include "matrix.h"　があるとコンパイルできないのでコメントアウトしたらコンパイルできました。
これは不要でいいのですか？
とりあえずwindowsでは動きました。（入力値は適当ですが）

それとソース読みにくいです。
※コードを貼り付ける場合は 〔ＣＯＤＥ〕と 〔／ＣＯＤＥ〕 で囲って下さい。でお願いします

[みけCAT]

>※コードを貼り付ける場合は 〔ＣＯＤＥ〕と 〔／ＣＯＤＥ〕 で囲って下さい。でお願いします
念のため補足します。
実際は半角の[code]～[/code]で囲ってください。
蛇足だったらすみません。