#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;
}
