uBLAS GMRES solver fails on first run

[Garth Wells]

(reported by Doug Arnold)

After cleaning the instant cache, the solver

from dolfin import *
parameters["linear_algebra_backend"] = "uBLAS"
mesh = UnitSquare(3,3)
V = FunctionSpace(mesh, 'CG', 1)
v = TestFunction(V)
u = TrialFunction(V)
a = u*v*dx
A = assemble(a)
L = v*dx
b = assemble(L)
x = Vector(b.size())
solve(A, x, b, "gmres")
print x.array()

crashes with the message:

Check failed in file /usr/include/boost/numeric/ublas/operation.hpp at line 260:
norm_1 (v - cv) <= 2 * std::numeric_limits<real_type>::epsilon () * verrorbound
Traceback (most recent call last):
  File "solve3.py", line 12, in <module>
    solve(A, x, b, 'gmres')
  File "/home/faculty/arnold/workshop/FEniCS/lib/python2.6/site-packages/dolfin/cpp.py", line 4169, in solve
    return _cpp.solve(*args)
StandardError: internal logic

Searching on the web for the error message came up with this:
  http://archives.free.net.ph/message/20081022.224821.1eddeaa9.da.html
where Gunter Winkler suggests "#define BOOST_UBLAS_TYPE_CHECK 0",
but I don't know how to do this from dolfin.

[Johannes Ring (johannr)]

I couldn't reproduce this problem either on my laptop or on any of the buildbots. However, on my desktop computer at home (running Lucid 64 bit), I ran into the same issue. It happens about half the times when I run the solver repeatedly. I tried to reconfigure DOLFIN like this:

  scons configure customCxxFlags="-DBOOST_UBLAS_TYPE_CHECK=0"

This got rid of the error message but sometimes the solver completely hangs and then I have to kill the process:

  johannr@ubuntu-htpc:/tmp$ python solve3.py
  Matrix of size 16 x 16 has 82 nonzero entries.
  Sorting sparsity pattern.
  Solving linear system of size 16 x 16 (uBLAS Krylov solver).
  Terminated
  johannr@ubuntu-htpc:/tmp$

Anyone else experiencing this behaviour when DOLFIN is configured with the flag above?

[Anders Logg (logg)]

I get the same error message. Will check what goes wrong.

[Anders Logg (logg)]

And it happens only occasionally, not every time. And it doesn't seem to correlate with instant-clean.

[Anders Logg (logg)]

The problem was that x = Vector(b.size()) created a vector which sometimes contained strange numbers. This has been changed now (in uBLASVector::resize()) so that the vector is always zero. Took me a while to work through all the 29 different axpy_prod that call each other in operation.hpp in uBLAS...

Should be fixed now.

[Mehdi Nikbakht (m-nikbakht)]

Is there any specific reason why we have this type of problems in
resizing uBLAS vectors not other type of vectors?

Note that resizing a vector becomes important in the adaptive methods.
Since when a function updates, we need to resize its underlying vector.
If we clear a vector, then we lose all data that we had there.

Mehdi

On Mon, 2010-06-28 at 15:01 +0000, Anders Logg wrote:
> The problem was that x = Vector(b.size()) created a vector which
> sometimes contained strange numbers. This has been changed now (in
> uBLASVector::resize()) so that the vector is always zero. Took me a
> while to work through all the 29 different axpy_prod that call each
> other in operation.hpp in uBLAS...
>
> Should be fixed now.
>
> ** Changed in: dolfin
> Status: Confirmed => Fix Committed
>

[Anders Logg (logg)]

On Tue, Jun 29, 2010 at 11:50:38AM +0200, Mehdi wrote:
> Is there any specific reason why we have this type of problems in
> resizing uBLAS vectors not other type of vectors?

Yes, because of a bug in the uBLAS resize function (it does not
initialize data).

> Note that resizing a vector becomes important in the adaptive methods.
> Since when a function updates, we need to resize its underlying vector.
> If we clear a vector, then we lose all data that we had there.

There's no point in keeping the old dofs as part of the vector since
the number of the mesh (and the old dofs) may change. What is needed
for adaptive methods is to interpolate the old solution to the new
mesh (which is already implemented in DOLFIN).

--
Anders

