# 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
# <../_static/ex4.py>`__.

# 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()