概述

本节介绍 cuDensityMat 库的概念和工作原理。 可能需要一些量子力学方面的知识。

通用的库使用工作流程如下：

  1. 创建库上下文句柄。 这将初始化库上下文。 库上下文与当前设置的 GPU 相关联。 所有后续 API 调用和数值运算都必须在同一 GPU 上执行。

  2. （可选）。 为了在多个 GPU 上并行运行，需要通过调用 C API 函数 cudensitymatResetDistributedConfiguration() 并为其提供指向有效 MPI 通信器的指针来显式激活分布式执行。 每个并行进程将被分配其在库上下文创建期间使用的当前 GPU。

  3. 在 GPU 内存中定义必要的基本张量算符的元素。 相应的数组必须使用广义列优先存储布局进行存储。

  4. 定义必要的基本张量算符、算符项和完整量子多体算符。

  5. 创建必要的量子态，查询所需的存储空间，并将 GPU 缓冲区附加到量子态，使其准备好进行计算。 附加的 GPU 缓冲区的原始内容定义了量子态的初始值。 对于多 GPU 运行，库提供了特殊的 API 函数来查询哪个特定的张量切片存储在当前并行进程中。

  6. 准备期望的量子算符对期望种类的量子态（形状、纯度和数据类型）的作用。 这只需要完成一次。

  7. 计算期望的量子算符对具体量子态的作用，该作用将被累积到另一种相同种类的量子态中。

  8. 还可以计算给定量子算符在给定量子态上的期望值。 此计算也需要首先准备（一次），然后根据需要执行多次。

  9. 如果需要，可以对量子态执行不同的基本操作，例如初始化为零、缩放、计算范数、迹、将量子态累积到另一个量子态中以及计算内积。

详细信息如下所示。

张量空间和量子态

张量空间 是一个复合线性空间，构建为一个或多个向量空间的张量积。 每个构成向量空间代表一个量子自由度，其维度等于向量空间的维度。 复合张量空间的形状是构成向量空间的维度元组，按照其原始顺序排列。 请注意，无论我们处理纯态（存在于张量空间中的向量）还是混合态（与同一张量空间关联的矩阵），张量空间的形状都保持不变。

纯量子态 由来自指定张量空间的张量表示，因此纯量子态的形状与张量空间形状一致。 混合量子态 由属于张量空间的张量表示，该张量空间构建为指定张量空间（主空间或右矢空间）及其对偶张量空间（对偶空间或左矢空间）的张量积。 因此，混合量子态的形状是两个张量空间形状的串联，其中张量模式的前半部分对应于右矢空间，而后半部分对应于左矢空间。 例如，如果张量 \(V[i0, i1, i2]\) 表示来自由三个向量空间（量子自由度）构建的张量空间的纯量子态，则张量 \(D[i0, i1, i2; j0, j1, j2]\) 表示来自同一张量空间的混合量子态，其中模式 \(i0, i1, i2\) 是右矢模式，而模式 \(j0, j1, j2\) 是左矢模式。 这两种情况下的对应张量空间形状均为 \([range(i0), range(i1), range(i2)]\)。

在 C API 中，纯/混合量子态由不透明的 cudensitymatState_t 类型表示。 cudensitymatCreateState() API 函数用于从期望形状的张量空间定义纯量子态或混合量子态。 可选地，可以通过指定大于 1 的批次大小来定义一批相同形状的量子态（批处理量子态）。 请注意，cudensitymatCreateState() 期望张量空间形状作为输入参数，而不是量子态形状（对于纯态相同，对于混合态不同）。 创建的量子态尚无关联的物理存储。 为了将物理存储附加到量子态，必须首先通过调用 cudensitymatStateGetNumComponents() 查询量子态表示所需的存储组件数量。 存储组件 是用于存储量子态表示中的特定张量的内存缓冲区。 默认情况下，纯量子态和混合量子态使用单个存储组件存储，该组件使用广义列优先存储布局存储表示量子态的密集张量。 批次维度（模式），如果大于 1，则始终是最后一个（最重要）模式。 相比之下，量子态的因子化表示可以使用多个存储组件，其中每个存储组件存储与请求的张量因子化方案关联的一个（可能批处理的）张量。 当前版本尚不支持因子化张量表示。

