Merge pull request #1020 from jprhyne/orm2r

Implementation of dorm2r and dlarf1f

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

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 \

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

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

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

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
 *

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

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

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

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
 *

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.

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

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 ),