外观
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