This is a pointer type to an opaque cuSolverDN context, which the user must initialize by calling prior to calling cusolverDnCreate() any other library function. The handle created and returned by cusolverDnCreate() must be passed to every cuSolverDN function.

The type indicates which part (lower or upper) of the dense matrix was filled and consequently should be used by the function. Its values correspond to Fortran characters ‘L’ or ‘l’ (lower) and ‘U’ or ‘u’ (upper) that are often used as parameters to legacy BLAS implementations.

The cublasOperation_t type indicates which operation needs to be performed with the dense matrix. Its values correspond to Fortran characters ‘N’ or ‘n’ (non-transpose), ‘T’ or ‘t’ (transpose) and ‘C’ or ‘c’ (conjugate transpose) that are often used as parameters to legacy BLAS implementations.

This is the same as cusolverStatus_t in the sparse LAPACK section.

This is a pointer type to an opaque cuSolverSP context, which the user must initialize by calling prior to calling cusolverSpCreate() any other library function. The handle created and returned by cusolverSpCreate() must be passed to every cuSolverSP function.

We have chosen to keep the same structure as exists in cuSparse to describe the shape and properties of a matrix. This enables calls to either cuSparse or cuSolver using the same matrix description.

typedef struct {
    cusparseMatrixType_t MatrixType;
    cusparseFillMode_t FillMode;
    cusparseDiagType_t DiagType;
    cusparseIndexBase_t IndexBase;
} cusparseMatDescr_t;

Please read documenation of CUSPARSE Library to understand each field of cusparseMatDescr_t.

This is a status type returned by the library functions and it can have the following values.

The cusolverRfHandle_t is a pointer to an opaque data structure that contains the cuSolverRF library handle. The user must initialize the handle by calling cusolverRfCreate() prior to any other cuSolverRF library calls. The handle is passed to all other cuSolverRF library calls.

The cusolverRfMatrixFormat_t is an enum that indicates the input/output matrix format assumed by the cusolverRfSetup(), cusolverRfSetupHost(), cusolverRfResetValues(), cusolveRfExtractBundledFactorsHost() and cusolverRfExtractSplitFactorsHost() routines.

The cusolverRfNumericBoostReport_t is an enum that indicates whether numeric boosting (of the pivot) was used during the cusolverRfRefactor() and cusolverRfSolve() routines. The numeric boosting is disabled by default.

The cusolverRfResetValuesFastMode_t is an enum that indicates the mode used for the cusolverRfResetValues() routine. The fast mode requires extra memory and is recommended only if very fast calls to cusolverRfResetValues() are needed.

The cusolverRfFactorization_t is an enum that indicates which (internal) algorithm is used for refactorization in the cusolverRfRefactor() routine.

The cusolverRfTriangularSolve_t is an enum that indicates which (internal) algorithm is used for triangular solve in the cusolverRfSolve() routine.

The cusolverRfUnitDiagonal_t is an enum that indicates whether and where the unit diagonal is stored in the input/output triangular factors in the cusolverRfSetup(), cusolverRfSetupHost() and cusolverRfExtractSplitFactorsHost() routines.

The cusolverStatus_t is an enum that indicates success or failure of the cuSolverRF library call. It is returned by all the cuSolver library routines, and it uses the same enumerated values as the sparse and dense Lapack routines.

cusolverStatus_t 
cusolverDnCreate(cusolverDnHandle_t *handle);

This function initializes the cuSolverDN library and creates a handle on the cuSolverDN context. It must be called before any other cuSolverDN API function is invoked. It allocates hardware resources necessary for accessing the GPU.

cusolverStatus_t 
cusolverDnDestroy(cusolverDnHandle_t handle);

This function releases CPU-side resources used by the cuSolverDN library.

cusolverStatus_t 
cusolverDnSetStream(cusolverDnHandle_t handle, cudaStream_t streamId)

This function sets the stream to be used by the cuSolverDN library to execute its routines.

cusolverStatus_t 
cusolverDnGetStream(cusolverDnHandle_t handle, cudaStream_t *streamId)

This function sets the stream to be used by the cuSolverDN library to execute its routines.

cusolverStatus_t
cusolverDnSpotrf_bufferSize(cusolverDnHandle_t handle,
                 cublasFillMode_t uplo,
                 int n,
                 float *A,
                 int lda,
                 int *Lwork );

cusolverStatus_t
cusolverDnDpotrf_bufferSize(cusolveDnHandle_t handle,
                 cublasFillMode_t uplo,
                 int n,
                 double *A,
                 int lda,
                 int *Lwork );

cusolverStatus_t 
cusolverDnCpotrf_bufferSize(cusolverDnHandle_t handle,
                 cublasFillMode_t uplo,
                 int n,
                 cuComplex *A,
                 int lda,
                 int *Lwork );

cusolverStatus_t 
cusolverDnZpotrf_bufferSize(cusolverDnHandle_t handle,
                 cublasFillMode_t uplo,
                 int n,
                 cuDoubleComplex *A,
                 int lda,
                 int *Lwork);

cusolverStatus_t 
cusolverDnSpotrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           float *A,
           int lda,
           float *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnDpotrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           double *A,
           int lda,
           double *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnCpotrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           cuComplex *A,
           int lda,
           cuComplex *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnZpotrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           cuDoubleComplex *A,
           int lda,
           cuDoubleComplex *Workspace,
           int Lwork,
           int *devInfo );

This function computes the Cholesky factorization of a Hermitian positive-definite matrix.

A is a n×n Hermitian matrix, only lower or upper part is meaningful. The input parameter uplo indicates which part of the matrix is used. The function would leave other part untouched.

If input parameter uplo is CUBLAS_FILL_MODE_LOWER, only lower triangular part of A is processed, and replaced by lower triangular Cholesky factor L.

If input parameter uplo is CUSBLAS_FILL_MODE_UPPER, only upper triangular part of A is processed, and replaced by upper triangular Cholesky factor U.

The user has to provide working space which is pointed by input parameter Workspace. The input parameter Lwork is size of the working space, and it is returned by potrf_bufferSize().

If Cholesky factorization failed, i.e. some leading minor of A is not positive definite, or equivalently some diagonal elements of L or U is not a real number. The output parameter devInfo would indicate smallest leading minor of A which is not positive definite.

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

cusolverStatus_t 
cusolverDnSpotrs(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           int nrhs,
           const float *A,
           int lda,
           float *B,
           int ldb,
           int *devInfo);

cusolverStatus_t 
cusolverDnDpotrs(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           int nrhs,
           const double *A,
           int lda,
           double *B,
           int ldb,
           int *devInfo);

