高性能可扩展线性求解算法库Hypre
目录
https://nixbit.com/software/hypre-review/
1 简介
Hypre是由劳伦斯利弗莫尔国家实验室开发并进行开源。劳伦斯利弗莫尔国家实验室(LLNL)成立于1952年冷战高峰时期，旨在通过推进核武器科学和技术来满足紧迫的国家安全需求。在其历史上，实验室通过其才华横溢和敬业的员工队伍和世界一流的研究能力，以科学和技术创新的传统，预测、开发和提供解决国家最具挑战性的问题的解决方案，加强了国家安全。65年来，LLNL一直在创造历史
Livermore的HYPRE线性求解器库可以比传统方法更快地解决问题，从而实现更大，更详细的仿真。 它提供了一套全面的可扩展求解器套件，用于大规模科学仿真，具有针对结构化和非结构化网格问题的并行多网格方法。 HYPRE库具有高度的可移植性，并支持多种语言。
HYPRE的工作始于1990年代后期。 此后，它已被研究机构和私人公司用来模拟地下水流动，托卡马克和恒星中的磁聚变能等离子体，流经心脏的血液，用于核电站的蒸汽发生器中的流体流动以及油库中的泵送活动。 只是几个领域。 2007年，HYPRE荣获R＆D Magazine颁发的R＆D 100奖，这是该年度最重大的技术突破之一。
HYPRE团队是最早为极端规模的并行超级计算机开发代数多重网格算法和软件的团队之一。 该团队在多网格研究社区中保持着积极的作用，并因其在算法和软件开发方面的领导地位而受到认可。
Hypre可扩展线性求解器项目的目标是开发可扩展算法和软件，用于在并行计算机上求解大型稀疏线性方程组。
主要软件产品是hypre，这是一个高性能的预处理器库，具有针对结构化和非结构化网格问题的并行多网格方法。
在LLNL和其他地方正在开发用于研究国防，环境，能源和生物科学中的物理现象的仿真代码时，会引起关注的问题。
尽管对于这些问题的数值解决方案，必须进行并行处理，但仅靠并行处理是不够的。 还需要可扩展的数值算法。 “可伸缩”通常是指有效使用附加计算资源来解决日益严重的问题的能力。 许多因素都会影响可伸缩性，包括并行计算机的体系结构和算法的并行实现。 但是，一个重要的问题经常被忽略：算法本身的可伸缩性。 在这里，可伸缩性是对总计算工作需求如何随问题大小增长的描述，可以独立于计算平台进行讨论。
今天的仿真代码中使用的许多算法都是基于昨天不可扩展的技术。 这意味着解决越来越大的问题所需的工作比线性增长（最佳速率）快得多。 使用可伸缩算法可以将仿真时间减少几个数量级，从而将MPP上的两天运行时间减少到30分钟。 此外，使用该技术的代码仅受机器内存大小的限制，因为它们能够有效地利用其他计算机资源来解决巨大的问题。
可伸缩的算法使应用科学家能够提出并回答新问题。 例如，如果给定的模拟（具有特定的分辨率）需要花费几天的时间来运行，而精炼的（即更准确的）模型将花费更长的时间，那么应用科学家可能会放弃更大，更高保真度的模拟。 由于每次运行都花费太长时间，因此他或她还可能被迫缩小参数研究的范围。 通过减少执行时间，可扩展算法允许科学家以更高的分辨率进行更多的模拟。
此版本中的新增功能：此版本添加了AME，这是一个基于AMS和LOBPCG的新的Maxwell本征求解器。
Python支持已添加到Babel接口。
内部标头中的命名空间问题已修复。
MG求解器模块和BoomerANG模块中的插值算法都得到了加速。
2 安装 hypre¶
如前所述，在大多数系统上，只需在顶级源目录中键入后键入，即可构建 hypre。或者，可以使用 CMake系统 [CMakeWeb]，并且是特别是在 Windows 系统上构建 hypre 的最佳方法。有关更多详细信息，请阅读与 hypre 分发一起提供的文件，或参阅本手册中的最后一章。请注意以下要求：configuremakeINSTALL
要并行运行，hypre 需要安装 MPI。
使用线程配置 hypre 需要实现 OpenMP。目前，只有 hypre 的子集是线程的。
hypre 库目前不直接支持复杂值的系统。
通过vcpkg可以进行安装
vcpkg install hpyre:x64-windows
3 使用Hypre
3.1 选择概念接口
在编写任何代码之前，要做出一项重要的决定是选择适当的概念接口。这些概念接口旨在表示应用程序开发人员自然地思考其线性问题的方式，并为它们提供自然接口，以便将定义其线性系统的数据传递到 hypre 中。从本质上讲，这些概念接口可以被视为帮助用户为 hypre 求解器和预置器构建矩阵数据结构的方便实用程序。图1 的顶行说明了许多概念性接口。通常，概念接口由不同类型的计算网格表示，但也可能使用其他应用程序功能，如几何信息。例如，使用结构化网格的应用程序（如图1中最左侧的接口）通常按差分和网格来查看其线性问题。另一方面，使用非结构化网格和有限元素的应用程序通常从元素和元件刚度矩阵的角度查看其线性问题。最后，最右边的接口是标准线性代数（矩阵行/列）查看线性问题的方式。
hypre 库目前支持四个概念接口，通常给定问题的适当选择相当明显，例如，结构化网格接口显然不适合非结构化网格应用程序。
- 结构化网格系统接口（结构）：此接口适用于其网格由逻辑矩形网格组成的联合，每个网格点都有固定的非零点差分模式。此接口仅支持每个网格点一个未知点。有关详细信息，请参阅结构化网格系统接口（结构）一章。
- 半结构化网格系统接口（结构）：此接口适用于网格大部分为结构化但具有一些非结构化功能的应用程序。示例包括块结构网格、结构化自适应网格优化 （AMR） 应用程序中的复合网格和过套网格。此接口支持每个单元格的多个未知值。有关详细信息，请参阅半结构化网格系统接口（结构）一章。
- 有限元接口 （FEI）：这适用于从有限元离散化形成线性系统的用户。接口反映典型的有限元数据结构，包括元刚度矩阵。虽然这个接口在hypre中提供，但其定义是在其他地方确定的（请发送电子邮件给艾伦·威廉姆斯·威廉·桑迪亚。政府了解更多信息）。有关详细信息，请参阅有限元接口章节。
- 线性代数系统接口 （IJ）：这是传统的线性代数接口。对于其他基于网格的接口不合适的用户，它可以作为最后的手段。它需要更多的工作在用户的一部分，虽然仍然小于建立并行稀疏数据结构。一般解算器和预置器可通过此接口获得，但不需要更多信息的专用解算器。我们的经验是，拥有旧代码的用户，在其中已经具有用于构建特定格式的矩阵的代码，发现 IJ 接口相对容易使用。有关详细信息，请参阅线性代数系统接口 （IJ）一章。
图1 说明概念线性系统接口概念的图形
Hpyre具有三层架构，数据层，线性求解器，然后是线性系统接口。通常，用户应选择与应用程序匹配的最具体的接口，因为这将允许他们使用专用和更有效的解算器和预置器，而不会丢失对更通用解算器的访问权限。例如，图 1 的第二行是一组线性解算器算法。每个线性解算器组都需要用户通过概念接口提供不同的信息。因此，例如，最左侧框中列出的几何多网格算法 （GMG） 只能与最左侧概念接口一起使用。另一方面，最右侧框中的 ILU 算法可与任何概念接口一起使用。每个解算器和预置器的矩阵要求在章节求解器和预置器以及 hypre 参考手册中提供。您所需的解算器策略可能会影响您选择的概念接口。典型用户将选择单个 Krylov 方法和单个预置器来解决其系统。
图1 的第三行是数据布局或矩阵/矢量存储方案的列表。线性解算器与存储方案的关系与概念接口和线性解算器的关系相似。请注意，hypre 中某些接口当前仅支持一个矩阵/矢量存储方案选择。概念接口、所需的解算器和预置器以及矩阵存储类必须全部兼容。
3.2 编写代码
如上一节所述，在编写任何代码之前，应做出以下决定：
- 选择概念接口。
- 选择所需的解算器策略。
- 查找每个解算器和预置器的矩阵要求。
- 选择与解算器、预置器以及概念接口兼容的矩阵存储类。
做出以前的决策后，是时候对应用程序进行编码以调用 hypre 了。此时，查看前面提到的与 hypre 库一起提供的示例代码可能非常有用。示例代码演示以下应用程序调用 hypre 的一般结构：
- 为所选概念界面构建任何必要的辅助结构。例如，如果使用结构化网格接口，包括网格和模具结构。
- 通过所选的概念界面构建矩阵、解决方案向量和右侧向量。每个概念界面都提供一系列调用，用于将有关问题的信息输入到 hypre 中。
- 生成解算器和预置器以及设置解算器参数（可选）。某些参数（如收敛容差）在求解器之间是相同的，而其他参数是特定于解算器的。
- 调用解算器的求解函数。
- 从解算器检索所需信息。根据您的应用程序，您可能想要对解决方案向量进行其他更改。此外，性能信息（如迭代次数）通常可用，但可能与解算器到解算器不同。
本用户手册的后续章节提供了更充分地理解每个概念接口和每个解算器的功能所需的详细信息。请记住，hypre 参考手册中提供了所有可用函数的完整列表，并且提供的示例代码可能证明有助于作为特定应用程序的模板。
4 实例测试
Hpyre的github源代码中提供了18个实例
https://github.com/hypre-space/hypre/tree/master/src/examples
选择其中的第一个进行测试。
/******************************************************************************
* Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
* HYPRE Project Developers. See the top-level COPYRIGHT file for details.
*
* SPDX-License-Identifier: (Apache-2.0 OR MIT)
******************************************************************************/
/*
Example 1
Interface: Structured interface (Struct)
Compile with: make ex1 (may need to edit HYPRE_DIR in Makefile)
Sample run: mpirun -np 2 ex1
Description: This is a two processor example. Each processor owns one
box in the grid. For reference, the two grid boxes are those
in the example diagram in the struct interface chapter
of the User's Manual. Note that in this example code, we have
used the two boxes shown in the diagram as belonging
to processor 0 (and given one box to each processor). The
solver is PCG with no preconditioner.
We recommend viewing examples 1-4 sequentially for
a nice overview/tutorial of the struct interface.
*/
#include <stdio.h>
#include <string.h>
/* Struct linear solvers header */
#include "HYPRE_struct_ls.h"
#ifdef HYPRE_EXVIS
#include "vis.c"
#endif
int main(int argc, char* argv[])
{
int i, j, myid, num_procs;
int vis = 0;
HYPRE_StructGrid grid;
HYPRE_StructStencil stencil;
HYPRE_StructMatrix A;
HYPRE_StructVector b;
HYPRE_StructVector x;
HYPRE_StructSolver solver;
/* Initialize MPI */
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
/* Initialize HYPRE */
HYPRE_Init();
if (num_procs != 2)
{
if (myid == 0) printf("Must run with 2 processors!\n");
MPI_Finalize();
return(0);
}
/* Parse command line */
{
int arg_index = 0;
int print_usage = 0;
while (arg_index < argc)
{
if (strcmp(argv[arg_index], "-vis") == 0)
{
arg_index++;
vis = 1;
}
else if (strcmp(argv[arg_index], "-help") == 0)
{
print_usage = 1;
break;
}
else
{
arg_index++;
}
}
if ((print_usage) && (myid == 0))
{
printf("\n");
printf("Usage: %s [<options>]\n", argv[0]);
printf("\n");
printf(" -vis : save the solution for GLVis visualization\n");
printf("\n");
}
if (print_usage)
{
MPI_Finalize();
return (0);
}
}
/* 1. Set up a grid. Each processor describes the piece
of the grid that it owns. */
{
/* Create an empty 2D grid object */
HYPRE_StructGridCreate(MPI_COMM_WORLD, 2, &grid);
/* Add boxes to the grid */
if (myid == 0)
{
int ilower[2] = { -3,1 }, iupper[2] = { -1,2 };
HYPRE_StructGridSetExtents(grid, ilower, iupper);
}
else if (myid == 1)
{
int ilower[2] = { 0,1 }, iupper[2] = { 2,4 };
HYPRE_StructGridSetExtents(grid, ilower, iupper);
}
/* This is a collective call finalizing the grid assembly.
The grid is now ``ready to be used'' */
HYPRE_StructGridAssemble(grid);
}
/* 2. Define the discretization stencil */
{
/* Create an empty 2D, 5-pt stencil object */
HYPRE_StructStencilCreate(2, 5, &stencil);
/* Define the geometry of the stencil. Each represents a
relative offset (in the index space). */
{
int entry;
int offsets[5][2] = { {0,0}, {-1,0}, {1,0}, {0,-1}, {0,1} };
/* Assign each of the 5 stencil entries */
for (entry = 0; entry < 5; entry++)
HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
}
}
/* 3. Set up a Struct Matrix */
{
/* Create an empty matrix object */
HYPRE_StructMatrixCreate(MPI_COMM_WORLD, grid, stencil, &A);
/* Indicate that the matrix coefficients are ready to be set */
HYPRE_StructMatrixInitialize(A);
/* Set the matrix coefficients. Each processor assigns coefficients
for the boxes in the grid that it owns. Note that the coefficients
associated with each stencil entry may vary from grid point to grid
point if desired. Here, we first set the same stencil entries for
each grid point. Then we make modifications to grid points near
the boundary. */
if (myid == 0)
{
int ilower[2] = { -3,1 }, iupper[2] = { -1,2 };
int stencil_indices[5] = { 0,1,2,3,4 }; /* labels for the stencil entries -
these correspond to the offsets
defined above */
int nentries = 5;
int nvalues = 30; /* 6 grid points, each with 5 stencil entries */
double values[30];
/* We have 6 grid points, each with 5 stencil entries */
for (i = 0; i < nvalues; i += nentries)
{
values[i] = 4.0;
for (j = 1; j < nentries; j++)
values[i + j] = -1.0;
}
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
stencil_indices, values);
}
else if (myid == 1)
{
int ilower[2] = { 0,1 }, iupper[2] = { 2,4 };
int stencil_indices[5] = { 0,1,2,3,4 };
int nentries = 5;
int nvalues = 60; /* 12 grid points, each with 5 stencil entries */
double values[60];
for (i = 0; i < nvalues; i += nentries)
{
values[i] = 4.0;
for (j = 1; j < nentries; j++)
values[i + j] = -1.0;
}
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
stencil_indices, values);
}
/* Set the coefficients reaching outside of the boundary to 0 */
if (myid == 0)
{
double values[3];
for (i = 0; i < 3; i++)
values[i] = 0.0;
{
/* values below our box */
int ilower[2] = { -3,1 }, iupper[2] = { -1,1 };
int stencil_indices[1] = { 3 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
{
/* values to the left of our box */
int ilower[2] = { -3,1 }, iupper[2] = { -3,2 };
int stencil_indices[1] = { 1 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
{
/* values above our box */
int ilower[2] = { -3,2 }, iupper[2] = { -1,2 };
int stencil_indices[1] = { 4 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
}
else if (myid == 1)
{
double values[4];
for (i = 0; i < 4; i++)
values[i] = 0.0;
{
/* values below our box */
int ilower[2] = { 0,1 }, iupper[2] = { 2,1 };
int stencil_indices[1] = { 3 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
{
/* values to the right of our box */
int ilower[2] = { 2,1 }, iupper[2] = { 2,4 };
int stencil_indices[1] = { 2 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
{
/* values above our box */
int ilower[2] = { 0,4 }, iupper[2] = { 2,4 };
int stencil_indices[1] = { 4 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
{
/* values to the left of our box
(that do not border the other box on proc. 0) */
int ilower[2] = { 0,3 }, iupper[2] = { 0,4 };
int stencil_indices[1] = { 1 };
HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
stencil_indices, values);
}
}
/* This is a collective call finalizing the matrix assembly.
The matrix is now ``ready to be used'' */
HYPRE_StructMatrixAssemble(A);
}
/* 4. Set up Struct Vectors for b and x. Each processor sets the vectors
corresponding to its boxes. */
{
/* Create an empty vector object */
HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &b);
HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &x);
/* Indicate that the vector coefficients are ready to be set */
HYPRE_StructVectorInitialize(b);
HYPRE_StructVectorInitialize(x);
/* Set the vector coefficients */
if (myid == 0)
{
int ilower[2] = { -3,1 }, iupper[2] = { -1,2 };
double values[6]; /* 6 grid points */
for (i = 0; i < 6; i++)
values[i] = 1.0;
HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);
for (i = 0; i < 6; i++)
values[i] = 0.0;
HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
}
else if (myid == 1)
{
int ilower[2] = { 0,1 }, iupper[2] = { 2,4 };
double values[12]; /* 12 grid points */
for (i = 0; i < 12; i++)
values[i] = 1.0;
HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);
for (i = 0; i < 12; i++)
values[i] = 0.0;
HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
}
/* This is a collective call finalizing the vector assembly.
The vectors are now ``ready to be used'' */
HYPRE_StructVectorAssemble(b);
HYPRE_StructVectorAssemble(x);
}
/* 5. Set up and use a solver (See the Reference Manual for descriptions
of all of the options.) */
{
/* Create an empty PCG Struct solver */
HYPRE_StructPCGCreate(MPI_COMM_WORLD, &solver);
/* Set some parameters */
HYPRE_StructPCGSetTol(solver, 1.0e-06); /* convergence tolerance */
HYPRE_StructPCGSetPrintLevel(solver, 2); /* amount of info. printed */
/* Setup and solve */
HYPRE_StructPCGSetup(solver, A, b, x);
HYPRE_StructPCGSolve(solver, A, b, x);
}
/* Save the solution for GLVis visualization, see vis/glvis-ex1.sh */
if (vis)
{
#ifdef HYPRE_EXVIS
GLVis_PrintStructGrid(grid, "vis/ex1.mesh", myid, NULL, NULL);
GLVis_PrintStructVector(x, "vis/ex1.sol", myid);
GLVis_PrintData("vis/ex1.data", myid, num_procs);
#endif
}
/* Free memory */
HYPRE_StructGridDestroy(grid);
HYPRE_StructStencilDestroy(stencil);
HYPRE_StructMatrixDestroy(A);
HYPRE_StructVectorDestroy(b);
HYPRE_StructVectorDestroy(x);
HYPRE_StructPCGDestroy(solver);
/* Finalize Hypre */
HYPRE_Finalize();
/* Finalize MPI */
MPI_Finalize();
return (0);
}
编译之后得到hypreTest.exe，通过mpi方式运行
mpiexec -n 2 hypreTest.exe
结果如下
<C*b,b>: 1.800000e+01
Iters ||r||_C conv.rate ||r||_C/||b||_C
----- ------------ --------- ------------
1 2.509980e+00 0.591608 5.916080e-01
2 9.888265e-01 0.393958 2.330686e-01
3 4.572262e-01 0.462393 1.077693e-01
4 1.706474e-01 0.373223 4.022197e-02
5 7.473022e-02 0.437922 1.761408e-02
6 3.402624e-02 0.455321 8.020061e-03
7 1.214929e-02 0.357057 2.863616e-03
8 3.533113e-03 0.290808 8.327628e-04
9 1.343893e-03 0.380371 3.167586e-04
10 2.968745e-04 0.220906 6.997400e-05
11 5.329671e-05 0.179526 1.256215e-05
12 7.308483e-06 0.137128 1.722626e-06
13 7.411552e-07 0.101410 1.746920e-07
5 文档说明
- Introduction
- Structured-Grid System Interface (Struct)
- Semi-Structured-Grid System Interface (SStruct)
- Finite Element Interface
- Linear-Algebraic System Interface (IJ)
- Solvers and Preconditioners
- General Information
- References
- API
