sy/hesv: Solves Ax = b using symmetric/Hermitian indefinite factorization (driver)

Functions
magma_int_t	magma_chesv (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, magmaFloatComplex *A, magma_int_t lda, magma_int_t *ipiv, magmaFloatComplex *B, magma_int_t ldb, magma_int_t *info)
	CHESV computes the solution to a complex system of linear equations A * X = B, where A is an n-by-n Hermitian matrix and X and B are n-by-nrhs matrices. More...
magma_int_t	magma_dssysv_gpu (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, magmaDouble_ptr dA, magma_int_t ldda, magmaDouble_ptr dB, magma_int_t lddb, magmaDouble_ptr dX, magma_int_t lddx, magmaDouble_ptr dworkd, magmaFloat_ptr dworks, magma_int_t *iter, magma_int_t *info)
	DSHESV computes the solution to a real system of linear equations A * X = B, where A is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices. More...
magma_int_t	magma_dsysv (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, double *A, magma_int_t lda, magma_int_t *ipiv, double *B, magma_int_t ldb, magma_int_t *info)
	DSYSV computes the solution to a real system of linear equations A * X = B, where A is an n-by-n symmetric matrix and X and B are n-by-nrhs matrices. More...
magma_int_t	magma_ssysv (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, float *A, magma_int_t lda, magma_int_t *ipiv, float *B, magma_int_t ldb, magma_int_t *info)
	SSYSV computes the solution to a real system of linear equations A * X = B, where A is an n-by-n symmetric matrix and X and B are n-by-nrhs matrices. More...
magma_int_t	magma_zchesv_gpu (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, magmaDoubleComplex_ptr dA, magma_int_t ldda, magmaDoubleComplex_ptr dB, magma_int_t lddb, magmaDoubleComplex_ptr dX, magma_int_t lddx, magmaDoubleComplex_ptr dworkd, magmaFloatComplex_ptr dworks, magma_int_t *iter, magma_int_t *info)
	ZCHESV computes the solution to a complex system of linear equations A * X = B, where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices. More...
magma_int_t	magma_zhesv (magma_uplo_t uplo, magma_int_t n, magma_int_t nrhs, magmaDoubleComplex *A, magma_int_t lda, magma_int_t *ipiv, magmaDoubleComplex *B, magma_int_t ldb, magma_int_t *info)
	ZHESV computes the solution to a complex system of linear equations A * X = B, where A is an n-by-n Hermitian matrix and X and B are n-by-nrhs matrices. More...

Detailed Description

Function Documentation

magma_int_t magma_chesv	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		magmaFloatComplex *	A,
		magma_int_t	lda,
		magma_int_t *	ipiv,
		magmaFloatComplex *	B,
		magma_int_t	ldb,
		magma_int_t *	info
	)

CHESV computes the solution to a complex system of linear equations A * X = B, where A is an n-by-n Hermitian matrix and X and B are n-by-nrhs matrices.

The diagonal pivoting method is used to factor A as A = U * D * U**H, if uplo = MagmaUpper, or A = L * D * L**H, if uplo = MagmaLower, where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. n >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. nrhs >= 0.
[in,out]	A	COMPLEX array, dimension (lda,n) On entry, the Hermitian matrix A. If uplo = MagmaUpper, 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 = MagmaLower, 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, if info = 0, the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**H or A = L*D*L**H as computed by CHETRF.

Parameters

[in]	lda	INTEGER The leading dimension of the array A. lda >= max(1,n).
[out]	ipiv	INTEGER array, dimension (n) Details of the interchanges and the block structure of D, as determined by CHETRF. If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block. If uplo = MagmaUpper and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If uplo = MagmaLower and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[in,out]	B	(input/output) COMPLEX array, dimension (ldb,nrhs) On entry, the n-by-nrhs right hand side matrix B. On exit, if info = 0, the n-by-nrhs solution matrix X.
[in]	ldb	INTEGER The leading dimension of the array B. ldb >= max(1,n).
[out]	info	INTEGER = 0: successful exit < 0: if info = -i, the i-th argument had an illegal value > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.

