From ead2c73f1a6dad1342bf32987c0b2f2eaf61f18a Mon Sep 17 00:00:00 2001 From: Julie Date: Tue, 15 Nov 2016 20:39:35 -0800 Subject: 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 --- SRC/csycon_3.f | 287 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 287 insertions(+) create mode 100644 SRC/csycon_3.f (limited to 'SRC/csycon_3.f') diff --git a/SRC/csycon_3.f b/SRC/csycon_3.f new file mode 100644 index 00000000..91aae29e --- /dev/null +++ b/SRC/csycon_3.f @@ -0,0 +1,287 @@ +*> \brief \b CSYCON_3 +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CSYCON_3 + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND, +* WORK, IWORK, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* REAL ANORM, RCOND +* .. +* .. Array Arguments .. +* INTEGER IPIV( * ), IWORK( * ) +* COMPLEX A( LDA, * ), E ( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> CSYCON_3 estimates the reciprocal of the condition number (in the +*> 1-norm) of a complex symmetric matrix A using the factorization +*> computed by CSYTRF_RK or CSYTRF_BK: +*> +*> A = P*U*D*(U**T)*(P**T) or A = P*L*D*(L**T)*(P**T), +*> +*> where U (or L) is unit upper (or lower) triangular matrix, +*> U**T (or L**T) is the transpose of U (or L), P is a permutation +*> matrix, P**T is the transpose of P, and D is symmetric and block +*> diagonal with 1-by-1 and 2-by-2 diagonal blocks. +*> +*> An estimate is obtained for norm(inv(A)), and the reciprocal of the +*> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))). +*> This routine uses BLAS3 solver CSYTRS_3. +*> \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**T)*(P**T); +*> = 'L': Lower triangular, form is A = P*L*D*(L**T)*(P**T). +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the matrix A. N >= 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 CSYTRF_RK and CSYTRF_BK: +*> a) ONLY diagonal elements of the symmetric 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 symmetric 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 CSYTRF_RK or CSYTRF_BK. +*> \endverbatim +*> +*> \param[in] ANORM +*> \verbatim +*> ANORM is REAL +*> The 1-norm of the original matrix A. +*> \endverbatim +*> +*> \param[out] RCOND +*> \verbatim +*> RCOND is REAL +*> The reciprocal of the condition number of the matrix A, +*> computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an +*> estimate of the 1-norm of inv(A) computed in this routine. +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (2*N) +*> \endverbatim +*> +*> \param[out] IWORK +*> \verbatim +*> IWORK is INTEGER array, dimension (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 complexSYcomputational +* +*> \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 CSYCON_3( UPLO, N, A, LDA, E, IPIV, ANORM, RCOND, + $ WORK, 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, N + REAL ANORM, RCOND +* .. +* .. Array Arguments .. + INTEGER IPIV( * ) + COMPLEX A( LDA, * ), E( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) + COMPLEX CZERO + PARAMETER ( CZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, KASE + REAL AINVNM +* .. +* .. Local Arrays .. + INTEGER ISAVE( 3 ) +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL CLACN2, CSYTRS_3, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input parameters. +* + 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( LDA.LT.MAX( 1, N ) ) THEN + INFO = -4 + ELSE IF( ANORM.LT.ZERO ) THEN + INFO = -7 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CSYCON_3', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + RCOND = ZERO + IF( N.EQ.0 ) THEN + RCOND = ONE + RETURN + ELSE IF( ANORM.LE.ZERO ) THEN + RETURN + END IF +* +* Check that the diagonal matrix D is nonsingular. +* + IF( UPPER ) THEN +* +* Upper triangular storage: examine D from bottom to top +* + DO I = N, 1, -1 + IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO ) + $ RETURN + END DO + ELSE +* +* Lower triangular storage: examine D from top to bottom. +* + DO I = 1, N + IF( IPIV( I ).GT.0 .AND. A( I, I ).EQ.CZERO ) + $ RETURN + END DO + END IF +* +* Estimate the 1-norm of the inverse. +* + KASE = 0 + 30 CONTINUE + CALL CLACN2( N, WORK( N+1 ), WORK, AINVNM, KASE, ISAVE ) + IF( KASE.NE.0 ) THEN +* +* Multiply by inv(L*D*L**T) or inv(U*D*U**T). +* + CALL CSYTRS_3( UPLO, N, 1, A, LDA, E, IPIV, WORK, N, INFO ) + GO TO 30 + END IF +* +* Compute the estimate of the reciprocal condition number. +* + IF( AINVNM.NE.ZERO ) + $ RCOND = ( ONE / AINVNM ) / ANORM +* + RETURN +* +* End of CSYCON_3 +* + END -- cgit v1.2.3