関数を呼び出して計算しているものと呼び出さないものの計算結果が異なる。
Q&A
解決したいこと
mainの外に関数を定義して計算した場合と、中で定義した計算が大きく異なります。関数の形に大きく違いはなく、代入する変数を変えて計算しています。
発生している問題・エラー
0.000000 0.000000 0.792215
0.000000 0.062210 3.106219
0.000000 0.124420 7.655335
0.000000 0.186629 15.129013
0.000000 0.248839 26.188469
0.000000 0.311049 41.464938
0.000000 0.373259 61.558238
0.000000 0.435468 87.035065
0.000000 0.497678 118.428168
0.000000 0.559888 156.237147
0.000000 0.622098 200.930657
0.000000 0.684307 252.949433
0.000000 0.746517 312.709557
0.000000 0.000000 1.000009
0.000000 0.062210 0.991301
0.000000 0.124420 0.966934
0.000000 0.186629 0.931030
0.000000 0.248839 0.888520
0.000000 0.311049 0.843716
0.000000 0.373259 0.799653
0.000000 0.435468 0.758121
0.000000 0.497678 0.719977
0.000000 0.559888 0.685483
0.000000 0.622098 0.654558
0.000000 0.684307 0.626953
0.000000 0.746517 0.602347
数値がかなり異なります。
該当するソースコード
mainの外
//グラフェンナノリボンエネルギーバンド 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
double b_pi;
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];
double kai[N][N]; //複素感受率
double B[N];
double KKK[P+1];
double static Omega_up[P+1][M][M], Omega_dw[P+1][M][M];
double static U_up_plus[N][P+1][M][M], U_up_minus[N][P+1][M][M];
double static U_dw_plus[N][P+1][M][M], U_dw_minus[N][P+1][M][M];
double static kai_1[N][N][P+1][M][M],kai_1_dw[N][N][P+1][M][M];
double static kai_1_up[N][N][P+1][M][M],kai_2[N][N][P+1][M],kai_3[N][N][P+1];
void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );
double f(double x, int j);
int main(int argc, char** argv){
FILE* fp_1;
FILE* fp_2;
FILE* fp_3;
FILE* fp_4;
FILE* fp_5;
FILE* fp_6;
FILE* fp_7;
FILE* fp_8;
FILE* fp_9;
FILE* fp_10;
b_pi=M_PI;
double U=1.0; //相互作用エネルギーの平均(eV) u/t
if ((fp_1 = fopen("kai.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_9 = fopen("kai_2.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("spinwave.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_3 = fopen("kai_omega0.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_8 = fopen("kaimatrix.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
fp_4 = fopen("u_dw.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_dw[l][i][j]) );
}
}
}
fp_5 = fopen("u_up.txt", "r");
if(fp_5 == 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_5, "%lf", &(u_up[l][i][j]) );
}
}
}
fp_6 = fopen("Energydw_1.txt", "r");
if(fp_6 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_6, "%lf", &(Energy_dw_1[l][i]) );
}
}
fp_7 = fopen("Energyup_1.txt", "r");
if(fp_7 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_7, "%lf", &(Energy_up_1[l][i]) );
}
}
if ((fp_10 = fopen("kai_pi0.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
//二分法
double x_start1 = 1.6;
double x_start2 = 0.0;
double err = 100.0;
double a = x_start1;
double b = x_start2;
double c;
/*while(y<=P){
while(err>0.00000001){
c = (a+b)/2.0;
double fc = f(c,y);
err = fc*fc;
if(fc > U){
b = c;
}
else{
a = c;
}
}
printf("result:x=%lf , y=%lf\n",c,y);
fprintf(fp_2, "%2f %2f\n", c ,2.0*b_pi/(P)*y);
y++;
}*/
for(int i=0;i<R;i++){
double e=x_start1/R*i;
for(int y=0;y<=P;y++){
// double fc=f(e,j);
f(e,y);
fprintf(fp_1, "%2f %2f %2f\n",e,2*b_pi/P*y,kai[0][0]/(P+1));
}
}
fclose(fp_1);
fclose(fp_2);
fclose(fp_3);
fclose(fp_4);
fclose(fp_5);
fclose(fp_6);
fclose(fp_7);
fclose(fp_8);
fclose(fp_9);
fclose(fp_10);
return 0;
}
double f(double x, int j){
for(int h=0;h<N;h++){ //site
for(int k=0;k<N;k++){ //site
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];
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];
kai_1_up[h][k][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[h][k][v][t][s]=U_dw_plus[k][v][t][s]*U_dw_minus[h][v][t][s]/Omega_dw[v][t][s];
kai_1[h][k][v][t][s]=-kai_1_up[h][k][v][t][s]+kai_1_dw[h][k][v][t][s];
kai_2[h][k][v][t]+=kai_1[h][k][v][t][s];
// printf("%2f %2f\n",kai_1[h][k][v][t][s],Omega_dw[v][t][s]);
}
kai_3[h][k][v]+=kai_2[h][k][v][t];
// printf("%2f \n",kai_2[h][k][v][t]);
}
kai[h][k]+=kai_3[h][k][v];
// printf("%2f \n",kai_3[h][k][v]);
}
// printf("%2f \n",kai[j][h][k]);
}
}
return kai[N][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[h][k];
// fprintf(fp_1, "%2f %2f \n",A[N*k+h],kai[i][j][k][h]);
}
}
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];
}
return B[N-1];
}*/
mainのなか
//グラフェンナノリボンエネルギーバンド 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 500
double b_pi;
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];
double static omega[R]; int static q[P+1]; //刻み
double static kai[R][P+1][N][N]; //複素感受率
double static B[R][P+1][N],B_2[R][P+1];
void zheev_( char*, char*, int*, double _Complex*, int*, double*, double _Complex*, int*, double*, int* );
int main(int argc, char** argv){
FILE* fp_1;
FILE* fp_2;
FILE* fp_3;
FILE* fp_4;
FILE* fp_5;
FILE* fp_6;
FILE* fp_7;
FILE* fp_8;
FILE* fp_9;
FILE* fp_10;
b_pi=M_PI;
double U=1.0; //相互作用エネルギーの平均(eV) u/t
if ((fp_1 = fopen("kai.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_9 = fopen("kai_2.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_2 = fopen("spinwave.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_3 = fopen("kai_omega0.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
if ((fp_8 = fopen("kaimatrix.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
fp_4 = fopen("u_dw.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_dw[l][i][j]) );
}
}
}
fp_5 = fopen("u_up.txt", "r");
if(fp_5 == 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_5, "%lf", &(u_up[l][i][j]) );
}
}
}
fp_6 = fopen("Energydw_1.txt", "r");
if(fp_6 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_6, "%lf", &(Energy_dw_1[l][i]) );
}
}
fp_7 = fopen("Energyup_1.txt", "r");
if(fp_7 == NULL) {
printf("ファイルを開くことが出来ませんでした．¥n");
}
for(int l=0; l<=P; l++){
for(int i=0;i<N;i++){
fscanf(fp_7, "%lf", &(Energy_up_1[l][i]) );
}
}
if ((fp_10 = fopen("kai_pi0.txt", "w")) == NULL)
{
printf("Cannot open the file\n");
exit(1);
}
for(int i=0;i<R;i++){
omega[i]=1.6/R*i;
}
for(int i=0;i<=P;i++){
q[i]=i;
}
//感受率
double static Omega_up[R][P+1][P+1][M][M], Omega_dw[R][P+1][P+1][M][M];
double static U_up_plus[P+1][N][P+1][M][M], U_up_minus[P+1][N][P+1][M][M];
double static U_dw_plus[P+1][N][P+1][M][M], U_dw_minus[P+1][N][P+1][M][M];
double static kai_1[R][(P+1)][N][N][P+1][M][M],kai_1_dw[R][(P+1)][N][N][P+1][M][M];
double static kai_1_up[R][(P+1)][N][N][P+1][M][M],kai_2[R][(P+1)][N][N][P+1][M],kai_3[R][(P+1)][N][N][P+1];
for(int i=0;i<R;i++){ //omega
for(int j=0;j<=P;j++){ //q
for(int h=0;h<N;h++){ //site
for(int k=0;k<N;k++){ //site
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[i][j][v][t][s]=omega[i]-Energy_dw_1[v][s]+Energy_up_1[(v+q[j])%(P+1)][t];
Omega_dw[i][j][v][t][s]=omega[i]+Energy_up_1[v][s]-Energy_dw_1[(v+q[j])%(P+1)][t];
U_up_plus[j][k][v][t][s]=conj(u_up[v][t][k])*u_dw[(v+q[j])%(P+1)][s][k];
U_up_minus[j][h][v][t][s]=conj(u_dw[(v+q[j])%(P+1)][s][h])*u_up[v][t][h];
U_dw_plus[j][k][v][t][s]=conj(u_dw[v][t][k])*u_up[(v+q[j])%(P+1)][s][k];
U_dw_minus[j][h][v][t][s]=conj(u_up[(v+q[j])%(P+1)][s][h])*u_dw[v][t][h];
/*U_up_plus[j][k][v][t][s]=(u_up[(v+q[j])%(P+1)][k][t])*u_dw[v][k][s];
U_up_minus[j][h][v][t][s]=(u_dw[v][h][s])*u_up[(v+q[j])%(P+1)][h][t];
U_dw_plus[j][k][v][t][s]=(u_dw[v][k][t])*u_up[(v+q[j])%(P+1)][k][s];
U_dw_minus[j][h][v][t][s]=(u_up[(v+q[j])%(P+1)][h][s])*u_dw[v][h][t];*/
kai_1_up[i][j][h][k][v][t][s]=U_up_plus[j][k][v][t][s]*U_up_minus[j][h][v][t][s]/Omega_up[i][j][v][t][s];
kai_1_dw[i][j][h][k][v][t][s]=U_dw_plus[j][k][v][t][s]*U_dw_minus[j][h][v][t][s]/Omega_dw[i][j][v][t][s];
kai_1[i][j][h][k][v][t][s]=-kai_1_up[i][j][h][k][v][t][s]+kai_1_dw[i][j][h][k][v][t][s];
kai_2[i][j][h][k][v][t]+=kai_1[i][j][h][k][v][t][s];
// printf("%2f %2f\n",kai_1[i][j][h][k][v][t][s],Omega_dw[i][j][v][t][s]);
}
kai_3[i][j][h][k][v]+=kai_2[i][j][h][k][v][t];
// printf("%2f \n",kai_2[i][j][h][k][v][t]);
}
kai[i][j][h][k]+=kai_3[i][j][h][k][v];
}
}
}
}
}
//感受率の対角化
for(int i=0;i<R;i++){
for(int j=0;j<=P;j++){
// 複素行列(エルミート)の対角化 (入力行列の配列は対角化後にユニタリ行列になる)
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*k+h]=kai[i][j][k][h];
// fprintf(fp_1, "%2f %2f \n",A[N*k+h],kai[i][j][k][h]);
}
}
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[i][j][m]=w[m];
}
}
}
for(int i=0;i<R;i++){
for(int j=0;j<(P+1)/2+1;j++){
fprintf(fp_1, "%2f %2f %2f\n",omega[i],2.0*b_pi/(P)*q[j],B[i][j][N-1]/(P+1));
for(int h=0;h<N;h++){
fprintf(fp_8, "%2f %2f %2f %2f\n",omega[i],2.0*b_pi/(P)*q[j],B[i][j][h]/(P+1),kai[i][j][h][h]/(P+1));
if(1/U-0.005<=1*B[i][j][h]/(P+1) && B[i][j][h]/(P+1)<=1/U+0.005){
// if(1/U-0.01<=1*B[i][j][h]/(P+1) && 1*B[i][j][h]/(P+1)<=1/U+0.01){
fprintf(fp_2, "%2f %2f %2f\n",omega[i],2.0*b_pi/(P)*q[j],B[i][j][h]/(P+1));
}
// fprintf(fp_10, "%2f %2f\n",2.0*b_pi/(P)*q[j],B[0][j][h]/N);
// fprintf(fp_3, "%2f %2f\n",omega[i],B[i][0][h]/N);
}
fprintf(fp_3, "%2f %2f\n",omega[i],B[i][0][N-1]/(P+1));
fprintf(fp_10, "%2f %2f\n",2.0*b_pi/(P)*q[j],B[0][j][N-1]/(P+1));
}
fprintf(fp_1, "\n");
fprintf(fp_8, "\n");
}
/*for(int i=0;i<(P+1)/2+4;i++){
for(int j=0;j<R;j++){
if(1/U<=1*B[j][i][N-1]/(P+1) && B[j][i][N-1]/(P+1)<=1/U+0.01){
// if(1/U-0.01<=1*B[i][j][h]/(P+1) && 1*B[i][j][h]/(P+1)<=1/U+0.01){
fprintf(fp_2, "%2f %2f %2f\n",omega[j],2.0*b_pi/(P)*q[i],B[j][i][N-1]/(P+1));
}
}
}*/
for(int i=0;i<R;i++){
for(int s=0;s<=P;s++){
for(int h=0;h<N;h++){
for(int t=0;t<N;t++){
fprintf(fp_9, "%2f %2f %2f\n",omega[i],2.0*b_pi/(P)*q[s],kai[i][s][h][t]/(P+1));
}
}
}
}
/*for(int i=0;i<R;i++){
for(int j=0;j<(P+1)/2+3;j++){
B_2[i][j]=abs(B[i][j][N-1]-1);
}
}
double tmp;
for(int h=0;h<R;h++){
for (int i=0; i<(P+1)/2+3; i++) {
for (int j=i+1; j<(P+1)/2+3; j++) {
if (B_2[h][i] > B_2[h][j]) {
tmp = B_2[h][i];
B_2[h][i] = B_2[h][j];
B_2[h][j] = tmp;
}
}
fprintf(fp_2, "%2f %2f %2f\n",omega[h],2.0*b_pi/(P)*q[i],B_2[h][i]/(P+1));
}
}*/
fclose(fp_1);
fclose(fp_2);
fclose(fp_3);
fclose(fp_4);
fclose(fp_5);
fclose(fp_6);
fclose(fp_7);
fclose(fp_8);
fclose(fp_9);
fclose(fp_10);
return 0;
}
自分で試したこと
現状、必要な計算は下のプログラムでもできていますが、最終的に二分法で評価したいことがあるので配列ではなく変数を連続パラメータのようにしたいです。
コンパイルでは icc -qmkl を用いていますが、コアダンプすることもあり、
icc -qmkl -mcmodel=medium -shared-intel で通しています
下のプログラムは対格化をしていますが、計算結果は対格化する前のもので比較しています。
環境はについて詳しくはわかりませんが
ubuntuでcmd上でコンパイルしています。
0