magma_int_t magma_dssysv_gpu	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		magmaDouble_ptr	dA,
		magma_int_t	ldda,
		magmaDouble_ptr	dB,
		magma_int_t	lddb,
		magmaDouble_ptr	dX,
		magma_int_t	lddx,
		magmaDouble_ptr	dworkd,
		magmaFloat_ptr	dworks,
		magma_int_t *	iter,
		magma_int_t *	info
	)

DSHESV computes the solution to a real system of linear equations A * X = B, where A is an N-by-N symmetric matrix and X and B are N-by-NRHS matrices.

DSHESV first attempts to factorize the matrix in real SINGLE PRECISION (without pivoting) and use this factorization within iterative refinements to produce a solution with real DOUBLE PRECISION norm-wise backward error quality (see below). If the approach fails the method switches to a real DOUBLE PRECISION factorization and solve.

The iterative refinement is not going to be a winning strategy if the ratio real SINGLE PRECISION performance over real DOUBLE PRECISION performance is too small or if there are many right-hand sides. A reasonable strategy should take the number of right-hand sides and the size of the matrix into account. This might be done with a call to ILAENV in the future. For now, we always try iterative refinement.

The iterative refinement process is stopped if ITER > ITERMAX or for all the RHS we have: RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX where o ITER is the number of the current iteration in the iterative refinement process o RNRM is the infinity-norm of the residual o XNRM is the infinity-norm of the solution o ANRM is the infinity-operator-norm of the matrix A o EPS is the machine epsilon returned by DLAMCH('Epsilon') The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
[in,out]	dA	DOUBLE PRECISION array on the GPU, dimension (LDDA,N) On entry, the symmetric matrix A. If UPLO = MagmaUpper, 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 = MagmaLower, 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, if iterative refinement has been successfully used (INFO.EQ.0 and ITER.GE.0, see description below), then A is unchanged, if double factorization has been used (INFO.EQ.0 and ITER.LT.0, see description below), then the array dA contains the factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T.
[in]	ldda	INTEGER The leading dimension of the array dA. LDDA >= max(1,N).
[in]	dB	DOUBLE PRECISION array on the GPU, dimension (LDDB,NRHS) The N-by-NRHS right hand side matrix B.
[in]	lddb	INTEGER The leading dimension of the array dB. LDDB >= max(1,N).
[out]	dX	DOUBLE PRECISION array on the GPU, dimension (LDDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.
[in]	lddx	INTEGER The leading dimension of the array dX. LDDX >= max(1,N).
	dworkd	(workspace) DOUBLE PRECISION array on the GPU, dimension (N*NRHS) This array is used to hold the residual vectors.
	dworks	(workspace) SINGLE PRECISION array on the GPU, dimension (N*(N+NRHS)) This array is used to store the real single precision matrix and the right-hand sides or solutions in single precision.
[out]	iter	INTEGER < 0: iterative refinement has failed, double precision factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of SPOTRF -31: stop the iterative refinement after the 30th iteration > 0: iterative refinement has been successfully used. Returns the number of iterations
[out]	info	INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i of (DOUBLE PRECISION) A is not positive definite, so the factorization could not be completed, and the solution has not been computed.

magma_int_t magma_dsysv	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		double *	A,
		magma_int_t	lda,
		magma_int_t *	ipiv,
		double *	B,
		magma_int_t	ldb,
		magma_int_t *	info
	)

DSYSV computes the solution to a real system of linear equations A * X = B, where A is an n-by-n symmetric matrix and X and B are n-by-nrhs matrices.

The diagonal pivoting method is used to factor A as A = U * D * U**H, if uplo = MagmaUpper, or A = L * D * L**H, if uplo = MagmaLower, where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. n >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. nrhs >= 0.
[in,out]	A	DOUBLE PRECISION array, dimension (lda,n) On entry, the symmetric matrix A. If uplo = MagmaUpper, 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 = MagmaLower, 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, if info = 0, the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**H or A = L*D*L**H as computed by DSYTRF.

Parameters

[in]	lda	INTEGER The leading dimension of the array A. lda >= max(1,n).
[out]	ipiv	INTEGER array, dimension (n) Details of the interchanges and the block structure of D, as determined by DSYTRF. If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block. If uplo = MagmaUpper and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If uplo = MagmaLower and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[in,out]	B	(input/output) DOUBLE PRECISION array, dimension (ldb,nrhs) On entry, the n-by-nrhs right hand side matrix B. On exit, if info = 0, the n-by-nrhs solution matrix X.
[in]	ldb	INTEGER The leading dimension of the array B. ldb >= max(1,n).
[out]	info	INTEGER = 0: successful exit < 0: if info = -i, the i-th argument had an illegal value > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.

