# Dolfin GMRES vs. CBC.Block LGMRES

This question is maybe way too general, but perhaps you can give me some hints:

I want to solve a nonsymmetric problem with CBC.Block.
The standard dolfin solve method with Trilinos backend solves the system in 5 iterations:

A, b = assemble_system(a, L)
solve(A, x.vector(), b, "gmres", "ilu")

However trying to solve the system with CBC.Block fails:
Ap = ILU(A)
Ainv = LGMRES(A, precond=Ap)
x = Ainv * b

with the following:

File "/usr/local/lib/python2.7/dist-packages/block/iterative/lgmres.py", line 190, in lgmres
y, resids, rank, s = lstsq(hess, e1)
File "/usr/lib/python2.7/dist-packages/scipy/linalg/basic.py", line 409, in lstsq
a1, b1 = map(asarray_chkfinite, (a, b))
File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 590, in asarray_chkfinite
"array must not contain infs or NaNs")
ValueError: array must not contain infs or NaNs

I also tried BiCGStab as a solver (3 iterations with builtin Trilinos solver, NaN residuum in case of the CBC.Block implementation).

As far as I understand, both the Trilinos GMRES and the CBC.Block LGMRES use the IFPACK ILU preconditioner. So does this mean the only difference between the two libraries is the implementation of the iterative methods?

I want to use the CBC.Block solver, because I want to implement a Schur complement method where I have to solve another linear system per GMRES step.

Do you have any hints what I could do to make things work with CBC.Block? I know this question is very general. I will gladly provide more information if needed.

## Question information

Language:
English Edit question
Status:
Solved
For:
CBC.Block Edit question
Assignee:
No assignee Edit question
Solved by:
Claas Abert
Solved:
Last query:
 Revision history for this message Joachim Haga (jobh) said on 2013-01-28: #1

I do not have any ideas off-hand. I have used LGMRES without seing this
problem before. Which dolfin and cbc.block version are you running? Also,
if you could provide a runnable example that shows this error (complete,
not necessarily minimal) then I can have a look. Send to <email address hidden>.

On 28 January 2013 12:06, Claas Abert

> New question #220333 on CBC.Block:
>
> This question is maybe way too general, but perhaps you can give me some
> hints:
>
> I want to solve a nonsymmetric problem with CBC.Block.
> The standard dolfin solve method with Trilinos backend solves the system
> in 5 iterations:
>
> A, b = assemble_system(a, L)
> solve(A, x.vector(), b, "gmres", "ilu")
>
> However trying to solve the system with CBC.Block fails:
> Ap = ILU(A)
> Ainv = LGMRES(A, precond=Ap)
> x = Ainv * b
>
> with the following:
>
> File "/usr/local/lib/python2.7/dist-packages/block/iterative/lgmres.py",
> line 190, in lgmres
> y, resids, rank, s = lstsq(hess, e1)
> File "/usr/lib/python2.7/dist-packages/scipy/linalg/basic.py", line 409,
> in lstsq
> a1, b1 = map(asarray_chkfinite, (a, b))
> File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line
> 590, in asarray_chkfinite
> "array must not contain infs or NaNs")
> ValueError: array must not contain infs or NaNs
>
> I also tried BiCGStab as a solver (3 iterations with builtin Trilinos
> solver, NaN residuum in case of the CBC.Block implementation).
>
> As far as I understand, both the Trilinos GMRES and the CBC.Block LGMRES
> use the IFPACK ILU preconditioner. So does this mean the only difference
> between the two libraries is the implementation of the iterative methods?
>
> I want to use the CBC.Block solver, because I want to implement a Schur
> complement method where I have to solve another linear system per GMRES
> step.
>
> Do you have any hints what I could do to make things work with CBC.Block?
> I know this question is very general. I will gladly provide more
> information if needed.
>
>
> --
> You received this question notification because you are a member of
> cbc.block maintainers, which is an answer contact for CBC.Block.
>
> --
> Post to : <email address hidden>
>

 Revision history for this message Kent-Andre Mardal (kent-and) said on 2013-01-28: #2

I can also look at it, so please post to cbc.block if the code
is not unreasonably large

Kent

On 28 January 2013 12:15, Joachim Haga <<email address hidden>
> wrote:

