[petsc-dev] Mat/DM dependency
Sorry for jumping in late. PDE discretizations and solvers many geometric decompositions of the linear space they ultimately operate on. However, many of these decompositions can be expressed rather naturally in linear algebraic terms (e.g., factorization of the stiffness matrix into a sum of element matrices in FEM). The gather-scatter matrices that glue together the element matrices are largely independent of the operator and encode the decomposition of the global function space into elementwise pieces. This is merely a linear algebraic form of the underlying partition of unity and, again, can be cast in purely algebraic terms as a decomposition of the identity operator: I = \sum_i R_i P_i where P_i is the projector onto the ith piece and R_i is the injection back into the larger space. This is formalism that a lot of BDDC is built around (e.g., J. Mandel, C. Dohrmann, Numer. Linear Algebra Appl. 2003; 10:639?659). Automating decomposition/assembly of operators relative to a decomposition of identity (doi) would, in my opinion, go a long way towards decoupling Mat from DM. DM would then convert their geometric data into a corresponding doi and Mat wouldn't need to know about DMs directly. That's what I'm trying to do with MatFwk -- a framework for general block decompositions of matrices, but that's been moving rather slowly. One particularly interesting aspect is the need to deal with dual-primal decompositions that require elimination of constraints or introduction of Lagrange multipliers. Mortar elements are a relatively simple example of that. Incidentally, there are similar Mat/DM dependence problems with wavelet(-like) transforms, including USFFT, that we've been trying to introduce into PETSc. Dmitry. On Sat, Jan 2, 2010 at 12:48 PM, Matthew Knepley knepley at gmail.com wrote: On Sat, Jan 2, 2010 at 12:44 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Jan 1, 2010, at 7:12 PM, Jed Brown wrote: On Fri, 1 Jan 2010 14:50:33 -0600, Matthew Knepley knepley at gmail.com wrote: That is definitely a problem. ?I am not convinced yet that this is a problem. The dependency is only in one specific Mat implementation that should not require pulling in the DM library unless that specific Mat implementation is specifically used. Since that specific Mat implementation is not used in the Mat examples it should not require pulling in the DM library. If it is not a practical problem, I still believe it is a conceptual problem. ? Matt ? Barry It makes me think we need an R^N part just like we have now, and then a more general part with operators on fiber bundles. The question is whether it is important to separate that stuff (which depends on both Vec and DM) from DM (which already depends on Vec). Jed -- 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
[petsc-dev] Mat/DM dependency
On Tue, 5 Jan 2010 10:52:36 -0600, Dmitry Karpeev karpeev at mcs.anl.gov wrote: Automating decomposition/assembly of operators relative to a decomposition of identity (doi) would, in my opinion, go a long way towards decoupling Mat from DM. DM would then convert their geometric data into a corresponding doi and Mat wouldn't need to know about DMs directly. I like keeping Mat purely algebraic. When assembly has nontrivial semantics, I write a MyDMMatSetValuesBlockedWhatever(MyDM,Mat,...) so that Mat doesn't see DM (since the DM is explicitly available when this function is called, I prefer using it explicitly over patching the Mat to retain a reference to the DM). That's what I'm trying to do with MatFwk -- a framework for general block decompositions of matrices What granularity are you targeting with MatFwk? I think it could be very difficult to have good performance for very fine grains (e.g. unassembled P1 elements). Jed
[petsc-dev] Mat/DM dependency
On Tue, Jan 5, 2010 at 11:42 AM, Jed Brown jed at 59a2.org wrote:
On Tue, 5 Jan 2010 10:52:36 -0600, Dmitry Karpeev karpeev at mcs.anl.gov
wrote:
Automating decomposition/assembly of operators relative to a
decomposition of identity (doi) would, in my opinion, go a long way
towards decoupling Mat from DM. ?DM would then convert their geometric
data into a corresponding doi and Mat wouldn't need to know about DMs
directly.
I like keeping Mat purely algebraic. ?When assembly has nontrivial
semantics, I write a MyDMMatSetValuesBlockedWhatever(MyDM,Mat,...) so
that Mat doesn't see DM (since the DM is explicitly available when this
function is called, I prefer using it explicitly over patching the Mat
to retain a reference to the DM).
That's what I'm trying to do with MatFwk -- a framework for general
block decompositions of matrices
What granularity are you targeting with MatFwk? ?I think it could be
very difficult to have good performance for very fine grains (e.g.
unassembled P1 elements).
I definitely see that problem.
My current view is that MatFwk could operate in two modes: (1)
unassembled and (2) assembled.
The assembled mode (below) may be better suited for fine granularity
decompositions. Even then
it might be challenging to make it perform well.
1) In the unassembled mode MatFwk would hold relatively few Mat
objects {M_k}, implementing the matrix blocks
as well as the Mat objects implementing the corresponding
projection/injection operators {P_k, R_k}
(P_k and R_k can frequently, but not always, implemented with
VecScatters wrapped as Mats).
MatFwk would then call these operators to effect a MatMult operation:
Mu = \sum_k R_k M_k P_k u.
In a sense, this will unify various other ways to compose matrices in
PETSc (MatComposite, etc).
There would definitely be overhead associated with the internally-held
Mat objects and the additional
Vec storage for the intermediate results. This is nothing new,
however, relative to MatComposite,
for example. As long as the granularity is relatively coarse (few
matrix blocks), the performance penalty is
probably acceptable.
Such unassembled matrices would definitely not be usable with many (if
not all) of the standard preconditioners.
However, this would allow for almost any Mat type to be used to
implement the matrix blocks, including using
MatFwk recursively. Matrix blocks could be extracted and manipulated
(as long as the layout stays to same) to
set their values, etc.
MatFwk, in a sense, would define the block sparsity or block
connectivity structure for blocks of Vec dofs:
for each k, the input block consists of those dofs that span the
preimage of P_k (loosely speaking, but it can
be made precise). Likewise, the output block will be the dofs that
span the image of R_k. This way MatFwk
is akin to an aij structure for blocks of dofs.
2) In the assembled mode, the matrix blocks would be defined only
conceptually. The dof block sparsity structure still needs
to be defined as well as the ways to inject blocks of local dofs
into the larger Vec and vice versa (i.e., R_ks and P_ks),
but there would be no per block overhead in the form of separate
objects implementing each M_k, P_k, R_k triple.
We still need the {P_k, R_k}, but they could be stored in a single
structure, amortizing the overhead. Then the corresponding
M_ks can be manipulated by setting their values, which are then
injected into M using the R_ks. The advantage is that
the user need only know how to form these local values, which are then
correctly accumulated within M.
This is the model that FEM uses.
In fact, PETSc already provides for this via MatSetValues, except (a)
all R_ks have the form of ADD_VALUES and (b)
the storage of the mapping from local to global degrees of freedom is
the responsibility of the user (the block index arrays supplied
to MatSetValues). So an assembled MatFwk would store these block
index arrays internally (and, possibly, R_ks matrix entries,
if we need more than just ADD_VALUES, as is the case for BDDC, mortar
elements, etc) implying overhead.
HOWEVER, this information is already stored as mesh data:
element-to-vertex lists and vertex coordinates, in the P1 case.
Distilling it to P_ks/R_ks would only (i) make it algebraic, (ii) make
it internal to MatFwk as a resolution of identity.
This is definitely a solver-centric viewpoint, which subordinates the
mesh/geometry to Mat, but without making references to
the corresponding software abstractions at the expense of algebrizing them.
DM implementations could then be built on top of MatFwk, and they
would be responsible for converting geometry
to algebra (i.e., connectivity and continuity into
projection/injection operators).
Dmitry.
Jed
[petsc-dev] Mat/DM dependency
On Tue, 5 Jan 2010 14:24:36 -0600, Dmitry Karpeev karpeev at mcs.anl.gov wrote: So an assembled MatFwk would store these block index arrays internally (and, possibly, R_ks matrix entries, if we need more than just ADD_VALUES, as is the case for BDDC How is something other than ADD_VALUES necessary for BDDC? Once the primal space is identified, each entry in the global numbering becomes split between a local unassembled part and a global assembled part. So P_k,R_k are not 1-1, but we still use ADD_VALUES. Jed
[petsc-dev] Mat/DM dependency
On Sun, 3 Jan 2010 16:52:18 -0600, Barry Smith bsmith at mcs.anl.gov wrote: I have moved the hyprestruct matrix definition from mat directory to the dm/da directory. This should resolve the problem for you. Great, thanks. Jed
[petsc-dev] Mat/DM dependency
I have moved the hyprestruct matrix definition from mat directory to the dm/da directory. This should resolve the problem for you. Barry On Jan 2, 2010, at 11:28 PM, Jed Brown wrote: On Sat, 02 Jan 2010 21:39:13 -0700, Jed Brown jed at 59A2.org wrote: This is standard behavior, and I believe it's a good thing To elaborate on this, suppose libA provides the symbol AFunc() which has some weird dependency. Suppose libB does not reference AFunc() and some application appC depends on libB. Since AFunc() is unreachable when we link appC, a permissive linker will let the silent dependencies slip by unresolved. Now we update libB (without changing it's ABI) so that AFunc() is referenced. Since appC has not changed, we should be able to just update libB and appC should not need to be touched. But since the linker let the dependencies for the formerly unreachable AFunc() slip by when we linked appC, we now have a broken executable. A stricter linker would have forced us to resolve the dependencies up front so that we can't silently end up with this broken executable. Also, I'm missing something else if the present isn't very brittle with static libs due to the link order. That is, -lpetscdm -lpetscmat could leave DA symbols in libpetscmat.a unresolved, requiring something like -lpetscdm -lpetscmat -lpetscdm. (This only applies if the application didn't reference these symbols, which any reasonable application probably needs to do in this case, explaining why it hasn't caused a problem yet.) Note that we already have MatSchurComplement in libpetscksp. Jed
[petsc-dev] Mat/DM dependency
Does this mean that Linux shared libraries require that EVERYTHING in a shared library be resolved even if the executable does not have a chain of calls that would ever reach the functions in the shared library that cannot be resolved? This is certainly not true for static libraries. Please confirm that this DOES happen for static libraries? This problem is present with every build I have where Hypre is enabled. It shouldn't happen for static libraries so something is going on that I do not understand. Perhaps there is some linux specific link flag that tells the linker to stop being stupid and not resolve things that don't need to be resolved? Barry On Jan 1, 2010, at 7:09 PM, Jed Brown wrote: On Fri, 1 Jan 2010 16:15:54 -0600, Barry Smith bsmith at mcs.anl.gov wrote: On Jan 1, 2010, at 1:19 PM, Jed Brown wrote: Mat does not normally depend on DM, but a couple implementations do, $ grep -rl DAGet petsc/src/mat/impls petsc/src/mat/impls/hypre/mhyp.c petsc/src/mat/impls/adic/nladic.c petsc/src/mat/impls/adic/matadic.c So when these are enabled, all the Mat examples fail to link because PETSC_MAT_LIB does not include PETSC_DM_LIB. Please explain in detail when this occurs. It does not happen to me on the Apple with static, shared or dynamic libraries. Does it happen when linking the example or running it? Please cut and past all error messages that appear. This problem is present with every build I have where Hypre is enabled. I normally use shared, but not dynamic so the errors show up at link time. $ make info |grep '^Using configure Options' Using configure Options: --with-zoltan-dir=/usr --download-ml --with- blas-lapack-dir=/usr --download-blacs --download-chaco --with-mpi- dir=/usr --download-mumps --download-superlu_dist --download-spooles --with-parmetis-dir=/usr --download-hypre --with-c2html --with-hdf5- dir=/usr --with-umfpack-dir=/usr --download-sundials --download- scalapack --with-lgrind --with-shared --with-sowing -PETSC_ARCH=ompi --download-superlu --download-spai $ cd src/mat/examples/tutorials make ex1 /usr/bin/mpicc -o ex1.o -c -fPIC -Wall -Wwrite-strings -Wno-strict- aliasing -g3 -I/home/jed/petsc/ompi/include -I/home/jed/petsc/ include -I/usr/include -I/usr/lib/openmpi -I/home/jed/petsc/ompi/ include -D__SDIR__=src/mat/examples/tutorials/ ex1.c /usr/bin/mpicc -fPIC -Wall -Wwrite-strings -Wno-strict-aliasing - g3 -o ex1 ex1.o -Wl,-rpath,/home/jed/petsc/ompi/lib -L/home/jed/ petsc/ompi/lib -lpetscmat -lpetscvec -lpetsc-Wl,-rpath,/usr/lib - L/usr/lib -lzoltan -lX11 -Wl,-rpath,/home/jed/petsc/ompi/lib -L/home/ jed/petsc/ompi/lib -lcmumps -ldmumps -lsmumps -lzmumps - lmumps_common -lpord -lscalapack -lblacs -lchaco -lspooles -lHYPRE - Wl,-rpath,/usr/lib/openmpi -Wl,-rpath,/usr/lib/gcc/x86_64-unknown- linux-gnu/4.4.2 -lmpi_cxx -lstdc++ -lhdf5 -lz -lsuperlu_dist_2.3 - lparmetis -lmetis -lml -lmpi_cxx -lstdc++ -lspai -lsundials_cvode - lsundials_nvecserial -lsundials_nvecparallel -lsuperlu_4.0 -lumfpack -lamd -Wl,-rpath,/usr -L/usr -llapack -lcblas -lf77blas -latlas -Wl,- rpath,/usr/lib/openmpi -L/usr/lib/openmpi -Wl,-rpath,/usr/lib/gcc/ x86_64-unknown-linux-gnu/4.4.2 -L/usr/lib/gcc/x86_64-unknown-linux- gnu/4.4.2 -ldl -lmpi -lopen-rte -lopen-pal -lnsl -lutil -lgcc_s - lpthread -lmpi_f90 -lmpi_f77 -lgfortran -lm -lm -Wl,-rpath,/usr/lib/ gcc/x86_64-unknown-linux-gnu -L/usr/lib/gcc/x86_64-unknown-linux-gnu -lm -lm -lmpi_cxx -lstdc++ -lmpi_cxx -lstdc++ -ldl -lmpi -lopen-rte - lopen-pal -lnsl -lutil -lgcc_s -lpthread -ldl /home/jed/petsc/ompi/lib/libpetscmat.so: undefined reference to `DAGetInfo' /home/jed/petsc/ompi/lib/libpetscmat.so: undefined reference to `DAGetCorners' /home/jed/petsc/ompi/lib/libpetscmat.so: undefined reference to `DAGetGhostCorners' /home/jed/petsc/ompi/lib/libpetscmat.so: undefined reference to `DAGetGlobalIndices' collect2: ld returned 1 exit status make: [ex1] Error 1 (ignored) /bin/rm -f ex1.o Let's not get carried away with redesigning all of PETSc because of a simple link issue (or crappy shared library design on Linux :-). The stuff in mhyp.c uses DA, this is a dependency loop because (1) this file resides in src/mat and (2) it is linked into libpetscmat which we claim can be used without libpetscdm. Jed
[petsc-dev] Mat/DM dependency
On Jan 1, 2010, at 7:12 PM, Jed Brown wrote: On Fri, 1 Jan 2010 14:50:33 -0600, Matthew Knepley knepley at gmail.com wrote: That is definitely a problem. I am not convinced yet that this is a problem. The dependency is only in one specific Mat implementation that should not require pulling in the DM library unless that specific Mat implementation is specifically used. Since that specific Mat implementation is not used in the Mat examples it should not require pulling in the DM library. Barry It makes me think we need an R^N part just like we have now, and then a more general part with operators on fiber bundles. The question is whether it is important to separate that stuff (which depends on both Vec and DM) from DM (which already depends on Vec). Jed
[petsc-dev] Mat/DM dependency
On Sat, Jan 2, 2010 at 12:44 PM, Barry Smith bsmith at mcs.anl.gov wrote: On Jan 1, 2010, at 7:12 PM, Jed Brown wrote: On Fri, 1 Jan 2010 14:50:33 -0600, Matthew Knepley knepley at gmail.com wrote: That is definitely a problem. I am not convinced yet that this is a problem. The dependency is only in one specific Mat implementation that should not require pulling in the DM library unless that specific Mat implementation is specifically used. Since that specific Mat implementation is not used in the Mat examples it should not require pulling in the DM library. If it is not a practical problem, I still believe it is a conceptual problem. Matt Barry It makes me think we need an R^N part just like we have now, and then a more general part with operators on fiber bundles. The question is whether it is important to separate that stuff (which depends on both Vec and DM) from DM (which already depends on Vec). Jed -- 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 -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20100102/b65677cf/attachment.html
[petsc-dev] Mat/DM dependency
On Sat, 2 Jan 2010 12:41:33 -0600, Barry Smith bsmith at mcs.anl.gov wrote: Does this mean that Linux shared libraries require that EVERYTHING in a shared library be resolved even if the executable does not have a chain of calls that would ever reach the functions in the shared library that cannot be resolved? This is standard behavior, and I believe it's a good thing because everything in a library should be accessible through it's interface. You seem to be arguing that it's a good thing to have hidden stuff in a library that is only accessible by a magic extern declaration, but I don't see the benefit of this. Why not move mhyp.c to src/dm? Please confirm that this DOES happen for static libraries? It does not happen with static libs, the builds I use I keep around all use shared libs. Jed
[petsc-dev] Mat/DM dependency
On Sat, 02 Jan 2010 21:39:13 -0700, Jed Brown jed at 59A2.org wrote: This is standard behavior, and I believe it's a good thing To elaborate on this, suppose libA provides the symbol AFunc() which has some weird dependency. Suppose libB does not reference AFunc() and some application appC depends on libB. Since AFunc() is unreachable when we link appC, a permissive linker will let the silent dependencies slip by unresolved. Now we update libB (without changing it's ABI) so that AFunc() is referenced. Since appC has not changed, we should be able to just update libB and appC should not need to be touched. But since the linker let the dependencies for the formerly unreachable AFunc() slip by when we linked appC, we now have a broken executable. A stricter linker would have forced us to resolve the dependencies up front so that we can't silently end up with this broken executable. Also, I'm missing something else if the present isn't very brittle with static libs due to the link order. That is, -lpetscdm -lpetscmat could leave DA symbols in libpetscmat.a unresolved, requiring something like -lpetscdm -lpetscmat -lpetscdm. (This only applies if the application didn't reference these symbols, which any reasonable application probably needs to do in this case, explaining why it hasn't caused a problem yet.) Note that we already have MatSchurComplement in libpetscksp. Jed
[petsc-dev] Mat/DM dependency
I would rather see then hierarchy change. I think it is natural for the operator (Mat) to depend on the space (DM) on which it is discretized. Just because shortsightedness in the past has confined us to really simple spaces (R^N) does not mean we can't change that. Matt On Fri, Jan 1, 2010 at 1:19 PM, Jed Brown jed at 59a2.org wrote: Mat does not normally depend on DM, but a couple implementations do, $ grep -rl DAGet petsc/src/mat/impls petsc/src/mat/impls/hypre/mhyp.c petsc/src/mat/impls/adic/nladic.c petsc/src/mat/impls/adic/matadic.c So when these are enabled, all the Mat examples fail to link because PETSC_MAT_LIB does not include PETSC_DM_LIB. This is especially annoying to me because the PETSc builds I use have Hypre enabled, meaning that I need to keep more special testing builds (without Hypre) just to be able to run the Mat tests. Is the HYPREStruct stuff destined to eventually move out of Mat? Jed -- 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 -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20100101/958d9e6e/attachment.html
[petsc-dev] Mat/DM dependency
On Fri, 1 Jan 2010 13:32:08 -0600, Matthew Knepley knepley at gmail.com wrote: I would rather see then hierarchy change. Maybe, but many things in the DM interface depend on Mat and Vec. Where would DMGetInterpolation and DMGetMatrix end up? I think it is natural for the operator (Mat) to depend on the space (DM) on which it is discretized. Agreed. Just because shortsightedness in the past has confined us to really simple spaces (R^N) does not mean we can't change that. KSP operates in a finite dimensional linear space, i.e. isometric to R^n, so at least up to that level, we're not really exploiting the more general interpretation. Jed
[petsc-dev] Mat/DM dependency
On Fri, Jan 1, 2010 at 1:58 PM, Jed Brown jed at 59a2.org wrote: On Fri, 1 Jan 2010 13:32:08 -0600, Matthew Knepley knepley at gmail.com wrote: I would rather see then hierarchy change. Maybe, but many things in the DM interface depend on Mat and Vec. Where would DMGetInterpolation and DMGetMatrix end up? That is definitely a problem. It makes me think we need an R^N part just like we have now, and then a more general part with operators on fiber bundles. I think it is natural for the operator (Mat) to depend on the space (DM) on which it is discretized. Agreed. Just because shortsightedness in the past has confined us to really simple spaces (R^N) does not mean we can't change that. KSP operates in a finite dimensional linear space, i.e. isometric to R^n, so at least up to that level, we're not really exploiting the more general interpretation. It depends how you see it. KSP can operate just as well on a fiber bundle. Matt Jed -- 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 -- next part -- An HTML attachment was scrubbed... URL: http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20100101/ed5b376c/attachment.html
[petsc-dev] Mat/DM dependency
On Jan 1, 2010, at 1:19 PM, Jed Brown wrote: Mat does not normally depend on DM, but a couple implementations do, $ grep -rl DAGet petsc/src/mat/impls petsc/src/mat/impls/hypre/mhyp.c petsc/src/mat/impls/adic/nladic.c petsc/src/mat/impls/adic/matadic.c So when these are enabled, all the Mat examples fail to link because PETSC_MAT_LIB does not include PETSC_DM_LIB. Please explain in detail when this occurs. It does not happen to me on the Apple with static, shared or dynamic libraries. Does it happen when linking the example or running it? Please cut and past all error messages that appear. Barry Let's not get carried away with redesigning all of PETSc because of a simple link issue (or crappy shared library design on Linux :-). This is especially annoying to me because the PETSc builds I use have Hypre enabled, meaning that I need to keep more special testing builds (without Hypre) just to be able to run the Mat tests. Is the HYPREStruct stuff destined to eventually move out of Mat? Jed
[petsc-dev] Mat/DM dependency
On Fri, 1 Jan 2010 14:50:33 -0600, Matthew Knepley knepley at gmail.com wrote: That is definitely a problem. It makes me think we need an R^N part just like we have now, and then a more general part with operators on fiber bundles. The question is whether it is important to separate that stuff (which depends on both Vec and DM) from DM (which already depends on Vec). Jed
