1. Introduction

The cuSolver library is a high-level package based on the cuBLAS and cuSPARSE libraries. It combines three separate libraries under a single umbrella, each of which can be used independently or in concert with other toolkit libraries.

The intent of cuSolver is to provide useful LAPACK-like features, such as common matrix factorization and triangular solve routines for dense matrices, a sparse least-squares solver and an eigenvalue solver. In addition cuSolver provides a new refactorization library useful for solving sequences of matrices with a shared sparsity pattern.

The first part of cuSolver is called cuSolverDN, and deals with dense matrix factorization and solve routines such as LU, QR, SVD and LDLT, as well as useful utilities such as matrix and vector permutations.

Next, cuSolverSP provides a new set of sparse routines based on a sparse QR factorization. Not all matrices have a good sparsity pattern for parallelism in factorization, so the cuSolverSP library also provides a CPU path to handle those sequential-like matrices. For those matrices with abundant parallelism, the GPU path will deliver higher performance. The library is designed to be called from C and C++.

The final part is cuSolverRF, a sparse re-factorization package that can provide very good performance when solving a sequence of matrices where only the coefficients are changed but the sparsity pattern remains the same.

The GPU path of the cuSolver library assumes data is already in the device memory. It is the responsibility of the developer to allocate memory and to copy data between GPU memory and CPU memory using standard CUDA runtime API routines, such as cudaMalloc(), cudaFree(), cudaMemcpy(), and cudaMemcpyAsync().

Note: The cuSolver library requires hardware with a CUDA compute capability (CC) of at least 2.0 or higher. Please see the NVIDIA CUDA C Programming Guide, Appendix A for a list of the compute capabilities corresponding to all NVIDIA GPUs.

cuSolverDN: Dense LAPACK

The cuSolverDN library was designed to solve dense linear systems of the form

A x = b

where the coefficient matrix A R nxn , right-hand-side vector b R n and solution vector x R n

The cuSolverDN library provides QR factorization and LU with partial pivoting to handle a general matrix A, which may be non-symmetric. Cholesky factorization is also provided for symmetric/Hermitian matrices. For symmetric indefinite matrices, we provide Bunch-Kaufman (LDL) factorization.

The cuSolverDN library also provides a helpful bidiagonalization routine and singular value decomposition (SVD).

The cuSolverDN library targets computationally-intensive and popular routines in LAPACK, and provides an API compatible with LAPACK. The user can accelerate these time-consuming routines with cuSolverDN and keep others in LAPACK without a major change to existing code.

cuSolverSP: Sparse LAPACK

The cuSolverSP library was mainly designed to a solve sparse linear system

A x = b

and the least-squares problem

x = argmin || A * z - b ||

where sparse matrix A R mxn , right-hand-side vector b R m and solution vector x R n . For a linear system, we require m=n.

The core algorithm is based on sparse QR factorization. The matrix A is accepted in CSR format. If matrix A is symmetric/Hermitian, the user has to provide a full matrix, ie fill missing lower or upper part.

If matrix A is symmetric positive definite and the user only needs to solve A x = b , Cholesky factorization can work and the user only needs to provide the lower triangular part of A.

On top of the linear and least-squares solvers, the cuSolverSP library provides a simple eigenvalue solver based on shift-inverse power method, and a function to count the number of eigenvalues contained in a box in the complex plane.

cuSolverRF: Refactorization

The cuSolverRF library was designed to accelerate solution of sets of linear systems by fast re-factorization when given new coefficients in the same sparsity pattern

A i x i = f i

where a sequence of coefficient matrices A i R nxn , right-hand-sides f i R n and solutions x i R n are given for i=1,...,k.

The cuSolverRF library is applicable when the sparsity pattern of the coefficient matrices A i as well as the reordering to minimize fill-in and the pivoting used during the LU factorization remain the same across these linear systems. In that case, the first linear system (i=1) requires a full LU factorization, while the subsequent linear systems (i=2,...,k) require only the LU re-factorization. The later can be performed using the cuSolverRF library.

Notice that because the sparsity pattern of the coefficient matrices, the reordering and pivoting remain the same, the sparsity pattern of the resulting triangular factors L i and U i also remains the same. Therefore, the real difference between the full LU factorization and LU re-factorization is that the required memory is known ahead of time.

1.4. Naming Conventions

The cuSolverDN library functions are available for data types float, double, cuComplex, and cuDoubleComplex. The naming convention is as follows:

cusolverDn<t><operation>

where <t> can be S, D, C, Z, or X, corresponding to the data types float, double, cuComplex, cuDoubleComplex, and the generic type, respectively. <operation> can be Cholesky factorization (potrf), LU with partial pivoting (getrf), QR factorization (geqrf) and Bunch-Kaufman factorization (sytrf).

The cuSolverSP library functions are available for data types float, double, cuComplex, and cuDoubleComplex. The naming convention is as follows:

cusolverSp[Host]<t>[<matrix data format>]<operation>[<output matrix data format>]<based on>

where cuSolverSp is the GPU path and cusolverSpHost is the corresponding CPU path. <t> can be S, D, C, Z, or X, corresponding to the data types float, double, cuComplex, cuDoubleComplex, and the generic type, respectively.

The <matrix data format> is csr, compressed sparse row format.

The <operation> can be ls, lsq, eig, eigs, corresponding to linear solver, least-square solver, eigenvalue solver and number of eigenvalues in a box, respectively.

The <output matrix data format> can be v or m, corresponding to a vector or a matrix.

<based on> describes which algorithm is used. For example, qr (sparse QR factorization) is used in linear solver and least-square solver.

All of the functions have the return type cusolverStatus_t and are explained in more detail in the chapters that follow.

cuSolverSP API
routine data format operation output format based on
csrlsvlu csr linear solver (ls) vector (v) LU (lu) with partial pivoting
csrlsvqr csr linear solver (ls) vector (v) QR factorization (qr)
csrlsvchol csr linear solver (ls) vector (v) Cholesky factorization (chol)
csrlsqvqr csr least-square solver (lsq) vector (v) QR factorization (qr)
csreigvsi csr eigenvalue solver (eig) vector (v) shift-inverse
csreigs csr number of eigenvalues in a box (eigs)
csrsymrcm csr Symmetric Reverse Cuthill-McKee (symrcm)

The cuSolverRF library routines are available for data type double. Most of the routines follow the naming convention:

cusolverRf_<operation>_[[Host]](...)

where the trailing optional Host qualifier indicates the data is accessed on the host versus on the device, which is the default. The <operation> can be Setup, Analyze, Refactor, Solve, ResetValues, AccessBundledFactors and ExtractSplitFactors.

Finally, the return type of the cuSolverRF library routines is cusolverStatus_t.

1.5. Asynchronous Execution

The cuSolver library functions prefer to keep asynchronous execution as much as possible. Developers can always use the cudaDeviceSynchronize() function to ensure that the execution of a particular cuSolver library routine has completed.

A developer can also use the cudaMemcpy() routine to copy data from the device to the host and vice versa, using the cudaMemcpyDeviceToHost and cudaMemcpyHostToDevice parameters, respectively. In this case there is no need to add a call to cudaDeviceSynchronize() because the call to cudaMemcpy() with the above parameters is blocking and completes only when the results are ready on the host.

1.6. Library Property

The libraryPropertyType data type is an enumeration of library property types. (ie. CUDA version X.Y.Z would yield MAJOR_VERSION=X, MINOR_VERSION=Y, PATCH_LEVEL=Z)

typedef enum libraryPropertyType_t
{
        MAJOR_VERSION,
        MINOR_VERSION,
        PATCH_LEVEL
} libraryPropertyType;

The following code can show the version of cusolver library.

    int major=-1,minor=-1,patch=-1;
    cusolverGetProperty(MAJOR_VERSION, &major);
    cusolverGetProperty(MINOR_VERSION, &minor);
    cusolverGetProperty(PATCH_LEVEL, &patch);
    printf("CUSOLVER Version (Major,Minor,PatchLevel): %d.%d.%d\n", major,minor,patch);

2. Using the cuSolver API

This chapter describes how to use the cuSolver library API. It is not a reference for the cuSolver API data types and functions; that is provided in subsequent chapters.

2.1. Thread Safety

The library is thread safe and its functions can be called from multiple host threads.

2.2. Scalar Parameters

In the cuSolver API, the scalar parameters can be passed by reference on the host.

2.3. Parallelism with Streams

If the application performs several small independent computations, or if it makes data transfers in parallel with the computation, CUDA streams can be used to overlap these tasks.

The application can conceptually associate a stream with each task. To achieve the overlap of computation between the tasks, the developer should create CUDA streams using the function cudaStreamCreate() and set the stream to be used by each individual cuSolver library routine by calling for example cusolverDnSetStream() just before calling the actual cuSolverDN routine. Then, computations performed in separate streams would be overlapped automatically on the GPU, when possible. This approach is especially useful when the computation performed by a single task is relatively small and is not enough to fill the GPU with work, or when there is a data transfer that can be performed in parallel with the computation.

3. cuSolver Types Reference

3.1. cuSolverDN Types

The float, double, cuComplex, and cuDoubleComplex data types are supported. The first two are standard C data types, while the last two are exported from cuComplex.h. In addition, cuSolverDN uses some familiar types from cuBlas.

3.1.1. cusolverDnHandle_t

This is a pointer type to an opaque cuSolverDN context, which the user must initialize by calling cusolverDnCreate() prior to calling any other library function. An un-initialized Handle object will lead to unexpected behavior, including crashes of cuSolverDN. The handle created and returned by cusolverDnCreate() must be passed to every cuSolverDN function.

3.1.2. cublasFillMode_t

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.

