コアダンプする箇所がわからない。

Segmentation fault (コアダンプ)

//グラフェンナノリボンエネルギーバンド　相互作用 反強磁性
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>

#define M  10	// 存在するサイト数 
#define N  20	//サイト数×２　A+Bの個数
#define P  100	//リボン長　全電子数　P*L
#define Q  1 //回数

double b_pi,N_c;
double complex t_1, t_2, t_3, t_4, t_5, t_6, t_7, t_8, T_1, T_2;	//t_方向
double k_x[P+1];
double n_2[2*N];		//代入する値
double meow_updw;		//up&dwのフェルミ

double Energy[P+1][2*N];
double complex u[P+1][2*N][2*N];	//固有ベクトル
double complex fermi_1[P+1][2*N];
double gap[P+1];       //Eg,フェルミエネルギーバンド

double u_5_1[P+1][2*N][2*N],u_6_1[P+1][2*N],u_7_1[2*N],u_8_1[2*N]; 

void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );

int main(void)
{
	FILE* fp;
	FILE* fp_1;
	FILE* fp_2;

	if ((fp = fopen("ziginteneantiferromagnetism M=4.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}
	if ((fp_1 = fopen("zigintantiferromagnetism.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}
	if ((fp_2 = fopen("zigintantiferromagnetism_1.txt", "w")) == NULL)
	{
		printf("Cannot open the file\n");
		exit(1);
	}

	b_pi=M_PI;
	N_c=(P+1)*N;
    double k_xx = -b_pi; 
	double kt=0.00001*8.6171*pow(10,-5);		//-273℃＊ボルツマン定数(eV/K)
	double U=4.32;		//相互作用エネルギーの平均(eV)
	double t=2.75;		//飛び移り積分(eV)
	double ribbon_length;		//リボン幅
	double Energy_all[(P+1)*2*N];		//全エネルギー

	//刻み幅
    for(int i=0;i<=P;i++){
        k_x[i]=k_xx+2.0*b_pi/P*i;
    }
	//リボン幅
	ribbon_length=0.246/pow(3,0.5)/2*(3*M-2);		//0.0710140

	//初期条件
	for(int i=0;i<2*N;i++){
		n_2[i]=0.4; 
	} 

	int z=0;
    while(z<Q){

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

        t_4 = 1;
	    t_5 = cexp(-1 * I * -1/ 2 * k_x[r]);			
        t_6 = cexp(-1 * I * 1/2 * k_x[r]);		
        t_1 =1;
    	t_2 = cexp(I * -1 / 2 * k_x[r]);
	    t_3 = cexp(I * 1 / 2 * k_x[r]);
    	T_1 = t_5 + t_6;
     	T_2 = t_2 + t_3;

		for(int i=0;i<2*N*2*N;i++){
			A[i]=0;
		}
		for(int i=0;i<2*N;i++){
			if(i<N){
				A[i*2*N+N+i]=U/t*n_2[i];
			}
			else{
				A[2*N*(i+N)+i]=U/t*n_2[i];
			}
		}
		for(int j=1;j<M-1;j++){
			A[2*N*j+M+j-1]=t_4;   A[j*2*N+M+j]=T_1;
			A[(j+M)*2*N+j]=T_2; A[(j+M)*2*N+j+1]=t_1;

			A[2*N*(j+N)+N+M+j-1]=t_4;   A[2*N*(j+N)+N+M+j]=T_1;
			A[(j+N+M)*2*N+N+j]=T_2; A[(j+N+M)*2*N+N+j+1]=t_1;
		}
		A[M]=T_1;
		A[2*N*(M-1)+N-2]=t_4; A[2*N*(M-1)+N-1]=T_1;
		A[2*N*(N+1)-M]=T_1; 
		A[2*N*(N+M)-2]=t_4; A[2*N*(M+N)-1]=T_1;

		A[2*N*M]=T_2;  A[2*N*M+1]=t_1; 
		A[2*N*(N-1)+M-1]=T_2;
		A[2*N*(N+M)+N]=T_2; A[2*N*(N+M)+N+1]=t_1; 
		A[2*N*2*N-M-1]=T_2;

		/*for(int i=0;i<2*N*2*N;i++){	
			fprintf(fp_2,"%2f\n",A[i]);	
		}*/

  		double w[2*N]              	;// 実固有値が小さい順に入る
	
 		int    lda  = 2*N	            ;// 対角化する正方行列のサイズ
  		double _Complex work[6*2*N]	;// 対角化する際に使用するメモリ
  		int    lwork = 6*2*N		    ;// workの次元
  		double rwork[3*2*N-2]      	;// 3*SIZE-2で固定
  		int    info                 ;// 成功すれば0、失敗すれば0以外を返す

  		zheev_( &jobz, &uplo, &n, A, &lda, w, work, &lwork, rwork, &info );

    	for(int i=0;i<2*N;i++){
			for(int j=0;j<2*N;j++){
				u[r][i][j]=0.0;
				u_5_1[r][i][j]=0.0; 
			}
			u_6_1[r][i]=0.0;
			u_7_1[i]=0.0;
			u_8_1[i]=0.0; 
		}
		for(int j=0;j<2*N;j++){
			Energy[r][j]=w[j];
		}
		for(int i=0;i<2*N;i++){
			for(int j=0;j<2*N;j++){
				u[r][i][j]=A[i*2*N+j];
			}	
		}
	}
	//フェルミエネルギー up&dw
	for(int i=0;i<=P;i++){
			for(int j=0;j<2*N;j++){
				Energy_all[i*2*N+j]=Energy[i][j];
			}
		}
		double tmp_ene;
		for (int i=0; i<=P*2*N; i++) {
   		 	for (int j=i+1; j<=P*2*N; j++) {
				if (Energy_all[i] > Energy_all[j]) {
	    			tmp_ene =  Energy_all[i];
    	    		Energy_all[i] = Energy_all[j];
        			Energy_all[j] = tmp_ene;
				}
			}
    	}
		meow_updw=(Energy_all[(P+1)*N]+Energy_all[(P+1)*N-2])/2;
		/*for(int i=0;i<=P*N*2;i++){
			printf("%f\n",Energy_all[i]);
		}
		printf("%f\n",meow_updw);*/
			
	//Eg
	double tmp;
	for (int i=0; i<=P; i++) {
		gap[i]=Energy[i][N]-Energy[i][N-1];
	}
 	for (int i=0; i<=P; ++i) {
    	for (int j=i+1; j<=P; ++j) {
      		if (gap[i] > gap[j]) {
        	tmp =  gap[i];
        	gap[i] = gap[j];
        	gap[j] = tmp;
    	  	}
    	}
  	}
	for(int l=0;l<=P;l++){
		for(int i=0;i<2*N;i++){
			fermi_1[l][i]=1/(1+exp((Energy[l][i]-meow_updw)/kt)); 
		}
	}
	for(int i=0;i<2*N;i++){
		for(int l=0;l<=P;l++){
			for(int j=0;j<2*N;j++){
				u_5_1[l][i][j]=pow(cabs(u[l][j][i]),2)*fermi_1[l][j];
				u_6_1[l][i]+=u_5_1[l][i][j];
			}
			u_7_1[i]+=u_6_1[l][i]; 
		}
		u_8_1[i]=u_7_1[i]/(P+1);
	}	
	for(int i=0;i<2*N;i++){
		n_2[i]=u_8_1[i];
	}

    z++;
    }

	for(int i=0;i<=P;i++){
    	for(int j=0;j<2*N;j++){
			fprintf(fp, "%2f %lf %lf %lf\n",k_x[i],Energy[i][j]-meow_updw,Energy[i][N]-meow_updw,Energy[i][N-1]-meow_updw);
			// fprintf(fp, "%2f %lf %lf\n",k_x[i],Energy[i][N+1]-meow_updw,Energy[i][N-2]-meow_updw);
		}   
	}
	for(int i=0;i<2*N;i++){	
		fprintf(fp_2,"%2f\n",n_2[i]);	
	}
	for(int i=0;i<M;i++){	
		fprintf(fp_1,"%2f\n",n_2[i]);
		fprintf(fp_1,"%2f\n",n_2[i+M]);
	}
	printf("リボン長：%2f\n", ribbon_length);
    printf("最小ギャップ：%f\n",gap[0]/(U/t));
    
    
 	fclose(fp);
	fclose(fp_1);
	fclose(fp_2);
	return 0;
}