PETSc - 在本地操作 BAIJ 矩阵

计算科学 C++ 宠物
2021-12-26 07:44:32

我的程序加载了一个并行的 PETSc 矩阵A在几个 MPI 进程上，每个进程都有一个块子矩阵Ai.

我想检索本地子矩阵Ai, 对应于当前排名的那个。

我不想改变矩阵Ai，我只需要阅读它并计算一些信息（行列式、特征值等）。

int main(int argc, char **args) {
    PetscErrorCode  ierr;
    PetscViewer     fd;
    Mat A;

    ierr = SlepcInitialize(&argc, &args, (char *) nullptr, nullptr);
    if (ierr) return ierr;

    // Open file
    PetscViewerBinaryOpen(PETSC_COMM_WORLD, "matrix.mat", FILE_MODE_READ, &fd);

    // Load into Petsc object
    MatCreate(PETSC_COMM_WORLD, &A);
    MatSetType(A, MATMPIBAIJ);
    MatLoad(A, fd);

    PetscViewerDestroy(&fd);

    // Here I would like to manipulate the local submatrix of A
    // For example to show information with MatView


    // Clear workspace
    MatDestroy(&A);
    SlepcFinalize();
    return 0;
}

检索当前等级的本地子矩阵的正确方法是什么？

注意：我SlepcInitialize使用 SLEPc 是为了在之后使用 SLEPc 作为特征值。

编辑 ：

我有指向的指针VecGetArray，我认为我需要矩阵的等价物，但MatGetArray不应该根据文档（和另一篇文章）使用

1个回答

我为 PETSc 3.7 找到的解决方案是使用MatGetSubMatrices。

警告：此方法在最后一个 API 中不存在。

我发现这BAIJ实际上不是用于块矩阵，而是在将矩阵分布到处理器上时按组保留行。

因此我使用AIJ了格式。（MPIAIJ得到一个平行矩阵）。

首先有权限行和列索引（这与 PETSc 无关）。

我的方法split_indices_rows并split_indices_rows拆分索引，以便每个 MPI 等级知道通过调用indices_rows_rank[mpi_rank]和处理矩阵的哪一部分indices_cols_rank[mpi_rank]。

// indices_rows contain all rows indices of original matrix : {O , .. , height}
vector<int> indices_rows;
for (int i = 0; i < matrix_height; i++)
    indices_rows.push_back(i);

// indices_columns contain all columns indices of original matrix : {0 , .. , width}
vector<int> indices_columns;
for (int i = 0; i < matrix_width; i++)
    indices_columns.push_back(i);

vector<vector<int> > indices_rows_rank = split_indices_rows(indices_rows);
vector<vector<int> > indices_cols_rank = split_indices_cols(indices_cols);

然后您可以创建索引集 (IS)。

int mpi_rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &mpi_rank);

// Create index sets
vector<int> local_rows = indices_rows_rank[mpi_rank];
IS local_row_indexes;
ISCreateGeneral(PETSC_COMM_WORLD, local_rows.size(), &local_rows[0], PETSC_COPY_VALUES, &local_row_indexes);

vector<int> local_cols = indices_cols_rank[mpi_rank];
IS local_col_indexes;
ISCreateGeneral(PETSC_COMM_WORLD, local_cols.size(), &local_cols[0], PETSC_COPY_VALUES, &local_col_indexes);

最后打电话MatGetSubMatrices

Mat *local_submatrices;
MatGetSubMatrices(A, 1, &local_row_indexes, &local_col_indexes, MAT_INITIAL_MATRIX, &local_submatrices);
其它你可能感兴趣的问题
上一篇具有涉及矩阵逆的约束的凸优化 下一篇是否可以使用 gmsh 重新网格化？