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

[Jared Crean]

     Yes, I meant MATAIJ. Sorry for the typo.

     Jared Crean

On 02/05/2018 12:46 PM, Stefano Zampini wrote:
> You mean MATAIJ?
>
> Note that you can run the same code in sequential or in parallel with MATAIJ, i.e.
>
> MatCreate(comm,&A);
> MatSetType(A,MATAIJ);
> MatSeqAIJSetPreallocation(A,dnz,ddnz); // dummy call if comm.size > 1
> MatMPIAIJSetPreallocation(A,dnz,ddnz,onz,oonz); // dummy call if comm.size == 1
>
> or by using MatXAIJSetPreallocation
>
>> On Feb 5, 2018, at 8:43 PM, Jared Crean <jcrean01 at gmail.com> wrote:
>>
>>      Switching to MPIAIJ fixed the problem.
>>
>>      Thanks for looking into this,
>>          Jared Crean
>>
>> On 02/02/2018 01:56 PM, Smith, Barry F. wrote:
>>>> On Feb 2, 2018, at 12:27 PM, Stefano Zampini <stefano.zampini at gmail.com> wrote:
>>>>
>>>> Jared,
>>>>
>>>> the problem arises using MPIAIJ for the uniprocessor case and it is related with BJacobi, nothing to do with your code.
>>>> If you set the type to MATAIJ and call MatSeqAIJSetPreallocation, it works fine.
>>>>
>>>> Barry, should we change the default pc type for MPIAIJ if the size is 1? MatGetFactor and MatSolverTypeGet cannot take this decision.
>>>    I think we need to improve the documentation and examples to emphasis the use of MATAIJ and not the particular Seq or MPI version so that people always think in terms of MATAIJ and so naturally write code that works for all cases.
>>>
>>>     Barry
>>>
>>> We tried once a very long time ago to have MPIAIJ become SeqAIJ on a single process but it ended up causing its own headaches.
>>>
>>>
>>>> Stefano
>>>>
>>>>> On Feb 1, 2018, at 8:48 PM, Jared Crean <jcrean01 at gmail.com> wrote:
>>>>>
>>>>> 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
>>>>>
>>>>>
>>>>> <testcase.c>
>>