magma_int_t magma_ssysv	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		float *	A,
		magma_int_t	lda,
		magma_int_t *	ipiv,
		float *	B,
		magma_int_t	ldb,
		magma_int_t *	info
	)

SSYSV computes the solution to a real system of linear equations A * X = B, where A is an n-by-n symmetric matrix and X and B are n-by-nrhs matrices.

The diagonal pivoting method is used to factor A as A = U * D * U**H, if uplo = MagmaUpper, or A = L * D * L**H, if uplo = MagmaLower, where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is symmetric and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. n >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. nrhs >= 0.
[in,out]	A	REAL array, dimension (lda,n) On entry, the symmetric matrix A. If uplo = MagmaUpper, 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 = MagmaLower, 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, if info = 0, the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**H or A = L*D*L**H as computed by SSYTRF.

Parameters

[in]	lda	INTEGER The leading dimension of the array A. lda >= max(1,n).
[out]	ipiv	INTEGER array, dimension (n) Details of the interchanges and the block structure of D, as determined by SSYTRF. If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block. If uplo = MagmaUpper and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If uplo = MagmaLower and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[in,out]	B	(input/output) REAL array, dimension (ldb,nrhs) On entry, the n-by-nrhs right hand side matrix B. On exit, if info = 0, the n-by-nrhs solution matrix X.
[in]	ldb	INTEGER The leading dimension of the array B. ldb >= max(1,n).
[out]	info	INTEGER = 0: successful exit < 0: if info = -i, the i-th argument had an illegal value > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.

magma_int_t magma_zchesv_gpu	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		magmaDoubleComplex_ptr	dA,
		magma_int_t	ldda,
		magmaDoubleComplex_ptr	dB,
		magma_int_t	lddb,
		magmaDoubleComplex_ptr	dX,
		magma_int_t	lddx,
		magmaDoubleComplex_ptr	dworkd,
		magmaFloatComplex_ptr	dworks,
		magma_int_t *	iter,
		magma_int_t *	info
	)

ZCHESV computes the solution to a complex system of linear equations A * X = B, where A is an N-by-N Hermitian matrix and X and B are N-by-NRHS matrices.

ZCHESV first attempts to factorize the matrix in complex SINGLE PRECISION (without pivoting) and use this factorization within iterative refinements to produce a solution with complex DOUBLE PRECISION norm-wise backward error quality (see below). If the approach fails the method switches to a complex DOUBLE PRECISION factorization and solve.

The iterative refinement is not going to be a winning strategy if the ratio complex SINGLE PRECISION performance over complex DOUBLE PRECISION performance is too small or if there are many right-hand sides. A reasonable strategy should take the number of right-hand sides and the size of the matrix into account. This might be done with a call to ILAENV in the future. For now, we always try iterative refinement.

The iterative refinement process is stopped if ITER > ITERMAX or for all the RHS we have: RNRM < SQRT(N)*XNRM*ANRM*EPS*BWDMAX where o ITER is the number of the current iteration in the iterative refinement process o RNRM is the infinity-norm of the residual o XNRM is the infinity-norm of the solution o ANRM is the infinity-operator-norm of the matrix A o EPS is the machine epsilon returned by DLAMCH('Epsilon') The value ITERMAX and BWDMAX are fixed to 30 and 1.0D+00 respectively.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. N >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. NRHS >= 0.
[in,out]	dA	COMPLEX_16 array on the GPU, dimension (LDDA,N) On entry, the Hermitian matrix A. If UPLO = MagmaUpper, 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 = MagmaLower, 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, if iterative refinement has been successfully used (INFO.EQ.0 and ITER.GE.0, see description below), then A is unchanged, if double factorization has been used (INFO.EQ.0 and ITER.LT.0, see description below), then the array dA contains the factor U or L from the Cholesky factorization A = U**T*U or A = L*L**T.
[in]	ldda	INTEGER The leading dimension of the array dA. LDDA >= max(1,N).
[in]	dB	COMPLEX_16 array on the GPU, dimension (LDDB,NRHS) The N-by-NRHS right hand side matrix B.
[in]	lddb	INTEGER The leading dimension of the array dB. LDDB >= max(1,N).
[out]	dX	COMPLEX_16 array on the GPU, dimension (LDDX,NRHS) If INFO = 0, the N-by-NRHS solution matrix X.
[in]	lddx	INTEGER The leading dimension of the array dX. LDDX >= max(1,N).
	dworkd	(workspace) COMPLEX_16 array on the GPU, dimension (N*NRHS) This array is used to hold the residual vectors.
	dworks	(workspace) COMPLEX array on the GPU, dimension (N*(N+NRHS)) This array is used to store the complex single precision matrix and the right-hand sides or solutions in single precision.
