外观
PETSc 中的 Mat 与 Vec
PETSc(Portable, Extensible Toolkit for Scientific Computation)是高性能科学计算的核心库,其中 Vec(向量)和 Mat(矩阵)是最基础也是最重要的数据类型。本文介绍这两种类型的基本用法,以及如何将它们导出为 MATLAB 可读格式,方便调试和后处理。
Vec:并行向量
Vec 是 PETSc 中表示向量的抽象类型。它支持多种后端实现(MPI 并行、顺序、CUDA 等),用户代码通过统一接口操作,无需关心底层实现。
在 MPI 并行环境中,一个 Vec 被分布在多个进程上,每个进程持有一段连续的局部数据。
创建 Vec
Vec x; // 声明 Vec 变量,此时仅声明指针,未分配内存
PetscInt n = 100; // 全局大小
// 创建标准 MPI 向量(自动均匀分布)
VecCreate(PETSC_COMM_WORLD, &x); // 创建空向量对象,绑定到所有进程;大小和类型尚未确定,x 会指向新创建的对象
VecSetSizes(x, PETSC_DECIDE, n); // 设置大小:本地大小由 PETSc 自动均匀分配,全局大小为 n
VecSetFromOptions(x); // 根据命令行选项确定向量类型(如 VECMPI/VECCUDA 等),完成后向量才可用基本操作
// 设置所有元素为同一值
VecSet(x, 1.0);
// 设置单个元素(全局索引)
PetscInt idx = 5;
PetscScalar val = 3.14;
VecSetValue(x, idx, val, INSERT_VALUES);
// 批量设置元素
PetscInt indices[] = {0, 1, 2};
PetscScalar values[] = {1.0, 2.0, 3.0};
VecSetValues(x, 3, indices, values, INSERT_VALUES);
// 完成赋值(必须调用)
VecAssemblyBegin(x);
VecAssemblyEnd(x);
// 向量加法:y = a*x + y
VecAXPY(y, alpha, x);
// 点积
PetscScalar dot;
VecDot(x, y, &dot);
// 2-范数
PetscReal norm;
VecNorm(x, NORM_2, &norm);
// 释放内存
VecDestroy(&x);访问本地数据
PetscScalar *array; // 声明指针,用于指向向量的底层数据数组
VecGetArray(x, &array); // 获取本进程局部数据的指针;若数据在 GPU 上,PETSc 会先自动拷贝到 CPU
for (PetscInt i = 0; i < local_size; i++) { // local_size 是本进程的元素个数,不是全局大小 n
array[i] = (PetscScalar)i; // 直接操作底层数组,比 VecSetValues 逐个设置更高效
}
VecRestoreArray(x, &array); // 归还指针(必须调用):通知 PETSc 数据修改完毕,同步内部状态;array 被置为 NULLMat:并行矩阵
Mat 是 PETSc 中表示矩阵的抽象类型,支持多种存储格式:
| 格式 | 说明 |
|---|---|
MATAIJ | 压缩稀疏行(CSR),最常用 |
MATBAIJ | 块稀疏行,适合块结构矩阵 |
MATSBAIJ | 对称块稀疏行 |
MATDENSE | 稠密矩阵 |
MATIS | 用于 FETI/DD 方法 |
在并行环境中,矩阵按行分布:每个进程负责若干连续行。
创建 Mat
Mat A;
PetscInt M = 100, N = 100; // 全局行数、列数
// 通用方式(推荐)
MatCreate(PETSC_COMM_WORLD, &A);
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N);
MatSetFromOptions(A);
MatSetUp(A);
// 直接创建 AIJ 格式的 MPI 矩阵
MatCreateAIJ(PETSC_COMM_WORLD,
PETSC_DECIDE, PETSC_DECIDE, // 本地行/列数(自动分配)
M, N, // 全局行/列数
5, NULL, // 对角块每行非零数(估计值)
2, NULL, // 非对角块每行非零数(估计值)
&A);基本操作
// 插入元素(全局行列索引)
PetscInt row = 0, col = 1;
PetscScalar val = 2.5;
MatSetValue(A, row, col, val, INSERT_VALUES);
// 批量插入
PetscInt rows[] = {0, 1}, cols[] = {0, 1};
PetscScalar vals[] = {1.0, 2.0, 3.0, 4.0}; // 按行优先排列
MatSetValues(A, 2, rows, 2, cols, vals, INSERT_VALUES);
// 完成组装(必须调用)
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
// 矩阵-向量乘法:y = A*x
MatMult(A, x, y);
// 矩阵-向量乘法(转置):y = A^T * x
MatMultTranspose(A, x, y);
// 释放内存
MatDestroy(&A);预分配非零结构
稀疏矩阵性能的关键在于预分配。PETSc 提供逐行指定非零数的接口:
// 假设每行最多 5 个非零元素(均匀估计)
MatMPIAIJSetPreallocation(A, 5, NULL, 2, NULL);
// 或者按行精确指定(d_nnz: 对角块,o_nnz: 非对角块)
PetscInt d_nnz[] = {3, 4, 3, ...}; // 每行对角块非零数
PetscInt o_nnz[] = {1, 1, 2, ...}; // 每行非对角块非零数
MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz);导出为 MATLAB 格式
调试大型稀疏线性系统时,将矩阵和向量导出到 MATLAB 非常实用。PETSc 提供了 PetscViewer 机制来实现这一功能。
方法一:导出为 MATLAB 二进制格式(推荐)
PETSc 的二进制格式可直接被 MATLAB 的 PetscBinaryRead 读取。
C 代码:
#include <petscviewer.h>
PetscViewer viewer;
// 导出向量
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vec_output.dat",
FILE_MODE_WRITE, &viewer);
VecView(x, viewer);
PetscViewerDestroy(&viewer);
// 导出矩阵
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mat_output.dat",
FILE_MODE_WRITE, &viewer);
MatView(A, viewer);
PetscViewerDestroy(&viewer);MATLAB 读取:
需要 PETSc 提供的 MATLAB 工具函数,位于 $PETSC_DIR/share/petsc/matlab/,将该目录加入 MATLAB 路径:
addpath('/path/to/petsc/share/petsc/matlab');
% 读取向量
x = PetscBinaryRead('vec_output.dat');
% 读取矩阵(自动识别为稀疏矩阵)
A = PetscBinaryRead('mat_output.dat');
% 验证
spy(A); % 可视化稀疏结构
norm(x) % 查看向量范数方法二:导出为 MATLAB ASCII 格式
对于小矩阵,可以直接生成 .m 脚本文件:
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, "matrix.m", &viewer);
PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
PetscObjectSetName((PetscObject)A, "A"); // 指定导出后在 MATLAB 中的变量名
PetscObjectSetName((PetscObject)x, "x"); // 若不设置,PETSc 会使用默认名称
MatView(A, viewer);
VecView(x, viewer);
PetscViewerPopFormat(viewer);
PetscViewerDestroy(&viewer);生成的 matrix.m 文件可直接在 MATLAB 中运行:
% MATLAB 中执行
run('matrix.m');
% 变量 A 和 x 自动加载到工作区注意: ASCII MATLAB 格式不适合大矩阵,建议仅用于调试小规模问题。
方法三:使用 MATLAB PetscViewer(并行直接写入)
对于并行程序,也可以使用专用的 MATLAB viewer:
PetscViewer viewer;
PetscViewerMatlabOpen(PETSC_COMM_WORLD, "data.mat",
FILE_MODE_WRITE, &viewer);
MatView(A, viewer);
VecView(x, viewer);
PetscViewerDestroy(&viewer);在 MATLAB 中直接用 load 读取:
data = load('data.mat');完整示例
以下是一个组装 1D 泊松方程刚度矩阵并导出的完整示例:
#include <petscmat.h>
#include <petscvec.h>
#include <petscviewer.h>
int main(int argc, char **argv) {
PetscInitialize(&argc, &argv, NULL, NULL);
PetscInt n = 10;
Mat A;
Vec b;
// 创建矩阵和向量
MatCreate(PETSC_COMM_WORLD, &A);
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
MatSetFromOptions(A);
MatSetUp(A);
VecCreate(PETSC_COMM_WORLD, &b);
VecSetSizes(b, PETSC_DECIDE, n);
VecSetFromOptions(b);
// 组装三对角矩阵(1D Laplacian)
for (PetscInt i = 0; i < n; i++) {
PetscScalar val = 2.0;
MatSetValue(A, i, i, val, INSERT_VALUES);
if (i > 0) {
val = -1.0;
MatSetValue(A, i, i - 1, val, INSERT_VALUES);
}
if (i < n - 1) {
val = -1.0;
MatSetValue(A, i, i + 1, val, INSERT_VALUES);
}
}
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
// 设置右端项
VecSet(b, 1.0);
// 导出为 MATLAB ASCII 格式
PetscViewer viewer;
PetscViewerASCIIOpen(PETSC_COMM_WORLD, "poisson.m", &viewer);
PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
PetscObjectSetName((PetscObject)A, "A");
PetscObjectSetName((PetscObject)b, "b");
MatView(A, viewer);
VecView(b, viewer);
PetscViewerPopFormat(viewer);
PetscViewerDestroy(&viewer);
MatDestroy(&A);
VecDestroy(&b);
PetscFinalize();
return 0;
}在 MATLAB 中读取:
run('poisson.m'); % 执行后 A 和 b 自动加载到工作区
% 验证
full(A)
spy(A)总结
| 需求 | 推荐方法 |
|---|---|
| 调试小矩阵 | ASCII MATLAB 格式 |
| 通用导出/读取 | 二进制格式 + PetscBinaryRead |
直接用 load 读取 | PetscViewerMatlabOpen |
PETSc 的 PetscViewer 机制非常灵活,掌握二进制导出方式后,可以方便地在 MATLAB、Python(scipy.io)等环境中对矩阵和向量进一步分析。
版权所有
版权归属:Guisong Wu