--- rpl/lapack/lapack/zhetrs_aa.f 2017/06/17 11:06:49 1.2 +++ rpl/lapack/lapack/zhetrs_aa.f 2020/05/21 21:46:06 1.5 @@ -38,8 +38,8 @@ *> \verbatim *> *> ZHETRS_AA solves a system of linear equations A*X = B with a complex -*> hermitian matrix A using the factorization A = U*T*U**H or -*> A = L*T*L**T computed by ZHETRF_AA. +*> hermitian matrix A using the factorization A = U**H*T*U or +*> A = L*T*L**H computed by ZHETRF_AA. *> \endverbatim * * Arguments: @@ -50,7 +50,7 @@ *> 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 = U*T*U**H; +*> = 'U': Upper triangular, form is A = U**H*T*U; *> = 'L': Lower triangular, form is A = L*T*L**H. *> \endverbatim *> @@ -67,7 +67,7 @@ *> of the matrix B. NRHS >= 0. *> \endverbatim *> -*> \param[in,out] A +*> \param[in] A *> \verbatim *> A is COMPLEX*16 array, dimension (LDA,N) *> Details of factors computed by ZHETRF_AA. @@ -98,14 +98,16 @@ *> The leading dimension of the array B. LDB >= max(1,N). *> \endverbatim *> -*> \param[in] WORK +*> \param[out] WORK *> \verbatim -*> WORK is DOUBLE array, dimension (MAX(1,LWORK)) +*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK)) *> \endverbatim *> *> \param[in] LWORK *> \verbatim -*> LWORK is INTEGER, LWORK >= MAX(1,3*N-2). +*> LWORK is INTEGER +*> The dimension of the array WORK. LWORK >= max(1,3*N-2). +*> \endverbatim *> *> \param[out] INFO *> \verbatim @@ -122,7 +124,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date December 2016 +*> \date November 2017 * *> \ingroup complex16HEcomputational * @@ -130,10 +132,10 @@ SUBROUTINE ZHETRS_AA( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, $ WORK, LWORK, INFO ) * -* -- LAPACK computational routine (version 3.7.0) -- +* -- LAPACK computational routine (version 3.8.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* December 2016 +* November 2017 * IMPLICIT NONE * @@ -160,7 +162,7 @@ EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL ZGTSV, ZSWAP, ZTRSM, XERBLA + EXTERNAL ZGTSV, ZSWAP, ZTRSM, ZLACGV, ZLACPY, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -199,61 +201,80 @@ * IF( UPPER ) THEN * -* Solve A*X = B, where A = U*T*U**T. +* Solve A*X = B, where A = U**H*T*U. * -* Pivot, P**T * B +* 1) Forward substitution with U**H * - DO K = 1, N - KP = IPIV( K ) - IF( KP.NE.K ) - $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) - END DO + IF( N.GT.1 ) THEN * -* Compute (U \P**T * B) -> B [ (U \P**T * B) ] +* Pivot, P**T * B -> B * - CALL ZTRSM('L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, - $ B( 2, 1 ), LDB) + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO * -* Compute T \ B -> B [ T \ (U \P**T * B) ] +* Compute U**H \ B -> B [ (U**H \P**T * B) ] * - CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1) + CALL ZTRSM( 'L', 'U', 'C', 'U', N-1, NRHS, ONE, A( 1, 2 ), + $ LDA, B( 2, 1 ), LDB ) + END IF +* +* 2) Solve with triangular matrix T +* +* Compute T \ B -> B [ T \ (U**H \P**T * B) ] +* + CALL ZLACPY( 'F', 1, N, A(1, 1), LDA+1, WORK(N), 1 ) IF( N.GT.1 ) THEN CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 2*N ), 1) - CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1) + CALL ZLACPY( 'F', 1, N-1, A( 1, 2 ), LDA+1, WORK( 1 ), 1 ) CALL ZLACGV( N-1, WORK( 1 ), 1 ) END IF - CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, - $ INFO) + CALL ZGTSV( N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, + $ INFO ) +* +* 3) Backward substitution with U +* + IF( N.GT.1 ) THEN * -* Compute (U**T \ B) -> B [ U**T \ (T \ (U \P**T * B) ) ] +* Compute U \ B -> B [ U \ (T \ (U**H \P**T * B) ) ] * - CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), LDA, - $ B(2, 1), LDB) + CALL ZTRSM( 'L', 'U', 'N', 'U', N-1, NRHS, ONE, A( 1, 2 ), + $ LDA, B(2, 1), LDB) * -* Pivot, P * B [ P * (U**T \ (T \ (U \P**T * B) )) ] +* Pivot, P * B [ P * (U**H \ (T \ (U \P**T * B) )) ] * - DO K = N, 1, -1 - KP = IPIV( K ) - IF( KP.NE.K ) - $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) - END DO + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO + END IF * ELSE * -* Solve A*X = B, where A = L*T*L**T. +* Solve A*X = B, where A = L*T*L**H. +* +* 1) Forward substitution with L * -* Pivot, P**T * B + IF( N.GT.1 ) THEN +* +* Pivot, P**T * B -> B * - DO K = 1, N - KP = IPIV( K ) - IF( KP.NE.K ) - $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) - END DO + DO K = 1, N + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO * -* Compute (L \P**T * B) -> B [ (L \P**T * B) ] +* Compute L \ B -> B [ (L \P**T * B) ] * - CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, - $ B(2, 1), LDB) + CALL ZTRSM( 'L', 'L', 'N', 'U', N-1, NRHS, ONE, A( 2, 1 ), + $ LDA, B(2, 1), LDB) + END IF +* +* 2) Solve with triangular matrix T * * Compute T \ B -> B [ T \ (L \P**T * B) ] * @@ -266,18 +287,23 @@ CALL ZGTSV(N, NRHS, WORK(1), WORK(N), WORK(2*N), B, LDB, $ INFO) * -* Compute (L**T \ B) -> B [ L**T \ (T \ (L \P**T * B) ) ] +* 3) Backward substitution with L**H +* + IF( N.GT.1 ) THEN * - CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), LDA, - $ B( 2, 1 ), LDB) +* Compute L**H \ B -> B [ L**H \ (T \ (L \P**T * B) ) ] * -* Pivot, P * B [ P * (L**T \ (T \ (L \P**T * B) )) ] + CALL ZTRSM( 'L', 'L', 'C', 'U', N-1, NRHS, ONE, A( 2, 1 ), + $ LDA, B( 2, 1 ), LDB) * - DO K = N, 1, -1 - KP = IPIV( K ) - IF( KP.NE.K ) - $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) - END DO +* Pivot, P * B [ P * (L**H \ (T \ (L \P**T * B) )) ] +* + DO K = N, 1, -1 + KP = IPIV( K ) + IF( KP.NE.K ) + $ CALL ZSWAP( NRHS, B( K, 1 ), LDB, B( KP, 1 ), LDB ) + END DO + END IF * END IF *