コアダンプする原因がわからない

Segmentation fault (コアダンプ)

//グラフェンナノリボンエネルギーバンド　spinwave susceptibility
#define _USE_MATH_DEFINES
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>

#define M  2	// 存在するサイト数   //M=24で5nmくらい
#define N  4	//サイト数×２　A+Bの個数
#define P  101	//リボン長　全電子数　P*
#define R  10

void zheev_( char*, char*, int*, double _Complex*,  int*, double*, double _Complex*, int*, double*, int* );
double f(double x, int j);

double u_up[P+1][N][N]={0.0}, u_dw[P+1][N][N]={0.0};		//固有ベクトル
double u_up_1[P+1][N][N]={0.0}, u_dw_1[P+1][N][N]={0.0};
double Energy_dw_1[P+1][N]={0.0},Energy_up_1[P+1][N]={0.0};  

int main(int argc, char** argv){ 
	FILE* fp;	
	FILE* fp_1;
    FILE* fp_2;
	FILE* fp_3;
	FILE* fp_4;

    double b_pi=M_PI;
	double U=1.0;		//相互作用エネルギーの平均(eV) u/t

    if ((fp = fopen("kai.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}
	fp_1 = fopen("u_dw.txt", "r");   
    if(fp_1 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				fscanf(fp_1, "%lf", &(u_dw[l][i][j]) );   	
			}
		} 
	}
	fp_2 = fopen("Energydw_1.txt", "r");   
    if(fp_2 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			fscanf(fp_2, "%lf", &(Energy_dw_1[l][i]) );   	
		} 
	}
    fp_3 = fopen("Energyup_1.txt", "r");   
    if(fp_3 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			fscanf(fp_3, "%lf", &(Energy_up_1[l][i]) );   	
		} 
	}
	fp_4 = fopen("u_up.txt", "r");   
    if(fp_4 == NULL) {
        printf("ファイルを開くことが出来ませんでした．¥n");
    }
    for(int l=0; l<=P; l++){
		for(int i=0;i<N;i++){
			for(int j=0;j<N;j++){
				fscanf(fp_4, "%lf", &(u_up[l][i][j]) );   	
			}
		} 
    }

	double x_start1=1.8;
	double x_start2=0.0;
	/*double c=0.0;
		
	for(int j=0;j<(P+1)/2;j++){
	    for(int i=0;i<R;i++){
			c=x_start1/R*i; 
			double fc=0.0;
			fc=f(c,j)/(P+1)-U;
			// printf("x=%lf,q=%lf ,%lf\n",c,2*b_pi/(P+1)*j,fc);
			fprintf(fp, "%lf %lf %lf\n", c,2.0*b_pi/(P+1)*j,fc);
		}
	}*/
	double err = 100.0;
  	double a = x_start1;
  	double b = x_start2;
  	double c;

	for(int j=0;j<(P+1)/2;j++){
	  	while(err>0.00000001){
    		c = (a+b)/2.0;
      		double fc = f(c,j)/(P+1)-U;
      		err = fc*fc;
      		if(fc > 0){
      			b = c;
    	  	}
    	  	else{
    	    	a = c;
    	  	}
    	  	printf("%lf\n",c);
			fprintf(fp, "%lf %lf %lf\n", c,2.0*b_pi/(P+1)*j,fc);
 		}
	}

    fclose(fp);	
 	fclose(fp_1);
    fclose(fp_2);
 	fclose(fp_3);	
	fclose(fp_4);
    return 0;
}

