diff --git a/SRC/CMakeLists.txt b/SRC/CMakeLists.txt index a6b42f552..be426cecd 100644 --- a/SRC/CMakeLists.txt +++ b/SRC/CMakeLists.txt @@ -307,7 +307,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 + dlarf.f dlarfb.f dlarfb_gett.f dlarfg.f dlarfgp.f dlarft.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 @@ -418,7 +418,7 @@ set(ZLASRC zlaqhb.f zlaqhe.f zlaqhp.f zlaqp2.f zlaqps.f zlaqp2rk.f zlaqp3rk.f zlaqsb.f 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 + zlarcm.f zlarf.f zlarfb.f zlarfb_gett.f zlarf1f.f zlarf1l.f zlarfg.f zlarfgp.f zlarft.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 diff --git a/SRC/Makefile b/SRC/Makefile index 7fb2a5670..0191626f0 100644 --- a/SRC/Makefile +++ b/SRC/Makefile @@ -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 \ + dlarf.o dlarfb.o dlarfb_gett.o dlarfg.o dlarfgp.o dlarft.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 \ @@ -453,7 +453,7 @@ ZLASRC = \ zlaqhb.o zlaqhe.o zlaqhp.o zlaqp2.o zlaqps.o zlaqp2rk.o zlaqp3rk.o zlaqsb.o \ 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 \ + zlarcm.o zlarf.o zlarfb.o zlarfb_gett.o zlarf1f.o zlarf1l.o \ zlarfg.o zlarft.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 \ diff --git a/SRC/dgebd2.f b/SRC/dgebd2.f index 98559fade..4994d7480 100644 --- a/SRC/dgebd2.f +++ b/SRC/dgebd2.f @@ -202,14 +202,14 @@ SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ZERO, ONE - PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D+0 ) * .. * .. Local Scalars .. INTEGER I * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1F, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -242,15 +242,13 @@ SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1, $ TAUQ( I ) ) D( I ) = A( I, I ) - A( I, I ) = ONE * * Apply H(i) to A(i:m,i+1:n) from the left * IF( I.LT.N ) - $ CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ CALL DLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, $ TAUQ( I ), $ A( I, I+1 ), LDA, WORK ) - A( I, I ) = D( I ) * IF( I.LT.N ) THEN * @@ -260,13 +258,11 @@ SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL DLARFG( N-I, A( I, I+1 ), A( I, MIN( I+2, N ) ), $ LDA, TAUP( I ) ) E( I ) = A( I, I+1 ) - A( I, I+1 ) = ONE * * Apply G(i) to A(i+1:m,i+1:n) from the right * - CALL DLARF( 'Right', M-I, N-I, A( I, I+1 ), LDA, + CALL DLARF1F( 'Right', M-I, N-I, A( I, I+1 ), LDA, $ TAUP( I ), A( I+1, I+1 ), LDA, WORK ) - A( I, I+1 ) = E( I ) ELSE TAUP( I ) = ZERO END IF @@ -283,14 +279,12 @@ SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) $ LDA, $ TAUP( I ) ) D( I ) = A( I, I ) - A( I, I ) = ONE * * Apply G(i) to A(i+1:m,i:n) from the right * IF( I.LT.M ) - $ CALL DLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, + $ CALL DLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, $ TAUP( I ), A( I+1, I ), LDA, WORK ) - A( I, I ) = D( I ) * IF( I.LT.M ) THEN * @@ -301,14 +295,12 @@ SUBROUTINE DGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) $ 1, $ TAUQ( I ) ) E( I ) = A( I+1, I ) - A( I+1, I ) = ONE * * Apply H(i) to A(i+1:m,i+1:n) from the left * - CALL DLARF( 'Left', M-I, N-I, A( I+1, I ), 1, + CALL DLARF1F( 'Left', M-I, N-I, A( I+1, I ), 1, $ TAUQ( I ), $ A( I+1, I+1 ), LDA, WORK ) - A( I+1, I ) = E( I ) ELSE TAUQ( I ) = ZERO END IF diff --git a/SRC/dgehd2.f b/SRC/dgehd2.f index 62417aa8c..4b18455a5 100644 --- a/SRC/dgehd2.f +++ b/SRC/dgehd2.f @@ -166,10 +166,9 @@ SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1F, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -199,20 +198,17 @@ SUBROUTINE DGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) * CALL DLARFG( IHI-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, $ TAU( I ) ) - AII = A( I+1, I ) - A( I+1, I ) = ONE * * Apply H(i) to A(1:ihi,i+1:ihi) from the right * - CALL DLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ), + CALL DLARF1F( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ), $ A( 1, I+1 ), LDA, WORK ) * * Apply H(i) to A(i+1:ihi,i+1:n) from the left * - CALL DLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1, TAU( I ), + CALL DLARF1F( 'Left', IHI-I, N-I, A( I+1, I ), 1, TAU( I ), $ A( I+1, I+1 ), LDA, WORK ) * - A( I+1, I ) = AII 10 CONTINUE * RETURN diff --git a/SRC/dgelq2.f b/SRC/dgelq2.f index 31dfc07a1..79af37b4e 100644 --- a/SRC/dgelq2.f +++ b/SRC/dgelq2.f @@ -146,10 +146,9 @@ SUBROUTINE DGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1F, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -183,12 +182,9 @@ SUBROUTINE DGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(i+1:m,i:n) from the right * - AII = A( I, I ) - A( I, I ) = ONE - CALL DLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, + CALL DLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, $ TAU( I ), $ A( I+1, I ), LDA, WORK ) - A( I, I ) = AII END IF 10 CONTINUE RETURN diff --git a/SRC/dgeql2.f b/SRC/dgeql2.f index dfb08ff31..18c14618f 100644 --- a/SRC/dgeql2.f +++ b/SRC/dgeql2.f @@ -140,10 +140,9 @@ SUBROUTINE DGEQL2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1L, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -177,12 +176,9 @@ SUBROUTINE DGEQL2( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(1:m-k+i,1:n-k+i-1) from the left * - AII = A( M-K+I, N-K+I ) - A( M-K+I, N-K+I ) = ONE - CALL DLARF( 'Left', M-K+I, N-K+I-1, A( 1, N-K+I ), 1, + CALL DLARF1L( 'Left', M-K+I, N-K+I-1, A( 1, N-K+I ), 1, $ TAU( I ), $ A, LDA, WORK ) - A( M-K+I, N-K+I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/dgeqp3rk.f b/SRC/dgeqp3rk.f index e14ea95c0..2bd82ca9b 100644 --- a/SRC/dgeqp3rk.f +++ b/SRC/dgeqp3rk.f @@ -670,7 +670,7 @@ SUBROUTINE DGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA, * 1) DGEQP3RK and DLAQP2RK: 2*N to store full and partial * column 2-norms. * 2) DLAQP2RK: N+NRHS-1 to use in WORK array that is used -* in DLARF subroutine inside DLAQP2RK to apply an +* in DLARF1F subroutine inside DLAQP2RK to apply an * elementary reflector from the left. * TOTAL_WORK_SIZE = 3*N + NRHS - 1 * @@ -686,7 +686,7 @@ SUBROUTINE DGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA, * 1) DGEQP3RK, DLAQP2RK, DLAQP3RK: 2*N to store full and * partial column 2-norms. * 2) DLAQP2RK: N+NRHS-1 to use in WORK array that is used -* in DLARF subroutine to apply an elementary reflector +* in DLARF1F subroutine to apply an elementary reflector * from the left. * 3) DLAQP3RK: NB*(N+NRHS) to use in the work array F that * is used to apply a block reflector from diff --git a/SRC/dgeqr2.f b/SRC/dgeqr2.f index bd4facfce..f42f8decb 100644 --- a/SRC/dgeqr2.f +++ b/SRC/dgeqr2.f @@ -147,10 +147,9 @@ SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1F, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -184,11 +183,8 @@ SUBROUTINE DGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(i:m,i+1:n) from the left * - AII = A( I, I ) - A( I, I ) = ONE - CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + CALL DLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), $ A( I, I+1 ), LDA, WORK ) - A( I, I ) = AII END IF 10 CONTINUE RETURN diff --git a/SRC/dgeqr2p.f b/SRC/dgeqr2p.f index b2f3188f3..b6f910f50 100644 --- a/SRC/dgeqr2p.f +++ b/SRC/dgeqr2p.f @@ -151,10 +151,9 @@ SUBROUTINE DGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFGP, XERBLA + EXTERNAL DLARF1F, DLARFGP, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -188,11 +187,8 @@ SUBROUTINE DGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(i:m,i+1:n) from the left * - AII = A( I, I ) - A( I, I ) = ONE - CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + CALL DLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), $ A( I, I+1 ), LDA, WORK ) - A( I, I ) = AII END IF 10 CONTINUE RETURN diff --git a/SRC/dgerq2.f b/SRC/dgerq2.f index a4ef46d85..7106fe029 100644 --- a/SRC/dgerq2.f +++ b/SRC/dgerq2.f @@ -140,10 +140,9 @@ SUBROUTINE DGERQ2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - DOUBLE PRECISION AII * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, XERBLA + EXTERNAL DLARF1L, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -177,11 +176,8 @@ SUBROUTINE DGERQ2( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right * - AII = A( M-K+I, N-K+I ) - A( M-K+I, N-K+I ) = ONE - CALL DLARF( 'Right', M-K+I-1, N-K+I, A( M-K+I, 1 ), LDA, + CALL DLARF1L( 'Right', M-K+I-1, N-K+I, A( M-K+I, 1 ), LDA, $ TAU( I ), A, LDA, WORK ) - A( M-K+I, N-K+I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/dlaqp2.f b/SRC/dlaqp2.f index 5bfa967ee..a4edebf64 100644 --- a/SRC/dlaqp2.f +++ b/SRC/dlaqp2.f @@ -168,7 +168,7 @@ SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MN, OFFPI, PVT - DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z + DOUBLE PRECISION TEMP, TEMP2, TOL3Z * .. * .. External Subroutines .. EXTERNAL DLARF, DLARFG, DSWAP @@ -219,11 +219,8 @@ SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, * * Apply H(i)**T to A(offset+i:m,i+1:n) from the left. * - AII = A( OFFPI, I ) - A( OFFPI, I ) = ONE - CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, + CALL DLARF1F( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, $ TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) ) - A( OFFPI, I ) = AII END IF * * Update partial column norms. diff --git a/SRC/dlaqp2rk.f b/SRC/dlaqp2rk.f index aecd6bb69..fd6d7ba63 100644 --- a/SRC/dlaqp2rk.f +++ b/SRC/dlaqp2rk.f @@ -253,7 +253,7 @@ *> \param[out] WORK *> \verbatim *> WORK is DOUBLE PRECISION array, dimension (N-1) -*> Used in DLARF subroutine to apply an elementary +*> Used in DLARF1F subroutine to apply an elementary *> reflector from the left. *> \endverbatim *> @@ -367,10 +367,10 @@ SUBROUTINE DLAQP2RK( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL, * .. Local Scalars .. INTEGER I, ITEMP, J, JMAXC2NRM, KK, KP, MINMNFACT, $ MINMNUPDT - DOUBLE PRECISION AIKK, HUGEVAL, TEMP, TEMP2, TOL3Z + DOUBLE PRECISION HUGEVAL, TEMP, TEMP2, TOL3Z * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFG, DSWAP + EXTERNAL DLARF1F, DLARFG, DSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT @@ -621,11 +621,8 @@ SUBROUTINE DLAQP2RK( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL, * condition is satisfied, not only KK < N+NRHS ) * IF( KK.LT.MINMNUPDT ) THEN - AIKK = A( I, KK ) - A( I, KK ) = ONE - CALL DLARF( 'Left', M-I+1, N+NRHS-KK, A( I, KK ), 1, + CALL DLARF1F( 'Left', M-I+1, N+NRHS-KK, A( I, KK ), 1, $ TAU( KK ), A( I, KK+1 ), LDA, WORK( 1 ) ) - A( I, KK ) = AIKK END IF * IF( KK.LT.MINMNFACT ) THEN diff --git a/SRC/dlaqr2.f b/SRC/dlaqr2.f index 8591d5d3b..a7a8392d1 100644 --- a/SRC/dlaqr2.f +++ b/SRC/dlaqr2.f @@ -312,7 +312,7 @@ SUBROUTINE DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, * .. External Subroutines .. EXTERNAL DCOPY, DGEHRD, DGEMM, DLACPY, $ DLAHQR, - $ DLANV2, DLARF, DLARFG, DLASET, DORMHR, DTREXC + $ DLANV2, DLARF1F, DLARFG, DLASET, DORMHR, DTREXC * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT @@ -597,16 +597,15 @@ SUBROUTINE DLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, CALL DCOPY( NS, V, LDV, WORK, 1 ) BETA = WORK( 1 ) CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU ) - WORK( 1 ) = ONE * CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), $ LDT ) * - CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, + CALL DLARF1F( 'L', NS, JW, WORK, 1, TAU, T, LDT, $ WORK( JW+1 ) ) - CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, + CALL DLARF1F( 'R', NS, NS, WORK, 1, TAU, T, LDT, $ WORK( JW+1 ) ) - CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, + CALL DLARF1F( 'R', JW, NS, WORK, 1, TAU, V, LDV, $ WORK( JW+1 ) ) * CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), diff --git a/SRC/dlaqr3.f b/SRC/dlaqr3.f index f8b0bc266..5ec235134 100644 --- a/SRC/dlaqr3.f +++ b/SRC/dlaqr3.f @@ -310,7 +310,7 @@ SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, * .. External Subroutines .. EXTERNAL DCOPY, DGEHRD, DGEMM, DLACPY, DLAHQR, $ DLANV2, - $ DLAQR4, DLARF, DLARFG, DLASET, DORMHR, DTREXC + $ DLAQR4, DLARF1F, DLARFG, DLASET, DORMHR, DTREXC * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, INT, MAX, MIN, SQRT @@ -608,16 +608,15 @@ SUBROUTINE DLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, CALL DCOPY( NS, V, LDV, WORK, 1 ) BETA = WORK( 1 ) CALL DLARFG( NS, BETA, WORK( 2 ), 1, TAU ) - WORK( 1 ) = ONE * CALL DLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), $ LDT ) * - CALL DLARF( 'L', NS, JW, WORK, 1, TAU, T, LDT, + CALL DLARF1F( 'L', NS, JW, WORK, 1, TAU, T, LDT, $ WORK( JW+1 ) ) - CALL DLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, + CALL DLARF1F( 'R', NS, NS, WORK, 1, TAU, T, LDT, $ WORK( JW+1 ) ) - CALL DLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, + CALL DLARF1F( 'R', JW, NS, WORK, 1, TAU, V, LDV, $ WORK( JW+1 ) ) * CALL DGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), diff --git a/SRC/dlarf1f.f b/SRC/dlarf1f.f new file mode 100644 index 000000000..104122999 --- /dev/null +++ b/SRC/dlarf1f.f @@ -0,0 +1,292 @@ +*> \brief \b DLARF1F applies an elementary reflector to a general rectangular +* matrix assuming v(1) = 1. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARF1F + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARF1F( 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 +*> +*> DLARF1F 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. V(1) is not referenced or modified. +*> \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 +* +* To take advantage of the fact that v(1) = 1, we do the following +* v = [ 1 v_2 ]**T +* If SIDE='L' +* |-----| +* | C_1 | +* C =| C_2 | +* |-----| +* C_1\in\mathbb{R}^{1\times n}, C_2\in\mathbb{R}^{m-1\times n} +* So we compute: +* C = HC = (I - \tau vv**T)C +* = C - \tau vv**T C +* w = C**T v = [ C_1**T C_2**T ] [ 1 v_2 ]**T +* = C_1**T + C_2**T v ( DGEMM then DAXPY ) +* C = C - \tau vv**T C +* = C - \tau vw**T +* Giving us C_1 = C_1 - \tau w**T ( DAXPY ) +* and +* C_2 = C_2 - \tau v_2w**T ( DGER ) +* If SIDE='R' +* +* C = [ C_1 C_2 ] +* C_1\in\mathbb{R}^{m\times 1}, C_2\in\mathbb{R}^{m\times n-1} +* So we compute: +* C = CH = C(I - \tau vv**T) +* = C - \tau Cvv**T +* +* w = Cv = [ C_1 C_2 ] [ 1 v_2 ]**T +* = C_1 + C_2v_2 ( DGEMM then DAXPY ) +* C = C - \tau Cvv**T +* = C - \tau wv**T +* Giving us C_1 = C_1 - \tau w ( DAXPY ) +* and +* C_2 = C_2 - \tau wv_2**T ( DGER ) +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup larf +* +* ===================================================================== + SUBROUTINE DLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- 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 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, DAXPY, DSCAL +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILADLR, ILADLC + EXTERNAL LSAME, ILADLR, ILADLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 1 + 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. +! Since we are assuming that V(1) = 1, and it is not stored, so we +! shouldn't access it. + DO WHILE( LASTV.GT.1 .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 + IF( LASTC.EQ.0 ) THEN + RETURN + END IF + IF( APPLYLEFT ) THEN +* +* Form H * C +* + ! Check if lastv = 1. This means v = 1, So we just need to compute + ! C := HC = (1-\tau)C. + IF( LASTV.EQ.1 ) THEN +* +* C(1,1:lastc) := ( 1 - tau ) * C(1,1:lastc) +* + CALL DSCAL(LASTC, ONE - TAU, C, LDC) + ELSE +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1) +* + ! w(1:lastc,1) := C(2:lastv,1:lastc)**T * v(2:lastv,1) + CALL DGEMV( 'Transpose', LASTV-1, LASTC, ONE, C(1+1,1), + $ LDC, V(1+INCV), INCV, ZERO, WORK, 1) + ! w(1:lastc,1) += C(1,1:lastc)**T * v(1,1) = C(1,1:lastc)**T + CALL DAXPY(LASTC, ONE, C, LDC, WORK, 1) +* +* C(1:lastv,1:lastc) := C(...) - tau * v(1:lastv,1) * w(1:lastc,1)**T +* + ! C(1, 1:lastc) := C(...) - tau * v(1,1) * w(1:lastc,1)**T + ! = C(...) - tau * w(1:lastc,1)**T + CALL DAXPY(LASTC, -TAU, WORK, 1, C, LDC) + ! C(2:lastv,1:lastc) := C(...) - tau * v(2:lastv,1)*w(1:lastc,1)**T + CALL DGER(LASTV-1, LASTC, -TAU, V(1+INCV), INCV, WORK, 1, + $ C(1+1,1), LDC) + END IF + ELSE +* +* Form C * H +* + ! Check if n = 1. This means v = 1, so we just need to compute + ! C := CH = C(1-\tau). + IF( LASTV.EQ.1 ) THEN +* +* C(1:lastc,1) := ( 1 - tau ) * C(1:lastc,1) +* + CALL DSCAL(LASTC, ONE - TAU, C, 1) + ELSE +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + ! w(1:lastc,1) := C(1:lastc,2:lastv) * v(2:lastv,1) + CALL DGEMV( 'No transpose', LASTC, LASTV-1, ONE, + $ C(1,1+1), LDC, V(1+INCV), INCV, ZERO, WORK, 1 ) + ! w(1:lastc,1) += C(1:lastc,1) v(1,1) = C(1:lastc,1) + CALL DAXPY(LASTC, ONE, C, 1, WORK, 1) +* +* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T +* + ! C(1:lastc,1) := C(...) - tau * w(1:lastc,1) * v(1,1)**T + ! = C(...) - tau * w(1:lastc,1) + CALL DAXPY(LASTC, -TAU, WORK, 1, C, 1) + ! C(1:lastc,2:lastv) := C(...) - tau * w(1:lastc,1) * v(2:lastv)**T + CALL DGER( LASTC, LASTV-1, -TAU, WORK, 1, V(1+INCV), + $ INCV, C(1,1+1), LDC ) + END IF + END IF + RETURN +* +* End of DLARF1F +* + END diff --git a/SRC/dlarf1l.f b/SRC/dlarf1l.f new file mode 100644 index 000000000..80a486f79 --- /dev/null +++ b/SRC/dlarf1l.f @@ -0,0 +1,252 @@ +*> \brief \b DLARF1L applies an elementary reflector to a general rectangular +* matrix assuming v(lastv) = 1 where lastv is the last non-zero +* element +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download DLARF1L + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE DLARF1L( 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 +*> +*> DLARF1L 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. +* +*> \ingroup larf +* +* ===================================================================== + SUBROUTINE DLARF1L( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- 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 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, FIRSTV, LASTV, LASTC +* .. +* .. External Subroutines .. + EXTERNAL DGEMV, DGER +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILADLR, ILADLC + EXTERNAL LSAME, ILADLR, ILADLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + FIRSTV = 1 + 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 + I = 1 +! Look for the last non-zero row in V. + DO WHILE( LASTV.GT.FIRSTV .AND. V( I ).EQ.ZERO ) + FIRSTV = FIRSTV + 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 + IF( LASTC.EQ.0 ) THEN + RETURN + END IF + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.GT.0 ) THEN + ! Check if m = 1. This means v = 1, So we just need to compute + ! C := HC = (1-\tau)C. + IF( LASTV.EQ.FIRSTV ) THEN + CALL DSCAL(LASTC, ONE - TAU, C( FIRSTV, 1), LDC) + ELSE +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1) +* + ! w(1:lastc,1) := C(1:lastv-1,1:lastc)**T * v(1:lastv-1,1) + CALL DGEMV( 'Transpose', LASTV-FIRSTV, LASTC, ONE, + $ C(FIRSTV,1), LDC, V(I), INCV, ZERO, + $ WORK, 1) + ! w(1:lastc,1) += C(lastv,1:lastc)**T * v(lastv,1) = C(lastv,1:lastc)**T + CALL DAXPY(LASTC, ONE, C(LASTV,1), LDC, WORK, 1) +* +* C(1:lastv,1:lastc) := C(...) - tau * v(1:lastv,1) * w(1:lastc,1)**T +* + ! C(lastv, 1:lastc) := C(...) - tau * v(lastv,1) * w(1:lastc,1)**T + ! = C(...) - tau * w(1:lastc,1)**T + CALL DAXPY(LASTC, -TAU, WORK, 1, C(LASTV,1), LDC) + ! C(1:lastv-1,1:lastc) := C(...) - tau * v(1:lastv-1,1)*w(1:lastc,1)**T + CALL DGER(LASTV-FIRSTV, LASTC, -TAU, V(I), INCV, + $ WORK, 1, C(FIRSTV,1), LDC) + END IF + END IF + ELSE +* +* Form C * H +* + IF( LASTV.GT.0 ) THEN + ! Check if n = 1. This means v = 1, so we just need to compute + ! C := CH = C(1-\tau). + IF( LASTV.EQ.FIRSTV ) THEN + CALL DSCAL(LASTC, ONE - TAU, C, 1) + ELSE +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + ! w(1:lastc,1) := C(1:lastc,1:lastv-1) * v(1:lastv-1,1) + CALL DGEMV( 'No transpose', LASTC, LASTV-FIRSTV, + $ ONE, C(1,FIRSTV), LDC, V(I), INCV, ZERO, WORK, 1 ) + ! w(1:lastc,1) += C(1:lastc,lastv) * v(lastv,1) = C(1:lastc,lastv) + CALL DAXPY(LASTC, ONE, C(1,LASTV), 1, WORK, 1) +* +* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T +* + ! C(1:lastc,lastv) := C(...) - tau * w(1:lastc,1) * v(lastv,1)**T + ! = C(...) - tau * w(1:lastc,1) + CALL DAXPY(LASTC, -TAU, WORK, 1, C(1,LASTV), 1) + ! C(1:lastc,1:lastv-1) := C(...) - tau * w(1:lastc,1) * v(1:lastv-1)**T + CALL DGER( LASTC, LASTV-FIRSTV, -TAU, WORK, 1, V(I), + $ INCV, C(1,FIRSTV), LDC ) + END IF + END IF + END IF + RETURN +* +* End of DLARF1L +* + END diff --git a/SRC/dopmtr.f b/SRC/dopmtr.f index fd2fd9c23..6289a2dae 100644 --- a/SRC/dopmtr.f +++ b/SRC/dopmtr.f @@ -261,12 +261,9 @@ SUBROUTINE DOPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, * * Apply H(i) * - AII = AP( II ) - AP( II ) = ONE - CALL DLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, + CALL DLARF1L( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, $ LDC, $ WORK ) - AP( II ) = AII * IF( FORWRD ) THEN II = II + I + 2 diff --git a/SRC/dorbdb.f b/SRC/dorbdb.f index f70813bdb..c796db13f 100644 --- a/SRC/dorbdb.f +++ b/SRC/dorbdb.f @@ -316,7 +316,7 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, DOUBLE PRECISION Z1, Z2, Z3, Z4 * .. * .. External Subroutines .. - EXTERNAL DAXPY, DLARF, DLARFGP, DSCAL, + EXTERNAL DAXPY, DLARF1F, DLARFGP, DSCAL, $ XERBLA * .. * .. External Functions .. @@ -422,7 +422,6 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, ELSE IF( P .EQ. I ) THEN CALL DLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) ) END IF - X11(I,I) = ONE IF ( M-P .GT. I ) THEN CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, $ TAUP2(I) ) @@ -430,25 +429,23 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1, $ TAUP2(I) ) END IF - X21(I,I) = ONE * IF ( Q .GT. I ) THEN - CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), + CALL DLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), $ X11(I,I+1), LDX11, WORK ) END IF IF ( M-Q+1 .GT. I ) THEN - CALL DLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, + CALL DLARF1F( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, $ TAUP1(I), $ X12(I,I), LDX12, WORK ) END IF IF ( Q .GT. I ) THEN - CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), - $ X21(I,I+1), LDX21, WORK ) + CALL DLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, + $ TAUP2(I), X21(I,I+1), LDX21, WORK ) END IF IF ( M-Q+1 .GT. I ) THEN - CALL DLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, - $ TAUP2(I), - $ X22(I,I), LDX22, WORK ) + CALL DLARF1F( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, + $ TAUP2(I), X22(I,I), LDX22, WORK ) END IF * IF( I .LT. Q ) THEN @@ -476,7 +473,6 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, $ TAUQ1(I) ) END IF - X11(I,I+1) = ONE END IF IF ( Q+I-1 .LT. M ) THEN IF ( M-Q .EQ. I ) THEN @@ -487,23 +483,22 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, $ TAUQ2(I) ) END IF END IF - X12(I,I) = ONE * IF( I .LT. Q ) THEN - CALL DLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, + CALL DLARF1F( 'R', P-I, Q-I, X11(I,I+1), LDX11, $ TAUQ1(I), $ X11(I+1,I+1), LDX11, WORK ) - CALL DLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, + CALL DLARF1F( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, $ TAUQ1(I), $ X21(I+1,I+1), LDX21, WORK ) END IF IF ( P .GT. I ) THEN - CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, + CALL DLARF1F( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, $ TAUQ2(I), $ X12(I+1,I), LDX12, WORK ) END IF IF ( M-P .GT. I ) THEN - CALL DLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, + CALL DLARF1F( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, $ TAUQ2(I), X22(I+1,I), LDX22, WORK ) END IF * @@ -521,15 +516,14 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, $ TAUQ2(I) ) END IF - X12(I,I) = ONE * IF ( P .GT. I ) THEN - CALL DLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, + CALL DLARF1F( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, $ TAUQ2(I), $ X12(I+1,I), LDX12, WORK ) END IF IF( M-P-Q .GE. 1 ) - $ CALL DLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, + $ CALL DLARF1F( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) * END DO @@ -546,9 +540,8 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), $ LDX22, TAUQ2(P+I) ) END IF - X22(Q+I,P+I) = ONE IF ( I .LT. M-P-Q ) THEN - CALL DLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), + CALL DLARF1F( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), $ LDX22, $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) END IF @@ -584,7 +577,6 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, * CALL DLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, $ TAUP1(I) ) - X11(I,I) = ONE IF ( I .EQ. M-P ) THEN CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21, $ TAUP2(I) ) @@ -592,24 +584,23 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, $ TAUP2(I) ) END IF - X21(I,I) = ONE * IF ( Q .GT. I ) THEN - CALL DLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, + CALL DLARF1F( 'R', Q-I, P-I+1, X11(I,I), LDX11, $ TAUP1(I), $ X11(I+1,I), LDX11, WORK ) END IF IF ( M-Q+1 .GT. I ) THEN - CALL DLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, + CALL DLARF1F( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, $ TAUP1(I), X12(I,I), LDX12, WORK ) END IF IF ( Q .GT. I ) THEN - CALL DLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, + CALL DLARF1F( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, $ TAUP2(I), $ X21(I+1,I), LDX21, WORK ) END IF IF ( M-Q+1 .GT. I ) THEN - CALL DLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, + CALL DLARF1F( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, $ TAUP2(I), X22(I,I), LDX22, WORK ) END IF * @@ -634,7 +625,6 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, $ TAUQ1(I) ) END IF - X11(I+1,I) = ONE END IF IF ( M-Q .GT. I ) THEN CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, @@ -643,20 +633,18 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I,I), 1, $ TAUQ2(I) ) END IF - X12(I,I) = ONE * IF( I .LT. Q ) THEN - CALL DLARF( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I), + CALL DLARF1F( 'L', Q-I, P-I, X11(I+1,I), 1, TAUQ1(I), $ X11(I+1,I+1), LDX11, WORK ) - CALL DLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, TAUQ1(I), - $ X21(I+1,I+1), LDX21, WORK ) + CALL DLARF1F( 'L', Q-I, M-P-I, X11(I+1,I), 1, + $ TAUQ1(I), X21(I+1,I+1), LDX21, WORK ) END IF - CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), - $ X12(I,I+1), LDX12, WORK ) + CALL DLARF1F( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), + $ X12(I,I+1), LDX12, WORK ) IF ( M-P-I .GT. 0 ) THEN - CALL DLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, - $ TAUQ2(I), - $ X22(I,I+1), LDX22, WORK ) + CALL DLARF1F( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, + $ TAUQ2(I), X22(I,I+1), LDX22, WORK ) END IF * END DO @@ -668,16 +656,14 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DSCAL( M-Q-I+1, -Z1*Z4, X12(I,I), 1 ) CALL DLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, $ TAUQ2(I) ) - X12(I,I) = ONE * IF ( P .GT. I ) THEN - CALL DLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, TAUQ2(I), - $ X12(I,I+1), LDX12, WORK ) + CALL DLARF1F( 'L', M-Q-I+1, P-I, X12(I,I), 1, + $ TAUQ2(I), X12(I,I+1), LDX12, WORK ) END IF IF( M-P-Q .GE. 1 ) - $ CALL DLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, - $ TAUQ2(I), - $ X22(I,Q+1), LDX22, WORK ) + $ CALL DLARF1F( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, + $ TAUQ2(I), X22(I,Q+1), LDX22, WORK ) * END DO * @@ -694,10 +680,10 @@ SUBROUTINE DORBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL DLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), $ 1, $ TAUQ2(P+I) ) - CALL DLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, - $ TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, WORK ) + CALL DLARF1F( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), + $ 1, TAUQ2(P+I), X22(P+I,Q+I+1), LDX22, + $ WORK ) END IF - X22(P+I,Q+I) = ONE * END DO * diff --git a/SRC/dorbdb1.f b/SRC/dorbdb1.f index 7dceca9bc..1972ef4bc 100644 --- a/SRC/dorbdb1.f +++ b/SRC/dorbdb1.f @@ -228,7 +228,7 @@ SUBROUTINE DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFGP, DORBDB5, DROT, + EXTERNAL DLARF1F, DLARFGP, DORBDB5, DROT, $ XERBLA * .. * .. External Functions .. @@ -287,12 +287,10 @@ SUBROUTINE DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( X21(I,I), X11(I,I) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - X11(I,I) = ONE - X21(I,I) = ONE - CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, + CALL DLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, $ I+1), $ LDX11, WORK(ILARF) ) - CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), + CALL DLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), $ X21(I,I+1), LDX21, WORK(ILARF) ) * IF( I .LT. Q ) THEN @@ -301,11 +299,11 @@ SUBROUTINE DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, CALL DLARFGP( Q-I, X21(I,I+1), X21(I,I+2), LDX21, $ TAUQ1(I) ) S = X21(I,I+1) - X21(I,I+1) = ONE - CALL DLARF( 'R', P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), + CALL DLARF1F( 'R', P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), $ X11(I+1,I+1), LDX11, WORK(ILARF) ) - CALL DLARF( 'R', M-P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), - $ X21(I+1,I+1), LDX21, WORK(ILARF) ) + CALL DLARF1F( 'R', M-P-I, Q-I, X21(I,I+1), LDX21, + $ TAUQ1(I), X21(I+1,I+1), LDX21, + $ WORK(ILARF) ) C = SQRT( DNRM2( P-I, X11(I+1,I+1), 1 )**2 $ + DNRM2( M-P-I, X21(I+1,I+1), 1 )**2 ) PHI(I) = ATAN2( S, C ) diff --git a/SRC/dorbdb2.f b/SRC/dorbdb2.f index e6b1b9710..e66a40d6b 100644 --- a/SRC/dorbdb2.f +++ b/SRC/dorbdb2.f @@ -227,7 +227,7 @@ SUBROUTINE DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFGP, DORBDB5, DROT, DSCAL, + EXTERNAL DLARF1F, DLARFGP, DORBDB5, DROT, DSCAL, $ XERBLA * .. * .. External Functions .. @@ -287,11 +287,10 @@ SUBROUTINE DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, END IF CALL DLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) ) C = X11(I,I) - X11(I,I) = ONE - CALL DLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + CALL DLARF1F( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL DLARF( 'R', M-P-I+1, Q-I+1, X11(I,I), LDX11, TAUQ1(I), - $ X21(I,I), LDX21, WORK(ILARF) ) + CALL DLARF1F( 'R', M-P-I+1, Q-I+1, X11(I,I), LDX11, + $ TAUQ1(I), X21(I,I), LDX21, WORK(ILARF) ) S = SQRT( DNRM2( P-I, X11(I+1,I), 1 )**2 $ + DNRM2( M-P-I+1, X21(I,I), 1 )**2 ) THETA(I) = ATAN2( S, C ) @@ -306,12 +305,10 @@ SUBROUTINE DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI(I) = ATAN2( X11(I+1,I), X21(I,I) ) C = COS( PHI(I) ) S = SIN( PHI(I) ) - X11(I+1,I) = ONE - CALL DLARF( 'L', P-I, Q-I, X11(I+1,I), 1, TAUP1(I), + CALL DLARF1F( 'L', P-I, Q-I, X11(I+1,I), 1, TAUP1(I), $ X11(I+1,I+1), LDX11, WORK(ILARF) ) END IF - X21(I,I) = ONE - CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), + CALL DLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), $ X21(I,I+1), LDX21, WORK(ILARF) ) * END DO @@ -320,8 +317,7 @@ SUBROUTINE DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, * DO I = P + 1, Q CALL DLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) - X21(I,I) = ONE - CALL DLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), + CALL DLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, TAUP2(I), $ X21(I,I+1), LDX21, WORK(ILARF) ) END DO * diff --git a/SRC/dorbdb3.f b/SRC/dorbdb3.f index 1914ce496..3c64d5226 100644 --- a/SRC/dorbdb3.f +++ b/SRC/dorbdb3.f @@ -226,7 +226,7 @@ SUBROUTINE DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL DLARF, DLARFGP, DORBDB5, DROT, + EXTERNAL DLARF1F, DLARFGP, DORBDB5, DROT, $ XERBLA * .. * .. External Functions .. @@ -287,10 +287,9 @@ SUBROUTINE DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, * CALL DLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) ) S = X21(I,I) - X21(I,I) = ONE - CALL DLARF( 'R', P-I+1, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + CALL DLARF1F( 'R', P-I+1, Q-I+1, X21(I,I), LDX21, TAUQ1(I), $ X11(I,I), LDX11, WORK(ILARF) ) - CALL DLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + CALL DLARF1F( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), $ X21(I+1,I), LDX21, WORK(ILARF) ) C = SQRT( DNRM2( P-I+1, X11(I,I), 1 )**2 $ + DNRM2( M-P-I, X21(I+1,I), 1 )**2 ) @@ -306,12 +305,10 @@ SUBROUTINE DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI(I) = ATAN2( X21(I+1,I), X11(I,I) ) C = COS( PHI(I) ) S = SIN( PHI(I) ) - X21(I+1,I) = ONE - CALL DLARF( 'L', M-P-I, Q-I, X21(I+1,I), 1, TAUP2(I), + CALL DLARF1F( 'L', M-P-I, Q-I, X21(I+1,I), 1, TAUP2(I), $ X21(I+1,I+1), LDX21, WORK(ILARF) ) END IF - X11(I,I) = ONE - CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, + CALL DLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, $ I+1), $ LDX11, WORK(ILARF) ) * @@ -321,8 +318,7 @@ SUBROUTINE DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, * DO I = M-P + 1, Q CALL DLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) - X11(I,I) = ONE - CALL DLARF( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, + CALL DLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, TAUP1(I), X11(I, $ I+1), $ LDX11, WORK(ILARF) ) END DO diff --git a/SRC/dorbdb4.f b/SRC/dorbdb4.f index c0150cb96..446ccc686 100644 --- a/SRC/dorbdb4.f +++ b/SRC/dorbdb4.f @@ -307,13 +307,10 @@ SUBROUTINE DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( PHANTOM(1), PHANTOM(P+1) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - PHANTOM(1) = ONE - PHANTOM(P+1) = ONE - CALL DLARF( 'L', P, Q, PHANTOM(1), 1, TAUP1(1), X11, - $ LDX11, - $ WORK(ILARF) ) - CALL DLARF( 'L', M-P, Q, PHANTOM(P+1), 1, TAUP2(1), X21, - $ LDX21, WORK(ILARF) ) + CALL DLARF1F( 'L', P, Q, PHANTOM(1), 1, TAUP1(1), X11, + $ LDX11, WORK(ILARF) ) + CALL DLARF1F( 'L', M-P, Q, PHANTOM(P+1), 1, TAUP2(1), + $ X21, LDX21, WORK(ILARF) ) ELSE CALL DORBDB5( P-I+1, M-P-I+1, Q-I+1, X11(I,I-1), 1, $ X21(I,I-1), 1, X11(I,I), LDX11, X21(I,I), @@ -326,21 +323,18 @@ SUBROUTINE DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( X11(I,I-1), X21(I,I-1) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - X11(I,I-1) = ONE - X21(I,I-1) = ONE - CALL DLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, TAUP1(I), - $ X11(I,I), LDX11, WORK(ILARF) ) - CALL DLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, TAUP2(I), - $ X21(I,I), LDX21, WORK(ILARF) ) + CALL DLARF1F( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, TAUP1(I), + $ X11(I,I), LDX11, WORK(ILARF) ) + CALL DLARF1F( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, + $ TAUP2(I), X21(I,I), LDX21, WORK(ILARF) ) END IF * CALL DROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C ) CALL DLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) ) C = X21(I,I) - X21(I,I) = ONE - CALL DLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + CALL DLARF1F( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL DLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + CALL DLARF1F( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), $ X21(I+1,I), LDX21, WORK(ILARF) ) IF( I .LT. M-Q ) THEN S = SQRT( DNRM2( P-I, X11(I+1,I), 1 )**2 @@ -354,10 +348,9 @@ SUBROUTINE DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, * DO I = M - Q + 1, P CALL DLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) ) - X11(I,I) = ONE - CALL DLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + CALL DLARF1F( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL DLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + CALL DLARF1F( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I), $ X21(M-Q+1,I), LDX21, WORK(ILARF) ) END DO * @@ -367,8 +360,7 @@ SUBROUTINE DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, CALL DLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), $ LDX21, $ TAUQ1(I) ) - X21(M-Q+I-P,I) = ONE - CALL DLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, + CALL DLARF1F( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, $ TAUQ1(I), $ X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) ) END DO diff --git a/SRC/dorg2l.f b/SRC/dorg2l.f index af8831085..b95430239 100644 --- a/SRC/dorg2l.f +++ b/SRC/dorg2l.f @@ -133,7 +133,7 @@ SUBROUTINE DORG2L( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, II, J, L * .. * .. External Subroutines .. - EXTERNAL DLARF, DSCAL, XERBLA + EXTERNAL DLARF1L, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -176,8 +176,8 @@ SUBROUTINE DORG2L( M, N, K, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(1:m-k+i,1:n-k+i) from the left * - A( M-N+II, II ) = ONE - CALL DLARF( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), + !A(M-N+II, II) = ONE + CALL DLARF1L( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), $ A, $ LDA, WORK ) CALL DSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 ) diff --git a/SRC/dorg2r.f b/SRC/dorg2r.f index 221b52bb8..0d293486f 100644 --- a/SRC/dorg2r.f +++ b/SRC/dorg2r.f @@ -133,7 +133,7 @@ SUBROUTINE DORG2R( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, J, L * .. * .. External Subroutines .. - EXTERNAL DLARF, DSCAL, XERBLA + EXTERNAL DLARF1F, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -176,8 +176,7 @@ SUBROUTINE DORG2R( M, N, K, A, LDA, TAU, WORK, INFO ) * Apply H(i) to A(i:m,i:n) from the left * IF( I.LT.N ) THEN - A( I, I ) = ONE - CALL DLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + CALL DLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), $ A( I, I+1 ), LDA, WORK ) END IF IF( I.LT.M ) diff --git a/SRC/dorgl2.f b/SRC/dorgl2.f index 98128b25d..3d71e3326 100644 --- a/SRC/dorgl2.f +++ b/SRC/dorgl2.f @@ -132,7 +132,7 @@ SUBROUTINE DORGL2( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, J, L * .. * .. External Subroutines .. - EXTERNAL DLARF, DSCAL, XERBLA + EXTERNAL DLARF1F, DSCAL, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -180,8 +180,7 @@ SUBROUTINE DORGL2( M, N, K, A, LDA, TAU, WORK, INFO ) * IF( I.LT.N ) THEN IF( I.LT.M ) THEN - A( I, I ) = ONE - CALL DLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, + CALL DLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, $ TAU( I ), A( I+1, I ), LDA, WORK ) END IF CALL DSCAL( N-I, -TAU( I ), A( I, I+1 ), LDA ) diff --git a/SRC/dorgr2.f b/SRC/dorgr2.f index 211db4f8d..ae317489b 100644 --- a/SRC/dorgr2.f +++ b/SRC/dorgr2.f @@ -180,8 +180,8 @@ SUBROUTINE DORGR2( M, N, K, A, LDA, TAU, WORK, INFO ) * * Apply H(i) to A(1:m-k+i,1:n-k+i) from the right * - A( II, N-M+II ) = ONE - CALL DLARF( 'Right', II-1, N-M+II, A( II, 1 ), LDA, + !A( II, N-M+II ) = ONE + CALL DLARF1L( 'Right', II-1, N-M+II, A( II, 1 ), LDA, $ TAU( I ), $ A, LDA, WORK ) CALL DSCAL( N-M+II-1, -TAU( I ), A( II, 1 ), LDA ) diff --git a/SRC/dorm2l.f b/SRC/dorm2l.f index b1a27ab21..37866ed0c 100644 --- a/SRC/dorm2l.f +++ b/SRC/dorm2l.f @@ -98,7 +98,6 @@ *> The i-th column must contain the vector which defines the *> elementary reflector H(i), for i = 1,2,...,k, as returned by *> DGEQLF in the last k columns of its array argument A. -*> A is modified by the routine but restored on exit. *> \endverbatim *> *> \param[in] LDA @@ -178,14 +177,13 @@ SUBROUTINE DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, MI, NI, NQ - DOUBLE PRECISION AII * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL DLARF, XERBLA + EXTERNAL DLARF1L, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -262,11 +260,8 @@ SUBROUTINE DORM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * * Apply H(i) * - AII = A( NQ-K+I, I ) - A( NQ-K+I, I ) = ONE - CALL DLARF( SIDE, MI, NI, A( 1, I ), 1, TAU( I ), C, LDC, + CALL DLARF1L( SIDE, MI, NI, A( 1, I ), 1, TAU( I ), C, LDC, $ WORK ) - A( NQ-K+I, I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/dorm2r.f b/SRC/dorm2r.f index d894a806c..e6e40a718 100644 --- a/SRC/dorm2r.f +++ b/SRC/dorm2r.f @@ -98,7 +98,6 @@ *> The i-th column must contain the vector which defines the *> elementary reflector H(i), for i = 1,2,...,k, as returned by *> DGEQRF in the first k columns of its array argument A. -*> A is modified by the routine but restored on exit. *> \endverbatim *> *> \param[in] LDA @@ -178,14 +177,13 @@ SUBROUTINE DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ - DOUBLE PRECISION AII * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL DLARF, XERBLA + EXTERNAL XERBLA, DLARF1F * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -266,12 +264,9 @@ SUBROUTINE DORM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * * Apply H(i) * - AII = A( I, I ) - A( I, I ) = ONE - CALL DLARF( SIDE, MI, NI, A( I, I ), 1, TAU( I ), C( IC, + CALL DLARF1F( SIDE, MI, NI, A( I, I ), 1, TAU( I ), C( IC, $ JC ), $ LDC, WORK ) - A( I, I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/dorml2.f b/SRC/dorml2.f index fcdf5b1b1..dfdab973c 100644 --- a/SRC/dorml2.f +++ b/SRC/dorml2.f @@ -100,7 +100,6 @@ *> The i-th row must contain the vector which defines the *> elementary reflector H(i), for i = 1,2,...,k, as returned by *> DGELQF in the first k rows of its array argument A. -*> A is modified by the routine but restored on exit. *> \endverbatim *> *> \param[in] LDA @@ -178,14 +177,13 @@ SUBROUTINE DORML2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ - DOUBLE PRECISION AII * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL DLARF, XERBLA + EXTERNAL DLARF1F, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -266,11 +264,8 @@ SUBROUTINE DORML2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * * Apply H(i) * - AII = A( I, I ) - A( I, I ) = ONE - CALL DLARF( SIDE, MI, NI, A( I, I ), LDA, TAU( I ), + CALL DLARF1F( SIDE, MI, NI, A( I, I ), LDA, TAU( I ), $ C( IC, JC ), LDC, WORK ) - A( I, I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/lapack_64.h b/SRC/lapack_64.h index 24e4fd148..4d7318d97 100644 --- a/SRC/lapack_64.h +++ b/SRC/lapack_64.h @@ -796,6 +796,8 @@ #define DLAR1V DLAR1V_64 #define DLAR2V DLAR2V_64 #define DLARF DLARF_64 +#define DLARF1F DLARF1F_64 +#define DLARF1L DLARF1L_64 #define DLARFB DLARFB_64 #define DLARFB_GETT DLARFB_GETT_64 #define DLARFG DLARFG_64 @@ -2029,6 +2031,8 @@ #define ZLAR2V ZLAR2V_64 #define ZLARCM ZLARCM_64 #define ZLARF ZLARF_64 +#define ZLARF1F ZLARF1F_64 +#define ZLARF1L ZLARF1L_64 #define ZLARFB ZLARFB_64 #define ZLARFB_GETT ZLARFB_GETT_64 #define ZLARFG ZLARFG_64 diff --git a/SRC/zgebd2.f b/SRC/zgebd2.f index ec1142954..510ab849b 100644 --- a/SRC/zgebd2.f +++ b/SRC/zgebd2.f @@ -202,16 +202,14 @@ SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) * ===================================================================== * * .. Parameters .. - COMPLEX*16 ZERO, ONE - PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ), - $ ONE = ( 1.0D+0, 0.0D+0 ) ) -* .. + COMPLEX*16 ZERO + PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. Local Scalars .. INTEGER I COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLACGV, ZLARF1F, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN @@ -245,12 +243,11 @@ SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL ZLARFG( M-I+1, ALPHA, A( MIN( I+1, M ), I ), 1, $ TAUQ( I ) ) D( I ) = DBLE( ALPHA ) - A( I, I ) = ONE * * Apply H(i)**H to A(i:m,i+1:n) from the left * IF( I.LT.N ) - $ CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, + $ CALL ZLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, $ DCONJG( TAUQ( I ) ), A( I, I+1 ), LDA, WORK ) A( I, I ) = D( I ) * @@ -264,11 +261,10 @@ SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL ZLARFG( N-I, ALPHA, A( I, MIN( I+2, N ) ), LDA, $ TAUP( I ) ) E( I ) = DBLE( ALPHA ) - A( I, I+1 ) = ONE * * Apply G(i) to A(i+1:m,i+1:n) from the right * - CALL ZLARF( 'Right', M-I, N-I, A( I, I+1 ), LDA, + CALL ZLARF1F( 'Right', M-I, N-I, A( I, I+1 ), LDA, $ TAUP( I ), A( I+1, I+1 ), LDA, WORK ) CALL ZLACGV( N-I, A( I, I+1 ), LDA ) A( I, I+1 ) = E( I ) @@ -289,12 +285,11 @@ SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL ZLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA, $ TAUP( I ) ) D( I ) = DBLE( ALPHA ) - A( I, I ) = ONE * * Apply G(i) to A(i+1:m,i:n) from the right * IF( I.LT.M ) - $ CALL ZLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, + $ CALL ZLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, $ TAUP( I ), A( I+1, I ), LDA, WORK ) CALL ZLACGV( N-I+1, A( I, I ), LDA ) A( I, I ) = D( I ) @@ -308,11 +303,10 @@ SUBROUTINE ZGEBD2( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO ) CALL ZLARFG( M-I, ALPHA, A( MIN( I+2, M ), I ), 1, $ TAUQ( I ) ) E( I ) = DBLE( ALPHA ) - A( I+1, I ) = ONE * * Apply H(i)**H to A(i+1:m,i+1:n) from the left * - CALL ZLARF( 'Left', M-I, N-I, A( I+1, I ), 1, + CALL ZLARF1F( 'Left', M-I, N-I, A( I+1, I ), 1, $ DCONJG( TAUQ( I ) ), A( I+1, I+1 ), LDA, $ WORK ) A( I+1, I ) = E( I ) diff --git a/SRC/zgehd2.f b/SRC/zgehd2.f index 63c9fce1c..22c03157b 100644 --- a/SRC/zgehd2.f +++ b/SRC/zgehd2.f @@ -166,10 +166,9 @@ SUBROUTINE ZGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLARF1F, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN @@ -197,22 +196,19 @@ SUBROUTINE ZGEHD2( N, ILO, IHI, A, LDA, TAU, WORK, INFO ) * * Compute elementary reflector H(i) to annihilate A(i+2:ihi,i) * - ALPHA = A( I+1, I ) - CALL ZLARFG( IHI-I, ALPHA, A( MIN( I+2, N ), I ), 1, + CALL ZLARFG( IHI-I, A( I+1, I ), A( MIN( I+2, N ), I ), 1, $ TAU( I ) ) - A( I+1, I ) = ONE * * Apply H(i) to A(1:ihi,i+1:ihi) from the right * - CALL ZLARF( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ), - $ A( 1, I+1 ), LDA, WORK ) + CALL ZLARF1F( 'Right', IHI, IHI-I, A( I+1, I ), 1, TAU( I ), + $ A( 1, I+1 ), LDA, WORK ) * * Apply H(i)**H to A(i+1:ihi,i+1:n) from the left * - CALL ZLARF( 'Left', IHI-I, N-I, A( I+1, I ), 1, - $ DCONJG( TAU( I ) ), A( I+1, I+1 ), LDA, WORK ) + CALL ZLARF1F( 'Left', IHI-I, N-I, A( I+1, I ), 1, + $ CONJG( TAU( I ) ), A( I+1, I+1 ), LDA, WORK ) * - A( I+1, I ) = ALPHA 10 CONTINUE * RETURN diff --git a/SRC/zgelq2.f b/SRC/zgelq2.f index bd3521caa..f84756167 100644 --- a/SRC/zgelq2.f +++ b/SRC/zgelq2.f @@ -146,10 +146,9 @@ SUBROUTINE ZGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLACGV, ZLARF1F, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -178,19 +177,16 @@ SUBROUTINE ZGELQ2( M, N, A, LDA, TAU, WORK, INFO ) * Generate elementary reflector H(i) to annihilate A(i,i+1:n) * CALL ZLACGV( N-I+1, A( I, I ), LDA ) - ALPHA = A( I, I ) - CALL ZLARFG( N-I+1, ALPHA, A( I, MIN( I+1, N ) ), LDA, + CALL ZLARFG( N-I+1, A( I, I ), A( I, MIN( I+1, N ) ), LDA, $ TAU( I ) ) IF( I.LT.M ) THEN * * Apply H(i) to A(i+1:m,i:n) from the right * - A( I, I ) = ONE - CALL ZLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, - $ TAU( I ), - $ A( I+1, I ), LDA, WORK ) + CALL ZLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, + $ TAU( I ), + $ A( I+1, I ), LDA, WORK ) END IF - A( I, I ) = ALPHA CALL ZLACGV( N-I+1, A( I, I ), LDA ) 10 CONTINUE RETURN diff --git a/SRC/zgeql2.f b/SRC/zgeql2.f index cdac186e9..d972c4ba1 100644 --- a/SRC/zgeql2.f +++ b/SRC/zgeql2.f @@ -140,10 +140,9 @@ SUBROUTINE ZGEQL2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLARF1L, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN @@ -172,15 +171,13 @@ SUBROUTINE ZGEQL2( M, N, A, LDA, TAU, WORK, INFO ) * Generate elementary reflector H(i) to annihilate * A(1:m-k+i-1,n-k+i) * - ALPHA = A( M-K+I, N-K+I ) - CALL ZLARFG( M-K+I, ALPHA, A( 1, N-K+I ), 1, TAU( I ) ) + CALL ZLARFG( M-K+I, A( M-K+I, N-K+I ), A( 1, N-K+I ), 1, + $ TAU( I ) ) * * Apply H(i)**H to A(1:m-k+i,1:n-k+i-1) from the left * - A( M-K+I, N-K+I ) = ONE - CALL ZLARF( 'Left', M-K+I, N-K+I-1, A( 1, N-K+I ), 1, - $ DCONJG( TAU( I ) ), A, LDA, WORK ) - A( M-K+I, N-K+I ) = ALPHA + CALL ZLARF1L( 'Left', M-K+I, N-K+I-1, A( 1, N-K+I ), 1, + $ CONJG( TAU( I ) ), A, LDA, WORK ) 10 CONTINUE RETURN * diff --git a/SRC/zgeqp3rk.f b/SRC/zgeqp3rk.f index 654093e31..322741b21 100644 --- a/SRC/zgeqp3rk.f +++ b/SRC/zgeqp3rk.f @@ -677,7 +677,7 @@ SUBROUTINE ZGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA, * Minimal workspace size in case of using only unblocked * BLAS 2 code in ZLAQP2RK. * 1) ZLAQP2RK: N+NRHS-1 to use in WORK array that is used -* in ZLARF subroutine inside ZLAQP2RK to apply an +* in ZLARF1F subroutine inside ZLAQP2RK to apply an * elementary reflector from the left. * TOTAL_WORK_SIZE = 3*N + NRHS - 1 * @@ -693,7 +693,7 @@ SUBROUTINE ZGEQP3RK( M, N, NRHS, KMAX, ABSTOL, RELTOL, A, LDA, * 1) ZGEQP3RK, ZLAQP2RK, ZLAQP3RK: 2*N to store full and * partial column 2-norms. * 2) ZLAQP2RK: N+NRHS-1 to use in WORK array that is used -* in ZLARF subroutine to apply an elementary reflector +* in ZLARF1F subroutine to apply an elementary reflector * from the left. * 3) ZLAQP3RK: NB*(N+NRHS) to use in the work array F that * is used to apply a block reflector from diff --git a/SRC/zgeqr2.f b/SRC/zgeqr2.f index 457404ad9..784d76e61 100644 --- a/SRC/zgeqr2.f +++ b/SRC/zgeqr2.f @@ -147,10 +147,9 @@ SUBROUTINE ZGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLARF1F, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN @@ -184,11 +183,8 @@ SUBROUTINE ZGEQR2( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i)**H to A(i:m,i+1:n) from the left * - ALPHA = A( I, I ) - A( I, I ) = ONE - CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, - $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) - A( I, I ) = ALPHA + CALL ZLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, + $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) END IF 10 CONTINUE RETURN diff --git a/SRC/zgeqr2p.f b/SRC/zgeqr2p.f index 93451faec..b5827494a 100644 --- a/SRC/zgeqr2p.f +++ b/SRC/zgeqr2p.f @@ -151,10 +151,9 @@ SUBROUTINE ZGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZLARFGP + EXTERNAL XERBLA, ZLARF1F, ZLARFGP * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX, MIN @@ -188,11 +187,8 @@ SUBROUTINE ZGEQR2P( M, N, A, LDA, TAU, WORK, INFO ) * * Apply H(i)**H to A(i:m,i+1:n) from the left * - ALPHA = A( I, I ) - A( I, I ) = ONE - CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, - $ DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) - A( I, I ) = ALPHA + CALL ZLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, + $ CONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK ) END IF 10 CONTINUE RETURN diff --git a/SRC/zgerq2.f b/SRC/zgerq2.f index f7bac6d57..0df8cc320 100644 --- a/SRC/zgerq2.f +++ b/SRC/zgerq2.f @@ -140,10 +140,9 @@ SUBROUTINE ZGERQ2( M, N, A, LDA, TAU, WORK, INFO ) * .. * .. Local Scalars .. INTEGER I, K - COMPLEX*16 ALPHA * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF, ZLARFG + EXTERNAL XERBLA, ZLACGV, ZLARF1L, ZLARFG * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN @@ -173,15 +172,13 @@ SUBROUTINE ZGERQ2( M, N, A, LDA, TAU, WORK, INFO ) * A(m-k+i,1:n-k+i-1) * CALL ZLACGV( N-K+I, A( M-K+I, 1 ), LDA ) - ALPHA = A( M-K+I, N-K+I ) - CALL ZLARFG( N-K+I, ALPHA, A( M-K+I, 1 ), LDA, TAU( I ) ) + CALL ZLARFG( N-K+I, A( M-K+I, N-K+I ), A( M-K+I, 1 ), LDA, + $ TAU( I ) ) * * Apply H(i) to A(1:m-k+i-1,1:n-k+i) from the right * - A( M-K+I, N-K+I ) = ONE - CALL ZLARF( 'Right', M-K+I-1, N-K+I, A( M-K+I, 1 ), LDA, - $ TAU( I ), A, LDA, WORK ) - A( M-K+I, N-K+I ) = ALPHA + CALL ZLARF1L( 'Right', M-K+I-1, N-K+I, A( M-K+I, 1 ), LDA, + $ TAU( I ), A, LDA, WORK ) CALL ZLACGV( N-K+I-1, A( M-K+I, 1 ), LDA ) 10 CONTINUE RETURN diff --git a/SRC/zlaqp2.f b/SRC/zlaqp2.f index ee47f5119..7189307bb 100644 --- a/SRC/zlaqp2.f +++ b/SRC/zlaqp2.f @@ -171,10 +171,9 @@ SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, * .. Local Scalars .. INTEGER I, ITEMP, J, MN, OFFPI, PVT DOUBLE PRECISION TEMP, TEMP2, TOL3Z - COMPLEX*16 AII * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFG, ZSWAP + EXTERNAL ZLARF1F, ZLARFG, ZSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, DCONJG, MAX, MIN, SQRT @@ -222,12 +221,9 @@ SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, * * Apply H(i)**H to A(offset+i:m,i+1:n) from the left. * - AII = A( OFFPI, I ) - A( OFFPI, I ) = CONE - CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, - $ DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA, - $ WORK( 1 ) ) - A( OFFPI, I ) = AII + CALL ZLARF1F( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, + $ CONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA, + $ WORK( 1 ) ) END IF * * Update partial column norms. diff --git a/SRC/zlaqp2rk.f b/SRC/zlaqp2rk.f index f6bf555c2..9fe92fc5a 100644 --- a/SRC/zlaqp2rk.f +++ b/SRC/zlaqp2rk.f @@ -254,7 +254,7 @@ *> \param[out] WORK *> \verbatim *> WORK is COMPLEX*16 array, dimension (N-1) -*> Used in ZLARF subroutine to apply an elementary +*> Used in ZLARF1F subroutine to apply an elementary *> reflector from the left. *> \endverbatim *> @@ -372,10 +372,9 @@ SUBROUTINE ZLAQP2RK( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL, INTEGER I, ITEMP, J, JMAXC2NRM, KK, KP, MINMNFACT, $ MINMNUPDT DOUBLE PRECISION HUGEVAL, TAUNAN, TEMP, TEMP2, TOL3Z - COMPLEX*16 AIKK * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFG, ZSWAP + EXTERNAL ZLARF1F, ZLARFG, ZSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT @@ -633,12 +632,9 @@ SUBROUTINE ZLAQP2RK( M, N, NRHS, IOFFSET, KMAX, ABSTOL, RELTOL, * condition is satisfied, not only KK < N+NRHS ) * IF( KK.LT.MINMNUPDT ) THEN - AIKK = A( I, KK ) - A( I, KK ) = CONE - CALL ZLARF( 'Left', M-I+1, N+NRHS-KK, A( I, KK ), 1, - $ DCONJG( TAU( KK ) ), A( I, KK+1 ), LDA, - $ WORK( 1 ) ) - A( I, KK ) = AIKK + CALL ZLARF1F( 'Left', M-I+1, N+NRHS-KK, A( I, KK ), 1, + $ CONJG( TAU( KK ) ), A( I, KK+1 ), LDA, + $ WORK( 1 ) ) END IF * IF( KK.LT.MINMNFACT ) THEN diff --git a/SRC/zlaqr2.f b/SRC/zlaqr2.f index e29c3875a..a107e2a41 100644 --- a/SRC/zlaqr2.f +++ b/SRC/zlaqr2.f @@ -293,7 +293,7 @@ SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 ) * .. * .. Local Scalars .. - COMPLEX*16 BETA, CDUM, S, TAU + COMPLEX*16 CDUM, S, TAU DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN, $ KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWKOPT @@ -305,7 +305,7 @@ SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, * .. External Subroutines .. EXTERNAL ZCOPY, ZGEHRD, ZGEMM, ZLACPY, $ ZLAHQR, - $ ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR + $ ZLARF1F, ZLARFG, ZLASET, ZTREXC, ZUNMHR * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN @@ -476,19 +476,17 @@ SUBROUTINE ZLAQR2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, DO 50 I = 1, NS WORK( I ) = DCONJG( WORK( I ) ) 50 CONTINUE - BETA = WORK( 1 ) - CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU ) - WORK( 1 ) = ONE + CALL ZLARFG( NS, WORK( 1 ), WORK( 2 ), 1, TAU ) * CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), $ LDT ) * - CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT, - $ WORK( JW+1 ) ) - CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, - $ WORK( JW+1 ) ) - CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, - $ WORK( JW+1 ) ) + CALL ZLARF1F( 'L', NS, JW, WORK, 1, CONJG( TAU ), T, LDT, + $ WORK( JW+1 ) ) + CALL ZLARF1F( 'R', NS, NS, WORK, 1, TAU, T, LDT, + $ WORK( JW+1 ) ) + CALL ZLARF1F( 'R', JW, NS, WORK, 1, TAU, V, LDV, + $ WORK( JW+1 ) ) * CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), $ LWORK-JW, INFO ) diff --git a/SRC/zlaqr3.f b/SRC/zlaqr3.f index a6f962611..b217cf0a6 100644 --- a/SRC/zlaqr3.f +++ b/SRC/zlaqr3.f @@ -290,7 +290,7 @@ SUBROUTINE ZLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, PARAMETER ( RZERO = 0.0d0, RONE = 1.0d0 ) * .. * .. Local Scalars .. - COMPLEX*16 BETA, CDUM, S, TAU + COMPLEX*16 CDUM, S, TAU DOUBLE PRECISION FOO, SAFMAX, SAFMIN, SMLNUM, ULP INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN, $ KNT, KROW, KWTOP, LTOP, LWK1, LWK2, LWK3, @@ -304,7 +304,7 @@ SUBROUTINE ZLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, * .. External Subroutines .. EXTERNAL ZCOPY, ZGEHRD, ZGEMM, ZLACPY, ZLAHQR, $ ZLAQR4, - $ ZLARF, ZLARFG, ZLASET, ZTREXC, ZUNMHR + $ ZLARF1F, ZLARFG, ZLASET, ZTREXC, ZUNMHR * .. * .. Intrinsic Functions .. INTRINSIC ABS, DBLE, DCMPLX, DCONJG, DIMAG, INT, MAX, MIN @@ -490,19 +490,17 @@ SUBROUTINE ZLAQR3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, DO 50 I = 1, NS WORK( I ) = DCONJG( WORK( I ) ) 50 CONTINUE - BETA = WORK( 1 ) - CALL ZLARFG( NS, BETA, WORK( 2 ), 1, TAU ) - WORK( 1 ) = ONE + CALL ZLARFG( NS, WORK( 1 ), WORK( 2 ), 1, TAU ) * CALL ZLASET( 'L', JW-2, JW-2, ZERO, ZERO, T( 3, 1 ), $ LDT ) * - CALL ZLARF( 'L', NS, JW, WORK, 1, DCONJG( TAU ), T, LDT, - $ WORK( JW+1 ) ) - CALL ZLARF( 'R', NS, NS, WORK, 1, TAU, T, LDT, - $ WORK( JW+1 ) ) - CALL ZLARF( 'R', JW, NS, WORK, 1, TAU, V, LDV, - $ WORK( JW+1 ) ) + CALL ZLARF1F( 'L', NS, JW, WORK, 1, CONJG( TAU ), T, LDT, + $ WORK( JW+1 ) ) + CALL ZLARF1F( 'R', NS, NS, WORK, 1, TAU, T, LDT, + $ WORK( JW+1 ) ) + CALL ZLARF1F( 'R', JW, NS, WORK, 1, TAU, V, LDV, + $ WORK( JW+1 ) ) * CALL ZGEHRD( JW, 1, NS, T, LDT, WORK, WORK( JW+1 ), $ LWORK-JW, INFO ) diff --git a/SRC/zlarf1f.f b/SRC/zlarf1f.f new file mode 100644 index 000000000..8203d4f6c --- /dev/null +++ b/SRC/zlarf1f.f @@ -0,0 +1,303 @@ +*> \brief \b ZLARF1F applies an elementary reflector to a general rectangular +* matrix assuming v(1) = 1. +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARF1F + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARF1F( 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 +*> +*> ZLARF1F applies a complex 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**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 +*> +*> \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. V(1) is not referenced or modified. +*> \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 +* To take advantage of the fact that v(1) = 1, we do the following +* v = [ 1 v_2 ]**T +* If SIDE='L' +* |-----| +* | C_1 | +* C =| C_2 | +* |-----| +* C_1\in\mathbb{C}^{1\times n}, C_2\in\mathbb{C}^{m-1\times n} +* So we compute: +* C = HC = (I - \tau vv**T)C +* = C - \tau vv**T C +* w = C**T v = [ C_1**T C_2**T ] [ 1 v_2 ]**T +* = C_1**T + C_2**T v ( ZGEMM then ZAXPYC-like ) +* C = C - \tau vv**T C +* = C - \tau vw**T +* Giving us C_1 = C_1 - \tau w**T ( ZAXPYC-like ) +* and +* C_2 = C_2 - \tau v_2w**T ( ZGERC ) +* If SIDE='R' +* +* C = [ C_1 C_2 ] +* C_1\in\mathbb{C}^{m\times 1}, C_2\in\mathbb{C}^{m\times n-1} +* So we compute: +* C = CH = C(I - \tau vv**T) +* = C - \tau Cvv**T +* +* w = Cv = [ C_1 C_2 ] [ 1 v_2 ]**T +* = C_1 + C_2v_2 ( ZGEMM then ZAXPYC-like ) +* C = C - \tau Cvv**T +* = C - \tau wv**T +* Giving us C_1 = C_1 - \tau w ( ZAXPYC-like ) +* and +* C_2 = C_2 - \tau wv_2**T ( ZGERC ) +* +* Authors: +* ======== +* +*> \author Univ. of Tennessee +*> \author Univ. of California Berkeley +*> \author Univ. of Colorado Denver +*> \author NAG Ltd. +* +*> \ingroup larf +* +* ===================================================================== + SUBROUTINE ZLARF1F( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- 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 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, J +* .. +* .. External Subroutines .. + EXTERNAL ZGEMV, ZGERC, ZSCAL +* .. Intrinsic Functions .. + INTRINSIC DCONJG +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAZLR, ILAZLC + EXTERNAL LSAME, ILAZLR, ILAZLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + LASTV = 1 + 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. +! Since we are assuming that V(1) = 1, and it is not stored, so we +! shouldn't access it. + DO WHILE( LASTV.GT.1 .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 + IF( LASTC.EQ.0 ) THEN + RETURN + END IF + IF( APPLYLEFT ) THEN +* +* Form H * C +* + ! Check if m = 1. This means v = 1, So we just need to compute + ! C := HC = (1-\tau)C. + IF( LASTV.EQ.1 ) THEN + CALL ZSCAL(LASTC, ONE - TAU, C, LDC) + ELSE +* +* w(1:lastc,1) := C(1:lastv,1:lastc)**H * v(1:lastv,1) +* + ! (I - tvv**H)C = C - tvv**H C + ! First compute w**H = v**H c -> w = C**H v + ! C = [ C_1 C_2 ]**T, v = [1 v_2]**T + ! w = C_1**H + C_2**Hv_2 + ! w = C_2**Hv_2 + CALL ZGEMV( 'Conjugate transpose', LASTV - 1, + $ LASTC, ONE, C( 1+1, 1 ), LDC, V( 1 + INCV ), + $ INCV, ZERO, WORK, 1 ) +* +* w(1:lastc,1) += v(1,1) * C(1,1:lastc)**H +* + DO I = 1, LASTC + WORK( I ) = WORK( I ) + DCONJG( C( 1, I ) ) + END DO +* +* C(1:lastv,1:lastc) := C(...) - tau * v(1:lastv,1) * w(1:lastc,1)**H +* + ! C(1, 1:lastc) := C(...) - tau * v(1,1) * w(1:lastc,1)**H + ! = C(...) - tau * Conj(w(1:lastc,1)) + ! This is essentially a zaxpyc + DO I = 1, LASTC + C( 1, I ) = C( 1, I ) - TAU * DCONJG( WORK( I ) ) + END DO +* +* C(2:lastv,1:lastc) += - tau * v(2:lastv,1) * w(1:lastc,1)**H +* + CALL ZGERC( LASTV - 1, LASTC, -TAU, V( 1 + INCV ), + $ INCV, WORK, 1, C( 1+1, 1 ), LDC ) + END IF + ELSE +* +* Form C * H +* + ! Check if n = 1. This means v = 1, so we just need to compute + ! C := CH = C(1-\tau). + IF( LASTV.EQ.1 ) THEN + CALL ZSCAL(LASTC, ONE - TAU, C, 1) + ELSE +* +* w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1) +* + ! w(1:lastc,1) := C(1:lastc,2:lastv) * v(2:lastv,1) + CALL ZGEMV( 'No transpose', LASTC, LASTV-1, ONE, + $ C(1,1+1), LDC, V(1+INCV), INCV, ZERO, WORK, 1 ) + ! w(1:lastc,1) += C(1:lastc,1) v(1,1) = C(1:lastc,1) + CALL ZAXPY(LASTC, ONE, C, 1, WORK, 1) +* +* C(1:lastc,1:lastv) := C(...) - tau * w(1:lastc,1) * v(1:lastv,1)**T +* + ! C(1:lastc,1) := C(...) - tau * w(1:lastc,1) * v(1,1)**T + ! = C(...) - tau * w(1:lastc,1) + CALL ZAXPY(LASTC, -TAU, WORK, 1, C, 1) + ! C(1:lastc,2:lastv) := C(...) - tau * w(1:lastc,1) * v(2:lastv)**T + CALL ZGERC( LASTC, LASTV-1, -TAU, WORK, 1, V(1+INCV), + $ INCV, C(1,1+1), LDC ) + END IF + END IF + RETURN +* +* End of ZLARF1F +* + END diff --git a/SRC/zlarf1l.f b/SRC/zlarf1l.f new file mode 100644 index 000000000..5bb6eac15 --- /dev/null +++ b/SRC/zlarf1l.f @@ -0,0 +1,268 @@ +*> \brief \b ZLARF1L applies an elementary reflector to a general rectangular +* matrix assuming v(lastv) = 1, where lastv is the last non-zero +* +* =========== DOCUMENTATION =========== +* +* Online html documentation available at +* http://www.netlib.org/lapack/explore-html/ +* +*> \htmlonly +*> Download ZLARF1L + dependencies +*> +*> [TGZ] +*> +*> [ZIP] +*> +*> [TXT] +*> \endhtmlonly +* +* Definition: +* =========== +* +* SUBROUTINE ZLARF1L( 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 +*> +*> ZLARF1L 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 real scalar and v is a real vector assuming v(lastv) = 1, +*> where lastv is the last non-zero element. +*> +*> 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*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. +* +*> \ingroup larf1f +* +* ===================================================================== + SUBROUTINE ZLARF1L( SIDE, M, N, V, INCV, TAU, C, LDC, WORK ) +* +* -- 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 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, J, LASTV, LASTC, FIRSTV +* .. +* .. External Subroutines .. + EXTERNAL ZGEMV, ZGERC, ZSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC DCONJG +* .. +* .. External Functions .. + LOGICAL LSAME + INTEGER ILAZLR, ILAZLC + EXTERNAL LSAME, ILAZLR, ILAZLC +* .. +* .. Executable Statements .. +* + APPLYLEFT = LSAME( SIDE, 'L' ) + FIRSTV = 1 + LASTC = 0 + IF( TAU.NE.ZERO ) THEN +! Set up variables for scanning V. LASTV begins pointing to the end +! of V up to V(1). + IF( APPLYLEFT ) THEN + LASTV = M + ELSE + LASTV = N + END IF + I = 1 +! Look for the last non-zero row in V. + DO WHILE( LASTV.GT.FIRSTV .AND. V( I ).EQ.ZERO ) + FIRSTV = FIRSTV + 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 + IF( LASTC.EQ.0 ) THEN + RETURN + END IF + IF( APPLYLEFT ) THEN +* +* Form H * C +* + IF( LASTV.EQ.FIRSTV ) THEN +* +* C(lastv,1:lastc) := ( 1 - tau ) * C(lastv,1:lastc) +* + CALL ZSCAL( LASTC, ONE - TAU, C( LASTV, 1 ), LDC ) + ELSE +* +* w(1:lastc,1) := C(firstv:lastv-1,1:lastc)**T * v(firstv:lastv-1,1) +* + CALL ZGEMV( 'Conjugate transpose', LASTV - FIRSTV, LASTC, + $ ONE, C( FIRSTV, 1 ), LDC, V( I ), INCV, ZERO, + $ WORK, 1 ) +* +* w(1:lastc,1) += C(lastv,1:lastc)**H * v(lastv,1) +* + DO J = 1, LASTC + WORK( J ) = WORK( J ) + CONJG( C( LASTV, J ) ) + END DO +* +* C(lastv,1:lastc) += - tau * v(lastv,1) * w(1:lastc,1)**H +* + DO J = 1, LASTC + C( LASTV, J ) = C( LASTV, J ) + $ - TAU * CONJG( WORK( J ) ) + END DO +* +* C(firstv:lastv-1,1:lastc) += - tau * v(firstv:lastv-1,1) * w(1:lastc,1)**H +* + CALL ZGERC( LASTV - FIRSTV, LASTC, -TAU, V( I ), INCV, + $ WORK, 1, C( FIRSTV, 1 ), LDC) + END IF + ELSE +* +* Form C * H +* + IF( LASTV.EQ.FIRSTV ) THEN +* +* C(1:lastc,lastv) := ( 1 - tau ) * C(1:lastc,lastv) +* + CALL ZSCAL( LASTC, ONE - TAU, C( 1, LASTV ), 1 ) + ELSE +* +* w(1:lastc,1) := C(1:lastc,firstv:lastv-1) * v(firstv:lastv-1,1) +* + CALL ZGEMV( 'No transpose', LASTC, LASTV - FIRSTV, ONE, + $ C( 1, FIRSTV ), LDC, V( I ), INCV, ZERO, + $ WORK, 1 ) +* +* w(1:lastc,1) += C(1:lastc,lastv) * v(lastv,1) +* + CALL ZAXPY( LASTC, ONE, C( 1, LASTV ), 1, WORK, 1 ) +* +* C(1:lastc,lastv) += - tau * v(lastv,1) * w(1:lastc,1) +* + CALL ZAXPY( LASTC, -TAU, WORK, 1, C( 1, LASTV ), 1 ) +* +* C(1:lastc,firstv:lastv-1) += - tau * w(1:lastc,1) * v(firstv:lastv-1)**H +* + CALL ZGERC( LASTC, LASTV - FIRSTV, -TAU, WORK, 1, V( I ), + $ INCV, C( 1, FIRSTV ), LDC ) + END IF + END IF + RETURN +* +* End of ZLARF1L +* + END diff --git a/SRC/zunbdb.f b/SRC/zunbdb.f index f05e46e6d..250ed59ad 100644 --- a/SRC/zunbdb.f +++ b/SRC/zunbdb.f @@ -316,7 +316,7 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, DOUBLE PRECISION Z1, Z2, Z3, Z4 * .. * .. External Subroutines .. - EXTERNAL ZAXPY, ZLARF, ZLARFGP, ZSCAL, + EXTERNAL ZAXPY, ZLARF1F, ZLARFGP, ZSCAL, $ XERBLA EXTERNAL ZLACGV * @@ -427,7 +427,6 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, ELSE IF ( P .EQ. I ) THEN CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I), 1, TAUP1(I) ) END IF - X11(I,I) = ONE IF ( M-P .GT. I ) THEN CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, $ TAUP2(I) ) @@ -435,19 +434,20 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I), 1, $ TAUP2(I) ) END IF - X21(I,I) = ONE * IF ( Q .GT. I ) THEN - CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, - $ DCONJG(TAUP1(I)), X11(I,I+1), LDX11, WORK ) - CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, - $ DCONJG(TAUP2(I)), X21(I,I+1), LDX21, WORK ) + CALL ZLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, + $ CONJG(TAUP1(I)), X11(I,I+1), LDX11, + $ WORK ) + CALL ZLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, + $ CONJG(TAUP2(I)), X21(I,I+1), LDX21, + $ WORK ) END IF IF ( M-Q+1 .GT. I ) THEN - CALL ZLARF( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, - $ DCONJG(TAUP1(I)), X12(I,I), LDX12, WORK ) - CALL ZLARF( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, - $ DCONJG(TAUP2(I)), X22(I,I), LDX22, WORK ) + CALL ZLARF1F( 'L', P-I+1, M-Q-I+1, X11(I,I), 1, + $ CONJG(TAUP1(I)), X12(I,I), LDX12, WORK ) + CALL ZLARF1F( 'L', M-P-I+1, M-Q-I+1, X21(I,I), 1, + $ CONJG(TAUP2(I)), X22(I,I), LDX22, WORK ) END IF * IF( I .LT. Q ) THEN @@ -477,7 +477,6 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL ZLARFGP( Q-I, X11(I,I+1), X11(I,I+2), LDX11, $ TAUQ1(I) ) END IF - X11(I,I+1) = ONE END IF IF ( M-Q+1 .GT. I ) THEN CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) @@ -489,24 +488,23 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, $ TAUQ2(I) ) END IF END IF - X12(I,I) = ONE * IF( I .LT. Q ) THEN - CALL ZLARF( 'R', P-I, Q-I, X11(I,I+1), LDX11, - $ TAUQ1(I), - $ X11(I+1,I+1), LDX11, WORK ) - CALL ZLARF( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, - $ TAUQ1(I), - $ X21(I+1,I+1), LDX21, WORK ) + CALL ZLARF1F( 'R', P-I, Q-I, X11(I,I+1), LDX11, + $ TAUQ1(I), + $ X11(I+1,I+1), LDX11, WORK ) + CALL ZLARF1F( 'R', M-P-I, Q-I, X11(I,I+1), LDX11, + $ TAUQ1(I), + $ X21(I+1,I+1), LDX21, WORK ) END IF IF ( P .GT. I ) THEN - CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, - $ TAUQ2(I), - $ X12(I+1,I), LDX12, WORK ) + CALL ZLARF1F( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, + $ TAUQ2(I), + $ X12(I+1,I), LDX12, WORK ) END IF IF ( M-P .GT. I ) THEN - CALL ZLARF( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, - $ TAUQ2(I), X22(I+1,I), LDX22, WORK ) + CALL ZLARF1F( 'R', M-P-I, M-Q-I+1, X12(I,I), LDX12, + $ TAUQ2(I), X22(I+1,I), LDX22, WORK ) END IF * IF( I .LT. Q ) @@ -529,16 +527,15 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I,I+1), LDX12, $ TAUQ2(I) ) END IF - X12(I,I) = ONE * IF ( P .GT. I ) THEN - CALL ZLARF( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, - $ TAUQ2(I), - $ X12(I+1,I), LDX12, WORK ) + CALL ZLARF1F( 'R', P-I, M-Q-I+1, X12(I,I), LDX12, + $ TAUQ2(I), + $ X12(I+1,I), LDX12, WORK ) END IF IF( M-P-Q .GE. 1 ) - $ CALL ZLARF( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, - $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) + $ CALL ZLARF1F( 'R', M-P-Q, M-Q-I+1, X12(I,I), LDX12, + $ TAUQ2(I), X22(Q+1,I), LDX22, WORK ) * CALL ZLACGV( M-Q-I+1, X12(I,I), LDX12 ) * @@ -553,9 +550,9 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) CALL ZLARFGP( M-P-Q-I+1, X22(Q+I,P+I), X22(Q+I,P+I+1), $ LDX22, TAUQ2(P+I) ) - X22(Q+I,P+I) = ONE - CALL ZLARF( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), LDX22, - $ TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, WORK ) + CALL ZLARF1F( 'R', M-P-Q-I, M-P-Q-I+1, X22(Q+I,P+I), + $ LDX22, TAUQ2(P+I), X22(Q+I+1,P+I), LDX22, + $ WORK ) * CALL ZLACGV( M-P-Q-I+1, X22(Q+I,P+I), LDX22 ) * @@ -595,7 +592,6 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, * CALL ZLARFGP( P-I+1, X11(I,I), X11(I,I+1), LDX11, $ TAUP1(I) ) - X11(I,I) = ONE IF ( I .EQ. M-P ) THEN CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I), LDX21, $ TAUP2(I) ) @@ -603,17 +599,16 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I,I+1), LDX21, $ TAUP2(I) ) END IF - X21(I,I) = ONE -* - CALL ZLARF( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), - $ X11(I+1,I), LDX11, WORK ) - CALL ZLARF( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, - $ TAUP1(I), - $ X12(I,I), LDX12, WORK ) - CALL ZLARF( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, TAUP2(I), - $ X21(I+1,I), LDX21, WORK ) - CALL ZLARF( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, - $ TAUP2(I), X22(I,I), LDX22, WORK ) +* + CALL ZLARF1F( 'R', Q-I, P-I+1, X11(I,I), LDX11, TAUP1(I), + $ X11(I+1,I), LDX11, WORK ) + CALL ZLARF1F( 'R', M-Q-I+1, P-I+1, X11(I,I), LDX11, + $ TAUP1(I), + $ X12(I,I), LDX12, WORK ) + CALL ZLARF1F( 'R', Q-I, M-P-I+1, X21(I,I), LDX21, + $ TAUP2(I), X21(I+1,I), LDX21, WORK ) + CALL ZLARF1F( 'R', M-Q-I+1, M-P-I+1, X21(I,I), LDX21, + $ TAUP2(I), X22(I,I), LDX22, WORK ) * CALL ZLACGV( P-I+1, X11(I,I), LDX11 ) CALL ZLACGV( M-P-I+1, X21(I,I), LDX21 ) @@ -639,23 +634,26 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, IF( I .LT. Q ) THEN CALL ZLARFGP( Q-I, X11(I+1,I), X11(I+2,I), 1, $ TAUQ1(I) ) - X11(I+1,I) = ONE END IF CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, $ TAUQ2(I) ) - X12(I,I) = ONE * IF( I .LT. Q ) THEN - CALL ZLARF( 'L', Q-I, P-I, X11(I+1,I), 1, - $ DCONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, WORK ) - CALL ZLARF( 'L', Q-I, M-P-I, X11(I+1,I), 1, - $ DCONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, WORK ) + CALL ZLARF1F( 'L', Q-I, P-I, X11(I+1,I), 1, + $ CONJG(TAUQ1(I)), X11(I+1,I+1), LDX11, + $ WORK ) + CALL ZLARF1F( 'L', Q-I, M-P-I, X11(I+1,I), 1, + $ CONJG(TAUQ1(I)), X21(I+1,I+1), LDX21, + $ WORK ) END IF - CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, - $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) + CALL ZLARF1F( 'L', M-Q-I+1, P-I, X12(I,I), 1, + $ CONJG(TAUQ2(I)), + $ X12(I,I+1), LDX12, WORK ) + IF ( M-P .GT. I ) THEN - CALL ZLARF( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, - $ DCONJG(TAUQ2(I)), X22(I,I+1), LDX22, WORK ) + CALL ZLARF1F( 'L', M-Q-I+1, M-P-I, X12(I,I), 1, + $ CONJG(TAUQ2(I)), X22(I,I+1), LDX22, + $ WORK ) END IF * END DO @@ -668,15 +666,16 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, $ 1 ) CALL ZLARFGP( M-Q-I+1, X12(I,I), X12(I+1,I), 1, $ TAUQ2(I) ) - X12(I,I) = ONE * IF ( P .GT. I ) THEN - CALL ZLARF( 'L', M-Q-I+1, P-I, X12(I,I), 1, - $ DCONJG(TAUQ2(I)), X12(I,I+1), LDX12, WORK ) + CALL ZLARF1F( 'L', M-Q-I+1, P-I, X12(I,I), 1, + $ CONJG(TAUQ2(I)), X12(I,I+1), LDX12, + $ WORK ) END IF IF( M-P-Q .GE. 1 ) - $ CALL ZLARF( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, - $ DCONJG(TAUQ2(I)), X22(I,Q+1), LDX22, WORK ) + $ CALL ZLARF1F( 'L', M-Q-I+1, M-P-Q, X12(I,I), 1, + $ CONJG(TAUQ2(I)), X22(I,Q+1), LDX22, + $ WORK ) * END DO * @@ -688,12 +687,10 @@ SUBROUTINE ZUNBDB( TRANS, SIGNS, M, P, Q, X11, LDX11, X12, $ X22(P+I,Q+I), 1 ) CALL ZLARFGP( M-P-Q-I+1, X22(P+I,Q+I), X22(P+I+1,Q+I), 1, $ TAUQ2(P+I) ) - X22(P+I,Q+I) = ONE -* IF ( M-P-Q .NE. I ) THEN - CALL ZLARF( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), 1, - $ DCONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), LDX22, - $ WORK ) + CALL ZLARF1F( 'L', M-P-Q-I+1, M-P-Q-I, X22(P+I,Q+I), + $ 1, CONJG(TAUQ2(P+I)), X22(P+I,Q+I+1), + $ LDX22, WORK ) END IF * END DO diff --git a/SRC/zunbdb1.f b/SRC/zunbdb1.f index b96c49993..af0eeaea7 100644 --- a/SRC/zunbdb1.f +++ b/SRC/zunbdb1.f @@ -228,7 +228,7 @@ SUBROUTINE ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFGP, ZUNBDB5, ZDROT, + EXTERNAL ZLARF1F, ZLARFGP, ZUNBDB5, ZDROT, $ XERBLA EXTERNAL ZLACGV * .. @@ -288,13 +288,13 @@ SUBROUTINE ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( DBLE( X21(I,I) ), DBLE( X11(I,I) ) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - X11(I,I) = ONE - X21(I,I) = ONE - CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)), - $ X11(I,I+1), LDX11, WORK(ILARF) ) - CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, - $ DCONJG(TAUP2(I)), - $ X21(I,I+1), LDX21, WORK(ILARF) ) + C = COS( THETA(I) ) + S = SIN( THETA(I) ) + CALL ZLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, CONJG(TAUP1(I)), + $ X11(I,I+1), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, + $ CONJG(TAUP2(I)), X21(I,I+1), LDX21, + $ WORK(ILARF) ) * IF( I .LT. Q ) THEN CALL ZDROT( Q-I, X11(I,I+1), LDX11, X21(I,I+1), LDX21, C, @@ -303,14 +303,14 @@ SUBROUTINE ZUNBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA, CALL ZLARFGP( Q-I, X21(I,I+1), X21(I,I+2), LDX21, $ TAUQ1(I) ) S = DBLE( X21(I,I+1) ) - X21(I,I+1) = ONE - CALL ZLARF( 'R', P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), - $ X11(I+1,I+1), LDX11, WORK(ILARF) ) - CALL ZLARF( 'R', M-P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), - $ X21(I+1,I+1), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', P-I, Q-I, X21(I,I+1), LDX21, TAUQ1(I), + $ X11(I+1,I+1), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'R', M-P-I, Q-I, X21(I,I+1), LDX21, + $ TAUQ1(I), X21(I+1,I+1), LDX21, + $ WORK(ILARF) ) CALL ZLACGV( Q-I, X21(I,I+1), LDX21 ) C = SQRT( DZNRM2( P-I, X11(I+1,I+1), 1 )**2 - $ + DZNRM2( M-P-I, X21(I+1,I+1), 1 )**2 ) + $ + DZNRM2( M-P-I, X21(I+1,I+1), 1 )**2 ) PHI(I) = ATAN2( S, C ) CALL ZUNBDB5( P-I, M-P-I, Q-I-1, X11(I+1,I+1), 1, $ X21(I+1,I+1), 1, X11(I+1,I+2), LDX11, diff --git a/SRC/zunbdb2.f b/SRC/zunbdb2.f index 245391982..3d0dc4a59 100644 --- a/SRC/zunbdb2.f +++ b/SRC/zunbdb2.f @@ -227,7 +227,7 @@ SUBROUTINE ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFGP, ZUNBDB5, ZDROT, ZSCAL, + EXTERNAL ZLARF1F, ZLARFGP, ZUNBDB5, ZDROT, ZSCAL, $ ZLACGV, $ XERBLA * .. @@ -289,11 +289,10 @@ SUBROUTINE ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, CALL ZLACGV( Q-I+1, X11(I,I), LDX11 ) CALL ZLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) ) C = DBLE( X11(I,I) ) - X11(I,I) = ONE - CALL ZLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), - $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL ZLARF( 'R', M-P-I+1, Q-I+1, X11(I,I), LDX11, TAUQ1(I), - $ X21(I,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + $ X11(I+1,I), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'R', M-P-I+1, Q-I+1, X11(I,I), LDX11, + $ TAUQ1(I), X21(I,I), LDX21, WORK(ILARF) ) CALL ZLACGV( Q-I+1, X11(I,I), LDX11 ) S = SQRT( DZNRM2( P-I, X11(I+1,I), 1 )**2 $ + DZNRM2( M-P-I+1, X21(I,I), 1 )**2 ) @@ -309,15 +308,13 @@ SUBROUTINE ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI(I) = ATAN2( DBLE( X11(I+1,I) ), DBLE( X21(I,I) ) ) C = COS( PHI(I) ) S = SIN( PHI(I) ) - X11(I+1,I) = ONE - CALL ZLARF( 'L', P-I, Q-I, X11(I+1,I), 1, - $ DCONJG(TAUP1(I)), - $ X11(I+1,I+1), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'L', P-I, Q-I, X11(I+1,I), 1, + $ CONJG(TAUP1(I)), + $ X11(I+1,I+1), LDX11, WORK(ILARF) ) END IF - X21(I,I) = ONE - CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, - $ DCONJG(TAUP2(I)), - $ X21(I,I+1), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, + $ CONJG(TAUP2(I)), X21(I,I+1), LDX21, + $ WORK(ILARF) ) * END DO * @@ -325,10 +322,9 @@ SUBROUTINE ZUNBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA, * DO I = P + 1, Q CALL ZLARFGP( M-P-I+1, X21(I,I), X21(I+1,I), 1, TAUP2(I) ) - X21(I,I) = ONE - CALL ZLARF( 'L', M-P-I+1, Q-I, X21(I,I), 1, - $ DCONJG(TAUP2(I)), - $ X21(I,I+1), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P-I+1, Q-I, X21(I,I), 1, + $ CONJG(TAUP2(I)), X21(I,I+1), LDX21, + $ WORK(ILARF) ) END DO * RETURN diff --git a/SRC/zunbdb3.f b/SRC/zunbdb3.f index 67b3eeedc..d7aa8fadb 100644 --- a/SRC/zunbdb3.f +++ b/SRC/zunbdb3.f @@ -226,7 +226,7 @@ SUBROUTINE ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFGP, ZUNBDB5, ZDROT, ZLACGV, + EXTERNAL ZLARF1F, ZLARFGP, ZUNBDB5, ZDROT, ZLACGV, $ XERBLA * .. * .. External Functions .. @@ -285,14 +285,12 @@ SUBROUTINE ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, $ S ) END IF * - CALL ZLACGV( Q-I+1, X21(I,I), LDX21 ) CALL ZLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) ) S = DBLE( X21(I,I) ) - X21(I,I) = ONE - CALL ZLARF( 'R', P-I+1, Q-I+1, X21(I,I), LDX21, TAUQ1(I), - $ X11(I,I), LDX11, WORK(ILARF) ) - CALL ZLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), - $ X21(I+1,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', P-I+1, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + $ X11(I,I), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + $ X21(I+1,I), LDX21, WORK(ILARF) ) CALL ZLACGV( Q-I+1, X21(I,I), LDX21 ) C = SQRT( DZNRM2( P-I+1, X11(I,I), 1 )**2 $ + DZNRM2( M-P-I, X21(I+1,I), 1 )**2 ) @@ -308,24 +306,20 @@ SUBROUTINE ZUNBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI(I) = ATAN2( DBLE( X21(I+1,I) ), DBLE( X11(I,I) ) ) C = COS( PHI(I) ) S = SIN( PHI(I) ) - X21(I+1,I) = ONE - CALL ZLARF( 'L', M-P-I, Q-I, X21(I+1,I), 1, - $ DCONJG(TAUP2(I)), X21(I+1,I+1), LDX21, - $ WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P-I, Q-I, X21(I+1,I), 1, + $ CONJG(TAUP2(I)), + $ X21(I+1,I+1), LDX21, WORK(ILARF) ) END IF - X11(I,I) = ONE - CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)), - $ X11(I,I+1), LDX11, WORK(ILARF) ) -* + CALL ZLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, CONJG(TAUP1(I)), + $ X11(I,I+1), LDX11, WORK(ILARF) ) END DO * * Reduce the bottom-right portion of X11 to the identity matrix * DO I = M-P + 1, Q CALL ZLARFGP( P-I+1, X11(I,I), X11(I+1,I), 1, TAUP1(I) ) - X11(I,I) = ONE - CALL ZLARF( 'L', P-I+1, Q-I, X11(I,I), 1, DCONJG(TAUP1(I)), - $ X11(I,I+1), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'L', P-I+1, Q-I, X11(I,I), 1, CONJG(TAUP1(I)), + $ X11(I,I+1), LDX11, WORK(ILARF) ) END DO * RETURN diff --git a/SRC/zunbdb4.f b/SRC/zunbdb4.f index a242d956d..90d70168f 100644 --- a/SRC/zunbdb4.f +++ b/SRC/zunbdb4.f @@ -239,7 +239,7 @@ SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, LOGICAL LQUERY * .. * .. External Subroutines .. - EXTERNAL ZLARF, ZLARFGP, ZUNBDB5, ZDROT, ZSCAL, + EXTERNAL ZLARF1F, ZLARFGP, ZUNBDB5, ZDROT, ZSCAL, $ ZLACGV, $ XERBLA * .. @@ -309,14 +309,12 @@ SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( DBLE( PHANTOM(1) ), DBLE( PHANTOM(P+1) ) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - PHANTOM(1) = ONE - PHANTOM(P+1) = ONE - CALL ZLARF( 'L', P, Q, PHANTOM(1), 1, DCONJG(TAUP1(1)), - $ X11, - $ LDX11, WORK(ILARF) ) - CALL ZLARF( 'L', M-P, Q, PHANTOM(P+1), 1, - $ DCONJG(TAUP2(1)), - $ X21, LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'L', P, Q, PHANTOM(1), 1, CONJG(TAUP1(1)), + $ X11, + $ LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P, Q, PHANTOM(P+1), 1, + $ CONJG(TAUP2(1)), + $ X21, LDX21, WORK(ILARF) ) ELSE CALL ZUNBDB5( P-I+1, M-P-I+1, Q-I+1, X11(I,I-1), 1, $ X21(I,I-1), 1, X11(I,I), LDX11, X21(I,I), @@ -329,23 +327,22 @@ SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, THETA(I) = ATAN2( DBLE( X11(I,I-1) ), DBLE( X21(I,I-1) ) ) C = COS( THETA(I) ) S = SIN( THETA(I) ) - X11(I,I-1) = ONE - X21(I,I-1) = ONE - CALL ZLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, - $ DCONJG(TAUP1(I)), X11(I,I), LDX11, WORK(ILARF) ) - CALL ZLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, - $ DCONJG(TAUP2(I)), X21(I,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'L', P-I+1, Q-I+1, X11(I,I-1), 1, + $ CONJG(TAUP1(I)), X11(I,I), LDX11, + $ WORK(ILARF) ) + CALL ZLARF1F( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1, + $ CONJG(TAUP2(I)), X21(I,I), LDX21, + $ WORK(ILARF) ) END IF * CALL ZDROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C ) CALL ZLACGV( Q-I+1, X21(I,I), LDX21 ) CALL ZLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) ) C = DBLE( X21(I,I) ) - X21(I,I) = ONE - CALL ZLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), - $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL ZLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), - $ X21(I+1,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + $ X11(I+1,I), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I), + $ X21(I+1,I), LDX21, WORK(ILARF) ) CALL ZLACGV( Q-I+1, X21(I,I), LDX21 ) IF( I .LT. M-Q ) THEN S = SQRT( DZNRM2( P-I, X11(I+1,I), 1 )**2 @@ -360,11 +357,10 @@ SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, DO I = M - Q + 1, P CALL ZLACGV( Q-I+1, X11(I,I), LDX11 ) CALL ZLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) ) - X11(I,I) = ONE - CALL ZLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), - $ X11(I+1,I), LDX11, WORK(ILARF) ) - CALL ZLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I), - $ X21(M-Q+1,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + $ X11(I+1,I), LDX11, WORK(ILARF) ) + CALL ZLARF1F( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I), + $ X21(M-Q+1,I), LDX21, WORK(ILARF) ) CALL ZLACGV( Q-I+1, X11(I,I), LDX11 ) END DO * @@ -375,10 +371,9 @@ SUBROUTINE ZUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, CALL ZLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), $ LDX21, $ TAUQ1(I) ) - X21(M-Q+I-P,I) = ONE - CALL ZLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, - $ TAUQ1(I), - $ X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) ) + CALL ZLARF1F( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, + $ TAUQ1(I), + $ X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) ) CALL ZLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 ) END DO * diff --git a/SRC/zung2l.f b/SRC/zung2l.f index 28854861b..814034e18 100644 --- a/SRC/zung2l.f +++ b/SRC/zung2l.f @@ -134,7 +134,7 @@ SUBROUTINE ZUNG2L( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, II, J, L * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZSCAL + EXTERNAL XERBLA, ZLARF1L, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -178,9 +178,9 @@ SUBROUTINE ZUNG2L( M, N, K, A, LDA, TAU, WORK, INFO ) * Apply H(i) to A(1:m-k+i,1:n-k+i) from the left * A( M-N+II, II ) = ONE - CALL ZLARF( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), - $ A, - $ LDA, WORK ) + CALL ZLARF1L( 'Left', M-N+II, II-1, A( 1, II ), 1, TAU( I ), + $ A, + $ LDA, WORK ) CALL ZSCAL( M-N+II-1, -TAU( I ), A( 1, II ), 1 ) A( M-N+II, II ) = ONE - TAU( I ) * diff --git a/SRC/zung2r.f b/SRC/zung2r.f index b73246b2b..80237cf31 100644 --- a/SRC/zung2r.f +++ b/SRC/zung2r.f @@ -134,7 +134,7 @@ SUBROUTINE ZUNG2R( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, J, L * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF, ZSCAL + EXTERNAL XERBLA, ZLARF1F, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC MAX @@ -177,9 +177,8 @@ SUBROUTINE ZUNG2R( M, N, K, A, LDA, TAU, WORK, INFO ) * Apply H(i) to A(i:m,i:n) from the left * IF( I.LT.N ) THEN - A( I, I ) = ONE - CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), - $ A( I, I+1 ), LDA, WORK ) + CALL ZLARF1F( 'Left', M-I+1, N-I, A( I, I ), 1, TAU( I ), + $ A( I, I+1 ), LDA, WORK ) END IF IF( I.LT.M ) $ CALL ZSCAL( M-I, -TAU( I ), A( I+1, I ), 1 ) diff --git a/SRC/zungl2.f b/SRC/zungl2.f index 83308c59b..4647c93f2 100644 --- a/SRC/zungl2.f +++ b/SRC/zungl2.f @@ -133,7 +133,7 @@ SUBROUTINE ZUNGL2( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, J, L * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF, ZSCAL + EXTERNAL XERBLA, ZLACGV, ZLARF1F, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -182,9 +182,9 @@ SUBROUTINE ZUNGL2( M, N, K, A, LDA, TAU, WORK, INFO ) IF( I.LT.N ) THEN CALL ZLACGV( N-I, A( I, I+1 ), LDA ) IF( I.LT.M ) THEN - A( I, I ) = ONE - CALL ZLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, - $ DCONJG( TAU( I ) ), A( I+1, I ), LDA, WORK ) + CALL ZLARF1F( 'Right', M-I, N-I+1, A( I, I ), LDA, + $ CONJG( TAU( I ) ), A( I+1, I ), LDA, + $ WORK ) END IF CALL ZSCAL( N-I, -TAU( I ), A( I, I+1 ), LDA ) CALL ZLACGV( N-I, A( I, I+1 ), LDA ) diff --git a/SRC/zungr2.f b/SRC/zungr2.f index 05c5fc74e..33f6be35f 100644 --- a/SRC/zungr2.f +++ b/SRC/zungr2.f @@ -134,7 +134,7 @@ SUBROUTINE ZUNGR2( M, N, K, A, LDA, TAU, WORK, INFO ) INTEGER I, II, J, L * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF, ZSCAL + EXTERNAL XERBLA, ZLACGV, ZLARF1L, ZSCAL * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -182,9 +182,8 @@ SUBROUTINE ZUNGR2( M, N, K, A, LDA, TAU, WORK, INFO ) * Apply H(i)**H to A(1:m-k+i,1:n-k+i) from the right * CALL ZLACGV( N-M+II-1, A( II, 1 ), LDA ) - A( II, N-M+II ) = ONE - CALL ZLARF( 'Right', II-1, N-M+II, A( II, 1 ), LDA, - $ DCONJG( TAU( I ) ), A, LDA, WORK ) + CALL ZLARF1L( 'Right', II-1, N-M+II, A( II, 1 ), LDA, + $ CONJG( TAU( I ) ), A, LDA, WORK ) CALL ZSCAL( N-M+II-1, -TAU( I ), A( II, 1 ), LDA ) CALL ZLACGV( N-M+II-1, A( II, 1 ), LDA ) A( II, N-M+II ) = ONE - DCONJG( TAU( I ) ) diff --git a/SRC/zunm2l.f b/SRC/zunm2l.f index 0e0ed1c06..08dac348c 100644 --- a/SRC/zunm2l.f +++ b/SRC/zunm2l.f @@ -178,14 +178,14 @@ SUBROUTINE ZUNM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, MI, NI, NQ - COMPLEX*16 AII, TAUI + COMPLEX*16 TAUI * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF + EXTERNAL XERBLA, ZLARF1L * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -266,10 +266,8 @@ SUBROUTINE ZUNM2L( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, ELSE TAUI = DCONJG( TAU( I ) ) END IF - AII = A( NQ-K+I, I ) - A( NQ-K+I, I ) = ONE - CALL ZLARF( SIDE, MI, NI, A( 1, I ), 1, TAUI, C, LDC, WORK ) - A( NQ-K+I, I ) = AII + CALL ZLARF1L( SIDE, MI, NI, A( 1, I ), 1, TAUI, C, LDC, + $ WORK ) 10 CONTINUE RETURN * diff --git a/SRC/zunm2r.f b/SRC/zunm2r.f index 6d6c802a6..ce7d53557 100644 --- a/SRC/zunm2r.f +++ b/SRC/zunm2r.f @@ -178,14 +178,14 @@ SUBROUTINE ZUNM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ - COMPLEX*16 AII, TAUI + COMPLEX*16 TAUI * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF + EXTERNAL XERBLA, ZLARF1F * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -270,12 +270,9 @@ SUBROUTINE ZUNM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, ELSE TAUI = DCONJG( TAU( I ) ) END IF - AII = A( I, I ) - A( I, I ) = ONE - CALL ZLARF( SIDE, MI, NI, A( I, I ), 1, TAUI, C( IC, JC ), + CALL ZLARF1F( SIDE, MI, NI, A( I, I ), 1, TAUI, C( IC, JC ), $ LDC, $ WORK ) - A( I, I ) = AII 10 CONTINUE RETURN * diff --git a/SRC/zunml2.f b/SRC/zunml2.f index 00385dc61..dace8cce0 100644 --- a/SRC/zunml2.f +++ b/SRC/zunml2.f @@ -178,14 +178,14 @@ SUBROUTINE ZUNML2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ - COMPLEX*16 AII, TAUI + COMPLEX*16 TAUI * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF + EXTERNAL XERBLA, ZLACGV, ZLARF1F * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -272,11 +272,8 @@ SUBROUTINE ZUNML2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, END IF IF( I.LT.NQ ) $ CALL ZLACGV( NQ-I, A( I, I+1 ), LDA ) - AII = A( I, I ) - A( I, I ) = ONE - CALL ZLARF( SIDE, MI, NI, A( I, I ), LDA, TAUI, C( IC, JC ), - $ LDC, WORK ) - A( I, I ) = AII + CALL ZLARF1F( SIDE, MI, NI, A( I, I ), LDA, TAUI, C( IC, + $ JC ), LDC, WORK ) IF( I.LT.NQ ) $ CALL ZLACGV( NQ-I, A( I, I+1 ), LDA ) 10 CONTINUE diff --git a/SRC/zunmr2.f b/SRC/zunmr2.f index 654752217..15768d23e 100644 --- a/SRC/zunmr2.f +++ b/SRC/zunmr2.f @@ -178,14 +178,14 @@ SUBROUTINE ZUNMR2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, * .. Local Scalars .. LOGICAL LEFT, NOTRAN INTEGER I, I1, I2, I3, MI, NI, NQ - COMPLEX*16 AII, TAUI + COMPLEX*16 TAUI * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLACGV, ZLARF + EXTERNAL XERBLA, ZLACGV, ZLARF1L * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -267,11 +267,8 @@ SUBROUTINE ZUNMR2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, TAUI = TAU( I ) END IF CALL ZLACGV( NQ-K+I-1, A( I, 1 ), LDA ) - AII = A( I, NQ-K+I ) - A( I, NQ-K+I ) = ONE - CALL ZLARF( SIDE, MI, NI, A( I, 1 ), LDA, TAUI, C, LDC, - $ WORK ) - A( I, NQ-K+I ) = AII + CALL ZLARF1L( SIDE, MI, NI, A( I, 1 ), LDA, TAUI, C, LDC, + $ WORK ) CALL ZLACGV( NQ-K+I-1, A( I, 1 ), LDA ) 10 CONTINUE RETURN diff --git a/SRC/zupmtr.f b/SRC/zupmtr.f index acf922f6d..b3b8b8eb1 100644 --- a/SRC/zupmtr.f +++ b/SRC/zupmtr.f @@ -170,14 +170,14 @@ SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, * .. Local Scalars .. LOGICAL FORWRD, LEFT, NOTRAN, UPPER INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ - COMPLEX*16 AII, TAUI + COMPLEX*16 TAUI * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZLARF + EXTERNAL XERBLA, ZLARF1, ZLARF1F * .. * .. Intrinsic Functions .. INTRINSIC DCONJG, MAX @@ -266,11 +266,8 @@ SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, ELSE TAUI = DCONJG( TAU( I ) ) END IF - AII = AP( II ) - AP( II ) = ONE - CALL ZLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, LDC, - $ WORK ) - AP( II ) = AII + CALL ZLARF1L( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, + $ LDC, WORK ) * IF( FORWRD ) THEN II = II + I + 2 @@ -306,8 +303,6 @@ SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, END IF * DO 20 I = I1, I2, I3 - AII = AP( II ) - AP( II ) = ONE IF( LEFT ) THEN * * H(i) or H(i)**H is applied to C(i+1:m,1:n) @@ -329,9 +324,8 @@ SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, ELSE TAUI = DCONJG( TAU( I ) ) END IF - CALL ZLARF( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, JC ), - $ LDC, WORK ) - AP( II ) = AII + CALL ZLARF1F( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, + $ JC ), LDC, WORK ) * IF( FORWRD ) THEN II = II + NQ - I + 1