[petsc-dev] factortype and -pc_type lu

[Pierre Jolivet]

> On Aug 17, 2016, at 5:21 PM, Matthew Knepley <knepley at gmail.com> wrote:
> 
> On Wed, Aug 17, 2016 at 10:15 AM, Pierre Jolivet <pierre.jolivet at enseeiht.fr <mailto:pierre.jolivet at enseeiht.fr>> wrote:
> Dear Hong,
> Thank you for your answer.
> I understand the error message, but I don’t understand why setting the preconditioner to be of type LU, i.e., factorizing the matrix, does not set the factortype to MAT_FACTOR_LU.
> 
> It sounds like you were using the wrong matrix. If the factorization is out-of-place, then A is not the factored matrix.

Ah, I didn’t think about in-place/out-of-place factorization. That makes perfect sense, thank you!

>    Matt
>  
> More puzzling to me is the fact that in the output, one can see:
> 
>       Factored matrix follows:
>         Mat Object: 1 MPI processes
> 
> but that the error is then that the matrix is unfactored. Anyway, as you suggested, this works:
>   Mat F;
>   PCFactorGetMatrix(pc, &F);
>   MatSolve(F, b, x);
> 
> Thank you,
> Pierre
> 
>> On Aug 17, 2016, at 4:34 PM, Hong <hzhang at mcs.anl.gov <mailto:hzhang at mcs.anl.gov>> wrote:
>> 
>> Pierre:
>> 
>> [0]PETSC ERROR: Object is in wrong state
>> [0]PETSC ERROR: Unfactored matrix
>> 
>> indicates that your input matrix is an  Unfactored matrix.
>> 
>> If you want to replace KSPSolve() with your own MatMatSolve, you need use following blocks
>> to replace KSPSolve() :
>> Mat F;
>>   ierr = KSPGetPC(ksp,&pc);
>>   ierr = PCFactorSetMatSolverPackage(pc,MATSOLVERPETSC);
>>   ierr = PCFactorSetUpMatSolverPackage(pc); /* call MatGetFactor() to create F */
>>   ierr = PCFactorGetMatrix(pc,&F);
>>   ierr = KSPSetUp(ksp);
>>   ierr = MatSolve(F,b,x);
>> 
>> Hong
>> 
>>  
>> [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html <http://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting.
>> [0]PETSC ERROR: Petsc Development GIT revision: v3.7.3-1165-gfeaa1dd  GIT Date: 2016-08-16 11:58:28 -0500
>> 
>> Is this the expected behavior? How should I call MatSolve?
>> The rationale for calling MatSolve instead of KSPSolve (for which there is obviously no error) is that I need to use MatMatSolve in my application—because there is no KSPSolve(KSP, Mat, Mat), and the behavior is the same for both MatSolve and MatMatSolve.
>> 
>> Thank you,
>> Pierre
>> 
> 
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20160818/c62e1805/attachment.html>