PETSc中mpiaij矩阵调用MatLUFactor报错 配置SuperLU仍无法求逆求助

问题根因 #

• 原生MatLUFactor接口本身不支持分布式mpiaij类型矩阵，即使安装了SuperLU，也需要通过PETSc的求解器框架调用分布式SuperLU（SuperLU_DIST），而不是直接调用该接口
• 你当前创建单位矩阵Ones和逆矩阵inv_matrix时用了PETSC_COMM_SELF的顺序稠密矩阵，和分布式mpiaij矩阵的通信域、存储结构不匹配，也是报错诱因

解决步骤 #

1. 先确认PETSc编译时已开启SuperLU_DIST支持，在终端执行命令查看：
   petscinfo | grep superlu_dist
   如果无输出，重新编译PETSc时添加--download-superlu_dist参数开启分布式SuperLU支持。
2. 放弃直接调用MatLUFactor，改用KSP线性求解器框架配置求解参数，通过求解AX=I（I为单位矩阵）的方式得到逆矩阵，这是PETSc官方推荐的分布式矩阵求逆实现方式。
3. 单位矩阵和逆矩阵要使用和原输入矩阵同通信域的并行稠密矩阵，不要用顺序矩阵。

修正后实现代码 #

void inv_PETSC_matrix(Mat& matrix,  Mat &inv_matrix)
{
  Mat             Ones;
  PetscInt        petsc_m=0, petsc_n=0,i, rstart, rend;
  PetscErrorCode  ierr=0;
  PetscScalar     v = 1.0;
  KSP             ksp;
  PC              pc;

  ierr = MatGetSize(matrix, &petsc_m, &petsc_n);CHKERRVOID(ierr);
  // 检查矩阵是否为方阵，只有方阵可求逆
  if(petsc_m != petsc_n) SETERRQ(PetscObjectComm((PetscObject)matrix), PETSC_ERR_ARG_WRONG, "Only square matrix can be inverted");

  // 创建同通信域的并行稠密单位矩阵和逆矩阵
  ierr = MatCreateDense(PetscObjectComm((PetscObject)matrix), PETSC_DECIDE, PETSC_DECIDE, petsc_m, petsc_n, NULL, &Ones);CHKERRVOID(ierr);
  ierr = MatCreateDense(PetscObjectComm((PetscObject)matrix), PETSC_DECIDE, PETSC_DECIDE, petsc_m, petsc_n, NULL, &inv_matrix);CHKERRVOID(ierr);

  // 并行填充单位矩阵，只填充当前进程负责的行
  ierr = MatGetOwnershipRange(Ones, &rstart, &rend);CHKERRVOID(ierr);
  for(i = rstart; i < rend; i++){
    v = 1.0;
    ierr = MatSetValues(Ones, 1, &i, 1, &i, &v, INSERT_VALUES);CHKERRVOID(ierr);
  }
  ierr = MatAssemblyBegin(Ones, MAT_FINAL_ASSEMBLY);CHKERRVOID(ierr);
  ierr = MatAssemblyEnd(Ones, MAT_FINAL_ASSEMBLY);CHKERRVOID(ierr);

  // 配置KSP求解器，使用LU分解，求解器指定为SuperLU_DIST
  ierr = KSPCreate(PetscObjectComm((PetscObject)matrix), &ksp);CHKERRVOID(ierr);
  ierr = KSPSetOperators(ksp, matrix, matrix);CHKERRVOID(ierr);
  ierr = KSPGetPC(ksp, &pc);CHKERRVOID(ierr);
  ierr = PCSetType(pc, PCLU);CHKERRVOID(ierr);
  ierr = PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU_DIST);CHKERRVOID(ierr);
  ierr = KSPSetTolerances(ksp, 1e-15, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);CHKERRVOID(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRVOID(ierr); // 允许运行时参数覆盖配置

  // 求解AX=I，结果就是逆矩阵
  ierr = KSPSolveMat(ksp, Ones, inv_matrix);CHKERRVOID(ierr);

  // 释放资源
  ierr = KSPDestroy(&ksp);CHKERRVOID(ierr);
  ierr = MatDestroy(&Ones);CHKERRVOID(ierr);
}

运行注意事项 #

• 运行程序时可添加命令行参数-ksp_view查看求解器配置，确认是否正确加载了SuperLU_DIST
• 小规模矩阵求逆可以直接转成顺序矩阵计算，分布式求逆更适合大尺度稀疏矩阵场景

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