[petsc-users] KSPSolve with MPIAIJ with non-square 'diagonal parts'

[Olivier Jamond]

Hello,

I am sorry because I surely miss something, but I cannot manage to solve 
a problem with a MPIAIJ matrix which has non-square 'diagonal parts'.

I copy/paste at the bottom of this message a very simple piece of code 
which causes me troubles. I this code I try to do 'x=KSP(A)*b' (with 
gmres/jacobi), but this fails whereas a matrix multiplication 'b=A*x' 
seems to work. Is ksp with such a matrix supposed to work (I can't find 
anything in the documentation about that, so I guess that it is...) ?

Many thanks,
Olivier

NB: this code should be launched with exactly 3 procs

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

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

static char help[] = "Trying to solve a linear system on a sub-block of 
a matrix\n\n";
int main(int argc, char **argv)
{
   MPI_Init(NULL, NULL);
   PetscErrorCode ierr;
   ierr = PetscInitialize(&argc, &argv, NULL, help);
   CHKERRQ(ierr);

   // clang-format off
   std::vector<std::vector<double>> AA = {
       { 1,  2,  0, /**/ 0,  3, /**/ 0,  0,  4},
       { 0,  5,  6, /**/ 7,  0, /**/ 0,  8,  0},
       { 9,  0, 10, /**/11,  0, /**/ 0, 12,  0},
       //---------------------------------------
       {13,  0, 14, /**/15, 16, /**/17,  0,  0},
       { 0, 18,  0, /**/19, 20, /**/21,  0,  0},
       { 0,  0,  0, /**/22, 23, /**/ 1, 24,  0},
       //--------------------------------------
       {25, 26, 27, /**/ 0,  0, /**/28, 29,  0},
       {30,  0,  0, /**/31, 32, /**/33,  0, 34},
   };


   std::vector<double> bb = {1.,
                             1.,
                             1.,
                             //
                             1.,
                             1.,
                             1.,
                             //
                             1.,
                             1.};


   std::vector<int> nDofsRow = {3, 3, 2};
   std::vector<int> nDofsCol = {3, 2, 3};
   // clang-format on

   int NDofs = std::accumulate(nDofsRow.begin(), nDofsRow.end(), 0);

   int pRank, nProc;
   MPI_Comm_rank(PETSC_COMM_WORLD, &pRank);
   MPI_Comm_size(PETSC_COMM_WORLD, &nProc);

   if (nProc != 3) {
     std::cerr << "THIS TEST MUST BE LAUNCHED WITH EXACTLY 3 PROCS\n";
     abort();
   }

   Mat A;
   MatCreate(PETSC_COMM_WORLD, &A);
   MatSetType(A, MATMPIAIJ);
   MatSetSizes(A, nDofsRow[pRank], nDofsCol[pRank], PETSC_DETERMINE, 
PETSC_DETERMINE);
   MatMPIAIJSetPreallocation(A, NDofs, NULL, NDofs, NULL);

   Vec b;
   VecCreate(PETSC_COMM_WORLD, &b);
   VecSetType(b, VECMPI);
   VecSetSizes(b, nDofsRow[pRank], PETSC_DECIDE);

   if (pRank == 0) {
     for (int i = 0; i < NDofs; ++i) {
       for (int j = 0; j < NDofs; ++j) {
         if (AA[i][j] != 0.) {
           MatSetValue(A, i, j, AA[i][j], ADD_VALUES);
         }
       }
       VecSetValue(b, i, bb[i], ADD_VALUES);
     }
   }

   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
   VecAssemblyBegin(b);
   VecAssemblyEnd(b);

   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, 
PETSC_VIEWER_ASCII_DENSE);
   MatView(A, PETSC_VIEWER_STDOUT_WORLD);
   VecView(b, PETSC_VIEWER_STDOUT_WORLD);

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

   PC pc;
   KSPGetPC(ksp, &pc);
   PCSetFromOptions(pc);

   Vec x;
   MatCreateVecs(A, &x, NULL);
   ierr = KSPSolve(ksp, b, x);    // this fails
   MatMult(A, x, b);              // whereas the seems to be ok...

   VecView(x, PETSC_VIEWER_STDOUT_WORLD);

   MPI_Finalize();

   return 0;
}