Value Meaning
CUBLAS_FILL_MODE_LOWER the lower part of the matrix is filled
CUBLAS_FILL_MODE_UPPER the upper part of the matrix is filled

3.1.3. cublasOperation_t

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.

Value Meaning
CUBLAS_OP_N the non-transpose operation is selected
CUBLAS_OP_T the transpose operation is selected
CUBLAS_OP_C the conjugate transpose operation is selected

3.1.4. cusolverEigType_t

The cusolverEigType_t type indicates which type of eigenvalue solver is. Its values correspond to Fortran integer 1 (A*x = lambda*B*x), 2 (A*B*x = lambda*x), 3 (B*A*x = lambda*x), used as parameters to legacy LAPACK implementations.

Value Meaning
CUSOLVER_EIG_TYPE_1 A*x = lambda*B*x
CUSOLVER_EIG_TYPE_2 A*B*x = lambda*x
CUSOLVER_EIG_TYPE_3 B*A*x = lambda*x

3.1.5. cusolverEigMode_t

The cusolverEigMode_t type indicates whether or not eigenvectors are computed. Its values correspond to Fortran character 'N' (only eigenvalues are computed), 'V' (both eigenvalues and eigenvectors are computed) used as parameters to legacy LAPACK implementations.

Value Meaning
CUSOLVER_EIG_MODE_NOVECTOR only eigenvalues are computed
CUSOLVER_EIG_MODE_VECTOR both eigenvalues and eigenvectors are computed

3.1.6. cusolverStatus_t

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

3.2. cuSolverSP Types

The float, double, cuComplex, and cuDoubleComplex data types are supported. The first two are standard C data types, while the last two are exported from cuComplex.h.

3.2.1. cusolverSpHandle_t

This is a pointer type to an opaque cuSolverSP context, which the user must initialize by calling cusolverSpCreate() prior to calling any other library function. An un-initialized Handle object will lead to unexpected behavior, including crashes of cuSolverSP. The handle created and returned by cusolverSpCreate() must be passed to every cuSolverSP function.

3.2.2. cusparseMatDescr_t

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.

3.2.3. cusolverStatus_t

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

CUSOLVER_STATUS_SUCCESS

The operation completed successfully.

CUSOLVER_STATUS_NOT_INITIALIZED

The cuSolver library was not initialized. This is usually caused by the lack of a prior call, an error in the CUDA Runtime API called by the cuSolver routine, or an error in the hardware setup.

To correct: call cusolverCreate() prior to the function call; and check that the hardware, an appropriate version of the driver, and the cuSolver library are correctly installed.

CUSOLVER_STATUS_ALLOC_FAILED

Resource allocation failed inside the cuSolver library. This is usually caused by a cudaMalloc() failure.

To correct: prior to the function call, deallocate previously allocated memory as much as possible.

CUSOLVER_STATUS_INVALID_VALUE

An unsupported value or parameter was passed to the function (a negative vector size, for example).

To correct: ensure that all the parameters being passed have valid values.

CUSOLVER_STATUS_ARCH_MISMATCH

The function requires a feature absent from the device architecture; usually caused by the lack of support for atomic operations or double precision.

To correct: compile and run the application on a device with compute capability 2.0 or above.

CUSOLVER_STATUS_EXECUTION_FAILED

The GPU program failed to execute. This is often caused by a launch failure of the kernel on the GPU, which can be caused by multiple reasons.

To correct: check that the hardware, an appropriate version of the driver, and the cuSolver library are correctly installed.

CUSOLVER_STATUS_INTERNAL_ERROR

An internal cuSolver operation failed. This error is usually caused by a cudaMemcpyAsync() failure.

To correct: check that the hardware, an appropriate version of the driver, and the cuSolver library are correctly installed. Also, check that the memory passed as a parameter to the routine is not being deallocated prior to the routine’s completion.

CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED

The matrix type is not supported by this function. This is usually caused by passing an invalid matrix descriptor to the function.

To correct: check that the fields in descrA were set correctly.

3.3. cuSolverRF Types

cuSolverRF only supports double.

cusolverRfHandle_t

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.

cusolverRfMatrixFormat_t

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

Value Meaning
CUSOLVER_MATRIX_FORMAT_CSR matrix format CSR is assumed. (default)
CUSOLVER_MATRIX_FORMAT_CSC matrix format CSC is assumed.

cusolverRfNumericBoostReport_t

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.

Value Meaning
CUSOLVER_NUMERIC_BOOST_NOT_USED numeric boosting not used. (default)
CUSOLVER_NUMERIC_BOOST_USED numeric boosting used.

cusolverRfResetValuesFastMode_t

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.

Value Meaning
CUSOLVER_RESET_VALUES_FAST_MODE_OFF fast mode disabled. (default)
CUSOLVER_RESET_VALUES_FAST_MODE_ON fast mode enabled.

cusolverRfFactorization_t

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

Value Meaning
CUSOLVER_FACTORIZATION_ALG0 algorithm 0. (default)
CUSOLVER_FACTORIZATION_ALG1 algorithm 1.
CUSOLVER_FACTORIZATION_ALG2 algorithm 2. Domino-based scheme.

cusolverRfTriangularSolve_t

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

Value Meaning
CUSOLVER_TRIANGULAR_SOLVE_ALG0 algorithm 0.
CUSOLVER_TRIANGULAR_SOLVE_ALG1 algorithm 1. (default)
CUSOLVER_TRIANGULAR_SOLVE_ALG2 algorithm 2. Domino-based scheme.
CUSOLVER_TRIANGULAR_SOLVE_ALG3 algorithm 3. Domino-based scheme.

cusolverRfUnitDiagonal_t

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

Value Meaning
CUSOLVER_UNIT_DIAGONAL_STORED_L unit diagonal is stored in lower triangular factor. (default)
CUSOLVER_UNIT_DIAGONAL_STORED_U unit diagonal is stored in upper triangular factor.
CUSOLVER_UNIT_DIAGONAL_ASSUMED_L unit diagonal is assumed in lower triangular factor.
CUSOLVER_UNIT_DIAGONAL_ASSUMED_U unit diagonal is assumed in upper triangular factor.

cusolverStatus_t

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.

4. cuSolver Formats Reference

4.1. Index Base Format

The CSR or CSC format requires either zero-based or one-based index for a sparse matrix A. The GLU library supports only zero-based indexing. Otherwise, both one-based and zero-based indexing are supported in cuSolver.

4.2. Vector (Dense) Format

The vectors are assumed to be stored linearly in memory. For example, the vector

x = x 1 x 2 x n

is represented as

x 1 x 2 x n

4.3. Matrix (Dense) Format

The dense matrices are assumed to be stored in column-major order in memory. The sub-matrix can be accessed using the leading dimension of the original matrix. For examle, the m*n (sub-)matrix

a 1 , 1 a 1 , n a 2 , 1 a 2 , n a m , 1 a m , n

is represented as

a 1 , 1 a 1 , n a 2 , 1 a 2 , n a m , 1 a m , n a lda , 1 a lda , n

with its elements arranged linearly in memory as

a 1 , 1 a 2 , 1 a m , 1 a lda , 1 a 1 , n a 2 , n a m , n a lda , n

where ldam is the leading dimension of A.

4.4. Matrix (CSR) Format

In CSR format the matrix is represented by the following parameters

parameter type size Meaning
n (int) the number of rows (and columns) in the matrix.
nnz (int) the number of non-zero elements in the matrix.
csrRowPtr (int *) n+1 the array of offsets corresponding to the start of each row in the arrays csrColInd and csrVal. This array has also an extra entry at the end that stores the number of non-zero elements in the matrix.
csrColInd (int *) nnz the array of column indices corresponding to the non-zero elements in the matrix. It is assumed that this array is sorted by row and by column within each row.
csrVal (S|D|C|Z)* nnz the array of values corresponding to the non-zero elements in the matrix. It is assumed that this array is sorted by row and by column within each row.

Note that in our CSR format sparse matrices are assumed to be stored in row-major order, in other words, the index arrays are first sorted by row indices and then within each row by column indices. Also it is assumed that each pair of row and column indices appears only once.

For example, the 4x4 matrix

A = 1.0 3.0 0.0 0.0 0.0 4.0 6.0 0.0 2.0 5.0 7.0 8.0 0.0 0.0 0.0 9.0

is represented as

csrRowPtr = 0 2 4 8 9

csrColInd = 0 1 1 2 0 1 2 3 3

csrVal = 1.0 3.0 4.0 6.0 2.0 5.0 7.0 8.0 9.0

4.5. Matrix (CSC) Format

In CSC format the matrix is represented by the following parameters

parameter type size Meaning
n (int) the number of rows (and columns) in the matrix.
nnz (int) the number of non-zero elements in the matrix.
cscColPtr (int *) n+1 the array of offsets corresponding to the start of each column in the arrays cscRowInd and cscVal. This array has also an extra entry at the end that stores the number of non-zero elements in the matrix.
cscRowInd (int *) nnz the array of row indices corresponding to the non-zero elements in the matrix. It is assumed that this array is sorted by column and by row within each column.
cscVal (S|D|C|Z)* nnz the array of values corresponding to the non-zero elements in the matrix. It is assumed that this array is sorted by column and by row within each column.

Note that in our CSC format sparse matrices are assumed to be stored in column-major order, in other words, the index arrays are first sorted by column indices and then within each column by row indices. Also it is assumed that each pair of row and column indices appears only once.

For example, the 4x4 matrix

A = 1.0 3.0 0.0 0.0 0.0 4.0 6.0 0.0 2.0 5.0 7.0 8.0 0.0 0.0 0.0 9.0

is represented as

cscColPtr = 0 2 5 7 9

