diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index 1a1813c1d..2d2edb73c 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -113,7 +113,7 @@ set(SLASRC slarz.f slarzb.f slarzt.f slasy2.f slasyf.f slasyf_rook.f slasyf_rk.f slasyf_aa.f slatbs.f slatdf.f slatps.f slatrd.f slatrs.f slatrs3.f slatrz.f - slauu2.f slauum.f sopgtr.f sopmtr.f sorg2l.f sorg2r.f + slauu2.f slauum.f slauum_blocked.f slauum_recursive.f sopgtr.f sopmtr.f sorg2l.f sorg2r.f sorgbr.f sorghr.f sorgl2.f sorglq.f sorgql.f sorgqr.f sorgr2.f sorgrq.f sorgtr.f sorgtsqr.f sorgtsqr_row.f sorm2l.f sorm2r.f sorm22.f sormbr.f sormhr.f sorml2.f sormlq.f sormql.f sormqr.f sormr2.f @@ -225,7 +225,7 @@ set(CLASRC clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f90 claswp.f clasyf.f clasyf_rook.f clasyf_rk.f clasyf_aa.f clatbs.f clatdf.f clatps.f clatrd.f clatrs.f clatrs3.f clatrz.f - clauu2.f clauum.f cpbcon.f cpbequ.f cpbrfs.f cpbstf.f cpbsv.f + clauu2.f clauum.f clauum_blocked.f clauum_recursive.f cpbcon.f cpbequ.f cpbrfs.f cpbstf.f cpbsv.f cpbsvx.f cpbtf2.f cpbtrf.f cpbtrs.f cpocon.f cpoequ.f cporfs.f cposv.f cposvx.f cpotf2.f cpotrf2.f cpotri.f cpstrf.f cpstf2.f cppcon.f cppequ.f cpprfs.f cppsv.f cppsvx.f cpptrf.f cpptri.f cpptrs.f @@ -314,7 +314,7 @@ set(DLASRC dlarz.f dlarzb.f dlarzt.f dlaswp.f dlasy2.f dlasyf.f dlasyf_rook.f dlasyf_rk.f dlasyf_aa.f dlatbs.f dlatdf.f dlatps.f dlatrd.f dlatrs.f dlatrs3.f dlatrz.f dlauu2.f - dlauum.f dopgtr.f dopmtr.f dorg2l.f dorg2r.f + dlauum.f dlauum_blocked.f dlauum_recursive.f dopgtr.f dopmtr.f dorg2l.f dorg2r.f dorgbr.f dorghr.f dorgl2.f dorglq.f dorgql.f dorgqr.f dorgr2.f dorgrq.f dorgtr.f dorgtsqr.f dorgtsqr_row.f dorm2l.f dorm2r.f dorm22.f dormbr.f dormhr.f dorml2.f dormlq.f dormql.f dormqr.f dormr2.f @@ -426,7 +426,7 @@ set(ZLASRC zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.f zlassq.f90 zlaswp.f zlasyf.f zlasyf_rook.f zlasyf_rk.f zlasyf_aa.f zlatbs.f zlatdf.f zlatps.f zlatrd.f zlatrs.f zlatrs3.f zlatrz.f zlauu2.f - zlauum.f zpbcon.f zpbequ.f zpbrfs.f zpbstf.f zpbsv.f + zlauum.f zlauum_blocked.f zlauum_recursive.f zpbcon.f zpbequ.f zpbrfs.f zpbstf.f zpbsv.f zpbsvx.f zpbtf2.f zpbtrf.f zpbtrs.f zpocon.f zpoequ.f zporfs.f zposv.f zposvx.f zpotf2.f zpotrf.f zpotrf2.f zpotri.f zpotrs.f zpstrf.f zpstf2.f zppcon.f zppequ.f zpprfs.f zppsv.f zppsvx.f zpptrf.f zpptri.f zpptrs.f diff --git a/SRC/Makefile b/SRC/Makefile index 39d4a0c07..f3e0ca01e 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -142,7 +142,7 @@ SLASRC = \ slarz.o slarzb.o slarzt.o slaswp.o slasy2.o slasyf.o slasyf_rook.o \ slasyf_rk.o \ slatbs.o slatdf.o slatps.o slatrd.o slatrs.o slatrs3.o slatrz.o \ - slauu2.o slauum.o sopgtr.o sopmtr.o sorg2l.o sorg2r.o \ + slauu2.o slauum.o slauum_blocked.o slauum_recursive.o sopgtr.o sopmtr.o sorg2l.o sorg2r.o \ sorgbr.o sorghr.o sorgl2.o sorglq.o sorgql.o sorgqr.o sorgr2.o \ sorgrq.o sorgtr.o sorgtsqr.o sorgtsqr_row.o sorm2l.o sorm2r.o sorm22.o \ sormbr.o sormhr.o sorml2.o sormlq.o sormql.o sormqr.o sormr2.o \ @@ -254,7 +254,7 @@ CLASRC = \ clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \ claswp.o clasyf.o clasyf_rook.o clasyf_rk.o clasyf_aa.o \ clatbs.o clatdf.o clatps.o clatrd.o clatrs.o clatrs3.o clatrz.o \ - clauu2.o clauum.o cpbcon.o cpbequ.o cpbrfs.o cpbstf.o cpbsv.o \ + clauu2.o clauum.o clauum_blocked.o clauum_recursive.o cpbcon.o cpbequ.o cpbrfs.o cpbstf.o cpbsv.o \ cpbsvx.o cpbtf2.o cpbtrf.o cpbtrs.o cpocon.o cpoequ.o cporfs.o \ cposv.o cposvx.o cpotf2.o cpotri.o cpstrf.o cpstf2.o \ cppcon.o cppequ.o cpprfs.o cppsv.o cppsvx.o cpptrf.o cpptri.o cpptrs.o \ @@ -344,7 +344,7 @@ DLASRC = \ dlarz.o dlarzb.o dlarzt.o dlaswp.o dlasy2.o \ dlasyf.o dlasyf_rook.o dlasyf_rk.o \ dlatbs.o dlatdf.o dlatps.o dlatrd.o dlatrs.o dlatrs3.o dlatrz.o dlauu2.o \ - dlauum.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \ + dlauum.o dlauum_blocked.o dlauum_recursive.o dlauum_recursive.o dopgtr.o dopmtr.o dorg2l.o dorg2r.o \ dorgbr.o dorghr.o dorgl2.o dorglq.o dorgql.o dorgqr.o dorgr2.o \ dorgrq.o dorgtr.o dorgtsqr.o dorgtsqr_row.o dorm2l.o dorm2r.o dorm22.o \ dormbr.o dormhr.o dorml2.o dormlq.o dormql.o dormqr.o dormr2.o \ @@ -459,7 +459,7 @@ ZLASRC = \ zlarz.o zlarzb.o zlarzt.o zlascl.o zlaset.o zlasr.o \ zlassq.o zlaswp.o zlasyf.o zlasyf_rook.o zlasyf_rk.o zlasyf_aa.o \ zlatbs.o zlatdf.o zlatps.o zlatrd.o zlatrs.o zlatrs3.o zlatrz.o zlauu2.o \ - zlauum.o zpbcon.o zpbequ.o zpbrfs.o zpbstf.o zpbsv.o \ + zlauum.o zlauum_blocked.o zlauum_recursive.o zpbcon.o zpbequ.o zpbrfs.o zpbstf.o zpbsv.o \ zpbsvx.o zpbtf2.o zpbtrf.o zpbtrs.o zpocon.o zpoequ.o zporfs.o \ zposv.o zposvx.o zpotf2.o zpotrf.o zpotri.o zpotrs.o zpstrf.o zpstf2.o \ zppcon.o zppequ.o zpprfs.o zppsv.o zppsvx.o zpptrf.o zpptri.o zpptrs.o \ diff --git a/SRC/clauum.f b/SRC/clauum.f index 49f4410c5..336a6f458 100644 --- a/SRC/clauum.f +++ b/SRC/clauum.f @@ -157,62 +157,13 @@ SUBROUTINE CLAUUM( UPLO, N, A, LDA, INFO ) IF( N.EQ.0 ) $ RETURN * -* Determine the block size for this environment. -* - NB = ILAENV( 1, 'CLAUUM', UPLO, N, -1, -1, -1 ) -* - IF( NB.LE.1 .OR. NB.GE.N ) THEN -* -* Use unblocked code -* - CALL CLAUU2( UPLO, N, A, LDA, INFO ) - ELSE -* -* Use blocked code -* - IF( UPPER ) THEN -* -* Compute the product U * U**H. -* - DO 10 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', - $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, - $ A( 1, I ), LDA ) - CALL CLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL CGEMM( 'No transpose', 'Conjugate transpose', - $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), - $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), - $ LDA ) - CALL CHERK( 'Upper', 'No transpose', IB, N-I-IB+1, - $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), - $ LDA ) - END IF - 10 CONTINUE - ELSE -* -* Compute the product L**H * L. -* - DO 20 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL CTRMM( 'Left', 'Lower', 'Conjugate transpose', - $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, - $ A( I, 1 ), LDA ) - CALL CLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL CGEMM( 'Conjugate transpose', 'No transpose', - $ IB, - $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, - $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) - CALL CHERK( 'Lower', 'Conjugate transpose', IB, - $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, - $ A( I, I ), LDA ) - END IF - 20 CONTINUE - END IF - END IF +* Here we dispatch to whatever is more efficient in a particular environment +* We are defaulting to recursive, but if you want to use the blocked variant +* Comment out the line starting with `CALL CLAUUM_RECURSIVE...` +* and uncomment the line starting with `CALL CLAUUM_BLOCKED...` * + CALL CLAUUM_RECURSIVE(UPLO, N, A, LDA, INFO) +* CALL CLAUUM_BLOCKED(UPLO, N, A, LDA, INFO) RETURN * * End of CLAUUM diff --git a/SRC/clauum_blocked.f b/SRC/clauum_blocked.f new file mode 100644 index 000000000..80f3cb85d --- /dev/null +++ b/SRC/clauum_blocked.f @@ -0,0 +1,220 @@ +*> \brief \b CLAUUM_BLOCKED computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download CLAUUM_BLOCKED + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE CLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLAUUM_BLOCKED computes the product U * U**H or L**H * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**H; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**H * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE CLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) + COMPLEX CONE + PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL CGEMM, CHERK, CLAUU2, CTRMM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CLAUUM_BLOCKED', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'CLAUUM_BLOCKED', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL CLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U**H. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, + $ A( 1, I ), LDA ) + CALL CLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), + $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), + $ LDA ) + CALL CHERK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L**H * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL CTRMM( 'Left', 'Lower', 'Conjugate transpose', + $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, + $ A( I, 1 ), LDA ) + CALL CLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL CGEMM( 'Conjugate transpose', 'No transpose', + $ IB, + $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) + CALL CHERK( 'Lower', 'Conjugate transpose', IB, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, + $ A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of CLAUUM_BLOCKED +* + END diff --git a/SRC/clauum_recursive.f b/SRC/clauum_recursive.f new file mode 100644 index 000000000..2e819e910 --- /dev/null +++ b/SRC/clauum_recursive.f @@ -0,0 +1,281 @@ +*> \brief \b CLAUUM_RECURSIVE computes the product UUH or LHL, where U and L are upper or lower triangular matrices (recursive algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download CLAUUM_RECURSIVE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE CLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLAUUM_RECURSIVE computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE CLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER K, NX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL CHERK, CTRMM, CLAUU2 +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAUUM_RECURSIVE', -INFO ) + RETURN + END IF +* +* Early termination criteria +* + IF( N.EQ.0 ) THEN + RETURN + END IF +* +* Base Case +* + IF( N.EQ.1 ) THEN + A(1,1) = A(1,1) * A(1,1) + RETURN + END IF +* +* Determine crossover point for when to bail to level2 +* + NX = ILAENV(3, "DLAUUM_RECURSIVE", UPLO, N, -1, -1, -1) + IF( K.LT.NX ) THEN + CALL CLAUU2(UPLO, N, A, LDA, INFO) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + K = N/2 + IF( UPPER ) THEN +* +* We are computing A = ut(U*U**H). +* +* Break apart U as follows +* |-----------------| +* U = | U_{11} U_{12} | +* | 0 U_{22} | +* |-----------------| +* +* Where +* U_{11}\in\R^{k\times k} U_{12}\in\R^{ k\times n-k} +* U_{22}\in\R^{n-k\times n-k} +* +* and U_{11},U_{22} are upper triangular and U_{12} is rectangular +* +* This gives us our operations as +* |--------------------| |------------------------| +* ut(U * U**H) = | U_{11} U_{12} | | U_{11}**H 0 | +* | 0 U_{22} | | U_{12}**H U_{22}**H | +* |--------------------| |------------------------| +* +* Thus we get +* +* U_{11} = U_{11}U_{11}**H + U_{12}U_{12}**H +* U_{12} = U_{12}U_{22}**H +* +* U_{22} = U_{22}U_{22}**H +* +* We break these operations apart as follows +* +* U_{11} = U_{11}U_{11}**H (This subroutine) +* U_{11} = U_{12}U_{12}**H + U_{11} (SYRK) +* +* U_{12} = U_{12}U_{22}**H (TRMM) +* +* U_{22} = U_{22}U_{22}**H (This subroutine) +* +* +* Compute U_{11} +* + CALL CLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL CHERK('Upper', 'No Transpose', K, N-K, + $ ONE, A(1,K+1), LDA, ONE, A, LDA) +* +* Compute U_{12} +* + CALL CTRMM('Right', 'Upper', 'Conjugate', 'Non-unit', + $ K, N-K, ONE, A(K+1,K+1), LDA, A(1,K+1), LDA) +* +* Compute U_{22} +* + CALL CLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + ELSE +* +* We are computing A = lt(L**H*L). +* +* Break apart L as follows +* |-----------------| +* L = | L_{11} 0 | +* | L_{21} L_{22} | +* |-----------------| +* +* Where +* L_{11}\in\R^{ k\times k} +* L_{21}\in\R^{n-k\times k} l_{22}\in\R^{n-k\times n-k} +* +* and L_{11},L_{22} are lower triangular and L_{21} is rectangular +* +* This gives us our operations as +* |--------------------------| |-----------------| +* lt(L**H * L) = | L_{11}**H L_{21}**H | | L_{11} 0 | +* | 0 L_{22}**H | | L_{21} L_{22} | +* |--------------------------| |-----------------| +* +* Thus we get +* +* L_{11} = L_{11}**H L_{11} + L_{21}**H L_{21} +* L_{21} = L_{22}**H L_{21} +* +* L_{22} = L_{22}**H L_{22} +* +* We break these operations apart as follows +* +* L_{11} = L_{11}**H L_{11} (This subroutine) +* L_{11} = L_{21}**H L_{21} + L_{11} (SYRK) +* +* L_{21} = L_{22}**H L_{21} (TRMM) +* +* L_{22} = L_{22}**H L_{22} (This subroutine) +* +* +* Compute L_{11} +* + CALL CLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL CHERK('Lower', 'Conjugate', K, N-K, + $ ONE, A(K+1,1), LDA, ONE, A, LDA) +* +* Compute L_{21} +* + CALL CTRMM('Left', 'Lower', 'Conjugate', 'Non-Unit', + $ N-K, K, ONE, A(K+1,K+1), LDA, A(K+1,1), LDA) +* +* Compute L_{22} +* + CALL CLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + END IF + END SUBROUTINE diff --git a/SRC/dlauum.f b/SRC/dlauum.f index 2321f238f..014358d8f 100644 --- a/SRC/dlauum.f +++ b/SRC/dlauum.f @@ -127,7 +127,7 @@ SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO ) EXTERNAL LSAME, ILAENV * .. * .. External Subroutines .. - EXTERNAL DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA + EXTERNAL DLAUUM_BLOCKED, DLAUUM_RECURSIVE * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -155,60 +155,13 @@ SUBROUTINE DLAUUM( UPLO, N, A, LDA, INFO ) IF( N.EQ.0 ) $ RETURN * -* Determine the block size for this environment. -* - NB = ILAENV( 1, 'DLAUUM', UPLO, N, -1, -1, -1 ) -* - IF( NB.LE.1 .OR. NB.GE.N ) THEN -* -* Use unblocked code -* - CALL DLAUU2( UPLO, N, A, LDA, INFO ) - ELSE -* -* Use blocked code -* - IF( UPPER ) THEN -* -* Compute the product U * U**T. -* - DO 10 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', - $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), - $ LDA ) - CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL DGEMM( 'No transpose', 'Transpose', I-1, IB, - $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, - $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) - CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, - $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), - $ LDA ) - END IF - 10 CONTINUE - ELSE -* -* Compute the product L**T * L. -* - DO 20 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', - $ IB, - $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) - CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL DGEMM( 'Transpose', 'No transpose', IB, I-1, - $ N-I-IB+1, ONE, A( I+IB, I ), LDA, - $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) - CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, - $ ONE, - $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) - END IF - 20 CONTINUE - END IF - END IF +* Here we dispatch to whatever is more efficient in a particular environment +* We are defaulting to recursive, but if you want to use the blocked variant +* Comment out the line starting with `CALL DLAUUM_RECURSIVE...` +* and uncomment the line starting with `CALL DLAUUM_BLOCKED...` * + CALL DLAUUM_RECURSIVE(UPLO, N, A, LDA, INFO) +* CALL DLAUUM_BLOCKED(UPLO, N, A, LDA, INFO) RETURN * * End of DLAUUM diff --git a/SRC/dlauum_blocked.f b/SRC/dlauum_blocked.f new file mode 100644 index 000000000..7a97c3652 --- /dev/null +++ b/SRC/dlauum_blocked.f @@ -0,0 +1,216 @@ +*> \brief \b DLAUUM_BLOCKED computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download DLAUUM_BLOCKED + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE DLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAUUM_BLOCKED computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE DLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DGEMM, DLAUU2, DSYRK, DTRMM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAUUM_BLOCKED', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'DLAUUM_BLOCKED', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL DLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U**T. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', + $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), + $ LDA ) + CALL DLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL DGEMM( 'No transpose', 'Transpose', I-1, IB, + $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, + $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) + CALL DSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L**T * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', + $ IB, + $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) + CALL DLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL DGEMM( 'Transpose', 'No transpose', IB, I-1, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) + CALL DSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, + $ ONE, + $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of DLAUUM_BLOCKED +* + END diff --git a/SRC/dlauum_recursive.f b/SRC/dlauum_recursive.f new file mode 100644 index 000000000..80cca24a2 --- /dev/null +++ b/SRC/dlauum_recursive.f @@ -0,0 +1,281 @@ +*> \brief \b DLAUUM_RECURSIVE computes the product UUH or LHL, where U and L are upper or lower triangular matrices (recursive algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download DLAUUM_RECURSIVE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE DLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLAUUM_RECURSIVE computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is DOUBLE PRECISION array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE DLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER K, NX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL DSYRK, DTRMM, DLAUU2 +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAUUM_RECURSIVE', -INFO ) + RETURN + END IF +* +* Early termination criteria +* + IF( N.EQ.0 ) THEN + RETURN + END IF +* +* Base Case +* + IF( N.EQ.1 ) THEN + A(1,1) = A(1,1) * A(1,1) + RETURN + END IF +* +* Determine crossover point for when to bail to level2 +* + NX = ILAENV(3, "DLAUUM_RECURSIVE", UPLO, N, -1, -1, -1) + IF( K.LT.NX ) THEN + CALL DLAUU2(UPLO, N, A, LDA, INFO) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + K = N/2 + IF( UPPER ) THEN +* +* We are computing A = ut(U*U**H). +* +* Break apart U as follows +* |-----------------| +* U = | U_{11} U_{12} | +* | 0 U_{22} | +* |-----------------| +* +* Where +* U_{11}\in\R^{k\times k} U_{12}\in\R^{ k\times n-k} +* U_{22}\in\R^{n-k\times n-k} +* +* and U_{11},U_{22} are upper triangular and U_{12} is rectangular +* +* This gives us our operations as +* |--------------------| |------------------------| +* ut(U * U**H) = | U_{11} U_{12} | | U_{11}**H 0 | +* | 0 U_{22} | | U_{12}**H U_{22}**H | +* |--------------------| |------------------------| +* +* Thus we get +* +* U_{11} = U_{11}U_{11}**H + U_{12}U_{12}**H +* U_{12} = U_{12}U_{22}**H +* +* U_{22} = U_{22}U_{22}**H +* +* We break these operations apart as follows +* +* U_{11} = U_{11}U_{11}**H (This subroutine) +* U_{11} = U_{12}U_{12}**H + U_{11} (SYRK) +* +* U_{12} = U_{12}U_{22}**H (TRMM) +* +* U_{22} = U_{22}U_{22}**H (This subroutine) +* +* +* Compute U_{11} +* + CALL DLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL DSYRK('Upper', 'No Transpose', K, N-K, + $ ONE, A(1,K+1), LDA, ONE, A, LDA) +* +* Compute U_{12} +* + CALL DTRMM('Right', 'Upper', 'Transpose', 'Non-unit', + $ K, N-K, ONE, A(K+1,K+1), LDA, A(1,K+1), LDA) +* +* Compute U_{22} +* + CALL DLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + ELSE +* +* We are computing A = lt(L**H*L). +* +* Break apart L as follows +* |-----------------| +* L = | L_{11} 0 | +* | L_{21} L_{22} | +* |-----------------| +* +* Where +* L_{11}\in\R^{ k\times k} +* L_{21}\in\R^{n-k\times k} l_{22}\in\R^{n-k\times n-k} +* +* and L_{11},L_{22} are lower triangular and L_{21} is rectangular +* +* This gives us our operations as +* |--------------------------| |-----------------| +* lt(L**H * L) = | L_{11}**H L_{21}**H | | L_{11} 0 | +* | 0 L_{22}**H | | L_{21} L_{22} | +* |--------------------------| |-----------------| +* +* Thus we get +* +* L_{11} = L_{11}**H L_{11} + L_{21}**H L_{21} +* L_{21} = L_{22}**H L_{21} +* +* L_{22} = L_{22}**H L_{22} +* +* We break these operations apart as follows +* +* L_{11} = L_{11}**H L_{11} (This subroutine) +* L_{11} = L_{21}**H L_{21} + L_{11} (SYRK) +* +* L_{21} = L_{22}**H L_{21} (TRMM) +* +* L_{22} = L_{22}**H L_{22} (This subroutine) +* +* +* Compute L_{11} +* + CALL DLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL DSYRK('Lower', 'Transpose', K, N-K, + $ ONE, A(K+1,1), LDA, ONE, A, LDA) +* +* Compute L_{21} +* + CALL DTRMM('Left', 'Lower', 'Transpose', 'Non-Unit', + $ N-K, K, ONE, A(K+1,K+1), LDA, A(K+1,1), LDA) +* +* Compute L_{22} +* + CALL DLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + END IF + END SUBROUTINE diff --git a/SRC/ilaenv.f b/SRC/ilaenv.f index c0371c0b0..2622d44cc 100644 --- a/SRC/ilaenv.f +++ b/SRC/ilaenv.f @@ -679,6 +679,12 @@ INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 ) ELSE IF( C2.EQ.'LA' ) THEN IF( C3.EQ.'RFT' ) THEN NX = 64 + ELSE IF( C3.EQ.'UUM' ) THEN +* +* Note that this matches on *lauum_recursive, which is the +* only place we currently use this +* + NX = 64 END IF END IF ILAENV = NX diff --git a/SRC/slauum.f b/SRC/slauum.f index f4e1d9df6..d92c68c3d 100644 --- a/SRC/slauum.f +++ b/SRC/slauum.f @@ -155,60 +155,13 @@ SUBROUTINE SLAUUM( UPLO, N, A, LDA, INFO ) IF( N.EQ.0 ) $ RETURN * -* Determine the block size for this environment. -* - NB = ILAENV( 1, 'SLAUUM', UPLO, N, -1, -1, -1 ) -* - IF( NB.LE.1 .OR. NB.GE.N ) THEN -* -* Use unblocked code -* - CALL SLAUU2( UPLO, N, A, LDA, INFO ) - ELSE -* -* Use blocked code -* - IF( UPPER ) THEN -* -* Compute the product U * U**T. -* - DO 10 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL STRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', - $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), - $ LDA ) - CALL SLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL SGEMM( 'No transpose', 'Transpose', I-1, IB, - $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, - $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) - CALL SSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, - $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), - $ LDA ) - END IF - 10 CONTINUE - ELSE -* -* Compute the product L**T * L. -* - DO 20 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL STRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', - $ IB, - $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) - CALL SLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL SGEMM( 'Transpose', 'No transpose', IB, I-1, - $ N-I-IB+1, ONE, A( I+IB, I ), LDA, - $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) - CALL SSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, - $ ONE, - $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) - END IF - 20 CONTINUE - END IF - END IF +* Here we dispatch to whatever is more efficient in a particular environment +* We are defaulting to recursive, but if you want to use the blocked variant +* Comment out the line starting with `CALL SLAUUM_RECURSIVE...` +* and uncomment the line starting with `CALL SLAUUM_BLOCKED...` * + CALL SLAUUM_RECURSIVE(UPLO, N, A, LDA, INFO) +* CALL SLAUUM_BLOCKED(UPLO, N, A, LDA, INFO) RETURN * * End of SLAUUM diff --git a/SRC/slauum_blocked.f b/SRC/slauum_blocked.f new file mode 100644 index 000000000..a66a3d0dd --- /dev/null +++ b/SRC/slauum_blocked.f @@ -0,0 +1,216 @@ +*> \brief \b SLAUUM_BLOCKED computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download SLAUUM_BLOCKED + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE SLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* REAL A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLAUUM_BLOCKED computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE SLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL SGEMM, SLAUU2, SSYRK, STRMM, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SLAUUM_BLOCKED', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'SLAUUM_BLOCKED', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL SLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U**T. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL STRMM( 'Right', 'Upper', 'Transpose', 'Non-unit', + $ I-1, IB, ONE, A( I, I ), LDA, A( 1, I ), + $ LDA ) + CALL SLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL SGEMM( 'No transpose', 'Transpose', I-1, IB, + $ N-I-IB+1, ONE, A( 1, I+IB ), LDA, + $ A( I, I+IB ), LDA, ONE, A( 1, I ), LDA ) + CALL SSYRK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L**T * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL STRMM( 'Left', 'Lower', 'Transpose', 'Non-unit', + $ IB, + $ I-1, ONE, A( I, I ), LDA, A( I, 1 ), LDA ) + CALL SLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL SGEMM( 'Transpose', 'No transpose', IB, I-1, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, ONE, A( I, 1 ), LDA ) + CALL SSYRK( 'Lower', 'Transpose', IB, N-I-IB+1, + $ ONE, + $ A( I+IB, I ), LDA, ONE, A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of SLAUUM_BLOCKED +* + END diff --git a/SRC/slauum_recursive.f b/SRC/slauum_recursive.f new file mode 100644 index 000000000..6ac7e218b --- /dev/null +++ b/SRC/slauum_recursive.f @@ -0,0 +1,281 @@ +*> \brief \b SLAUUM_RECURSIVE computes the product UUH or LHL, where U and L are upper or lower triangular matrices (recursive algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download SLAUUM_RECURSIVE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE SLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* REAL A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLAUUM_RECURSIVE computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is REAL array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE SLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + REAL A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER K, NX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL SSYRK, STRMM, SLAUU2 +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAUUM_RECURSIVE', -INFO ) + RETURN + END IF +* +* Early termination criteria +* + IF( N.EQ.0 ) THEN + RETURN + END IF +* +* Base Case +* + IF( N.EQ.1 ) THEN + A(1,1) = A(1,1) * A(1,1) + RETURN + END IF +* +* Determine crossover point for when to bail to level2 +* + NX = ILAENV(3, "DLAUUM_RECURSIVE", UPLO, N, -1, -1, -1) + IF( K.LT.NX ) THEN + CALL SLAUU2(UPLO, N, A, LDA, INFO) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + K = N/2 + IF( UPPER ) THEN +* +* We are computing A = ut(U*U**H). +* +* Break apart U as follows +* |-----------------| +* U = | U_{11} U_{12} | +* | 0 U_{22} | +* |-----------------| +* +* Where +* U_{11}\in\R^{k\times k} U_{12}\in\R^{ k\times n-k} +* U_{22}\in\R^{n-k\times n-k} +* +* and U_{11},U_{22} are upper triangular and U_{12} is rectangular +* +* This gives us our operations as +* |--------------------| |------------------------| +* ut(U * U**H) = | U_{11} U_{12} | | U_{11}**H 0 | +* | 0 U_{22} | | U_{12}**H U_{22}**H | +* |--------------------| |------------------------| +* +* Thus we get +* +* U_{11} = U_{11}U_{11}**H + U_{12}U_{12}**H +* U_{12} = U_{12}U_{22}**H +* +* U_{22} = U_{22}U_{22}**H +* +* We break these operations apart as follows +* +* U_{11} = U_{11}U_{11}**H (This subroutine) +* U_{11} = U_{12}U_{12}**H + U_{11} (SYRK) +* +* U_{12} = U_{12}U_{22}**H (TRMM) +* +* U_{22} = U_{22}U_{22}**H (This subroutine) +* +* +* Compute U_{11} +* + CALL SLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL SSYRK('Upper', 'No Transpose', K, N-K, + $ ONE, A(1,K+1), LDA, ONE, A, LDA) +* +* Compute U_{12} +* + CALL STRMM('Right', 'Upper', 'Transpose', 'Non-unit', + $ K, N-K, ONE, A(K+1,K+1), LDA, A(1,K+1), LDA) +* +* Compute U_{22} +* + CALL SLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + ELSE +* +* We are computing A = lt(L**H*L). +* +* Break apart L as follows +* |-----------------| +* L = | L_{11} 0 | +* | L_{21} L_{22} | +* |-----------------| +* +* Where +* L_{11}\in\R^{ k\times k} +* L_{21}\in\R^{n-k\times k} l_{22}\in\R^{n-k\times n-k} +* +* and L_{11},L_{22} are lower triangular and L_{21} is rectangular +* +* This gives us our operations as +* |--------------------------| |-----------------| +* lt(L**H * L) = | L_{11}**H L_{21}**H | | L_{11} 0 | +* | 0 L_{22}**H | | L_{21} L_{22} | +* |--------------------------| |-----------------| +* +* Thus we get +* +* L_{11} = L_{11}**H L_{11} + L_{21}**H L_{21} +* L_{21} = L_{22}**H L_{21} +* +* L_{22} = L_{22}**H L_{22} +* +* We break these operations apart as follows +* +* L_{11} = L_{11}**H L_{11} (This subroutine) +* L_{11} = L_{21}**H L_{21} + L_{11} (SYRK) +* +* L_{21} = L_{22}**H L_{21} (TRMM) +* +* L_{22} = L_{22}**H L_{22} (This subroutine) +* +* +* Compute L_{11} +* + CALL SLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL SSYRK('Lower', 'Transpose', K, N-K, + $ ONE, A(K+1,1), LDA, ONE, A, LDA) +* +* Compute L_{21} +* + CALL STRMM('Left', 'Lower', 'Transpose', 'Non-Unit', + $ N-K, K, ONE, A(K+1,K+1), LDA, A(K+1,1), LDA) +* +* Compute L_{22} +* + CALL SLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + END IF + END SUBROUTINE diff --git a/SRC/zlauum.f b/SRC/zlauum.f index d7c13fc37..6eeab4ef5 100644 --- a/SRC/zlauum.f +++ b/SRC/zlauum.f @@ -157,62 +157,13 @@ SUBROUTINE ZLAUUM( UPLO, N, A, LDA, INFO ) IF( N.EQ.0 ) $ RETURN * -* Determine the block size for this environment. -* - NB = ILAENV( 1, 'ZLAUUM', UPLO, N, -1, -1, -1 ) -* - IF( NB.LE.1 .OR. NB.GE.N ) THEN -* -* Use unblocked code -* - CALL ZLAUU2( UPLO, N, A, LDA, INFO ) - ELSE -* -* Use blocked code -* - IF( UPPER ) THEN -* -* Compute the product U * U**H. -* - DO 10 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', - $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, - $ A( 1, I ), LDA ) - CALL ZLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL ZGEMM( 'No transpose', 'Conjugate transpose', - $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), - $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), - $ LDA ) - CALL ZHERK( 'Upper', 'No transpose', IB, N-I-IB+1, - $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), - $ LDA ) - END IF - 10 CONTINUE - ELSE -* -* Compute the product L**H * L. -* - DO 20 I = 1, N, NB - IB = MIN( NB, N-I+1 ) - CALL ZTRMM( 'Left', 'Lower', 'Conjugate transpose', - $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, - $ A( I, 1 ), LDA ) - CALL ZLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) - IF( I+IB.LE.N ) THEN - CALL ZGEMM( 'Conjugate transpose', 'No transpose', - $ IB, - $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, - $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) - CALL ZHERK( 'Lower', 'Conjugate transpose', IB, - $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, - $ A( I, I ), LDA ) - END IF - 20 CONTINUE - END IF - END IF +* Here we dispatch to whatever is more efficient in a particular environment +* We are defaulting to recursive, but if you want to use the blocked variant +* Comment out the line starting with `CALL ZLAUUM_RECURSIVE...` +* and uncomment the line starting with `CALL ZLAUUM_BLOCKED...` * + CALL ZLAUUM_RECURSIVE(UPLO, N, A, LDA, INFO) +* CALL ZLAUUM_BLOCKED(UPLO, N, A, LDA, INFO) RETURN * * End of ZLAUUM diff --git a/SRC/zlauum_blocked.f b/SRC/zlauum_blocked.f new file mode 100644 index 000000000..28fa7cbb5 --- /dev/null +++ b/SRC/zlauum_blocked.f @@ -0,0 +1,220 @@ +*> \brief \b ZLAUUM_BLOCKED computes the product UUH or LHL, where U and L are upper or lower triangular matrices (blocked algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download ZLAUUM_BLOCKED + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE ZLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX*16 A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLAUUM_BLOCKED computes the product U * U**H or L**H * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**H; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**H * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE ZLAUUM_BLOCKED( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) + COMPLEX*16 CONE + PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER I, IB, NB +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL XERBLA, ZGEMM, ZHERK, ZLAUU2, ZTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX, MIN +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'ZLAUUM_BLOCKED', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* +* Determine the block size for this environment. +* + NB = ILAENV( 1, 'ZLAUUM_BLOCKED', UPLO, N, -1, -1, -1 ) +* + IF( NB.LE.1 .OR. NB.GE.N ) THEN +* +* Use unblocked code +* + CALL ZLAUU2( UPLO, N, A, LDA, INFO ) + ELSE +* +* Use blocked code +* + IF( UPPER ) THEN +* +* Compute the product U * U**H. +* + DO 10 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Non-unit', I-1, IB, CONE, A( I, I ), LDA, + $ A( 1, I ), LDA ) + CALL ZLAUU2( 'Upper', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ I-1, IB, N-I-IB+1, CONE, A( 1, I+IB ), + $ LDA, A( I, I+IB ), LDA, CONE, A( 1, I ), + $ LDA ) + CALL ZHERK( 'Upper', 'No transpose', IB, N-I-IB+1, + $ ONE, A( I, I+IB ), LDA, ONE, A( I, I ), + $ LDA ) + END IF + 10 CONTINUE + ELSE +* +* Compute the product L**H * L. +* + DO 20 I = 1, N, NB + IB = MIN( NB, N-I+1 ) + CALL ZTRMM( 'Left', 'Lower', 'Conjugate transpose', + $ 'Non-unit', IB, I-1, CONE, A( I, I ), LDA, + $ A( I, 1 ), LDA ) + CALL ZLAUU2( 'Lower', IB, A( I, I ), LDA, INFO ) + IF( I+IB.LE.N ) THEN + CALL ZGEMM( 'Conjugate transpose', 'No transpose', + $ IB, + $ I-1, N-I-IB+1, CONE, A( I+IB, I ), LDA, + $ A( I+IB, 1 ), LDA, CONE, A( I, 1 ), LDA ) + CALL ZHERK( 'Lower', 'Conjugate transpose', IB, + $ N-I-IB+1, ONE, A( I+IB, I ), LDA, ONE, + $ A( I, I ), LDA ) + END IF + 20 CONTINUE + END IF + END IF +* + RETURN +* +* End of ZLAUUM_BLOCKED +* + END diff --git a/SRC/zlauum_recursive.f b/SRC/zlauum_recursive.f new file mode 100644 index 000000000..e3d58295f --- /dev/null +++ b/SRC/zlauum_recursive.f @@ -0,0 +1,281 @@ +*> \brief \b ZLAUUM_RECURSIVE computes the product UUH or LHL, where U and L are upper or lower triangular matrices (recursive algorithm). +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> Download ZLAUUM_RECURSIVE + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +* +* Definition: +* =========== +* +* SUBROUTINE ZLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) +* +* .. Scalar Arguments .. +* CHARACTER UPLO +* INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. +* COMPLEX*16 A( LDA, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLAUUM_RECURSIVE computes the product U * U**T or L**T * L, where the triangular +*> factor U or L is stored in the upper or lower triangular part of +*> the array A. +*> +*> If UPLO = 'U' or 'u' then the upper triangle of the result is stored, +*> overwriting the factor U in A. +*> If UPLO = 'L' or 'l' then the lower triangle of the result is stored, +*> overwriting the factor L in A. +*> +*> This is the blocked form of the algorithm, calling Level 3 BLAS. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] UPLO +*> \verbatim +*> UPLO is CHARACTER*1 +*> Specifies whether the triangular factor stored in the array A +*> is upper or lower triangular: +*> = 'U': Upper triangular +*> = 'L': Lower triangular +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the triangular factor U or L. N >= 0. +*> \endverbatim +*> +*> \param[in,out] A +*> \verbatim +*> A is COMPLEX*16 array, dimension (LDA,N) +*> On entry, the triangular factor U or L. +*> On exit, if UPLO = 'U', the upper triangle of A is +*> overwritten with the upper triangle of the product U * U**T; +*> if UPLO = 'L', the lower triangle of A is overwritten with +*> the lower triangle of the product L**T * L. +*> \endverbatim +*> +*> \param[in] LDA +*> \verbatim +*> LDA is INTEGER +*> The leading dimension of the array A. LDA >= max(1,N). +*> \endverbatim +*> +*> \param[out] INFO +*> \verbatim +*> INFO is INTEGER +*> = 0: successful exit +*> < 0: if INFO = -k, the k-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. +* +*> \ingroup lauum +* +* ===================================================================== + SUBROUTINE ZLAUUM_RECURSIVE( UPLO, N, A, LDA, INFO ) + IMPLICIT NONE +* +* -- LAPACK auxiliary routine -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* +* .. Scalar Arguments .. + CHARACTER UPLO + INTEGER INFO, LDA, N +* .. +* .. Array Arguments .. + COMPLEX*16 A( LDA, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL UPPER + INTEGER K, NX +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAENV + EXTERNAL LSAME, ILAENV +* .. +* .. External Subroutines .. + EXTERNAL ZHERK, ZTRMM, ZLAUU2 +* .. +* .. 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 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DLAUUM_RECURSIVE', -INFO ) + RETURN + END IF +* +* Early termination criteria +* + IF( N.EQ.0 ) THEN + RETURN + END IF +* +* Base Case +* + IF( N.EQ.1 ) THEN + A(1,1) = A(1,1) * A(1,1) + RETURN + END IF +* +* Determine crossover point for when to bail to level2 +* + NX = ILAENV(3, "DLAUUM_RECURSIVE", UPLO, N, -1, -1, -1) + IF( K.LT.NX ) THEN + CALL ZLAUU2(UPLO, N, A, LDA, INFO) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + K = N/2 + IF( UPPER ) THEN +* +* We are computing A = ut(U*U**H). +* +* Break apart U as follows +* |-----------------| +* U = | U_{11} U_{12} | +* | 0 U_{22} | +* |-----------------| +* +* Where +* U_{11}\in\R^{k\times k} U_{12}\in\R^{ k\times n-k} +* U_{22}\in\R^{n-k\times n-k} +* +* and U_{11},U_{22} are upper triangular and U_{12} is rectangular +* +* This gives us our operations as +* |--------------------| |------------------------| +* ut(U * U**H) = | U_{11} U_{12} | | U_{11}**H 0 | +* | 0 U_{22} | | U_{12}**H U_{22}**H | +* |--------------------| |------------------------| +* +* Thus we get +* +* U_{11} = U_{11}U_{11}**H + U_{12}U_{12}**H +* U_{12} = U_{12}U_{22}**H +* +* U_{22} = U_{22}U_{22}**H +* +* We break these operations apart as follows +* +* U_{11} = U_{11}U_{11}**H (This subroutine) +* U_{11} = U_{12}U_{12}**H + U_{11} (SYRK) +* +* U_{12} = U_{12}U_{22}**H (TRMM) +* +* U_{22} = U_{22}U_{22}**H (This subroutine) +* +* +* Compute U_{11} +* + CALL ZLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL ZHERK('Upper', 'No Transpose', K, N-K, + $ ONE, A(1,K+1), LDA, ONE, A, LDA) +* +* Compute U_{12} +* + CALL ZTRMM('Right', 'Upper', 'Conjugate', 'Non-unit', + $ K, N-K, ONE, A(K+1,K+1), LDA, A(1,K+1), LDA) +* +* Compute U_{22} +* + CALL ZLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + ELSE +* +* We are computing A = lt(L**H*L). +* +* Break apart L as follows +* |-----------------| +* L = | L_{11} 0 | +* | L_{21} L_{22} | +* |-----------------| +* +* Where +* L_{11}\in\R^{ k\times k} +* L_{21}\in\R^{n-k\times k} l_{22}\in\R^{n-k\times n-k} +* +* and L_{11},L_{22} are lower triangular and L_{21} is rectangular +* +* This gives us our operations as +* |--------------------------| |-----------------| +* lt(L**H * L) = | L_{11}**H L_{21}**H | | L_{11} 0 | +* | 0 L_{22}**H | | L_{21} L_{22} | +* |--------------------------| |-----------------| +* +* Thus we get +* +* L_{11} = L_{11}**H L_{11} + L_{21}**H L_{21} +* L_{21} = L_{22}**H L_{21} +* +* L_{22} = L_{22}**H L_{22} +* +* We break these operations apart as follows +* +* L_{11} = L_{11}**H L_{11} (This subroutine) +* L_{11} = L_{21}**H L_{21} + L_{11} (SYRK) +* +* L_{21} = L_{22}**H L_{21} (TRMM) +* +* L_{22} = L_{22}**H L_{22} (This subroutine) +* +* +* Compute L_{11} +* + CALL ZLAUUM_RECURSIVE(UPLO, K, A, LDA, INFO) + CALL ZHERK('Lower', 'Conjugate', K, N-K, + $ ONE, A(K+1,1), LDA, ONE, A, LDA) +* +* Compute L_{21} +* + CALL ZTRMM('Left', 'Lower', 'Conjugate', 'Non-Unit', + $ N-K, K, ONE, A(K+1,K+1), LDA, A(K+1,1), LDA) +* +* Compute L_{22} +* + CALL ZLAUUM_RECURSIVE(UPLO, N-K, A(K+1,K+1), LDA, INFO) + END IF + END SUBROUTINE