# Solving saddle point problem having non-invertible top-left block with a PETSc nested matrix

My system is a symmetric FE problem with lagrange multipliers:

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

The matrix $$A$$ is positive semi-definite, non-invertible. The whole matrix is invertible.

I am working on the development of a finite element software where I want to solve such problems with PETSc in parallel. At this time I would like the $$Z$$ matrix to be a MATNEST in which the $$A$$ and $$C$$ matrices are assembled independently.

With the representative simple code at the end of this message (1D laplacian with Lagrange multipliers), and using a PCFIELDSPLIT preconditioner, the solver diverges. If I artificially force the matrix $$A$$ to be invertible (adding the command-line option -force_invertible), then I get the right solution.

So, what can I do to solve such a system with a nested matrix?

Another question: is it possible to use direct methods (Mumps) with nested matrices in parallel?

#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)) {
}
else {
}
};

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) {
}

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) {
}

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);
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;
}

• As for direct methods, you definitely can create a monolithic matrix for Z by simply joining the sparse representations. The memory overhead even if you duplicate the storage should be negligible compared to the cost of the factorization. Jan 18, 2020 at 2:03