From 9cf77ce1d80cf17aa79c5da95b578ee2a4490152 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Mon, 12 Nov 2012 15:20:37 +0100 Subject: Add support for Sparse QR factorization --- lapack/CMakeLists.txt | 8 +- lapack/clarf.f | 232 +++++++++++++++ lapack/clarfb.f | 771 +++++++++++++++++++++++++++++++++++++++++++++++++ lapack/clarfg.f | 203 +++++++++++++ lapack/clarft.f | 328 +++++++++++++++++++++ lapack/dlarf.f | 227 +++++++++++++++ lapack/dlarfb.f | 762 +++++++++++++++++++++++++++++++++++++++++++++++++ lapack/dlarfg.f | 196 +++++++++++++ lapack/dlarft.f | 326 +++++++++++++++++++++ lapack/slarf.f | 227 +++++++++++++++ lapack/slarfb.f | 763 +++++++++++++++++++++++++++++++++++++++++++++++++ lapack/slarfg.f | 196 +++++++++++++ lapack/slarft.f | 326 +++++++++++++++++++++ lapack/zlarf.f | 232 +++++++++++++++ lapack/zlarfb.f | 774 ++++++++++++++++++++++++++++++++++++++++++++++++++ lapack/zlarfg.f | 203 +++++++++++++ lapack/zlarft.f | 327 +++++++++++++++++++++ 17 files changed, 6100 insertions(+), 1 deletion(-) create mode 100644 lapack/clarf.f create mode 100644 lapack/clarfb.f create mode 100644 lapack/clarfg.f create mode 100644 lapack/clarft.f create mode 100644 lapack/dlarf.f create mode 100644 lapack/dlarfb.f create mode 100644 lapack/dlarfg.f create mode 100644 lapack/dlarft.f create mode 100644 lapack/slarf.f create mode 100644 lapack/slarfb.f create mode 100644 lapack/slarfg.f create mode 100644 lapack/slarft.f create mode 100644 lapack/zlarf.f create mode 100644 lapack/zlarfb.f create mode 100644 lapack/zlarfg.f create mode 100644 lapack/zlarft.f (limited to 'lapack') diff --git a/lapack/CMakeLists.txt b/lapack/CMakeLists.txt index 062845a3f..174e77d8d 100644 --- a/lapack/CMakeLists.txt +++ b/lapack/CMakeLists.txt @@ -18,6 +18,13 @@ single.cpp double.cpp complex_single.cpp complex_double.cpp ../blas/xerbla.cpp if(EIGEN_Fortran_COMPILER_WORKS) +set(EigenLapack_SRCS +slarft.f dlarft.f clarft.f zlarft.f +slarfb.f dlarfb.f clarfb.f zlarfb.f +slarfg.f dlarfg.f clarfg.f zlarfg.f +slarf.f dlarf.f clarf.f zlarf.f +) + get_filename_component(eigen_full_path_to_reference_to_reference_lapack "./reference/" ABSOLUTE) if(EXISTS ${eigen_full_path_to_reference_to_reference_lapack}) set(EigenLapack_SRCS ${EigenLapack_SRCS} @@ -26,7 +33,6 @@ set(EigenLapack_SRCS ${EigenLapack_SRCS} # reference/dgetrf.f reference/cgetrf.f reference/sgetrf.f reference/zgetrf.f # reference/cgetrs.f reference/dgetrs.f reference/sgetrs.f reference/zgetrs.f # reference/dsyev.f reference/ssyev.f - reference/dlamch.f reference/ilaver.f reference/lsame.f reference/slamch.f reference/second_NONE.f reference/dsecnd_NONE.f reference/cbdsqr.f reference/ctbrfs.f reference/dorml2.f reference/sla_porfsx_extended.f reference/zggglm.f reference/cgbbrd.f reference/ctbtrs.f reference/dormlq.f reference/sla_porpvgrw.f reference/zgghrd.f diff --git a/lapack/clarf.f b/lapack/clarf.f new file mode 100644 index 000000000..ca0328fb5 --- /dev/null +++ b/lapack/clarf.f @@ -0,0 +1,232 @@ +*> \brief \b CLARF +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLARF + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* COMPLEX TAU +* .. +* .. Array Arguments .. +* COMPLEX C( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLARF applies a complex elementary reflector H to a complex M-by-N +*> matrix C, from either the left or the right. H is represented in the +*> form +*> +*> H = I - tau * v * v**H +*> +*> where tau is a complex scalar and v is a complex vector. +*> +*> If tau = 0, then H is taken to be the unit matrix. +*> +*> To apply H**H (the conjugate transpose of H), supply conjg(tau) instead +*> tau. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form H * C +*> = 'R': form C * H +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of H. V is not used if +*> TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is COMPLEX +*> The value tau in the representation of H. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is COMPLEX array, dimension (LDC,N) +*> On entry, the M-by-N matrix C. +*> On exit, C is overwritten by the matrix H * C if SIDE = 'L', +*> or C * H if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension +*> (N) if SIDE = 'L' +*> or (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexOTHERauxiliary +* +* ===================================================================== + SUBROUTINE CLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + COMPLEX TAU +* .. +* .. Array Arguments .. + COMPLEX C( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE, ZERO + PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), + $ ZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL APPLYLEFT + INTEGER I, LASTV, LASTC +* .. +* .. External Subroutines .. + EXTERNAL CGEMV, CGERC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILACLR, ILACLC + EXTERNAL LSAME, ILACLR, ILACLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 0 + LASTC = 0 + IF( TAU.NE.ZERO ) THEN +! Set up variables for scanning V. LASTV begins pointing to the end +! of V. + IF( APPLYLEFT ) THEN + LASTV = M + ELSE + LASTV = N + END IF + IF( INCV.GT.0 ) THEN + I = 1 + (LASTV-1) * INCV + ELSE + I = 1 + END IF +! Look for the last non-zero row in V. + DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) + LASTV = LASTV - 1 + I = I - INCV + END DO + IF( APPLYLEFT ) THEN +! Scan for the last non-zero column in C(1:lastv,:). + LASTC = ILACLC(LASTV, N, C, LDC) + ELSE +! Scan for the last non-zero row in C(:,1:lastv). + LASTC = ILACLR(M, LASTV, C, LDC) + END IF + END IF +! Note that lastc.eq.0 renders the BLAS operations null; no special +! case is needed at this level. + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**H * v(1:lastv,1) +* + CALL CGEMV( 'Conjugate transpose', LASTV, LASTC, ONE, + $ C, LDC, V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**H +* + CALL CGERC( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC ) + END IF + ELSE +* +* Form C * H +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + CALL CGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC, + $ V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**H +* + CALL CGERC( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC ) + END IF + END IF + RETURN +* +* End of CLARF +* + END diff --git a/lapack/clarfb.f b/lapack/clarfb.f new file mode 100644 index 000000000..40bbdf487 --- /dev/null +++ b/lapack/clarfb.f @@ -0,0 +1,771 @@ +*> \brief \b CLARFB +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLARFB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, +* T, LDT, C, LDC, WORK, LDWORK ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, SIDE, STOREV, TRANS +* INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. +* COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), +* $ WORK( LDWORK, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLARFB applies a complex block reflector H or its transpose H**H to a +*> complex M-by-N matrix C, from either the left or the right. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': apply H or H**H from the Left +*> = 'R': apply H or H**H from the Right +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> = 'N': apply H (No transpose) +*> = 'C': apply H**H (Conjugate transpose) +*> \endverbatim +*> +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Indicates how H is formed from a product of elementary +*> reflectors +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Indicates how the vectors which define the elementary +*> reflectors are stored: +*> = 'C': Columnwise +*> = 'R': Rowwise +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the matrix T (= the number of elementary +*> reflectors whose product defines the block reflector). +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,M) if STOREV = 'R' and SIDE = 'L' +*> (LDV,N) if STOREV = 'R' and SIDE = 'R' +*> The matrix V. See Further Details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +*> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +*> if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] T +*> \verbatim +*> T is COMPLEX array, dimension (LDT,K) +*> The triangular K-by-K matrix T in the representation of the +*> block reflector. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is COMPLEX array, dimension (LDC,N) +*> On entry, the M-by-N matrix C. +*> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX array, dimension (LDWORK,K) +*> \endverbatim +*> +*> \param[in] LDWORK +*> \verbatim +*> LDWORK is INTEGER +*> The leading dimension of the array WORK. +*> If SIDE = 'L', LDWORK >= max(1,N); +*> if SIDE = 'R', LDWORK >= max(1,M). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored; the corresponding +*> array elements are modified but restored on exit. The rest of the +*> array is not used. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE CLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, + $ T, LDT, C, LDC, WORK, LDWORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX C( LDC, * ), T( LDT, * ), V( LDV, * ), + $ WORK( LDWORK, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE + PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + CHARACTER TRANST + INTEGER I, J, LASTV, LASTC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILACLR, ILACLC + EXTERNAL LSAME, ILACLR, ILACLC +* .. +* .. External Subroutines .. + EXTERNAL CCOPY, CGEMM, CLACGV, CTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 ) + $ RETURN +* + IF( LSAME( TRANS, 'N' ) ) THEN + TRANST = 'C' + ELSE + TRANST = 'N' + END IF +* + IF( LSAME( STOREV, 'C' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 ) (first K rows) +* ( V2 ) +* where V1 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) + LASTC = ILACLC( LASTV, N, C, LDC ) +* +* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) +* +* W := C1**H +* + DO 10 J = 1, K + CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + CALL CLACGV( LASTC, WORK( 1, J ), 1 ) + 10 CONTINUE +* +* W := W * V1 +* + CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**H *V2 +* + CALL CGEMM( 'Conjugate transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC, + $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**H +* + IF( M.GT.K ) THEN +* +* C2 := C2 - V2 * W**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTV-K, LASTC, K, -ONE, V( K+1, 1 ), LDV, + $ WORK, LDWORK, ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1**H +* + CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**H +* + DO 30 J = 1, K + DO 20 I = 1, LASTC + C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) + 20 CONTINUE + 30 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) + LASTC = ILACLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C1 +* + DO 40 J = 1, K + CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 40 CONTINUE +* +* W := W * V1 +* + CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2 +* + CALL CGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**H +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1**H +* + CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 60 J = 1, K + DO 50 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 50 CONTINUE + 60 CONTINUE + END IF +* + ELSE +* +* Let V = ( V1 ) +* ( V2 ) (last K rows) +* where V2 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILACLR( M, K, V, LDV ) ) + LASTC = ILACLC( LASTV, N, C, LDC ) +* +* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) +* +* W := C2**H +* + DO 70 J = 1, K + CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + CALL CLACGV( LASTC, WORK( 1, J ), 1 ) + 70 CONTINUE +* +* W := W * V2 +* + CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**H*V1 +* + CALL CGEMM( 'Conjugate transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1 * W**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**H +* + CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**H +* + DO 90 J = 1, K + DO 80 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - + $ CONJG( WORK( I, J ) ) + 80 CONTINUE + 90 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILACLR( N, K, V, LDV ) ) + LASTC = ILACLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C2 +* + DO 100 J = 1, K + CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 100 CONTINUE +* +* W := W * V2 +* + CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1 +* + CALL CGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**H +* + CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W +* + DO 120 J = 1, K + DO 110 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) + $ - WORK( I, J ) + 110 CONTINUE + 120 CONTINUE + END IF + END IF +* + ELSE IF( LSAME( STOREV, 'R' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 V2 ) (V1: first K columns) +* where V1 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) + LASTC = ILACLC( LASTV, N, C, LDC ) +* +* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) +* +* W := C1**H +* + DO 130 J = 1, K + CALL CCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + CALL CLACGV( LASTC, WORK( 1, J ), 1 ) + 130 CONTINUE +* +* W := W * V1**H +* + CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**H*V2**H +* + CALL CGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL CTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**H * W**H +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2**H * W**H +* + CALL CGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTV-K, LASTC, K, + $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, + $ ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**H +* + DO 150 J = 1, K + DO 140 I = 1, LASTC + C( J, I ) = C( J, I ) - CONJG( WORK( I, J ) ) + 140 CONTINUE + 150 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) + LASTC = ILACLR( M, LASTV, C, LDC ) +* +* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) +* +* W := C1 +* + DO 160 J = 1, K + CALL CCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 160 CONTINUE +* +* W := W * V1**H +* + CALL CTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC, + $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL CTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2 +* + CALL CGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL CTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 180 J = 1, K + DO 170 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 170 CONTINUE + 180 CONTINUE +* + END IF +* + ELSE +* +* Let V = ( V1 V2 ) (V2: last K columns) +* where V2 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILACLC( K, M, V, LDV ) ) + LASTC = ILACLC( LASTV, N, C, LDC ) +* +* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) +* +* W := C2**H +* + DO 190 J = 1, K + CALL CCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + CALL CLACGV( LASTC, WORK( 1, J ), 1 ) + 190 CONTINUE +* +* W := W * V2**H +* + CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**H * V1**H +* + CALL CGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTC, K, LASTV-K, + $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL CTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**H * W**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1**H * W**H +* + CALL CGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTV-K, LASTC, K, + $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**H +* + DO 210 J = 1, K + DO 200 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - + $ CONJG( WORK( I, J ) ) + 200 CONTINUE + 210 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILACLC( K, N, V, LDV ) ) + LASTC = ILACLR( M, LASTV, C, LDC ) +* +* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) +* +* W := C2 +* + DO 220 J = 1, K + CALL CCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 220 CONTINUE +* +* W := W * V2**H +* + CALL CTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1**H +* + CALL CGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE, + $ WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL CTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1 +* + CALL CGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL CTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 240 J = 1, K + DO 230 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) + $ - WORK( I, J ) + 230 CONTINUE + 240 CONTINUE +* + END IF +* + END IF + END IF +* + RETURN +* +* End of CLARFB +* + END diff --git a/lapack/clarfg.f b/lapack/clarfg.f new file mode 100644 index 000000000..d64f396c3 --- /dev/null +++ b/lapack/clarfg.f @@ -0,0 +1,203 @@ +*> \brief \b CLARFG +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLARFG + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* COMPLEX ALPHA, TAU +* .. +* .. Array Arguments .. +* COMPLEX X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLARFG generates a complex elementary reflector H of order n, such +*> that +*> +*> H**H * ( alpha ) = ( beta ), H**H * H = I. +*> ( x ) ( 0 ) +*> +*> where alpha and beta are scalars, with beta real, and x is an +*> (n-1)-element complex vector. H is represented in the form +*> +*> H = I - tau * ( 1 ) * ( 1 v**H ) , +*> ( v ) +*> +*> where tau is a complex scalar and v is a complex (n-1)-element +*> vector. Note that H is not hermitian. +*> +*> If the elements of x are all zero and alpha is real, then tau = 0 +*> and H is taken to be the unit matrix. +*> +*> Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 . +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the elementary reflector. +*> \endverbatim +*> +*> \param[in,out] ALPHA +*> \verbatim +*> ALPHA is COMPLEX +*> On entry, the value alpha. +*> On exit, it is overwritten with the value beta. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is COMPLEX array, dimension +*> (1+(N-2)*abs(INCX)) +*> On entry, the vector x. +*> On exit, it is overwritten with the vector v. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between elements of X. INCX > 0. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX +*> The value tau. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complexOTHERauxiliary +* +* ===================================================================== + SUBROUTINE CLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + INTEGER INCX, N + COMPLEX ALPHA, TAU +* .. +* .. Array Arguments .. + COMPLEX X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + REAL ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + REAL SCNRM2, SLAMCH, SLAPY3 + COMPLEX CLADIV + EXTERNAL SCNRM2, SLAMCH, SLAPY3, CLADIV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, AIMAG, CMPLX, REAL, SIGN +* .. +* .. External Subroutines .. + EXTERNAL CSCAL, CSSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.0 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = SCNRM2( N-1, X, INCX ) + ALPHR = REAL( ALPHA ) + ALPHI = AIMAG( ALPHA ) +* + IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' ) + RSAFMN = ONE / SAFMIN +* + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + 10 CONTINUE + KNT = KNT + 1 + CALL CSSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHI = ALPHI*RSAFMN + ALPHR = ALPHR*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = SCNRM2( N-1, X, INCX ) + ALPHA = CMPLX( ALPHR, ALPHI ) + BETA = -SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + END IF + TAU = CMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) + ALPHA = CLADIV( CMPLX( ONE ), ALPHA-BETA ) + CALL CSCAL( N-1, ALPHA, X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of CLARFG +* + END diff --git a/lapack/clarft.f b/lapack/clarft.f new file mode 100644 index 000000000..981447f77 --- /dev/null +++ b/lapack/clarft.f @@ -0,0 +1,328 @@ +*> \brief \b CLARFT +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download CLARFT + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, STOREV +* INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. +* COMPLEX T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> CLARFT forms the triangular factor T of a complex block reflector H +*> of order n, which is defined as a product of k elementary reflectors. +*> +*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; +*> +*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. +*> +*> If STOREV = 'C', the vector which defines the elementary reflector +*> H(i) is stored in the i-th column of the array V, and +*> +*> H = I - V * T * V**H +*> +*> If STOREV = 'R', the vector which defines the elementary reflector +*> H(i) is stored in the i-th row of the array V, and +*> +*> H = I - V**H * T * V +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Specifies the order in which the elementary reflectors are +*> multiplied to form the block reflector: +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Specifies how the vectors which define the elementary +*> reflectors are stored (see also Further Details): +*> = 'C': columnwise +*> = 'R': rowwise +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the block reflector H. N >= 0. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the triangular factor T (= the number of +*> elementary reflectors). K >= 1. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,N) if STOREV = 'R' +*> The matrix V. See further details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is COMPLEX array, dimension (K) +*> TAU(i) must contain the scalar factor of the elementary +*> reflector H(i). +*> \endverbatim +*> +*> \param[out] T +*> \verbatim +*> T is COMPLEX array, dimension (LDT,K) +*> The k by k triangular factor T of the block reflector. +*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is +*> lower triangular. The rest of the array is not used. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup complexOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* -- LAPACK auxiliary routine (version 3.4.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. + COMPLEX T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE, ZERO + PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), + $ ZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J, PREVLASTV, LASTV +* .. +* .. External Subroutines .. + EXTERNAL CGEMV, CLACGV, CTRMV +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( LSAME( DIRECT, 'F' ) ) THEN + PREVLASTV = N + DO I = 1, K + PREVLASTV = MAX( PREVLASTV, I ) + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = 1, I + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) +* + CALL CGEMV( 'Conjugate transpose', J-I, I-1, + $ -TAU( I ), V( I+1, 1 ), LDV, + $ V( I+1, I ), 1, + $ ONE, T( 1, I ), 1 ) + ELSE +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( J , I ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H +* + CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), + $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, + $ ONE, T( 1, I ), LDT ) + END IF +* +* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) +* + CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, + $ LDT, T( 1, I ), 1 ) + T( I, I ) = TAU( I ) + IF( I.GT.1 ) THEN + PREVLASTV = MAX( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + END DO + ELSE + PREVLASTV = 1 + DO I = K, 1, -1 + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = I, K + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( I.LT.K ) THEN + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) +* + CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I, + $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), + $ 1, ONE, T( I+1, I ), 1 ) + ELSE +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( J, N-K+I ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H +* + CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), + $ V( I+1, J ), LDV, V( I, J ), LDV, + $ ONE, T( I+1, I ), LDT ) + END IF +* +* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) +* + CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, + $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) + IF( I.GT.1 ) THEN + PREVLASTV = MIN( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + T( I, I ) = TAU( I ) + END IF + END DO + END IF + RETURN +* +* End of CLARFT +* + END diff --git a/lapack/dlarf.f b/lapack/dlarf.f new file mode 100644 index 000000000..2a82ff439 --- /dev/null +++ b/lapack/dlarf.f @@ -0,0 +1,227 @@ +*> \brief \b DLARF +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARF + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* DOUBLE PRECISION TAU +* .. +* .. Array Arguments .. +* DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLARF applies a real elementary reflector H to a real m by n matrix +*> C, from either the left or the right. H is represented in the form +*> +*> H = I - tau * v * v**T +*> +*> where tau is a real scalar and v is a real vector. +*> +*> If tau = 0, then H is taken to be the unit matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form H * C +*> = 'R': form C * H +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of H. V is not used if +*> TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION +*> The value tau in the representation of H. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is DOUBLE PRECISION array, dimension (LDC,N) +*> On entry, the m by n matrix C. +*> On exit, C is overwritten by the matrix H * C if SIDE = 'L', +*> or C * H if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension +*> (N) if SIDE = 'L' +*> or (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleOTHERauxiliary +* +* ===================================================================== + SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + DOUBLE PRECISION TAU +* .. +* .. Array Arguments .. + DOUBLE PRECISION C( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL APPLYLEFT + INTEGER I, LASTV, LASTC +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, DGER +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILADLR, ILADLC + EXTERNAL LSAME, ILADLR, ILADLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 0 + LASTC = 0 + IF( TAU.NE.ZERO ) THEN +! Set up variables for scanning V. LASTV begins pointing to the end +! of V. + IF( APPLYLEFT ) THEN + LASTV = M + ELSE + LASTV = N + END IF + IF( INCV.GT.0 ) THEN + I = 1 + (LASTV-1) * INCV + ELSE + I = 1 + END IF +! Look for the last non-zero row in V. + DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) + LASTV = LASTV - 1 + I = I - INCV + END DO + IF( APPLYLEFT ) THEN +! Scan for the last non-zero column in C(1:lastv,:). + LASTC = ILADLC(LASTV, N, C, LDC) + ELSE +! Scan for the last non-zero row in C(:,1:lastv). + LASTC = ILADLR(M, LASTV, C, LDC) + END IF + END IF +! Note that lastc.eq.0 renders the BLAS operations null; no special +! case is needed at this level. + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1) +* + CALL DGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV, + $ ZERO, WORK, 1 ) +* +* C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T +* + CALL DGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC ) + END IF + ELSE +* +* Form C * H +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + CALL DGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC, + $ V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T +* + CALL DGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC ) + END IF + END IF + RETURN +* +* End of DLARF +* + END diff --git a/lapack/dlarfb.f b/lapack/dlarfb.f new file mode 100644 index 000000000..206d3b268 --- /dev/null +++ b/lapack/dlarfb.f @@ -0,0 +1,762 @@ +*> \brief \b DLARFB +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARFB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, +* T, LDT, C, LDC, WORK, LDWORK ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, SIDE, STOREV, TRANS +* INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), +* $ WORK( LDWORK, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLARFB applies a real block reflector H or its transpose H**T to a +*> real m by n matrix C, from either the left or the right. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': apply H or H**T from the Left +*> = 'R': apply H or H**T from the Right +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> = 'N': apply H (No transpose) +*> = 'T': apply H**T (Transpose) +*> \endverbatim +*> +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Indicates how H is formed from a product of elementary +*> reflectors +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Indicates how the vectors which define the elementary +*> reflectors are stored: +*> = 'C': Columnwise +*> = 'R': Rowwise +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the matrix T (= the number of elementary +*> reflectors whose product defines the block reflector). +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,M) if STOREV = 'R' and SIDE = 'L' +*> (LDV,N) if STOREV = 'R' and SIDE = 'R' +*> The matrix V. See Further Details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +*> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +*> if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] T +*> \verbatim +*> T is DOUBLE PRECISION array, dimension (LDT,K) +*> The triangular k by k matrix T in the representation of the +*> block reflector. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is DOUBLE PRECISION array, dimension (LDC,N) +*> On entry, the m by n matrix C. +*> On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is DOUBLE PRECISION array, dimension (LDWORK,K) +*> \endverbatim +*> +*> \param[in] LDWORK +*> \verbatim +*> LDWORK is INTEGER +*> The leading dimension of the array WORK. +*> If SIDE = 'L', LDWORK >= max(1,N); +*> if SIDE = 'R', LDWORK >= max(1,M). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored; the corresponding +*> array elements are modified but restored on exit. The rest of the +*> array is not used. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, + $ T, LDT, C, LDC, WORK, LDWORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION C( LDC, * ), T( LDT, * ), V( LDV, * ), + $ WORK( LDWORK, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + CHARACTER TRANST + INTEGER I, J, LASTV, LASTC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILADLR, ILADLC + EXTERNAL LSAME, ILADLR, ILADLC +* .. +* .. External Subroutines .. + EXTERNAL DCOPY, DGEMM, DTRMM +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 ) + $ RETURN +* + IF( LSAME( TRANS, 'N' ) ) THEN + TRANST = 'T' + ELSE + TRANST = 'N' + END IF +* + IF( LSAME( STOREV, 'C' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 ) (first K rows) +* ( V2 ) +* where V1 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) + LASTC = ILADLC( LASTV, N, C, LDC ) +* +* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) +* +* W := C1**T +* + DO 10 J = 1, K + CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + 10 CONTINUE +* +* W := W * V1 +* + CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**T *V2 +* + CALL DGEMM( 'Transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2 * W**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, + $ C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1**T +* + CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**T +* + DO 30 J = 1, K + DO 20 I = 1, LASTC + C( J, I ) = C( J, I ) - WORK( I, J ) + 20 CONTINUE + 30 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) + LASTC = ILADLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C1 +* + DO 40 J = 1, K + CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 40 CONTINUE +* +* W := W * V1 +* + CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2 +* + CALL DGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, + $ C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1**T +* + CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 60 J = 1, K + DO 50 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 50 CONTINUE + 60 CONTINUE + END IF +* + ELSE +* +* Let V = ( V1 ) +* ( V2 ) (last K rows) +* where V2 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILADLR( M, K, V, LDV ) ) + LASTC = ILADLC( LASTV, N, C, LDC ) +* +* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) +* +* W := C2**T +* + DO 70 J = 1, K + CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + 70 CONTINUE +* +* W := W * V2 +* + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**T*V1 +* + CALL DGEMM( 'Transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1 * W**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**T +* + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**T +* + DO 90 J = 1, K + DO 80 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) + 80 CONTINUE + 90 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILADLR( N, K, V, LDV ) ) + LASTC = ILADLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C2 +* + DO 100 J = 1, K + CALL DCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) + 100 CONTINUE +* +* W := W * V2 +* + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1 +* + CALL DGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**T +* + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W +* + DO 120 J = 1, K + DO 110 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) + 110 CONTINUE + 120 CONTINUE + END IF + END IF +* + ELSE IF( LSAME( STOREV, 'R' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 V2 ) (V1: first K columns) +* where V1 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) + LASTC = ILADLC( LASTV, N, C, LDC ) +* +* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) +* +* W := C1**T +* + DO 130 J = 1, K + CALL DCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + 130 CONTINUE +* +* W := W * V1**T +* + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**T*V2**T +* + CALL DGEMM( 'Transpose', 'Transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL DTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**T * W**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2**T * W**T +* + CALL DGEMM( 'Transpose', 'Transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, + $ ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**T +* + DO 150 J = 1, K + DO 140 I = 1, LASTC + C( J, I ) = C( J, I ) - WORK( I, J ) + 140 CONTINUE + 150 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) + LASTC = ILADLR( M, LASTV, C, LDC ) +* +* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) +* +* W := C1 +* + DO 160 J = 1, K + CALL DCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 160 CONTINUE +* +* W := W * V1**T +* + CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL DTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2 +* + CALL DGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL DTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 180 J = 1, K + DO 170 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 170 CONTINUE + 180 CONTINUE +* + END IF +* + ELSE +* +* Let V = ( V1 V2 ) (V2: last K columns) +* where V2 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILADLC( K, M, V, LDV ) ) + LASTC = ILADLC( LASTV, N, C, LDC ) +* +* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) +* +* W := C2**T +* + DO 190 J = 1, K + CALL DCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + 190 CONTINUE +* +* W := W * V2**T +* + CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**T * V1**T +* + CALL DGEMM( 'Transpose', 'Transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL DTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**T * W**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1**T * W**T +* + CALL DGEMM( 'Transpose', 'Transpose', + $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**T +* + DO 210 J = 1, K + DO 200 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) + 200 CONTINUE + 210 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILADLC( K, N, V, LDV ) ) + LASTC = ILADLR( M, LASTV, C, LDC ) +* +* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) +* +* W := C2 +* + DO 220 J = 1, K + CALL DCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 220 CONTINUE +* +* W := W * V2**T +* + CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1**T +* + CALL DGEMM( 'No transpose', 'Transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL DTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1 +* + CALL DGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL DTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 240 J = 1, K + DO 230 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) + 230 CONTINUE + 240 CONTINUE +* + END IF +* + END IF + END IF +* + RETURN +* +* End of DLARFB +* + END diff --git a/lapack/dlarfg.f b/lapack/dlarfg.f new file mode 100644 index 000000000..458ad2e05 --- /dev/null +++ b/lapack/dlarfg.f @@ -0,0 +1,196 @@ +*> \brief \b DLARFG +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARFG + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* DOUBLE PRECISION ALPHA, TAU +* .. +* .. Array Arguments .. +* DOUBLE PRECISION X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLARFG generates a real elementary reflector H of order n, such +*> that +*> +*> H * ( alpha ) = ( beta ), H**T * H = I. +*> ( x ) ( 0 ) +*> +*> where alpha and beta are scalars, and x is an (n-1)-element real +*> vector. H is represented in the form +*> +*> H = I - tau * ( 1 ) * ( 1 v**T ) , +*> ( v ) +*> +*> where tau is a real scalar and v is a real (n-1)-element +*> vector. +*> +*> If the elements of x are all zero, then tau = 0 and H is taken to be +*> the unit matrix. +*> +*> Otherwise 1 <= tau <= 2. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the elementary reflector. +*> \endverbatim +*> +*> \param[in,out] ALPHA +*> \verbatim +*> ALPHA is DOUBLE PRECISION +*> On entry, the value alpha. +*> On exit, it is overwritten with the value beta. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is DOUBLE PRECISION array, dimension +*> (1+(N-2)*abs(INCX)) +*> On entry, the vector x. +*> On exit, it is overwritten with the vector v. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between elements of X. INCX > 0. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION +*> The value tau. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup doubleOTHERauxiliary +* +* ===================================================================== + SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + INTEGER INCX, N + DOUBLE PRECISION ALPHA, TAU +* .. +* .. Array Arguments .. + DOUBLE PRECISION X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2 + EXTERNAL DLAMCH, DLAPY2, DNRM2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, SIGN +* .. +* .. External Subroutines .. + EXTERNAL DSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.1 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = DNRM2( N-1, X, INCX ) +* + IF( XNORM.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + RSAFMN = ONE / SAFMIN + 10 CONTINUE + KNT = KNT + 1 + CALL DSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHA = ALPHA*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = DNRM2( N-1, X, INCX ) + BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA ) + END IF + TAU = ( BETA-ALPHA ) / BETA + CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of DLARFG +* + END diff --git a/lapack/dlarft.f b/lapack/dlarft.f new file mode 100644 index 000000000..4b7550403 --- /dev/null +++ b/lapack/dlarft.f @@ -0,0 +1,326 @@ +*> \brief \b DLARFT +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARFT + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, STOREV +* INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. +* DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> DLARFT forms the triangular factor T of a real block reflector H +*> of order n, which is defined as a product of k elementary reflectors. +*> +*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; +*> +*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. +*> +*> If STOREV = 'C', the vector which defines the elementary reflector +*> H(i) is stored in the i-th column of the array V, and +*> +*> H = I - V * T * V**T +*> +*> If STOREV = 'R', the vector which defines the elementary reflector +*> H(i) is stored in the i-th row of the array V, and +*> +*> H = I - V**T * T * V +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Specifies the order in which the elementary reflectors are +*> multiplied to form the block reflector: +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Specifies how the vectors which define the elementary +*> reflectors are stored (see also Further Details): +*> = 'C': columnwise +*> = 'R': rowwise +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the block reflector H. N >= 0. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the triangular factor T (= the number of +*> elementary reflectors). K >= 1. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is DOUBLE PRECISION array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,N) if STOREV = 'R' +*> The matrix V. See further details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is DOUBLE PRECISION array, dimension (K) +*> TAU(i) must contain the scalar factor of the elementary +*> reflector H(i). +*> \endverbatim +*> +*> \param[out] T +*> \verbatim +*> T is DOUBLE PRECISION array, dimension (LDT,K) +*> The k by k triangular factor T of the block reflector. +*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is +*> lower triangular. The rest of the array is not used. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup doubleOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* -- LAPACK auxiliary routine (version 3.4.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, PREVLASTV, LASTV +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, DTRMV +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( LSAME( DIRECT, 'F' ) ) THEN + PREVLASTV = N + DO I = 1, K + PREVLASTV = MAX( I, PREVLASTV ) + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = 1, I + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( I , J ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) +* + CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ), + $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, + $ T( 1, I ), 1 ) + ELSE +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( J , I ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T +* + CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ), + $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE, + $ T( 1, I ), 1 ) + END IF +* +* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) +* + CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, + $ LDT, T( 1, I ), 1 ) + T( I, I ) = TAU( I ) + IF( I.GT.1 ) THEN + PREVLASTV = MAX( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + END DO + ELSE + PREVLASTV = 1 + DO I = K, 1, -1 + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = I, K + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( I.LT.K ) THEN + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( N-K+I , J ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) +* + CALL DGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ), + $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, + $ T( I+1, I ), 1 ) + ELSE +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( J, N-K+I ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T +* + CALL DGEMV( 'No transpose', K-I, N-K+I-J, + $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, + $ ONE, T( I+1, I ), 1 ) + END IF +* +* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) +* + CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, + $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) + IF( I.GT.1 ) THEN + PREVLASTV = MIN( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + T( I, I ) = TAU( I ) + END IF + END DO + END IF + RETURN +* +* End of DLARFT +* + END diff --git a/lapack/slarf.f b/lapack/slarf.f new file mode 100644 index 000000000..8a8ff308e --- /dev/null +++ b/lapack/slarf.f @@ -0,0 +1,227 @@ +*> \brief \b SLARF +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SLARF + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* REAL TAU +* .. +* .. Array Arguments .. +* REAL C( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLARF applies a real elementary reflector H to a real m by n matrix +*> C, from either the left or the right. H is represented in the form +*> +*> H = I - tau * v * v**T +*> +*> where tau is a real scalar and v is a real vector. +*> +*> If tau = 0, then H is taken to be the unit matrix. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form H * C +*> = 'R': form C * H +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is REAL array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of H. V is not used if +*> TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is REAL +*> The value tau in the representation of H. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is REAL array, dimension (LDC,N) +*> On entry, the m by n matrix C. +*> On exit, C is overwritten by the matrix H * C if SIDE = 'L', +*> or C * H if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension +*> (N) if SIDE = 'L' +*> or (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup realOTHERauxiliary +* +* ===================================================================== + SUBROUTINE SLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + REAL TAU +* .. +* .. Array Arguments .. + REAL C( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL APPLYLEFT + INTEGER I, LASTV, LASTC +* .. +* .. External Subroutines .. + EXTERNAL SGEMV, SGER +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILASLR, ILASLC + EXTERNAL LSAME, ILASLR, ILASLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 0 + LASTC = 0 + IF( TAU.NE.ZERO ) THEN +! Set up variables for scanning V. LASTV begins pointing to the end +! of V. + IF( APPLYLEFT ) THEN + LASTV = M + ELSE + LASTV = N + END IF + IF( INCV.GT.0 ) THEN + I = 1 + (LASTV-1) * INCV + ELSE + I = 1 + END IF +! Look for the last non-zero row in V. + DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) + LASTV = LASTV - 1 + I = I - INCV + END DO + IF( APPLYLEFT ) THEN +! Scan for the last non-zero column in C(1:lastv,:). + LASTC = ILASLC(LASTV, N, C, LDC) + ELSE +! Scan for the last non-zero row in C(:,1:lastv). + LASTC = ILASLR(M, LASTV, C, LDC) + END IF + END IF +! Note that lastc.eq.0 renders the BLAS operations null; no special +! case is needed at this level. + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1) +* + CALL SGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV, + $ ZERO, WORK, 1 ) +* +* C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T +* + CALL SGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC ) + END IF + ELSE +* +* Form C * H +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + CALL SGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC, + $ V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T +* + CALL SGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC ) + END IF + END IF + RETURN +* +* End of SLARF +* + END diff --git a/lapack/slarfb.f b/lapack/slarfb.f new file mode 100644 index 000000000..eb95990b3 --- /dev/null +++ b/lapack/slarfb.f @@ -0,0 +1,763 @@ +*> \brief \b SLARFB +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SLARFB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, +* T, LDT, C, LDC, WORK, LDWORK ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, SIDE, STOREV, TRANS +* INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. +* REAL C( LDC, * ), T( LDT, * ), V( LDV, * ), +* $ WORK( LDWORK, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLARFB applies a real block reflector H or its transpose H**T to a +*> real m by n matrix C, from either the left or the right. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': apply H or H**T from the Left +*> = 'R': apply H or H**T from the Right +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> = 'N': apply H (No transpose) +*> = 'T': apply H**T (Transpose) +*> \endverbatim +*> +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Indicates how H is formed from a product of elementary +*> reflectors +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Indicates how the vectors which define the elementary +*> reflectors are stored: +*> = 'C': Columnwise +*> = 'R': Rowwise +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the matrix T (= the number of elementary +*> reflectors whose product defines the block reflector). +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is REAL array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,M) if STOREV = 'R' and SIDE = 'L' +*> (LDV,N) if STOREV = 'R' and SIDE = 'R' +*> The matrix V. See Further Details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +*> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +*> if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] T +*> \verbatim +*> T is REAL array, dimension (LDT,K) +*> The triangular k by k matrix T in the representation of the +*> block reflector. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is REAL array, dimension (LDC,N) +*> On entry, the m by n matrix C. +*> On exit, C is overwritten by H*C or H**T*C or C*H or C*H**T. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is REAL array, dimension (LDWORK,K) +*> \endverbatim +*> +*> \param[in] LDWORK +*> \verbatim +*> LDWORK is INTEGER +*> The leading dimension of the array WORK. +*> If SIDE = 'L', LDWORK >= max(1,N); +*> if SIDE = 'R', LDWORK >= max(1,M). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup realOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored; the corresponding +*> array elements are modified but restored on exit. The rest of the +*> array is not used. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE SLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, + $ T, LDT, C, LDC, WORK, LDWORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + REAL C( LDC, * ), T( LDT, * ), V( LDV, * ), + $ WORK( LDWORK, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + CHARACTER TRANST + INTEGER I, J, LASTV, LASTC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILASLR, ILASLC + EXTERNAL LSAME, ILASLR, ILASLC +* .. +* .. External Subroutines .. + EXTERNAL SCOPY, SGEMM, STRMM +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 ) + $ RETURN +* + IF( LSAME( TRANS, 'N' ) ) THEN + TRANST = 'T' + ELSE + TRANST = 'N' + END IF +* + IF( LSAME( STOREV, 'C' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 ) (first K rows) +* ( V2 ) +* where V1 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILASLR( M, K, V, LDV ) ) + LASTC = ILASLC( LASTV, N, C, LDC ) +* +* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) +* +* W := C1**T +* + DO 10 J = 1, K + CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + 10 CONTINUE +* +* W := W * V1 +* + CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**T *V2 +* + CALL SGEMM( 'Transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2 * W**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, ONE, + $ C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1**T +* + CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**T +* + DO 30 J = 1, K + DO 20 I = 1, LASTC + C( J, I ) = C( J, I ) - WORK( I, J ) + 20 CONTINUE + 30 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILASLR( N, K, V, LDV ) ) + LASTC = ILASLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C1 +* + DO 40 J = 1, K + CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 40 CONTINUE +* +* W := W * V1 +* + CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2 +* + CALL SGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, ONE, + $ C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1**T +* + CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 60 J = 1, K + DO 50 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 50 CONTINUE + 60 CONTINUE + END IF +* + ELSE +* +* Let V = ( V1 ) +* ( V2 ) (last K rows) +* where V2 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILASLR( M, K, V, LDV ) ) + LASTC = ILASLC( LASTV, N, C, LDC ) +* +* W := C**T * V = (C1**T * V1 + C2**T * V2) (stored in WORK) +* +* W := C2**T +* + DO 70 J = 1, K + CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + 70 CONTINUE +* +* W := W * V2 +* + CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**T*V1 +* + CALL SGEMM( 'Transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1 * W**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**T +* + CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**T +* + DO 90 J = 1, K + DO 80 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) + 80 CONTINUE + 90 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILASLR( N, K, V, LDV ) ) + LASTC = ILASLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C2 +* + DO 100 J = 1, K + CALL SCOPY( LASTC, C( 1, N-K+J ), 1, WORK( 1, J ), 1 ) + 100 CONTINUE +* +* W := W * V2 +* + CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1 +* + CALL SGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**T +* + CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W +* + DO 120 J = 1, K + DO 110 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) - WORK(I, J) + 110 CONTINUE + 120 CONTINUE + END IF + END IF +* + ELSE IF( LSAME( STOREV, 'R' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 V2 ) (V1: first K columns) +* where V1 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILASLC( K, M, V, LDV ) ) + LASTC = ILASLC( LASTV, N, C, LDC ) +* +* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) +* +* W := C1**T +* + DO 130 J = 1, K + CALL SCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + 130 CONTINUE +* +* W := W * V1**T +* + CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**T*V2**T +* + CALL SGEMM( 'Transpose', 'Transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL STRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**T * W**T +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2**T * W**T +* + CALL SGEMM( 'Transpose', 'Transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, + $ ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**T +* + DO 150 J = 1, K + DO 140 I = 1, LASTC + C( J, I ) = C( J, I ) - WORK( I, J ) + 140 CONTINUE + 150 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILASLC( K, N, V, LDV ) ) + LASTC = ILASLR( M, LASTV, C, LDC ) +* +* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) +* +* W := C1 +* + DO 160 J = 1, K + CALL SCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 160 CONTINUE +* +* W := W * V1**T +* + CALL STRMM( 'Right', 'Upper', 'Transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL STRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2 +* + CALL SGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL STRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 180 J = 1, K + DO 170 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 170 CONTINUE + 180 CONTINUE +* + END IF +* + ELSE +* +* Let V = ( V1 V2 ) (V2: last K columns) +* where V2 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**T * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILASLC( K, M, V, LDV ) ) + LASTC = ILASLC( LASTV, N, C, LDC ) +* +* W := C**T * V**T = (C1**T * V1**T + C2**T * V2**T) (stored in WORK) +* +* W := C2**T +* + DO 190 J = 1, K + CALL SCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + 190 CONTINUE +* +* W := W * V2**T +* + CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**T * V1**T +* + CALL SGEMM( 'Transpose', 'Transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**T or W * T +* + CALL STRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**T * W**T +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1**T * W**T +* + CALL SGEMM( 'Transpose', 'Transpose', + $ LASTV-K, LASTC, K, -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**T +* + DO 210 J = 1, K + DO 200 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - WORK(I, J) + 200 CONTINUE + 210 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**T where C = ( C1 C2 ) +* + LASTV = MAX( K, ILASLC( K, N, V, LDV ) ) + LASTC = ILASLR( M, LASTV, C, LDC ) +* +* W := C * V**T = (C1*V1**T + C2*V2**T) (stored in WORK) +* +* W := C2 +* + DO 220 J = 1, K + CALL SCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 220 CONTINUE +* +* W := W * V2**T +* + CALL STRMM( 'Right', 'Lower', 'Transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1**T +* + CALL SGEMM( 'No transpose', 'Transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**T +* + CALL STRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1 +* + CALL SGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL STRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 240 J = 1, K + DO 230 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) + $ - WORK( I, J ) + 230 CONTINUE + 240 CONTINUE +* + END IF +* + END IF + END IF +* + RETURN +* +* End of SLARFB +* + END diff --git a/lapack/slarfg.f b/lapack/slarfg.f new file mode 100644 index 000000000..4f10ffcaf --- /dev/null +++ b/lapack/slarfg.f @@ -0,0 +1,196 @@ +*> \brief \b SLARFG +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SLARFG + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* REAL ALPHA, TAU +* .. +* .. Array Arguments .. +* REAL X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLARFG generates a real elementary reflector H of order n, such +*> that +*> +*> H * ( alpha ) = ( beta ), H**T * H = I. +*> ( x ) ( 0 ) +*> +*> where alpha and beta are scalars, and x is an (n-1)-element real +*> vector. H is represented in the form +*> +*> H = I - tau * ( 1 ) * ( 1 v**T ) , +*> ( v ) +*> +*> where tau is a real scalar and v is a real (n-1)-element +*> vector. +*> +*> If the elements of x are all zero, then tau = 0 and H is taken to be +*> the unit matrix. +*> +*> Otherwise 1 <= tau <= 2. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the elementary reflector. +*> \endverbatim +*> +*> \param[in,out] ALPHA +*> \verbatim +*> ALPHA is REAL +*> On entry, the value alpha. +*> On exit, it is overwritten with the value beta. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is REAL array, dimension +*> (1+(N-2)*abs(INCX)) +*> On entry, the vector x. +*> On exit, it is overwritten with the vector v. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between elements of X. INCX > 0. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is REAL +*> The value tau. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup realOTHERauxiliary +* +* ===================================================================== + SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + INTEGER INCX, N + REAL ALPHA, TAU +* .. +* .. Array Arguments .. + REAL X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + REAL BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + REAL SLAMCH, SLAPY2, SNRM2 + EXTERNAL SLAMCH, SLAPY2, SNRM2 +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, SIGN +* .. +* .. External Subroutines .. + EXTERNAL SSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.1 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = SNRM2( N-1, X, INCX ) +* + IF( XNORM.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' ) + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + RSAFMN = ONE / SAFMIN + 10 CONTINUE + KNT = KNT + 1 + CALL SSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHA = ALPHA*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = SNRM2( N-1, X, INCX ) + BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA ) + END IF + TAU = ( BETA-ALPHA ) / BETA + CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of SLARFG +* + END diff --git a/lapack/slarft.f b/lapack/slarft.f new file mode 100644 index 000000000..30b0668e4 --- /dev/null +++ b/lapack/slarft.f @@ -0,0 +1,326 @@ +*> \brief \b SLARFT +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download SLARFT + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, STOREV +* INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. +* REAL T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> SLARFT forms the triangular factor T of a real block reflector H +*> of order n, which is defined as a product of k elementary reflectors. +*> +*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; +*> +*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. +*> +*> If STOREV = 'C', the vector which defines the elementary reflector +*> H(i) is stored in the i-th column of the array V, and +*> +*> H = I - V * T * V**T +*> +*> If STOREV = 'R', the vector which defines the elementary reflector +*> H(i) is stored in the i-th row of the array V, and +*> +*> H = I - V**T * T * V +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Specifies the order in which the elementary reflectors are +*> multiplied to form the block reflector: +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Specifies how the vectors which define the elementary +*> reflectors are stored (see also Further Details): +*> = 'C': columnwise +*> = 'R': rowwise +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the block reflector H. N >= 0. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the triangular factor T (= the number of +*> elementary reflectors). K >= 1. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is REAL array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,N) if STOREV = 'R' +*> The matrix V. See further details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is REAL array, dimension (K) +*> TAU(i) must contain the scalar factor of the elementary +*> reflector H(i). +*> \endverbatim +*> +*> \param[out] T +*> \verbatim +*> T is REAL array, dimension (LDT,K) +*> The k by k triangular factor T of the block reflector. +*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is +*> lower triangular. The rest of the array is not used. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup realOTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* -- LAPACK auxiliary routine (version 3.4.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. + REAL T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE, ZERO + PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) +* .. +* .. Local Scalars .. + INTEGER I, J, PREVLASTV, LASTV +* .. +* .. External Subroutines .. + EXTERNAL SGEMV, STRMV +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( LSAME( DIRECT, 'F' ) ) THEN + PREVLASTV = N + DO I = 1, K + PREVLASTV = MAX( I, PREVLASTV ) + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = 1, I + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( I , J ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) +* + CALL SGEMV( 'Transpose', J-I, I-1, -TAU( I ), + $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, + $ T( 1, I ), 1 ) + ELSE +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( J , I ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T +* + CALL SGEMV( 'No transpose', I-1, J-I, -TAU( I ), + $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, + $ ONE, T( 1, I ), 1 ) + END IF +* +* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) +* + CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, + $ LDT, T( 1, I ), 1 ) + T( I, I ) = TAU( I ) + IF( I.GT.1 ) THEN + PREVLASTV = MAX( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + END DO + ELSE + PREVLASTV = 1 + DO I = K, 1, -1 + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = I, K + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( I.LT.K ) THEN + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( N-K+I , J ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) +* + CALL SGEMV( 'Transpose', N-K+I-J, K-I, -TAU( I ), + $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, + $ T( I+1, I ), 1 ) + ELSE +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( J, N-K+I ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T +* + CALL SGEMV( 'No transpose', K-I, N-K+I-J, + $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, + $ ONE, T( I+1, I ), 1 ) + END IF +* +* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) +* + CALL STRMV( 'Lower', 'No transpose', 'Non-unit', K-I, + $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) + IF( I.GT.1 ) THEN + PREVLASTV = MIN( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + T( I, I ) = TAU( I ) + END IF + END DO + END IF + RETURN +* +* End of SLARFT +* + END diff --git a/lapack/zlarf.f b/lapack/zlarf.f new file mode 100644 index 000000000..53f314d64 --- /dev/null +++ b/lapack/zlarf.f @@ -0,0 +1,232 @@ +*> \brief \b ZLARF +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARF + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* .. Scalar Arguments .. +* CHARACTER SIDE +* INTEGER INCV, LDC, M, N +* COMPLEX*16 TAU +* .. +* .. Array Arguments .. +* COMPLEX*16 C( LDC, * ), V( * ), WORK( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLARF applies a complex elementary reflector H to a complex M-by-N +*> matrix C, from either the left or the right. H is represented in the +*> form +*> +*> H = I - tau * v * v**H +*> +*> where tau is a complex scalar and v is a complex vector. +*> +*> If tau = 0, then H is taken to be the unit matrix. +*> +*> To apply H**H, supply conjg(tau) instead +*> tau. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': form H * C +*> = 'R': form C * H +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX*16 array, dimension +*> (1 + (M-1)*abs(INCV)) if SIDE = 'L' +*> or (1 + (N-1)*abs(INCV)) if SIDE = 'R' +*> The vector v in the representation of H. V is not used if +*> TAU = 0. +*> \endverbatim +*> +*> \param[in] INCV +*> \verbatim +*> INCV is INTEGER +*> The increment between elements of v. INCV <> 0. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is COMPLEX*16 +*> The value tau in the representation of H. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is COMPLEX*16 array, dimension (LDC,N) +*> On entry, the M-by-N matrix C. +*> On exit, C is overwritten by the matrix H * C if SIDE = 'L', +*> or C * H if SIDE = 'R'. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension +*> (N) if SIDE = 'L' +*> or (M) if SIDE = 'R' +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERauxiliary +* +* ===================================================================== + SUBROUTINE ZLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER SIDE + INTEGER INCV, LDC, M, N + COMPLEX*16 TAU +* .. +* .. Array Arguments .. + COMPLEX*16 C( LDC, * ), V( * ), WORK( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + LOGICAL APPLYLEFT + INTEGER I, LASTV, LASTC +* .. +* .. External Subroutines .. + EXTERNAL ZGEMV, ZGERC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAZLR, ILAZLC + EXTERNAL LSAME, ILAZLR, ILAZLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 0 + LASTC = 0 + IF( TAU.NE.ZERO ) THEN +* Set up variables for scanning V. LASTV begins pointing to the end +* of V. + IF( APPLYLEFT ) THEN + LASTV = M + ELSE + LASTV = N + END IF + IF( INCV.GT.0 ) THEN + I = 1 + (LASTV-1) * INCV + ELSE + I = 1 + END IF +* Look for the last non-zero row in V. + DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO ) + LASTV = LASTV - 1 + I = I - INCV + END DO + IF( APPLYLEFT ) THEN +* Scan for the last non-zero column in C(1:lastv,:). + LASTC = ILAZLC(LASTV, N, C, LDC) + ELSE +* Scan for the last non-zero row in C(:,1:lastv). + LASTC = ILAZLR(M, LASTV, C, LDC) + END IF + END IF +* Note that lastc.eq.0 renders the BLAS operations null; no special +* case is needed at this level. + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**H * v(1:lastv,1) +* + CALL ZGEMV( 'Conjugate transpose', LASTV, LASTC, ONE, + $ C, LDC, V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**H +* + CALL ZGERC( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC ) + END IF + ELSE +* +* Form C * H +* + IF( LASTV.GT.0 ) THEN +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + CALL ZGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC, + $ V, INCV, ZERO, WORK, 1 ) +* +* C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**H +* + CALL ZGERC( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC ) + END IF + END IF + RETURN +* +* End of ZLARF +* + END diff --git a/lapack/zlarfb.f b/lapack/zlarfb.f new file mode 100644 index 000000000..30fc4b940 --- /dev/null +++ b/lapack/zlarfb.f @@ -0,0 +1,774 @@ +*> \brief \b ZLARFB +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARFB + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, +* T, LDT, C, LDC, WORK, LDWORK ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, SIDE, STOREV, TRANS +* INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. +* COMPLEX*16 C( LDC, * ), T( LDT, * ), V( LDV, * ), +* $ WORK( LDWORK, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLARFB applies a complex block reflector H or its transpose H**H to a +*> complex M-by-N matrix C, from either the left or the right. +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] SIDE +*> \verbatim +*> SIDE is CHARACTER*1 +*> = 'L': apply H or H**H from the Left +*> = 'R': apply H or H**H from the Right +*> \endverbatim +*> +*> \param[in] TRANS +*> \verbatim +*> TRANS is CHARACTER*1 +*> = 'N': apply H (No transpose) +*> = 'C': apply H**H (Conjugate transpose) +*> \endverbatim +*> +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Indicates how H is formed from a product of elementary +*> reflectors +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Indicates how the vectors which define the elementary +*> reflectors are stored: +*> = 'C': Columnwise +*> = 'R': Rowwise +*> \endverbatim +*> +*> \param[in] M +*> \verbatim +*> M is INTEGER +*> The number of rows of the matrix C. +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The number of columns of the matrix C. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the matrix T (= the number of elementary +*> reflectors whose product defines the block reflector). +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX*16 array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,M) if STOREV = 'R' and SIDE = 'L' +*> (LDV,N) if STOREV = 'R' and SIDE = 'R' +*> See Further Details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M); +*> if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N); +*> if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] T +*> \verbatim +*> T is COMPLEX*16 array, dimension (LDT,K) +*> The triangular K-by-K matrix T in the representation of the +*> block reflector. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +*> +*> \param[in,out] C +*> \verbatim +*> C is COMPLEX*16 array, dimension (LDC,N) +*> On entry, the M-by-N matrix C. +*> On exit, C is overwritten by H*C or H**H*C or C*H or C*H**H. +*> \endverbatim +*> +*> \param[in] LDC +*> \verbatim +*> LDC is INTEGER +*> The leading dimension of the array C. LDC >= max(1,M). +*> \endverbatim +*> +*> \param[out] WORK +*> \verbatim +*> WORK is COMPLEX*16 array, dimension (LDWORK,K) +*> \endverbatim +*> +*> \param[in] LDWORK +*> \verbatim +*> LDWORK is INTEGER +*> The leading dimension of the array WORK. +*> If SIDE = 'L', LDWORK >= max(1,N); +*> if SIDE = 'R', LDWORK >= max(1,M). +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored; the corresponding +*> array elements are modified but restored on exit. The rest of the +*> array is not used. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE ZLARFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, V, LDV, + $ T, LDT, C, LDC, WORK, LDWORK ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, SIDE, STOREV, TRANS + INTEGER K, LDC, LDT, LDV, LDWORK, M, N +* .. +* .. Array Arguments .. + COMPLEX*16 C( LDC, * ), T( LDT, * ), V( LDV, * ), + $ WORK( LDWORK, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + CHARACTER TRANST + INTEGER I, J, LASTV, LASTC +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAZLR, ILAZLC + EXTERNAL LSAME, ILAZLR, ILAZLC +* .. +* .. External Subroutines .. + EXTERNAL ZCOPY, ZGEMM, ZLACGV, ZTRMM +* .. +* .. Intrinsic Functions .. + INTRINSIC DCONJG +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( M.LE.0 .OR. N.LE.0 ) + $ RETURN +* + IF( LSAME( TRANS, 'N' ) ) THEN + TRANST = 'C' + ELSE + TRANST = 'N' + END IF +* + IF( LSAME( STOREV, 'C' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 ) (first K rows) +* ( V2 ) +* where V1 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILAZLR( M, K, V, LDV ) ) + LASTC = ILAZLC( LASTV, N, C, LDC ) +* +* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) +* +* W := C1**H +* + DO 10 J = 1, K + CALL ZCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + CALL ZLACGV( LASTC, WORK( 1, J ), 1 ) + 10 CONTINUE +* +* W := W * V1 +* + CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**H *V2 +* + CALL ZGEMM( 'Conjugate transpose', 'No transpose', + $ LASTC, K, LASTV-K, ONE, C( K+1, 1 ), LDC, + $ V( K+1, 1 ), LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**H +* + IF( M.GT.K ) THEN +* +* C2 := C2 - V2 * W**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V( K+1, 1 ), LDV, WORK, LDWORK, + $ ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1**H +* + CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**H +* + DO 30 J = 1, K + DO 20 I = 1, LASTC + C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) ) + 20 CONTINUE + 30 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILAZLR( N, K, V, LDV ) ) + LASTC = ILAZLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C1 +* + DO 40 J = 1, K + CALL ZCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 40 CONTINUE +* +* W := W * V1 +* + CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2 +* + CALL ZGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C( 1, K+1 ), LDC, V( K+1, 1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**H +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( K+1, 1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1**H +* + CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 60 J = 1, K + DO 50 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 50 CONTINUE + 60 CONTINUE + END IF +* + ELSE +* +* Let V = ( V1 ) +* ( V2 ) (last K rows) +* where V2 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILAZLR( M, K, V, LDV ) ) + LASTC = ILAZLC( LASTV, N, C, LDC ) +* +* W := C**H * V = (C1**H * V1 + C2**H * V2) (stored in WORK) +* +* W := C2**H +* + DO 70 J = 1, K + CALL ZCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + CALL ZLACGV( LASTC, WORK( 1, J ), 1 ) + 70 CONTINUE +* +* W := W * V2 +* + CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**H*V1 +* + CALL ZGEMM( 'Conjugate transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C, LDC, V, LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V * W**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1 * W**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTV-K, LASTC, K, + $ -ONE, V, LDV, WORK, LDWORK, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**H +* + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**H +* + DO 90 J = 1, K + DO 80 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - + $ DCONJG( WORK( I, J ) ) + 80 CONTINUE + 90 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILAZLR( N, K, V, LDV ) ) + LASTC = ILAZLR( M, LASTV, C, LDC ) +* +* W := C * V = (C1*V1 + C2*V2) (stored in WORK) +* +* W := C2 +* + DO 100 J = 1, K + CALL ZCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 100 CONTINUE +* +* W := W * V2 +* + CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1 +* + CALL ZGEMM( 'No transpose', 'No transpose', + $ LASTC, K, LASTV-K, + $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2**H +* + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( LASTV-K+1, 1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W +* + DO 120 J = 1, K + DO 110 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) + $ - WORK( I, J ) + 110 CONTINUE + 120 CONTINUE + END IF + END IF +* + ELSE IF( LSAME( STOREV, 'R' ) ) THEN +* + IF( LSAME( DIRECT, 'F' ) ) THEN +* +* Let V = ( V1 V2 ) (V1: first K columns) +* where V1 is unit upper triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILAZLC( K, M, V, LDV ) ) + LASTC = ILAZLC( LASTV, N, C, LDC ) +* +* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) +* +* W := C1**H +* + DO 130 J = 1, K + CALL ZCOPY( LASTC, C( J, 1 ), LDC, WORK( 1, J ), 1 ) + CALL ZLACGV( LASTC, WORK( 1, J ), 1 ) + 130 CONTINUE +* +* W := W * V1**H +* + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2**H*V2**H +* + CALL ZGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTC, K, LASTV-K, + $ ONE, C( K+1, 1 ), LDC, V( 1, K+1 ), LDV, + $ ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL ZTRMM( 'Right', 'Upper', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**H * W**H +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - V2**H * W**H +* + CALL ZGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTV-K, LASTC, K, + $ -ONE, V( 1, K+1 ), LDV, WORK, LDWORK, + $ ONE, C( K+1, 1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W**H +* + DO 150 J = 1, K + DO 140 I = 1, LASTC + C( J, I ) = C( J, I ) - DCONJG( WORK( I, J ) ) + 140 CONTINUE + 150 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILAZLC( K, N, V, LDV ) ) + LASTC = ILAZLR( M, LASTV, C, LDC ) +* +* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) +* +* W := C1 +* + DO 160 J = 1, K + CALL ZCOPY( LASTC, C( 1, J ), 1, WORK( 1, J ), 1 ) + 160 CONTINUE +* +* W := W * V1**H +* + CALL ZTRMM( 'Right', 'Upper', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V, LDV, WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C2 * V2**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, K, LASTV-K, ONE, C( 1, K+1 ), LDC, + $ V( 1, K+1 ), LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL ZTRMM( 'Right', 'Upper', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C2 := C2 - W * V2 +* + CALL ZGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, + $ -ONE, WORK, LDWORK, V( 1, K+1 ), LDV, + $ ONE, C( 1, K+1 ), LDC ) + END IF +* +* W := W * V1 +* + CALL ZTRMM( 'Right', 'Upper', 'No transpose', 'Unit', + $ LASTC, K, ONE, V, LDV, WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 180 J = 1, K + DO 170 I = 1, LASTC + C( I, J ) = C( I, J ) - WORK( I, J ) + 170 CONTINUE + 180 CONTINUE +* + END IF +* + ELSE +* +* Let V = ( V1 V2 ) (V2: last K columns) +* where V2 is unit lower triangular. +* + IF( LSAME( SIDE, 'L' ) ) THEN +* +* Form H * C or H**H * C where C = ( C1 ) +* ( C2 ) +* + LASTV = MAX( K, ILAZLC( K, M, V, LDV ) ) + LASTC = ILAZLC( LASTV, N, C, LDC ) +* +* W := C**H * V**H = (C1**H * V1**H + C2**H * V2**H) (stored in WORK) +* +* W := C2**H +* + DO 190 J = 1, K + CALL ZCOPY( LASTC, C( LASTV-K+J, 1 ), LDC, + $ WORK( 1, J ), 1 ) + CALL ZLACGV( LASTC, WORK( 1, J ), 1 ) + 190 CONTINUE +* +* W := W * V2**H +* + CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1**H * V1**H +* + CALL ZGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTC, K, LASTV-K, + $ ONE, C, LDC, V, LDV, ONE, WORK, LDWORK ) + END IF +* +* W := W * T**H or W * T +* + CALL ZTRMM( 'Right', 'Lower', TRANST, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - V**H * W**H +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - V1**H * W**H +* + CALL ZGEMM( 'Conjugate transpose', + $ 'Conjugate transpose', LASTV-K, LASTC, K, + $ -ONE, V, LDV, WORK, LDWORK, ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C2 := C2 - W**H +* + DO 210 J = 1, K + DO 200 I = 1, LASTC + C( LASTV-K+J, I ) = C( LASTV-K+J, I ) - + $ DCONJG( WORK( I, J ) ) + 200 CONTINUE + 210 CONTINUE +* + ELSE IF( LSAME( SIDE, 'R' ) ) THEN +* +* Form C * H or C * H**H where C = ( C1 C2 ) +* + LASTV = MAX( K, ILAZLC( K, N, V, LDV ) ) + LASTC = ILAZLR( M, LASTV, C, LDC ) +* +* W := C * V**H = (C1*V1**H + C2*V2**H) (stored in WORK) +* +* W := C2 +* + DO 220 J = 1, K + CALL ZCOPY( LASTC, C( 1, LASTV-K+J ), 1, + $ WORK( 1, J ), 1 ) + 220 CONTINUE +* +* W := W * V2**H +* + CALL ZTRMM( 'Right', 'Lower', 'Conjugate transpose', + $ 'Unit', LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) + IF( LASTV.GT.K ) THEN +* +* W := W + C1 * V1**H +* + CALL ZGEMM( 'No transpose', 'Conjugate transpose', + $ LASTC, K, LASTV-K, ONE, C, LDC, V, LDV, ONE, + $ WORK, LDWORK ) + END IF +* +* W := W * T or W * T**H +* + CALL ZTRMM( 'Right', 'Lower', TRANS, 'Non-unit', + $ LASTC, K, ONE, T, LDT, WORK, LDWORK ) +* +* C := C - W * V +* + IF( LASTV.GT.K ) THEN +* +* C1 := C1 - W * V1 +* + CALL ZGEMM( 'No transpose', 'No transpose', + $ LASTC, LASTV-K, K, -ONE, WORK, LDWORK, V, LDV, + $ ONE, C, LDC ) + END IF +* +* W := W * V2 +* + CALL ZTRMM( 'Right', 'Lower', 'No transpose', 'Unit', + $ LASTC, K, ONE, V( 1, LASTV-K+1 ), LDV, + $ WORK, LDWORK ) +* +* C1 := C1 - W +* + DO 240 J = 1, K + DO 230 I = 1, LASTC + C( I, LASTV-K+J ) = C( I, LASTV-K+J ) + $ - WORK( I, J ) + 230 CONTINUE + 240 CONTINUE +* + END IF +* + END IF + END IF +* + RETURN +* +* End of ZLARFB +* + END diff --git a/lapack/zlarfg.f b/lapack/zlarfg.f new file mode 100644 index 000000000..a90ae9f74 --- /dev/null +++ b/lapack/zlarfg.f @@ -0,0 +1,203 @@ +*> \brief \b ZLARFG +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARFG + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) +* +* .. Scalar Arguments .. +* INTEGER INCX, N +* COMPLEX*16 ALPHA, TAU +* .. +* .. Array Arguments .. +* COMPLEX*16 X( * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLARFG generates a complex elementary reflector H of order n, such +*> that +*> +*> H**H * ( alpha ) = ( beta ), H**H * H = I. +*> ( x ) ( 0 ) +*> +*> where alpha and beta are scalars, with beta real, and x is an +*> (n-1)-element complex vector. H is represented in the form +*> +*> H = I - tau * ( 1 ) * ( 1 v**H ) , +*> ( v ) +*> +*> where tau is a complex scalar and v is a complex (n-1)-element +*> vector. Note that H is not hermitian. +*> +*> If the elements of x are all zero and alpha is real, then tau = 0 +*> and H is taken to be the unit matrix. +*> +*> Otherwise 1 <= real(tau) <= 2 and abs(tau-1) <= 1 . +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the elementary reflector. +*> \endverbatim +*> +*> \param[in,out] ALPHA +*> \verbatim +*> ALPHA is COMPLEX*16 +*> On entry, the value alpha. +*> On exit, it is overwritten with the value beta. +*> \endverbatim +*> +*> \param[in,out] X +*> \verbatim +*> X is COMPLEX*16 array, dimension +*> (1+(N-2)*abs(INCX)) +*> On entry, the vector x. +*> On exit, it is overwritten with the vector v. +*> \endverbatim +*> +*> \param[in] INCX +*> \verbatim +*> INCX is INTEGER +*> The increment between elements of X. INCX > 0. +*> \endverbatim +*> +*> \param[out] TAU +*> \verbatim +*> TAU is COMPLEX*16 +*> The value tau. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date November 2011 +* +*> \ingroup complex16OTHERauxiliary +* +* ===================================================================== + SUBROUTINE ZLARFG( N, ALPHA, X, INCX, TAU ) +* +* -- LAPACK auxiliary routine (version 3.4.0) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* November 2011 +* +* .. Scalar Arguments .. + INTEGER INCX, N + COMPLEX*16 ALPHA, TAU +* .. +* .. Array Arguments .. + COMPLEX*16 X( * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ONE, ZERO + PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER J, KNT + DOUBLE PRECISION ALPHI, ALPHR, BETA, RSAFMN, SAFMIN, XNORM +* .. +* .. External Functions .. + DOUBLE PRECISION DLAMCH, DLAPY3, DZNRM2 + COMPLEX*16 ZLADIV + EXTERNAL DLAMCH, DLAPY3, DZNRM2, ZLADIV +* .. +* .. Intrinsic Functions .. + INTRINSIC ABS, DBLE, DCMPLX, DIMAG, SIGN +* .. +* .. External Subroutines .. + EXTERNAL ZDSCAL, ZSCAL +* .. +* .. Executable Statements .. +* + IF( N.LE.0 ) THEN + TAU = ZERO + RETURN + END IF +* + XNORM = DZNRM2( N-1, X, INCX ) + ALPHR = DBLE( ALPHA ) + ALPHI = DIMAG( ALPHA ) +* + IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN +* +* H = I +* + TAU = ZERO + ELSE +* +* general case +* + BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' ) + RSAFMN = ONE / SAFMIN +* + KNT = 0 + IF( ABS( BETA ).LT.SAFMIN ) THEN +* +* XNORM, BETA may be inaccurate; scale X and recompute them +* + 10 CONTINUE + KNT = KNT + 1 + CALL ZDSCAL( N-1, RSAFMN, X, INCX ) + BETA = BETA*RSAFMN + ALPHI = ALPHI*RSAFMN + ALPHR = ALPHR*RSAFMN + IF( ABS( BETA ).LT.SAFMIN ) + $ GO TO 10 +* +* New BETA is at most 1, at least SAFMIN +* + XNORM = DZNRM2( N-1, X, INCX ) + ALPHA = DCMPLX( ALPHR, ALPHI ) + BETA = -SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR ) + END IF + TAU = DCMPLX( ( BETA-ALPHR ) / BETA, -ALPHI / BETA ) + ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA-BETA ) + CALL ZSCAL( N-1, ALPHA, X, INCX ) +* +* If ALPHA is subnormal, it may lose relative accuracy +* + DO 20 J = 1, KNT + BETA = BETA*SAFMIN + 20 CONTINUE + ALPHA = BETA + END IF +* + RETURN +* +* End of ZLARFG +* + END diff --git a/lapack/zlarft.f b/lapack/zlarft.f new file mode 100644 index 000000000..6a6151fd0 --- /dev/null +++ b/lapack/zlarft.f @@ -0,0 +1,327 @@ +*> \brief \b ZLARFT +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARFT + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* .. Scalar Arguments .. +* CHARACTER DIRECT, STOREV +* INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. +* COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* +*> \par Purpose: +* ============= +*> +*> \verbatim +*> +*> ZLARFT forms the triangular factor T of a complex block reflector H +*> of order n, which is defined as a product of k elementary reflectors. +*> +*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular; +*> +*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular. +*> +*> If STOREV = 'C', the vector which defines the elementary reflector +*> H(i) is stored in the i-th column of the array V, and +*> +*> H = I - V * T * V**H +*> +*> If STOREV = 'R', the vector which defines the elementary reflector +*> H(i) is stored in the i-th row of the array V, and +*> +*> H = I - V**H * T * V +*> \endverbatim +* +* Arguments: +* ========== +* +*> \param[in] DIRECT +*> \verbatim +*> DIRECT is CHARACTER*1 +*> Specifies the order in which the elementary reflectors are +*> multiplied to form the block reflector: +*> = 'F': H = H(1) H(2) . . . H(k) (Forward) +*> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> \endverbatim +*> +*> \param[in] STOREV +*> \verbatim +*> STOREV is CHARACTER*1 +*> Specifies how the vectors which define the elementary +*> reflectors are stored (see also Further Details): +*> = 'C': columnwise +*> = 'R': rowwise +*> \endverbatim +*> +*> \param[in] N +*> \verbatim +*> N is INTEGER +*> The order of the block reflector H. N >= 0. +*> \endverbatim +*> +*> \param[in] K +*> \verbatim +*> K is INTEGER +*> The order of the triangular factor T (= the number of +*> elementary reflectors). K >= 1. +*> \endverbatim +*> +*> \param[in] V +*> \verbatim +*> V is COMPLEX*16 array, dimension +*> (LDV,K) if STOREV = 'C' +*> (LDV,N) if STOREV = 'R' +*> The matrix V. See further details. +*> \endverbatim +*> +*> \param[in] LDV +*> \verbatim +*> LDV is INTEGER +*> The leading dimension of the array V. +*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K. +*> \endverbatim +*> +*> \param[in] TAU +*> \verbatim +*> TAU is COMPLEX*16 array, dimension (K) +*> TAU(i) must contain the scalar factor of the elementary +*> reflector H(i). +*> \endverbatim +*> +*> \param[out] T +*> \verbatim +*> T is COMPLEX*16 array, dimension (LDT,K) +*> The k by k triangular factor T of the block reflector. +*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is +*> lower triangular. The rest of the array is not used. +*> \endverbatim +*> +*> \param[in] LDT +*> \verbatim +*> LDT is INTEGER +*> The leading dimension of the array T. LDT >= K. +*> \endverbatim +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \date April 2012 +* +*> \ingroup complex16OTHERauxiliary +* +*> \par Further Details: +* ===================== +*> +*> \verbatim +*> +*> The shape of the matrix V and the storage of the vectors which define +*> the H(i) is best illustrated by the following example with n = 5 and +*> k = 3. The elements equal to 1 are not stored. +*> +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> +*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) +*> ( v1 1 ) ( 1 v2 v2 v2 ) +*> ( v1 v2 1 ) ( 1 v3 v3 ) +*> ( v1 v2 v3 ) +*> ( v1 v2 v3 ) +*> +*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> +*> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) +*> ( v1 v2 v3 ) ( v2 v2 v2 1 ) +*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) +*> ( 1 v3 ) +*> ( 1 ) +*> \endverbatim +*> +* ===================================================================== + SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) +* +* -- LAPACK auxiliary routine (version 3.4.1) -- +* -- LAPACK is a software package provided by Univ. of Tennessee, -- +* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- +* April 2012 +* +* .. Scalar Arguments .. + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N +* .. +* .. Array Arguments .. + COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) +* .. +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX*16 ONE, ZERO + PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), + $ ZERO = ( 0.0D+0, 0.0D+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J, PREVLASTV, LASTV +* .. +* .. External Subroutines .. + EXTERNAL ZGEMV, ZLACGV, ZTRMV +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.EQ.0 ) + $ RETURN +* + IF( LSAME( DIRECT, 'F' ) ) THEN + PREVLASTV = N + DO I = 1, K + PREVLASTV = MAX( PREVLASTV, I ) + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = 1, I + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) +* + CALL ZGEMV( 'Conjugate transpose', J-I, I-1, + $ -TAU( I ), V( I+1, 1 ), LDV, + $ V( I+1, I ), 1, ONE, T( 1, I ), 1 ) + ELSE +* Skip any trailing zeros. + DO LASTV = N, I+1, -1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = 1, I-1 + T( J, I ) = -TAU( I ) * V( J , I ) + END DO + J = MIN( LASTV, PREVLASTV ) +* +* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H +* + CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), + $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, + $ ONE, T( 1, I ), LDT ) + END IF +* +* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) +* + CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, + $ LDT, T( 1, I ), 1 ) + T( I, I ) = TAU( I ) + IF( I.GT.1 ) THEN + PREVLASTV = MAX( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + END DO + ELSE + PREVLASTV = 1 + DO I = K, 1, -1 + IF( TAU( I ).EQ.ZERO ) THEN +* +* H(i) = I +* + DO J = I, K + T( J, I ) = ZERO + END DO + ELSE +* +* general case +* + IF( I.LT.K ) THEN + IF( LSAME( STOREV, 'C' ) ) THEN +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( LASTV, I ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) +* + CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I, + $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), + $ 1, ONE, T( I+1, I ), 1 ) + ELSE +* Skip any leading zeros. + DO LASTV = 1, I-1 + IF( V( I, LASTV ).NE.ZERO ) EXIT + END DO + DO J = I+1, K + T( J, I ) = -TAU( I ) * V( J, N-K+I ) + END DO + J = MAX( LASTV, PREVLASTV ) +* +* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H +* + CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ), + $ V( I+1, J ), LDV, V( I, J ), LDV, + $ ONE, T( I+1, I ), LDT ) + END IF +* +* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) +* + CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I, + $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) + IF( I.GT.1 ) THEN + PREVLASTV = MIN( PREVLASTV, LASTV ) + ELSE + PREVLASTV = LASTV + END IF + END IF + T( I, I ) = TAU( I ) + END IF + END DO + END IF + RETURN +* +* End of ZLARFT +* + END -- cgit v1.2.3