cscRowInd = 0 2 0 1 2 1 2 2 3

cscVal = 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

cuSolverDN: dense LAPACK Function Reference

This chapter describes the API of cuSolverDN, which provides a subset of dense LAPACK functions.

cuSolverDN Helper Function Reference

The cuSolverDN helper functions are described in this section.

5.1.1. cusolverDnCreate()

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.

parameter Memory In/out Meaning
handle host output the pointer to the handle to the cuSolverDN context.
Status Returned
CUSOLVER_STATUS_SUCCESS the initialization succeeded.
CUSOLVER_STATUS_NOT_INITIALIZED the CUDA Runtime initialization failed.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.

5.1.2. cusolverDnDestroy()

cusolverStatus_t 
cusolverDnDestroy(cusolverDnHandle_t handle);

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

parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
Status Returned
CUSOLVER_STATUS_SUCCESS the shutdown succeeded.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.

cusolverDnSetStream()

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.

parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
streamId host input the stream to be used by the library.
Status Returned
CUSOLVER_STATUS_SUCCESS the stream was set successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.

cusolverDnGetStream()

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.

parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
streamId host output the stream to be used by the library.
Status Returned
CUSOLVER_STATUS_SUCCESS the stream was set successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.

Dense Linear Solver Reference

This chapter describes linear solver API of cuSolverDN, including Cholesky factorization, LU with partial pivoting, QR factorization and Bunch-Kaufman (LDLT) factorization.

cusolverDn<t>potrf()

These helper functions calculate the necessary size of work buffers.
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);
The S and D data types are real valued single and double precision, respectively.
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 );
The C and Z data types are complex valued single and double precision, respectively.
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.

A = L * L H

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.

A = U * U H

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.

API of potrf
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
uplo host input indicates if matrix A lower or upper part is stored, the other part is not referenced.
n host input number of rows and columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n).
lda host input leading dimension of two-dimensional array used to store matrix A.
Workspace device in/out working space, <type> array of size Lwork.
Lwork host input size of Workspace, returned by potrf_bufferSize.
devInfo device output if devInfo = 0, the Cholesky factorization is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo = i, the leading minor of order i is not positive definite.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0 or lda<max(1,n)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>potrs()

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

A * X = B

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.

API of potrs
parameter Memory In/out Meaning
handle host input handle to the cuSolveDN library context.
uplo host input indicates if matrix A lower or upper part is stored, the other part is not referenced.
n host input number of rows and columns of matrix A.
nrhs host input number of columns of matrix X and B.
A device input <type> array of dimension lda * n with lda is not less than max(1,n). A is either lower cholesky factor L or upper Cholesky factor U.
lda host input leading dimension of two-dimensional array used to store matrix A.
B device in/out <type> array of dimension ldb * nrhs. ldb is not less than max(1,n). As an input, B is right hand side matrix. As an output, B is the solution matrix.
devInfo device output if devInfo = 0, the Cholesky factorization is successful. if devInfo = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0, nrhs<0, lda<max(1,n) or ldb<max(1,n)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>getrf()

These helper functions calculate the size of work buffers needed.
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 );
The S and D data types are real single and double precision, respectively.
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 );
The C and Z data types are complex valued single and double precision, respectively.
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

P * A = L * U

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).

API of getrf
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
m host input number of rows of matrix A.
n host input number of columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,m).
lda host input leading dimension of two-dimensional array used to store matrix A.
Workspace device in/out working space, <type> array of size Lwork.
devIpiv device output array of size at least min(m,n), containing pivot indices.
devInfo device output if devInfo = 0, the LU factorization is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo = i, the U(i,i) = 0.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or lda<max(1,m)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>getrs()

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

op(A) * X = B

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

op ( A ) = A if  trans == CUBLAS_OP_N A T if  trans == CUBLAS_OP_T A H if  trans == CUBLAS_OP_C

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.

parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
trans host input operation op(A) that is non- or (conj.) transpose.
n host input number of rows and columns of matrix A.
nrhs host input number of right-hand sides.
A device input <type> array of dimension lda * n with lda is not less than max(1,n).
lda host input leading dimension of two-dimensional array used to store matrix A.
devIpiv device input array of size at least n, containing pivot indices.
B device output <type> array of dimension ldb * nrhs with ldb is not less than max(1,n).
ldb host input leading dimension of two-dimensional array used to store matrix B.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0 or lda<max(1,n) or ldb<max(1,n)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>geqrf()

These helper functions calculate the size of work buffers needed.
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 );
The S and D data types are real valued single and double precision, respectively.
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 );
The C and Z data types are complex valued single and double precision, respectively.
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

A = Q * R

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

I - 2 * v * v H = I - τ * q * q H

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

API of geqrf
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
m host input number of rows of matrix A.
n host input number of columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,m).
lda host input leading dimension of two-dimensional array used to store matrix A.
TAU device output <type> array of dimension at least min(m,n).
Workspace device in/out working space, <type> array of size Lwork.
Lwork host input size of working array Workspace.
devInfo device output if info = 0, the LU factorization is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or lda<max(1,m)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>ormqr()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSormqr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasOperation_t trans,
    int m,
    int n,
    int k,
    const float *A,
    int lda,
    const float *C,
    int ldc,
    int *lwork);

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

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

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

The S and D data types are real valued single and double precision, respectively.
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);
The C and Z data types are complex valued single and double precision, respectively.
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

C = op ( Q ) * C if  side == CUBLAS_SIDE_LEFT C * op ( Q ) if  side == CUBLAS_SIDE_RIGHT

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

op ( Q ) = Q if  transa == CUBLAS_OP_N Q T if  transa == CUBLAS_OP_T Q H if  transa == CUBLAS_OP_C

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() or ormqr_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.

API of ormqr
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
side host input indicates if matrix Q is on the left or right of C.
trans host input operation op(Q) that is non- or (conj.) transpose.
m host input number of rows of matrix A.
n host input number of columns of matrix A.
k host input number of elementary relfections.
A device in/out <type> array of dimension lda * k with lda is not less than max(1,m). The matrix A is from geqrf, so i-th column contains elementary reflection vector.
lda host input leading dimension of two-dimensional array used to store matrix A. if side is CUBLAS_SIDE_LEFT, lda >= max(1,m); if side is CUBLAS_SIDE_RIGHT, lda >= max(1,n).
tau device output <type> array of dimension at least min(m,n). The vector tau is from geqrf, so tau(i) is the scalar of i-th elementary reflection vector.
C device in/out <type> array of size ldc * n. On exit, C is overwritten by op(Q)*C.
ldc host input leading dimension of two-dimensional array of matrix C. ldc >= max(1,m).
work device in/out working space, <type> array of size lwork.
lwork host input size of working array work.
devInfo device output if info = 0, the ormqr is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or wrong lda or ldc).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>orgqr()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSorgqr_bufferSize(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    const float *A,
    int lda,
    int *lwork);

cusolverStatus_t 
cusolverDnDorgqr_bufferSize(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    const double *A,
    int lda,
    int *lwork);

cusolverStatus_t 
cusolverDnCungqr_bufferSize(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    const cuComplex *A,
    int lda,
    int *lwork);

cusolverStatus_t 
cusolverDnZungqr_bufferSize(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    const cuDoubleComplex *A,
    int lda,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSorgqr(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    float *A,
    int lda,
    const float *tau,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDorgqr(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    double *A,
    int lda,
    const double *tau,
    double *work,
    int lwork,
    int *devInfo);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCungqr(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    cuComplex *A,
    int lda,
    const cuComplex *tau,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t
cusolverDnZungqr(
    cusolverDnHandle_t handle,
    int m,
    int n,
    int k,
    cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function overwrites m×n matrix A by

Q = H(1) * H(2) * ... * H(k)

where Q is a unitary matrix formed by a sequence of elementary reflection vectors stored in 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 orgqr_bufferSize().

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

The user can combine geqrf, orgqr to complete orthogonalization. Please refer to appendix C.2.

API of ormqr
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
m host input number of rows of matrix Q. m >= 0;
n host input number of columns of matrix Q. m >= n >= 0;
k host input number of elementary relfections whose product defines the matrix Q. n >= k >= 0;
A device in/out <type> array of dimension lda * n with lda is not less than max(1,m). i-th column of A contains elementary reflection vector.
lda host input leading dimension of two-dimensional array used to store matrix A. lda >= max(1,m).
tau device output <type> array of dimension k. tau(i) is the scalar of i-th elementary reflection vector.
work device in/out working space, <type> array of size lwork.
lwork host input size of working array work.
devInfo device output if info = 0, the orgqr is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n,k<0, n>m, k>n or lda<m).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>sytrf()

These helper functions calculate the size of the needed buffers.
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 );
The S and D data types are real valued single and double precision, respectively.
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 );
The C and Z data types are complex valued single and double precision, respectively.
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.

P * A * P T = L * D * L T

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.

P * A * P T = U * D * U T

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.

API of sytrf
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
uplo host input indicates if matrix A lower or upper part is stored, the other part is not referenced.
n host input number of rows and columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n).
lda host input leading dimension of two-dimensional array used to store matrix A.
ipiv device output array of size at least n, containing pivot indices.
work device in/out working space, <type> array of size lwork.
lwork host input size of working space work.
devInfo device output if devInfo = 0, the LU factorization is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo = i, the D(i,i) = 0.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0 or lda<max(1,n)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

Dense Eigenvalue Solver Reference

This chapter describes eigenvalue solver API of cuSolverDN, including bidiagonalization and SVD.

cusolverDn<t>gebrd()

These helper functions calculate the size of work buffers needed.
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 );
The S and D data types are real valued single and double precision, respectively.
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 );
The C and Z data types are complex valued single and double precision, respectively.
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 m×n matrix A to a real 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.