[out]	iter	INTEGER < 0: iterative refinement has failed, double precision factorization has been performed -1 : the routine fell back to full precision for implementation- or machine-specific reasons -2 : narrowing the precision induced an overflow, the routine fell back to full precision -3 : failure of SPOTRF -31: stop the iterative refinement after the 30th iteration > 0: iterative refinement has been successfully used. Returns the number of iterations
[out]	info	INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value > 0: if INFO = i, the leading minor of order i of (DOUBLE PRECISION) A is not positive definite, so the factorization could not be completed, and the solution has not been computed.

magma_int_t magma_zhesv	(	magma_uplo_t	uplo,
		magma_int_t	n,
		magma_int_t	nrhs,
		magmaDoubleComplex *	A,
		magma_int_t	lda,
		magma_int_t *	ipiv,
		magmaDoubleComplex *	B,
		magma_int_t	ldb,
		magma_int_t *	info
	)

ZHESV computes the solution to a complex system of linear equations A * X = B, where A is an n-by-n Hermitian matrix and X and B are n-by-nrhs matrices.

The diagonal pivoting method is used to factor A as A = U * D * U**H, if uplo = MagmaUpper, or A = L * D * L**H, if uplo = MagmaLower, where U (or L) is a product of permutation and unit upper (lower) triangular matrices, and D is Hermitian and block diagonal with 1-by-1 and 2-by-2 diagonal blocks. The factored form of A is then used to solve the system of equations A * X = B.

Parameters

[in]	uplo	magma_uplo_t = MagmaUpper: Upper triangle of A is stored; = MagmaLower: Lower triangle of A is stored.
[in]	n	INTEGER The number of linear equations, i.e., the order of the matrix A. n >= 0.
[in]	nrhs	INTEGER The number of right hand sides, i.e., the number of columns of the matrix B. nrhs >= 0.
[in,out]	A	COMPLEX*16 array, dimension (lda,n) On entry, the Hermitian matrix A. If uplo = MagmaUpper, 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 = MagmaLower, 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, if info = 0, the block diagonal matrix D and the multipliers used to obtain the factor U or L from the factorization A = U*D*U**H or A = L*D*L**H as computed by ZHETRF.

Parameters

[in]	lda	INTEGER The leading dimension of the array A. lda >= max(1,n).
[out]	ipiv	INTEGER array, dimension (n) Details of the interchanges and the block structure of D, as determined by ZHETRF. If ipiv(k) > 0, then rows and columns k and ipiv(k) were interchanged, and D(k,k) is a 1-by-1 diagonal block. If uplo = MagmaUpper and ipiv(k) = ipiv(k-1) < 0, then rows and columns k-1 and -ipiv(k) were interchanged and D(k-1:k,k-1:k) is a 2-by-2 diagonal block. If uplo = MagmaLower and ipiv(k) = ipiv(k+1) < 0, then rows and columns k+1 and -ipiv(k) were interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
[in,out]	B	(input/output) COMPLEX*16 array, dimension (ldb,nrhs) On entry, the n-by-nrhs right hand side matrix B. On exit, if info = 0, the n-by-nrhs solution matrix X.
[in]	ldb	INTEGER The leading dimension of the array B. ldb >= max(1,n).
[out]	info	INTEGER = 0: successful exit < 0: if info = -i, the i-th argument had an illegal value > 0: if info = i, D(i,i) is exactly zero. The factorization has been completed, but the block diagonal matrix D is exactly singular, so the solution could not be computed.