segmentation fault してしまう。
解決したいこと
プログラムを動かすと、セグメントしてしまう。エラーをデバッグで見るとファイル（fp）に例外があると出るが具体的にどう対処すべきかわからない。
発生している問題・エラー
Segmentation fault (コアダンプ)
例外が発生しました
Segmentation fault
上がcmdでの表示。
下がVSでのデバッグでの表示
該当するソースコード
//グラフェンナノリボンエネルギーバンド 相互作用 強磁性
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <stdlib.h>
#define M 6 // 存在するサイト数 //M=24で5nmくらい
#define N 12 //サイト数×２ A+Bの個数
#define P 101 //リボン長 全電子数 P*L
#define Q 200 //回数
#define R 200
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 e[P+1],k_x[P+1];
double n_up[N],n_dw[N]; //代入する値
double n_1_up[M],n_1_dw[M],n_updw;
double meow_updw; //up&dwのフェルミ
double Energy_dw[P+1][N],Energy_up[P+1][N];
double Energy_sum,Energy_sum_1;
double complex u_up[P+1][N][N], u_dw[P+1][N][N]; //固有ベクトル
double fermi_up_1[P+1][N],fermi_dw_1[P+1][N];
double gap[P+1],gap_up[P+1],gap_dw[P+1]; //Eg,フェルミエネルギーバンド
double complex u_5_up_1[P+1][N][N],u_6_up_1[P+1][N],u_7_up_1[N],u_8_up_1[N];
double complex u_5_dw_1[P+1][N][N],u_6_dw_1[P+1][N],u_7_dw_1[N],u_8_dw_1[N];
double DOS[P+1][P+1][N];
double ddd[P+1][P+1],dd[P+1];
void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );
int main(void)
{
FILE* fp;
if ((fp = fopen("kai.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 ene=-4.0;
double kt=0.00001*8.6171*pow(10,-5); //-273℃＊ボルツマン定数(eV/K)
// double kt=300*8.6171*pow(10,-5); //27℃＊ボルツマン定数(eV/K)
double U=1.0; //相互作用エネルギーの平均(eV) u/t
double t=1; //2.75; //飛び移り積分(eV)
double delta_E=2.75*pow(10,-2); //エネルギー分解能
double ribbon_length; //リボン幅
double Energy_all[(P+1)*2*N*2]; //全エネルギー
double meu_B= 5.788*pow(10,-5); //ボーア磁子 [/eV/T]
//刻み幅
for(int i=0;i<=P;i++){
e[i]=ene+8.0/P*i;
k_x[i]=k_xx+2.0*b_pi/(P+1)*i;
}
//リボン幅
ribbon_length=0.246/pow(3,0.5)/2*(3*M-2); //0.0710140
//初期条件
//端が違う
n_up[0]=0.01;
n_up[1]=0.01;
n_up[M]=0.6;
n_up[M+1]=0.6;
n_up[N-2]=0.99;
n_up[N-1]=0.99;
n_dw[0]=0.99;
n_dw[1]=0.99;
n_dw[M]=0.4;
n_dw[M+1]=0.4;
n_dw[N-2]=0.01;
n_dw[N-1]=0.01;
int z=0;
while(z<Q){
//行列計算
//dw
for(int r=0;r<=P;r++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = N ;// 対角化する正方行列のサイズ
double _Complex A[N*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<N*N;i++){
A[i]=0;
}
for(int i=0;i<N;i++){
A[i*N+i]=U/t*n_up[i];
}
for(int j=1;j<M-1;j++){
A[j*N+M+j]=t_4; A[j*N+M+j+1]=T_1;
A[(j+M)*N+j]=T_2; A[(j+M)*N+j+1]=t_1;
}
A[M]=T_1; A[M+1]=0; A[N*M-2]=t_4; A[N*M-1]=T_1;
A[N*M]=T_2; A[N*M+1]=t_1; A[N*N-1-M]=T_2; A[N*N-M]=0;
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 i=0;i<N;i++){
for(int j=0;j<N;j++){
u_dw[r][i][j]=0.0;
u_5_dw_1[r][i][j]=0.0;
}
u_6_dw_1[r][i]=0.0;
u_7_dw_1[i]=0.0;
u_8_dw_1[i]=0.0;
}
for(int j=0;j<N;j++){
Energy_dw[r][j]=w[j];
}
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u_dw[r][i][j]=A[i*N+j];
}
}
}
//up
for(int r=0;r<=P;r++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = N ;// 対角化する正方行列のサイズ
double _Complex A[N*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<N*N;i++){
A[i]=0;
}
for(int i=0;i<N;i++){
A[i*N+i]=U/t*n_dw[i];
}
for(int j=1;j<M-1;j++){
A[j*N+M+j]=t_4; A[j*N+M+j+1]=T_1;
A[(j+M)*N+j]=T_2; A[(j+M)*N+j+1]=t_1;
}
A[M]=T_1; A[M+1]=0; A[N*M-2]=t_4; A[N*M-1]=T_1;
A[N*M]=T_2; A[N*M+1]=t_1; A[N*N-1-M]=T_2; A[N*N-M]=0;
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 i=0;i<N;i++){
for(int j=0;j<N;j++){
u_up[r][i][j]=0.0;
u_5_up_1[r][i][j]=0.0;
}
u_6_up_1[r][i]=0.0;
u_7_up_1[i]=0.0;
u_8_up_1[i]=0.0;
}
for(int j=0;j<N;j++){
Energy_up[r][j]=w[j];
}
for(int i=0;i<N;i++){
for(int j=0;j<N;j++){
u_up[r][i][j]=A[i*N+j];
}
}
}
//フェルミエネルギー up&dw
for(int i=0;i<=P;i++){
for(int j=0;j<N;j++){
Energy_all[i*N+j]=Energy_dw[i][j];
Energy_all[i*N+j+(P+1)*N]=Energy_up[i][j];
}
}
double tmp_ene;
for (int i=0; i<(P+1)*N*2; i++) {
for (int j=i+1; j<(P+1)*N*2; 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-1])/2;
/*for(int i=0;i<(P+1)*N*2;i++){
printf("%f\n",Energy_all[i]);
}
printf("%f\n",meow_updw);*/
for(int l=0;l<=P;l++){
for(int i=0;i<N;i++){
fermi_dw_1[l][i]=1/(1+exp((Energy_dw[l][i]-meow_updw)/kt));
fermi_up_1[l][i]=1/(1+exp((Energy_up[l][i]-meow_updw)/kt));
}
}
for(int i=0;i<N;i++){ //site
for(int l=0;l<=P;l++){ //wavenumber
for(int j=0;j<N;j++){ //state?
u_5_dw_1[l][i][j]=pow(cabs(u_dw[l][j][i]),2)*fermi_dw_1[l][j];
u_6_dw_1[l][i]+=u_5_dw_1[l][i][j];
u_5_up_1[l][i][j]=pow(cabs(u_up[l][j][i]),2)*fermi_up_1[l][j];
u_6_up_1[l][i]+=u_5_up_1[l][i][j];
}
u_7_dw_1[i]+=u_6_dw_1[l][i];
u_7_up_1[i]+=u_6_up_1[l][i];
}
u_8_dw_1[i]=u_7_dw_1[i]/(P+1);
u_8_up_1[i]=u_7_up_1[i]/(P+1);
}
for(int i=0;i<N;i++){
n_up[i]=u_8_up_1[i];
n_dw[i]=u_8_dw_1[i];
}
z++;
}
//Eg
double tmp,tmp_dw, tmp_up;
for (int i=0; i<=P; i++) {
gap_up[i]=Energy_up[i][M]-Energy_up[i][M-1];
gap_dw[i]=Energy_dw[i][M]-Energy_dw[i][M-1];
gap[i]=Energy_up[i][M]-Energy_dw[i][M-1];
// gap[i]=Energy_up[i][M]-Energy_up[i][M-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 i=0;i<N;i++){
n_updw+=n_up[i]*n_dw[i];
}
//状態密度
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
for(int k=0;k<N;k++){
DOS[i][j][k]=1/pow(b_pi*2,0.5)/delta_E/b_pi*exp(-1*pow(e[i]-Energy_dw[j][k]+meow_updw,2)/2/delta_E/delta_E);
// printf("%f %f %f %f\n",creal(u_up[i][j][k]),cimag(u_up[i][j][k]),creal(u_dw[i][j][k]),cimag(u_dw[i][j][k]));
}
}
}
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
for(int k=0;k<N;k++){
ddd[i][j]+=DOS[i][j][k];
}
}
}
for(int i=0;i<=P;i++){
for(int j=0;j<=P;j++){
dd[i]+=ddd[i][j];
}
}
for(int i=0;i<(P+1)*N*2;i++){
Energy_sum+=Energy_all[i];
}
double omega[R]; int q[(P+1)]; //刻み
double kai[R][(P+1)], kai_1[R][(P+1)][P+1][M][M],kai_2[R][(P+1)][P+1][M],kai_3[R][(P+1)][P+1]; //複素感受率
for(int i=0;i<R;i++){
omega[i]=2.0/R*i;
}
for(int i=0;i<=P;i++){
q[i]=i;
}
for(int i=0;i<R;i++){ //omega
for(int s=0;s<=P;s++){ //q
for(int j=0;j<=P;j++){ //k
for(int k=0;k<M;k++){ //valence
for(int h=M;h<N;h++){ //conductive
kai_1[i][s][j][k][h]=-(conj(u_6_up_1[(j+q[s])%(P+1)][k])*u_6_dw_1[j][h]*conj(u_6_dw_1[(j-q[s]+P+1)%(P+1)][h])*u_6_up_1[j][k])/(omega[i]+Energy_up[(j+q[s])%(P+1)][h]-Energy_dw[j][k])+(conj(u_6_dw_1[(j+q[s])%(P+1)][k])*u_6_up_1[j][h]*conj(u_6_up_1[(j-q[s]+P+1)%(P+1)][h])*u_6_dw_1[j][k])/(omega[i]+Energy_dw[(j+q[s])%(P+1)][h]-Energy_up[j][k]);
kai_2[i][s][j][k]+=kai_1[i][s][j][k][h];
}
kai_3[i][s][j]+=kai_2[i][s][j][k];
}
kai[i][s]+=kai_3[i][s][j];
}
}
}
for(int i=0;i<R;i++){
for(int j=0;j<=P;j++){
fprintf(fp, "%2f %2f %2f\n",omega[i],-k_xx+2.0*b_pi/(P+1)*q[j],kai[i][j]/(P+1)/N);
}
}
printf("リボン長：%lf\n", ribbon_length);
printf("最小ギャップ：%lf\n",gap[0]);
printf("All Energy sum：%f /N_c:%lf \n",Energy_sum-U*n_updw*(P+1),(Energy_sum-U*n_updw*(P+1))/N_c);
printf("fermi Energy：%f\n",meow_updw);
fclose(fp);
return 0;
}
自分で試したこと
下のomega,kai.qの配列を用いて計算しているのですが、ここの行列サイズによるエラーなのかと思いましたが、for文の回数などを確認しましたが、あまりうまく行かなかったです。
環境はlinux cmdでコンパイルしています。
0