API of gebrd
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
m host input number of rows of matrix A.
n host input number of columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n).
lda host input leading dimension of two-dimensional array used to store matrix A.
D device output real array of dimension min(m,n). The diagonal elements of the bidiagonal matrix B: D(i) = A(i,i).
E device output real array of dimension min(m,n). The off-diagonal elements of the bidiagonal matrix B: if m>=n, E(i) = A(i,i+1) for i = 1,2,...,n-1; if m<n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
TAUQ device output <type> array of dimension min(m,n). The scalar factors of the elementary reflectors which represent the orthogonal matrix Q.
TAUP device output <type> array of dimension min(m,n). The scalar factors of the elementary reflectors which represent the orthogonal matrix P.
Work device in/out working space, <type> array of size Lwork.
Lwork host input size of Work, returned by gebrd_bufferSize.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0, or lda<max(1,m)).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>orgbr()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSorgbr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    const float *A,
    int lda,
    const float *tau,
    int *lwork);

cusolverStatus_t
cusolverDnDorgbr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    const double *A,
    int lda,
    const double *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnCungbr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    const cuComplex *A,
    int lda,
    const cuComplex *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnZungbr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    const cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSorgbr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    float *A,
    int lda,
    const float *tau,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDorgbr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    double *A,
    int lda,
    const double *tau,
    double *work,
    int lwork,
    int *devInfo);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCungbr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    cuComplex *A,
    int lda,
    const cuComplex *tau,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnZungbr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    int m,
    int n,
    int k,
    cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function generates one of the unitary matrices Q or P**H determined by gebrd when reducing a matrix A to bidiagonal form: Q H * A * P = B

Q and P**H are defined as products of elementary reflectors H(i) or G(i) respectively.

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 orgbr_bufferSize().

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

API of orgbr
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
side host input if side = CUBLAS_SIDE_LEFT, generate Q. if side = CUBLAS_SIDE_RIGHT, generate P**T.
m host input number of rows of matrix Q or P**T.
n host input if side = CUBLAS_SIDE_LEFT, m>= n>= min(m,k). if side = CUBLAS_SIDE_RIGHT, n>= m>= min(n,k).
k host input if side = CUBLAS_SIDE_LEFT, the number of columns in the original m-by-k matrix reduced by gebrd. if side = CUBLAS_SIDE_RIGHT, the number of rows in the original k-by-n matrix reduced by gebrd.
A device in/out <type> array of dimension lda * n On entry, the vectors which define the elementary reflectors, as returned by gebrd. On exit, the m-by-n matrix Q or P**T.
lda host input leading dimension of two-dimensional array used to store matrix A. lda >= max(1,m);
tau device output <type> array of dimension min(m,k) if side is CUBLAS_SIDE_LEFT; of dimension min(n,k) if side is CUBLAS_SIDE_RIGHT; tau(i) must contain the scalar factor of the elementary reflector H(i) or G(i), which determines Q or P**T, as returned by gebrd in its array argument TAUQ or TAUP.
work device in/out working space, <type> array of size lwork.
lwork host input size of working array work.
devInfo device output if info = 0, the ormqr is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or wrong lda ).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>sytrd()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSsytrd_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const float *A,
    int lda,
    const float *d,
    const float *e,
    const float *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnDsytrd_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const double *A,
    int lda,
    const double *d,
    const double *e,
    const double *tau,
    int *lwork);

cusolverStatus_t
cusolverDnChetrd_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const cuComplex *A,
    int lda,
    const float *d,
    const float *e,
    const cuComplex *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnZhetrd_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const cuDoubleComplex *A,
    int lda,
    const double *d,
    const double *e,
    const cuDoubleComplex *tau,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSsytrd(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    float *A,
    int lda,
    float *d,
    float *e,
    float *tau,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDsytrd(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    double *A,
    int lda,
    double *d,
    double *e,
    double *tau,
    double *work,
    int lwork,
    int *devInfo);
The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnChetrd(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    cuComplex *A,
    int lda,
    float *d,
    float *e,
    cuComplex *tau,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t CUDENSEAPI cusolverDnZhetrd(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    cuDoubleComplex *A,
    int lda,
    double *d,
    double *e,
    cuDoubleComplex *tau,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function reduces a general symmetric (Hermitian) n×n matrix A to real symmetric tridiagonal form T by an orthogonal transformation: Q H * A * Q = T

As an output, A contains T and householder reflection vectors. If uplo = CUBLAS_FILL_MODE_UPPER, the diagonal and first superdiagonal of A are overwritten by the corresponding elements of the tridiagonal matrix T, and the elements above the first superdiagonal, with the array tau, represent the orthogonal matrix Q as a product of elementary reflectors; If uplo = CUBLAS_FILL_MODE_LOWER, the diagonal and first subdiagonal of A are overwritten by the corresponding elements of the tridiagonal matrix T, and the elements below the first subdiagonal, with the array tau, represent the orthogonal matrix Q 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 sytrd_bufferSize().

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

API of sytrd
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
uplo host input specifies which part of A is stored. uplo = CUBLAS_FILL_MODE_LOWER: Lower triangle of A is stored. uplo = CUBLAS_FILL_MODE_UPPER: Upper triangle of A is stored.
n host input number of rows (columns) of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n). If uplo = CUBLAS_FILL_MODE_UPPER, the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If uplo = CUBLAS_FILL_MODE_LOWER, the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, A is overwritten by T and householder reflection vectors.
lda host input leading dimension of two-dimensional array used to store matrix A. lda >= max(1,n).
D device output real array of dimension n. The diagonal elements of the tridiagonal matrix T: D(i) = A(i,i).
E device output real array of dimension (n-1). The off-diagonal elements of the tridiagonal matrix T: if uplo = CUBLAS_FILL_MODE_UPPER, E(i) = A(i,i+1). if uplo = CUBLAS_FILL_MODE_LOWERE(i) = A(i+1,i).
tau device output <type> array of dimension (n-1). The scalar factors of the elementary reflectors which represent the orthogonal matrix Q.
work device in/out working space, <type> array of size lwork.
lwork host input size of work, returned by sytrd_bufferSize.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0, or lda<max(1,n), or uplo is not CUBLAS_FILL_MODE_LOWER or CUBLAS_FILL_MODE_UPPER).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>ormtr()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSormtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    const float *A,
    int lda,
    const float *tau,
    const float *C,
    int ldc,
    int *lwork);

cusolverStatus_t 
cusolverDnDormtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    const double *A,
    int lda,
    const double *tau,
    const double *C,
    int ldc,
    int *lwork);

cusolverStatus_t 
cusolverDnCunmtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    const cuComplex *A,
    int lda,
    const cuComplex *tau,
    const cuComplex *C,
    int ldc,
    int *lwork);

cusolverStatus_t 
cusolverDnZunmtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    const cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    const cuDoubleComplex *C,
    int ldc,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSormtr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    float *A,
    int lda,
    float *tau,
    float *C,
    int ldc,
    float *work,
    int lwork,
    int *info);

cusolverStatus_t 
cusolverDnDormtr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    double *A,
    int lda,
    double *tau,
    double *C,
    int ldc,
    double *work,
    int lwork,
    int *info);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCunmtr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    cuComplex *A,
    int lda,
    cuComplex *tau,
    cuComplex *C,
    int ldc,
    cuComplex *work,
    int lwork,
    int *info);

cusolverStatus_t 
cusolverDnZunmtr(
    cusolverDnHandle_t handle,
    cublasSideMode_t side,
    cublasFillMode_t uplo,
    cublasOperation_t trans,
    int m,
    int n,
    cuDoubleComplex *A,
    int lda,
    cuDoubleComplex *tau,
    cuDoubleComplex *C,
    int ldc,
    cuDoubleComplex *work,
    int lwork,
    int *info);

This function overwrites m×n matrix C by

C = op ( Q ) * C if  side == CUBLAS_SIDE_LEFT C * op ( Q ) if  side == CUBLAS_SIDE_RIGHT

where Q is a unitary matrix formed by a sequence of elementary reflection vectors from sytrd.

The operation on Q is defined by

op ( Q ) = Q if  transa == CUBLAS_OP_N Q T if  transa == CUBLAS_OP_T Q H if  transa == CUBLAS_OP_C

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 ormtr_bufferSize().

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

API of ormtr
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
side host input side = CUBLAS_SIDE_LEFT, apply Q or Q**T from the Left; side = CUBLAS_SIDE_RIGHT, apply Q or Q**T from the Right.
uplo host input uplo = CUBLAS_FILL_MODE_LOWER: Lower triangle of A contains elementary reflectors from sytrd. uplo = CUBLAS_FILL_MODE_UPPER: Upper triangle of A contains elementary reflectors from sytrd.
trans host input operation op(Q) that is non- or (conj.) transpose.
m host input number of rows of matrix C.
n host input number of columns of matrix C.
A device in/out <type> array of dimension lda * m if side = CUBLAS_SIDE_LEFT; lda * n if side = CUBLAS_SIDE_RIGHT. The matrix A from sytrd contains the elementary reflectors.
lda host input leading dimension of two-dimensional array used to store matrix A. if side is CUBLAS_SIDE_LEFT, lda >= max(1,m); if side is CUBLAS_SIDE_RIGHT, lda >= max(1,n).
tau device output <type> array of dimension (m-1) if side is CUBLAS_SIDE_LEFT; of dimension (n-1) if side is CUBLAS_SIDE_RIGHT; The vector tau is from sytrd, so tau(i) is the scalar of i-th elementary reflection vector.
C device in/out <type> array of size ldc * n. On exit, C is overwritten by op(Q)*C or C*op(Q).
ldc host input leading dimension of two-dimensional array of matrix C. ldc >= max(1,m).
work device in/out working space, <type> array of size lwork.
lwork host input size of working array work.
devInfo device output if info = 0, the ormqr is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or wrong lda or ldc).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>orgtr()

