我的系统是一个带有拉格朗日乘数的对称有限元问题：

$Z=\begin{pmatrix}A & C^T \\ C & 0\end{pmatrix}$

矩阵$A$是半正定的，不可逆的。整个矩阵是可逆的。

我正在开发一个有限元软件，我想用 PETSc 并行解决这些问题。这时候我想$Z$矩阵是MATNEST其中的$A$和$C$矩阵是独立组装的。

使用此消息末尾的代表性简单代码（具有拉格朗日乘数的一维拉普拉斯算子），并使用PCFIELDSPLIT预条件器，求解器发散。如果我人为地强制矩阵$A$是可逆的（添加命令行选项-force_invertible），那么我得到了正确的解决方案。

那么，我能做些什么来解决这样一个带有嵌套矩阵的系统呢？

另一个问题：是否可以并行使用带有嵌套矩阵的直接方法（Mumps）？

#include "petscsys.h" /* framework routines */
#include "petscvec.h" /* vectors */
#include "petscmat.h" /* matrices */
#include "petscksp.h"

#include <vector>
#include <string>
#include <iostream>
#include <numeric>

// Try to solve with Petsc a simple 1D laplacian problem (or a series of springs) 
// o-////-o-////-o-////-o-////-o-////-o-////-o-////-o-////-o ...
//
// The boundary conditions (Dirichlet) are imposed using Lagrange multipliers.
// The system to be solved is of the form:
//
//      Z     x  =  y
//    --^--  -^-   -^-
//    [A Ct] [u] = [b]
//    [C 0 ] [l]   [d]
//
// A is positive semi-definite (1 eigenvalue is 0)
// Z is indefinite invertible

