```001:       SUBROUTINE DLARF( SIDE, M, N, V, INCV, TAU, C, LDC, WORK )
002:       IMPLICIT NONE
003: *
004: *  -- LAPACK auxiliary routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          SIDE
011:       INTEGER            INCV, LDC, M, N
012:       DOUBLE PRECISION   TAU
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   C( LDC, * ), V( * ), WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DLARF applies a real elementary reflector H to a real m by n matrix
022: *  C, from either the left or the right. H is represented in the form
023: *
024: *        H = I - tau * v * v'
025: *
026: *  where tau is a real scalar and v is a real vector.
027: *
028: *  If tau = 0, then H is taken to be the unit matrix.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  SIDE    (input) CHARACTER*1
034: *          = 'L': form  H * C
035: *          = 'R': form  C * H
036: *
037: *  M       (input) INTEGER
038: *          The number of rows of the matrix C.
039: *
040: *  N       (input) INTEGER
041: *          The number of columns of the matrix C.
042: *
043: *  V       (input) DOUBLE PRECISION array, dimension
044: *                     (1 + (M-1)*abs(INCV)) if SIDE = 'L'
045: *                  or (1 + (N-1)*abs(INCV)) if SIDE = 'R'
046: *          The vector v in the representation of H. V is not used if
047: *          TAU = 0.
048: *
049: *  INCV    (input) INTEGER
050: *          The increment between elements of v. INCV <> 0.
051: *
052: *  TAU     (input) DOUBLE PRECISION
053: *          The value tau in the representation of H.
054: *
055: *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
056: *          On entry, the m by n matrix C.
057: *          On exit, C is overwritten by the matrix H * C if SIDE = 'L',
058: *          or C * H if SIDE = 'R'.
059: *
060: *  LDC     (input) INTEGER
061: *          The leading dimension of the array C. LDC >= max(1,M).
062: *
063: *  WORK    (workspace) DOUBLE PRECISION array, dimension
064: *                         (N) if SIDE = 'L'
065: *                      or (M) if SIDE = 'R'
066: *
067: *  =====================================================================
068: *
069: *     .. Parameters ..
070:       DOUBLE PRECISION   ONE, ZERO
071:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
072: *     ..
073: *     .. Local Scalars ..
074:       LOGICAL            APPLYLEFT
075:       INTEGER            I, LASTV, LASTC
076: *     ..
077: *     .. External Subroutines ..
078:       EXTERNAL           DGEMV, DGER
079: *     ..
080: *     .. External Functions ..
081:       LOGICAL            LSAME
084: *     ..
085: *     .. Executable Statements ..
086: *
087:       APPLYLEFT = LSAME( SIDE, 'L' )
088:       LASTV = 0
089:       LASTC = 0
090:       IF( TAU.NE.ZERO ) THEN
091: !     Set up variables for scanning V.  LASTV begins pointing to the end
092: !     of V.
093:          IF( APPLYLEFT ) THEN
094:             LASTV = M
095:          ELSE
096:             LASTV = N
097:          END IF
098:          IF( INCV.GT.0 ) THEN
099:             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)' * 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)'
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)'
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
154: ```