[petsc-users] solving multiple linear systems with same matrix (sequentially, not simultaneously)

[Daniel Goldberg]

Hello

In the code I am developing, there is a point where I will need to solve a
linear system with the same matrix but different right hand sides. It is a
time-dependent model and the matrix in question will change each timestep,
but for a single timestep ~30-40 linear solves will be required, with the
same matrix but differing right hand sides. Each right hand side depends on
the previous, so they cannot all be solved simultaneously.

Generally I solve the matrix via CG with  Block Jacobi / ILU(0)
preconditioning. I don't persist the matrix in between solves (i.e. I
destroy the Mat, KSP and Vec objects after each linear solve) because I do
not get the feeling this makes a large difference to performance and it is
easier in dealing with the rest of the model code, which does not use
PETSc. The matrix is positive definite with sparsity of 18 nonzero elements
per row, and generally the linear system is smaller than 1 million degrees
of freedom.

If my matrix was dense, then likely retaining LU factors for
foward/backward substitution would be a good strategy (though I'm not sure
if this is easily doable using PETSc direct solvers). But given the
sparsity I am unsure of whether I can take advantage of using the same
matrix 30-40 times in a row.

The comments in ex16.c state that when a KSP object is used multiple times
for a linear solve, the preconditioner operations are not done each time --
and so I figured that if I changed my preconditioning parameters I might be
able to make subsequent solves with the same KSP object faster. I tried an
experment in which I increased pc_factor_levels, and then created a number
of random vectors, e.g. vec_1, vec_2, vec_3... and called

KSPSolve(ksp, vec_i, solution_i)

sequentially with the same ksp object. I did see a decrease in the number
of CG iterations required as pc_factor_levels was increased, as well as an
increase in time for the first linear solve, but no decrease in time
required for subsequent solves. *Should* there be a timing decrease? But
more generally is there a way to optimize the fact I'm solving 40 linear
systems with the same matrix?

Apologies for not providing a copy of the relevant bits of code -- this is
my first email to the group, and a rather long one already -- happy to be
more concrete about anything I've said above.

Thanks
Dan

-- 

Daniel Goldberg, PhD
Lecturer in Glaciology
School of Geosciences, University of Edinburgh
Geography Building, Drummond Street, Edinburgh EH8 9XP


em: D <dgoldber at mit.edu>an.Goldberg at ed.ac.uk
web: http://ocean.mit.edu/~dgoldberg
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150222/d366f05b/attachment.html>