These helper functions calculate the size of work buffers needed.
cusolverStatus_t 
cusolverDnSorgtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const float *A,
    int lda,
    const float *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnDorgtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const double *A,
    int lda,
    const double *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnCungtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const cuComplex *A,
    int lda,
    const cuComplex *tau,
    int *lwork);

cusolverStatus_t 
cusolverDnZungtr_bufferSize(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    const cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSorgtr(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    float *A,
    int lda,
    const float *tau,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDorgtr(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    double *A,
    int lda,
    const double *tau,
    double *work,
    int lwork,
    int *devInfo);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCungtr(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    cuComplex *A,
    int lda,
    const cuComplex *tau,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnZungtr(
    cusolverDnHandle_t handle,
    cublasFillMode_t uplo,
    int n,
    cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *tau,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function generates a unitary matrix Q which is defined as the product of n-1 elementary reflectors of order n, as returned by sytrd:

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 orgtr_bufferSize().

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

API of orgtr
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
uplo host input uplo = CUBLAS_FILL_MODE_LOWER: Lower triangle of A contains elementary reflectors from sytrd. uplo = CUBLAS_FILL_MODE_UPPER: Upper triangle of A contains elementary reflectors from sytrd.
n host input number of rows (columns) of matrix Q.
A device in/out <type> array of dimension lda * n On entry, matrix A from sytrd contains the elementary reflectors. On exit, matrix A contains the n-by-n orthogonal matrix Q.
lda host input leading dimension of two-dimensional array used to store matrix A. lda >= max(1,n).
tau device output <type> array of dimension (n-1)tau(i) is the scalar of i-th elementary reflection vector.
work device in/out working space, <type> array of size lwork.
lwork host input size of working array work.
devInfo device output if info = 0, the orgtr is successful. if info = -i, the i-th parameter is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0 or wrong lda ).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>gesvd()

The helper functions below can calculate the sizes needed for pre-allocated buffer.
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 );

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSgesvd (
    cusolverDnHandle_t handle,
    signed char jobu,
    signed 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,
    signed char jobu,
    signed 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);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCgesvd (
    cusolverDnHandle_t handle,
    signed char jobu,
    signed 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,
    signed char jobu,
    signed 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

A = U * Σ * V H

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 did not converge to zero.

The rwork is real array of dimension (min(m,n)-1). If devInfo>0 and rwork is not nil, rwork contains the unconverged superdiagonal elements of an upper bidiagonal matrix. This is slightly different from LAPACK which puts unconverged superdiagonal elements in work if type is real; in rwork if type is complex. rwork can be a NULL pointer if the user does not want the information from supperdiagonal.

Appendix E.1 provides a simple example of gesvd.

Remark 1: gesvd only supports m>=n.

Remark 2: the routine returns V H , not V.

API of gesvd
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
jobu host input specifies options for computing all or part of the matrix U: = 'A': all m columns of U are returned in array U: = 'S': the first min(m,n) columns of U (the left singular vectors) are returned in the array U; = 'O': the first min(m,n) columns of U (the left singular vectors) are overwritten on the array A; = 'N': no columns of U (no left singular vectors) are computed.
jobvt host input specifies options for computing all or part of the matrix V**T: = 'A': all N rows of V**T are returned in the array VT; = 'S': the first min(m,n) rows of V**T (the right singular vectors) are returned in the array VT; = 'O': the first min(m,n) rows of V**T (the right singular vectors) are overwritten on the array A; = 'N': no rows of V**T (no right singular vectors) are computed.
m host input number of rows of matrix A.
n host input number of columns of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,m). On exit, the contents of A are destroyed.
lda host input leading dimension of two-dimensional array used to store matrix A.
S device output real array of dimension min(m,n). The singular values of A, sorted so that S(i) >= S(i+1).
U device output <type> array of dimension ldu * m with ldu is not less than max(1,m). U contains the m×m unitary matrix U.
ldu host input leading dimension of two-dimensional array used to store matrix U.
VT device output <type> array of dimension ldvt * n with ldvt is not less than max(1,n). VT contains the n×n unitary matrix V**T.
ldvt host input leading dimension of two-dimensional array used to store matrix Vt.
work device in/out working space, <type> array of size lwork.
lwork host input size of work, returned by gesvd_bufferSize.
rwork device input real array of dimension min(m,n)-1. It contains the unconverged superdiagonal elements of an upper bidiagonal matrix if devInfo > 0.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo > 0, devInfo indicates how many superdiagonals of an intermediate bidiagonal form did not converge to zero.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n<0 or lda<max(1,m) or ldu<max(1,m) or ldvt<max(1,n) ).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>syevd()

The helper functions below can calculate the sizes needed for pre-allocated buffer.
cusolverStatus_t 
cusolverDnSsyevd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const float *A,
    int lda,
    const float *W,
    int *lwork);

cusolverStatus_t 
cusolverDnDsyevd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const double *A,
    int lda,
    const double *W,
    int *lwork);

cusolverStatus_t 
cusolverDnCheevd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const cuComplex *A,
    int lda,
    const float *W,
    int *lwork);

cusolverStatus_t 
cusolverDnZheevd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const cuDoubleComplex *A,
    int lda,
    const double *W,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnSsyevd(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    float *A,
    int lda,
    float *W,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDsyevd(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    double *A,
    int lda,
    double *W,
    double *work,
    int lwork,
    int *devInfo);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnCheevd(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    cuComplex *A,
    int lda,
    float *W,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnZheevd(
    cusolverDnHandle_t handle,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    cuDoubleComplex *A,
    int lda,
    double *W,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function computes eigenvalues and eigenvectors of a symmetric (Hermitian) n×n matrix A. The standard symmetric eigenvalue problem is

A * V = V * Λ

where Λ is a real n×n diagonal matrix. V is an n×n unitary matrix. The diagonal elements of Λ are the eigenvalues of A in ascending order.

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 syevd_bufferSize().

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong. If devInfo = i (greater than zero), i off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

if jobz = CUSOLVER_EIG_MODE_VECTOR, A contains the orthonormal eigenvectors of the matrix A. The eigenvectors are computed by a divide and conquer algorithm.

Appendix D.1 provides a simple example of syevd.

API of syevd
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
jobz host input specifies options to either compute eigenvalue only or compute eigen-pair: jobz = CUSOLVER_EIG_MODE_NOVECTOR: Compute eigenvalues only; jobz = CUSOLVER_EIG_MODE_VECTOR: Compute eigenvalues and eigenvectors.
uplo host input specifies which part of A is stored. uplo = CUBLAS_FILL_MODE_LOWER: Lower triangle of A is stored. uplo = CUBLAS_FILL_MODE_UPPER: Upper triangle of A is stored.
n host input number of rows (or columns) of matrix A.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n). If uplo = CUBLAS_FILL_MODE_UPPER, the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A. If uplo = CUBLAS_FILL_MODE_LOWER, the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A. On exit, if jobz = 'V', and devInfo = 0, A contains the orthonormal eigenvectors of the matrix A. If jobz = 'N', the contents of A are destroyed.
lda host input leading dimension of two-dimensional array used to store matrix A.
W device output a real array of dimension n. The eigenvalue values of A, in ascending order ie, sorted so that W(i) <= W(i+1).
work device in/out working space, <type> array of size lwork.
Lwork host input size of work, returned by syevd_bufferSize.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo = i (> 0), devInfo indicates i off-diagonal elements of an intermediate tridiagonal form did not converge to zero;
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0, or lda<max(1,n), or jobz is not 'N' or 'V', or uplo is not CUBLAS_FILL_MODE_LOWER or CUBLAS_FILL_MODE_UPPER).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cusolverDn<t>sygvd()

The helper functions below can calculate the sizes needed for pre-allocated buffer.
cusolverStatus_t 
cusolverDnSsygvd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const float *A,
    int lda,
    const float *B,
    int ldb,
    const float *W,
    int *lwork);

cusolverStatus_t
cusolverDnDsygvd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const double *A,
    int lda,
    const double *B,
    int ldb,
    const double *W,
    int *lwork);

cusolverStatus_t 
cusolverDnChegvd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const cuComplex *A,
    int lda,
    const cuComplex *B,
    int ldb,
    const float *W,
    int *lwork);