double f(double x, int j){  	
    double kai[N][N] = {0.0};
    double kai_2[N][N][P + 1][M] = {0.0};
    double kai_1[N][N][P + 1][M][M] = {0.0};
    double kai_1_dw[N][N][P + 1][M][M] = {0.0}, kai_1_up[N][N][P + 1][M][M] = {0.0};
    double kai_3[N][N][P + 1] = {0.0};
    double B[N] = {0.0};
	double Omega_up[P+1][M][M] = {0.0}, Omega_dw[P+1][M][M] = {0.0};
	double U_up_plus[N][P+1][M][M] = {0.0}, U_up_minus[N][P+1][M][M] = {0.0};
	double U_dw_plus[N][P+1][M][M] = {0.0}, U_dw_minus[N][P+1][M][M] = {0.0};	
	
	for(int v=0;v<=P;v++){	//wavenumber
		for (int t=0;t<M;t++){	//valence
			for(int s=M;s<N;s++){	//conductive
					Omega_up[v][t][s]=x-Energy_dw_1[v][s]+Energy_up_1[(v+j)%(P+1)][t];
					Omega_dw[v][t][s]=x+Energy_up_1[v][s]-Energy_dw_1[(v+j)%(P+1)][t];
			}		
		}
	}
	for(int k=0;k<N;k++){
		for(int h=0;h<N;h++){
			for(int v=0;v<=P;v++){	//wavenumber
				for (int t=0;t<M;t++){	//valence
					for(int s=M;s<N;s++){	//conductive
						U_up_plus[k][v][t][s]=conj(u_up[v][t][k])*u_dw[(v+j)%(P+1)][s][k];
						U_up_minus[h][v][t][s]=conj(u_dw[(v+j)%(P+1)][s][h])*u_up[v][t][h];
						U_dw_plus[k][v][t][s]=conj(u_dw[v][t][k])*u_up[(v+j)%(P+1)][s][k];
						U_dw_minus[h][v][t][s]=conj(u_up[(v+j)%(P+1)][s][h])*u_dw[v][t][h];
					}
				}
			}
		}
	}
	for(int k=0;k<N;k++){
		for(int h=0;h<N;h++){
			for(int v=0;v<=P;v++){	//wavenumber
				for (int t=0;t<M;t++){	//valence
					for(int s=M;s<N;s++){	//conductive
						kai_1_up[k][h][v][t][s]=U_up_plus[k][v][t][s]*U_up_minus[h][v][t][s]/Omega_up[v][t][s];
						kai_1_dw[k][h][v][t][s]=U_dw_plus[k][v][t][s]*U_dw_minus[h][v][t][s]/Omega_dw[v][t][s];
						kai_1[k][h][v][t][s]=-kai_1_up[k][h][v][t][s]+kai_1_dw[k][h][v][t][s];

						kai_2[k][h][v][t]+=kai_1[k][h][v][t][s];
						// printf("%2f\n",kai_1[k][h][v][t][s]);
					}
					kai_3[k][h][v]+=kai_2[k][h][v][t];
					// printf("%2f \n",kai_2[k][h][v][t]);
				}	
				kai[k][h]+=kai_3[k][h][v];
				// printf("%2f \n",kai_3[k][h][v]);
			}			
				// printf("%2f \n",kai[k][h]);
		}
	}								

    //感受率の対角化
                    // 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
                    char   jobz = 'V'           ;// 固有ベクトルを計算する
                    char   uplo = 'U'           ;// Aに上三角行列を格納
                    int    n    = N             ;// 対角化する正方行列のサイズ
                    double _Complex A[N*N]		;// 対角化する行列。対角化後は固有ベクトルが並ぶ

					for(int k=0;k<N;k++){
						for(int h=0;h<N;h++){
		                    A[N*k+h]=kai[k][h];
						}
					}
					// fprintf(fp_1, "%2f %2f \n",A[N*k+h],kai[i][j][k][h]);
					// printf("%2f %2f\n",A[N*h+k],kai[h][k]);

                    double w[N]              	;// 実固有値が小さい順に入る
            
                    int    lda  = N	            ;// 対角化する正方行列のサイズ
                    double _Complex work[6*N]	;// 対角化する際に使用するメモリ
                    int    lwork = 6*N		    ;// workの次元
                    double rwork[3*N-2]      	;// 3*SIZE-2で固定
                    int    info                 ;// 成功すれば0、失敗すれば0以外を返す
  		            
                    zheev_( &jobz, &uplo, &n, A, &lda, w, work, &lwork, rwork, &info );
                        
                    for(int m=0;m<N;m++){
		            	B[m]=w[m];
					// printf("%2f\n",B[m]);
		            }

                double kai_acc=0.0;
				kai_acc=B[N-1];

  return kai_acc;
}