SLEPc SVD solvers with dolfin

Question

I would like to use SVD factorization from SLEPc in rectangular sparse matrix.

The script below illustrates this type of matrices.

mesh0 = UnitSquareMesh(64, 64)

V0 = FunctionSpace(mesh0, "CG", 2)
V1 = FunctionSpace(mesh0, "CG", 1)

W = MixedFunctionSpace([V0, V1])

V0_dofs = W.sub(0).collapse(mesh0)[1].values()

(u0, u1) = TrialFunctions(W)
(v0, v1) = TestFunctions(W)

f = Expression("-10exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
g = Expression("sin(5x[0])")

a = inner(grad(u0 + u1), grad(v0)) * dx
L = f * (v0) * dx + g * (v0) * ds

A, b = assemble_system(A_form = a, b_form = L)

A0 = A[V0_dofs, :]
b0 = b[V0_dofs]

Can I access with dolfin this feature of SLEPc for computing the partial singular value decomposition (SVD) of a large sparse rectangular matrix?

( Are solvers such as Lanczos and thick-restart Lanczos bidiagonalization available? )

Thanks in advance!

Comments

Hi, you mentioned rectangular matrices but the matrix you assemble is square. Is the matrix
correct?

---

Thanks for the observation! I try to fix the imprecision. I would like to use the SVD solver for rectangular matrices or rank deficient square ones.

Answer 1

Hi, you could perhaps use slepc4py

import sys, slepc4py
slepc4py.init(sys.argv)
from slepc4py import SLEPc
from dolfin import *

mesh = UnitSquareMesh(20, 20)
CG = FunctionSpace(mesh, 'CG', 1)
CR = FunctionSpace(mesh, 'CR', 1)
cg = TrialFunction(CR)
cr = TestFunction(CG)

# A is a matrix from CR to CG, m=CG.dim() x n=CR.dim()
a = inner(cg, cr)*dx
A = PETScMatrix()
A = assemble(a, tensor=A)
A_mat = A.mat()

S = SLEPc.SVD()
S.create()
S.setOperator(A_mat)
S.setType(S.Type.LANCZOS)
S.solve()

nconv = S.getConverged()
if nconv > 0:
    v, u = A_mat.getVecs()
    for i in range(nconv):
        sigma = S.getSingularTriplet(i, u, v)
        print sigma

Comments

MiroK,

many thanks for your reply!

I was trying to do this directly as in the example

http://fenicsproject.org/documentation/dolfin/1.3.0/python/demo/documented/eigenvalue/python/documentation.html

u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx

A = PETScMatrix()
assemble(a, tensor=A)

eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()

r, c, rx, cx = eigensolver.get_eigenpair(0)

print "Largest eigenvalue: ", r

but I can't find some class like SLEPcSVDSolver.

---

Hi, I believe only eigenvalue solver from SLEPc is wrapped in DOLFIN and not the SVD routines. That's why I suggested slepc4py. The installation is much like petsc4py. For petsc4py there is a demo.

---

Thanks! I installed the petsc4py, slepc4py and tao4py through pip
on Ubuntu 12.04 LST.

Do you know some script demo of tao4py with FEniCS?