1*> \brief \b DLARZT forms the triangular factor T of a block reflector H = I - vtvH.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLARZT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarzt.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarzt.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarzt.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
22*
23*       .. Scalar Arguments ..
24*       CHARACTER          DIRECT, STOREV
25*       INTEGER            K, LDT, LDV, N
26*       ..
27*       .. Array Arguments ..
28*       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> DLARZT forms the triangular factor T of a real block reflector
38*> H of order > n, which is defined as a product of k elementary
39*> reflectors.
40*>
41*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
42*>
43*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
44*>
45*> If STOREV = 'C', the vector which defines the elementary reflector
46*> H(i) is stored in the i-th column of the array V, and
47*>
48*>    H  =  I - V * T * V**T
49*>
50*> If STOREV = 'R', the vector which defines the elementary reflector
51*> H(i) is stored in the i-th row of the array V, and
52*>
53*>    H  =  I - V**T * T * V
54*>
55*> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
56*> \endverbatim
57*
58*  Arguments:
59*  ==========
60*
61*> \param[in] DIRECT
62*> \verbatim
63*>          DIRECT is CHARACTER*1
64*>          Specifies the order in which the elementary reflectors are
65*>          multiplied to form the block reflector:
66*>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
67*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
68*> \endverbatim
69*>
70*> \param[in] STOREV
71*> \verbatim
72*>          STOREV is CHARACTER*1
73*>          Specifies how the vectors which define the elementary
74*>          reflectors are stored (see also Further Details):
75*>          = 'C': columnwise                        (not supported yet)
76*>          = 'R': rowwise
77*> \endverbatim
78*>
79*> \param[in] N
80*> \verbatim
81*>          N is INTEGER
82*>          The order of the block reflector H. N >= 0.
83*> \endverbatim
84*>
85*> \param[in] K
86*> \verbatim
87*>          K is INTEGER
88*>          The order of the triangular factor T (= the number of
89*>          elementary reflectors). K >= 1.
90*> \endverbatim
91*>
92*> \param[in,out] V
93*> \verbatim
94*>          V is DOUBLE PRECISION array, dimension
95*>                               (LDV,K) if STOREV = 'C'
96*>                               (LDV,N) if STOREV = 'R'
97*>          The matrix V. See further details.
98*> \endverbatim
99*>
100*> \param[in] LDV
101*> \verbatim
102*>          LDV is INTEGER
103*>          The leading dimension of the array V.
104*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
105*> \endverbatim
106*>
107*> \param[in] TAU
108*> \verbatim
109*>          TAU is DOUBLE PRECISION array, dimension (K)
110*>          TAU(i) must contain the scalar factor of the elementary
111*>          reflector H(i).
112*> \endverbatim
113*>
114*> \param[out] T
115*> \verbatim
116*>          T is DOUBLE PRECISION array, dimension (LDT,K)
117*>          The k by k triangular factor T of the block reflector.
118*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
119*>          lower triangular. The rest of the array is not used.
120*> \endverbatim
121*>
122*> \param[in] LDT
123*> \verbatim
124*>          LDT is INTEGER
125*>          The leading dimension of the array T. LDT >= K.
126*> \endverbatim
127*
128*  Authors:
129*  ========
130*
131*> \author Univ. of Tennessee
132*> \author Univ. of California Berkeley
133*> \author Univ. of Colorado Denver
134*> \author NAG Ltd.
135*
136*> \ingroup doubleOTHERcomputational
137*
138*> \par Contributors:
139*  ==================
140*>
141*>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
142*
143*> \par Further Details:
144*  =====================
145*>
146*> \verbatim
147*>
148*>  The shape of the matrix V and the storage of the vectors which define
149*>  the H(i) is best illustrated by the following example with n = 5 and
150*>  k = 3. The elements equal to 1 are not stored; the corresponding
151*>  array elements are modified but restored on exit. The rest of the
152*>  array is not used.
153*>
154*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
155*>
156*>                                              ______V_____
157*>         ( v1 v2 v3 )                        /            \
158*>         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )
159*>     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   )
160*>         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )
161*>         ( v1 v2 v3 )
162*>            .  .  .
163*>            .  .  .
164*>            1  .  .
165*>               1  .
166*>                  1
167*>
168*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
169*>
170*>                                                        ______V_____
171*>            1                                          /            \
172*>            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )
173*>            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )
174*>            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )
175*>            .  .  .
176*>         ( v1 v2 v3 )
177*>         ( v1 v2 v3 )
178*>     V = ( v1 v2 v3 )
179*>         ( v1 v2 v3 )
180*>         ( v1 v2 v3 )
181*> \endverbatim
182*>
183*  =====================================================================
184      SUBROUTINE DLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
185*
186*  -- LAPACK computational routine --
187*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
188*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190*     .. Scalar Arguments ..
191      CHARACTER          DIRECT, STOREV
192      INTEGER            K, LDT, LDV, N
193*     ..
194*     .. Array Arguments ..
195      DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
196*     ..
197*
198*  =====================================================================
199*
200*     .. Parameters ..
201      DOUBLE PRECISION   ZERO
202      PARAMETER          ( ZERO = 0.0D+0 )
203*     ..
204*     .. Local Scalars ..
205      INTEGER            I, INFO, J
206*     ..
207*     .. External Subroutines ..
208      EXTERNAL           DGEMV, DTRMV, XERBLA
209*     ..
210*     .. External Functions ..
211      LOGICAL            LSAME
212      EXTERNAL           LSAME
213*     ..
214*     .. Executable Statements ..
215*
216*     Check for currently supported options
217*
218      INFO = 0
219      IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
220         INFO = -1
221      ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
222         INFO = -2
223      END IF
224      IF( INFO.NE.0 ) THEN
225         CALL XERBLA( 'DLARZT', -INFO )
226         RETURN
227      END IF
228*
229      DO 20 I = K, 1, -1
230         IF( TAU( I ).EQ.ZERO ) THEN
231*
232*           H(i)  =  I
233*
234            DO 10 J = I, K
235               T( J, I ) = ZERO
236   10       CONTINUE
237         ELSE
238*
239*           general case
240*
241            IF( I.LT.K ) THEN
242*
243*              T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)**T
244*
245               CALL DGEMV( 'No transpose', K-I, N, -TAU( I ),
246     $                     V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO,
247     $                     T( I+1, I ), 1 )
248*
249*              T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i)
250*
251               CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
252     $                     T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
253            END IF
254            T( I, I ) = TAU( I )
255         END IF
256   20 CONTINUE
257      RETURN
258*
259*     End of DLARZT
260*
261      END
262