summaryrefslogtreecommitdiff
path: root/SRC/chetrs_3.f
diff options
context:
space:
mode:
authorJulie <julie@cs.utk.edu>2016-11-15 20:39:35 -0800
committerJulie <julie@cs.utk.edu>2016-11-15 20:39:35 -0800
commitead2c73f1a6dad1342bf32987c0b2f2eaf61f18a (patch)
treeb82e9ad49e12960ad410a418d03d68adc7e2e653 /SRC/chetrs_3.f
parent39698bc46ca55081ebd94c81c5c95771c9f125cd (diff)
downloadlapack-ead2c73f1a6dad1342bf32987c0b2f2eaf61f18a.tar.gz
lapack-ead2c73f1a6dad1342bf32987c0b2f2eaf61f18a.tar.bz2
lapack-ead2c73f1a6dad1342bf32987c0b2f2eaf61f18a.zip
Added (S,D,C,Z) (SY,HE) routines, drivers for new rook code
Close #82 Added routines for new factorization code for symmetric indefinite ( or Hermitian indefinite ) matrices with bounded Bunch-Kaufman ( rook ) pivoting algorithm. New more efficient storage format for factors U ( or L ), block-diagonal matrix D, and pivoting information stored in IPIV: factor L is stored explicitly in lower triangle of A; diagonal of D is stored on the diagonal of A; subdiagonal elements of D are stored in array E; IPIV format is the same as in *_ROOK routines, but differs from SY Bunch-Kaufman routines (e.g. *SYTRF). The factorization output of these new rook _RK routines is not compatible with the existing _ROOK routines and vice versa. This new factorization format is designed in such a way, that there is a possibility in the future to write new Bunch-Kaufman routines that conform to this new factorization format. Then the future Bunch-Kaufman routines could share solver *TRS_3,inversion *TRI_3 and condition estimator *CON_3. To convert between the factorization formats in both ways the following routines are developed: CONVERSION ROUTINES BETWEEN FACTORIZATION FORMATS DOUBLE PRECISION (symmetric indefinite matrices): new file: SRC/dsyconvf.f new file: SRC/dsyconvf_rook.f REAL (symmetric indefinite matrices): new file: SRC/csyconvf.f new file: SRC/csyconvf_rook.f COMPLEX*16 (symmetric indefinite and Hermitian indefinite matrices): new file: SRC/zsyconvf.f new file: SRC/zsyconvf_rook.f COMPLEX (symmetric indefinite and Hermitian indefinite matrices): new file: SRC/ssyconvf.f new file: SRC/ssyconvf_rook.f *SYCONVF routine converts between old Bunch-Kaufman storage format ( denote (L1,D1,IPIV1) ) that is used by *SYTRF and new rook storage format ( denote (L2,D2, IPIV2)) that is used by *SYTRF_RK *SYCONVF_ROOK routine between old rook storage format ( denote (L1,D1,IPIV2) ) that is used by *SYTRF_ROOK and new rook storage format ( denote (L2,D2, IPIV2)) that is used by *SYTRF_RK ROUTINES AND DRIVERS DOUBLE PRECISION (symmetric indefinite matrices): new file: SRC/dsytf2_rk.f BLAS2 unblocked factorization new file: SRC/dlasyf_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/dsytrf_rk.f BLAS3 blocked factorization new file: SRC/dsytrs_3.f BLAS3 solver new file: SRC/dsycon_3.f BLAS3 condition number estimator new file: SRC/dsytri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/dsytri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/dsysv_rk.f BLAS3 solver driver REAL (symmetric indefinite matrices): new file: SRC/ssytf2_rk.f BLAS2 unblocked factorization new file: SRC/slasyf_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/ssytrf_rk.f BLAS3 blocked factorization new file: SRC/ssytrs_3.f BLAS3 solver new file: SRC/ssycon_3.f BLAS3 condition number estimator new file: SRC/ssytri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/ssytri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/ssysv_rk.f BLAS3 solver driver COMPLEX*16 (symmetric indefinite matrices): new file: SRC/zsytf2_rk.f BLAS2 unblocked factorization new file: SRC/zlasyf_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/zsytrf_rk.f BLAS3 blocked factorization new file: SRC/zsytrs_3.f BLAS3 solver new file: SRC/zsycon_3.f BLAS3 condition number estimator new file: SRC/zsytri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/zsytri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/zsysv_rk.f BLAS3 solver driver COMPLEX*16 (Hermitian indefinite matrices): new file: SRC/zhetf2_rk.f BLAS2 unblocked factorization new file: SRC/zlahef_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/zhetrf_rk.f BLAS3 blocked factorization new file: SRC/zhetrs_3.f BLAS3 solver new file: SRC/zhecon_3.f BLAS3 condition number estimator new file: SRC/zhetri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/zhetri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/zhesv_rk.f BLAS3 solver driver COMPLEX (symmetric indefinite matrices): new file: SRC/csytf2_rk.f BLAS2 unblocked factorization new file: SRC/clasyf_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/csytrf_rk.f BLAS3 blocked factorization new file: SRC/csytrs_3.f BLAS3 solver new file: SRC/csycon_3.f BLAS3 condition number estimator new file: SRC/csytri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/csytri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/csysv_rk.f BLAS3 solver driver COMPLEX (Hermitian indefinite matrices): new file: SRC/chetf2_rk.f BLAS2 unblocked factorization new file: SRC/clahef_rk.f BLAS3 auxiliary blocked partial factorization new file: SRC/chetrf_rk.f BLAS3 blocked factorization new file: SRC/chetrs_3.f BLAS3 solver new file: SRC/checon_3.f BLAS3 condition number estimator new file: SRC/chetri_3.f BLAS3 inversion, sets the size of work array and calls *sytri_3x new file: SRC/chetri_3x.f BLAS3 auxiliary inversion, actually computes blocked inversion new file: SRC/chesv_rk.f BLAS3 solver driver MISC modified: SRC/CMakeLists.txt modified: SRC/Makefile TEST CODE modified: TESTING/LIN/CMakeLists.txt modified: TESTING/LIN/Makefile modified: TESTING/LIN/aladhd.f modified: TESTING/LIN/alaerh.f modified: TESTING/LIN/alahd.f DOUBLE PRECISION (symmetric indefinite matrices): modified: TESTING/LIN/dchkaa.f modified: TESTING/LIN/derrsy.f modified: TESTING/LIN/derrsyx.f modified: TESTING/LIN/derrvx.f modified: TESTING/LIN/derrvxx.f modified: TESTING/dtest.in new file: TESTING/LIN/dchksy_rk.f new file: TESTING/LIN/ddrvsy_rk.f new file: TESTING/LIN/dsyt01_3.f REAL (symmetric indefinite matrices): modified: TESTING/LIN/schkaa.f modified: TESTING/LIN/serrsy.f modified: TESTING/LIN/serrsyx.f modified: TESTING/LIN/serrvx.f modified: TESTING/LIN/serrvxx.f modified: TESTING/stest.in new file: TESTING/LIN/schksy_rk.f new file: TESTING/LIN/sdrvsy_rk.f new file: TESTING/LIN/ssyt01_3.f COMPLEX*16 (symmetric indefinite and Hermitian indefinite matrices): modified: TESTING/LIN/zchkaa.f modified: TESTING/LIN/zerrsy.f modified: TESTING/LIN/zerrsyx.f modified: TESTING/LIN/zerrhe.f modified: TESTING/LIN/zerrhex.f modified: TESTING/LIN/zerrvx.f modified: TESTING/LIN/zerrvxx.f modified: TESTING/ztest.in new file: TESTING/LIN/zchksy_rk.f new file: TESTING/LIN/zdrvsy_rk.f new file: TESTING/LIN/zsyt01_3.f new file: TESTING/LIN/zchkhe_rk.f new file: TESTING/LIN/zdrvhe_rk.f new file: TESTING/LIN/zhet01_3.f COMPLEX (symmetric indefinite and Hermitian indefinite matrices): modified: TESTING/LIN/cchkaa.f modified: TESTING/LIN/cerrsy.f modified: TESTING/LIN/cerrsyx.f modified: TESTING/LIN/cerrhe.f modified: TESTING/LIN/cerrhex.f modified: TESTING/LIN/cerrvx.f modified: TESTING/LIN/cerrvxx.f modified: TESTING/ctest.in new file: TESTING/LIN/cchksy_rk.f new file: TESTING/LIN/cdrvsy_rk.f new file: TESTING/LIN/csyt01_3.f new file: TESTING/LIN/cchkhe_rk.f new file: TESTING/LIN/cdrvhe_rk.f new file: TESTING/LIN/chet01_3.f
Diffstat (limited to 'SRC/chetrs_3.f')
-rw-r--r--SRC/chetrs_3.f374
1 files changed, 374 insertions, 0 deletions
diff --git a/SRC/chetrs_3.f b/SRC/chetrs_3.f
new file mode 100644
index 00000000..2799aa24
--- /dev/null
+++ b/SRC/chetrs_3.f
@@ -0,0 +1,374 @@
+*> \brief \b CHETRS_3
+*
+* =========== DOCUMENTATION ===========
+*
+* Online html documentation available at
+* http://www.netlib.org/lapack/explore-html/
+*
+*> \htmlonly
+*> Download CHETRS_3 + dependencies
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/chetrs_3.f">
+*> [TGZ]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/chetrs_3.f">
+*> [ZIP]</a>
+*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/chetrs_3.f">
+*> [TXT]</a>
+*> \endhtmlonly
+*
+* Definition:
+* ===========
+*
+* SUBROUTINE CHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
+* INFO )
+*
+* .. Scalar Arguments ..
+* CHARACTER UPLO
+* INTEGER INFO, LDA, LDB, N, NRHS
+* ..
+* .. Array Arguments ..
+* INTEGER IPIV( * )
+* COMPLEX A( LDA, * ), B( LDB, * ), E( * )
+* ..
+*
+*
+*> \par Purpose:
+* =============
+*>
+*> \verbatim
+*> CHETRS_3 solves a system of linear equations A * X = B with a complex
+*> Hermitian matrix A using the factorization computed
+*> by CHETRF_RK or CHETRF_BK:
+*>
+*> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T),
+*>
+*> where U (or L) is unit upper (or lower) triangular matrix,
+*> U**H (or L**H) is the conjugate of U (or L), P is a permutation
+*> matrix, P**T is the transpose of P, and D is Hermitian and block
+*> diagonal with 1-by-1 and 2-by-2 diagonal blocks.
+*>
+*> This algorithm is using Level 3 BLAS.
+*> \endverbatim
+*
+* Arguments:
+* ==========
+*
+*> \param[in] UPLO
+*> \verbatim
+*> UPLO is CHARACTER*1
+*> Specifies whether the details of the factorization are
+*> stored as an upper or lower triangular matrix:
+*> = 'U': Upper triangular, form is A = P*U*D*(U**H)*(P**T);
+*> = 'L': Lower triangular, form is A = P*L*D*(L**H)*(P**T).
+*> \endverbatim
+*>
+*> \param[in] N
+*> \verbatim
+*> N is INTEGER
+*> The order of the matrix A. N >= 0.
+*> \endverbatim
+*>
+*> \param[in] NRHS
+*> \verbatim
+*> NRHS is INTEGER
+*> The number of right hand sides, i.e., the number of columns
+*> of the matrix B. NRHS >= 0.
+*> \endverbatim
+*>
+*> \param[in] A
+*> \verbatim
+*> A is COMPLEX array, dimension (LDA,N)
+*> Diagonal of the block diagonal matrix D and factors U or L
+*> as computed by CHETRF_RK and CHETRF_BK:
+*> a) ONLY diagonal elements of the Hermitian block diagonal
+*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k);
+*> (superdiagonal (or subdiagonal) elements of D
+*> should be provided on entry in array E), and
+*> b) If UPLO = 'U': factor U in the superdiagonal part of A.
+*> If UPLO = 'L': factor L in the subdiagonal part of A.
+*> \endverbatim
+*>
+*> \param[in] LDA
+*> \verbatim
+*> LDA is INTEGER
+*> The leading dimension of the array A. LDA >= max(1,N).
+*> \endverbatim
+*>
+*> \param[in] E
+*> \verbatim
+*> E is COMPLEX array, dimension (N)
+*> On entry, contains the superdiagonal (or subdiagonal)
+*> elements of the Hermitian block diagonal matrix D
+*> with 1-by-1 or 2-by-2 diagonal blocks, where
+*> If UPLO = 'U': E(i) = D(i-1,i),i=2:N, E(1) not refernced;
+*> If UPLO = 'L': E(i) = D(i+1,i),i=1:N-1, E(N) not referenced.
+*>
+*> NOTE: For 1-by-1 diagonal block D(k), where
+*> 1 <= k <= N, the element E(k) is not referenced in both
+*> UPLO = 'U' or UPLO = 'L' cases.
+*> \endverbatim
+*>
+*> \param[in] IPIV
+*> \verbatim
+*> IPIV is INTEGER array, dimension (N)
+*> Details of the interchanges and the block structure of D
+*> as determined by CHETRF_RK or CHETRF_BK.
+*> \endverbatim
+*>
+*> \param[in,out] B
+*> \verbatim
+*> B is COMPLEX array, dimension (LDB,NRHS)
+*> On entry, the right hand side matrix B.
+*> On exit, the solution matrix X.
+*> \endverbatim
+*>
+*> \param[in] LDB
+*> \verbatim
+*> LDB is INTEGER
+*> The leading dimension of the array B. LDB >= max(1,N).
+*> \endverbatim
+*>
+*> \param[out] INFO
+*> \verbatim
+*> INFO is INTEGER
+*> = 0: successful exit
+*> < 0: if INFO = -i, the i-th argument had an illegal value
+*> \endverbatim
+*
+* Authors:
+* ========
+*
+*> \author Univ. of Tennessee
+*> \author Univ. of California Berkeley
+*> \author Univ. of Colorado Denver
+*> \author NAG Ltd.
+*
+*> \date November 2016
+*
+*> \ingroup complexHEcomputational
+*
+*> \par Contributors:
+* ==================
+*>
+*> \verbatim
+*>
+*> November 2016, Igor Kozachenko,
+*> Computer Science Division,
+*> University of California, Berkeley
+*>
+*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas,
+*> School of Mathematics,
+*> University of Manchester
+*>
+*> \endverbatim
+*
+* =====================================================================
+ SUBROUTINE CHETRS_3( UPLO, N, NRHS, A, LDA, E, IPIV, B, LDB,
+ $ INFO )
+*
+* -- LAPACK computational routine (version 3.7.0) --
+* -- LAPACK is a software package provided by Univ. of Tennessee, --
+* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
+* November 2016
+*
+* .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+* ..
+* .. Array Arguments ..
+ INTEGER IPIV( * )
+ COMPLEX A( LDA, * ), B( LDB, * ), E( * )
+* ..
+*
+* =====================================================================
+*
+* .. Parameters ..
+ COMPLEX ONE
+ PARAMETER ( ONE = ( 1.0E+0,0.0E+0 ) )
+* ..
+* .. Local Scalars ..
+ LOGICAL UPPER
+ INTEGER I, J, K, KP
+ REAL S
+ COMPLEX AK, AKM1, AKM1K, BK, BKM1, DENOM
+* ..
+* .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+* ..
+* .. External Subroutines ..
+ EXTERNAL CSSCAL, CSWAP, CTRSM, XERBLA
+* ..
+* .. Intrinsic Functions ..
+ INTRINSIC ABS, CONJG, MAX, REAL
+* ..
+* .. Executable Statements ..
+*
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'CHETRS_3', -INFO )
+ RETURN
+ END IF
+*
+* Quick return if possible
+*
+ IF( N.EQ.0 .OR. NRHS.EQ.0 )
+ $ RETURN
+*
+ IF( UPPER ) THEN
+*
+* Begin Upper
+*
+* Solve A*X = B, where A = U*D*U**H.
+*
+* P**T * B
+*
+* Interchange rows K and IPIV(K) of matrix B in the same order
+* that the formation order of IPIV(I) vector for Upper case.
+*
+* (We can do the simple loop over IPIV with decrement -1,
+* since the ABS value of IPIV(I) represents the row index
+* of the interchange with row i in both 1x1 and 2x2 pivot cases)
+*
+ DO K = N, 1, -1
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K ) THEN
+ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ END DO
+*
+* Compute (U \P**T * B) -> B [ (U \P**T * B) ]
+*
+ CALL CTRSM( 'L', 'U', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
+*
+* Compute D \ B -> B [ D \ (U \P**T * B) ]
+*
+ I = N
+ DO WHILE ( I.GE.1 )
+ IF( IPIV( I ).GT.0 ) THEN
+ S = REAL( ONE ) / REAL( A( I, I ) )
+ CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
+ ELSE IF ( I.GT.1 ) THEN
+ AKM1K = E( I )
+ AKM1 = A( I-1, I-1 ) / AKM1K
+ AK = A( I, I ) / CONJG( AKM1K )
+ DENOM = AKM1*AK - ONE
+ DO J = 1, NRHS
+ BKM1 = B( I-1, J ) / AKM1K
+ BK = B( I, J ) / CONJG( AKM1K )
+ B( I-1, J ) = ( AK*BKM1-BK ) / DENOM
+ B( I, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ END DO
+ I = I - 1
+ END IF
+ I = I - 1
+ END DO
+*
+* Compute (U**H \ B) -> B [ U**H \ (D \ (U \P**T * B) ) ]
+*
+ CALL CTRSM( 'L', 'U', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB )
+*
+* P * B [ P * (U**H \ (D \ (U \P**T * B) )) ]
+*
+* Interchange rows K and IPIV(K) of matrix B in reverse order
+* from the formation order of IPIV(I) vector for Upper case.
+*
+* (We can do the simple loop over IPIV with increment 1,
+* since the ABS value of IPIV(I) represents the row index
+* of the interchange with row i in both 1x1 and 2x2 pivot cases)
+*
+ DO K = 1, N, 1
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K ) THEN
+ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ END DO
+*
+ ELSE
+*
+* Begin Lower
+*
+* Solve A*X = B, where A = L*D*L**H.
+*
+* P**T * B
+* Interchange rows K and IPIV(K) of matrix B in the same order
+* that the formation order of IPIV(I) vector for Lower case.
+*
+* (We can do the simple loop over IPIV with increment 1,
+* since the ABS value of IPIV(I) represents the row index
+* of the interchange with row i in both 1x1 and 2x2 pivot cases)
+*
+ DO K = 1, N, 1
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K ) THEN
+ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ END DO
+*
+* Compute (L \P**T * B) -> B [ (L \P**T * B) ]
+*
+ CALL CTRSM( 'L', 'L', 'N', 'U', N, NRHS, ONE, A, LDA, B, LDB )
+*
+* Compute D \ B -> B [ D \ (L \P**T * B) ]
+*
+ I = 1
+ DO WHILE ( I.LE.N )
+ IF( IPIV( I ).GT.0 ) THEN
+ S = REAL( ONE ) / REAL( A( I, I ) )
+ CALL CSSCAL( NRHS, S, B( I, 1 ), LDB )
+ ELSE IF( I.LT.N ) THEN
+ AKM1K = E( I )
+ AKM1 = A( I, I ) / CONJG( AKM1K )
+ AK = A( I+1, I+1 ) / AKM1K
+ DENOM = AKM1*AK - ONE
+ DO J = 1, NRHS
+ BKM1 = B( I, J ) / CONJG( AKM1K )
+ BK = B( I+1, J ) / AKM1K
+ B( I, J ) = ( AK*BKM1-BK ) / DENOM
+ B( I+1, J ) = ( AKM1*BK-BKM1 ) / DENOM
+ END DO
+ I = I + 1
+ END IF
+ I = I + 1
+ END DO
+*
+* Compute (L**H \ B) -> B [ L**H \ (D \ (L \P**T * B) ) ]
+*
+ CALL CTRSM('L', 'L', 'C', 'U', N, NRHS, ONE, A, LDA, B, LDB )
+*
+* P * B [ P * (L**H \ (D \ (L \P**T * B) )) ]
+*
+* Interchange rows K and IPIV(K) of matrix B in reverse order
+* from the formation order of IPIV(I) vector for Lower case.
+*
+* (We can do the simple loop over IPIV with decrement -1,
+* since the ABS value of IPIV(I) represents the row index
+* of the interchange with row i in both 1x1 and 2x2 pivot cases)
+*
+ DO K = N, 1, -1
+ KP = ABS( IPIV( K ) )
+ IF( KP.NE.K ) THEN
+ CALL CSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB )
+ END IF
+ END DO
+*
+* END Lower
+*
+ END IF
+*
+ RETURN
+*
+* End of CHETRS_3
+*
+ END