一旦确定了所需存储组件的数量，用户可以通过调用 C API 函数 cudensitymatStateGetComponentStorageSize() 查询每个存储组件所需的大小。 了解给定量子态的每个存储组件所需的大小后，用户可以为每个存储组件分配一个 GPU 缓冲区，并通过调用 C API 函数 cudensitymatStateAttachComponentStorage() 将 GPU 缓冲区附加到量子态。 此时，量子态已完全设置，可以用于实际计算。 附加的 GPU 缓冲区必须在量子态的整个生命周期内保持分配状态。 请注意，用于存储量子态表示的 GPU 缓冲区必须是 256 字节对齐的，这在使用 cudaMalloc 或 cudaMallocManaged 时可以保证。 对于自定义 GPU 内存池，必须显式遵守此对齐方式，否则可能会导致未定义的行为。 存储在与给定存储组件关联的附加 GPU 缓冲区中的数据定义了与该存储组件关联的张量的初始值。

对于分布式多 GPU 运行，可以通过调用 C API 函数 cudensitymatResetDistributedConfiguration() 并为其提供 MPI 通信器的副本来激活，张量存储将分布在所有 GPU 上。 分布式执行配置假定每个并行进程一个 GPU，并行进程的总数是二的幂（1、2、4、8、16、32、64 等）。 在这种情况下，密集纯/混合量子态仍将使用单个存储组件（张量）存储，但是每个并行进程将在本地存储完整张量的自身切片。 要查询当前并行进程存储的特定张量切片，可以首先调用 C API 函数 cudensitymatStateGetComponentNumModes() 以检索张量切片模式的数量，然后调用 cudensitymatStateGetComponentInfo() 以检索张量切片信息（模式偏移量、模式范围）。

可以在量子态上执行许多基本操作。 C API 函数 cudensitymatStateInitializeZero() 将量子态初始化为零（将所有张量元素设置为零）。 C API 函数 cudensitymatStateComputeScaling() 将量子态乘以给定的标量（或批处理量子态的标量数组，每个批处理实例一个标量）。 C API 函数 cudensitymatStateComputeNorm() 分别计算纯/混合量子态的平方欧几里得/弗罗贝尼乌斯范数（取平方根将产生范数）。 对于批处理量子态，它将返回一个批次大小的数组。 C API 函数 cudensitymatStateComputeTrace() 计算混合量子态（由密度矩阵表示）的迹，或纯量子态的平方欧几里得范数。 对于批处理量子态，它将返回一个批次大小的数组。 C API 函数 cudensitymatStateComputeAccumulation() 将一个量子态累积（相加）到另一个量子态中，并将可选的缩放系数应用于前者。 对于批处理量子态，它接受批次大小的缩放系数数组。 C API 函数 cudensitymatStateComputeInnerProduct() 计算两个量子态的内积，对于纯量子态，它是常规欧几里得内积，对于混合量子态，它是弗罗贝尼乌斯内积。 对于批处理量子态，它将返回一个批次大小的数组。 在所有量子态为批处理且返回数组的情况下，为结果提供的 GPU 存储空间必须足以存储这种批次大小的数组。

请注意，所有数值运算都接受 CUDA 流参数，也就是说，这些运算的执行通常是异步的，并且需要通过 cudaStreamSynchronize 或 cudaDeviceSynchronize 进行显式同步。 一旦不再需要量子态并且对其执行的所有先前操作都已通过 CUDA 同步 API 同步，则可以通过调用 C API 函数 cudensitymatDestroyState() 安全地销毁它。

量子多体算符（超算符）

广泛的量子 Liouvillian 超算符或广义量子多体算符可以表示为基本张量算符的张量积之和，其中具有一些通常是时间相关的系数，其中基本张量算符是作用于特定数量的量子自由度的张量算符，可以从右侧或左侧或两侧作用。 每个基本张量算符具有相同数量的右矢和左矢模式，其中右矢模式在其张量中位于左矢模式之前。 基本张量算符的左矢模式作用于量子态的右矢模式（从左侧作用），而基本张量算符的右矢模式作用于量子态的左矢模式（从右侧作用）。 实际上，许多量子 Liouvillian 算符都用单体基本张量算符的张量积来表示，例如产生/湮灭算符、粒子数算符、坐标/动量算符、单个自旋算符等。 为了具体起见，让我们考虑以下量子 Liouvillian 算符作为示例：

\[L = \omega_{c} b^{\dagger} b + \frac{1}{2} \omega_{a} \sigma_{z} + [ \frac{1}{2} g(t, \alpha, \beta) b^{\dagger} \sigma_{-} + \frac{1}{2} g(t, \alpha, \beta) b \sigma_{+} ] + \gamma \sigma_{-]\]

