Direct LU solver

[Hong Zhang]

Amit,

> Thanks for your response. A couple of follow-up questions -
>
> I go over the steps
>
> MatCreate()
> MatSetType(MPIAIJ)
> MatMPIAIJSetPreallocation()
> MatSetType(MATSUPERLU_DIST).
>
> Even though I want to do the direct LU solves repeatedly (with the same
> matrix), I don't want the program to do LU factorization repeatedly. I hope
> I can get that functionality by using the "SAME_PRECONDITIONER" option
> (along with -ksppreonly), while defining the KSP. When the program does the
> factorization, does it do it in-place or does it allocate new storage ?

allocate new storage unless you call
MatLUFactor() explicitly.

It is done by
  KSPSetUp()
    which calls PCSetUp() -> PCSetUp_LU()
      -> MatLUFactorSymbolic() & MatLUFactorNumeric(

You can do following for reusing factored L and U:
   KSPCreate()
   KSPSetOperators()
   KSPSetFromOptions()
   KSPSetUp()
   while ( num_rhs-- ) {
     KSPSolve(ksp,b,x);
   }

See src/ksp/ksp/examples/tutorials/ex10.c

> After doing the LU factorization, are the other operations such as
> MatMult() preserved ?

Yes.

Hong

>
> Thanks
>
> Rgds,
> Amit
>
>
>
>
>
>             Satish Balay
>             <balay at mcs.anl.go
>             v>                                                         To
>             Sent by:                  petsc-users at mcs.anl.gov
>             owner-petsc-users                                          cc
>             @mcs.anl.gov
>             No Phone Info                                         Subject
>             Available                 Re: Direct LU solver
>
>
>             02/29/2008 02:06
>             PM
>
>
>             Please respond to
>             petsc-users at mcs.a
>                  nl.gov
>
>
>
>
>
>
> On Fri, 29 Feb 2008, Amit.Itagi at seagate.com wrote:
>
>>
>> My woes continue.  Based on the earlier discussions, I implemented the
>> matrix as
>>
>>
> //=========================================================================
>>
>>  //   Option 1
>>   ierr=MatCreate(PETSC_COMM_WORLD,&A); CHKERRQ(ierr);
>>   ierr=MatSetSizes(A,1,1,2,2); CHKERRQ(ierr);
>>
>>
>>   /*    Option 2
>>   PetscInt d_nnz=1, o_nnz=1;
>>   ierr=MatCreateMPIAIJ(PETSC_COMM_WORLD,1,1,2,2,0,&d_nnz,0,&o_nnz,&A);
>> CHKERRQ(ierr);
>>   */
>>
>>   /*   Option 3
>>
>>
> ierr=MatCreateMPIAIJ(PETSC_COMM_WORLD,1,1,2,2,0,PETSC_NULL,0,PETSC_NULL,&A);
>
>>  CHKERRQ(ierr);
>>   */
>>
>>   ierr=MatSetType(A,MATSUPERLU_DIST); CHKERRQ(ierr);
>>   ierr=MatSetFromOptions(A); CHKERRQ(ierr);
>>
>>   // (After this, I set the values and do the assembly). I then use the
>> direct LU solver.
>>
>>
> //============================================================================
>
>>
>> Note: I have a simple 2 by 2 matrix (with non-zero values in all 4
> places).
>> If I use "option 1" (based on Satish's email), the program executes
>> successfully. If instead of "option 1", I use "option 2" or "option 3", I
>> get a crash.
>> If I am not mistaken, options 1 and 3 are the same. Option 2,
> additionally,
>> does a pre-allocation. Am I correct ?
>
> Nope - Option 3 is same as:
>
> MatCreate()
> MatSetType(MPIAIJ)
> MatMPIAIJSetPreallocation()
> MatSetType(MATSUPERLU_DIST)
>
> [i.e first you are setting type as MPIAIJ, and then changing to
> MATSUPERLU_DIST]
>
> What you want is:
>
> MatCreate()
> MatSetType(MATSUPERLU_DIST)
> MatMPIAIJSetPreallocation()
>
> [Ideally you need MatSuerLU_DistSetPreallocation() - but that would be
> same as MatMPIAIJSetPreallocation()]
>
> Satish
>
>
>
>