Questions about petsc matrix and KSP

[Gang Cheng]

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.

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? 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.

Sorry for this long email.

-GC