在上述 Liouville 算符中，\(b\)、\(b^{\dagger}\)、\(\sigma_{z}\)、\(\sigma_{-}\)、\(\sigma_{+}\) 是基本张量算符，它们各自仅作用于单个量子自由度，由 2 阶张量 \(U[i0; j0]\) 表示，其中 \(i0\) 是右矢模式，\(j0\) 是左矢模式（k 阶基本张量算符由具有 \(k\) 个右矢模式和 \(k\) 个左矢模式的张量表示）。\(\omega_{c}\)、\(\omega_{a}\) 和 \(\gamma\) 是常数标量。\(g(t, \alpha, \beta)\) 是一个随时间变化的标量系数，它取决于时间和两个实参数 \(\alpha\) 和 \(\beta\)，这两个参数参数化了 Liouville 算符。

C API 函数 cudensitymatCreateElementaryOperator() 通过指定它作用的量子自由度（空间模式）的数量、这些空间模式的范围、数值数据类型（实数/复数，FP32/FP64）、包含基本张量算符的张量元素的 GPU 缓冲区（广义列优先布局）以及可选的用户定义的张量回调函数，来创建一个基本张量算符 (cudensitymatElementaryOperator_t)。通过注册用户提供的张量回调函数，我们允许 cuDensityMat 库在计算算符对量子态的作用期间，基于当前时间值以及任意一组实参数自动更新基本张量算符的张量元素（张量回调函数指针 cudensitymatTensorCallback_t 的签名在公共 C 头文件 cudensitymat.h 中提供）。请注意，用户提供的张量回调函数在传递给函数回调的临时 CPU 缓冲区内执行张量元素的更新，然后自动将其传播到在定义基本张量算符期间使用的原始 GPU 缓冲区。一旦不再需要基本张量算符并且任何算符项都不再使用它，就可以通过调用 C API 函数 cudensitymatDestroyElementaryOperator() 安全地销毁它。请注意，存储基本张量算符的张量元素的 GPU 缓冲区必须在基本张量算符的整个生命周期内保持有效。

通常方便地将给定的量子多体算符表示为不同项的总和，每项都有其自身的物理意义，例如系统哈密顿量、脉冲/驱动哈密顿量、耗散项等。例如，上面定义的 Liouville 算符可以被认为是四个算符项组成的（方括号中的算符项包含基本张量算符的两个张量积）。cuDensityMat 库为此目的提供了算符项的概念，由 cudensitymatOperatorTerm_t C 类型表示。C API 函数 cudensitymatCreateOperatorTerm() 仅通过指定算符项作用的量子态所驻留的张量空间形状来创建一个空的算符项。随后，用户可以通过逐步将基本张量算符的张量积与一些通常是随时间变化的系数附加到其中来构建所需的算符项。C API 函数 cudensitymatOperatorTermAppendElementaryProduct() 用于此目的。它接受一个构成所需张量积的基本张量算符数组、张量积作用的量子自由度（空间模式）列表（在其各自的顺序中跨所有构成基本张量算符连接），以及作用于每个空间模式的对偶性。如果对偶性设置为 0，则选择基本张量算符的左矢模式作用于量子态的右矢模式（从左侧作用）；否则，选择基本张量算符的右矢模式作用于量子态的左矢模式（从右侧作用）。最后，C API 函数 cudensitymatOperatorTermAppendElementaryProduct() 还接受与附加的张量积关联的静态系数以及可选的用户定义的系数回调函数 (cudensitymatScalarCallback_t)，该函数也将与附加的张量积关联。cuDensityMat 库将使用用户提供的系数回调函数来检索与附加的张量积关联的总系数的随时间变化的部分，其中总系数是静态系数和系数回调函数对于给定的时间和一组任意实参数在算符作用计算期间返回的值的乘积。例如，张量积 \(\frac{1}{2} g(t, \alpha, \beta) b \sigma_{+}\) 包含两个基本张量算符 \(b\) 和 \(\sigma_{+}\)，一个常数（静态）系数 \(\frac{1}{2}\)，以及一个随时间变化的参数化系数 \(g(t, \alpha, \beta)\)，它可以定义为用户提供的 CPU 回调函数。

一旦不再需要算符项并且任何算符都不再使用它，就可以通过调用 C API 函数 cudensitymatDestroyOperatorTerm() 安全地销毁它。请注意，构成算符项的基本张量算符必须在算符项的整个生命周期内保持有效。

