# Exercise 4: Singular Value Decomposition In this exercise we turn our attention to the Singular Value Decomposition (SVD). Remember that in positive-definite real symmetric (or complex Hermitian) matrices, singular values coincide with eigenvalues, but in general this is not the case. The SVD is defined for any matrix, even rectangular ones. Singular values are always non-negative real values. This example works also by reading a matrix from a file. In particular, the matrix to be used is related to a 2D reaction-diffusion model. More details about this problem can be found at [Matrix Market](https://math.nist.gov/MatrixMarket/data/NEP/brussel/brussel.html). ## Compiling Copy the file {{'[ex14.c](https://slepc.upv.es/{}/src/svd/tutorials/ex14.c.html)'.format(branch)}} to your directory and add these lines to the makefile ```{code} make ex14: ex14.o -${CLINKER} -o ex14 ex14.o ${SLEPC_SVD_LIB} ${RM} ex14.o ``` Build the executable with the command ```{code} console $ make ex14 ``` ## Running the Program In order to run the program, type the following: ```{code} console $ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc ``` You should get an output similar to this: ```{code} Singular value problem stored in file. Reading REAL matrix from a binary file... Number of iterations of the method: 3 Solution method: cross Number of requested singular values: 1 Stopping condition: tol=1e-08, maxit=100 SVD solve converged (1 singular triplet) due to CONVERGED_TOL; iterations 3 ---------------------- -------------------- sigma relative error ---------------------- -------------------- 35.007519 1.40809e-10 ---------------------- -------------------- ``` ## Source Code Details The way in which the `SVD` object works is very similar to that of `EPS`. However, some important differences exist. Examine the source code of the example program and pay attention to the differences with respect to `EPS`. After loading the matrix, the problem is solved by the following sequence of function calls: * `SVDCreate(MPI_Comm comm,SVD *svd);` * `SVDSetOperators`(`SVD` `svd`,{external:doc}`Mat` `A`,{external:doc}`Mat` `B`);` * `SVDSetFromOptions(SVD svd);` * `SVDSolve(SVD svd);` * `SVDGetConverged`(`SVD` `svd`, {external:doc}`PetscInt` `*nconv`); * `SVDGetSingularTriplet`(`SVD` `svd`,{external:doc}`PetscInt` `i`,{external:doc}`PetscReal` `*sigma`,{external:doc}`Vec` `u`,{external:doc}`Vec` `v`); * `SVDDestroy(SVD svd)`; First, the singular value solver (`SVD`) context is created and the matrix associated with the problem is specified. Then various options are set for customized solution. After that, the program solves the problem, retrieves the solution, and finally destroys the `SVD` context. Note that the singular value, `sigma`, is defined as a {external:doc}`PetscReal`, and that the singular vectors are simple {external:doc}`Vec`'s. ## SVD Options Most of the options available in the `EPS` object have their equivalent in `SVD`. A full list of command-line options can be obtained by running the example with the option `-help`. To show information about the `SVD` solver, add the `-svd_view` option: ```{code} console $ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc -svd_view ``` :::{note} All the command-line options related to the `SVD` object have the `-svd_` prefix. ::: Your output will include something like this: ```{code} SVD Object: 1 MPI processes type: cross implicit matrix EPS Object: (svd_cross_) 1 MPI processes type: krylovschur 50% of basis vectors kept after restart using the locking variant problem type: symmetric eigenvalue problem selected portion of the spectrum: largest real parts number of eigenvalues (nev): 1 number of column vectors (ncv): 16 maximum dimension of projected problem (mpd): 16 maximum number of iterations: 100 tolerance: 1e-09 convergence test: relative to the eigenvalue BV Object: (svd_cross_) 1 MPI processes type: svec 17 columns of global length 200 vector orthogonalization method: classical Gram-Schmidt orthogonalization refinement: if needed (eta: 0.7071) block orthogonalization method: GS doing matmult as a single matrix-matrix product DS Object: (svd_cross_) 1 MPI processes type: hep solving the problem with: Implicit QR method (_steqr) ST Object: (svd_cross_) 1 MPI processes type: shift shift: 0. number of matrices: 1 using a shell matrix transpose mode: explicit selected portion of the spectrum: largest number of singular values (nsv): 1 number of column vectors (ncv): 16 maximum dimension of projected problem (mpd): 16 maximum number of iterations: 100 tolerance: 1e-08 convergence test: relative to the singular value ``` The output shows all the options that are susceptible of being changed, either from the command line or from the source code of the program: the method, the portion of the spectrum (largest or smallest singular values), the number of singular values (`nsv`), etc. Try to change some of the values, for instance: ```{code} console $ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc -svd_nsv 10 -svd_ncv 40 -svd_smallest ``` The "transpose mode" refers to whether the transpose of matrix A is being built explicitly or not (see `SVDSetImplicitTranspose` for an explanation). Note that in the sample output above, the `SVD` object contains an `EPS` object. This only happens in some `SVD` solver types, as detailed below. ## Changing the Singular Value Solver SLEPc provides several solvers for computing the SVD, which can be selected in the source code with the function `SVDSetType`, or at run time: ```{code} console $ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc -svd_type trlanczos ``` The following table shows the list of SVD solvers available in SLEPc. Solver | Command-line Name | Parameter --- | --- | --- Cross Product | cross | SVDCROSS Cyclic Matrix | cyclic | SVDCYCLIC Lanczos with Explicit Restart | lanczos | SVDLANCZOS Lanczos with Thick Restart | trlanczos | SVDTRLANCZOS LAPACK | lapack | SVDLAPACK :::{note} The LAPACK solver is not really a full-featured singular value solver but simply an interface to some LAPACK routines. These routines operate sequentially in dense mode and therefore are suitable only for small size problems. This solver should be used only for debugging purposes. ::: :::{note} The default solver is `cross`. ::: The first two solvers, `cross` and `cyclic`, are not real methods implemented in the `SVD` module, but are two convenient ways of solving the SVD problem by making use of the eigensolvers available in the `EPS` module. In those two cases, the `SVD` object manages an `EPS` object internally, whose parameters can be set as desired (typically only the method). For example: ```{code} console $ ./ex14 -file $SLEPC_DIR/share/slepc/datafiles/matrices/rdb200.petsc -svd_type cyclic -svd_cyclic_eps_type lanczos -svd_cyclic_eps_lanczos_reorthog local ```