[petsc-users] Error Using Matrix-Free Linear Operator with Matrix-Explicit Preconditioner

[Jared Crean]

Hello,

    I am trying to use a matrix-free linear operator with a 
matrix-explicit preconditioner, but when I try to do the KSP solve it 
gives the error:

[0]PETSC ERROR: See 
http://www.mcs.anl.gov/petsc/documentation/linearsolvertable.html for 
possible LU and Cholesky solvers
[0]PETSC ERROR: Could not locate a solver package. Perhaps you must 
./configure with --download-<package>
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.8.3, Dec, 09, 2017
[0]PETSC ERROR: ./testcase on a arch-linux2-c-debug named jared-r15 by 
jared Thu Feb  1 12:40:57 2018
[0]PETSC ERROR: Configure options
[0]PETSC ERROR: #1 MatGetFactor() line 4346 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/mat/interface/matrix.c
[0]PETSC ERROR: #2 PCSetUp_ILU() line 142 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/pc/impls/factor/ilu/ilu.c
[0]PETSC ERROR: #3 PCSetUp() line 924 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #4 KSPSetUp() line 381 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #5 PCSetUpOnBlocks_BJacobi_Singleblock() line 618 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/pc/impls/bjacobi/bjacobi.c
[0]PETSC ERROR: #6 PCSetUpOnBlocks() line 955 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/pc/interface/precon.c
[0]PETSC ERROR: #7 KSPSetUpOnBlocks() line 213 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #8 KSPSolve() line 613 in 
/home/jared/.julia/v0.4/PETSc2/deps/petsc-3.8.3/src/ksp/ksp/interface/itfunc.c


  The code to reproduce this error is attached.  The error is present 
on Petsc 3.7.6 and 3.8.3.  I noticed two things while creating the test 
case: 1) using a jacobi preconditioner works (using block jacobi with 
ILU on each block does not), and 2) if I replace the shell matrix with 
the preconditioner matrix in KSPSetOperators(), there is no error (with 
the block jacobi ILU preconditioner).

  Is this a bug in Petsc or did I setup the preconditioner incorrectly?

    Jared Crean


static char help[] = "Test using a matrix-free linear operator with a matrix-explicit preconditioner";

#include <petscksp.h>

// use a global constant for the matrix values
const PetscInt N = 4;  // matrix is N x N
const PetscScalar Avals[] = { 1 + 2, 1,     1,      1,
                              1,     1 + 3, 1,      1,
                              1,     1,     1 + 4,  1,
                              1,     1,     1,      1 + 5};

// multiplication function for shell matrix
PetscErrorCode applyLinearOperator(Mat A, Vec x, Vec b);

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{

  PetscInitialize(&argc,&args,(char*)0,help);
  const char* option_keys[] = {"-ksp_monitor",
                               "-pc_type",
                               "-sub_pc_factor_levels",
                               "-sub_pc_type"
                        };
  const char* option_vals[] = {"",
                               "bjacobi",
                               "4",
                               "ilu"
                              };

  int noptions = sizeof(option_keys)/sizeof(option_keys[0]);
  for (int i=0; i < noptions; ++i)
    PetscOptionsSetValue(NULL, option_keys[i], option_vals[i]);


  // create preconditioner matrix
  Mat Ap;
  MatCreate(MPI_COMM_WORLD, &Ap);
  MatSetFromOptions(Ap);
  MatSetType(Ap, MATMPIAIJ);
  MatSetSizes(Ap, N, N, PETSC_DECIDE, PETSC_DECIDE);
  MatMPIAIJSetPreallocation(Ap, N, NULL, 0, NULL);

  // create the shell matrix
  Mat A;
  MatCreateShell(MPI_COMM_WORLD, N, N, PETSC_DECIDE, PETSC_DECIDE, NULL, &A);
  MatSetFromOptions(A);
  MatSetUp(A);
  MatShellSetOperation(A, MATOP_MULT, (void*)(&applyLinearOperator));

  // create the KSP
  KSP ksp;
  KSPCreate(MPI_COMM_WORLD, &ksp);
  KSPSetFromOptions(ksp);
  KSPSetOperators(ksp, A, Ap);

  // create vectors
  Vec x, b;
  VecCreate(MPI_COMM_WORLD, &b);
  VecSetType(b, VECMPI);
  VecSetSizes(b, N, PETSC_DECIDE);

  VecCreate(MPI_COMM_WORLD, &x);
  VecSetType(x, VECMPI);
  VecSetSizes(x, N, PETSC_DECIDE);


  // put values into Ap and b
  PetscInt idx[N];
  PetscScalar bvals[N];

  for (int i=0; i < N; ++i)
  {
    idx[i] = i;
    bvals[i] = 1;
  }

  MatSetValues(Ap, N, idx, N, idx, Avals, ADD_VALUES);
  VecSetValues(b, N, idx, bvals, ADD_VALUES);


  // assemble
  MatAssemblyBegin(Ap, MAT_FINAL_ASSEMBLY);
  VecAssemblyBegin(b);

  MatAssemblyEnd(Ap, MAT_FINAL_ASSEMBLY);
  VecAssemblyEnd(b);

  // solve - this triggers the error
  KSPSolve(ksp, b, x);
  VecView(x, PETSC_VIEWER_STDOUT_WORLD);


  // cleanup
  VecDestroy(&b);
  VecDestroy(&x);
  MatDestroy(&A);
  KSPDestroy(&ksp);

  PetscFinalize();
  return 0;
}


// multiplication function for shell matrix
PetscErrorCode applyLinearOperator(Mat A, Vec x, Vec b)
{

  PetscScalar* btmp;
  const PetscScalar* xtmp;

  VecGetArray(b, &btmp);
  VecGetArrayRead(x, &xtmp);

  // zero out b just in case
  for (int i=0; i < N; ++i)
    btmp[i] = 0;

  // multiply
  for (int i=0; i < N; ++i)
    for (int j=0; j < N; ++j)
      btmp[i] += Avals[j + i*N]*xtmp[j];

  VecRestoreArray(b, &btmp);
  VecRestoreArrayRead(x, &xtmp);

  return 0;
}