コアダンプする原因がわからない
Q&A
Closed
解決したいこと
コアダンプします。計算している変数をすべて調べたところkai_1_upとmainのなかのｊのループを切るとコアダンプしないことがわかります。
しかし、配列上では特に問題なく定義しているようにもおもえ。gdbの方で確認しても目立った問題の方は見られませんでした。
前の質問の方で初期化し忘れているということを多くアドバイスをされ、初期化についてはいろいろ気を使いました。
発生している問題・エラー
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;
}
自分で試したこと
printf やらデバッグやらはいろいろしましたが全く原因がわかりませんでした。
正直根本的なプログラム知識がないので、触ってはいけないメモリ領域に入っているのかもしれませんが、ご容赦ください。
環境についてはよくわかりませんが、
Ubuntu 20.04.5 LTS
でコンパイルには icc -qmkl -mcmodel=medium -shared-intel
を用いています。
0
