PETSc sparsity

[Barry Smith]

  Toby,

    For Seq matrices you can call MatGetRowIJ() and it will give
you exactly what you want. For MPI matrices the matrices are stored
in two parts on each process (the "diagonal" part and the "off-diagonal" part)
you can, of course, call the MatGetRowIJ() on each part but this is
not exactly what you are looking for.

  For algorithms that require dealing with the sparsity structure of
the matrix we generally just include the appropriate private include file
for the matrix format and access the data directly in the underlying format.
This allows the fastest performance.

  Other approaches to handling boundary conditions are to eliminated
fixed boundary values at finite element assembly time. That is to 
NEVER form the matrix with all unknowns, instead form the matrix without
the fixed values and adjust the right hand side load values at this time

Or use MatGetSubMatrix() to pull out the interior unknowns, another 
MatGetSubMatrix() to get the part of the matrix coupling the interior and
boundary unknowns and do a matrix vector product to undate the interior
unknown right hand side values based on the fixed unknowns. MatGetSubMatrix()
is actually fast and this is an efficient and easy way to handle the 
boundary conditions.

   Barry


On Fri, 15 Jun 2007, Toby Young wrote:

> 
> 
> Hello petsc-users and petsc-dev.
> 
> I am using PETSc to solve finite element problems. In particular I am
> working with the dealii package which has C++ wrappers designed to
> interface with PETSc. I am working on developing some of these wrappers.
> 
> The difficulty I face is that, after applying boundary conditions to the
> matrices, these matrices return as nonsymmetric; I really need them to be
> symmetric! I can achieve this by performing the so-called Gauss
> elimination step; however this requires access to the sparsity pattern of
> the PETSc matrix.
> 
> "The gauss elimination step "per se" is not a problem. The problem is
> that we don't have access to the sparsity of the matrix, so it is not
> so easy to find out which entries are not zero in the petsc matrices."
> -- dealii developer
> 
> Is there *any* possible way to obtain the sparsity of the PETSc matrix?
> Maybe something like GetMatrixSparsity? I have done my research, but I
> cannot find any information on this so I am pretty clueless on how
> to proceed (sorry if I have missed something obvious).
> 
> Any comments or suggests would be greatly appreciated.
> 
> Best,
> 	Toby
> 
> -----
> 
> Toby D. Young (Adiunkt)
> Department of Computational Science
> Institute of Fundamental Technological Research
> Polish Academy of Science
> Room 206, ul. Swietokrzyska 21
> 00-049 Warszawa, POLAND
> 
> 
> 
> 
> Luca.
> 
> --
> Luca Heltai <luca.heltai at gmail.com>
> http://www-dimat.unipv.it/heltai
> --
> There are no answers, only cross references.
> 
> 
> 
> On Jun 15, 2007, at 11:17 AM, Toby D. Young wrote:
> 
> > On Thu, 14 Jun 2007 10:55:20 -0400
> > Luca Heltai <luca.heltai at gmail.com> wrote:
> >
> >> Toby, if you do some serious development (like what you are
> >> doing... :) ...) I strongly suggest you switch to the subversion
> >> tree.
> >
> > Okay, fair suggestion, and not too tricky to do, I'm running of dealii
> > 5.3pre now.
> >
> >> In the subversion tree that is implemented.
> >
> > Yet it doesn't seem to be and the problem does not go away. I'll
> > have to
> > think of a way around this, either that or actually implement Guass
> > elimination for the Petsc SParseMatrix < scary thought :-) >.
> >
> > You can verify this yourself by running step-17, after first changing
> > the final flag false->true in apply_boundary_conditions.
> >
> > Here's what we find:
> >
> > Cycle 0:
> >    Number of active cells        : 1024
> >    Number of degrees of freedom  : 1089 (by partition: 1089)
> > --------------------------------------------------------
> > An error occurred in line <534> of file
> > <source/numerics/matrices.all_dimensions.cc>
> > in function void dealii::PETScWrappers::apply_boundary_values(const
> > std::map<unsigned int, double, std::less<unsigned int>,
> > std::allocator<std::pair<const unsigned int, double> > >&,
> > PETScMatrix&, PETScVector&, PETScVector&, bool) [with PETScMatrix =
> > dealii::PETScWrappers::MPI::SparseMatrix, PETScVector =
> > dealii::PETScWrappers::MPI::Vector]
> > The violated condition was:
> > preserve_symmetry == false
> > The name and call sequence of the exception was:
> > ExcNotImplemented()
> > Additional Information:
> > (none)
> >
> > Stacktrace:
> > -----------
> > #0  /home/tyoung/fem/deal.II/lib/libdeal_II_1d.g.so: void
> > dealii::PETScWrappers::apply_boundary_values<dealii::PETScWrappers::MP
> > I::SparseMatrix,
> > dealii::PETScWrappers::MPI::Vector>(std::map<unsigned int, double,
> > std::less<unsigned int>, std::allocator<std::pair<unsigned int const,
> > double> > > const&, dealii::PETScWrappers::MPI::SparseMatrix&,
> > double> > > dealii::PETScWrappers::MPI::Vector&,
> > double> > > dealii::PETScWrappers::MPI::Vector&, bool)
> > double> > > #1  /home/tyoung/fem/deal.II/lib/libdeal_II_1d.g.so:
> > double> > > dealii::MatrixTools::apply_boundary_values
> > (std::map<unsigned
> > double> > > int, double, std::less<unsigned int>,
> > double> > > std::allocator<std::pair<unsigned int const, double> > >
> > double> > > const&, dealii::PETScWrappers::MPI::SparseMatrix&,
> > double> > > dealii::PETScWrappers::MPI::Vector&,
> > double> > > dealii::PETScWrappers::MPI::Vector&, bool)
> > #2  ./step-75: double> > > InfiniteWellProblem<2>::assemble_system()
> > #3  ./step-75: double> > > InfiniteWellProblem<2>::run()
> > #4  ./step-75: main
> > --------------------------------------------------------
> >
> >
> >
> > --
> >
> > Toby D. Young ( Adiunkt )
> > Department of Computational Science
> > Institute of Fundamental Technical Research, Polish Academy of
> > Sciences
> > Room 206, Swietokrzyska 21, 00-049 Warsaw, POLAND
> 
>