一旦定义了所有必要的算符项，用户就可以继续定义完整的量子多体算符（或超算符），由 cudensitymatOperator_t C 类型表示，通过调用 C API 函数 cudensitymatCreateOperator()。与创建算符项类似，此 API 函数仅通过获取它将作用的量子态所驻留的张量空间的形状来创建一个空算符。随后，用户可以通过 C API 函数 cudensitymatOperatorAppendTerm() 将构成算符项附加到其中来构建所需的算符。每个附加的算符项都可以提供一个静态系数以及一个可选的用户定义的系数回调函数，这样与附加的算符项关联的总系数将是两者的乘积。此外，用户可以请求以其对偶形式附加算符项，其中其所有模式的对偶性将被翻转（右矢 <–> 左矢）。这对于定义主方程的冯·诺依曼部分非常方便，其中哈密顿量从密度矩阵的不同侧作用以形成对易式。一旦不再需要算符并且参与的所有计算都已同步，就可以通过调用 C API 函数 cudensitymatDestroyOperator() 安全地销毁它。请注意，构成算符的算符项必须在算符的整个生命周期内保持有效。

最后，为了将定义的量子多体算符应用于量子态，用户需要首先准备其作用（仅一次），然后根据需要多次计算其作用。C API 函数 cudensitymatOperatorPrepareAction() 准备算符对给定种类的量子态的作用。它将返回实际计算所需的工作区缓冲区大小。然后，用户可以分配所需大小（或更大）的 GPU 工作区缓冲区，并将其附加到调用中使用的工作区描述符 (cudensitymatWorkspaceDescriptor_t)。一旦工作区准备就绪，调用 C API 函数 cudensitymatOperatorComputeAction() 将把算符应用于给定的（输入）量子态，并将其作用累积到另一个（输出）量子态中。此 API 调用采用当前时间值以及一组任意实参数，这些参数将在开始计算之前触发通过注册的回调函数更新算符的值。请求的计算可以插入到 CUDA 流中，并相对于主机异步执行，在这种情况下，用户有责任通过 cudaStreamSynchronize 同步计算。

工作区描述符

cuDensityMat 库支持的每个非平凡数值运算通常都需要在 GPU 内存中分配工作区缓冲区。C API 函数 cudensitymatCreateWorkspace() 创建一个空的工作区描述符。所需数值计算的准备阶段在工作区描述符内设置所需的工作区缓冲区大小。然后，用户可以通过 cudensitymatWorkspaceGetMemorySize() 查询所需大小，并分配至少该大小的 GPU 工作区缓冲区。一旦工作区缓冲区被分配（或从预分配的内存池中获取），就可以通过调用 C API 函数 cudensitymatWorkspaceSetMemory() 将其附加到工作区描述符。此时，工作区描述符已准备好传递给实际的计算调用。在传递到的计算调用已同步后，工作区描述符可以重复使用，前提是工作区缓冲区大小对于下一次计算调用足够。工作区描述符使用的 GPU 缓冲区必须在工作区描述符使用期间保持分配状态。一旦不再需要工作区描述符，就可以通过 cudensitymatDestroyWorkspace() 安全地销毁它。

请注意，附加到工作区描述符的任何 GPU 缓冲区都必须与 256 字节边界对齐，这在使用 cudaMalloc 时得到保证。对于自定义 GPU 内存池，必须显式遵守此对齐方式。

耦合量子动力学主方程（常微分方程组）

虽然 Liouville 超算符及其对量子态的作用足以指定广泛使用的主方程（如 Lindblad 方程）的右侧，但在更一般的设置中，可能需要能够指定耦合常微分方程组 (ODE)，以便在各自的 Liouville 算符下同时演化多个（可能批量处理的）量子态，就像分层运动方程 (HEOM) 方法中所做的那样。在这种情况下，用户需要定义所谓的复合算符作用，由一组 cudensitymatOperatorAction_t 类型的 C 对象表示。它们可以通过调用 C API 函数 cudensitymatCreateOperatorAction() 构建。复合算符作用假设由 \(M\) 个耦合 ODE 系统描述的 \(M\) 个量子态的同时演化。对于 ODE 系统中的每个耦合 ODE，都需要通过调用 cudensitymatCreateOperatorAction() 来构建单独的算符作用对象 (cudensitymatOperatorAction_t)，因此对于完整的 ODE 系统规范需要 \(M\) 个算符作用对象。构建表示耦合 ODE 系统中特定 ODE 的每个算符作用对象接受 \(M\) 个算符，这些算符对应于 \(M\) 个演化的量子态，其中每个算符将作用于其对应的量子态，并累积到该特定 ODE 的右侧。反过来，该特定 ODE 的右侧确定了正在随时间演化的一组 \(M\) 个量子态中对应量子态的时间导数。如果耦合 ODE 系统中的特定 ODE 中不存在特定量子态，则相应的算符在形式上为零，这通过在相应位置提供 NULL 指针来指定。一旦构建了 \(M\) 个算符作用 (cudensitymatOperatorAction_t) 对象，即耦合 ODE 系统已完全指定，就可以通过调用 C API 函数 cudensitymatOperatorActionPrepare() 来准备每个算符作用对象的计算，该函数将在工作区描述符内返回所需工作区缓冲区的大小。为了准备每个算符作用对象，此 API 函数需要知道算符将应用于哪种量子态。准备好后，用户可以查询所需的工作区缓冲区大小，分配所需大小（或更大）的 GPU 工作区缓冲区，并将其附加到调用中使用的工作区描述符 (cudensitymatWorkspaceDescriptor_t)。一旦工作区准备就绪，调用 C API 函数 cudensitymatOperatorActionCompute() 将把所有算符应用于给定的（输入）量子态，并将其作用累积到输出量子态中。此 API 调用采用当前时间值以及一组任意实参数，这些参数将在开始计算之前触发通过注册的回调函数更新算符的值。请求的计算可以插入到 CUDA 流中，并相对于主机异步执行，在这种情况下，用户有责任通过 cudaStreamSynchronize 同步计算。请注意，如果所有构建和准备的算符作用 (cudensitymatOperatorAction_t) 对象都放置在同一个 CUDA 流中，它们将按顺序计算，在这种情况下，只要缓冲区大小对于它们中的任何一个都足够大，就可以为所有这些准备的算符作用重用相同的工作区缓冲区。