> Question #220333 on CBC.Block changed:
>
>
> Joachim Haga proposed the following answer:
> I do not have any ideas off-hand. I have used LGMRES without seing this
> problem before. Which dolfin and cbc.block version are you running? Also,
> if you could provide a runnable example that shows this error (complete,
> not necessarily minimal) then I can have a look. Send to <email address hidden>.
>
>
> On 28 January 2013 12:06, Claas Abert
>
> > New question #220333 on CBC.Block:
> >
> > This question is maybe way too general, but perhaps you can give me some
> > hints:
> >
> > I want to solve a nonsymmetric problem with CBC.Block.
> > The standard dolfin solve method with Trilinos backend solves the system
> > in 5 iterations:
> >
> > A, b = assemble_system(a, L)
> > solve(A, x.vector(), b, "gmres", "ilu")
> >
> > However trying to solve the system with CBC.Block fails:
> > Ap = ILU(A)
> > Ainv = LGMRES(A, precond=Ap)
> > x = Ainv * b
> >
> > with the following:
> >
> > File "/usr/local/lib/python2.7/dist-packages/block/iterative/lgmres.py",
> > line 190, in lgmres
> > y, resids, rank, s = lstsq(hess, e1)
> > File "/usr/lib/python2.7/dist-packages/scipy/linalg/basic.py", line 409,
> > in lstsq
> > a1, b1 = map(asarray_chkfinite, (a, b))
> > File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line
> > 590, in asarray_chkfinite
> > "array must not contain infs or NaNs")
> > ValueError: array must not contain infs or NaNs
> >
> > I also tried BiCGStab as a solver (3 iterations with builtin Trilinos
> > solver, NaN residuum in case of the CBC.Block implementation).
> >
> > As far as I understand, both the Trilinos GMRES and the CBC.Block LGMRES
> > use the IFPACK ILU preconditioner. So does this mean the only difference
> > between the two libraries is the implementation of the iterative methods?
> >
> > I want to use the CBC.Block solver, because I want to implement a Schur
> > complement method where I have to solve another linear system per GMRES
> > step.
> >
> > Do you have any hints what I could do to make things work with CBC.Block?
> > I know this question is very general. I will gladly provide more
> > information if needed.
> >
> > Thanks in advance, Claas
> >
> > --
> > You received this question notification because you are a member of
> > cbc.block maintainers, which is an answer contact for CBC.Block.
> >
> > --
> > Post to : <email address hidden>
> > More help : https://help.launchpad.net/ListHelp
> >
>
> --
> You received this question notification because you are a member of
> cbc.block maintainers, which is an answer contact for CBC.Block.
>
> --
> Post to : <email address hidden>
>

 Revision history for this message Claas Abert (cabert) said on 2013-01-28: #3

Thanks Joachim Haga and Kent-Andre Mardal for the quick response.
I'm using the dolfin-dev rev 7250 (from the beginning of 2013) and the most recent dev-version of cbc.block (rev 121).
Here is a (not very minimal) example showing the problem.

The problem itself is a saddle point problem and I know that I should probably deal with this block structure properly using CBC.Block. However this is just the starting point for a bigger block system and I wanted to keep things simple in the beginning...

from dolfin import *
from block import *
from block.iterative import *
from block.algebraic.trilinos import *

# Set backend to Epetra
parameters['linear_algebra_backend'] = 'Epetra'

# Mesh and Function Spaces
mesh = BoxMesh(-50, -50, -5, 50, 50, 5, 10, 10, 1)
VV = VectorFunctionSpace(mesh, "CG", 1)
VS = FunctionSpace(mesh, "CG", 1)
V = VV * VS

# Test and Trial Functions
(v, sigma) = TrialFunctions(V)
(w2, w3) = TestFunctions(V)

# Define some constants and coeffecients
arg = "sqrt((3.141592*(x[0]/1e1))*(3.141592*(x[0]/1e1)))"
m_expr = Expression(("cos(%s)" % arg, "sin(%s)" % arg, "0.0"))
m = interpolate(m_expr, VV)

f_ex = -5.71823816178e+12
alpha = 0.1
dt = 1e-12

# Define weak forms
a = alpha * dot(v, w2) * dx + dot(cross(m, v), w2) * dx
a += - 0.5 * dt * f_ex * Dx(v[i],j) * Dx(w2[i],j) * dx
a += sigma * inner(m, w2) * dx
a += inner(m, v) * w3 * dx

L = f_ex * Dx(m[i],j) * Dx(w2[i],j) * dx

# Assemble
A = assemble(a)
b = assemble(L)

# Solve with Trilinos solvers
x = Function(V)
solve(A, x.vector(), b, "gmres", "ilu")
solve(A, x.vector(), b, "bicgstab", "ilu")

# Solve with CBC.Block solvers
Ap = ILU(A)
Ainv = LGMRES(A, precond=Ap)
#Ainv = BiCGStab(A, precond=Ap)
x = Ainv * b

 Revision history for this message Joachim Haga (jobh) said on 2013-01-29: #4

I have looked briefly. The convergence of builtin Trilinos solvers is illusory; try

# Solve with Trilinos solvers
x = Function(V)
solve(A, x.vector(), b, "gmres", "ilu")
plot(x[0])
solve(A, x.vector(), b, "bicgstab", "ilu")
plot(x[0])
solve(A, x.vector(), b, "bicgstab", "ilu")
plot(x[0])
interactive()

--- you'll see three very different "solutions". I suspect that the problem is that the full matrix A is indefinite; the preconditioner should probably be positive. The built-in solver probably bails out when encountering division-by-zero (which generates NaN/Inf), while cbc.block doesn't check for this and hence throws an exception.

 Revision history for this message Joachim Haga (jobh) said on 2013-01-29: #5

Oops, that last solve() should be just

solve(A, x.vector(), b)

to use a direct solver.

Anyway, it looks like simple doesn't work here. Go blockwise ;)

 Revision history for this message Claas Abert (cabert) said on 2013-01-30: #6

Thanks a lot Joachim!
I'm sorry that I didn't check that before posting.
I will go for the blocks now :)