Function cusolverDnZgesvdj 

pub unsafe extern "C" fn cusolverDnZgesvdj(
    handle: cusolverDnHandle_t,
    jobz: cusolverEigMode_t,
    econ: c_int,
    m: c_int,
    n: c_int,
    A: *mut cuDoubleComplex,
    lda: c_int,
    S: *mut f64,
    U: *mut cuDoubleComplex,
    ldu: c_int,
    V: *mut cuDoubleComplex,
    ldv: c_int,
    work: *mut cuDoubleComplex,
    lwork: c_int,
    info: *mut c_int,
    params: gesvdjInfo_t,
) -> cusolverStatus_t

The helper functions below can calculate the sizes needed for pre-allocated buffer.

The S and D data types are real valued single and double precision, respectively.

The C and Z data types are complex valued single and double precision, respectively.

This function computes the singular value decomposition (SVD) of an $m \times n$ matrix A and corresponding the left and/or right singular vectors. The SVD is written: $$ A = U\*\Sigma\*V^{H} $$

where $\Sigma$ is an $m \times n$ matrix which is zero except for its min(m,n) diagonal elements, U is an $m \times m$ unitary matrix, and V is an $n \times n$ unitary matrix. The diagonal elements of $\Sigma$ are the singular values of A; they are real and non-negative, and are returned in descending order. The first min(m,n) columns of U and V are the left and right singular vectors of A.

gesvdj has the same functionality as gesvd. The difference is that gesvd uses QR algorithm and gesvdj uses Jacobi method. The parallelism of Jacobi method gives GPU better performance on small and medium size matrices. Moreover the user can configure gesvdj to perform approximation up to certain accuracy.

gesvdj iteratively generates a sequence of unitary matrices to transform matrix A to the following form: $$ U^{H}\*A\*V = S + E $$

where S is diagonal and diagonal of E is zero.

During the iterations, the Frobenius norm of E decreases monotonically. As E goes down to zero, S is the set of singular values. In practice, Jacobi method stops if: $$ {\|E\|}{F}\leq\operatorname{eps}\*{\|A\|}{F} $$

where eps is given tolerance. Note that if the real residual norm: $$ {\|{S} - {U}^{H}\{A}\{V}\|}_{F} $$

is computed, it will differ from ${\|{E}\|}_{F}$ up to roundoff errors of order $N = max(m, n)$, to still have the standard SVD accuracy expectation: $$ \frac{\|S - U^{H} \* A \* V\|_F}{O(N) \* \|A\|_F} \leq \frac{\|E\|_F}{\|A\|_F} \leq \operatorname{eps} $$

$O(N)$ is typically $N$, but the constant depends on the number of sweeps, which gives an upper roundoff error bound of $sweeps \* N$.

gesvdj has two parameters to control the accuracy. First parameter is tolerance (eps). The default value is machine accuracy but The user can use function cusolverDnXgesvdjSetTolerance to set a priori tolerance. The second parameter is maximum number of sweeps which controls number of iterations of Jacobi method. The default value is 100 but the user can use function cusolverDnXgesvdjSetMaxSweeps to set a proper bound. The experiments show 15 sweeps are good enough to converge to machine accuracy. gesvdj stops either tolerance is met or maximum number of sweeps is met.

Jacobi method has quadratic convergence, so the accuracy is not proportional to number of sweeps. To guarantee certain accuracy, the user should configure tolerance only.

The user has to provide working space which is pointed by input parameter work. The input parameter lwork is the size of the working space, and it is returned by gesvdj_bufferSize(). Please note that the size in bytes of the working space is equal to sizeof(<type>) * lwork.

If output parameter info = -i (less than zero), the i-th parameter is wrong (not counting handle). If info = min(m,n)+1, gesvdj does not converge under given tolerance and maximum sweeps.

If the user sets an improper tolerance, gesvdj may not converge. For example, tolerance should not be smaller than machine accuracy.

Please visit cuSOLVER Library Samples - gesvdj for a code example.

Remark 1: gesvdj supports any combination of m and n.

Remark 2: the routine returns V, not $V^{H}$. This is different from gesvd.