> Mehdi
>
> On Mon, 2010-06-28 at 15:01 +0000, Anders Logg wrote:
> > The problem was that x = Vector(b.size()) created a vector which
> > sometimes contained strange numbers. This has been changed now (in
> > uBLASVector::resize()) so that the vector is always zero. Took me a
> > while to work through all the 29 different axpy_prod that call each
> > other in operation.hpp in uBLAS...
> >
> > Should be fixed now.
> >
> > ** Changed in: dolfin
> > Status: Confirmed => Fix Committed
> >
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp

[Mehdi Nikbakht (m-nikbakht)]

On Tue, 2010-06-29 at 12:13 +0200, Anders Logg wrote:
> On Tue, Jun 29, 2010 at 11:50:38AM +0200, Mehdi wrote:
> > Is there any specific reason why we have this type of problems in
> > resizing uBLAS vectors not other type of vectors?
>
> Yes, because of a bug in the uBLAS resize function (it does not
> initialize data).
>
> > Note that resizing a vector becomes important in the adaptive methods.
> > Since when a function updates, we need to resize its underlying vector.
> > If we clear a vector, then we lose all data that we had there.
>
> There's no point in keeping the old dofs as part of the vector since
> the number of the mesh (and the old dofs) may change. What is needed
> for adaptive methods is to interpolate the old solution to the new
> mesh (which is already implemented in DOLFIN).

This is not the case in XFEM where we have a fixed mesh but the number
of degrees of freedom is changing. Although this is not a big challenge,
since we can define some local variables to transfer data from the old
vector to new one.

Mehdi
>
> --
> Anders
>
>
> > Mehdi
> >
> > On Mon, 2010-06-28 at 15:01 +0000, Anders Logg wrote:
> > > The problem was that x = Vector(b.size()) created a vector which
> > > sometimes contained strange numbers. This has been changed now (in
> > > uBLASVector::resize()) so that the vector is always zero. Took me a
> > > while to work through all the 29 different axpy_prod that call each
> > > other in operation.hpp in uBLAS...
> > >
> > > Should be fixed now.
> > >
> > > ** Changed in: dolfin
> > > Status: Confirmed => Fix Committed
> > >
> >
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~dolfin
> > Post to : <email address hidden>
> > Unsubscribe : https://launchpad.net/~dolfin
> > More help : https://help.launchpad.net/ListHelp

[Garth Wells (garth-wells)]

On 29/06/10 14:33, Mehdi wrote:
> On Tue, 2010-06-29 at 12:13 +0200, Anders Logg wrote:
>> On Tue, Jun 29, 2010 at 11:50:38AM +0200, Mehdi wrote:
>>> Is there any specific reason why we have this type of problems in
>>> resizing uBLAS vectors not other type of vectors?
>>
>> Yes, because of a bug in the uBLAS resize function (it does not
>> initialize data).
>>
>>> Note that resizing a vector becomes important in the adaptive methods.
>>> Since when a function updates, we need to resize its underlying vector.
>>> If we clear a vector, then we lose all data that we had there.
>>
>> There's no point in keeping the old dofs as part of the vector since
>> the number of the mesh (and the old dofs) may change. What is needed
>> for adaptive methods is to interpolate the old solution to the new
>> mesh (which is already implemented in DOLFIN).
>
> This is not the case in XFEM where we have a fixed mesh but the number
> of degrees of freedom is changing. Although this is not a big challenge,
> since we can define some local variables to transfer data from the old
> vector to new one.
>

No all backends support preservation of data when resizing, so it's not
an option to preserve data when calling GenericVector::resize.

Garth

> Mehdi
>>
>> --
>> Anders
>>
>>
>>> Mehdi
>>>
>>> On Mon, 2010-06-28 at 15:01 +0000, Anders Logg wrote:
>>>> The problem was that x = Vector(b.size()) created a vector which
>>>> sometimes contained strange numbers. This has been changed now (in
>>>> uBLASVector::resize()) so that the vector is always zero. Took me a
>>>> while to work through all the 29 different axpy_prod that call each
>>>> other in operation.hpp in uBLAS...
>>>>
>>>> Should be fixed now.
>>>>
>>>> ** Changed in: dolfin
>>>> Status: Confirmed => Fix Committed
>>>>
>>>
>>>
>>> _______________________________________________
>>> Mailing list: https://launchpad.net/~dolfin
>>> Post to : <email address hidden>
>>> Unsubscribe : https://launchpad.net/~dolfin
>>> More help : https://help.launchpad.net/ListHelp
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~dolfin
> More help : https://help.launchpad.net/ListHelp