Reuse matrix and vector

[Jed Brown]

jarunan at ascomp.ch wrote:
> 
>> I suspect you are using a stale preconditioner, but
>> MatMPIAIJSetPreallocationCSR should not be called every iteration.
> 
> What do you mean 'stale preconditioner' ? I use Additive Schwarz.

That the preconditioner was not being updated when you change the
matrix.  If you reset everything inside the loop, then that isn't the
problem.

>> It's better to call this once at the beginning and update with
>> MatSetValues (insert one row at a time, it's very little code).
>> MatMPIAIJSetPreallocationCSR reallocates memory because it doesn't check
>> to see if the nonzero pattern has changed (because that's not what it's
>> for).
>>
> 
> Yes, thank you for the advice. I will modify this. But MatSetValues() is
> not efficient with a big problem. It takes much time.

No, either you are inserting values that have not been preallocated
(check with -info | grep mallocs) or you are inserting single values.
You should insert a full row every time you call MatSetValues.

>> Where are you calling KSPSetOperators?  You have to call this every time
>> the matrix changes.
> 
> I did shortcut of the code, actually...just to put it here. After
> creating and setting Matrix and vector. In each loop I create the KSP:
> 
>     call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
>     call KSPSetOperators(ksp,Ap,Ap,DIFFERENT_NONZERO_PATTERN,ierr)
>     call KSPGetPC(ksp,pc,ierr)
>     call PCSetType(pc,pct,ierr)
> 
>     if (pct == PCASM) then
>        call PCASMSetTotalSubdomains(pc,glob_nblocks_psc,PETSC_NULL_OBJECT,ierr)
>        call PCASMSetLocalSubdomains(pc,nblocks_psc,PETSC_NULL_OBJECT,ierr)

Only call one of these.

>        call PetscOptionsSetValue('-pc_asm_overlap','1',ierr)
>        call PetscOptionsSetValue('-sub_pc_type','lu',ierr)
>        call PetscOptionsSetValue('-sub_pc_factor_zeropivot','0.0',ierr)
>     endif
> 
>     call KSPSetTolerances(ksp,resin,1.e-20,  &
>              PETSC_DEFAULT_DOUBLE_PRECISION,nswp_psc,ierr)
> 
>     call KSPSetType(ksp,kspt,ierr)
>     call KSPSetFromOptions(ksp,ierr)

Of the above, only KSPSetOperators() should be called inside the loop,
everything else is setup that should happen before your loop.

>     call KSPSolve(ksp,rhs,sol,ierr)
>
>     call KSPDestroy(ksp,ierr)


How was your code, and the convergence different before?

Jed

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 261 bytes
Desc: OpenPGP digital signature
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20091105/0ea1ef14/attachment.pgp>