Skip to main content

dgesvdj

Function dgesvdj 

Source
pub fn dgesvdj(
    ctx: &Context,
    jobz: EigenMode,
    econ: bool,
    m: usize,
    n: usize,
    a: MatrixMut<'_, f64>,
    s: &mut DeviceMemory<f64>,
    u: Option<MatrixMut<'_, f64>>,
    v: Option<MatrixMut<'_, f64>>,
    workspace: &mut DeviceMemory<f64>,
    dev_info: &mut DeviceMemory<i32>,
    params: &GesvdjInfo,
) -> Result<()>
Expand description

Use the matching buffer-size helper to calculate the sizes needed for pre-allocated workspace.

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.

Computes the singular value decomposition (SVD) of an $m \times n$ matrix A and the corresponding left and/or right singular vectors. The SVD is written as $A = U \Sigma V^{H}$, where Σ 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 Σ 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 computes the same decomposition as gesvd. The difference is that gesvd uses a QR algorithm and gesvdj uses the Jacobi method. The Jacobi method gives GPUs better parallelism on small and medium-size matrices. Callers can configure gesvdj to target a chosen accuracy.

gesvdj iteratively generates a sequence of unitary matrices that transform A toward $A = U(S + E)V^{H}$, where S is diagonal and the 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, the Jacobi method stops when the off-diagonal residual is below the configured tolerance eps. If the real residual norm is computed, it differs from ${\|{E}\|}_{F}$ by roundoff errors of order $N = max(m, n)$ while still meeting the standard SVD accuracy expectation.

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

gesvdj has two parameters to control the accuracy. The first parameter is the tolerance (eps). The default value is machine accuracy, but GesvdjInfo::set_tolerance can set an a priori tolerance. The maximum-sweep parameter is the maximum number of sweeps, which controls the number of Jacobi iterations. The default value is 100, but GesvdjInfo::set_max_sweeps can set a different bound. Experiments show that 15 sweeps are enough to converge to machine accuracy. gesvdj stops when either the tolerance or the maximum number of sweeps is reached.

The Jacobi method has quadratic convergence, so the accuracy is not proportional to the number of sweeps. To guarantee a target accuracy, configure only the tolerance.

Provide workspace through workspace. Use the corresponding *_buffer_size helper to query the required workspace length. The workspace size in bytes is size_of::<T>() * lwork.

If the reported info value is -i, the ith parameter is invalid. If info == min(m, n) + 1, gesvdj did not converge within the given tolerance and maximum sweep count.

If the tolerance is too small, gesvdj may not converge. Use a tolerance no smaller than machine accuracy.

  • gesvdj supports any combination of m and n.

  • Returns V, not $V^{H}$. This is different from gesvd.

§Errors

Returns an error if cuSOLVER has not been initialized, if the matrix dimensions, leading dimensions, or vector-computation mode are invalid, or if cuSOLVER reports an internal failure.