# YH-DIAG

**Repository Path**: e-level-parallel-algorithm/yh-diag

## Basic Information

- **Project Name**: YH-DIAG
- **Description**: YH-DIAG是基于Hubbard模型的大尺寸量子多体系统并行精确对角化计算软件，由国防科技大学的李彪负责研发。YH-DIAG的目的是为了解决增加量子多体的系统尺寸时，随着所需的内存暴增，导致程序设计困难的问题。YH-DIAG仅从Hubbard-Hamiltonian矩阵表示自身出发，利用其稀疏性和对称性，提出了一种快速查找Hamiltonian矩阵非零元的算法，通过即时计算获取矩阵非零元素的策略来减少内存占用，提出了负载均衡的分布式Lanczos过程中矩阵向量乘法的并行化方案。截止到目前，YH-DIAG可以并行扩展到36格点730-billion-dimensional 一维Hamiltonian矩阵精确对角化求解，在资源允许的条件下此计算尺寸还可以增加。YH-DIAG软件采用Fortran77+MPI编程实现，支持双精度浮点运算。
- **Primary Language**: FORTRAN
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2021-02-24
- **Last Updated**: 2021-06-19

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# YH-DIAG

#### 介绍
YH-Diag是基于Hubbard模型的大尺寸量子多体系统并行精确对角化计算软件，由国防科技大学的李彪负责研发。YH-Diag的目的是为了解决增加量子多体的系统尺寸时，随着所需的内存暴增，导致程序设计困难的问题。YH-Diag仅从Hubbard-Hamiltonian矩阵表示自身出发，利用其稀疏性和对称性，提出了一种快速查找Hamiltonian矩阵非零元的算法，通过即时计算获取矩阵非零元素的策略来减少内存占用，提出了负载均衡的分布式Lanczos过程中矩阵向量乘法的并行化方案。截止到目前，YH-Diag可以并行扩展到36格点730-billion-dimensional 一维Hamiltonian矩阵精确对角化求解，在资源允许的条件下此计算尺寸还可以增加。YH-Diag软件采用Fortran77+MPI编程实现，支持双精度浮点运算。


#### 一、软件安装及使用


大规模量子多体系统精确对角化软件采用MPI + fortran 77语言编写。

1.2.1编译和安装

mpif77/mpiifort -O2 -o diag diag.f 

程序采用-O2编译选项当前工作目录下生成可执行程序diag。

1.2.2 运行

在当前工作目录下直接运行命令：

$mpirun -np N ./diag   # N表示并行CPU核数。


 二、软件结构与组成

本软件包含的主要文件，相关说明如下：

uvj.h          头文件包含系统规模以及定义全局数组变量

iuvj_32_08     32格点8个电子系统规模的输入文件

inbond_32      32格点规模的全局数组变量的输入文件

diag.f          源文件

iuvj_32_08和inbond_32两个输入文件可以通过修改规模参数来改变系统规模。具体参数说明如图所示：


