为了提高准确性,我一直在尝试将程序从 double 升级到 long double。但是,我收到以下错误。
haread.cpp:178:43: 错误: '2.0e+0 * std::operator+ [with _Tp = long double](( (((const std::complex )KEp) + ( (long unsigned int)(((long unsigned int)i) * 32ul)))), ( (const std::complex )(& energia)
这条线是
phicheb[1][i] = (2.0*(KEp[i] + energia)/dE)-phistate[is][i] ;
数组是使用“新”动态定义的。所有变量都是 long double。当某个数字与复数数组相乘时,所有情况都会出现问题。因此,我将 2.0 更改为 2.0L 以及所有类似的情况。所以,现在编译成功了。但是,我的程序不再正常工作。我是否需要为我定义的所有长双打后缀 L?我将它们定义如下:
长双 xmin = 1.6 ; 长双 xmax = 7.0 ;
我什至不确定从哪里开始。有什么建议么?编译器是 gcc 版本 4.6.2 (SUSE Linux)。编辑:
for ( int is = 0 ; is < nstates ; is++ ){
for ( int i = 0 ; i < xgrid ; i++ ) phi[i] = phistate[is][i];
fft(&phi[0], &kphi[0], -1);
for ( int i = 0 ; i < xgrid ; i++ ){
kphi[i] = akx2[i]*kphi[i]*ixgrid; //calculates the KE using Fourier transform
}
fft(&kphi[0],&KEp[0],1);
for ( int i = 0 ; i < xgrid ; i++ ){
x = xmin + (i*dx) ;
energia = complex<long double>(0.0,0.0);
energia = hmatrix[is][is][i]*phistate[is][i] ; //Potential energy
phicheb[0][i] = phistate[is][i];
phicheb[1][i] = (2.0L*(KEp[i] + energia)/dE)-phistate[is][i] ;
}
for ( int j = 2 ; j < ncheb ; j++ ){ //this is a recursion relation
for ( int i = 0 ; i < xgrid ; i++ ) phi[i] = phicheb[j-1][i];
fft(&phi[0], &kphi[0], -1);
for ( int i = 0 ; i < xgrid ; i++ ){
kphi[i] = akx2[i]*kphi[i]*ixgrid;
}
fft(&kphi[0],&KEp[0],1) ;
for ( int i = 0 ; i < xgrid ; i++ ){
x = xmin + (i*dx) ;
energia = complex<long double>(0.0,0.0);
energia = hmatrix[is][is][i]*phi[i];
phicheb[j][i] = (2.0L*((2.0L*(KEp[i] + energia)/dE) - phicheb[j-1][i])) - phicheb[j-2][i]; //recursion
}
}
for ( int i = 0 ; i < xgrid ; i++){
for ( int j = 0 ; j < ncheb ; j++ ){
phistate_new[is][i] += chebc[j]*phicheb[j][i] ;
}
}
}
这是代码的核心。它用于在网格上传播波函数。我将尝试尽可能少地解释代码。我有 phistate[0][i] 的值,它是对应于基态的波函数(其余状态均为零)。然后我计算系统的能量(使用傅里叶变换的动能和势能已经给出)。现在,使用一种特殊的方法(切比雪夫方法),我及时传播波函数。该方法基于递归关系。在此递归中,使用多项式展开计算传播。Chebc 是 long double 系数,其值定义明确。所有未显示定义的变量都是long double。
当我只使用双而不是长双时,我得到结果没有问题。当使用 long double 时,我的结果 phistate_new 具有 NaN 给出的所有值。