Subject: [FEniCS-users] Computing Eigenvalues for a symmetric matrix using slepc

From: Jehanzeb Hameed <[email-masked]>
Date: Wed, 21 Oct 2009 22:33:52 -0500

Hello,

I am trying to use slepc to compute eigenvalues of a symmetic matrix
using SLEPc eigensolver.  However I get the following error:

[0]PETSC ERROR: Requested method is only available for Hermitian problems!

I think what is happening is that PETSC matrix is set to non-symmetric
by default when it is downcast from a Matrix. Anyone has any ideas?
Is there a way to  set PETSc matrix's symmetric flag?

Thanks,
-Jehanzeb

---

From: Chong Luo <[email-masked]>
Date: Thu, 22 Oct 2009 07:43:00 -0700 (PDT)

I use the c++ version of ALGLIB (http://www.alglib.net/) to do the eigenvalue decomposition and SVD, QR decomposition etc.. It's very easy to use.

Best,
Chong

________________________________
From: Jehanzeb Hameed <[email-masked]>
To: [email-masked]
Sent: Wed, October 21, 2009 10:33:52 PM
Subject: [FEniCS-users] Computing Eigenvalues for a symmetric matrix using slepc

Hello,

Thanks,
-Jehanzeb
_______________________________________________
FEniCS-users mailing list
[email-masked]
http://fenics.org/mailman/listinfo/fenics-users

---

From: Marie Rognes <[email-masked]>
Date: Thu, 22 Oct 2009 16:44:26 +0200

> Is there a way to  set PETSc matrix's symmetric flag?

Not from the DOLFIN interface at the moment. This could probably be added.

Which solver do you want to use? Have you tried krylov-schur?

--
Marie

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 22 Oct 2009 15:48:12 +0100

> Is there a way to  set PETSc matrix's symmetric flag?

No changes are made to the PETSc object when down casting. Are you sure 
that your matrix is symmetric? One ways of applying boundary conditions 
in DOLFIN will destroy symmetry. In PETSc's simple interface it is not 
possible say that a matrix is symmetric, but there is a function that 
tests for symmetry which SLEPc is probably using.

Can you post a very short piece of code that reproduces the problem?

Garth

---

From: Jehanzeb Hameed <[email-masked]>
Date: Thu, 22 Oct 2009 10:25:09 -0500

I think the matrix is numerically symmetric. I verified this using
matlab (max difference of A and A^T is 1e-15). So the only loss of
symmetry is coming from assembly of ufl form, not from applying
boundary conditions.  Maybe this is the problem?

Here is the piece of code to reproduce the bug:

PETScMatrix& AA(A.down_cast<PETScMatrix>()); // A is symmetric
unsigned int n = 1;
SLEPcEigenSolver esolver;
esolver.parameters["spectrum"] = "smallest magnitude";
esolver.set("eigenvalue solver","lanczos");
esolver.solve(AA, n);

I finally managed to solve this using lanczos by directly using SLEPc.
The relevant line I had to set was,

EPSSetProblemType( eps, EPS_HEP );

Thanks,
-Jehanzeb

On Thu, Oct 22, 2009 at 9:48 AM, Garth N. Wells <[email-masked]> wrote:
> No changes are made to the PETSc object when down casting. Are you sure that
> your matrix is symmetric? One ways of applying boundary conditions in DOLFIN
> will destroy symmetry. In PETSc's simple interface it is not possible say
> that a matrix is symmetric, but there is a function that tests for symmetry
> which SLEPc is probably using.
>

---

From: Jehanzeb Hameed <[email-masked]>
Date: Thu, 22 Oct 2009 10:27:23 -0500

I want to use lanczos. I think I tried Kyrlov-shur a long time ago,
and couldnt get it to work.

-Jehanzeb

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 22 Oct 2009 16:32:01 +0100

> esolver.set("eigenvalue solver","lanczos");

But how is the matrix constructed?

Garth

---

From: Jehanzeb Hameed <[email-masked]>
Date: Thu, 22 Oct 2009 10:37:48 -0500

by

assemble_system(A, b, a, L, bc1);

On Thu, Oct 22, 2009 at 10:32 AM, Garth N. Wells <[email-masked]> wrote:
>>> that

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 22 Oct 2009 16:42:58 +0100

> assemble_system(A, b, a, L, bc1);

Doesn't help too much since you haven't indicated what the form a is.

There was a bug in assemble_system which was fixed in 0.94. I can't 
recall exactly whether or not it would impact symmetry.

To debug, you can access the PETSc object associated with the Matrix and 
use the PETSc function to test for symmetry.

Garth

---

From: Jehanzeb Hameed <[email-masked]>
Date: Thu, 22 Oct 2009 10:59:07 -0500

The form is quite complicated, I can still post it if you wish.

I did use PETSc function to check for symmetry, and gave it a
tolerance of 1e-10. It return TRUE. Using 0 for tolerance gives false.

On installing 0.9.4. I havent figured out how to do these upgrades.
Will apt-get work for doing from 0.9.3 to 0.9.4? If I try it, it
installs 0.9.3 along with PETSc 2.

The experiments I ran were using dorsal-9cdbff7b6a6b. I am not sure
which version of dolfin it installs, but it uses PETSc 3. What is your
recommended way of updating/installing dolfin?

Thanks,
-Jehanzeb

On Thu, Oct 22, 2009 at 10:42 AM, Garth N. Wells <[email-masked]> wrote:
> There was a bug in assemble_system which was fixed in 0.94. I can't recall
> exactly whether or not it would impact symmetry.
>
> To debug, you can access the PETSc object associated with the Matrix and use
> the PETSc function to test for symmetry.
>
> Garth
>

---

From: Marie Rognes <[email-masked]>
Date: Thu, 22 Oct 2009 18:21:29 +0200

The SLEPc solvers that are restricted to Hermittian problems
actually require the

to be set (as Jehanzeb noted). If not, they fail even if the matrix is 
symmetric.

--
Marie

---

From: Marie Rognes <[email-masked]>
Date: Thu, 22 Oct 2009 18:22:51 +0200

> and couldnt get it to work.

Ok.  (I asked the SLEPc developers about this at one point, and they 
recommended using
krylov-schur instead of lanczos.  You might want to try to use a 
shift-and-invert
spectral transform if krylov-schur gives zero pivot error.)

--
Marie

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 22 Oct 2009 17:26:42 +0100

> to be set (as Jehanzeb noted). If not, they fail even if the matrix is 

OK. I just registered a bug so that it doesn't get forgotten,

    https://bugs.launchpad.net/dolfin/+bug/458309

Garth

>> To debug, you can access the PETSc object associated with the Matrix 
>> and use the PETSc function to test for symmetry.
>>
>> Garth
>>
>>

---

From: Jehanzeb Hameed <[email-masked]>
Date: Thu, 22 Oct 2009 13:58:17 -0500

On a related note, is there a way to change the values in a matrix
after assembly? I tried using setrow and setitem, but I get a PETSc
error saying

[0]PETSC ERROR: Not for unassembled matrix!

-Jehanzeb

---

From: Anders Logg <[email-masked]>
Date: Thu, 22 Oct 2009 21:10:22 +0200

If you modify the entries of a PETScMatrix, you must do

  A.apply()

afterwards.

This corresponds to calling MatAssemblyBegin/MatAssemblyEnd in PETSc.

--
Anders

---

From: Murtazo Nazarov <[email-masked]>
Date: Fri, 23 Oct 2009 01:12:14 +0200

Hi,

I am trying to implement an artificial viscosity of form 
nu[i]*Dx(v[i],j)*Dx(u[i],j)*dx, where nu[i] is constant vector on a 
cell, in FFC 0.7.0. The UFL form of the expression is as the following:

cell = triangle
vector = VectorElement("Lagrange", cell, 1)
constant_vector = VectorElement("Discontinuous Lagrange", cell, 0)

v  = TestFunction(vector)
nu = Function(constant_vector)
u  = Function(vector)

L = nu[i]*Dx(v[i],j)*Dx(u[i],j)*dx

It gives:

...
ufl.log.UFLException: Trying to integrate expression of rank 0 with free 
indices (Index(0),).

However, nu[0]*Dx(v[i],j)*Dx(u[i],j)*dx works.

Thanks,
murtazo

---

From: Anders Logg <[email-masked]>
Date: Fri, 23 Oct 2009 01:18:08 +0200

> L = nu[i]*Dx(v[i],j)*Dx(u[i],j)*dx

Indices must be repeated exactly twice in summation expressions.

Are you sure you want to have a vector-valued viscosity?

It would be more natural to do something like

L = nu*inner(grad(v), grad(u))*dx

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 23 Oct 2009 00:26:37 +0100

> Are you sure you want to have a vector-valued viscosity?

Indeed, a vector-valued viscosity doesn't make sense. It must be a 
scalar or a second-order tensor.

Garth

>> ------------------------------------------------------------------------
>>

---

From: Murtazo Nazarov <[email-masked]>
Date: Fri, 23 Oct 2009 09:07:17 +0200

>> Are you sure you want to have a vector-valued viscosity?

Yes, I want to have a vector valued viscosity. For example for u[0] and 
u[1]. Here Nu[i] is an artificial viscosity which are residual based for 
each component.

>> L = nu*inner(grad(v), grad(u))*dx

I have tried it already.

/murtazo

---

From: Anders Logg <[email-masked]>
Date: Fri, 23 Oct 2009 10:36:53 +0200

On Fri, Oct 23, 2009 at 09:07:17AM +0200, Murtazo Nazarov wrote:
> Yes, I want to have a vector valued viscosity. For example for u[0]
> and u[1]. Here Nu[i] is an artificial viscosity which are residual
> based for each component.

It seems to me you really have a matrix-valued viscosity (which makes
sense) but it is diagonal.

Try something like this:

NU = as_matrix([[nu[0], 0, 0], [0, nu[1], 0], [0, 0, nu[2]]])

L = inner(grad(v), NU*grad(u))*dx

Note that you may need to use grad(v).T and grad(u).T since grad() is
defined as du_j/du_i in UFL. (I'm getting more and more confused with
this and will propose we change this at some point.)

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Fri, 23 Oct 2009 09:49:12 +0100

> this and will propose we change this at some point.)

