ex4.py: Singular value decomposition of the Lauchli matrix#

This example illustrates the use of the SVD solver in slepc4py. It computes singular values and vectors of the Lauchli matrix, whose condition number depends on a parameter mu.

The full source code for this demo can be downloaded here.

Initialization is similar to previous examples.

import sys
import slepc4py

slepc4py.init(sys.argv)

from petsc4py import PETSc
from slepc4py import SLEPc

Print = PETSc.Sys.Print

This example takes two command-line arguments, the matrix size n and the mu parameter.

opts = PETSc.Options()
n = opts.getInt('n', 30)
mu = opts.getReal('mu', 1e-6)

Print(f'Lauchli singular value decomposition, ({n + 1} x {n}) mu={mu}\n')

Create the matrix and fill its nonzero entries. Every MPI process will insert its locally owned part only.

A = PETSc.Mat().create()
A.setSizes([n + 1, n])
A.setFromOptions()

rstart, rend = A.getOwnershipRange()

for i in range(rstart, rend):
    if i == 0:
        for j in range(n):
            A[0, j] = 1.0
    else:
        A[i, i - 1] = mu

A.assemble()

The singular value solver is similar to the eigensolver used in previous examples. In this case, we select the thick-restart Lanczos bidiagonalization method.

S = SLEPc.SVD().create()

S.setOperator(A)
S.setType(S.Type.TRLANCZOS)
S.setFromOptions()

S.solve()

After solve, we print some informative data and extract the computed solution, showing the list of singular values and the corresponding residual errors.

Print('******************************')
Print('*** SLEPc Solution Results ***')
Print('******************************\n')

svd_type = S.getType()
Print(f'Solution method: {svd_type}')

its = S.getIterationNumber()
Print(f'Number of iterations of the method: {its}')

nsv, _ncv, _mpd = S.getDimensions()
Print(f'Number of requested singular values: {nsv}')

tol, maxit = S.getTolerances()
Print(f'Stopping condition: tol={tol:.4g}, maxit={maxit}')

nconv = S.getConverged()
Print(f'Number of converged approximate singular triplets {nconv}')

if nconv > 0:
    v, u = A.createVecs()
    Print()
    Print('    sigma       residual norm ')
    Print('-------------  ---------------')
    for i in range(nconv):
        sigma = S.getSingularTriplet(i, u, v)
        error = S.computeError(i)
        Print(f'   {sigma:6f}     {error:12g}')
    Print()

This Page