1       SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
  2       IMPLICIT NONE
  3 *
  4 *  -- LAPACK auxiliary routine (version 3.3.1) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *  -- April 2011                                                      --
  8 *
  9 *     .. Scalar Arguments ..
 10       CHARACTER          SIDE
 11       INTEGER            INCV, LDC, M, N
 12       DOUBLE PRECISION   TAU
 13 *     ..
 14 *     .. Array Arguments ..
 15       DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  DLARF applies a real elementary reflector H to a real m by n matrix
 22 *  C, from either the left or the right. H is represented in the form
 23 *
 24 *        H = I - tau * v * v**T
 25 *
 26 *  where tau is a real scalar and v is a real vector.
 27 *
 28 *  If tau = 0, then H is taken to be the unit matrix.
 29 *
 30 *  Arguments
 31 *  =========
 32 *
 33 *  SIDE    (input) CHARACTER*1
 34 *          = 'L': form  H * C
 35 *          = 'R': form  C * H
 36 *
 37 *  M       (input) INTEGER
 38 *          The number of rows of the matrix C.
 39 *
 40 *  N       (input) INTEGER
 41 *          The number of columns of the matrix C.
 42 *
 43 *  V       (input) DOUBLE PRECISION array, dimension
 44 *                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
 45 *                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
 46 *          The vector v in the representation of H. V is not used if
 47 *          TAU = 0.
 48 *
 49 *  INCV    (input) INTEGER
 50 *          The increment between elements of v. INCV <> 0.
 51 *
 52 *  TAU     (input) DOUBLE PRECISION
 53 *          The value tau in the representation of H.
 54 *
 55 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 56 *          On entry, the m by n matrix C.
 57 *          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
 58 *          or C * H if SIDE = 'R'.
 59 *
 60 *  LDC     (input) INTEGER
 61 *          The leading dimension of the array C. LDC >= max(1,M).
 62 *
 63 *  WORK    (workspace) DOUBLE PRECISION array, dimension
 64 *                         (N) if SIDE = 'L'
 65 *                      or (M) if SIDE = 'R'
 66 *
 67 *  =====================================================================
 68 *
 69 *     .. Parameters ..
 70       DOUBLE PRECISION   ONE, ZERO
 71       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 72 *     ..
 73 *     .. Local Scalars ..
 74       LOGICAL            APPLYLEFT
 75       INTEGER            I, LASTV, LASTC
 76 *     ..
 77 *     .. External Subroutines ..
 78       EXTERNAL           DGEMV, DGER
 79 *     ..
 80 *     .. External Functions ..
 81       LOGICAL            LSAME
 82       INTEGER            ILADLR, ILADLC
 83       EXTERNAL           LSAME, ILADLR, ILADLC
 84 *     ..
 85 *     .. Executable Statements ..
 86 *
 87       APPLYLEFT = LSAME( SIDE, 'L' )
 88       LASTV = 0
 89       LASTC = 0
 90       IF( TAU.NE.ZERO ) THEN
 91 !     Set up variables for scanning V.  LASTV begins pointing to the end
 92 !     of V.
 93          IF( APPLYLEFT ) THEN
 94             LASTV = M
 95          ELSE
 96             LASTV = N
 97          END IF
 98          IF( INCV.GT.0 ) THEN
 99             I = 1 + (LASTV-1* INCV
100          ELSE
101             I = 1
102          END IF
103 !     Look for the last non-zero row in V.
104          DO WHILE( LASTV.GT.0 .AND. V( I ).EQ.ZERO )
105             LASTV = LASTV - 1
106             I = I - INCV
107          END DO
108          IF( APPLYLEFT ) THEN
109 !     Scan for the last non-zero column in C(1:lastv,:).
110             LASTC = ILADLC(LASTV, N, C, LDC)
111          ELSE
112 !     Scan for the last non-zero row in C(:,1:lastv).
113             LASTC = ILADLR(M, LASTV, C, LDC)
114          END IF
115       END IF
116 !     Note that lastc.eq.0 renders the BLAS operations null; no special
117 !     case is needed at this level.
118       IF( APPLYLEFT ) THEN
119 *
120 *        Form  H * C
121 *
122          IF( LASTV.GT.0 ) THEN
123 *
124 *           w(1:lastc,1) := C(1:lastv,1:lastc)**T * v(1:lastv,1)
125 *
126             CALL DGEMV( 'Transpose', LASTV, LASTC, ONE, C, LDC, V, INCV,
127      $           ZERO, WORK, 1 )
128 *
129 *           C(1:lastv,1:lastc) := C(...) - v(1:lastv,1) * w(1:lastc,1)**T
130 *
131             CALL DGER( LASTV, LASTC, -TAU, V, INCV, WORK, 1, C, LDC )
132          END IF
133       ELSE
134 *
135 *        Form  C * H
136 *
137          IF( LASTV.GT.0 ) THEN
138 *
139 *           w(1:lastc,1) := C(1:lastc,1:lastv) * v(1:lastv,1)
140 *
141             CALL DGEMV( 'No transpose', LASTC, LASTV, ONE, C, LDC,
142      $           V, INCV, ZERO, WORK, 1 )
143 *
144 *           C(1:lastc,1:lastv) := C(...) - w(1:lastc,1) * v(1:lastv,1)**T
145 *
146             CALL DGER( LASTC, LASTV, -TAU, WORK, 1, V, INCV, C, LDC )
147          END IF
148       END IF
149       RETURN
150 *
151 *     End of DLARF
152 *
153       END