Questions about petsc matrix and KSP

[Barry Smith]

On Fri, 3 Feb 2006, Gang Cheng wrote:

> Hi, guys,
>
> I have some questions about the usage of petsc matrix and KSP. The
> parallel Fortran code I'm working on simulates a biological process in
> which a diffusion-reaction PDE needs to be solved at each time step in a
> loop for many steps. Moreover, the coefficient matrix (H) for the
> discretized (using finite difference) PDE changes as the system evolves.
> So it is necessary to update H at each time step. My current algorithm is:
>
> Before the loop
> ---------------
> 1A. Create the parallel matrix H using MatCreateMPIAIJ and save it in a
> common block.
> 2A. According to the initial condition of the domain, fill non-zero values
> in H row by row using MatSetValues. H is then assembed by calling
> MatAssemblyBegin and MatAssemblyEnd.
>
> At each time step within the loop
> ---------------------------------
> 1B. Each processor scans local domain and ONLY updates (using
> MatSetValues) the rows in H that need to be updated (e.g., a certain grid
> has changed and a new boundary condition needs to be applied to its
> interface with other grids). MatAssemblyBegin and MatAssemblyEnd are
> called even if no update is done in the local processor to avoid deadlock.
> 2B. Create a KSP from H, calculate right-hand-side vector and solve the
> linear system using KSPSolve.
> 3B. Destroy the KSP.

You should not create a new KSP each time, reuse the old one and just
call KSPSetOperators() again with the matrix.

>
> Here is the problem I'm having with the above algorithm. In step 1B, I
> didn't update the entire H with the intention to save some CPU time (I
> found that the MatSetValues operiaiton is very slow for parallel
> matrices). However, I found that the updated values were not inserted in
> the supposed postion in H, causing errors in the solution of the PDE (I'm
> certain that the input parameters for MatSetValues were right). If I
> instead update the entire H while keeping everything else the same, the
> PDE solution is correct. So my specific question is, if a matrix needs to
> be updated, do we need to update ALL its rows before assembling it? Can we
> update just some of them?

  You are absolutely suppose to be able to just change some of the values
(those that are not changed should stay where they are), but if you call
MatZeroEntries() then obviously you need to reset everything.

  I would do a MatView() on the matrix (the second time through the loop)
for both approaches (1) change only some values and (2) set all the values
and try to see what is different about the two matrices. In theory they
should be the same.

> Also, I'm sure some of you guys are using petsc
> for problems similar to mine. If you have any suggestions on how I may
> improve the efficiency of my algorithm, I'd greatly appreciate them.

After you have it working well you can try to lag the preconditioner.
You can call KSPSetOperators() with the flag SAME_PRECONDITIONER
for several iterations, this may speed things up since it will not
have to spend the time rebuilding the preconditioner each step.
As the matrix changes likely the preconditioner will work less and less
well and so you have to call KSPSetOperators() without that flag
to get it to build a new preconditioner. Of course, this approach
need not always work; if the matrix spectrum is changing rapidly
you will need to rebuild the preconditioner each time. But it is
worth trying.

   Barry

>
> Sorry for this long email.
>
> -GC
>
>
>
>
>
>
>