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