cusolverStatus_t 
cusolverDnCpotrs(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           int nrhs,
           const cuComplex *A,
           int lda,
           cuComplex *B,
           int ldb,
           int *devInfo);

cusolverStatus_t 
cusolverDnZpotrs(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           int nrhs,
           const cuDoubleComplex *A,
           int lda,
           cuDoubleComplex *B,
           int ldb,
           int *devInfo);

This function solves a system of linear equations

where A is a n×n Hermitian matrix, only lower or upper part is meaningful. The input parameter uplo indicates which part of the matrix is used. The function would leave other part untouched.

The user has to call potrf first to factorize matrix A. If input parameter uplo is CUBLAS_FILL_MODE_LOWER, A is lower triangular Cholesky factor L correspoding to $A=L*L^H$ . If input parameter uplo is CUSBLAS_FILL_MODE_UPPER, A is upper triangular Cholesky factor U corresponding to $A=U*U^H$ .

The operation is in-place, i.e. matrix X overwrites matrix B with the same leading dimension ldb.

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

cusolverStatus_t 
cusolverDnSgetrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      float *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnDgetrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      double *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnCgetrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      cuComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnZgetrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      cuDoubleComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnSgetrf(cusolverDnHandle_t handle,
           int m,
           int n,
           float *A,
           int lda,
           float *Workspace,
           int *devIpiv,
           int *devInfo );

cusolverStatus_t 
cusolverDnDgetrf(cusolverDnHandle_t handle,
           int m,
           int n,
           double *A,
           int lda,
           double *Workspace,
           int *devIpiv,
           int *devInfo );

cusolverStatus_t 
cusolverDnCgetrf(cusolverDnHandle_t handle,
           int m,
           int n,
           cuComplex *A,
           int lda,
           cuComplex *Workspace,
           int *devIpiv,
           int *devInfo );

cusolverStatus_t 
cusolverDnZgetrf(cusolverDnHandle_t handle,
           int m,
           int n,
           cuDoubleComplex *A,
           int lda,
           cuDoubleComplex *Workspace,
           int *devIpiv,
           int *devInfo );

This function computes the LU factorization of a m×n matrix

where A is a m×n matrix, P is a permutation matrix, L is a lower triangular matrix with unit diagonal, and U is an upper triangular matrix.

The user has to provide working space which is pointed by input parameter Workspace. The input parameter Lwork is size of the working space, and it is returned by getrf_bufferSize().

If LU factorization failed, i.e. matrix A (U) is singular, The output parameter devInfo=i indicates U(i,i) = 0.

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

No matter LU factorization failed or not, the output parameter devIpiv contains pivoting sequence, row i is interchanged with row devIpiv(i).

cusolverStatus_t 
cusolverDnSgetrs(cusolverDnHandle_t handle,
           cublasOperation_t trans,
           int n,
           int nrhs,
           const float *A,
           int lda,
           const int *devIpiv,
           float *B,
           int ldb,
           int *devInfo );

cusolverStatus_t 
cusolverDnDgetrs(cusolverDnHandle_t handle,
           cublasOperation_t trans,
           int n,
           int nrhs,
           const double *A,
           int lda,
           const int *devIpiv,
           double *B,
           int ldb,
           int *devInfo );

cusolverStatus_t 
cusolverDnCgetrs(cusolverDnHandle_t handle,
           cublasOperation_t trans,
           int n,
           int nrhs,
           const cuComplex *A,
           int lda,
           const int *devIpiv,
           cuComplex *B,
           int ldb,
           int *devInfo );

cusolverStatus_t 
cusolverDnZgetrs(cusolverDnHandle_t handle,
           cublasOperation_t trans,
           int n,
           int nrhs,
           const cuDoubleComplex *A,
           int lda,
           const int *devIpiv,
           cuDoubleComplex *B,
           int ldb,
           int *devInfo );

This function solves a linear system of multiple right-hand sides

where A is a n×n matrix, and was LU-factored by getrf, that is, lower trianular part of A is L, and upper triangular part (including diagonal elements) of A is U. B is a n×nrhs right-hand side matrix.

The input parameter trans is defined by

