H-66

分数：200

提交次数限制：0/25

题面

提交

附件

（如果你急着做题，请跳到“任务描述”部分）

背景

室温超导在推动能量的无损传输、磁场应用以及人类在材料和能源方面的研究等方面有着非常重要的前景。H星球上的居民将其视为通往无限能源时代最重要的钥匙。你的朋友小金对这个话题也非常感兴趣，每天都泡在实验室里，试图合成室温超导材料。

有一天，小金通过某种奇妙的方法合成出了一种似乎有抗磁性的材料。他非常激动，向你炫耀这一成果，并表示，如果你投资1亿H币用来和他一起大规模生产这种材料，他会用你的名字小何命名这种材料，叫它“H-66”。

作为理论物理的学渣，你当然不会按他的方法进行材料合成。不过小金给了你一份他测出来的样品的结构，你或许可以理论上算看看这个结构是否是超导体。你并不会什么复杂的方法，于是查资料找到了一种叫Hubbard Model的模型。这似乎是一种简单的方法，能帮你获取样品的一些性质。

计算背景

Hubbard Model 指的是固体物理学中的一种简化模型，用于研究电子在晶格中的行为，特别是在强电子关联系统中。这个模型最初是由物理学家John Hubbard在20世纪60年代提出的，主要用于理解固体中的磁性和超导性。

在Hubbard模型中，电子可以在一个简化的原子格点系统中跃迁，同时受到两个主要因素的影响：

跃迁项：描述电子从一个格点跃迁到另一个格点的能力，通常与电子的动能相关。

相互作用项：描述当两个电子占据同一格点时的库仑排斥能。

这个模型的重要之处在于其简单性和灵活性。虽然它是高度简化的，但仍然能够捕捉到复杂电子系统的一些关键特性，如绝缘体与金属的转变、磁性的形成，以及高温超导体的某些性质。

例如，可以建立一个一维的哈伯德模型，并计算携带不同动量的电子（黑色）或空穴（白色）的激发谱，如下图所示。电子和空穴之间的能隙告诉我们它是一种绝缘体而非导体。

在量子多粒子系统的理论和光谱测量的理论描述中，如角分辨光电子能谱（ARPES）和X射线吸收谱（XAS），常常会出现像 ⟨i∣1ω−H∣i⟩\bra i \frac 1{\omega-H} \ket i⟨i∣ω−H1​∣i⟩ 这样的表达式。ARPES 直接与单粒子激发的谱函数相关，它简要描述了电子结构。而 XAS 则涉及到核心层到价带的激发，并且对元素的价态和轨道占据敏感1^11。现在我们关注一种特定类型的 A(ω)=−Im(⟨x∣1ω+iΓ−H∣x⟩)A(\omega)=-Im(\bra x \frac 1{\omega+i\Gamma-H} \ket x)A(ω)=−Im(⟨x∣ω+iΓ−H1​∣x⟩)，其中 H 是一个厄米算符。如果我们用常见的基来表示这个系统，向量 x 和厄米矩阵 H 的维数会非常大，以至于不能直接计算这个表达式。另一种方法是使用迭代方法为每个 ω\omegaω 解线性方程，但这仍然太慢。一个更实用的方法是尝试将 H 分解为几个特征值和特征向量，H≈∑vieivi′H \approx \sum v_ie_iv_i'H≈∑vi​ei​vi′​，那么 A(ω)≈∑−Im(1ω+iΓ−ei)∣vi′x∣2A(\omega)\approx\sum -Im(\frac 1{\omega+i\Gamma-e_i}) {\lvert v_i'x \rvert}^2A(ω)≈∑−Im(ω+iΓ−ei​1​)∣vi′​x∣2。

我们希望优化对角化计算。这里将使用一个具有随机跃迁项的一维Hubbard Model， H 是一个实对称矩阵。

任务描述

附件中包含样例程序base.cpp，你需要优化样例程序从开始到结束的完整运行时间。该程序读取一个厄米算符（实数，每项的最大长度为 4）和一个向量，并输出一个由该向量展开的 Krylov 子空间中的对称三对角矩阵。程序主要包含以下部分。并非每个部分都值得优化，你需要进行程序热点分析。

getss：读取占据数配置并确定哪些本征态将被包含在基中。

readop：以二次量子化形式读取算符。

act：生成算符的稀疏矩阵。

getsp：实对称稀疏矩阵的 Lanczos 过程

其中，给定一个对称矩阵和一个起始向量，Lanczos 过程构建一个嵌套基和在这个基中的三对角矩阵。迭代步骤如下所定义。Lanczos 过程具有固有的数值不稳定性，不同求和顺序得到的三对角矩阵将完全不同。好在我们真正关心的能谱是稳定的，至少在这个模型中。

inputs: mmm, xxx

x=x/norm(x)x=x/norm(x)x=x/norm(x)

v0=0v_0=0v0​=0

v1=xv_1=xv1​=x

for i=1,2,...\quad i=1,2,...i=1,2,...

x=m∗x\quad x=m*xx=m∗x

ai=vi′∗x\quad a_i=v_i'*xai​=vi′​∗x

x=x−ai∗vi−bi−1∗vi−1\quad x=x-a_i*v_i-b_{i-1}*v_{i-1}x=x−ai​∗vi​−bi−1​∗vi−1​

bi=norm(x)\quad b_i=norm(x)bi​=norm(x)

x=x/bi\quad x=x/b_ix=x/bi​

vi+1=x\quad v_{i+1}=xvi+1​=x

输入输出

输入输出函数在base.cpp中已经给出，请不要改动。这里说明输入输出数据的格式，不过你并不需要关心。

输入：二进制文件 conf.data，通过readss和readop函数进行读取。

输出：二进制文件out.data，通过fwrite(result.data(), 1, 16 * itn, file)输出。

同时，评测系统的输入输出使用非常高速文件系统(> 10GB/s)，不需要优化这一部分。

样例

本题提供两个样例输入文件conf_12.data和conf_14.data以及对应的参考能谱sp_12.data和sp_14.data。后者与测评时的输入具有一样的规模，前者只有约1/12的规模。

为了大家方便验证程序输出是否正确，我们提供一个python脚本，check.ipynb，可以计算最终能谱并与参考能谱比较，这个脚本中的算法也会应用在正式测评中，是目前唯一可行的验证方式。仅通过减少迭代次数加快计算的方式将视为作弊，尽管其得到的结果有概率通过验证。

评测细节

将优化后的单个cpp文件提交。提交的代码会以最简单的方式编译并执行，编译选项为-O3 -fopenmp -march=native，通常情况下编译器会尝试avx2向量化，但你仍可以使用avx512内置函数。如果需要链接特定系统或编译器自带库，请联系工作人员并说明理由，编译条件的更改会对所有选手可见。一共两个规模一样的测试点，各100分。程序运行在一整台intel Xeon 8358 * 2的机器上，样例程序大约需要8分钟。每个测试点优化到2分钟内能够拿到10分基本分，优化到20秒内能够拿到100分满分。分数随运行时间的倒数线性变化。

提示

除了Lanczos算法的核心逻辑外的每一个部分都可以自行修改，包括基组顺序，双向映射实现方式，稀疏矩阵的生成方式，稀疏矩阵与向量的储存形式。

有大量现成的类似算法，可以自行搜索exact diagonalization。

1: arXiv:cond-mat/0505214