It must be the symmetric gradient,

     L = inner(grad(v), NU*sym_grad(u))*dx

in which case the definition of grad won't matter.

Garth

---

From: Anders Logg <[email-masked]>
Date: Fri, 23 Oct 2009 11:46:43 +0200

Hi Murtazo,

Your emails don't seem to make it to the list since you are posting
from [email-masked] which is not registered on the list. The
address you have registered is [email-masked].

--
Anders

On Fri, Oct 23, 2009 at 11:08:59AM +0200, Murtazo Nazarov wrote:
> I think the diagonal matrix makes sense. Thank you, I will try this.
>
> /murtazo
>
>>>
>>>
>>>

---

From: Murtazo Nazarov <[email-masked]>
Date: Thu, 29 Oct 2009 19:01:23 +0100

>>> L = inner(grad(v), NU*grad(u))*dx

It worked.

>>      L = inner(grad(v), NU*sym_grad(u))*dx

It did not work:
NameError: name 'sym_grad' is not defined

Thanks!
Murtazo

---

From: "Garth N. Wells" <[email-masked]>
Date: Thu, 29 Oct 2009 19:05:08 +0000

>>>> L = inner(grad(v), NU*grad(u))*dx

This approach is no physically correct. The viscous stress cannot depend 
on the plain gradient.

> NameError: name 'sym_grad' is not defined

You need to provide the sym_grad function yourself.

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 02 Nov 2009 09:11:21 +0100

> You need to provide the sym_grad function yourself.

I think sym(grad(u)) should do the trick.

Kristian

---

From: Murtazo Nazarov <[email-masked]>
Date: Mon, 02 Nov 2009 12:03:32 +0100

> I think sym(grad(u)) should do the trick.

Great!

murtazo