(Fast?) assembly of upper triangle and MatZeroRows
Hi, I want my FE-program to be able to use LU- and LDLt-decomp (MUMPS). LU works fine and now I want to implement LDLt. My Questions: 1.) How can I speed up the assembly of the upper triangle matrix? * '-mat_ignore_lower_triangular' makes it terrible slow (1000 elements -> 1s explode to 18s * MatSetValues with Blocks of 6x6 and single values at the diagonal is even worse. 2.) SBAIJ doesn't like MatZeroRows. What are the alternatives to apply boundary conditions to my matrix? Thx for your help - Clemens
On Mon, Aug 8, 2011 at 08:40, Clemens Domanig <[email protected]>wrote:
Hi,
I want my FE-program to be able to use LU- and LDLt-decomp (MUMPS). LU works fine and now I want to implement LDLt. My Questions: 1.) How can I speed up the assembly of the upper triangle matrix? * '-mat_ignore_lower_triangular' makes it terrible slow (1000 elements -> 1s explode to 18s * MatSetValues with Blocks of 6x6 and single values at the diagonal is even worse.
SBAIJ needs different preallocation. If you set it with Mat{Seq,MPI}SBAIJSetPreallocation(), assembly will be fast.
2.) SBAIJ doesn't like MatZeroRows. What are the alternatives to apply boundary conditions to my matrix?
That function would make the matrix nonsymmetric because the columns aren't also zeroed, so you couldn't use it with a symmetric solver anyway. I believe in imposing boundary conditions symmetrically at the local level (either by removing them or by imposing them while computing residuals for those nodes), but many people like to do it differently. To do boundary conditions by zeroing rows and columns after assembly, you would need to modify the right hand side in a compatible way. We could automate this by having the Mat "remember" the entries in the columns that it zeroed so that each right hand side could automatically be modified before starting the solve, but then solving A*x = b would either yield an x that doesn't satisfy this equation or modify b, both of which I think would be confusing. Suggestions for imposing boundary conditions globally without ruining symmetry or having confusing semantics would be welcome.
You can use petsc-dev AIJ matrix format for mumps' LDLt-decomp. See http://www.mcs.anl.gov/petsc/petsc-as/developers/index.html on how to get petsc-dev. Hong On Mon, Aug 8, 2011 at 9:10 AM, Jed Brown <[email protected]> wrote:
On Mon, Aug 8, 2011 at 08:40, Clemens Domanig <[email protected]> wrote:
Hi,
I want my FE-program to be able to use LU- and LDLt-decomp (MUMPS). LU works fine and now I want to implement LDLt. My Questions: 1.) How can I speed up the assembly of the upper triangle matrix? * '-mat_ignore_lower_triangular' makes it terrible slow (1000 elements -> 1s explode to 18s * MatSetValues with Blocks of 6x6 and single values at the diagonal is even worse.
SBAIJ needs different preallocation. If you set it with Mat{Seq,MPI}SBAIJSetPreallocation(), assembly will be fast.
2.) SBAIJ doesn't like MatZeroRows. What are the alternatives to apply boundary conditions to my matrix?
That function would make the matrix nonsymmetric because the columns aren't also zeroed, so you couldn't use it with a symmetric solver anyway. I believe in imposing boundary conditions symmetrically at the local level (either by removing them or by imposing them while computing residuals for those nodes), but many people like to do it differently. To do boundary conditions by zeroing rows and columns after assembly, you would need to modify the right hand side in a compatible way. We could automate this by having the Mat "remember" the entries in the columns that it zeroed so that each right hand side could automatically be modified before starting the solve, but then solving A*x = b would either yield an x that doesn't satisfy this equation or modify b, both of which I think would be confusing. Suggestions for imposing boundary conditions globally without ruining symmetry or having confusing semantics would be welcome.
PETSc developers, We have a MatZeroRowsColumns() for MPIAIJ and SeqSBAIJ, what about having one for MPISBAIJ? Anyone want to code it up? Barry On Aug 8, 2011, at 9:10 AM, Jed Brown wrote:
On Mon, Aug 8, 2011 at 08:40, Clemens Domanig <[email protected]> wrote: Hi,
I want my FE-program to be able to use LU- and LDLt-decomp (MUMPS). LU works fine and now I want to implement LDLt. My Questions: 1.) How can I speed up the assembly of the upper triangle matrix? * '-mat_ignore_lower_triangular' makes it terrible slow (1000 elements -> 1s explode to 18s * MatSetValues with Blocks of 6x6 and single values at the diagonal is even worse.
SBAIJ needs different preallocation. If you set it with Mat{Seq,MPI}SBAIJSetPreallocation(), assembly will be fast.
2.) SBAIJ doesn't like MatZeroRows. What are the alternatives to apply boundary conditions to my matrix?
That function would make the matrix nonsymmetric because the columns aren't also zeroed, so you couldn't use it with a symmetric solver anyway. I believe in imposing boundary conditions symmetrically at the local level (either by removing them or by imposing them while computing residuals for those nodes), but many people like to do it differently.
To do boundary conditions by zeroing rows and columns after assembly, you would need to modify the right hand side in a compatible way. We could automate this by having the Mat "remember" the entries in the columns that it zeroed so that each right hand side could automatically be modified before starting the solve, but then solving A*x = b would either yield an x that doesn't satisfy this equation or modify b, both of which I think would be confusing.
Suggestions for imposing boundary conditions globally without ruining symmetry or having confusing semantics would be welcome.
On Mon, Aug 8, 2011 at 14:23, Barry Smith <[email protected]> wrote:
We have a MatZeroRowsColumns() for MPIAIJ and SeqSBAIJ, what about having one for MPISBAIJ?
Of course this makes sense to have, but I don't think it solves the problem. How do you solve more than one system? How would you use this interface with SNES?
On Aug 8, 2011, at 2:37 PM, Jed Brown wrote:
On Mon, Aug 8, 2011 at 14:23, Barry Smith <[email protected]> wrote: We have a MatZeroRowsColumns() for MPIAIJ and SeqSBAIJ, what about having one for MPISBAIJ?
Of course this makes sense to have, but I don't think it solves the problem. How do you solve more than one system?
You don't
How would you use this interface with SNES?
The user calls this as part the FormJacobian() (in the same way as if they explicitly remove those rows and columns (what Matt likes) in the matrix assembly. Barry
On Mon, Aug 8, 2011 at 14:44, Barry Smith <[email protected]> wrote:
How would you use this interface with SNES?
The user calls this as part the FormJacobian() (in the same way as if they explicitly remove those rows and columns (what Matt likes) in the matrix assembly.
But they need access to the right hand side at the same time. It seems too fragile to reach into the SNES and modify the residual vector during Jacobian evaluation.
On Aug 8, 2011, at 2:54 PM, Jed Brown wrote:
On Mon, Aug 8, 2011 at 14:44, Barry Smith <[email protected]> wrote:
How would you use this interface with SNES?
The user calls this as part the FormJacobian() (in the same way as if they explicitly remove those rows and columns (what Matt likes) in the matrix assembly.
But they need access to the right hand side at the same time. It seems too fragile to reach into the SNES and modify the residual vector during Jacobian evaluation.
Add a KSP/SNES./TSSetDirichlet(ksp,snes,ts,IS listofnodes,Vec valuesfornodes)? Barry
On Mon, Aug 8, 2011 at 17:16, Barry Smith <[email protected]> wrote:
Add a KSP/SNES./TSSetDirichlet(ksp,snes,ts,IS listofnodes,Vec valuesfornodes)?
This might be worthwhile. It needs a new Mat primitive which zeros the rows and columns while also "remembering" the column values that it zeroed (as another assembled matrix). I think it's best to apply the transformation immediately after residuals are calculated so that norms are correct (not in KSP at all). Matrix values are not available at this time so the Dirichlet equations may not be scaled correctly for GMG. We have to be careful that the residual evaluation passed to MatMFFD has the Dirichlet projection inside it.
On Mon, Aug 8, 2011 at 10:31 PM, Jed Brown <[email protected]> wrote:
On Mon, Aug 8, 2011 at 17:16, Barry Smith <[email protected]> wrote:
Add a KSP/SNES./TSSetDirichlet(ksp,snes,ts,IS listofnodes,Vec valuesfornodes)?
This might be worthwhile. It needs a new Mat primitive which zeros the rows and columns while also "remembering" the column values that it zeroed (as another assembled matrix). I think it's best to apply the transformation immediately after residuals are calculated so that norms are correct (not in KSP at all). Matrix values are not available at this time so the Dirichlet equations may not be scaled correctly for GMG. We have to be careful that the residual evaluation passed to MatMFFD has the Dirichlet projection inside it.
This is exactly the reason that I abandoned this approach when I got to nonlinear equations. It is complex and fragile, and requires that too many things be coordinated out of sight of the user. Matt -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
On Mon, Aug 8, 2011 at 17:31, Jed Brown <[email protected]> wrote:
It needs a new Mat primitive which zeros the rows and columns while also "remembering" the column values that it zeroed (as another assembled matrix).
I spoke too soon. If you apply the transformation while evaluating residuals and during assembly, there is no need to have this after-the-fact modification. Is there a reason this information should not be held by the DM? In Dohp, I use explicit projections into homogeneous and inhomogeneous spaces. Both are needed, so it's not clear to me that this can cleanly be placed in SNES or KSP without the weird semantics (modification of force vector or matrix, or having MatMult tell us that we didn't solve the equations. The concept of Dirichlet values (with homogeneous and inhomogeneous spaces) could be pushed all the way down into Mat and Vec (some of which I did for Dohp, in order to use VecGhost).
participants (5)
-
Barry Smith -
Clemens Domanig -
Hong Zhang -
Jed Brown -
Matthew Knepley