请注意，算符作用对象也可以为单个（可能批量处理的）量子态构建，该量子态由单个 Liouville 算符作用，在这种情况下，它等效于直接通过算符 API 准备和计算算符作用。当不再需要算符作用对象时，可以通过调用 cudensitymatDestroyOperatorAction() 安全地销毁它。

量子态的属性

为了计算量子态的不同属性，用户可以创建、准备和计算相应的属性计算对象。C API 函数 cudensitymatCreateExpectation() 接受任意算符 (cudensitymatOperator_t) 并为其创建期望值计算对象 (cudensitymatExpectation_t)。C API 函数 cudensitymatExpectationPrepare()，只需要调用一次，为给定种类的量子态准备期望值计算对象以进行计算。与其他准备函数类似，它在工作区描述符内设置所需的工作区缓冲区大小。然后，用户将需要分配一个该大小（或更大）的 GPU 缓冲区，并将其附加到工作区描述符。准备好后，可以通过 cudensitymatExpectationCompute() 对相同种类（相同形状、纯度和数据类型）的相同或不同量子态多次计算算符的期望值。如果量子态是批量的，cudensitymatExpectationCompute() 将计算批次中所有实例的期望值，从而返回一个期望值数组，该数组大小预计至少为批次大小。

多 GPU 多节点执行

cuDensityMat 库 C API 在单 GPU 和多 GPU 上均可工作。在多 GPU 上运行并行作业需要在下载库存档后构建 MPI 插件。解压缩存档后，可以在 distributed_interfaces 文件夹内找到名为 activate_mpi_cudm.sh 的脚本。通过在该文件夹内运行该脚本，将构建 MPI 接口插件库 libcudensitymat_distributed_interface_mpi.so。该脚本需要环境变量 CUDA_PATH 和 MPI_PATH 分别指向 CUDA 和 MPI 安装目录。提供的 MPI 安装必须指向 CUDA 感知的 MPI 实现。该脚本还将设置环境变量 CUDENSITYMAT_COMM_LIB 以指向 MPI 接口插件库位置。应全局导出此变量，以确保 cuDensityMat 库可以在运行时找到 MPI 接口插件库。

一旦设置了 MPI 接口插件，通过调用 C API 函数 cudensitymatResetDistributedConfiguration() 来激活在多 GPU 上的并行执行。此函数接受 MPI 通信器的副本（以类型擦除的形式），例如通过 MPI_Comm_duplicate 调用获得的 MPI_COMM_WORLD 通信器的副本（请查阅 MPI 库）。在那之后，所有数值计算将在提供的 MPI 通信器内所有可用的 MPI 进程之间并行运行。用户可以通过调用 cudensitymatGetNumRanks() 和 cudensitymatGetProcRank() 分别查询并行进程的总数和当前进程的秩。cuDensityMat 库仅支持并行进程数为 2 的幂次方时的并行执行，每个进程一个 GPU。

引用 cuQuantum

    1. Bayraktar 等人，“cuQuantum SDK：用于加速量子科学的高性能库”，2023 年 IEEE 国际量子计算与工程会议 (QCE)，美国华盛顿州贝尔维尤，2023 年，第 1050-1061 页，doi: 10.1109/QCE57702.2023.00119。