cusolverStatus_t 
cusolverDnZhegvd_bufferSize(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    const cuDoubleComplex *A,
    int lda,
    const cuDoubleComplex *B,
    int ldb,
    const double *W,
    int *lwork);

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t
cusolverDnSsygvd(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    float *A,
    int lda,
    float *B,
    int ldb,
    float *W,
    float *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnDsygvd(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    double *A,
    int lda,
    double *B,
    int ldb,
    double *W,
    double *work,
    int lwork,
    int *devInfo);

The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverDnChegvd(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    cuComplex *A,
    int lda,
    cuComplex *B,
    int ldb,
    float *W,
    cuComplex *work,
    int lwork,
    int *devInfo);

cusolverStatus_t 
cusolverDnZhegvd(
    cusolverDnHandle_t handle,
    cusolverEigType_t itype,
    cusolverEigMode_t jobz,
    cublasFillMode_t uplo,
    int n,
    cuDoubleComplex *A,
    int lda,
    cuDoubleComplex *B,
    int ldb,
    double *W,
    cuDoubleComplex *work,
    int lwork,
    int *devInfo);

This function computes eigenvalues and eigenvectors of a symmetric (Hermitian) n×n matrix-pair (A,B). The generalized symmetric-definite eigenvalue problem is

eig(A,B) = A * V = B * V * Λ if  itype = CUSOLVER_EIG_TYPE_1 A * B * V = V * Λ if  itype = CUSOLVER_EIG_TYPE_2 B * A * V = V * Λ if  itype = CUSOLVER_EIG_TYPE_3

where the matrix B is positive definite. Λ is a real n×n diagonal matrix. The diagonal elements of Λ are the eigenvalues of (A, B) in ascending order. V is an n×n orthogonal matrix. The eigenvectors are normalized as follows:

V H * B * V = I if  itype = CUSOLVER_EIG_TYPE_1, CUSOLVER_EIG_TYPE_2 V H * inv(B) * V = I if  itype = CUSOLVER_EIG_TYPE_3

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 sygvd_bufferSize().

If output parameter devInfo = -i (less than zero), the i-th parameter is wrong. If devInfo = i (i > 0 and i<=n) and jobz = CUSOLVER_EIG_MODE_NOVECTOR, i off-diagonal elements of an intermediate tridiagonal form did not converge to zero. If devInfo = N + i (i > 0), then the leading minor of order i of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed.

if jobz = CUSOLVER_EIG_MODE_VECTOR, A contains the orthogonal eigenvectors of the matrix A. The eigenvectors are computed by divide and conquer algorithm.

Appendix D.2 provides a simple example of sygvd.

API of sygvd
parameter Memory In/out Meaning
handle host input handle to the cuSolverDN library context.
itype host input Specifies the problem type to be solved: itype=CUSOLVER_EIG_TYPE_1: A*x = (lambda)*B*x. itype=CUSOLVER_EIG_TYPE_2: A*B*x = (lambda)*x. itype=CUSOLVER_EIG_TYPE_3: B*A*x = (lambda)*x.
jobz host input specifies options to either compute eigenvalue only or compute eigen-pair: jobz = CUSOLVER_EIG_MODE_NOVECTOR: Compute eigenvalues only; jobz = CUSOLVER_EIG_MODE_VECTOR: Compute eigenvalues and eigenvectors.
uplo host input specifies which part of A and B are stored. uplo = CUBLAS_FILL_MODE_LOWER: Lower triangle of A and B are stored. uplo = CUBLAS_FILL_MODE_UPPER: Upper triangle of A and B are stored.
n host input number of rows (or columns) of matrix A and B.
A device in/out <type> array of dimension lda * n with lda is not less than max(1,n). If uplo = CUBLAS_FILL_MODE_UPPER, the leading n-by-n upper triangular part of A contains the upper triangular part of the matrix A. If uplo = CUBLAS_FILL_MODE_LOWER, the leading n-by-n lower triangular part of A contains the lower triangular part of the matrix A. On exit, if jobz = 'V', and devInfo = 0, A contains the orthonormal eigenvectors of the matrix A. If jobz = 'N', the contents of A are destroyed.
lda host input leading dimension of two-dimensional array used to store matrix A. lda is not less than max(1,n).
B device in/out <type> array of dimension ldb * n. If uplo = CUBLAS_FILL_MODE_UPPER, the leading n-by-n upper triangular part of B contains the upper triangular part of the matrix B. If uplo = CUBLAS_FILL_MODE_LOWER, the leading n-by-n lower triangular part of B contains the lower triangular part of the matrix B. On exit, if devInfo is less than n, B is overwritten by triangular factor U or L from the Cholesky factorization of B.
ldb host input leading dimension of two-dimensional array used to store matrix B. ldb is not less than max(1,n).
W device output a real array of dimension n. The eigenvalue values of A, sorted so that W(i) >= W(i+1).
work device in/out working space, <type> array of size lwork.
Lwork host input size of work, returned by sygvd_bufferSize.
devInfo device output if devInfo = 0, the operation is successful. if devInfo = -i, the i-th parameter is wrong. if devInfo = i (> 0), devInfo indicates either potrf or syevd is wrong.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n<0, or lda<max(1,n), or ldb<max(1,n), or itype is not 1, 2 or 3, or jobz is not 'N' or 'V', or uplo is not CUBLAS_FILL_MODE_LOWER or CUBLAS_FILL_MODE_UPPER).
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.

cuSolverSP: sparse LAPACK Function Reference

This chapter describes the API of cuSolverSP, which provides a subset of LAPACK funtions for sparse matrices in CSR or CSC format.

Helper Function Reference

cusolverSpCreate()

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.

Output
handle the pointer to the handle to the cuSolverSP context.
Status Returned
CUSOLVER_STATUS_SUCCESS the initialization succeeded.
CUSOLVER_STATUS_NOT_INITIALIZED the CUDA Runtime initialization failed.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.

cusolverSpDestroy()

cusolverStatus_t
cusolverSpDestroy(cusolverSpHandle_t handle)

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

Input
handle the handle to the cuSolverSP context.
Status Returned
CUSOLVER_STATUS_SUCCESS the shutdown succeeded.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.

cusolverSpSetStream()

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.

Input
handle the handle to the cuSolverSP context.
streamId the stream to be used by the library.
Status Returned
CUSOLVER_STATUS_SUCCESS the stream was set successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.

cusolverSpXcsrissym()


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.

Input
parameter MemorySpace description
handle host handle to the cuSolverSP library context.
m host number of rows and columns of matrix A.
nnzA host number of nonzeros of matrix A. It is the size of csrValA and csrColIndA.
descrA host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA host integer array of m elements that contains the start of every row.
csrEndPtrA host integer array of m elements that contains the end of the last row plus one.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix A.
Output
parameter MemorySpace description
issym host 1 if A is symmetric; 0 otherwise.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

High Level Function Reference

This section describes high level API of cuSolverSP, including linear solver, least-square solver and eigenvalue solver. The high-level API is designed for ease-of-use, so it allocates any required memory under the hood automatically. If the host or GPU system memory is not enough, an error is returned.

6.2.1. cusolverSp<t>csrlsvlu()

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 * x = b

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 would be wrong.

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

P * A = L * U

cusolver library provides two reordering schemes, symrcm and symamd, to reduce zero fill-in which dramactically affects the performance of LU factorization. The input parameter reorder can enable symrcm (or symamd) if reorder is 1 (or 2), otherwise, no reordering is performed.

If reorder is nonzero, csrlsvlu does

P * A * Q T = L * U

where Q = symrcm ( A + A T ) .

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

|U(j,j)| < tol for some j

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 performs 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. 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 propose 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. The idea to avoid computing A * A T in [7] is adopted. 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: approximate minimum degree ordering (symamd) is a well-known technique to reduce zero fill-in of QR factorization. However in most cases, symrcm still performs well.

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

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

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
n host host number of rows and columns of matrix A.
nnzA host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnzA ( = csrRowPtrA(n) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of n + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnzA ( = csrRowPtrA(n) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
b device host right hand side vector of size n.
tol host host tolerance to decide if singular or not.
reorder host host no ordering if reorder=0. Otherwise, symrcm is used to reduce zero fill-in.
Output
parameter cusolverSp MemSpace *Host MemSpace description
x device host solution vector of size n, x = inv(A)*b.
singularity host host -1 if A is invertible. Otherwise, first index j such that U(j,j)≈0
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.2.2. cusolverSp<t>csrlsvqr()

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 * x = b

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 would be wrong.

The linear system is solved by sparse QR factorization,

A = Q*R

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

|R(j,j)| < tol for some j

The output parameter singularity is the smallest index of such j. If A is non-singular, singularity is -1. The singularity 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 R(1,1)≈0.

cusolver library provides two reordering schemes, symrcm and symamd, to reduce zero fill-in which dramactically affects the performance of LU factorization. The input parameter reorder can enable symrcm (or symamd) if reorder is 1 (or 2), otherwise, no reordering is performed.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
m host host number of rows and columns of matrix A.
nnz host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of m + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
b device host right hand side vector of size m.
tol host host tolerance to decide if singular or not.
reorder host host no effect.
Output
parameter cusolverSp MemSpace *Host MemSpace description
x device host solution vector of size m, x = inv(A)*b.
singularity host host -1 if A is invertible. Otherwise, first index j such that R(j,j)≈0
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,nnz<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.2.3. cusolverSp<t>csrlsvchol()

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 * x = b

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 (if parameter reorder is zero). 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 . If parameter reorder is nonzero, the user has to extend A to a full matrix, otherwise the solution would be wrong.

The linear system is solved by sparse Cholesky factorization,

A = G * G H

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 k,k <= tol . singularity is the minimum of such k.

singularity is base-0. If A is positive definite and not near singular under tolerance, singularity is -1. If the user wants to know if A is postive definite or not, tol=0 is enough.

cusolver library provides two reordering schemes, symrcm and symamd, to reduce zero fill-in which dramactically affects the performance of LU factorization. The input parameter reorder can enable symrcm (or symamd) if reorder is 1 (or 2), otherwise, no reordering is performed.

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

Remark 2: 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.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
m host host number of rows and columns of matrix A.
nnz host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of m + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
b device host right hand side vector of size m.
tol host host tolerance to decide singularity.
reorder host host no effect.
Output
parameter cusolverSp MemSpace *Host MemSpace description
x device host solution vector of size m, x = inv(A)*b.
singularity host host -1 if A is symmetric postive definite.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,nnz<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.2.4. cusolverSp<t>csrlsqvqr()

The S and D data types are real valued single and double precision, respectively.
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);
The C and Z data types are complex valued single and double precision, respectively.
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

x = argmin || A * z - b ||

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,

A * P T = Q * R

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:

A * P T = A 1 A 2 = Q 1 Q 2 R 1 1 R 1 2 R 2 2

where R 1 1 and A have the same rank, but R 2 2 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 2 2 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

min || A * x - b || = min || R * y - c ||

or in matrix form

R 1 1 R 1 2 R 2 2 y 1 y 2 = c 1 c 2

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

min || y || subject to R 1 1 * y 1 + R 1 2 * y 2 = c 1

Or equivalently another least-square problem

min|| R 1 1 \ R 1 2 I * y 2 - R 1 1 \ c 1 O ||

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.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolver library context.
m host host number of rows of matrix A.
n host host number of columns of matrix A.
nnz host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of m + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
b device host right hand side vector of size m.
tol host host tolerance to decide rank of A.
Output
parameter cusolverSp MemSpace *Host MemSpace description
rankA host host numerical rank of A.
x device host solution vector of size n, x=pinv(A)*b.
p device host a vector of size n, which represents the permuation matrix P satisfying A*P^T=Q*R.
min_norm host host ||A*x-b||, x=pinv(A)*b.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n,nnz<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.2.5. cusolverSp<t>csreigvsi()

The S and D data types are real valued single and double precision, respectively.
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);
The C and Z data types are complex valued single and double precision, respectively.
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 = λ * 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.