static char help[] = "Saddle point problem: scalar 1D laplacian with lagrange multipliers.\n\n";
int main(int argc,char **argv)
{

  // Initialization ------------------------------------------------ //
  MPI_Init(NULL, NULL);  
  PetscErrorCode ierr;
  ierr = PetscInitialize(&argc,&argv,NULL,help);CHKERRQ(ierr);

  int rank, nproc;
  MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
  MPI_Comm_size(PETSC_COMM_WORLD, &nproc);

  PetscViewer viewer = PETSC_VIEWER_STDOUT_(PETSC_COMM_WORLD);
  PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_DENSE);

  PetscBool forceInvertible=PETSC_FALSE;
  PetscOptionsGetBool(NULL,NULL, "-force_invertible", &forceInvertible, NULL);

  // Problem definition and dof spliting among processes ----------- //
  double k=1.;

  std::vector<int> globalDofs={0, 1, 2, 3, 4, 5, 6, 7, 8};

  std::vector<std::pair<int, double>> globalDirichlet={{0, 0.},
                                                       {8, 1.}}; // first: global dof number, second: imposed value

  auto start=[&](int rank){
    int q=(globalDofs.size()-1)/nproc;
    int r=(globalDofs.size()-1)%nproc;
    return q*rank + std::min(rank, r);
  };

  std::vector<int> dofs;
  for (int i=start(rank); i<=start(rank+1); ++i) {
    dofs.push_back(globalDofs[i]);
  }

  std::vector<int> isLocal(dofs.size(), 1);
  if (rank!=0) isLocal[0]=0;
  int nLocalDofs=std::accumulate(isLocal.begin(), isLocal.end(), 0);  

  std::vector<int> dirichletIdx;
  for (int j=0; j<globalDirichlet.size(); ++j) {
    for (int i=0; i<dofs.size(); ++i) {
      if (isLocal[i] && globalDirichlet[j].first==dofs[i]) {
        dirichletIdx.push_back(j);
      }
    }
  }

  for (int p=0; p<nproc; ++p) {
    MPI_Barrier(PETSC_COMM_WORLD);
    if (rank==p) {
      std::cout << "Rank: " << rank << std::endl;

      std::cout << "  dofs   : ";
      for (int d: dofs) {
        std::cout << d << " ";
      }
      std::cout << std::endl;

      std::cout << "  isLocal: ";
      for (int d: isLocal) {
        std::cout << d << " ";
      }
      std::cout << std::endl;

      std::cout << "  nLocalDofs: " << nLocalDofs;
      std::cout << std::endl;

      std::cout << "  dirichlet: ";
      for (int i: dirichletIdx) {
        std::cout << "{"
                  << globalDirichlet[i].first << ", "
                  << globalDirichlet[i].second
                  << "} ";
      }
      std::cout << std::endl;

    }
  }

  // Matrix A ------------------------------------------------------ //
  Mat A;
  MatCreate(PETSC_COMM_WORLD, &A);
  MatSetType(A, MATMPIAIJ);
  MatSetSizes(A, nLocalDofs, nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE);

  MatMPIAIJSetPreallocation(A, 5, NULL, 1, NULL);

  auto setValue=[&](int i, int j, double v){
    if (forceInvertible && (i==2 || j==2)) {
      MatSetValue(A, 2, 2, 1., ADD_VALUES);
    }
    else {
      MatSetValue(A, i, j, v, ADD_VALUES);
    }
  };

  for (int i=0; i<dofs.size()-1; ++i) {
    // k * [ -1  1 ]
    //     [  1 -1 ]
    setValue(dofs[i]  , dofs[i]  , -k);
    setValue(dofs[i+1], dofs[i+1], -k);
    setValue(dofs[i]  , dofs[i+1],  k);
    setValue(dofs[i+1], dofs[i]  ,  k);

  }

  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

  MatView(A, viewer);


  // Matrix C ------------------------------------------------------ //  
  Mat C;
  MatCreate(PETSC_COMM_WORLD, &C);
  MatSetType(C, MATMPIAIJ);

  MatSetSizes(C, dirichletIdx.size(), nLocalDofs, PETSC_DETERMINE, PETSC_DETERMINE);
  MatMPIAIJSetPreallocation(C, 1, NULL, 0, NULL);
  for (int i: dirichletIdx) {
    MatSetValue(C, i, globalDirichlet[i].first, 1., ADD_VALUES);
  }

  MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);

  // Zero Matrix --------------------------------------------------- //
  Mat N;
  MatCreate(PETSC_COMM_WORLD, &N);
  MatSetType(N, MATMPIAIJ);
  MatSetSizes(N, dirichletIdx.size(), dirichletIdx.size(), PETSC_DETERMINE, PETSC_DETERMINE);

  MatMPIAIJSetPreallocation(N, 0, NULL, 0, NULL);
  MatAssemblyBegin(N, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(N, MAT_FINAL_ASSEMBLY);  

  // Matrix Z ------------------------------------------------------ //
  Mat Z;
  Mat subZ[4];

  subZ[0]=A;
  MatCreateTranspose(C, &subZ[1]);
  subZ[2]=C;
  subZ[3]=N;

  MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, subZ, &Z);
  MatView(Z, viewer);

  // Vector b ------------------------------------------------------ //
  Vec b;
  VecCreate(PETSC_COMM_WORLD, &b);
  VecSetType(b, VECMPI);
  VecSetSizes(b, nLocalDofs, PETSC_DECIDE);

  VecAssemblyBegin(b);
  VecAssemblyEnd(b);

  // Vector d ------------------------------------------------------ //
  Vec d;
  VecCreate(PETSC_COMM_WORLD, &d);
  VecSetType(d, VECMPI);
  VecSetSizes(d, dirichletIdx.size(), PETSC_DECIDE);

  for (int i: dirichletIdx) {
    VecSetValue(d, i, globalDirichlet[i].second, ADD_VALUES);
  }

  VecAssemblyBegin(d);
  VecAssemblyEnd(d);

  // Vector y ------------------------------------------------------ //
  Vec y;
  Vec suby[2];

  suby[0]=b;
  suby[1]=d;

  VecCreateNest(PETSC_COMM_WORLD, 2, NULL, suby, &y);

  VecView(y, viewer);

  // KSP ----------------------------------------------------------- //

  KSP ksp;
  KSPCreate(PETSC_COMM_WORLD,&ksp);
  KSPSetOperators(ksp, Z, Z);
  KSPSetFromOptions(ksp);

  PC pc;
  KSPGetPC(ksp,&pc);
  PCSetType(pc, PCFIELDSPLIT);
  PCFieldSplitSetDetectSaddlePoint(pc, PETSC_TRUE);
  IS isg[2];
  MatNestGetISs(Z, isg, NULL);
  PCFieldSplitSetIS(pc, "u", isg[0]);
  PCFieldSplitSetIS(pc, "l", isg[1]);

  // PC pc;
  // KSPSetType(ksp, KSPPREONLY);
  // KSPGetPC(ksp,&pc);
  // PCSetType(pc, PCLU);
  // PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);  

  Vec x;
  ierr = VecDuplicate(y, &x);
  KSPSolve(ksp, y, x);

  // View ---------------------------------------------------------- //
  MatView(Z, viewer);
  VecView(y, viewer);
  KSPView(ksp, viewer);
  VecView(x, viewer);
  // --------------------------------------------------------------- //

  ierr = PetscFinalize();

  MPI_Finalize();
  return ierr;
}