Skip to main content
`Ctrl`+`K`

# 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.

    try: range = xrange
    except: pass

    import sys, slepc4py
    slepc4py.init(sys.argv)

    from petsc4py import PETSc
    from slepc4py import SLEPc

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)

    PETSc.Sys.Print( "Lauchli singular value decomposition, (%d x %d) mu=%g\n" % (n+1,n,mu) )

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

    A = PETSc.Mat(); A.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(); S.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 = PETSc.Sys.Print

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

    svd_type = S.getType()
    Print( "Solution method: %s" % svd_type )

    its = S.getIterationNumber()
    Print( "Number of iterations of the method: %d" % its )

    nsv, ncv, mpd = S.getDimensions()
    Print( "Number of requested singular values: %d" % nsv )

    tol, maxit = S.getTolerances()
    Print( "Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit) )

    nconv = S.getConverged()
    Print( "Number of converged approximate singular triplets %d" % 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( "   %6f     %12g" % (sigma, error) )
      Print()