$\text{op}(A)=\left\{\begin{array}{ll}A& \text{if }\mathsf{\text{trans == CUBLAS\_OP\_N}}\\ A^T& \text{if }\mathsf{\text{trans == CUBLAS\_OP\_T}}\\ A^H& \text{if }\mathsf{\text{trans == CUBLAS\_OP\_C}}\end{array}\right.$

The input parameter devIpiv is an output of getrf. It contains pivot indices, which are used to permutate right-hand sides.

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

cusolverStatus_t 
cusolverDnSgeqrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      float *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnDgeqrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      double *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnCgeqrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      cuComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnZgeqrf_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      cuDoubleComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnSgeqrf(cusolverDnHandle_t handle,
           int m,
           int n,
           float *A,
           int lda,
           float *TAU,
           float *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnDgeqrf(cusolverDnHandle_t handle,
           int m,
           int n,
           double *A,
           int lda,
           double *TAU,
           double *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnCgeqrf(cusolverDnHandle_t handle,
           int m,
           int n,
           cuComplex *A,
           int lda,
           cuComplex *TAU,
           cuComplex *Workspace,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnZgeqrf(cusolverDnHandle_t handle,
           int m,
           int n,
           cuDoubleComplex *A,
           int lda,
           cuDoubleComplex *TAU,
           cuDoubleComplex *Workspace,
           int Lwork,
           int *devInfo );

This function computes the QR factorization of a m×n matrix

where A is a m×n matrix, Q is a m×n matrix, and R is a n×n upper triangular matrix.

The user has to provide working space which is pointed by input parameter Workspace. The input parameter Lwork is size of the working space, and it is returned by geqrf_bufferSize().

The matrix R is overwritten in upper triangular part of A, including diagonal elements.

The matrix Q is not formed explicitly, instead, a sequence of householder vectors are stored in lower triangular part of A. The leading nonzero element of householder vector is assumed to be 1 such that output parameter TAU contains the scaling factor τ. If v is original householder vector, q is the new householder vector corresponding to τ, satisying the following relation

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

cusolverStatus_t 
cusolverDnSormqr(cusolverDnHandle_t handle,
           cublasSideMode_t side,
           cublasOperation_t trans,
           int m,
           int n,
           int k,
           const float *A,
           int lda,
           const float *tau,
           float *C,
           int ldc,
           float *work,
           int lwork,
           int *devInfo);

cusolverStatus_t
cusolverDnDormqr(cusolverDnHandle_t handle,
           cublasSideMode_t side,
           cublasOperation_t trans,
           int m,
           int n,
           int k,
           const double *A,
           int lda,
           const double *tau,
           double *C,
           int ldc,
           double *work,
           int lwork,
           int *devInfo);

cusolverStatus_t 
cusolverDnCunmqr(cusolverDnHandle_t handle,
           cublasSideMode_t side,
           cublasOperation_t trans,
           int m,
           int n,
           int k,
           const cuComplex *A,
           int lda,
           const cuComplex *tau,
           cuComplex *C,
           int ldc,
           cuComplex *work,
           int lwork,
           int *devInfo);

cusolverStatus_t 
cusolverDnZunmqr(cusolverDnHandle_t handle,
           cublasSideMode_t side,
           cublasOperation_t trans,
           int m,
           int n,
           int k,
           const cuDoubleComplex *A,
           int lda,
           const cuDoubleComplex *tau,
           cuDoubleComplex *C,
           int ldc,
           cuDoubleComplex *work,
           int lwork,
           int *devInfo);

This function overwrites m×n matrix C by

where Q is a unitary matrix formed by a sequence of elementary reflection vectors from QR factorization of A. Also for Q

The user has to provide working space which is pointed by input parameter work. The input parameter lwork is size of the working space, and it is returned by geqrf_bufferSize().

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

The user can combine geqrf, ormqr and trsm to complete a linear solver or a least-square solver. Please refer to appendix C.1.

cusolverStatus_t 
cusolverDnSsytrf_bufferSize(cusolverDnHandle_t handle,
                      int n,
                      float *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnDsytrf_bufferSize(cusolverDnHandle_t handle,
                      int n,
                      double *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnCsytrf_bufferSize(cusolverDnHandle_t handle,
                      int n,
                      cuComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnZsytrf_bufferSize(cusolverDnHandle_t handle,
                      int n,
                      cuDoubleComplex *A,
                      int lda,
                      int *Lwork );

cusolverStatus_t 
cusolverDnSsytrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           float *A,
           int lda,
           int *ipiv,
           float *work,
           int lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnDsytrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           double *A,
           int lda,
           int *ipiv,
           double *work,
           int lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnCsytrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           cuComplex *A,
           int lda,
           int *ipiv,
           cuComplex *work,
           int lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnZsytrf(cusolverDnHandle_t handle,
           cublasFillMode_t uplo,
           int n,
           cuDoubleComplex *A,
           int lda,
           int *ipiv,
           cuDoubleComplex *work,
           int lwork,
           int *devInfo );

This function computes the Bunch-Kaufman factorization of a n×n symmetric indefinite matrix

A is a n×n symmetric matrix, only lower or upper part is meaningful. The input parameter uplo which part of the matrix is used. The function would leave other part untouched.

If input parameter uplo is CUBLAS_FILL_MODE_LOWER, only lower triangular part of A is processed, and replaced by lower triangular factor L and block diagonal matrix D. Each block of D is either 1x1 or 2x2 block, depending on pivoting.

If input parameter uplo is CUBLAS_FILL_MODE_UPPER, only upper triangular part of A is processed, and replaced by upper triangular factor U and block diagonal matrix D.

The user has to provide working space which is pointed by input parameter work. The input parameter lwork is size of the working space, and it is returned by sytrf_bufferSize().

If Bunch-Kaufman factorization failed, i.e. A is singular. The output parameter devInfo = i would indicate D(i,i)=0.

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

The output parameter devIpiv contains pivoting sequence. If devIpiv(i) = k > 0, D(i,i) is 1x1 block, and i-th row/column of A is interchanged with k-th row/column of A. If uplo is CUSBLAS_FILL_MODE_UPPER and devIpiv(i-1) = devIpiv(i) = -m < 0, D(i-1:i,i-1:i) is a 2x2 block, and (i-1)-th row/column is interchanged with m-th row/column. If uplo is CUSBLAS_FILL_MODE_LOWER and devIpiv(i+1) = devIpiv(i) = -m < 0, D(i:i+1,i:i+1) is a 2x2 block, and (i+1)-th row/column is interchanged with m-th row/column.

cusolverStatus_t 
cusolverDnSgebrd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnDgebrd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t
cusolverDnCgebrd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnZgebrd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnSgebrd(cusolverDnHandle_t handle,
           int m,
           int n,
           float *A,
           int lda,
           float *D,
           float *E,
           float *TAUQ,
           float *TAUP,
           float *Work,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnDgebrd(cusolverDnHandle_t handle,
           int m,
           int n,
           double *A,
           int lda,
           double *D,
           double *E,
           double *TAUQ,
           double *TAUP,
           double *Work,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnCgebrd(cusolverDnHandle_t handle,
           int m,
           int n,
           cuComplex *A,
           int lda,
           float *D,
           float *E,
           cuComplex *TAUQ,
           cuComplex *TAUP,
           cuComplex *Work,
           int Lwork,
           int *devInfo );

cusolverStatus_t 
cusolverDnZgebrd(cusolverDnHandle_t handle,
           int m,
           int n,
           cuDoubleComplex *A,
           int lda,
           double *D,
           double *E,
           cuDoubleComplex *TAUQ,
           cuDoubleComplex *TAUP,
           cuDoubleComplex *Work,
           int Lwork,
           int *devInfo );

This function reduces a general real m×n matrix A to upper or lower bidiagonal form B by an orthogonal transformation: $Q^H*A*P=B$

If m>=n, B is upper bidiagonal; if m<n, B is lower bidiagonal.

The matrix Q and P are overwritten into matrix A in the following sense:

if m>=n, the diagonal and the first superdiagonal are overwritten with the upper bidiagonal matrix B; the elements below the diagonal, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors, and the elements above the first superdiagonal, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors.

if m<n, the diagonal and the first subdiagonal are overwritten with the lower bidiagonal matrix B; the elements below the first subdiagonal, with the array TAUQ, represent the orthogonal matrix Q as a product of elementary reflectors, and the elements above the diagonal, with the array TAUP, represent the orthogonal matrix P as a product of elementary reflectors.

The user has to provide working space which is pointed by input parameter Work. The input parameter Lwork is size of the working space, and it is returned by gebrd_bufferSize().

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong.

Remark: gebrd only supports m>=n.

cusolverStatus_t 
cusolverDnSgesvd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnDgesvd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnCgesvd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t
cusolverDnZgesvd_bufferSize(cusolverDnHandle_t handle,
                      int m,
                      int n,
                      int *Lwork );

cusolverStatus_t 
cusolverDnSgesvd (cusolverDnHandle_t handle,
            char jobu,
            char jobvt,
            int m,
            int n,
            float *A,
            int lda,
            float *S,
            float *U,
            int ldu,
            float *VT,
            int ldvt,
            float *Work,
            int Lwork,
            float *rwork,
            int  *devInfo);

cusolverStatus_t 
cusolverDnDgesvd (cusolverDnHandle_t handle,
            char jobu,
            char jobvt,
            int m,
            int n,
            double *A,
            int lda,
            double *S,
            double *U,
            int ldu,
            double *VT,
            int ldvt,
            double *Work,
            int Lwork,
            double *rwork,
            int *devInfo);

cusolverStatus_t 
cusolverDnCgesvd (cusolverDnHandle_t handle,
            char jobu,
            char jobvt,
            int m,
            int n,
            cuComplex *A,
            int lda,
            float *S,
            cuComplex *U,
            int ldu,
            cuComplex *VT,
            int ldvt,
            cuComplex *Work,
            int Lwork,
            float *rwork,
            int *devInfo);

cusolverStatus_t 
cusolverDnZgesvd (cusolverDnHandle_t handle,
            char jobu,
            char jobvt,
            int m,
            int n,
            cuDoubleComplex *A,
            int lda,
            double *S,
            cuDoubleComplex *U,
            int ldu,
            cuDoubleComplex *VT,
            int ldvt,
            cuDoubleComplex *Work,
            int Lwork,
            double *rwork,
            int *devInfo);

This function computes the singular value decomposition (SVD) of a m×n matrix A and corresponding the left and/or right singular vectors. The SVD is written

where Σ is an m×n matrix which is zero except for its min(m,n) diagonal elements, U is an m×m unitary matrix, and V is an n×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.

The user has to provide working space which is pointed by input parameter Work. The input parameter Lwork is size of the working space, and it is returned by gesvd_bufferSize().

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong. if bdsqr did not converge, devInfo specifies how many superdiagonals of an intermediate bidiagonal form B did not converge to zero.

Note that the routine returns $V^H$ , not V.

Remark: gesvd only supports m>=n.

Remark: gesvd returns matrix $U$ and $V^H$ .

cusolverStatus_t
cusolverSpCreate(cusolverSpHandle_t *handle)

This function initializes the cuSolverSP library and creates a handle on the cuSolver context. It must be called before any other cuSolverSP API function is invoked. It allocates hardware resources necessary for accessing the GPU.

cusolverStatus_t
cusolverSpDestroy(cusolverSpHandle_t handle)

This function releases CPU-side resources used by the cuSolverSP library.

cusolverStatus_t
cusolverSpSetStream(cusolverSpHandle_t handle, cudaStream_t streamId)

This function sets the stream to be used by the cuSolverSP library to execute its routines.

cusolverStatus_t 
cusolverSpXcsrissymHost(cusolverSpHandle_t handle,
              int m,
              int nnzA,
              const cusparseMatDescr_t descrA,
              const int *csrRowPtrA,
              const int *csrEndPtrA,
              const int *csrColIndA,
              int *issym);

This function checks if A has symmetric pattern or not. The output parameter issym reports 1 if A is symmetric; otherwise, it reports 0.

The matrix A is an m×m sparse matrix that is defined in CSR storage format by the four arrays csrValA, csrRowPtrA, csrEndPtrA and csrColIndA.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL.

The csrlsvlu and csrlsvqr do not accept non-general matrix. the user has to extend the matrix into its missing upper/lower part, otherwise the result is not expected. The user can use csrissym to check if the matrix has symmetric pattern or not.

Remark 1: only CPU path is provided.

Remark 2: the user has to check returned status to get valid information. The function converts A to CSC format and compare CSR and CSC format. If the CSC failed because of insufficient resources, issym is undefined, and this state can only be detected by the return status code.

cusolverStatus_t 
cusolverSpScsrlsvlu[Host](cusolverSpHandle_t handle,
                 int n,
                 int nnzA,
                 const cusparseMatDescr_t descrA,
                 const float *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const float *b,
                 float tol,
                 int reorder,
                 float *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpDcsrlsvlu[Host](cusolverSpHandle_t handle,
                 int n,
                 int nnzA,
                 const cusparseMatDescr_t descrA,
                 const double *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const double *b,
                 double tol,
                 int reorder,
                 double *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpCcsrlsvlu[Host](cusolverSpHandle_t handle,
                 int n,
                 int nnzA,
                 const cusparseMatDescr_t descrA,
                 const cuComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const cuComplex *b,
                 float tol,
                 int reorder,
                 cuComplex *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpZcsrlsvlu[Host](cusolverSpHandle_t handle,
                 int n,
                 int nnzA,
                 const cusparseMatDescr_t descrA,
                 const cuDoubleComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const cuDoubleComplex *b,
                 double tol,
                 int reorder,
                 cuDoubleComplex *x,
                 int *singularity);

This function solves the linear system

A is an n×n sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA. b is the right-hand-side vector of size n, and x is the solution vector of size n.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If matrix A is symmetric/Hermitian and only lower/upper part is used or meaningful, the user has to extend the matrix into its missing upper/lower part, otherwise the result is wrong.

The linear system is solved by sparse LU with partial pivoting,

The performance of LU factorization mainly depends on zero fill-in. cusolver library provides symrcm to reduce zero fill-in. The input parameter reorder can enable symrcm if reorder is nonzero, otherwise, no reordering is performed.

If reorder is nonzero, csrlsvlu does

where $Q=\mathrm{symrcm}(A+A^T)$ .

If A is singular under given tolerance (max(tol,0)), then some diagonal elements of U is zero, i.e.

The output parameter singularity is the smallest index of such j. If A is non-singular, singularity is -1. The index is base-0, independent of base index of A. For example, if 2nd column of A is the same as first column, then A is singular and singularity = 1 which means U(1,1)≈0.

Remark 1: csrlsvlu does traditional LU with partial pivoting, the pivot of k-th column is determined dynamically based on the k-th column of intermediate matrix. csrlsvlu follows Gilbert and Peierls's algorithm [4] which uses depth-first-search and topological ordering to solve triangular system (Davis also describes this algorithm in detail in his book [1]). Before performing LU factorization, csrlsvlu over-estimates size of L and U, and allocates a buffer to contain factors L and U. There is no easy way to find a tight upper bound of size of LU, George and Ng [5] proves that sparsity pattern of cholesky factor of $A*A^T$ is a superset of sparsity pattern of L and U. Furthermore, they provides an algorithm to find sparisty pattern of QR factorization which is a superset of LU [6]. csrlsvlu uses QR factorization to estimate size of LU in the analysis phase. The cost of analysis phase is mainly on figuring out sparsity pattern of householder vectors in QR factorization. If system memory is insufficient to keep sparsity pattern of QR, csrlsvlu returns CUSOLVER_STATUS_ALLOC_FAILED. If the matrix is not banded, it is better to enable reordering to avoid CUSOLVER_STATUS_ALLOC_FAILED.

Remark 2: minimum degree ordering is the well-known technique to reduce zero fill-in of QR factorization. However in most cases, symrcm still performs well.

Remark 3: The user can reorder the matrix by other reordering scheme before calling csrlsvlu For example, symrcm is not a good choice on circuit simulation. The user can try column AMD.

Remark 4: if matrix A is singular and reorder is on, the output parameter singularity reports singularity of $P*A*P^T$ .

Remark 5: only CPU (Host) path is provided.

Remark 6: multithreaded csrlsvlu is not avaiable yet. If QR does not incur much zero fill-in, csrlsvqr would be faster than csrlsvlu.

cusolverStatus_t 
cusolverSpScsrlsvqr[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const float *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const float *b,
                 float tol,
                 int reorder,
                 float *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpDcsrlsvqr[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const double *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const double *b,
                 double tol,
                 int reorder,
                 double *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpCcsrlsvqr[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const cuComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const cuComplex *b,
                 float tol,
                 int reorder,
                 cuComplex *x,
                 int *singularity);

cusolverStatus_t 
cusolverSpZcsrlsvqr[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const cuDoubleComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 const cuDoubleComplex *b,
                 double tol,
                 int reorder,
                 cuDoubleComplex *x,
                 int *singularity);

This function solves the linear system

A is an m×m sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA. b is the right-hand-side vector of size m, and x is the solution vector of size m.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If matrix A is symmetric/Hermitian and only lower/upper part is used or meaningful, the user has to extend the matrix into its missing upper/lower part, otherwise the result is wrong.

The linear system is solved by sparse QR factorization,

If A is singular under given tolerance (max(tol,0)), then some diagonal elements of R is zero, i.e.

The output parameter singularity is the smallest index of such j. If A is non-singular, singularity is -1. The index is consistent with base index of A. For example, if 2nd column of A is the same as first column, then A is singular. If base index is 1, then singularity = 2 which means R(2,2)≈0. If base index is 0, then singularity = 1 which means R(1,1)≈0.

Remark: the parameter reorder has no effect because reordering is not implemented yet.

cusolverStatus_t 
cusolverSpScsrlsvchol[Host](cusolverSpHandle_t handle,
                     int m,
                     int nnz,
                     const cusparseMatDescr_t descrA,
                     const float *csrVal,
                     const int *csrRowPtr,
                     const int *csrColInd,
                     const float *b,
                     float tol,
                     int reorder,
                     float *x,
                     int *singularity);

cusolverStatus_t 
cusolverSpDcsrlsvchol[Host](cusolverSpHandle_t handle,
                     int m,
                     int nnz,
                     const cusparseMatDescr_t descrA,
                     const double *csrVal,
                     const int *csrRowPtr,
                     const int *csrColInd,
                     const double *b,
                     double tol,
                     int reorder,
                     double *x,
                     int *singularity);

cusolverStatus_t 
cusolverSpCcsrlsvchol[Host](cusolverSpHandle_t handle,
                     int m,
                     int nnz,
                     const cusparseMatDescr_t descrA,
                     const cuComplex *csrVal,
                     const int *csrRowPtr,
                     const int *csrColInd,
                     const cuComplex *b,
                     float tol,
                     int reorder,
                     cuComplex *x,
                     int *singularity);

cusolverStatus_t 
cusolverSpZcsrlsvchol[Host](cusolverSpHandle_t handle,
                     int m,
                     int nnz,
                     const cusparseMatDescr_t descrA,
                     const cuDoubleComplex *csrVal,
                     const int *csrRowPtr,
                     const int *csrColInd,
                     const cuDoubleComplex *b,
                     double tol,
                     int reorder,
                     cuDoubleComplex *x,
                     int *singularity);

This function solves the linear system

A is an m×m symmetric postive definite sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA. b is the right-hand-side vector of size m, and x is the solution vector of size m.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL and upper triangular part of A is ignored. In other words, suppose input matrix A is decomposed as $A=L+D+U$ , where L is lower triangular, D is diagonal and U is upper triangular. The function would ignore U and regard A as a symmetric matrix with the formula $A=L+D+L^H$ . Thereafter, A is assumed symmetric/Hermitian.

The linear system is solved by sparse Cholesky factorization,

where G is the Cholesky factor, a lower triangular matrix.

The output parameter singularity has two meanings:

If A is not postive definite, there exists some integer k such that A(0:k, 0:k) is not positive definite. singularity is the minimum of such k.

If A is postive definite but near singular under tolerance (max(tol,0)), i.e. there exists some integer k such that $G\left(\begin{array}{c}\mathrm{k,k}\end{array}\right)<=\mathrm{tol}$ . singularity is the minimum of such k.

Remark 1: singularity is base-0. If A is positive definite and not near singular under tolerance, singularity is -1.

Remark 2: if the user wants to know if A is postive definite or not, tol=0 is enough.

Remark 3: reordering can greatly affect zero fill-in and have a huge performance impact. cuSolver does not provide reordering scheme, the user has to apply some reordering scheme to have decent performance, for example, symrcm or symamd.

Remark 4: the function works for in-place (x and b point to the same memory block) and out-of-place.

Remark 5: the function only works on 32-bit index, if matrix G has large zero fill-in such that number of nonzeros is bigger than $2^{31}$ , then CUSOLVER_STATUS_ALLOC_FAILED is returned.

Remark 6: the parameter reorder has no effect because reordering is not implemented yet.

cusolverStatus_t 
cusolverSpScsrlsqvqr[Host](cusolverSpHandle_t handle,
                  int m,
                  int n,
                  int nnz,
                  const cusparseMatDescr_t descrA,
                  const float *csrValA,
                  const int *csrRowPtrA,
                  const int *csrColIndA,
                  const float *b,
                  float tol,
                  int *rankA,
                  float *x,
                  int *p,
                  float *min_norm);

cusolverStatus_t 
cusolverSpDcsrlsqvqr[Host](cusolverSpHandle_t handle,
                  int m,
                  int n,
                  int nnz,
                  const cusparseMatDescr_t descrA,
                  const double *csrValA,
                  const int *csrRowPtrA,
                  const int *csrColIndA,
                  const double *b,
                  double tol,
                  int *rankA,
                  double *x,
                  int *p,
                  double *min_norm);

cusolverStatus_t 
cusolverSpCcsrlsqvqr[Host](cusolverSpHandle_t handle,
                  int m,
                  int n,
                  int nnz,
                  const cusparseMatDescr_t descrA,
                  const cuComplex *csrValA,
                  const int *csrRowPtrA,
                  const int *csrColIndA,
                  const cuComplex *b,
                  float tol,
                  int *rankA,
                  cuComplex *x,
                  int *p,
                  float *min_norm);

cusolverStatus_t 
cusolverSpZcsrlsqvqr[Host](cusolverSpHandle_t handle,
                  int m,
                  int n,
                  int nnz,
                  const cusparseMatDescr_t descrA,
                  const cuDoubleComplex *csrValA,
                  const int *csrRowPtrA,
                  const int *csrColIndA,
                  const cuDoubleComplex *b,
                  double tol,
                  int *rankA,
                  cuDoubleComplex *x,
                  int *p,
                  double *min_norm);

This function solves the following least-square problem

A is an m×n sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA. b is the right-hand-side vector of size m, and x is the least-square solution vector of size n.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If A is square, symmetric/Hermitian and only lower/upper part is used or meaningful, the user has to extend the matrix into its missing upper/lower part, otherwise the result is wrong.

This function only works if m is greater or equal to n, in other words, A is a tall matrix.

The least-square problem is solved by sparse QR factorization with column pivoting,

If A is of full rank (i.e. all columns of A are linear independent), then matrix P is an identity. Suppose rank of A is k, less than n, the permutation matrix P reorders columns of A in the following sense:

where $R_{11}$ and A have the same rank, but $R_{22}$ is almost zero, i.e. every column of $A_2$ is linear combination of $A_1$ .

The input parameter tol decides numerical rank. The absolute value of every entry in $R_{22}$ is less than or equal to tolerance=max(tol,0).

The output parameter rankA denotes numerical rank of A.

Suppose $y=P*x$ and $c=Q^H*b$ , the least square problem can be reformed by

or in matrix form

The output parameter min_norm is $||c_2||$ , which is minimum value of least-square problem.

If A is not of full rank, above equation does not have a unique solution. The least-square problem is equivalent to

Or equivalently another least-square problem

The output parameter x is $P^T*y$ , the solution of least-square problem.

The output parameter p is a vector of size n. It corresponds to a permutation matrix P. p(i)=j means (P*x)(i) = x(j). If A is of full rank, p=0:n-1.

Remark 1: p is always base 0, independent of base index of A.

Remark 2: only CPU (Host) path is provided.

cusolverStatus_t 
cusolverSpScsreigvsi[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const float *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 float mu0,
                 const float *x0,
                 int maxite,
                 float tol,
                 float *mu,
                 float *x);

cusolverStatus_t
cusolverSpDcsreigvsi[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const double *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 double mu0,
                 const double *x0,
                 int maxite,
                 double tol,
                 double *mu,
                 double *x);

cusolverStatus_t 
cusolverSpCcsreigvsi[Host](cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const cuComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 cuComplex mu0,
                 const cuComplex *x0,
                 int maxite,
                 float tol,
                 cuComplex *mu,
                 cuComplex *x);

cusolverStatus_t 
cusolverSpZcsreigvsi(cusolverSpHandle_t handle,
                 int m,
                 int nnz,
                 const cusparseMatDescr_t descrA,
                 const cuDoubleComplex *csrValA,
                 const int *csrRowPtrA,
                 const int *csrColIndA,
                 cuDoubleComplex mu0,
                 const cuDoubleComplex *x0,
                 int maxite,
                 double tol,
                 cuDoubleComplex *mu,
                 cuDoubleComplex *x);

This function solves the simple eigenvalue problem $A*x=\lambda *x$ by shift-inverse method.

A is an m×m sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA. The output paramter x is the approximated eigenvector of size m,

The following shift-inverse method corrects eigenpair step-by-step until convergence.

It accepts several parameters:

mu0 is an initial guess of eigenvalue. The shift-inverse method will converge to the eigenvalue mu nearest mu0 if mu is a singleton. Otherwise, the shift-inverse method may not converge.

x0 is an initial eigenvector. If the user has no preference, just chose x0 randomly. x0 must be nonzero. It can be non-unit length.

tol is the tolerance to decide convergence. If tol is less than zero, it would be treated as zero.

maxite is maximum number of iterations. It is useful when shift-inverse method does not converge because the tolerance is too small or the desired eigenvalue is not a singleton.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If A is symmetric/Hermitian and only lower/upper part is used or meaningful, the user has to extend the matrix into its missing upper/lower part, otherwise the result is wrong.

Remark 1: [cu|h]solver[S|D]csreigvsi only allows mu0 as a real number. This works if A is symmetric. Otherwise, the non-real eigenvalue has a conjugate counterpart on the complex plan, and shift-inverse method would not converge to such eigevalue even the eigenvalue is a singleton. The user has to extend A to complex numbre and call [cu|h]solver[C|Z]csreigvsi with mu0 not on real axis.

Remark 2: the tolerance tol should not be smaller than |mu0|*eps, where eps is machine zero. Otherwise, shift-inverse may not converge because of small tolerance.

cusolverStatus_t 
solverspScsreigs[Host](cusolverSpHandle_t handle,
                int m,
                int nnz,
                const cusparseMatDescr_t descrA,
                const float *csrValA,
                const int *csrRowPtrA,
                const int *csrColIndA,
                cuComplex left_bottom_corner,
                cuComplex right_upper_corner,
                int *num_eigs);

cusolverStatus_t 
cusolverSpDcsreigs[Host](cusolverSpHandle_t handle,
                int m,
                int nnz,
                const cusparseMatDescr_t descrA,
                const double *csrValA,
                const int *csrRowPtrA,
                const int *csrColIndA,
                cuDoubleComplex left_bottom_corner,
                cuDoubleComplex right_upper_corner,
                int *num_eigs);

cusolverStatus_t 
cusolverSpCcsreigs[Host](cusolverSpHandle_t handle,
                int m,
                int nnz,
                const cusparseMatDescr_t descrA,
                const cuComplex *csrValA,
                const int *csrRowPtrA,
                const int *csrColIndA,
                cuComplex left_bottom_corner,
                cuComplex right_upper_corner,
                int *num_eigs);

cusolverStatus_t 
cusolverSpZcsreigs[Host](cusolverSpHandle_t handle,
                int m,
                int nnz,
                const cusparseMatDescr_t descrA,
                const cuDoubleComplex *csrValA,
                const int *csrRowPtrA,
                const int *csrColIndA,
                cuDoubleComplex left_bottom_corner,
                cuDoubleComplex right_upper_corner,
                int *num_eigs);

This function computes number of algebraic eigenvalues in a given box B by contour integral

where closed line C is boundary of the box B which is a rectangle specified by two points, one is left bottom corner (input parameter left_botoom_corner) and the other is right upper corner (input parameter right_upper_corner). P(z)=det(A - z*I) is the characteristic polynomial of A.

A is an m×m sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA.

The output parameter num_eigs is number of algebraic eigenvalues in the box B. This number may not be accurate due to several reasons:

1. the contour C is close to some eigenvalues or even passes through some eigenvalues.

2. the numerical integration is not accurate due to coarse grid size. The default resolution is 1200 grids along contour C uniformly.

Even though csreigs may not be accurate, it still can give the user some idea how many eigenvalues in a region where the resolution of disk theorem is bad. For example, standard 3-point stencil of finite difference of Laplacian operator is a tridiagonal matrix, and disk theorem would show "all eigenvalues are in the interval [0, 4*N^2]" where N is number of grids. In this case, csreigs is useful for any interval inside [0, 4*N^2].

Remark 1: if A is symmetric in real or hermitian in complex, all eigenvalues are real. The user still needs to specify a box, not an interval. The height of the box can be much smaller than the width.

Remark 2: only CPU (Host) path is provided.

cusolverStatus_t 
cusolverSpXcsrsymrcmHost(cusolverSpHandle_t handle,
             int n,
             int nnzA,
             const cusparseMatDescr_t descrA,
             const int *csrRowPtrA,
             const int *csrColIndA,
             int *p);

This function implements Symmetric Reverse Cuthill-McKee permutation. It returns a permutation vector p such that A(p,p) would concentrate nonzeros to diagonal. This is equivalent to symrcm in MATLAB, however the result may not be the same because of different heuristics in the pseudoperipheral finder. The cuSolverSP library implements symrcm based on the following two papers:

E. Chuthill and J. McKee, reducing the bandwidth of sparse symmetric matrices, ACM '69 Proceedings of the 1969 24th national conference, Pages 157-172

Alan George, Joseph W. H. Liu, An Implementation of a Pseudoperipheral Node Finder, ACM Transactions on Mathematical Software (TOMS) Volume 5 Issue 3, Sept. 1979, Pages 284-295

The output parameter p is an integer array of n elements. It represents a permutation array and it indexed using the base-0 convention. The permutation array p corresponds to a permutation matrix P, and satisfies the following relation:

A is an n×n sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA, and csrColIndA.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Internally rcm works on $A+A^T$ , the user does not need to extend the matrix if the matrix is not symmetric.

Remark 1: only CPU (Host) path is provided.

cusolverStatus_t 
cusolverSpXcsrperm_bufferSizeHost(cusolverSpHandle_t handle,
                          int m,
                          int n,
                          int nnzA,
                          const cusparseMatDescr_t descrA,
                          int *csrRowPtrA,
                          int *csrColIndA,
                          const int *p,
                          const int *q,
                          size_t *bufferSizeInBytes);

cusolverStatus_t 
cusolverSpXcsrpermHost(cusolverSpHandle_t handle,
                int m,
                int n,
                int nnzA,
                const cusparseMatDescr_t descrA,
                int *csrRowPtrA,
                int *csrColIndA,
                const int *p,
                const int *q,
                int *map,
                void *pBuffer);

Given a left permutation vector p which corresponds to permutation matrix P and a right permutation vector q which corresponds to permutation matrix Q, this function computes permutation of matrix A by

A is an m×n sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA and csrColIndA.

The operation is in-place, i.e. the matrix A is overwritten by B.

The permutation vector p and q are base 0. p performs row permutation while q performs column permutation. One can also use MATLAB command $B=\mathrm{A(p,q)}$ to permutate matrix A.

This function only computes sparsity pattern of B. The user can use parameter map to get csrValB as well. The parameter map is an input/output. If the user sets map=0:1:(nnzA-1) before calling csrperm, csrValB=csrValA(map).

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If A is symmetric and only lower/upper part is provided, the user has to pass $A+A^T$ into this function.

This function requires a buffer size returned by csrperm_bufferSize(). The address of pBuffer must be a multiple of 128 bytes. If it is not, CUSOLVER_STATUS_INVALID_VALUE is returned.

For example, if matrix A is

and left permutation vector p=(0,2,1), right permutation vector q=(2,1,0), then $P*A*Q^T$ is

Remark 1: only CPU (Host) path is provided.

Remark 2: the user can combine csrsymrcm and csrperm to get $P*A*P^T$ which has less zero fill-in during QR factorization.

cusolverStatus_t 
cusolverSpCreateCsrqrInfo(csrqrInfo_t *info);

cusolverStatus_t 
cusolverSpDestroyCsrqrInfo(csrqrInfo_t info);

cusolverStatus_t 
cusolverSpXcsrqrAnalysisBatched(cusolverSpHandle_t handle,
                           int m,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrqrInfo_t info);

cusolverStatus_t 
cusolverSpScsrqrBufferInfoBatched(cusolverSpHandle_t handle,
                           int m,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const float *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           int batchSize,
                           csrqrInfo_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpDcsrqrBufferInfoBatched(cusolverSpHandle_t handle,
                           int m,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const double *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           int batchSize,
                           csrqrInfo_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpCcsrqrBufferInfoBatched(cusolverSpHandle_t handle,
                           int m,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           int batchSize,
                           csrqrInfo_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpZcsrqrBufferInfoBatched(cusolverSpHandle_t handle,
                           int m,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuDoubleComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           int batchSize,
                           csrqrInfo_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpScsrqrsvBatched(cusolverSpHandle_t handle,
                        int m,
                        int n,
                        int nnzA,
                        const cusparseMatDescr_t descrA,
                        const float *csrValA,
                        const int *csrRowPtrA,
                        const int *csrColIndA,
                        const float *b,
                        float *x,
                        int batchSize,
                        csrqrInfo_t info,
                        void *pBuffer);

cusolverStatus_t 
cusolverSpDcsrqrsvBatched(cusolverSpHandle_t handle,
                        int m,
                        int n,
                        int nnz,
                        const cusparseMatDescr_t descrA,
                        const double *csrValA,
                        const int *csrRowPtrA,
                        const int *csrColIndA,
                        const double *b,
                        double *x,
                        int batchSize,
                        csrqrInfo_t info,
                        void *pBuffer);

cusolverStatus_t 
cusolverSpCcsrqrsvBatched(cusolverSpHandle_t handle,
                        int m,
                        int n,
                        int nnzA,
                        const cusparseMatDescr_t descrA,
                        const cuComplex *csrValA,
                        const int *csrRowPtrA,
                        const int *csrColIndA,
                        const cuComplex *b,
                        cuComplex *x,
                        int batchSize,
                        csrqrInfo_t info,
                        void *pBuffer);

cusolverStatus_t 
cusolverSpZcsrqrsvBatched(cusolverSpHandle_t handle,
                        int m,
                        int n,
                        int nnzA,
                        const cusparseMatDescr_t descrA,
                        const cuDoubleComplex *csrValA,
                        const int *csrRowPtrA,
                        const int *csrColIndA,
                        const cuDoubleComplex *b,
                        cuDoubleComplex *x,
                        int batchSize,
                        csrqrInfo_t info,
                        void *pBuffer);

The batched sparse QR factorization is used to solve either a set of least-squares problems

or a set of linear systems

where each $A_j$ is a m×n sparse matrix that is defined in CSR storage format by the four arrays csrValA, csrRowPtrA and csrColIndA.

The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. If A is symmetric and only lower/upper part is prvided, the user has to pass $A+A^H$ into this function.

The prerequisite to use batched sparse QR has two-folds. First all matrices $A_j$ must have the same sparsity pattern. Second, no column pivoting is used in least-square problem, so the solution is valid only if $A_j$ is of full rank for all j = 1,2,..., batchSize . All matrices have the same sparity pattern, so only one copy of csrRowPtrA and csrColIndA is used. But the array csrValA stores coefficients of $A_j$ one after another. In other words, csrValA[k*nnzA : (k+1)*nnzA] is the value of $A_k$ .

The batched QR uses opaque data structure csrqrInfo to keep intermediate data, for example, matrix Q and matrix R of QR factorization. The user needs to create csrqrInfo first by cusolverSpCreateCsrqrInfo before any function in batched QR operation. The csrqrInfo would not release internal data until cusolverSpDestroyCsrqrInfo is called.

There are three routines in batched sparse QR, cusolverSpXcsrqrAnalysisBatched, cusolverSp[S|D|C|Z]csrqrBufferInfoBatched and cusolverSp[S|D|C|Z]csrqrsvBatched.

First, cusolverSpXcsrqrAnalysisBatched is the analysis phase, used to analyze sparsity pattern of matrix Q and matrix R of QR factorization. Also parallelism is extracted during analysis phase. Once analysis phase is done, the size of working space to perform QR is known. However cusolverSpXcsrqrAnalysisBatched uses CPU to analyze the structure of matrix A, and this may consume a lot of memory. If host memory is not sufficient to finish the analysis, CUSOLVER_STATUS_ALLOC_FAILED is returned. The required memory for analysis is proportional to zero fill-in in QR factorization. The user may need to perform some kind of reordering to minimize zero fill-in, for example, colamd or symrcm in MATLAB. cuSolverSP library provides symrcm (cusolverSpXcsrsymrcm).

Second, the user needs to choose proper batchSize and to prepare working space for sparse QR. There are two memory blocks used in batched sparse QR. One is internal memory block used to store matrix Q and matrix R. The other is working space used to perform numerical factorization. The size of the former is proportional to batchSize, and the size is specified by returned parameter internalDataInBytes of cusolverSp[S|D|C|Z]csrqrBufferInfoBatched. while the size of the latter is almost independent of batchSize, and the size is specified by returned parameter workspaceInBytes of cusolverSp[S|D|C|Z]csrqrBufferInfoBatched. The internal memory block is allocated implicitly during first call of cusolverSp[S|D|C|Z]csrqrsvBatched. The user only needs to allocate working space for cusolverSp[S|D|C|Z]csrqrsvBatched.

Instead of trying all batched matrices, the user can find maximum batchSize by querying cusolverSp[S|D|C|Z]csrqrBufferInfoBatched. For example, the user can increase batchSize till summation of internalDataInBytes and workspaceInBytes is greater than size of available device memory.

Suppose that the user needs to perform 253 linear solvers and available device memory is 2GB. if cusolverSp[S|D|C|Z]csrqrsvBatched can only afford batchSize 100, the user has to call cusolverSp[S|D|C|Z]csrqrsvBatched three times to finish all. The user calls cusolverSp[S|D|C|Z]csrqrBufferInfoBatched with batchSize 100. The opaque info would remember this batchSize and any subsequent call of cusolverSp[S|D|C|Z]csrqrsvBatched cannot exceed this value. In this example, the first two calls of cusolverSp[S|D|C|Z]csrqrsvBatched will use batchSize 100, and last call of cusolverSp[S|D|C|Z]csrqrsvBatched will use batchSize 53.

Example: suppose that A0, A1, .., A9 have the same sparsity pattern, the following code solves 10 linear systems $A_j{x}_j=b_j\mathrm{,\; j\; =\; 0,2,...,\; 9}$ by batched sparse QR.

// Suppose that A0, A1, .., A9 are m x m sparse matrix represented by CSR format, 
// Each matrix Aj has nonzero nnzA, and shares the same csrRowPtrA and csrColIndA.
// csrValA is aggregation of A0, A1, ..., A9.
int m ; // number of rows and columns of each Aj 
int nnzA ; // number of nonzeros of each Aj
int *csrRowPtrA ; // each Aj has the same csrRowPtrA 
int *csrColIndA ; // each Aj has the same csrColIndA
double *csrValA ; // aggregation of A0,A1,...,A9
cont int batchSize = 10; // 10 linear systems 

cusolverSpHandle_t handle; // handle to cusolver library
csrqrInfo_t info = NULL;
cusparseMatDescr_t descrA = NULL;
void *pBuffer = NULL; // working space for numerical factorization
 
// step 1: create a descriptor
cusparseCreateMatDescr(&descrA);
cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ONE); // A is base-1
cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL); // A is a general matrix

// step 2: create empty info structure
cusolverSpCreateCsrqrInfo(&info);

// step 3: symbolic analysis
cusolverSpXcsrqrAnalysisBatched(
    handle, m, m, nnzA,
    descrA, csrRowPtrA, csrColIndA, info);

// step 4: allocate working space for Aj*xj=bj
cusolverSpDcsrqrBufferInfoBatched(
    handle, m, m, nnzA,
    descrA,
    csrValA, csrRowPtrA, csrColIndA,
    batchSize, 
    info,
    &internalDataInBytes,
    &workspaceInBytes);

cudaMalloc(&pBuffer, workspaceInBytes);

// step 5: solve Aj*xj = bj
cusolverSpDcsrqrsvBatched(
    handle, m, m, nnzA,
    descrA, csrValA, csrRowPtrA, csrColIndA,
    b,
    x,
    batchSize,
    info,
    pBuffer);
 
// step 7: destroy info
cusolverSpDestroyCsrqrInfo(info);

Please refer to Appendix B for detailed examples.

Remark 1: only GPU (device) path is provided.