使用SLEPc并行输出特征向量至文件的优化及跨语言读取方案问询

一、更优的特征向量输出方式 #

1. 使用二进制格式输出（推荐） #

二进制格式避免ASCII的并行输出混乱，效率更高，且PETSc原生支持并行读写。修改C代码如下：

PetscViewer viewer;
// 打开二进制文件写入器
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "./data/eigvecs.bin", FILE_MODE_WRITE, &viewer);
// 输出特征向量
EPSVectorsView(eps, viewer);
// 释放资源
PetscViewerDestroy(&viewer);

后续可以通过PETSc命令行工具petscvecview直接查看，也能通过petsc4py（Python）或Petsc.jl（Julia）的接口直接读取完整向量，无需手动处理并行分区。

2. 单进程收集后统一输出ASCII #

如果必须使用ASCII格式，可让主进程（如进程0）收集所有进程的向量数据后再写入文件，保证输出是完整连续的特征向量：

Vec *vecs;
PetscInt n, n_total, i;
// 获取特征向量数组和数量
EPSGetVectors(eps, &vecs);
EPSGetDimensions(eps, &n, NULL, NULL);
// 获取向量总长度
VecGetSize(vecs[0], &n_total);

for (i = 0; i < n; i++) {
    Vec full_vec;
    // 创建全局向量（仅进程0需要）
    if (PetscRank() == 0) {
        VecCreate(PETSC_COMM_SELF, &full_vec);
        VecSetSizes(full_vec, PETSC_DECIDE, n_total);
        VecSetFromOptions(full_vec);
    }
    // 将分布式向量收集到进程0的全局向量中
    VecScatter scatter;
    VecScatterCreateToAll(vecs[i], &scatter, &full_vec);
    VecScatterBegin(scatter, vecs[i], full_vec, INSERT_VALUES, SCATTER_FORWARD);
    VecScatterEnd(scatter, vecs[i], full_vec, INSERT_VALUES, SCATTER_FORWARD);
    VecScatterDestroy(&scatter);

    // 进程0写入文件
    if (PetscRank() == 0) {
        PetscViewerASCIIOpen(PETSC_COMM_SELF, "./data/eigvecs.txt", &viewer);
        VecView(full_vec, viewer);
        PetscViewerDestroy(&viewer);
        VecDestroy(&full_vec);
    }
}

二、现有并行ASCII文件的读取方法 #

Python实现 #

需提前知晓进程数n_procs、特征向量总数n_vecs、向量总长度n_total：

import numpy as np

def read_slepc_parallel_eigvecs(filename, n_procs, n_vecs, n_total):
    # 读取文件原始数据
    raw_data = np.loadtxt(filename)
    # 按全零行分割各个特征向量的块
    vec_blocks = []
    current_block = []
    for row in raw_data:
        if np.allclose(row, 0, atol=1e-12):
            if current_block:
                vec_blocks.append(np.array(current_block))
                current_block = []
        else:
            current_block.append(row)
    if current_block:
        vec_blocks.append(np.array(current_block))
    
    # 计算每个进程负责的向量分量长度
    local_lengths = [
        n_total // n_procs + (1 if i < n_total % n_procs else 0)
        for i in range(n_procs)
    ]
    split_points = np.cumsum([0] + local_lengths)
    
    # 拼接每个向量的完整分量
    full_eigvecs = []
    for block in vec_blocks:
        proc_parts = [
            block[split_points[i]:split_points[i+1]]
            for i in range(n_procs)
        ]
        full_vec = np.concatenate(proc_parts).flatten()
        full_eigvecs.append(full_vec)
    
    # 返回矩阵，每列对应一个特征向量
    return np.array(full_eigvecs).T

Julia实现 #

同样需提前知晓n_procs、n_vecs、n_total：

using DelimitedFiles, LinearAlgebra

function read_slepc_parallel_eigvecs(filename, n_procs, n_vecs, n_total)
    raw_data = readdlm(filename)
    vec_blocks = []
    current_block = []
    for row in eachrow(raw_data)
        if all(isapprox.(row, 0; atol=1e-12))
            if !isempty(current_block)
                push!(vec_blocks, hcat(current_block...))
                current_block = []
            end
        else
            push!(current_block, row)
        end
    end
    if !isempty(current_block)
        push!(vec_blocks, hcat(current_block...))
    end

    # 计算各进程本地分量长度
    local_lengths = [
        n_total ÷ n_procs + (i <= n_total % n_procs ? 1 : 0)
        for i in 1:n_procs
    ]
    split_points = cumsum([0; local_lengths])

    full_eigvecs = []
    for block in vec_blocks
        proc_parts = [
            block[split_points[i]+1:split_points[i+1], :]
            for i in 1:n_procs
        ]
        full_vec = vcat(proc_parts...)
        push!(full_eigvecs, full_vec)
    end

    # 返回矩阵，每列对应一个特征向量
    return hcat(full_eigvecs...)
end

内容的提问来源于stack exchange，提问作者devanshu shekhar