コアダンプの原因がわかりません。ご教授お願いします。
解決したいこと
現在大学でプログラミングを始めたのですが、セグメント違反の原因がわからないため質問しました。
テーマは物理で、線形ハミルトニアンの計算を行っています。
発生している問題・エラー
Segmentation fault
配列内に過剰に要素を入れたり、大きすぎる数を計算したりなどは行っていない(はず)なのですが、セグメしてしまいます。
printfでいろいろ確認したところ、matrix_complex()とhamiltonin_element_complex()に原因があることは分かったのですが、それ以上がわかりません。
gdbを使ったのですが、
Program received signal SIGSEGV, Segmentation fault.
0x0000000000408bdc in hamiltonin_element_complex ()
(gdb) bt
#0 0x0000000000408bdc in hamiltonin_element_complex ()
#1 0x3faed73b855cd58a in ?? ()
#2 0x0000000000000000 in ?? ()
(gdb)
ということだけがわかりました。
該当するソースコード
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <stdlib.h>
#define N 6
#define L 2*N
int ket_index_up = 0, ket_index_down = 0;
int NUM_KET_up = 0, NUM_KET_down = 0, NUM_KET_creation_up = 0, NUM_KET_creation_down = 0, NUM_KET_annihilation_up = 0, NUM_KET_annihilation_down = 0;
int size = 0, size_creation_up = 0, size_creation_down = 0, size_annihilation_up = 0, size_annihilation_down = 0;
double t_opt = 1.0;
void shokika_double(int length, double hairetu[length]){
for(int i = 0; i < length; i++){
hairetu[i] = 0.0;
}
}
void shokika_complex(int length, double complex hairetu[length]){
for(int i = 0; i < length; i++){
hairetu[i] = 0.0 + 0.0 * I;
}
}
void shokika_int(int length, int hairetu[length]){
for(int i = 0; i < length; i++){
hairetu[i] = 0.0;
}
}
void shokika_matrix(int length, double hairetu[length][length]){
for(int i = 0; i < length; i++){
for(int j = 0; j < length; j++){
hairetu[i][j] = 0.0;
}
}
}
void shokika_matrix_complex(int length, double complex hairetu[length][length]){
for(int i = 0; i < length; i++){
for(int j = 0; j < length; j++){
hairetu[i][j] = 0.0;
}
}
}
void shokika_vector_convert(int length, double eigenvector_spin[N][length]){
for(int i = 0; i < N; i++){
for(int j = 0; j < length; j++){
eigenvector_spin[i][j] = 0.0;
}
}
}
void shokika_renbunsuu_complex(double complex alpha[N][N][101], double beta[N][N][100]){
for(int i = 0; i < N; i++){
for(int j = 0; j < N; j++){
shokika_complex(101, alpha[i][j]); shokika_double(100, beta[i][j]);
}
}
}
void copy_double(int length, double hairetu_before[length], double hairetu_after[length]){
for(int i = 0; i < length; i++){
hairetu_after[i] = hairetu_before[i];
}
}
void copy_double_complex(int length, double complex hairetu_before[length], double complex hairetu_after[length]){
for(int i = 0; i < length; i++){
hairetu_after[i] = hairetu_before[i];
}
}
void copy_matrix_double_complex(int length, double complex hairetu_before[length][length], double complex hairetu_after[length][length]){
for(int i = 0; i < length; i++){
for(int j = 0; j < length; j++){
hairetu_after[i][j] = hairetu_before[i][j];
}
}
}
void copy_int(int length, int hairetu_before[length], int hairetu_after[length]){
for(int i = 0; i < length; i++){
hairetu_after[i] = hairetu_before[i];
}
}
void hyouji(int matrix_size, int search_size[matrix_size], int decimal, int sign, int index, int is){
// (is == 1) ? printf("%d %d %d %d \n", search_size[index], index+1, decimal, sign) : 0;
}
void hyouji_on_site(int matrix_size, int index, int counta){
// (counta != 0) ? printf("%d %d \n", index+1, counta) : 0;
}
int combination(int n, int r) {
int result = 1;
for (int i = 1; i <= r; i++) {
result *= (n - i + 1);
result /= i;
}
return result;
}
double dot_product(int length, double vector_1[length], double vector_2[length]){ //内積
double dot = 0.0;
for(int i = 0; i < length; i++){
dot += vector_1[i] * vector_2[i];
}
return dot;
}
double complex ket_complex(int length, double complex vector_1[length], double complex vector_2[length]){ //内積 複素数
double complex dot_complex = 0.0;
for(int i = 0; i < length; i++){
dot_complex += conj(vector_1[i]) * vector_2[i];
}
// printf("%lf %lf\n", creal(dot), cimag(dot));
return dot_complex ;
}
double norm(int length, double vector[length]){ //ノルム
double norm = 0.0;
norm = dot_product(length, vector, vector);
norm = pow(norm, 0.5);
return norm;
}
double norm_complex(int length, double complex vector[length]){ //ノルム
double norm = 0.0;
norm = ket_complex(length, vector, vector);
norm = pow(norm, 0.5);
return norm;
}
void normalization(int length, double vector[length]){
double n = 0.0;
n = norm(length, vector);
for(int i = 0; i < length; i++){
vector[i] = vector[i] / n;
}
}
void normalization_complex(int length, double complex vector[length]){
double n = 0.0;
n = norm_complex(length, vector);
for(int i = 0; i < length; i++){
vector[i] = vector[i] / n;
}
}
//random vector creation
void r_v_c(int size_length, double vector_q_0[size_length]){
// 乱数の初期化
srand(time(NULL));
// ランダムな成分からなるベクトルを生成
for (int i = 0; i < size_length; i++) {
vector_q_0[i] = ((double)rand() / RAND_MAX) * 2.0 - 1.0;
}
double norm = 0.0;
for (int i = 0; i < size_length; i++) {
norm += vector_q_0[i] * vector_q_0[i];
}
norm = sqrt(norm);
for (int i = 0; i < size_length; i++) {
vector_q_0[i] /= norm;
}
}
void r_v_c_c(int size_length, double complex vector_q_0[size_length]){
// 乱数の初期化
srand(time(NULL));
// ランダムな成分からなるベクトルを生成
for (int i = 0; i < size_length; i++) {
vector_q_0[i] = ((double)rand() / RAND_MAX) * 2.0 - 1.0;
vector_q_0[i] += I * (((double)rand() / RAND_MAX) * 9.90 / 7 - 1.0);
}
normalization_complex(size_length, vector_q_0);
}
int countCrossing_new(int binary_1[N], int binary_2[N], int index_left) { //operator
int andresult[N] = {0};
int count = 0, sign = 1;
for(int i = 0; i < index_left; i++){
andresult[i] = binary_1[i] & binary_2[i];
}
for(int i = 0; i < N; i++){
count += andresult[i];
}
sign = pow(-1, count);
return sign;
}
// 二進数から10進数への変換関数
int binaryToDecimal(int *binary, int length) {
int decimal = 0;
for (int i = 0; i < length; i++) {
decimal = (decimal << 1) | binary[i];
}
return decimal;
}
void decimalToBinary(int decimal, int *binary, int length) {
for (int i = length - 1; i >= 0; i--) {
binary[i] = decimal & 1;
decimal >>= 1;
}
}
// ケットベクトルを2進法でket配列に保存する関数
void fillKet_up(int ones, int *binary, int binary_length, int binary_index, int *ket_libra_junbiki_up, int num_ket_length_up){
if (ket_index_up >= num_ket_length_up) {
return;
}
if (binary_index >= binary_length) {
return;
}
if (ones > 0) {
// 1を配置してみる
binary[binary_index] = 1;
fillKet_up(ones - 1, binary, binary_length, binary_index + 1, ket_libra_junbiki_up, num_ket_length_up);
if (ones - 1 == 0) {
ket_libra_junbiki_up[num_ket_length_up - ket_index_up - 1] = binaryToDecimal(binary, binary_length) ;
ket_index_up++;
}
}
// 0を配置してみる (onesが0でも残りの桁を0で埋める必要がある)
binary[binary_index] = 0;
fillKet_up(ones, binary, binary_length, binary_index + 1, ket_libra_junbiki_up, num_ket_length_up);
}
void fillKet_down(int ones, int *binary, int binary_length, int binary_index, int *ket_libra_junbiki_down, int num_ket_length_down) {
if (ket_index_down >= num_ket_length_down) {
return;
}
if (binary_index >= binary_length) {
return;
}
if (ones > 0) {
// 1を配置してみる
binary[binary_index] = 1;
fillKet_down(ones - 1, binary, binary_length, binary_index + 1, ket_libra_junbiki_down, num_ket_length_down);
if (ones - 1 == 0) {
ket_libra_junbiki_down[num_ket_length_down - ket_index_down - 1] = binaryToDecimal(binary, binary_length) ;
ket_index_down++;
}
}
// 0を配置してみる (onesが0でも残りの桁を0で埋める必要がある)
binary[binary_index] = 0;
fillKet_down(ones, binary, binary_length, binary_index + 1, ket_libra_junbiki_down, num_ket_length_down);
// printf("ket_libra_junbiki_down:\n");
// for (int i = 0; i < num_ket_length_down; i++) {
// printf("%d\n", ket_libra_junbiki_down[i]);
// }
}
int ket_libra_gyakubiki_func(int size_length, int decimal, int search_size[size_length]) { // 辞書 [decimal]⇒index 逆引き
int index = 0;
while (index < size_length && search_size[index] != decimal) {
index++;
}
return (index < size_length) ? (index + 1) : 0;
}
void spin_separate(int binary[L], int binary_up[N], int binary_down[N]){
for(int i = 0; i < N; i++){
binary_up[i] = binary[i];
binary_down[i] = binary[i + N];
}
}
void spin_connect(int binary[L], int binary_up[N], int binary_down[N]){
for(int i = 0; i < N; i++){
binary[i] = binary_up[i];
binary[i + N] = binary_down[i];
}
}
void element_search(int *ket_libra_junbiki_up, int *ket_libra_junbiki_down, int *search_size, int num_ket_length_up, int num_ket_length_down){
int binary[L] = {0.0}, binary_up[N] = {0.0}, binary_down[N] = {0.0}; shokika_int(L, binary); shokika_int(N, binary_up); shokika_int(N, binary_down);
for(int i = 0; i < num_ket_length_up; i++){
for(int j = 0; j < num_ket_length_down; j++){
decimalToBinary(ket_libra_junbiki_up[i], binary_up, N);
decimalToBinary(ket_libra_junbiki_down[j], binary_down, N);
spin_connect(binary, binary_up, binary_down);
search_size[i * num_ket_length_down + j] = binaryToDecimal(binary, L);
// printf("%d\n", binaryToDecimal(binary, L));
}
// printf("%d\n",ket_libra_junbiki_down[i] );
}
}
void make_search(int num_ket_length_up, int num_ket_length_down, int *search_size, int one_up, int one_down, int SPIN_up, int SPIN_down){
int binary_spin[L] = {0}; shokika_int(L, binary_spin);
// ケットベクトルを2進法で保存する配列
int ket_libra_junbiki_up[num_ket_length_up]; int ket_libra_junbiki_down[num_ket_length_down]; // 辞書 [index]⇒decimal 順引き
shokika_int(num_ket_length_up, ket_libra_junbiki_up); shokika_int(num_ket_length_down, ket_libra_junbiki_down);
ket_index_up = 0, ket_index_down = 0;
fillKet_up(SPIN_up + one_up, binary_spin, N, 0, ket_libra_junbiki_up, num_ket_length_up);
fillKet_down(SPIN_down + one_down, binary_spin, N, 0, ket_libra_junbiki_down, num_ket_length_down);
element_search(ket_libra_junbiki_up, ket_libra_junbiki_down, search_size, num_ket_length_up, num_ket_length_down);
}
void search(int search_2D[size], int search_creation_up[size_creation_up], int search_creation_down[size_creation_down], int search_annihilation_up[size_annihilation_up], int search_annihilation_down[size_annihilation_down], int SPIN_up, int SPIN_down){
shokika_int(size, search_2D);
shokika_int(size_creation_up, search_creation_up); shokika_int(size_creation_down, search_creation_down);
shokika_int(size_annihilation_up, search_annihilation_up); shokika_int(size_annihilation_down, search_annihilation_down);
make_search(NUM_KET_creation_up, NUM_KET_down, search_creation_up, 1, 0, SPIN_up, SPIN_down); make_search(NUM_KET_up, NUM_KET_creation_down, search_creation_down, 0, 1, SPIN_up, SPIN_down);
make_search(NUM_KET_annihilation_up, NUM_KET_down, search_annihilation_up, -1, 0, SPIN_up, SPIN_down); make_search(NUM_KET_up, NUM_KET_annihilation_down, search_annihilation_down, 0, -1, SPIN_up, SPIN_down);
make_search(NUM_KET_up, NUM_KET_down, search_2D, 0, 0, SPIN_up, SPIN_down);
}
int flipBits01_ribbon(int *binary, int i, int is01) {
is01 = (binary[i] == 1) && (binary[i + 1] == 0);
(is01 == 1) ? (binary[i] = 0, binary[i + 1] = 1) : 0;
return is01;
}
int flipBits10_ribbon(int *binary, int i, int is10) {
is10 = (binary[i] == 0) && (binary[i + 1] == 1);
(is10 == 1) ? (binary[i] = 1, binary[i + 1] = 0) : 0;
// printf("%d %d %d\n", is10, binary[i], binary[i + 1]);
return is10;
}
void make_hamiltonian_complex(int matrix_size, double complex out_put_vector[matrix_size], int index, int search_size[matrix_size],
int decimal, int is, double complex vector_q[matrix_size]) {
(is != 0) ? out_put_vector[ket_libra_gyakubiki_func(matrix_size, decimal, search_size) - 1] += - t_opt * (is) * vector_q[index]: 0;
}
void make_hamiltonian_diag_complex(int matrix_size, double complex out_put_vector[matrix_size], int counta, int index, double U_t, double complex vector_q[matrix_size]) {
out_put_vector[index] += U_t * counta * vector_q[index];
}
void hopping_up_complex(int matrix_size, double complex out_put_vector[matrix_size], int binary_up[N], int binary_down[N], int search_size[matrix_size],
int index, double complex vector_q_0[matrix_size]){
int binary_up_hop[N], binary_convert[L], is01 = 0, is10 = 0, decimal = 0; shokika_int(N, binary_up_hop); shokika_int(L, binary_convert);
for(int i = 0; i < N; i++){
copy_int(N, binary_up, binary_up_hop); is01 = flipBits01_ribbon(binary_up_hop, i, is01);
(is01 == 1) ? (spin_connect(binary_convert, binary_up_hop, binary_down),
decimal = binaryToDecimal(binary_convert, L),
make_hamiltonian_complex(matrix_size, out_put_vector, index, search_size, decimal, is01, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is01);
}
for(int i = 0; i < N; i++){
copy_int(N, binary_up, binary_up_hop); is10 = flipBits10_ribbon(binary_up_hop, i, is10);
(is10 == 1) ? (spin_connect(binary_convert, binary_up_hop, binary_down),
decimal = binaryToDecimal(binary_convert, L),
make_hamiltonian_complex(matrix_size, out_put_vector, index, search_size, decimal, is10, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is10);
}
}
void hopping_down_complex(int matrix_size, double complex out_put_vector[matrix_size], int binary_up[N], int binary_down[N], int search_size[matrix_size],
int index, double complex vector_q_0[matrix_size]){
int binary_down_hop[N], binary_convert[L], is01 = 0, is10 = 0, decimal = 0; shokika_int(N, binary_down_hop); shokika_int(L, binary_convert);
for(int i = 0; i < N; i++){
copy_int(N, binary_down, binary_down_hop); is01 = flipBits01_ribbon(binary_down_hop, i, is01);
(is01 == 1) ? (spin_connect(binary_convert, binary_up, binary_down_hop),
decimal = binaryToDecimal(binary_convert, L),
make_hamiltonian_complex(matrix_size, out_put_vector, index, search_size, decimal, is01, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is01);
}
for(int i = 0; i < N; i++){
copy_int(N, binary_down, binary_down_hop); is10 = flipBits10_ribbon(binary_down_hop, i, is10);
(is10 == 1) ? (spin_connect(binary_convert, binary_up, binary_down_hop),
decimal = binaryToDecimal(binary_convert, L),
make_hamiltonian_complex(matrix_size, out_put_vector, index, search_size, decimal, is10, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is10);
}
}
void matrix_element_diag_on_site_complex(int matrix_size, double complex out_put_vector[matrix_size], int binary_up[N], int binary_down[N], double U_t,
int index, double complex vector_q_0[matrix_size]){
int counta = 0, binary_diag[N]; shokika_int(N, binary_diag);
for(int i = 0; i < N; i++){
binary_diag[i] = binary_up[i] & binary_down[i];
counta += binary_diag[i];
}
(counta != 0) ? make_hamiltonian_diag_complex(matrix_size, out_put_vector, counta, index, U_t, vector_q_0) : 0;
hyouji_on_site(matrix_size, index, counta);
}
void hamiltonin_element_complex(int matrix_size, double complex out_put_vector[matrix_size], int search_size[matrix_size], int index,
double U_t, double complex vector_q_0[matrix_size]){
int binary_spin[L], binary_up[N], binary_down[N]; shokika_int(L, binary_spin); shokika_int(N, binary_up); shokika_int(N, binary_down);
decimalToBinary(search_size[index], binary_spin, L); spin_separate(binary_spin, binary_up, binary_down);
hopping_up_complex(matrix_size, out_put_vector, binary_up, binary_down, search_size, index, vector_q_0);
hopping_down_complex(matrix_size, out_put_vector, binary_up, binary_down, search_size, index, vector_q_0);
matrix_element_diag_on_site_complex(matrix_size, out_put_vector, binary_up, binary_down, U_t, index, vector_q_0);
}
void matrix_complex(int matrix_size, double complex out_put_vector[matrix_size], int search_size[matrix_size], double U_t, double complex vector_q_0[matrix_size]){
for(int index = 0; index < matrix_size; index++){
hamiltonin_element_complex(matrix_size, out_put_vector, search_size, index, U_t, vector_q_0);
}
}
int main(){
double U_t = 0.0, mu_chemical = U_t / 2.0;
int SPIN_up = 2, SPIN_down = 2, filling = SPIN_up + SPIN_down;
// chemical_filling(mu_chemical, U_t, &SPIN_up, &SPIN_down);
NUM_KET_up = combination(N, SPIN_up); NUM_KET_down = combination(N, SPIN_down);
NUM_KET_creation_up = combination(N, SPIN_up + 1); NUM_KET_creation_down = combination(N, SPIN_down + 1);
NUM_KET_annihilation_up = combination(N, SPIN_up - 1); NUM_KET_annihilation_down = combination(N, SPIN_down - 1);
size = NUM_KET_up * NUM_KET_down;
size_creation_up = NUM_KET_creation_up * NUM_KET_down; size_creation_down = NUM_KET_up * NUM_KET_creation_down;
size_annihilation_up = NUM_KET_annihilation_up * NUM_KET_down; size_annihilation_down = NUM_KET_up * NUM_KET_annihilation_down;
//辞書 normal creation annihilatioin
int search_2D[size];
int search_creation_up[size_creation_up], search_creation_down[size_creation_down];
int search_annihilation_up[size_annihilation_up], search_annihilation_down[size_annihilation_down];
search(search_2D, search_creation_up, search_creation_down, search_annihilation_up, search_annihilation_down, SPIN_up, SPIN_down);
double complex vector_q_cm_1[size]; shokika_complex(size, vector_q_cm_1);
r_v_c_c(size, vector_q_cm_1);
double complex Aq[size]; shokika_complex(size, Aq);
// matrix_complex(size, Aq, search_2D, U_t, vector_q_cm_1);
// // printf("%d\n", size);
hamiltonin_element_complex(size, Aq, search_2D, 0, U_t, vector_q_cm_1);
return 0;
}
自分で試したこと
matrix_complex()とhamiltonin_element_complex()が異常を引き起こすことが分かったので、
1.matrix_complex()の計算について、関数内のindexの上限を変えて回るのか確認したところ、1回でもセグメしました。
2.hamiltonin_element_complex()で呼び出している関数をコメントアウトして回したところ、hopping_up_complex(), hopping_down_complex()のどちらかをコメントアウトすればセグメしなくなる、という感じでした。
正直お手上げです。何かデバッグオプションや岡目八目的な助言をいただけないでしょうか。
環境について、細かいことはわかりませんが、
メモリ:8Gb
os:windows11でWSL2
コマンド:icx -qmkl
です。足りなければ申し訳ございません。