PETSc矩阵求逆

PETSc简介

PETSc是Portable, Extensible Toolkits for Scientific Computing的缩写，由美国能源部(DOE)下属的Argonne国家实验室开发。PETSc基于MPI、BLAS、LAPACK等底层库，用于并行求解偏微分方程形成的大型线性方程组(SLES)，也包含非线性求解器(SNES)。自3.5版本开始，优化包Toolkit for Advanced Optimization (TAO) 也包含在PETSc中，现在一般称为PETSc/TAO。PETSc官网是: https://petsc.org。

在科学计算(高性能计算)领域，PETSc应用非常广泛。本人近期工作也基于PETSc展开，正在增加相关知识储备。

PETSc矩阵求逆

PETSc矩阵求逆的需求较为小众，官方不推荐直接进行求逆操作，而应该使用KSP相关函数(KSPCreate、KSPSolve等)求解线性系统。但有时候会遇到对矩阵的某个小矩阵块求逆的特殊需求，经过阅读官方文档和测试，本文给出PETSc矩阵求逆的源码供参考：

// 创建稠密单位矩阵
void MatCreateDenseDiag(int n, PetscScalar value, Mat* eye) {
  MatCreateDense(PETSC_COMM_SELF, n, n, n, n, NULL, eye);
  std::vector<PetscInt> rows(n);
  std::iota(rows.begin(), rows.end(), 0);
  MatZeroRows(*eye, n, rows.data(), value, NULL, NULL);
  MatAssemblyBegin(*eye, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(*eye, MAT_FINAL_ASSEMBLY);
}

// 求矩阵A的逆，输出稠密逆矩阵invA
void InverseMatrix(Mat A, Mat* invA) {
  PetscInt m, n;
  MatGetSize(A, &m, &n);
  assert(m == n);

  // 创建单元矩阵
  Mat eye;
  MatCreateDenseDiag(n, 1, &eye);

  // 一般来说，逆矩阵也都是稠密的
  MatDuplicate(eye, Mat_DO_NOT_COPY_VALUES, invA);

  // LU分解并求A的逆
  Mat factor;
  MatFactorInfo info;
  IS isrow, iscol;
  MatGetOrdering(A, MATORDERINGNATURAL, &isrow, &iscol);
  MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &factor);
  MatLUFactorSymbolic(factor, A, isrow, iscol, &info);
  MATLUFactorNumeric(factor, A, &info);
  MatMatSolve(factor, eye, *invA);

  // 销毁变量
  MatDestroy(&factor);
  MatDestroy(&eye);
  ISDestroy(&isrow);
  ISDestroy(&iscol);
}

一个简单的算例验证代码的正确性：

constexpr int n = 10;
PetscScalar value = 10;

Mat A, invA;
MatCreateConstantDiagonal(PETSC_COMM_SELF, n, n, n, n, value, &A);

InverseMatrix(A, &invA);

MatView(A, PETSC_VIEWER_STDOUT_SELF);
MatView(invA, PETSC_VIEWER_STDOUT_SELF);

MatDestroy(&A);
MatDestory(&invA);

参考

Matrix Factorization