只求完工而不求代码效率,现在终于写完这个程序了。接下来就要做代码的优化和并行化了——事实上需要优化的内容简直太多。
最后求解本征能量与本征函数的时候使用了 intel Math Kernel Library 提供的函数,即 lapack 中的 LAPACKE_dsyevd 函数。编译的时候需要正确调用 intel MKL 才可以编译完成。
/*
luozhen, 2015-05-25
Prof. Cgliu's Group, itcc, NJU
intel MKL linked as: -L${MKLROOT}/lib/intel64 -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -liomp5 -lpthread -lm
*/
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <string>
#include <cstring>
#include <time.h>
#include <cstdlib>
#include <vector>
#include <iomanip>
#include "mkl_lapacke.h"
using namespace std;
//debug
//定义class的方式效率并不高
class config_wave
{
public:
int *alpha;
int *beta;
};
//--end debug
//根据输入文件名确定输出文件名
string FileName(string name)
{
int site = name.find_last_of('.');
if (site == 0)
{
name = "output";
}
else if(site > 0)
{
name.erase(site, name.size() - site + 1);
}
return name;
}
//获取当前系统日期与时间
string _local_time()
{
char buf[64];
time_t curtime = time(NULL);
tm *ptm = localtime(&curtime);
sprintf(buf, "%d/%02d/%02d %02d:%02d:%02d", ptm->tm_year+1900, ptm->tm_mon+1,
ptm->tm_mday, ptm->tm_hour, ptm->tm_min, ptm->tm_sec);
string _curtime(buf);
return _curtime;
}
//获取当前主机的名称
string HostName()
{
string hostname;
//对于 Linux 系统
ifstream hostname_file("/etc/hostname", ifstream::in);
hostname_file >> hostname;
return hostname;
}
//计算组合数的函数,递归方法实现
//int类型时计算上限是C(32,16)
int combinatorial(int n, int k)
{
if(n < 0 || k < 0)
{
cerr<<"Error when calculating the combinatorial number of ("
<<n<<", "<<k<<")!"<<endl;
exit(10);
}
if(k > n)
{
return 0;
}
else if(k > n/2)
{
k = n - k;
}
else if(k == 1)
{
return n;
}
else if(k == 0)
{
return 1;
}
else
{
return combinatorial(n-1, k-1) + combinatorial(n-1, k);
}
}
//由索引值生成01序列
//因为更高效的bitset不能在运行时动态生成,所以只能退而求其次使用数组
int* GetBits(int index, int atoms)
{
int *bits;
bits = new int [atoms];
for(int i = 0; i != atoms; ++ i)
{
bits[i] = 0;
}
int n = atoms - 1;
int k = atoms / 2 - 1;
int flag = combinatorial(n, k);
for(int i = 0; i != atoms; ++ i)
{
if(index > flag)
{
index -= flag;
-- n;
}
else
{
bits[i] = 1;
-- n;
-- k;
}
if(k < 0)
{
break;
}
flag = combinatorial(n, k);
}
return bits;
}
//由组态生成索引,对单种自旋的电子有效
int GetIndex(int *_bits, int atoms)
{
vector<int> _site;
_site.push_back(0);
for(int ix = 0; ix !=