Shift-Inverse Method

Given a initial guess of eigenvalue μ0 and initial vector x0 x (0) = x0 of unit length for j = 0 : maxite solve ( A - μ0 * I ) * x (k+1) = x (k) normalize x (k+1) to unit length compute approx. eigenvalue μ = x H * A * x where x = x (k+1) if || A * x (k+1) - μ * x (k+1) || < tolerance, then stop endfor

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.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolver library context.
m host host number of rows and columns of matrix A.
nnz host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of m + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
mu0 host host initial guess of eigenvalue.
x0 device host initial guess of eigenvector, a vecotr of size m.
maxite host host maximum iterations in shift-inverse method.
tol host host tolerance for convergence.
Output
parameter cusolverSp MemSpace *Host MemSpace description
mu device host approximated eigenvalue nearest mu0 under tolerance.
x device host approximated eigenvector of size m.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,nnz<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.2.6. cusolverSp<t>csreigs()

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

number of algebraic eigenvalues in box B = 1 2 * π * - 1 C P (z) P(z) d z

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.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
m host host number of rows and columns of matrix A.
nnz host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) nonzero elements of matrix A.
csrRowPtrA device host integer array of m + 1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnz ( = csrRowPtrA(m) - csrRowPtrA(0) ) column indices of the nonzero elements of matrix A.
left_bottom_corner host host left bottom corner of the box.
right_upper_corner host host right upper corner of the box.
Output
parameter cusolverSp MemSpace *Host MemSpace description
num_eigs host host number of algebraic eigenvalues in a box.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,nnz<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

Low Level Function Reference

This section describes low level API of cuSolverSP, including symrcm and batched QR.

6.3.1. cusolverSpXcsrsymrcm()


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(p,p) = P * A * P T

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.

Input
parameter *Host MemSpace description
handle host handle to the cuSolverSP library context.
n host number of rows and columns of matrix A.
nnzA host number of nonzeros of matrix A. It is the size of csrValA and csrColIndA.
descrA host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix A.
Output
parameter hsolver description
p host permutation vector of size n.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.3.2. cusolverSpXcsrsymmdq()


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

This function implements Symmetric Minimum Degree Algorithm based on Quotient Graph. It returns a permutation vector p such that A(p,p) would have less zero fill-in during Cholesky factorization. The cuSolverSP library implements symmdq based on the following two papers:

Patrick R. Amestoy, Timothy A. Davis, Iain S. Duff, An Approximate Minimum Degree Ordering Algorithm, SIAM J. Matrix Analysis Applic. Vol 17, no 4, pp. 886-905, Dec. 1996.

Alan George, Joseph W. Liu, A Fast Implementation of the Minimum Degree Algorithm Using Quotient Graphs, ACM Transactions on Mathematical Software, Vol 6, No. 3, September 1980, page 337-358.

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

A(p,p) = P * A * P T

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 mdq 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.

Input
parameter *Host MemSpace description
handle host handle to the cuSolverSP library context.
n host number of rows and columns of matrix A.
nnzA host number of nonzeros of matrix A. It is the size of csrValA and csrColIndA.
descrA host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix A.
Output
parameter hsolver description
p host permutation vector of size n.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.3.3. cusolverSpXcsrsymamd()


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

This function implements Symmetric Approximate Minimum Degree Algorithm based on Quotient Graph. It returns a permutation vector p such that A(p,p) would have less zero fill-in during Cholesky factorization. The cuSolverSP library implements symamd based on the following paper:

Patrick R. Amestoy, Timothy A. Davis, Iain S. Duff, An Approximate Minimum Degree Ordering Algorithm, SIAM J. Matrix Analysis Applic. Vol 17, no 4, pp. 886-905, Dec. 1996.

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

A(p,p) = P * A * P T

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 amd 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.

Input
parameter *Host MemSpace description
handle host handle to the cuSolverSP library context.
n host number of rows and columns of matrix A.
nnzA host number of nonzeros of matrix A. It is the size of csrValA and csrColIndA.
descrA host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix A.
Output
parameter hsolver description
p host permutation vector of size n.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.3.4. cusolverSpXcsrperm()


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

B = P * A * Q T

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 = 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

A = 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0

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

P * A * Q T = 3.0 2.0 1.0 9.0 8.0 7.0 6.0 5.0 4.0

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.

Input
parameter cusolverSp MemSpace description
handle host handle to the cuSolver library context.
m host number of rows of matrix A.
n host number of columns of matrix A.
nnzA host number of nonzeros of matrix A. It is the size of csrValA and csrColIndA.
descrA host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA host integer array of m+1 elements that contains the start of every row and end of last row plus one of matrix A.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix A.
p host left permutation vector of size m.
q host right permutation vector of size n.
map host integer array of nnzA indices. If the user wants to get relationship between A and B, map must be set 0:1:(nnzA-1).
pBuffer host buffer allocated by the user, the size is returned by csrperm_bufferSize().
Output
parameter hsolver description
csrRowPtrA host integer array of m+1 elements that contains the start of every row and end of last row plus one of matrix B.
csrColIndA host integer array of nnzAcolumn indices of the nonzero elements of matrix B.
map host integer array of nnzA indices that maps matrix A to matrix B.
pBufferSizeInBytes host number of bytes of the buffer.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.3.5. cusolverSpXcsrqrBatched()

The create and destroy methods start and end the lifetime of a csrqrInfo object.
cusolverStatus_t 
cusolverSpCreateCsrqrInfo(csrqrInfo_t *info);

cusolverStatus_t 
cusolverSpDestroyCsrqrInfo(csrqrInfo_t info);
Analysis is the same for all data types, but each data type has a unique buffer size.
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);
Calculate buffer sizes for complex valued data types.
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);
The S and D data types are real valued single and double precision, respectively.
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);
The C and Z data types are complex valued single and double precision, respectively.
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

x j = argmin || A j * z - b j || , j = 1,2,..., batchSize

or a set of linear systems

A j x j = b j , j = 1,2,..., batchSize

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 , 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.

