ソースちゃんと読んではいませんが、kai_eleをループ毎に生成して毎回ひとつの要素にしか代入してないように見えるのですが大丈夫でしょうか？
Like!
二分法について、うまくいかない
解決したいこと
二分法を用いてgc＝1.00となるｃとiの組を見つけたい。
発生している問題・エラー
現状二分法はできている。
しかし、やりたいことはiごとにｃの値を計算してほしいのだが、特定のｃにおける値しか返してくれない。
該当するソースコード
//グラフェンナノリボンエネルギーバンド 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 100
void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );
double f(double x, int j, int h, int k);
double g(double kai_ele[N][N]);
double u_up[P+1][N][N], u_dw[P+1][N][N]; //固有ベクトル
double u_up_1[P+1][N][N], u_dw_1[P+1][N][N];
double Energy_dw_1[P+1][N],Energy_up_1[P+1][N];
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.6;
double x_start2 = 0.0;
for(int i=0;i<(P+1)/2+3;i++){
double err = 100.0;
double a = x_start1;
double b = x_start2;
double c;
// int i=8;
while(err>0.00000001){
c = (a+b)/2.0;
double gc=0.0;
for(int h=0;h<N;h++){
for(int k=0;k<N;k++){
double kai_ele[N][N];
kai_ele[h][k]=f(c,i,h,k);
gc=g(kai_ele)/(P+1);
}
}
err = gc*gc;
if(gc > U){
b = c;
}
else{
a = c;
}
// printf("result:x=%lf , q=%lf\n",c,2.0*b_pi/(P+1)*i);
// fprintf(fp, "%2f %2f\n", c ,2.0*b_pi/(P+1)*i);
}
// printf("result:x=%lf , q=%lf\n",c,2.0*b_pi/(P+1)*i);
fprintf(fp, "%lf %lf\n", c ,2.0*b_pi/(P+1)*i);
}
/*for(int i=0;i<R;i++){
double e=x_start1/R*i;
for(int j=0;j<(P+1)/2+3;j++){
for(int h=0;h<N;h++){
for(int k=0;k<N;k++){
// double fc=f(e,j);
double kai_ele[N][N];
kai_ele[h][k]=f(e,j,h,k);
gc=g(kai_ele);
}
}
// fprintf(fp, "%2f %2f %2f\n",e,2*b_pi/P*j,gc/(P+1));
if(1/U-0.005<=gc/(P+1) && gc/(P+1)<=1/U+0.005){
fprintf(fp, "%2f %2f %2f\n",e,2*b_pi/P*j,gc/(P+1));
}
}
// fprintf(fp, "%2f %2f\n",e,gc/(P+1));
}*/
fclose(fp);
fclose(fp_1);
fclose(fp_2);
fclose(fp_3);
fclose(fp_4);
return 0;
}
double f(double x, int j, int h, int k){
double kai=0;0;
double kai_2[P+1][M]={0.0};
double kai_1[P+1][M][M]={0.0},kai_1_dw[P+1][M][M]={0.0},kai_1_up[P+1][M][M]={0.0};
double kai_3[P+1]={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
double Omega_up[P+1][M][M], Omega_dw[P+1][M][M];
double U_up_plus[P+1][M][M], U_up_minus[P+1][M][M];
double U_dw_plus[P+1][M][M], U_dw_minus[P+1][M][M];
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];
U_up_plus[v][t][s]=conj(u_up[v][t][k])*u_dw[(v+j)%(P+1)][s][k];
U_up_minus[v][t][s]=conj(u_dw[(v+j)%(P+1)][s][h])*u_up[v][t][h];
U_dw_plus[v][t][s]=conj(u_dw[v][t][k])*u_up[(v+j)%(P+1)][s][k];
U_dw_minus[v][t][s]=conj(u_up[(v+j)%(P+1)][s][h])*u_dw[v][t][h];
kai_1_up[v][t][s]=U_up_plus[v][t][s]*U_up_minus[v][t][s]/Omega_up[v][t][s];
kai_1_dw[v][t][s]=U_dw_plus[v][t][s]*U_dw_minus[v][t][s]/Omega_dw[v][t][s];
kai_1[v][t][s]=-kai_1_up[v][t][s]+kai_1_dw[v][t][s];
kai_2[v][t]+=kai_1[v][t][s];
// printf("%2f\n",kai_1[v][t][s]);
}
kai_3[v]+=kai_2[v][t];
// printf("%2f \n",kai_2[v][t]);
}
kai+=kai_3[v];
// printf("%2f \n",kai_3[v]);
}
return kai;
}
double g(double kai_ele[N][N]){
double B[N];
//感受率の対角化
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
char jobz = 'V' ;// 固有ベクトルを計算する
char uplo = 'U' ;// Aに上三角行列を格納
int n = N ;// 対角化する正方行列のサイズ
double _Complex A[N*N] ;// 対角化する行列。対角化後は固有ベクトルが並ぶ
for(int h=0;h<N;h++){
for(int k=0;k<N;k++){
A[N*h+k]=kai_ele[h][k];
}
}
// 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=B[N-1];
return kai_acc;
}
自分で試したこと
出力する場所を変えたりしてみたがうまく行かない。
環境についてはよくわかりませんが、
Ubuntu 20.04.5 LTS
でコンパイルには icc -qmkl -mcmodel=medium -shared-intel
を用いています。
0
2Answer
1
Comments
@kisara11235
Questioner要素としてはkai_eke[N][N]で対格化してN個要素がありますが、実行的に必要なものは一つだけなので大丈夫だと思います。
配列を渡した先ではすべての要素を使用しているように見えます
@kisara11235
Questionerそうですね、関数ｆでの計算ではｈ，ｋの要素からなるN*N行列を関数gに渡しています。
そして関数gで体格化して一番大きい固有値を渡して計算するように組んでいるつもりです。ここで見当違いなことを言っているとしたら申し訳ないです。
最初の回答に戻っちゃうのですが、代入していない値を使用しているように見えたので不定値入ってたらまずくないですかと思ったまでです
double gc[N][N];
for(int i=0;i<(P+1)/2+3;i++){
//省略...
while(err>0.00000001){
c = (a+b)/2.0;
for(int h=0;h<N;h++){
for(int k=0;k<N;k++){
double kai_ele[N][N];
kai_ele[h][k]=f(c,i,h,k);
gc[h][k]=g(kai_ele)/(P+1);
}
}
//省略...
}
for(int h=0;h<N;h++){
for(int k=0;k<N;k++){
fprintf(fp, "%lf %lf %lf %lf\n", c ,2.0*b_pi/(P+1)*i, h, k, gc[h][k]);
}
}
}
この辺でしょうか。参考までに。
Like!
Comments
@kisara11235
Questioner考えているかつほしい値はｈ，ｋで作る行列全部ではないのでこのようなプログラムではどうしようもないです。
そもそも関数gは対格化を行っているので複数の値を計算しますが、N*N行列にはなりません。
また戻り値も行列では定義していません。