Skip to main content

cusolverDnXpolar

singe_cusolver_sys

Function cusolverDnXpolar 

Source 
pub unsafe extern "C" fn cusolverDnXpolar(
    handle: cusolverDnHandle_t,
    params: cusolverDnParams_t,
    uplo: cublasFillMode_t,
    M: i64,
    N: i64,
    dataTypeA: cudaDataType,
    A: *mut c_void,
    lda: i64,
    dataTypeH: cudaDataType,
    H: *mut c_void,
    ldh: i64,
    computeType: cudaDataType,
    bufferOnDevice: *mut c_void,
    workspaceInBytesOnDevice: size_t,
    bufferOnHost: *mut c_void,
    workspaceInBytesOnHost: size_t,
    d_res_nrm: *mut f64,
    d_A_nrmF: *mut f64,
    d_rcond: *mut f64,
    d_info: *mut c_int,
) -> cusolverStatus_t
Expand description

The following helper function calculates the required workspace sizes for cusolverDnXpolar.

The following routine computes the polar decomposition.

This function computes the polar decomposition of an $m \times n$ matrix A (where $m \geq n$). The polar decomposition is written: $$ A = U_p \cdot H $$

where $U_p$ is an $m \times n$ matrix with orthonormal columns and $H$ is an $n \times n$ Hermitian positive semidefinite matrix.

cusolverDnXpolar uses the QDWH (QR-based Dynamically Weighted Halley) iteration [14][16] to compute $U_p$. The QDWH iteration converges cubically and requires at most ~6 iterations. For ill-conditioned or rank-deficient matrices A, the routine adds a small perturbation $\xi$ to ensure the polar decomposition converges correctly. The output parameter d_res_nrm reports the magnitude of the residual $\|A - U_p \cdot H\|_2$.

The user has to provide device and host working spaces which are pointed by input parameters bufferOnDevice and bufferOnHost. The input parameters workspaceInBytesOnDevice (and workspaceInBytesOnHost) is size in bytes of the device (and host) working space, and it is returned by cusolverDnXpolar_bufferSize.

If output parameter info = -i (less than zero), the i-th parameter is wrong (not counting handle).

Currently, cusolverDnXpolar supports only the default algorithm.

Algorithms supported by cusolverDnXpolar
cusolverAlgMode_t::CUSOLVER_ALG_0 or NULLDefault algorithm.

Remark 1: If H is NULL, computation of H is skipped.

Remark 2: uplo = CUBLAS_FILL_MODE_UPPER treats A as upper triangular, saving one QR decomposition. cublasFillMode_t::CUBLAS_FILL_MODE_LOWER is not supported.

Remark 3: d_res_nrm serves dual purpose. On input, if *d_res_nrm >= 0, it overrides the perturbation $\xi$ (skipping automatic search). On output, it reports the residual norm $\|A - U_p \cdot H\|_2$.

List of input arguments for cusolverDnXpolar_bufferSize and cusolverDnXpolar:

The generic API has three different types, dataTypeA is data type of the matrix A, dataTypeH is data type of the matrix H, and computeType is compute type of the operation. cusolverDnXpolar only supports the following four combinations:

Valid combination of data type and compute type

DataTypeADataTypeHComputeTypeMeaning
CUDA_R_32FCUDA_R_32FCUDA_R_32FSPOLAR
CUDA_R_64FCUDA_R_64FCUDA_R_64FDPOLAR
CUDA_C_32FCUDA_C_32FCUDA_C_32FCPOLAR
CUDA_C_64FCUDA_C_64FCUDA_C_64FZPOLAR

§Parameters

  • handle: Handle to the cuSolverDN library context.
  • params: Structure with information collected by cusolverDnSetAdvOptions.
  • uplo: Specifies the structure of matrix A. uplo = cublasFillMode_t::CUBLAS_FILL_MODE_UPPER: A is upper triangular, saving one QR decomposition. uplo = cublasFillMode_t::CUBLAS_FILL_MODE_FULL: A is a general matrix. cublasFillMode_t::CUBLAS_FILL_MODE_LOWER is not supported.
  • dataTypeA: Data type of array A.
  • A: Array of dimension lda * n with lda is not less than max(1,m). On exit, A is overwritten with the unitary polar factor $U_p$ with orthonormal columns.
  • lda: Leading dimension of two-dimensional array used to store matrix A. lda >= max(1,m).
  • dataTypeH: Data type of array H. Must equal dataTypeA.
  • H: Array of dimension ldh * n with ldh is not less than max(1,n). $n \times n$ Hermitian positive semidefinite factor. If H is NULL, computation of H is skipped.
  • ldh: Leading dimension of two-dimensional array used to store matrix H. ldh >= max(1,n) when H is not NULL.
  • computeType: Data type of computation. Must equal dataTypeA.
  • bufferOnDevice: Device workspace. Array of type void of size workspaceInBytesOnDevice bytes.
  • workspaceInBytesOnDevice: Size in bytes of bufferOnDevice, returned by cusolverDnXpolar_bufferSize.
  • bufferOnHost: Host workspace. Array of type void of size workspaceInBytesOnHost bytes.
  • workspaceInBytesOnHost: Size in bytes of bufferOnHost, returned by cusolverDnXpolar_bufferSize.
  • d_res_nrm: On input, if *d_res_nrm >= 0, it overrides the perturbation $\xi$ (skipping automatic search). On output, reports $\Vert A - U_p \cdot H \Vert_2$.
  • d_A_nrmF: $\Vert A \Vert_F$ (Frobenius norm of A). If NULL, the computation is skipped.
  • d_rcond: Reciprocal condition number $1/\text{cond}(A, \text{fro})$. If NULL, the computation is skipped. Undefined if a user-supplied $\xi > 0$ was provided via d_res_nrm.
  • d_info: If info = 0, the operation is successful. If info = -i, the i-th parameter is wrong (not counting handle). If info = 1, the algorithm failed to find rcond > 1.e-15 by adding perturbation. If info = 10, the algorithm failed to converge within 10 iterations.

§Return value