diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt
index 676397e76..bfe747563 100644
--- a/SRC/CMakeLists.txt
+++ b/SRC/CMakeLists.txt
@@ -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
diff --git a/SRC/Makefile b/SRC/Makefile
index 0191626f0..13e47020d 100644
--- a/SRC/Makefile
+++ b/SRC/Makefile
@@ -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 \
diff --git a/SRC/VARIANTS/Makefile b/SRC/VARIANTS/Makefile
index 4b0575cc6..0fd8ebf73 100644
--- a/SRC/VARIANTS/Makefile
+++ b/SRC/VARIANTS/Makefile
@@ -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
diff --git a/SRC/VARIANTS/larft/LL-LVL2/clarft.f b/SRC/VARIANTS/larft/LL-LVL2/clarft.f
index 185a5638e..527b970d9 100644
--- a/SRC/VARIANTS/larft/LL-LVL2/clarft.f
+++ b/SRC/VARIANTS/larft/LL-LVL2/clarft.f
@@ -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
*>
*> [TGZ]
*>
@@ -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 @@ SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
END IF
RETURN
*
-* End of CLARFT
+* End of CLARFT_LVL2
*
END
diff --git a/SRC/VARIANTS/larft/LL-LVL2/dlarft.f b/SRC/VARIANTS/larft/LL-LVL2/dlarft.f
index de25fa3e1..6437e31dc 100644
--- a/SRC/VARIANTS/larft/LL-LVL2/dlarft.f
+++ b/SRC/VARIANTS/larft/LL-LVL2/dlarft.f
@@ -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
*>
*> [TGZ]
*>
@@ -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 @@ SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
END IF
RETURN
*
-* End of DLARFT
+* End of DLARFT_LVL2
*
END
diff --git a/SRC/VARIANTS/larft/LL-LVL2/slarft.f b/SRC/VARIANTS/larft/LL-LVL2/slarft.f
index 070f67f56..ffabdd16c 100644
--- a/SRC/VARIANTS/larft/LL-LVL2/slarft.f
+++ b/SRC/VARIANTS/larft/LL-LVL2/slarft.f
@@ -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
*>
*> [TGZ]
*>
@@ -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 @@ SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
END IF
RETURN
*
-* End of SLARFT
+* End of SLARFT_LVL2
*
END
diff --git a/SRC/VARIANTS/larft/LL-LVL2/zlarft.f b/SRC/VARIANTS/larft/LL-LVL2/zlarft.f
index 64cf8b72c..a2e0abc7b 100644
--- a/SRC/VARIANTS/larft/LL-LVL2/zlarft.f
+++ b/SRC/VARIANTS/larft/LL-LVL2/zlarft.f
@@ -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
*>
*> [TGZ]
*>
@@ -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 @@ SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
END IF
RETURN
*
-* End of ZLARFT
+* End of ZLARFT_LVL2
*
END
diff --git a/SRC/clarft.f b/SRC/clarft.f
index fd43bbbce..3bf2448ba 100644
--- a/SRC/clarft.f
+++ b/SRC/clarft.f
@@ -178,11 +178,12 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. 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 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. External Functions..
*
LOGICAL LSAME
- EXTERNAL LSAME
+ INTEGER ILAENV
+ EXTERNAL LSAME, ILAENV
*
* .. Intrinsic Functions..
*
@@ -218,6 +220,14 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV,
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
diff --git a/SRC/clarft_lvl2.f b/SRC/clarft_lvl2.f
new file mode 100644
index 000000000..3b0aea113
--- /dev/null
+++ b/SRC/clarft_lvl2.f
@@ -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
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*
+* 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
diff --git a/SRC/dlarft.f b/SRC/dlarft.f
index 5f87623cf..fa2ea0d13 100644
--- a/SRC/dlarft.f
+++ b/SRC/dlarft.f
@@ -182,7 +182,7 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV,
*
* .. Local Scalars ..
*
- INTEGER I,J,L
+ INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +192,8 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. 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 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV,
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
diff --git a/SRC/dlarft_lvl2.f b/SRC/dlarft_lvl2.f
new file mode 100644
index 000000000..9614df466
--- /dev/null
+++ b/SRC/dlarft_lvl2.f
@@ -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
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*
+* 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
diff --git a/SRC/ilaenv.f b/SRC/ilaenv.f
index 66215c1b3..e108108d7 100644
--- a/SRC/ilaenv.f
+++ b/SRC/ilaenv.f
@@ -666,6 +666,10 @@ INTEGER FUNCTION ILAENV( ISPEC, NAME, OPTS, N1, N2, N3, N4 )
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
diff --git a/SRC/lapack_64.h b/SRC/lapack_64.h
index e8000bf2c..3b9d4275a 100644
--- a/SRC/lapack_64.h
+++ b/SRC/lapack_64.h
@@ -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
diff --git a/SRC/slarft.f b/SRC/slarft.f
index 47fd34988..3e0eac751 100644
--- a/SRC/slarft.f
+++ b/SRC/slarft.f
@@ -182,7 +182,7 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV,
*
* .. Local Scalars ..
*
- INTEGER I,J,L
+ INTEGER I,J,L,NX
LOGICAL QR,LQ,QL,DIRF,COLV
*
* .. External Subroutines ..
@@ -192,7 +192,8 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. 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 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV,
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
diff --git a/SRC/slarft_lvl2.f b/SRC/slarft_lvl2.f
new file mode 100644
index 000000000..7107a91d5
--- /dev/null
+++ b/SRC/slarft_lvl2.f
@@ -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
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*
+* 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
diff --git a/SRC/zlarft.f b/SRC/zlarft.f
index 06cea3860..efd52037d 100644
--- a/SRC/zlarft.f
+++ b/SRC/zlarft.f
@@ -178,11 +178,12 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. 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 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV,
* .. External Functions..
*
LOGICAL LSAME
- EXTERNAL LSAME
+ INTEGER ILAENV
+ EXTERNAL LSAME, ILAENV
*
* .. Intrinsic Functions..
*
@@ -218,6 +220,14 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV,
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
diff --git a/SRC/zlarft_lvl2.f b/SRC/zlarft_lvl2.f
new file mode 100644
index 000000000..808c7fdb2
--- /dev/null
+++ b/SRC/zlarft_lvl2.f
@@ -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
+*>
+*> [TGZ]
+*>
+*> [ZIP]
+*>
+*> [TXT]
+*
+* 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