phynexis-deps/lapack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1163 from jprhyne/dlarftTermCase

Browse Source 
Add non-trivial terminating case to xLARFT
This commit is contained in:
langou
2025-10-23 10:20:45 -06:00
committed by GitHub
parent c08dcfa106 dd076164da
commit 7014148a5b
17 changed files with 1397 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SRC/CMakeLists.txt
+4 -4
View File
@@ -108,7 +108,7 @@ set(SLASRC
slaqgb.f slaqge.f slaqp2.f slaqps.f slaqp2rk.f slaqp3rk.f slaqsb.f slaqsp.f slaqsy.f
slaqr0.f slaqr1.f slaqr2.f slaqr3.f slaqr4.f slaqr5.f
slaqtr.f slar1v.f slar2v.f ilaslr.f ilaslc.f
slarf.f slarf1f.f slarf1l.f slarfb.f slarfb_gett.f slarfg.f slarfgp.f slarft.f slarfx.f slarfy.f
slarf.f slarf1f.f slarf1l.f slarfb.f slarfb_gett.f slarfg.f slarfgp.f slarft.f slarft_lvl2.f slarfx.f slarfy.f
slargv.f slarmm.f slarrv.f slartv.f
slarz.f slarzb.f slarzt.f slasy2.f
slasyf.f slasyf_rook.f slasyf_rk.f slasyf_aa.f
@@ -220,7 +220,7 @@ set(CLASRC
claqhb.f claqhe.f claqhp.f claqp2.f claqps.f claqp2rk.f claqp3rk.f claqsb.f
claqr0.f claqr1.f claqr2.f claqr3.f claqr4.f claqr5.f
claqsp.f claqsy.f clar1v.f clar2v.f ilaclr.f ilaclc.f
clarf.f clarf1f.f clarf1l.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f
clarf.f clarf1f.f clarf1l.f clarfb.f clarfb_gett.f clarfg.f clarfgp.f clarft.f clarft_lvl2.f
clarfx.f clarfy.f clargv.f clarnv.f clarrv.f clartg.f90 clartv.f
clarz.f clarzb.f clarzt.f clascl.f claset.f clasr.f classq.f90
claswp.f clasyf.f clasyf_rook.f clasyf_rk.f clasyf_aa.f
@@ -309,7 +309,7 @@ set(DLASRC
dlaqgb.f dlaqge.f dlaqp2.f dlaqps.f dlaqp2rk.f dlaqp3rk.f dlaqsb.f dlaqsp.f dlaqsy.f
dlaqr0.f dlaqr1.f dlaqr2.f dlaqr3.f dlaqr4.f dlaqr5.f
dlaqtr.f dlar1v.f dlar2v.f iladlr.f iladlc.f
dlarf.f dlarfb.f dlarfb_gett.f dlarfg.f dlarfgp.f dlarft.f dlarfx.f dlarfy.f dlarf1f.f dlarf1l.f
dlarf.f dlarfb.f dlarfb_gett.f dlarfg.f dlarfgp.f dlarft.f dlarft_lvl2.f dlarfx.f dlarfy.f dlarf1f.f dlarf1l.f
dlargv.f dlarmm.f dlarrv.f dlartv.f
dlarz.f dlarzb.f dlarzt.f dlaswp.f dlasy2.f
dlasyf.f dlasyf_rook.f dlasyf_rk.f dlasyf_aa.f
@@ -421,7 +421,7 @@ set(ZLASRC
zlaqr0.f zlaqr1.f zlaqr2.f zlaqr3.f zlaqr4.f zlaqr5.f
zlaqsp.f zlaqsy.f zlar1v.f zlar2v.f ilazlr.f ilazlc.f
zlarcm.f zlarf.f zlarfb.f zlarfb_gett.f zlarf1f.f zlarf1l.f
zlarfg.f zlarfgp.f zlarft.f
zlarfg.f zlarfgp.f zlarft.f zlarft_lvl2.f
zlarfx.f zlarfy.f zlargv.f zlarnv.f zlarrv.f zlartg.f90 zlartv.f
zlarz.f zlarzb.f zlarzt.f zlascl.f zlaset.f zlasr.f
zlassq.f90 zlaswp.f zlasyf.f zlasyf_rook.f zlasyf_rk.f zlasyf_aa.f
SRC/Makefile
+4 -4
View File
@@ -137,7 +137,7 @@ SLASRC = \
slaqgb.o slaqge.o slaqp2.o slaqps.o slaqp2rk.o slaqp3rk.o slaqsb.o slaqsp.o slaqsy.o \
slaqr0.o slaqr1.o slaqr2.o slaqr3.o slaqr4.o slaqr5.o \
slaqtr.o slar1v.o slar2v.o ilaslr.o ilaslc.o \
slarf.o slarf1f.o slarf1l.o slarfb.o slarfb_gett.o slarfg.o slarfgp.o slarft.o slarfx.o slarfy.o \
slarf.o slarf1f.o slarf1l.o slarfb.o slarfb_gett.o slarfg.o slarfgp.o slarft.o slarft_lvl2.o slarfx.o slarfy.o \
slargv.o slarmm.o slarrv.o slartv.o \
slarz.o slarzb.o slarzt.o slaswp.o slasy2.o slasyf.o slasyf_rook.o \
slasyf_rk.o \
@@ -249,7 +249,7 @@ CLASRC = \
claqhb.o claqhe.o claqhp.o claqp2.o claqps.o claqp2rk.o claqp3rk.o claqsb.o \
claqr0.o claqr1.o claqr2.o claqr3.o claqr4.o claqr5.o \
claqsp.o claqsy.o clar1v.o clar2v.o ilaclr.o ilaclc.o \
clarf.o clarf1f.o clarf1l.o clarfb.o clarfb_gett.o clarfg.o clarft.o clarfgp.o \
clarf.o clarf1f.o clarf1l.o clarfb.o clarfb_gett.o clarfg.o clarft.o clarft_lvl2.o clarfgp.o \
clarfx.o clarfy.o clargv.o clarnv.o clarrv.o clartg.o clartv.o \
clarz.o clarzb.o clarzt.o clascl.o claset.o clasr.o classq.o \
claswp.o clasyf.o clasyf_rook.o clasyf_rk.o clasyf_aa.o \
@@ -339,7 +339,7 @@ DLASRC = \
dlaqgb.o dlaqge.o dlaqp2.o dlaqps.o dlaqp2rk.o dlaqp3rk.o dlaqsb.o dlaqsp.o dlaqsy.o \
dlaqr0.o dlaqr1.o dlaqr2.o dlaqr3.o dlaqr4.o dlaqr5.o \
dlaqtr.o dlar1v.o dlar2v.o iladlr.o iladlc.o \
dlarf.o dlarfb.o dlarfb_gett.o dlarfg.o dlarfgp.o dlarft.o dlarfx.o dlarfy.o dlarf1f.o dlarf1l.o\
dlarf.o dlarfb.o dlarfb_gett.o dlarfg.o dlarfgp.o dlarft.o dlarft_lvl2.o dlarfx.o dlarfy.o dlarf1f.o dlarf1l.o\
dlargv.o dlarmm.o dlarrv.o dlartv.o \
dlarz.o dlarzb.o dlarzt.o dlaswp.o dlasy2.o \
dlasyf.o dlasyf_rook.o dlasyf_rk.o \
@@ -454,7 +454,7 @@ ZLASRC = \
zlaqr0.o zlaqr1.o zlaqr2.o zlaqr3.o zlaqr4.o zlaqr5.o \
zlaqsp.o zlaqsy.o zlar1v.o zlar2v.o ilazlr.o ilazlc.o \
zlarcm.o zlarf.o zlarfb.o zlarfb_gett.o zlarf1f.o zlarf1l.o \
zlarfg.o zlarft.o zlarfgp.o \
zlarfg.o zlarft.o zlarft_lvl2.o zlarfgp.o \
zlarfx.o zlarfy.o zlargv.o zlarnv.o zlarrv.o zlartg.o zlartv.o \
zlarz.o zlarzb.o zlarzt.o zlascl.o zlaset.o zlasr.o \
zlassq.o zlaswp.o zlasyf.o zlasyf_rook.o zlasyf_rk.o zlasyf_aa.o \
SRC/VARIANTS/Makefile
+2 -9
View File
@@ -30,11 +30,8 @@ LUREC = lu/REC/cgetrf.o lu/REC/dgetrf.o lu/REC/sgetrf.o lu/REC/zgetrf.o
QRLL = qr/LL/cgeqrf.o qr/LL/dgeqrf.o qr/LL/sgeqrf.o qr/LL/zgeqrf.o
LARFTL2 = larft/LL-LVL2/clarft.o larft/LL-LVL2/dlarft.o larft/LL-LVL2/slarft.o larft/LL-LVL2/zlarft.o
.PHONY: all
all: cholrl.a choltop.a lucr.a lull.a lurec.a qrll.a larftl2.a
all: cholrl.a choltop.a lucr.a lull.a lurec.a qrll.a
cholrl.a: $(CHOLRL)
$(AR) $(ARFLAGS) $@ $^
@@ -60,13 +57,9 @@ qrll.a: $(QRLL)
$(AR) $(ARFLAGS) $@ $^
$(RANLIB) $@
larftl2.a: $(LARFTL2)
$(AR) $(ARFLAGS) $@ $^
$(RANLIB) $@
.PHONY: clean cleanobj cleanlib
clean: cleanobj cleanlib
cleanobj:
rm -f $(CHOLRL) $(CHOLTOP) $(LUCR) $(LULL) $(LUREC) $(QRLL) $(LARFTL2)
rm -f $(CHOLRL) $(CHOLTOP) $(LUCR) $(LULL) $(LUREC) $(QRLL)
cleanlib:
rm -f *.a
SRC/VARIANTS/larft/LL-LVL2/clarft.f
+6 -5
View File
@@ -1,11 +1,11 @@
*> \brief \b CLARFT VARIANT: left-looking Level 2 BLAS version of the algorithm
*> \brief \b CLARFT_LVL2: Level 2 BLAS version for terminating case of CLARFT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download CLARFT + dependencies
*> Download CLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f">
@@ -16,7 +16,8 @@
* Definition:
* ===========
*
* SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
* SUBROUTINE CLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
* T, LDT )
*
* .. Scalar Arguments ..
* CHARACTER DIRECT, STOREV
@@ -32,7 +33,7 @@
*>
*> \verbatim
*>
*> CLARFT forms the triangular factor T of a complex block reflector H
*> CLARFT_LVL2 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;
@@ -322,6 +323,6 @@
END IF
RETURN
*
* End of CLARFT
* End of CLARFT_LVL2
*
END
SRC/VARIANTS/larft/LL-LVL2/dlarft.f
+6 -5
View File
@@ -1,11 +1,11 @@
*> \brief \b DLARFT VARIANT: left-looking Level 2 BLAS version of the algorithm
*> \brief \b DLARFT_LVL2: Level 2 BLAS version for terminating case of DLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download DLARFT + dependencies
*> Download DLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f">
@@ -16,7 +16,8 @@
* Definition:
* ===========
*
* SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
* SUBROUTINE DLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
* T, LDT )
*
* .. Scalar Arguments ..
* CHARACTER DIRECT, STOREV
@@ -32,7 +33,7 @@
*>
*> \verbatim
*>
*> DLARFT forms the triangular factor T of a real block reflector H
*> DLARFT_LVL2 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;
@@ -320,6 +321,6 @@
END IF
RETURN
*
* End of DLARFT
* End of DLARFT_LVL2
*
END
SRC/VARIANTS/larft/LL-LVL2/slarft.f
+6 -5
View File
@@ -1,11 +1,11 @@
*> \brief \b SLARFT VARIANT: left-looking Level 2 BLAS version of the algorithm.
*> \brief \b SLARFT_LVL2: Level 2 BLAS version for terminating case of SLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download SLARFT + dependencies
*> Download SLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f">
@@ -16,7 +16,8 @@
* Definition:
* ===========
*
* SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
* SUBROUTINE SLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
* T, LDT )
*
* .. Scalar Arguments ..
* CHARACTER DIRECT, STOREV
@@ -32,7 +33,7 @@
*>
*> \verbatim
*>
*> SLARFT forms the triangular factor T of a real block reflector H
*> SLARFT_LVL2 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;
@@ -320,6 +321,6 @@
END IF
RETURN
*
* End of SLARFT
* End of SLARFT_LVL2
*
END
SRC/VARIANTS/larft/LL-LVL2/zlarft.f
+6 -5
View File
@@ -1,11 +1,11 @@
*> \brief \b ZLARFT VARIANT: left-looking Level 2 BLAS version of the algorithm.
*> \brief \b ZLARFT_LVL2: Level 2 BLAS version for terminating case of ZLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download ZLARFT + dependencies
*> Download ZLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f">
@@ -16,7 +16,8 @@
* Definition:
* ===========
*
* SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
* SUBROUTINE ZLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
* T, LDT )
*
* .. Scalar Arguments ..
* CHARACTER DIRECT, STOREV
@@ -32,7 +33,7 @@
*>
*> \verbatim
*>
*> ZLARFT forms the triangular factor T of a complex block reflector H
*> ZLARFT_LVL2 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;
@@ -321,6 +322,6 @@
END IF
RETURN
*
* End of ZLARFT
* End of ZLARFT_LVL2
*
END
SRC/clarft.f
+13 -3
View File
@@ -178,11 +178,12 @@
* .. Parameters ..
*
COMPLEX ONE, NEG_ONE, ZERO
PARAMETER(ONE=1.0E+0, ZERO = 0.0E+0, NEG_ONE=-1.0E+0)
PARAMETER(ONE=(1.0E+0,0.0E+0), ZERO = (0.0E+0,0.0E+0),
$ NEG_ONE=(-1.0E+0,0.0E+0))
*
* .. Local Scalars ..
*
INTEGER I,J,L
INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +193,8 @@
* .. External Functions..
*
LOGICAL LSAME
EXTERNAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
*
* .. Intrinsic Functions..
*
@@ -218,6 +220,14 @@
RETURN
END IF
*
* Determine when to cross over into the level 2 based implementation
*
NX = ILAENV(3, "CLARFT", DIRECT // STOREV, N, K, -1, -1)
IF(K.LT.NX) THEN
CALL CLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
RETURN
END IF
*
* Beginning of executable statements
*
L = K / 2
SRC/clarft_lvl2.f
+328
View File
@@ -0,0 +1,328 @@
*> \brief \b CLARFT_LVL2: Level 2 BLAS version for terminating case of CLARFT
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download CLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clarft.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clarft.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE CLARFT_LVL2( 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_LVL2 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.
*
*> \ingroup larft
*
*> \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_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
$ T, LDT )
*
* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER 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 CGEMM, CGEMV, 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_LVL2
*
END
SRC/dlarft.f
+11 -2
View File
@@ -182,7 +182,7 @@
*
* .. Local Scalars ..
*
INTEGER I,J,L
INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +192,8 @@
* .. External Functions..
*
LOGICAL LSAME
EXTERNAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
*
* The general scheme used is inspired by the approach inside DGEQRT3
* which was (at the time of writing this code):
@@ -214,6 +215,14 @@
RETURN
END IF
*
* Determine when to cross over into the level 2 based implementation
*
NX = ILAENV(3, "DLARFT", DIRECT // STOREV, N, K, -1, -1)
IF(K.LT.NX) THEN
CALL DLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
RETURN
END IF
*
* Beginning of executable statements
*
L = K / 2
SRC/dlarft_lvl2.f
+326
View File
@@ -0,0 +1,326 @@
*> \brief \b DLARFT_LVL2: Level 2 BLAS version for terminating case of DLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download DLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarft.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarft.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE DLARFT_LVL2( 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_LVL2 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.
*
*> \ingroup larft
*
*> \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_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
$ T, LDT )
*
* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER 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_LVL2
*
END
SRC/ilaenv.f
+4
View File
@@ -666,6 +666,10 @@
IF( C3.EQ.'HD3' ) THEN
NX = 128
END IF
ELSE IF( C2.EQ.'LA' ) THEN
IF( C3.EQ.'RFT' ) THEN
NX = 64
END IF
END IF
ILAENV = NX
RETURN
SRC/lapack_64.h
+4
View File
@@ -335,6 +335,7 @@
#define CLARFG CLARFG_64
#define CLARFGP CLARFGP_64
#define CLARFT CLARFT_64
#define CLARFT_LVL2 CLARFT_LVL2_64
#define CLARFX CLARFX_64
#define CLARFY CLARFY_64
#define CLARGV CLARGV_64
@@ -809,6 +810,7 @@
#define DLARFG DLARFG_64
#define DLARFGP DLARFGP_64
#define DLARFT DLARFT_64
#define DLARFT_LVL2 DLARFT_LVL2_64
#define DLARFX DLARFX_64
#define DLARFY DLARFY_64
#define DLARGV DLARGV_64
@@ -1404,6 +1406,7 @@
#define SLARFG SLARFG_64
#define SLARFGP SLARFGP_64
#define SLARFT SLARFT_64
#define SLARFT_LVL2 SLARFT_LVL2_64
#define SLARFX SLARFX_64
#define SLARFY SLARFY_64
#define SLARGV SLARGV_64
@@ -2050,6 +2053,7 @@
#define ZLARFG ZLARFG_64
#define ZLARFGP ZLARFGP_64
#define ZLARFT ZLARFT_64
#define ZLARFT_LVL2 ZLARFT_LVL2_64
#define ZLARFX ZLARFX_64
#define ZLARFY ZLARFY_64
#define ZLARGV ZLARGV_64
SRC/slarft.f
+11 -2
View File
@@ -182,7 +182,7 @@
*
* .. Local Scalars ..
*
INTEGER I,J,L
INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +192,8 @@
* .. External Functions..
*
LOGICAL LSAME
EXTERNAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
*
* The general scheme used is inspired by the approach inside DGEQRT3
* which was (at the time of writing this code):
@@ -214,6 +215,14 @@
RETURN
END IF
*
* Determine when to cross over into the level 2 based implementation
*
NX = ILAENV(3, "SLARFT", DIRECT // STOREV, N, K, -1, -1)
IF(K.LT.NX) THEN
CALL SLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
RETURN
END IF
*
* Beginning of executable statements
*
L = K / 2
SRC/slarft_lvl2.f
+326
View File
@@ -0,0 +1,326 @@
*> \brief \b SLARFT_LVL2: Level 2 BLAS version for terminating case of SLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download SLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE SLARFT_LVL2( 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_LVL2 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.
*
*> \ingroup larft
*
*> \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_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
$ T, LDT )
*
* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER 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_LVL2
*
END
SRC/zlarft.f
+13 -3
View File
@@ -178,11 +178,12 @@
* .. Parameters ..
*
COMPLEX*16 ONE, NEG_ONE, ZERO
PARAMETER(ONE=1.0D+0, ZERO = 0.0D+0, NEG_ONE=-1.0D+0)
PARAMETER(ONE=(1.0D+0,0.0D+0), ZERO = (0.0D+0,0.0D+0),
$ NEG_ONE=(-1.0D+0,0.0D+0))
*
* .. Local Scalars ..
*
INTEGER I,J,L
INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +193,8 @@
* .. External Functions..
*
LOGICAL LSAME
EXTERNAL LSAME
INTEGER ILAENV
EXTERNAL LSAME, ILAENV
*
* .. Intrinsic Functions..
*
@@ -218,6 +220,14 @@
RETURN
END IF
*
* Determine when to cross over into the level 2 based implementation
*
NX = ILAENV(3, "ZLARFT", DIRECT // STOREV, N, K, -1, -1)
IF(K.LT.NX) THEN
CALL ZLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT)
RETURN
END IF
*
* Beginning of executable statements
*
L = K / 2
SRC/zlarft_lvl2.f
+327
View File
@@ -0,0 +1,327 @@
*> \brief \b ZLARFT_LVL2: Level 2 BLAS version for terminating case of ZLARFT.
*
* =========== DOCUMENTATION ===========
*
* Online html documentation available at
* http://www.netlib.org/lapack/explore-html/
*
*> Download ZLARFT_LVL2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f">
*> [TXT]</a>
*
* Definition:
* ===========
*
* SUBROUTINE ZLARFT_LVL2( 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_LVL2 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.
*
*> \ingroup larft
*
*> \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_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU,
$ T, LDT )
*
* -- LAPACK auxiliary routine --
* -- LAPACK is a software package provided by Univ. of Tennessee, --
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
* .. Scalar Arguments ..
CHARACTER 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, ZTRMV, ZGEMM
* ..
* .. 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_LVL2
*
END