Input
parameter cusolverSp MemSpace description
handle host handle to the cuSolverSP library context.
m host number of rows of each matrix Aj.
n host number of columns of each matrix Aj.
nnzA host number of nonzeros of each matrix Aj. It is the size csrColIndA.
descrA host the descriptor of each matrix Aj. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device <type> array of nnzA*batchSize nonzero elements of matrices A0, A1, .... All matrices are aggregated one after another.
csrRowPtrA device integer array of m+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device integer array of nnzAcolumn indices of the nonzero elements of each matrix Aj.
b device <type> array of m*batchSize of right-hand-side vectors b0, b1, .... All vectors are aggregated one after another.
batchSize host number of systems to be solved.
info host opaque structure for QR factorization.
pBuffer device buffer allocated by the user, the size is returned by cusolverSpXcsrqrBufferInfoBatched().
Output
parameter cusolverSp MemSpace description
x device <type> array of m*batchSize of solution vectors x0, x1, .... All vectors are aggregated one after another.
internalDataInBytes host number of bytes of the internal data.
workspaceInBytes host number of bytes of the buffer in numerical factorization.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (m,n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_ARCH_MISMATCH the device only supports compute capability 2.0 and above.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.4. cuda 7.5 Preview

This section describes new low level APIs of cuSolverSP in cuda 7.5. The low level APIs include sparse LU, sparse Cholesky and sparse QR. The user has to include header file cusolverSp_LOWLEVEL_PREVIEW.h.

LU, Cholesky and QR have the same flow, including

  • analysis phase to find sparsity pattern of numerical factor.
  • query size of buffer.
  • numerical factorization.
  • report singularity of numerical factorization.
  • numerical solve to complete linear solver or least-square solver.

The user has to follow the above sequence to perform either a linear solver or a least-square solver.

6.4.1. cusolverSpXcsrlu()

The sparse LU factorization is used to factorize matrix A in the following form

P * A * Q T = L * U

A is a n×n sparse matrix that is defined in CSR storage format by the three arrays csrValA, csrRowPtrA and csrColIndA. P is a left permutation matrix mainly on pivoting and Q is a right permutation matrix from postordering of the elimination tree. L is a lower triangular matrix with implicit diagonal one while U is a upper triangular matrix.

If A is symmetric, the user has to extend it to a full matrix and sets the matrix type as CUSPARSE_MATRIX_TYPE_GENERAL.

The low-level API does not reorder the matrix to minimize zero fill-in. The user can use cusolverSpXcsrsymrcm or cusolverSpXcsrsymamd to reorder the matrix to reduce zero fill-in.

cusolverSP LU can be first step of refactorization. Please refer SDK samples/7_CUDALibraries/cuSolverRf.

6.4.1.1. cusolverSpCreateCsrluInfo()

The create and destroy methods start and end the lifetime of a csrluInfo object.
cusolverStatus_t 
cusolverSpCreateCsrluInfo[Host](csrluInfo[Host]_t *info);

cusolverStatus_t 
cusolverSpDestroyCsrluInfo[Host](csrluInfo[Host]_t info);

The function cusolverSpCreateCsrluInfo creates and initializes the opaque structure of LU to default values.

The function cusolverSpDestroyCsrluInfo releases any memory required by the structure.

Remark 1: only CPU path is provided.

Output
parameter cusolverSp MemSpace *Host MemSpace description
info host host opaque structure for LU factorization.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.

6.4.1.2. cusolverSpXcsrluAnalysis()

cusolverStatus_t 
cusolverSpXcsrluAnalysis[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info);

This function analyzes sparsity pattern of matrix L and matrix U of LU factorization. The pivoting is determined at runtime, so only superset of L and U can be found. After analysis, the size of working space to perform LU can be retrieved from cusolverSpXcsrluBufferInfo.

The analysis phase needs working space to estimate sparsity pattern of L and U. If host memory is not sufficient to finish the analysis, CUSOLVER_STATUS_ALLOC_FAILED is returned.

Remark 1: only CPU path is provided.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
n host host number of rows and columns of matrix A.
nnzA host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrRowPtrA device host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnzAcolumn indices of the nonzero elements.
Output
parameter cusolverSp MemSpace *Host MemSpace description
info host host recording scheduling information used in numerical factorization.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.4.1.3. cusolverSpXcsrluBufferInfo()

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverSpScsrluBufferInfo[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const float *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpDcsrluBufferInfo[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const double *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);
The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverSpCcsrluBufferInfo[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

cusolverStatus_t 
cusolverSpZcsrluBufferInfo[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuDoubleComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           size_t *internalDataInBytes,
                           size_t *workspaceInBytes);

There are two memory blocks used in sparse LU. One is internal memory used to store matrix L and matrix U. The other is working space used to perform numerical factorization. The size of the former is specified by returned parameter internalDataInBytes; while the size of the latter is specified by returned parameter workspaceInBytes.

The first call of cusolverSpXcsrluFactor would allocate L and U whose size is bounded by internalDataInBytes. Once internal memory (of size internalDataInBytes bytes) is allocated by cusolverSpXcsrluFactor, the life time is the same as info. Such internal memory is different from working space of size workspaceInBytes bytes, whose life time starts at the beginning of the calling function and ends when the function returns.

Remark 1: only CPU path is provided.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
n host host number of rows and columns of matrix A.
nnzA host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnzA nonzero elements of matrix A.
csrRowPtrA device host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnzAcolumn indices of the nonzero elements.
info host host opaque structure for LU factorization.
Output
parameter cusolverSp MemSpace *Host MemSpace description
internalDataInBytes host host number of bytes of the internal data.
workspaceInBytes host host number of bytes of the buffer in numerical factorization.
info host host recording internal parameters for buffer.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.4.1.4. cusolverSpXcsrluFactor()

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverSpScsrluFactor[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const float *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           float pivot_threshold,
                           void *pBuffer);

cusolverSpDcsrluFactor[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const double *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           double pivot_threshold,
                           void *pBuffer);
The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t 
cusolverSpCcsrluFactor[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           float pivot_threshold,
                           void *pBuffer);

cusolverStatus_t 
cusolverSpZcsrluFactor[Host](cusolverSpHandle_t handle,
                           int n,
                           int nnzA,
                           const cusparseMatDescr_t descrA,
                           const cuDoubleComplex *csrValA,
                           const int *csrRowPtrA,
                           const int *csrColIndA,
                           csrluInfo[Host]_t info,
                           double pivot_threshold,
                           void *pBuffer);

This function performs numerical factorization

P * A * Q T = L * U

The first call to cusolverSpXcsrluFactor would allocate space for L and U. If the memory is insufficient, CUSOLVER_STATUS_ALLOC_FAILED is returned. The numerical factor L and U are kept in structure info and can be used in cusolverSpXcsrluSolve.

The parameter pivot_threshold is for diagonal pivoting. The value is between 0 and 1. If pivot_threshold is 0, then no pivoting is chosen; if pivot_threshold is 1, traditional pivoting is chosen. Assuming that first j-1 columns are done, A is updated, and ξ = max{|A(j:end,j)|} is the condition of traditional pivoting, the formula to choose diagonal A(j,j) as the pivot is

pivot_threshold * ξ <= | A j,j |

Remark 1: only CPU path is provided.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
n host host number of rows and columns of matrix A.
nnzA host host number of nonzeros of matrix A.
descrA host host the descriptor of matrix A. The supported matrix type is CUSPARSE_MATRIX_TYPE_GENERAL. Also, the supported index bases are CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_INDEX_BASE_ONE.
csrValA device host <type> array of nnzA nonzero elements of matrix A.
csrRowPtrA device host integer array of n+1 elements that contains the start of every row and the end of the last row plus one.
csrColIndA device host integer array of nnzAcolumn indices of the nonzero elements.
info host host opaque structure for LU factorization.
pivot_threshold host host a threshold to enable diagonal pivoting.
pBuffer device host buffer allocated by the user, the size is returned by cusolverSpXcsrluBufferInfo().
Output
parameter cusolverSp MemSpace *Host MemSpace description
info host host containing numerical factor L and Q.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_ALLOC_FAILED the resources could not be allocated.
CUSOLVER_STATUS_INVALID_VALUE invalid parameters were passed (n,nnzA<=0), base index is not 0 or 1.
CUSOLVER_STATUS_INTERNAL_ERROR an internal operation failed.
CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED the matrix type is not supported.

6.4.1.5. cusolverSpXcsrluZeroPivot()

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t
cusolverSpScsrluZeroPivot[Host](cusolverSpHandle_t handle,
                           csrluInfo[Host]_t info,
                           float tol,
                           int *position);

cusolverStatus_t
cusolverSpDcsrluZeroPivot[Host](cusolverSpHandle_t handle,
                           csrluInfo[Host]_t info,
                           double tol,
                           int *position);
The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t
cusolverSpCcsrluZeroPivot[Host](cusolverSpHandle_t handle,
                           csrluInfo[Host]_t info,
                           float tol,
                           int *position);
cusolverStatus_t
cusolverSpZcsrluZeroPivot[Host](cusolverSpHandle_t handle,
                           csrluInfo[Host]_t info,
                           double tol,
                           int *position);

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

|U(j,j)| < tol for some j

The output parameter position is the smallest index of such j. If A is non-singular, position 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 position = 1 which means U(1,1)≈0.

The numerical factorization must be done before calling this function, otherwise, CUSOLVER_STATUS_INVALID_VALUE is returned.

Remark 1: only a CPU path is provided.

Remark 2: This routine is not intended to prove that a matrix is singular or non-singular, but to show the need for pivoting. When the pivot threshold is set to 0.0 (no pivoting) this routine may return false positives, ie show a zero pivot when the matrix is not singular. When the pivoting threshold is 1.0, you can trust the output of the zero-pivot routine.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
info host host opaque structure for LU factorization.
tol host host tolerance to determine singularity.
Output
parameter cusolverSp MemSpace *Host MemSpace description
position host host -1 if A is non-singular; otherwise, first column that U(j,j) is zero under given tolerance.
Status Returned
CUSOLVER_STATUS_SUCCESS the operation completed successfully.
CUSOLVER_STATUS_NOT_INITIALIZED the library was not initialized.
CUSOLVER_STATUS_INVALID_VALUE invalid calling sequence.

6.4.1.6. cusolverSpXcsrluSolve()

The S and D data types are real valued single and double precision, respectively.
cusolverStatus_t 
cusolverSpScsrluSolve[Host](cusolverSpHandle_t handle,
                           int n,
                           const float *b,
                           float *x,
                           csrluInfo[Host]_t info,
                           void *pBuffer);

cusolverStatus_t 
cusolverSpDcsrluSolve[Host](cusolverSpHandle_t handle,
                           int n,
                           const double *b,
                           double *x,
                           csrluInfo[Host]_t info,
                           void *pBuffer);
The C and Z data types are complex valued single and double precision, respectively.
cusolverStatus_t
cusolverSpCcsrluSolve[Host](cusolverSpHandle_t handle,
                           int n,
                           const cuComplex *b,
                           cuComplex *x,
                           csrluInfo[Host]_t info,
                           void *pBuffer);

cusolverStatus_t
cusolverSpZcsrluSolve[Host](cusolverSpHandle_t handle,
                           int n,
                           const cuDoubleComplex *b,
                           cuDoubleComplex *x,
                           csrluInfo[Host]_t info,
                           void *pBuffer);

This function solves the linear system A * x = b by forward and backward substitution. The user has to complete numerical factorization before calling this function. If numerical factorization is not done, CUSOLVER_STATUS_INVALID_VALUE is returned.

The numerical factorization must be done before calling this function, otherwise, CUSOLVER_STATUS_INVALID_VALUE is returned.

Remark 1: only CPU path is provided.

Input
parameter cusolverSp MemSpace *Host MemSpace description
handle host host handle to the cuSolverSP library context.
n host host number of rows and columns of matrix A.
b device host <type> array of n of right-hand-side vectors b.
info host host opaque structure for LU factorization.
pBuffer device host buffer allocated by the user, the size is returned by cusolverSpXcsrluBufferInfo().
Output