セグメの原因がわからない。
解決したいこと
セグメの原因がわからない。
下のコードは物理のハミルトニアンを構築し、固有値問題をLANCZOS法を用いて計算するプログラムです。
現状、Nの値を8以下の偶数にし、spinをちょうど半分の数字にすると動くのだが、N＝10でspin=5にするとセグメしてしまいます。またspinが4以下ですと問題なく動きます。
発生している問題・エラー
Segmentation fault
以下はgdbした場合の結果です。
GNU gdb (Ubuntu 9.2-0ubuntu1~20.04.1) 9.2
Copyright (C) 2020 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.
Type "show copying" and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
Type "show configuration" for configuration details.
For bug reporting instructions, please see:
--Type <RET> for more, q to quit, c to continue without paging--
<http://www.gnu.org/software/gdb/bugs/>.
Find the GDB manual and other documentation resources online at:
<http://www.gnu.org/software/gdb/documentation/>.
For help, type "help".
Type "apropos word" to search for commands related to "word"...
Reading symbols from ./a.out...
(No debugging symbols found in ./a.out)
(gdb) run
Starting program: /home/sagyouspace/tiba_gakui/change/lanczos/L=10/a.out
warning: the debug information found in "/opt/intel/oneapi/compiler/2024.0/lib/libiomp5.dbg" does not match "/opt/intel/oneapi/compiler/2024.0/lib/libiomp5.so" (CRC mismatch).
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".
Program received signal SIGSEGV, Segmentation fault.
0x0000000000406712 in lanczos_first ()
(gdb) bt
#0 0x0000000000406712 in lanczos_first ()
#1 0x0000000000406ef8 in lanczos_step ()
#2 0x0000000000407e91 in lanczos_vector ()
#3 0x0000000000408b23 in make_eigenvector_green ()
#4 0x0000000000409043 in main ()
ここにあるようにlanczos_first ()に問題があるとのことですが、コメントアウトしてもあまり効果はなかったです。むしろ、make_eigenvector()の部分であったり、lanczos_vector ()で問題が起きているようですが、解決方法や根本的な問題を見つけることができていません。
該当するソースコード
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <stdlib.h>
#define grand_potential_convergence 10e-7
#define Matsubara_frequency 1024 * 4
#define Lanczos_stop 8
#define N 10
#define L 2*N
#define Cn 100
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, size_whole_system = 0, size_gyaku = 0;
int M = 0, limit = 0; double t_opt = 1.0, k_bT = 0.001;
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 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_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 norm(int length, double vector[length]){ //ノルム
double norm = 0.0;
norm = dot_product(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;
}
}
//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;
}
}
// 二進数から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;
}
}
//対格化
void dsyev_(char* JOBZ, char* UPLO, int* n, double* A, int* LDA, double* W, double* WORK, int* LWORK, int* INFO);
void diagonalizaton(int matrix_size, double hamiltonian[matrix_size][matrix_size], double Energy[matrix_size], double eigenvector[matrix_size][matrix_size]){
char JOBZ = 'V';
char UPLO = 'U';
int n = matrix_size;
double A[matrix_size * matrix_size];
for(int i = 0; i < matrix_size; i++){
for(int j = 0; j < matrix_size; j++){
A[i * matrix_size + j] = hamiltonian[i][j];
}
}
int LDA = matrix_size;
double W[matrix_size];
double WORK[6 * matrix_size];
int LWORK = 6 * matrix_size;
int INFO;
dsyev_(&JOBZ, &UPLO, &n, A, &LDA, W, WORK, &LWORK, &INFO);
for(int i = 0; i < matrix_size; i++){
// printf("%lf\n",W[i]);
Energy[i] = W[i];
}
for(int i = 0; i < matrix_size; i++){
for(int j = 0; j < matrix_size; j++){
eigenvector[i][j] = A[i * matrix_size + j]; //基底状態の固有ベクトル
// printf("%d %lf\n", j, lanczos_eigenvector[j]);
}
}
}
void diagonalizaton_eigenvalue(int matrix_size, double hamiltonian[matrix_size][matrix_size], double Energy[matrix_size]){
char JOBZ = 'V';
char UPLO = 'U';
int n = matrix_size;
double A[matrix_size * matrix_size];
for(int i = 0; i < matrix_size; i++){
for(int j = 0; j < matrix_size; j++){
A[i * matrix_size + j] = hamiltonian[i][j];
}
}
int LDA = matrix_size;
double W[matrix_size];
double WORK[6 * matrix_size];
int LWORK = 6 * matrix_size;
int INFO;
dsyev_(&JOBZ, &UPLO, &n, A, &LDA, W, WORK, &LWORK, &INFO);
for(int i = 0; i < matrix_size; i++){
// printf("%lf\n",W[i]);
Energy[i] = W[i];
}
}
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];
}
}
int get(int decimal, int SPIN_up, int SPIN_down){
int i_up = 0, i_down = 0, binary[L] = {0}; decimalToBinary(decimal, binary, L);
for(int i = 0; i < N; i++){
i_up += binary[i]; i_down += binary[i + N];
}
int check_spin = - 1; (i_up == SPIN_up && i_down == SPIN_down) ? check_spin = 0: (check_spin = - 1);
return check_spin;
}
void make_search(int size_length, int search_junbiki[size_length], int search_gyakubiki[size_gyaku], int SPIN_up, int SPIN_down){
int i_cnt = 0;
for(int i = 0; i < pow(2, L) && i_cnt < size_length; i++){
int check_spin = get(i, SPIN_up, SPIN_down);
(check_spin == 0) ? (search_junbiki[i_cnt] = i, search_gyakubiki[search_junbiki[i_cnt]] = i_cnt, i_cnt++): 0;
}
}
void search_chemical(int search_2D[size], int search_gyakubiki[size_gyaku], int SPIN_up, int SPIN_down){
shokika_int(size, search_2D); shokika_int(size_gyaku, search_gyakubiki); make_search(size, search_2D, search_gyakubiki, SPIN_up, SPIN_down);
}
void search(int search_2D[size], int search_gyakubiki[size_gyaku],
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_gyaku, search_gyakubiki);
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(size, search_2D, search_gyakubiki, SPIN_up, SPIN_down);
make_search(size_creation_up, search_creation_up, search_gyakubiki, SPIN_up + 1, SPIN_down);
make_search(size_annihilation_up, search_annihilation_up, search_gyakubiki, SPIN_up - 1, SPIN_down);
make_search(size_creation_down, search_creation_down, search_gyakubiki, SPIN_up, SPIN_down + 1);
make_search(size_annihilation_down, search_annihilation_down, search_gyakubiki, SPIN_up, SPIN_down - 1);
}
int flipBits01_ribbon(int binary[N], 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[N], 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(int matrix_size, double out_put_vector[matrix_size], int index, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
int decimal, int is, double vector_q[matrix_size]) {
(is != 0) ? out_put_vector[search_gyakubiki[decimal]] += - t_opt * (is) * vector_q[index]: 0;
}
void make_hamiltonian_diag(int matrix_size, double out_put_vector[matrix_size], int counta, int index, double U_t, double vector_q[matrix_size]) {
(counta != 0) ? out_put_vector[index] += U_t * counta * vector_q[index]: 0;
}
void hopping_up(int matrix_size, double out_put_vector[matrix_size], int binary_up[N], int binary_down[N],
int search_size[matrix_size], int search_gyakubiki[size_gyaku], int index, double 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 - 1; 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(matrix_size, out_put_vector, index, search_size, search_gyakubiki, decimal, is01, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is01);
}
for(int i = 0; i < N - 1; 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(matrix_size, out_put_vector, index, search_size, search_gyakubiki, decimal, is10, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is10);
}
}
void hopping_down(int matrix_size, double out_put_vector[matrix_size], int binary_up[N], int binary_down[N],
int search_size[matrix_size], int search_gyakubiki[size_gyaku], int index, double 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 - 1; 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(matrix_size, out_put_vector, index, search_size, search_gyakubiki, decimal, is01, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is01);
}
for(int i = 0; i < N - 1; 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(matrix_size, out_put_vector, index, search_size, search_gyakubiki, decimal, is10, vector_q_0)) : 0;
hyouji(matrix_size, search_size, decimal, 1, index, is10);
}
}
void matrix_element_diag_on_site(int matrix_size, double out_put_vector[matrix_size], int binary_up[N], int binary_down[N], double U_t, int index, double 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(matrix_size, out_put_vector, counta, index, U_t, vector_q_0) : 0;
hyouji_on_site(matrix_size, index, counta);
}
void matrix_element_diag_chemical(int matrix_size, double out_put_vector[matrix_size], int binary_up[N], int binary_down[N],
double mu_chemical, double variational_chemical, int index, double 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(matrix_size, out_put_vector, counta, index, mu_chemical + variational_chemical, vector_q_0) : 0;
}
void hamiltonin_element(int matrix_size, double out_put_vector[matrix_size], int search_size[matrix_size], int search_gyakubiki[size_gyaku],
int index, double U_t, double 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(matrix_size, out_put_vector, binary_up, binary_down, search_size, search_gyakubiki, index, vector_q_0);
hopping_down(matrix_size, out_put_vector, binary_up, binary_down, search_size, search_gyakubiki, index, vector_q_0);
matrix_element_diag_on_site(matrix_size, out_put_vector, binary_up, binary_down, U_t, index, vector_q_0);
}
void matrix(int matrix_size, double out_put_vector[matrix_size], int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double U_t, double vector_q_0[matrix_size]){
for(int index = 0; index < matrix_size; index++){
hamiltonin_element(matrix_size, out_put_vector, search_size, search_gyakubiki, index, U_t, vector_q_0);
}
}
double variance(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku], double eigenvalue[100], double eigenvector[matrix_size], double U_t){
double Hq[matrix_size], Eq[matrix_size]; shokika_double(matrix_size, Hq); shokika_double(matrix_size, Eq);
double dot_Hq = 0.0, dot_Eq = 0.0;
for(int i = 0; i < matrix_size; i++){
Eq[i] = eigenvalue[0] * eigenvector[i];
}
matrix(matrix_size, Hq, search_size, search_gyakubiki, U_t, eigenvector);
dot_Hq = dot_product(matrix_size, eigenvector, Hq); dot_Eq = dot_product(matrix_size, eigenvector, Eq);
// printf("%.10lf %.10lf\n", dot_Eq, dot_Hq);
double variance = fabs(dot_Hq - dot_Eq);
return variance;
}
double new_variance(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku], double eigenvector[matrix_size], double U_t){
double Hq[matrix_size]; shokika_double(matrix_size, Hq);
double qHHq = 0.0, qHq = 0.0;
matrix(matrix_size, Hq, search_size, search_gyakubiki, U_t, eigenvector);
qHHq = dot_product(matrix_size, Hq, Hq); qHq = dot_product(matrix_size, eigenvector, Hq);
double variance = fabs(qHHq - qHq * qHq);
return variance;
}
void make_eigenvector(int matrix_size, double eigenvector_re[matrix_size], double vector_q_0[matrix_size],
double lanczos_eigenvector[100], int index){
for(int i = 0; i < matrix_size; i++){
eigenvector_re[i] += lanczos_eigenvector[index] * vector_q_0[i];
// printf("%lf %lf \n", eigenvector_re[i], lanczos_eigenvector[index]);
}
}
void diago(double lanczos_matrix[M][M], double Energy[100], double lanczos_eigenvector[100]){
shokika_double(100, Energy), shokika_double(100, lanczos_eigenvector);
char JOBZ = 'V';
char UPLO = 'U';
int n = M;
double A[M * M];
for(int i = 0; i < M; i++){
for(int j = 0; j < M; j++){
A[i * M + j] = lanczos_matrix[i][j];
}
}
int LDA = M;
double W[M];
double WORK[6 * M];
int LWORK = 6 * M;
int INFO;
dsyev_(&JOBZ, &UPLO, &n, A, &LDA, W, WORK, &LWORK, &INFO);
for(int i = 0; i < M; i++){
Energy[i] = W[i];
}
for(int j = 0; j < M; j++){
lanczos_eigenvector[j] = A[0 * M + j]; //基底状態の固有ベクトル
}
}
void make_lanczos_matrix(double lanczos_matrix[M][M], double alpha[101], double beta[100]){
lanczos_matrix[0][0] = alpha[0]; lanczos_matrix[0][1] = beta[0]; lanczos_matrix[1][0] = beta[0];
for(int i = 1; i < M - 1; i++){
lanczos_matrix[i][i] = alpha[i], lanczos_matrix[i][i + 1] = beta[i], lanczos_matrix[i + 1][i] = beta[i];
}
lanczos_matrix[M - 1][M - 1] = alpha[M - 1];
}
void kakunin_diago(double alpha[101], double beta[100], double eigenvalue[100], double lanczos_eigenvector[100], double pre_eigenvalue[2]){
double lanczos_matrix[M][M]; shokika_matrix(M, lanczos_matrix);
make_lanczos_matrix(lanczos_matrix, alpha, beta);
diago(lanczos_matrix, eigenvalue, lanczos_eigenvector);
}
double break_point_kitei(double eigenvalue[100], double pre_eigenvalue[2]){
double differ = eigenvalue[0] - pre_eigenvalue[0];
pre_eigenvalue[0] = eigenvalue[0];
return fabs(differ);
}
double break_point_daiitireiki(double eigenvalue[100], double pre_eigenvalue[2]){
double differ = eigenvalue[1] - pre_eigenvalue[1];
pre_eigenvalue[1] = eigenvalue[1];
return fabs(differ);
}
//lanczos eigenvalue
void lanczos_first(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double vector_q_1[matrix_size], double alpha[101], double beta[100], double U_t){
double a = 0, b = 0, Aq[matrix_size], Bq[matrix_size]; shokika_double(matrix_size, Aq), shokika_double(matrix_size, Bq);
matrix(matrix_size, Aq, search_size, search_gyakubiki, U_t, vector_q_0);
a = dot_product(matrix_size, vector_q_0, Aq); alpha[0] = a;
for(int i = 0; i < matrix_size; i++){
Bq[i] = Aq[i] - a * vector_q_0[i];
}
b = norm(matrix_size, Bq); beta[0] = b;
for(int i = 0; i < matrix_size; i++){
vector_q_1[i] = Bq[i] / b;
}
// printf("%lf %lf\n", a, b);
}
void lanczos_middle(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double vector_q_1[matrix_size], double alpha[101], double beta[100], int index, double U_t){
double a = 0.0, b = 0.0, Aq[matrix_size], Bq[matrix_size]; shokika_double(matrix_size, Aq), shokika_double(matrix_size, Bq);
double vector_q_2[matrix_size]; shokika_double(matrix_size, vector_q_2);
matrix(matrix_size, Aq, search_size, search_gyakubiki, U_t, vector_q_1);
a = dot_product(matrix_size, vector_q_1, Aq); alpha[index] = a;
for(int i = 0; i < matrix_size; i++){
Bq[i] = Aq[i] - a * vector_q_1[i] - beta[index - 1] * vector_q_0[i];
}
b = norm(matrix_size, Bq); beta[index] = b;
for(int i = 0; i < matrix_size; i++){
vector_q_2[i] = Bq[i] / b;
vector_q_0[i] = vector_q_1[i];
vector_q_1[i] = vector_q_2[i];
}
// printf("%lf %lf\n", a, b);
}
void lanczos_step(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double alpha[101], double beta[100], double eigenvalue[100], double lanczos_eigenvector[100], double U_t){
double vector_q_1[matrix_size], vector_q_2[matrix_size]; shokika_double(matrix_size, vector_q_1); shokika_double(matrix_size, vector_q_2);
double pre_eigenvalue[2] = {0.0}, differ_kitei_value = 10.0, differ_daiitireiki_value = 10.0;
lanczos_first(matrix_size, search_size, search_gyakubiki, vector_q_0, vector_q_1, alpha, beta, U_t);
int lanczos_index = 1, should_break = 0; M = 0, limit = 0;
while(lanczos_index < 100 && should_break == 0){
M = lanczos_index;
lanczos_middle(matrix_size, search_size, search_gyakubiki, vector_q_0, vector_q_1, alpha, beta, lanczos_index, U_t);
(lanczos_index == Lanczos_stop) ? kakunin_diago(alpha, beta, eigenvalue, lanczos_eigenvector, pre_eigenvalue) : 0;
(lanczos_index > Lanczos_stop && lanczos_index % 2 == 0) ? kakunin_diago(alpha, beta, eigenvalue, lanczos_eigenvector, pre_eigenvalue) : 0;
(lanczos_index > Lanczos_stop && lanczos_index % 2 == 0) ? (differ_kitei_value = break_point_kitei(eigenvalue, pre_eigenvalue),
differ_daiitireiki_value = break_point_daiitireiki(eigenvalue, pre_eigenvalue)): 0;
should_break = (differ_kitei_value < 1e-6 && differ_daiitireiki_value < 1e-6);
// printf("%d %lf %lf\n", should_break, differ_kitei_value, differ_daiitireiki_value);
lanczos_index++;
}
limit = M;
}
//lanczos eigenvector
void lanczos_first_re(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double vector_q_1[matrix_size], double alpha[101], double beta[100],
double lanczos_eigenvector[100], double eigenvector_re[matrix_size], double U_t){
double a = 0, b = 0, Aq[matrix_size], Bq[matrix_size]; shokika_double(matrix_size, Aq), shokika_double(matrix_size, Bq);
matrix(matrix_size, Aq, search_size, search_gyakubiki, U_t, vector_q_0);
a = dot_product(matrix_size, vector_q_0, Aq); alpha[0] = a;
for(int i = 0; i < matrix_size; i++){
Bq[i] = Aq[i] - a * vector_q_0[i];
}
b = norm(matrix_size, Bq); beta[0] = b;
for(int i = 0; i < matrix_size; i++){
vector_q_1[i] = Bq[i] / b;
}
make_eigenvector(matrix_size, eigenvector_re, vector_q_0, lanczos_eigenvector, 0);
}
void lanczos_middle_re(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double vector_q_1[matrix_size],
double alpha[101], double beta[100], double lanczos_eigenvector[100], double eigenvector_re[matrix_size], int index, double U_t){
double a = 0, b = 0, Aq[matrix_size], Bq[matrix_size]; shokika_double(matrix_size, Aq), shokika_double(matrix_size, Bq);
double vector_q_2[matrix_size]; shokika_double(matrix_size, vector_q_2);
matrix(matrix_size, Aq, search_size, search_gyakubiki, U_t, vector_q_1);
a = dot_product(matrix_size, vector_q_1, Aq); alpha[index] = a;
for(int i = 0; i < matrix_size; i++){
Bq[i] = Aq[i] - a * vector_q_1[i] - beta[index - 1] * vector_q_0[i];
}
b = norm(matrix_size, Bq); beta[index] = b;
for(int i = 0; i < matrix_size; i++){
vector_q_2[i] = Bq[i] / b;
vector_q_0[i] = vector_q_1[i];
vector_q_1[i] = vector_q_2[i];
}
make_eigenvector(matrix_size, eigenvector_re, vector_q_0, lanczos_eigenvector, index);
}
void lanczos_last_re(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_1[matrix_size], double alpha[101], int index, double U_t){
double a = 0, Aq[matrix_size]; shokika_double(matrix_size, Aq);
matrix(matrix_size, Aq, search_size, search_gyakubiki, U_t, vector_q_1);
a = dot_product(matrix_size, vector_q_1, Aq); alpha[index] = a;
}
void lanczos_step_re(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku], double vector_q_0[matrix_size],
double alpha[101], double beta[100], double lanczos_eigenvector[100], double eigenvector_re[matrix_size], double U_t){
double vector_q_1[matrix_size], vector_q_2[matrix_size]; shokika_double(matrix_size, vector_q_1); shokika_double(matrix_size, vector_q_2);
lanczos_first_re(matrix_size, search_size, search_gyakubiki, vector_q_0, vector_q_1, alpha, beta, lanczos_eigenvector, eigenvector_re, U_t);
int lanczos_index_re = 1;
while(lanczos_index_re <= limit){
lanczos_middle_re(matrix_size, search_size, search_gyakubiki, vector_q_0, vector_q_1, alpha, beta, lanczos_eigenvector, eigenvector_re, lanczos_index_re, U_t);
lanczos_index_re++;
}
lanczos_last_re(matrix_size, search_size, search_gyakubiki, vector_q_1, alpha, lanczos_index_re, U_t);
M += 1;
}
void lanczos_vector(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku],
double vector_q_0[matrix_size], double alpha[101], double beta[100], double eigenvalue[100], double eigenvector[matrix_size], double U_t){
int z = 0; double v = 10.0, vector_Q[matrix_size], lanczos_eigenvector[100], eigenvector_re[matrix_size]; normalization(matrix_size, vector_q_0);
shokika_double(matrix_size, vector_Q), shokika_double(101, alpha), shokika_double(100, beta);
shokika_double(100, lanczos_eigenvector), shokika_double(matrix_size, eigenvector_re);
while(z < 2 || v > 1e-8){
shokika_double(matrix_size, vector_Q), shokika_double(101, alpha), shokika_double(100, beta);
shokika_double(100, lanczos_eigenvector), shokika_double(matrix_size, eigenvector_re); copy_double(matrix_size, vector_q_0, vector_Q);
lanczos_step(matrix_size, search_size, search_gyakubiki, vector_q_0, alpha, beta, eigenvalue, lanczos_eigenvector, U_t);
shokika_double(matrix_size, vector_q_0), shokika_double(101, alpha), shokika_double(100, beta);
// copy_double(matrix_size, vector_Q, vector_q_0); shokika_double(matrix_size, vector_Q);
lanczos_step_re(matrix_size, search_size, search_gyakubiki, vector_Q, alpha, beta, lanczos_eigenvector, eigenvector_re, U_t);
// copy_double(matrix_size, eigenvector_re, vector_q_0); shokika_double(matrix_size, eigenvector_re); normalization(matrix_size, vector_q_0);
v = variance(matrix_size, search_size, search_gyakubiki, eigenvalue, vector_q_0, U_t);
// v = new_variance(matrix_size, search_size, search_gyakubiki, vector_q_0, U_t);
printf("%d %d\n", M, matrix_size);
z++;
}
copy_double(matrix_size, vector_q_0, eigenvector); shokika_double(matrix_size, vector_q_0);
}
void lanczos(int matrix_size, int search_size[matrix_size], int search_gyakubiki[size_gyaku], double vector_q_0[matrix_size], double alpha[101], double beta[100],
double eigenvalue[100], double eigenvector[matrix_size], double U_t){
int z = 0; double v = 10.0, vector_Q[matrix_size], lanczos_eigenvector[100], eigenvector_re[matrix_size]; normalization(matrix_size, vector_q_0);
while(z < 10 && v > 1e-8){
shokika_double(matrix_size, vector_Q), shokika_double(101, alpha), shokika_double(100, beta), shokika_double(100, lanczos_eigenvector), shokika_double(matrix_size, eigenvector_re); copy_double(matrix_size, vector_q_0, vector_Q);
lanczos_step(matrix_size, search_size, search_gyakubiki, vector_q_0, alpha, beta, eigenvalue, lanczos_eigenvector, U_t);
copy_double(matrix_size, vector_Q, vector_q_0); shokika_double(matrix_size, vector_Q), shokika_double(101, alpha), shokika_double(100, beta);
lanczos_step_re(matrix_size, search_size, search_gyakubiki, vector_q_0, alpha, beta, lanczos_eigenvector, eigenvector_re, U_t);
copy_double(matrix_size, eigenvector_re, vector_q_0); shokika_double(matrix_size, eigenvector_re); normalization(matrix_size, vector_q_0);
v = variance(matrix_size, search_size, search_gyakubiki, eigenvalue, eigenvector, U_t);
// v = new_variance(matrix_size, search_size, search_gyakubiki, eigenvector, U_t);
// printf("%.20lf\n", v);
z++;
}
copy_double(matrix_size, vector_q_0, eigenvector); shokika_double(matrix_size, vector_q_0);
}
double make_eigenvector_green(int size_length, double eigenvector[size_length], int search[size_length], int search_gyakubiki[size_gyaku], double U_t){
double eigenvalue[100], alpha[101], beta[100], vector_q_0[size_length], kitei_Energy = 0.0;
shokika_double(100, eigenvalue); shokika_double(101, alpha); shokika_double(100, beta); shokika_double(size_length, vector_q_0); r_v_c(size_length, vector_q_0);
lanczos_vector(size_length, search, search_gyakubiki, vector_q_0, alpha, beta, eigenvalue, eigenvector, U_t); kitei_Energy = eigenvalue[0];
return kitei_Energy;
}
double make_kitei_energy_green(int size_length, int search[size_length], int search_gyakubiki[size_gyaku], double U_t){
double eigenvector[size_length], eigenvalue[100], alpha[101], beta[100], vector_q_0[size_length], kitei_Energy = 0.0;
shokika_double(size_length, eigenvector); shokika_double(100, eigenvalue);
shokika_double(101, alpha); shokika_double(100, beta); shokika_double(size_length, vector_q_0); r_v_c(size_length, vector_q_0); lanczos_vector(size_length, search, search_gyakubiki, vector_q_0, alpha, beta, eigenvalue, eigenvector, U_t); kitei_Energy = eigenvalue[0];
return kitei_Energy;
}
int main(){
double U_t = 3.0; int size_gyaku = (int)(pow(2, L));
int SPIN_up = 5, SPIN_down = 5;
// 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;
size_whole_system = combination(Cn * N, 1);
//辞書 normal creation annihilatioin
int search_2D[size], search_gyakubiki[size_gyaku]; shokika_int(size_gyaku, search_gyakubiki);
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_gyakubiki, search_creation_up, search_creation_down, search_annihilation_up, search_annihilation_down, SPIN_up, SPIN_down);
double eigenvector[size], eigenvector_0[size]; shokika_double(size, eigenvector_0); shokika_double(size, eigenvector);
double kitei_Energy_U = make_eigenvector_green(size, eigenvector, search_2D, search_gyakubiki, U_t);
double kitei_Energy_0 = make_eigenvector_green(size, eigenvector_0, search_2D, search_gyakubiki, 0.0);
printf("%lf %lf\n", kitei_Energy_U, kitei_Energy_0);
printf("%lf\n", kitei_Energy_U);
return 0;
}
自分で試したこと
該当コードに対し、コメントアウトや計算結果が配列の外に出ていないかなどを確認しました。
特に関数同士で渡しているMやindexの数値に問題がないかを確認しました。
正直、いろいろ確認しましたがよくわからなくなってしまったため、岡目八目的にアドバイスをいただきたく質問しました。
環境はメモリ８GB、windouws11,WSL上でoneAPIを用いて計算しております。
申し訳ございませんが、よろしければｊ
0
