ident_zeros fails with PETSc.

[Douglas Brinkerhoff]

Using the development version of dolfin, when calling ident_zeros on any matrix that actually has null rows, PETSc throws an "Argument out of range" error. For example

from dolfin import *
mesh = UnitSquare(10,10)
Q = FunctionSpace(mesh,"CG",1)
u = TrialFunction(Q)
phi = TestFunction(Q)
A = u*phi*ds
a = assemble(A)
a.ident_zeros()

yields:

[0]PETSC ERROR: Configure options COPTFLAGS=-O2 --with-debugging=0 --with-shared-libraries=1 --with-clanguage=cxx --with-parmetis=1 --with-parmetis-dir=/usr --download-umfpack=1 --download-hypre=1 --download-mumps=1 --download-scalapack=1 --download-blacs=1 --download-ptscotch=1 --download-scotch=1 --prefix=/home/brinkerhoff/Source/FEniCS
[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: MatSetValues_SeqAIJ() line 331 in /home/brinkerhoff/Source/FEniCS/src/petsc-3.2-p6/src/mat/impls/aij/seq/aij.c
[0]PETSC ERROR: MatSetValues() line 1115 in /home/brinkerhoff/Source/FEniCS/src/petsc-3.2-p6/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Argument out of range!
[0]PETSC ERROR: New nonzero at (94,94) caused a malloc!
[0]PETSC ERROR: ------------------------------------------------------------------------

for every non-zero row. The above code worked with the previous development version.

[Timothée Lecomte (timothee-lecomte-5)]

I'a

[Timothée Lecomte (timothee-lecomte-5)]

I am also affected by that same bug. It is very annoying since I am dolfin to solve a boundary optimisation problem, so that the matrices are actually mostly zeros on the diagonal. The error reported is displayed thousands of times, so many times that the console is flooded and rendered completely useless...

[Timothée Lecomte (timothee-lecomte-5)]

I forgot to say that I see the problem when running on Ubuntu 12.10 x86, with dolfin from the launchpad fenics-dev PPA. Dolfin is currently at rev 7238, and petsc version from ubuntu is 3.2.0, patch 5.

[Jan Blechta (blechta)]

Present also on
Petsc Release Version 3.2.0, Patch 7, Thu Mar 15 09:30:51 CDT 2012
DOLFIN rev 7352

[Jan Blechta (blechta)]

Fix in lp:~blechta/dolfin/ident at rev 7435

[Cian Wilson (cwilson)]

Is there not already a fix for this in the trunk that doesn't require turning automatically turning off the petsc warning/error (which I find quite useful as it is flagging a significant cost if I hit it accidentally)... the keep_diagonal member of the AssemblerBase class? This ensures that the diagonal entries are not collapsed out of the sparsity pattern during initialization of the tensor.

[Jan Blechta (blechta)]

You're right, keep_diagonal fixes the problem. So I exposed keep_diagonal to python and push to lp:~blechta/dolfin/keep_diagonal

lp:~blechta/dolfin/ident could be merged to trunk also, because it simply makes PETScMatrix::ident() and ident_zeros() running on arbitrary PETScMatrix

[Jan Blechta (blechta)]

Moreover lp:~blechta/dolfin/ident fix runs on 1.1.0

[Jan Blechta (blechta)]

Should be fixed now: responsible PETSc error should be reported by DOLFIN with a hint to use keep_diagonal=True assemble option.

[Garth Wells (garth-wells)]

The sample code is flawed. Its attempts to modify the sparse matrix after it has been finalised.

[Jan Blechta (blechta)]

On Sat, 30 Mar 2013 08:08:22 -0000
Garth Wells <email address hidden> wrote:
> The sample code is flawed. Its attempts to modify the sparse matrix
> after it has been finalised.
>
> ** Changed in: dolfin
> Status: Fix Committed => New
>

Where's the problem? Which sample code is flawed?

[Garth Wells (garth-wells)]

On Saturday, 30 March 2013, Jan Blechta wrote:

> On Sat, 30 Mar 2013 08:08:22 -0000
> Garth Wells <<email address hidden> <javascript:;>> wrote:
> > The sample code is flawed. Its attempts to modify the sparse matrix
> > after it has been finalised.
> >
> > ** Changed in: dolfin
> > Status: Fix Committed => New
> >
>
> Where's the problem? Which sample code is flawed?

The code in the original post.

The problems with the merge request code are outlined in response to the
merge request.

Garth

>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1033611
>
> Title:
> ident_zeros fails with PETSc.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1033611/+subscriptions
>

--
Garth N. Wells
Department of Engineering, University of Cambridge
http://www.eng.cam.ac.uk/~gnw20

[Jan Blechta (blechta)]

Changed status to Fix Released as PETScMatrix::ident_zeros() works fine for square matrices in DOLFIN 1.2.0.

It only does not have unit test and assembler with keep_diagonal does not work for matrices with more rows than columns which is another, not reported bug (now being fixed by lp:~blechta/dolfin/keep_diagonal)

[Garth Wells (garth-wells)]

This is not fixed. The sample code gives an error:

[0]PETSC ERROR: --------------------- Error Message ------------------------------------
[0]PETSC ERROR: Object is in wrong state!
[0]PETSC ERROR: Matrix is missing diagonal entry in row 12!
[0]PETSC ERROR: ------------------------------------------------------------------------

It may work with older version of PETSc (<3.3), but performance will be poor.

[Jan Blechta (blechta)]

On Sun, 31 Mar 2013 13:21:42 -0000
Garth Wells <email address hidden> wrote:
> This is not fixed. The sample code gives an error:
>
>
> [0]PETSC ERROR: --------------------- Error Message
> ------------------------------------ [0]PETSC ERROR: Object is in
> wrong state! [0]PETSC ERROR: Matrix is missing diagonal entry in row
> 12! [0]PETSC ERROR:
> ------------------------------------------------------------------------
>
> It may work with older version of PETSc (<3.3), but performance will
> be poor.
>
> ** Changed in: dolfin
> Status: Fix Released => New
>

Garth, you did not understand the fix! It does not insert non-allocated
diagonal entries. It raises dolfin_error which suggest to allocate
diagonal entries (for example by keep_diagnal=True). This works with
PETSc 3.2 with no performance overhead! Please read my merge proposal
comment - inserting non-alocated entries is NOT allowed by DOLFIN.

Does this raise dolfin_error in addition to PETSc error with PETSc 3.3?
_______________________________
from dolfin import *
mesh = UnitSquare(10,10)
Q = FunctionSpace(mesh,"CG",1)
u = TrialFunction(Q)
phi = TestFunction(Q)
A = u*phi*ds
a = assemble(A)
a.ident_zeros()
_______________________________

Does this pass with PETSc 3.3?
_______________________________
from dolfin import *
mesh = UnitSquare(10,10)
Q = FunctionSpace(mesh,"CG",1)
u = TrialFunction(Q)
phi = TestFunction(Q)
A = u*phi*ds
a = assemble(A, keep_diagonal=True)
a.ident_zeros()
_______________________________

If answers are both yes then it is probably fixed. Definite
answer gives new test/unit/la/python/Matrix.py

Jan

[Garth Wells (garth-wells)]

On 31 March 2013 21:45, Jan Blechta <email address hidden> wrote:
> On Sun, 31 Mar 2013 13:21:42 -0000
> Garth Wells <email address hidden> wrote:
>> This is not fixed. The sample code gives an error:
>>
>>
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------ [0]PETSC ERROR: Object is in
>> wrong state! [0]PETSC ERROR: Matrix is missing diagonal entry in row
>> 12! [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>>
>> It may work with older version of PETSc (<3.3), but performance will
>> be poor.
>>
>> ** Changed in: dolfin
>> Status: Fix Released => New
>>
>
> Garth, you did not understand the fix!

OK, I see the error message now. When I test my terminal is filled
with PETSc errors which obscured the DOLFIN error.

It would be nice, if possible, to check the allocation and throw a
DOLFIN error before getting the screens full of PETSc error messages.

Garth

> It does not insert non-allocated
> diagonal entries. It raises dolfin_error which suggest to allocate
> diagonal entries (for example by keep_diagnal=True). This works with
> PETSc 3.2 with no performance overhead! Please read my merge proposal
> comment - inserting non-alocated entries is NOT allowed by DOLFIN.
>
> Does this raise dolfin_error in addition to PETSc error with PETSc 3.3?
> _______________________________
> from dolfin import *
> mesh = UnitSquare(10,10)
> Q = FunctionSpace(mesh,"CG",1)
> u = TrialFunction(Q)
> phi = TestFunction(Q)
> A = u*phi*ds
> a = assemble(A)
> a.ident_zeros()
> _______________________________
>
> Does this pass with PETSc 3.3?
> _______________________________
> from dolfin import *
> mesh = UnitSquare(10,10)
> Q = FunctionSpace(mesh,"CG",1)
> u = TrialFunction(Q)
> phi = TestFunction(Q)
> A = u*phi*ds
> a = assemble(A, keep_diagonal=True)
> a.ident_zeros()
> _______________________________
>
> If answers are both yes then it is probably fixed. Definite
> answer gives new test/unit/la/python/Matrix.py
>
> Jan

[Jan Blechta (blechta)]

On Sun, 31 Mar 2013 21:56:00 +0800
"Garth N. Wells" <email address hidden> wrote:
> On 31 March 2013 21:45, Jan Blechta <email address hidden>
> wrote:
> > On Sun, 31 Mar 2013 13:21:42 -0000
> > Garth Wells <email address hidden> wrote:
> >> This is not fixed. The sample code gives an error:
> >>
> >>
> >> [0]PETSC ERROR: --------------------- Error Message
> >> ------------------------------------ [0]PETSC ERROR: Object is in
> >> wrong state! [0]PETSC ERROR: Matrix is missing diagonal entry in
> >> row 12! [0]PETSC ERROR:
> >> ------------------------------------------------------------------------
> >>
> >> It may work with older version of PETSc (<3.3), but performance
> >> will be poor.
> >>
> >> ** Changed in: dolfin
> >> Status: Fix Released => New
> >>
> >
> > Garth, you did not understand the fix!
>
> OK, I see the error message now. When I test my terminal is filled
> with PETSc errors which obscured the DOLFIN error.

dolfin_error is not always last on stdout?

>
> It would be nice, if possible, to check the allocation and throw a
> DOLFIN error before getting the screens full of PETSc error messages.

Sure, it would be nice if PETSc error disappeared. But checking is
unnecessary overhead. It seems that stderr of PETSc could be somehow
redirected conditionally (for error catched by dolfin)
http://www.mcs.anl.gov/petsc/documentation/faq.html#redirectstdout
but I don't have a time now:(

Jan

>
> Garth
>
>
> > It does not insert non-allocated
> > diagonal entries. It raises dolfin_error which suggest to allocate
> > diagonal entries (for example by keep_diagnal=True). This works with
> > PETSc 3.2 with no performance overhead! Please read my merge
> > proposal comment - inserting non-alocated entries is NOT allowed by
> > DOLFIN.
> >
> > Does this raise dolfin_error in addition to PETSc error with PETSc
> > 3.3? _______________________________
> > from dolfin import *
> > mesh = UnitSquare(10,10)
> > Q = FunctionSpace(mesh,"CG",1)
> > u = TrialFunction(Q)
> > phi = TestFunction(Q)
> > A = u*phi*ds
> > a = assemble(A)
> > a.ident_zeros()
> > _______________________________
> >
> > Does this pass with PETSc 3.3?
> > _______________________________
> > from dolfin import *
> > mesh = UnitSquare(10,10)
> > Q = FunctionSpace(mesh,"CG",1)
> > u = TrialFunction(Q)
> > phi = TestFunction(Q)
> > A = u*phi*ds
> > a = assemble(A, keep_diagonal=True)
> > a.ident_zeros()
> > _______________________________
> >
> > If answers are both yes then it is probably fixed. Definite
> > answer gives new test/unit/la/python/Matrix.py
> >
> > Jan