![图1 iuvj_32_08输入文件格式与说明](https://images.gitee.com/uploads/images/2021/0619/150132_14c593c4_8717528.png "屏幕截图.png")


 三、软件功能介绍

3.1 Lanczos算法

由于希尔伯特空间维数会随着格点数和粒子数呈指数增长，完整对角化哈密顿矩阵的方法只适用于格点数和粒子数都非常小的系统。通常情况下， 我们只需要求解哈密顿矩阵的最小的特征值和特征向量，所以我们可以通过数值方法来对角化哈密顿矩阵，其中最常用的就是Lanczos算法。Lanczos算法的思想是将Hamiltonian投射到Krylov子空间的一组正交基上：

![输入图片说明](https://images.gitee.com/uploads/images/2021/0619/150302_156c5793_8717528.png "屏幕截图.png")

其中b是一个随机的起始向量，m是子空间的维度，Lanczos过程的结果是得到一个三对角矩阵。Krylov空间的维度不需要事先知道，相反，我们可以在每次迭代后检查收敛性，例如通过计算上一次迭代的基态能量差。在m≪dimH的情况下，最低(和最高)特征值已经准确地逼近H的相应特征值。

![图2 Lanczos算法](https://images.gitee.com/uploads/images/2021/0619/150320_c8234094_8717528.png "屏幕截图.png")

这就是Lanczos算法的处理过程，通过对Krylov子空间的基向量进行正交化处理生成一组Lanczos正交基q1, q2, ..., qm。然后把H投影到这组正交基上生成三对角矩阵 

![输入图片说明](https://images.gitee.com/uploads/images/2021/0619/150343_2891c629_8717528.png "屏幕截图.png")

随着规模增加，对内存的需求急剧增大，我们已经不能单个服务器来存储算法1中的向量q和b，以及哈密顿矩阵H。为了解决此问题，我们需要引入分布式存储并行lanczos算法。

3.2 分布式矩阵向量乘 

当我们考虑更大的规模时，例如32或36格点规模，此时Lanczos算法实现过程成中我们遇到的最大的麻烦是算法1中第8行的矩阵乘向量运算。当规模很小时，我们可以完全存储Hamiltonian矩阵，并按照通用的矩阵向量乘法进行计算，然而当格点数规模超过16时，我们的哈密顿矩阵将变得非常大，单个的服务器已经不足以完全存储矩阵。我们需要考虑通过分布式的方式来存储哈密顿矩阵，但是当格点数超过16格点时，我们也不能轻易的完全存储哈密顿矩阵。我们需要重新设计Hamiltonian矩阵的存储形式，哈密顿矩阵可以某种稀疏矩阵格式存储在内存中[7]，也可以在每次执行矩阵向量乘法时生成，我们采用后一种方法，在每次执行矩阵向量乘法计算之前生成对应的矩阵非零元素。

假设我们在计算Hb的过程中使用了P个进程，进程编号分别为0,1, ... , P − 1。因为进程数可以由用户任意选取，我们考虑选取P满足整除dimH，相应的对矩阵H和向量b，q采取按行的方式均匀的划分为P个块矩阵和P个向量

![输入图片说明](https://images.gitee.com/uploads/images/2021/0619/150434_668f2fb4_8717528.png "屏幕截图.png")

则每个块矩阵和块向量的行数都为dim(H)/P. 对矩阵我们并不进行存储，上述对矩阵按行划分给对应的进程只是形式上的划分，为了使方法阐述的更清晰。全局向量b、q分布式存储到对应进程，进程p存储局部向量b(p)、q(p)。由哈密顿算法的性质可得到哈密顿矩阵的每行的非零元个数不超过2N（N是格点个数数）因此我们对格点进行遍历，对于确定的格点位置，分别考虑自旋向上和自旋向下两种情况.

![图2简略版本的并行矩阵向量乘算法](https://images.gitee.com/uploads/images/2021/0619/150452_af58a5a4_8717528.png "屏幕截图.png")


算法2给出了形式上简略的并行矩阵向量乘Hb的计算流程，首先对格点位置从1到N进行遍历，每次的外循环迭代中当前进程遍历的计算本地存储的行所对应的非零元，对于固定的格点位置k，每一行的非零元不超过两个，分别对自旋向上和自旋向下两个情形，我们分开进行计算。由于是并行算法，导致算法2需要克服一个难点，那就是bjk, b'jk可能不存储在本地进程，此时需要涉及到通信，如何合理设计通信模式。

3.3 并行通信模式设计

并行通信主要是算法2的第4和6行，虽然可以在本地计算出H(p)i,jk和H(p)i,j'k，但是bjk和bj'k可能不存储在本地进程，需要存储它们的进程把数据发送给当前进程。不过此通信存在一个困难，由于jk, j'k是数据接收方计算得到的，这就导致存储bjk和bj'k的发送方进程并不掌握任何通信信息包括数据索引和目的进程。所以无法简单的给数据接收和发送方建立点对点通信。所以为了解决此困难，我们必须要考虑哈密顿矩阵的对称性，由于自旋向上和自旋向下两种情形的计算是类似的，所以只需考虑算法2的第4行。

![图3 哈密顿矩阵的对称性](https://images.gitee.com/uploads/images/2021/0619/150831_aa6926d6_8717528.png "屏幕截图.png")

![图4 分层通信模式](https://images.gitee.com/uploads/images/2021/0619/150848_22f41935_8717528.png "屏幕截图.png")