• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1*> \brief \b ZLARFT
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLARFT + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlarft.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlarft.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlarft.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE ZLARFT( 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*       COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
29*       ..
30*
31*
32*> \par Purpose:
33*  =============
34*>
35*> \verbatim
36*>
37*> ZLARFT forms the triangular factor T of a complex block reflector H
38*> of order n, which is defined as a product of k elementary reflectors.
39*>
40*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
41*>
42*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
43*>
44*> If STOREV = 'C', the vector which defines the elementary reflector
45*> H(i) is stored in the i-th column of the array V, and
46*>
47*>    H  =  I - V * T * V**H
48*>
49*> If STOREV = 'R', the vector which defines the elementary reflector
50*> H(i) is stored in the i-th row of the array V, and
51*>
52*>    H  =  I - V**H * T * V
53*> \endverbatim
54*
55*  Arguments:
56*  ==========
57*
58*> \param[in] DIRECT
59*> \verbatim
60*>          DIRECT is CHARACTER*1
61*>          Specifies the order in which the elementary reflectors are
62*>          multiplied to form the block reflector:
63*>          = 'F': H = H(1) H(2) . . . H(k) (Forward)
64*>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
65*> \endverbatim
66*>
67*> \param[in] STOREV
68*> \verbatim
69*>          STOREV is CHARACTER*1
70*>          Specifies how the vectors which define the elementary
71*>          reflectors are stored (see also Further Details):
72*>          = 'C': columnwise
73*>          = 'R': rowwise
74*> \endverbatim
75*>
76*> \param[in] N
77*> \verbatim
78*>          N is INTEGER
79*>          The order of the block reflector H. N >= 0.
80*> \endverbatim
81*>
82*> \param[in] K
83*> \verbatim
84*>          K is INTEGER
85*>          The order of the triangular factor T (= the number of
86*>          elementary reflectors). K >= 1.
87*> \endverbatim
88*>
89*> \param[in] V
90*> \verbatim
91*>          V is COMPLEX*16 array, dimension
92*>                               (LDV,K) if STOREV = 'C'
93*>                               (LDV,N) if STOREV = 'R'
94*>          The matrix V. See further details.
95*> \endverbatim
96*>
97*> \param[in] LDV
98*> \verbatim
99*>          LDV is INTEGER
100*>          The leading dimension of the array V.
101*>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
102*> \endverbatim
103*>
104*> \param[in] TAU
105*> \verbatim
106*>          TAU is COMPLEX*16 array, dimension (K)
107*>          TAU(i) must contain the scalar factor of the elementary
108*>          reflector H(i).
109*> \endverbatim
110*>
111*> \param[out] T
112*> \verbatim
113*>          T is COMPLEX*16 array, dimension (LDT,K)
114*>          The k by k triangular factor T of the block reflector.
115*>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
116*>          lower triangular. The rest of the array is not used.
117*> \endverbatim
118*>
119*> \param[in] LDT
120*> \verbatim
121*>          LDT is INTEGER
122*>          The leading dimension of the array T. LDT >= K.
123*> \endverbatim
124*
125*  Authors:
126*  ========
127*
128*> \author Univ. of Tennessee
129*> \author Univ. of California Berkeley
130*> \author Univ. of Colorado Denver
131*> \author NAG Ltd.
132*
133*> \date April 2012
134*
135*> \ingroup complex16OTHERauxiliary
136*
137*> \par Further Details:
138*  =====================
139*>
140*> \verbatim
141*>
142*>  The shape of the matrix V and the storage of the vectors which define
143*>  the H(i) is best illustrated by the following example with n = 5 and
144*>  k = 3. The elements equal to 1 are not stored.
145*>
146*>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
147*>
148*>               V = (  1       )                 V = (  1 v1 v1 v1 v1 )
149*>                   ( v1  1    )                     (     1 v2 v2 v2 )
150*>                   ( v1 v2  1 )                     (        1 v3 v3 )
151*>                   ( v1 v2 v3 )
152*>                   ( v1 v2 v3 )
153*>
154*>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
155*>
156*>               V = ( v1 v2 v3 )                 V = ( v1 v1  1       )
157*>                   ( v1 v2 v3 )                     ( v2 v2 v2  1    )
158*>                   (  1 v2 v3 )                     ( v3 v3 v3 v3  1 )
159*>                   (     1 v3 )
160*>                   (        1 )
161*> \endverbatim
162*>
163*  =====================================================================
164      SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
165*
166*  -- LAPACK auxiliary routine (version 3.4.1) --
167*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
168*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169*     April 2012
170*
171*     .. Scalar Arguments ..
172      CHARACTER          DIRECT, STOREV
173      INTEGER            K, LDT, LDV, N
174*     ..
175*     .. Array Arguments ..
176      COMPLEX*16         T( LDT, * ), TAU( * ), V( LDV, * )
177*     ..
178*
179*  =====================================================================
180*
181*     .. Parameters ..
182      COMPLEX*16         ONE, ZERO
183      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
184     $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
185*     ..
186*     .. Local Scalars ..
187      INTEGER            I, J, PREVLASTV, LASTV
188*     ..
189*     .. External Subroutines ..
190      EXTERNAL           ZGEMV, ZLACGV, ZTRMV
191*     ..
192*     .. External Functions ..
193      LOGICAL            LSAME
194      EXTERNAL           LSAME
195*     ..
196*     .. Executable Statements ..
197*
198*     Quick return if possible
199*
200      IF( N.EQ.0 )
201     $   RETURN
202*
203      IF( LSAME( DIRECT, 'F' ) ) THEN
204         PREVLASTV = N
205         DO I = 1, K
206            PREVLASTV = MAX( PREVLASTV, I )
207            IF( TAU( I ).EQ.ZERO ) THEN
208*
209*              H(i)  =  I
210*
211               DO J = 1, I
212                  T( J, I ) = ZERO
213               END DO
214            ELSE
215*
216*              general case
217*
218               IF( LSAME( STOREV, 'C' ) ) THEN
219*                 Skip any trailing zeros.
220                  DO LASTV = N, I+1, -1
221                     IF( V( LASTV, I ).NE.ZERO ) EXIT
222                  END DO
223                  DO J = 1, I-1
224                     T( J, I ) = -TAU( I ) * CONJG( V( I , J ) )
225                  END DO
226                  J = MIN( LASTV, PREVLASTV )
227*
228*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i)
229*
230                  CALL ZGEMV( 'Conjugate transpose', J-I, I-1,
231     $                        -TAU( I ), V( I+1, 1 ), LDV,
232     $                        V( I+1, I ), 1, ONE, T( 1, I ), 1 )
233               ELSE
234*                 Skip any trailing zeros.
235                  DO LASTV = N, I+1, -1
236                     IF( V( I, LASTV ).NE.ZERO ) EXIT
237                  END DO
238                  DO J = 1, I-1
239                     T( J, I ) = -TAU( I ) * V( J , I )
240                  END DO
241                  J = MIN( LASTV, PREVLASTV )
242*
243*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H
244*
245                  CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ),
246     $                        V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
247     $                        ONE, T( 1, I ), LDT )
248               END IF
249*
250*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
251*
252               CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T,
253     $                     LDT, T( 1, I ), 1 )
254               T( I, I ) = TAU( I )
255               IF( I.GT.1 ) THEN
256                  PREVLASTV = MAX( PREVLASTV, LASTV )
257               ELSE
258                  PREVLASTV = LASTV
259               END IF
260             END IF
261         END DO
262      ELSE
263         PREVLASTV = 1
264         DO I = K, 1, -1
265            IF( TAU( I ).EQ.ZERO ) THEN
266*
267*              H(i)  =  I
268*
269               DO J = I, K
270                  T( J, I ) = ZERO
271               END DO
272            ELSE
273*
274*              general case
275*
276               IF( I.LT.K ) THEN
277                  IF( LSAME( STOREV, 'C' ) ) THEN
278*                    Skip any leading zeros.
279                     DO LASTV = 1, I-1
280                        IF( V( LASTV, I ).NE.ZERO ) EXIT
281                     END DO
282                     DO J = I+1, K
283                        T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) )
284                     END DO
285                     J = MAX( LASTV, PREVLASTV )
286*
287*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i)
288*
289                     CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I,
290     $                           -TAU( I ), V( J, I+1 ), LDV, V( J, I ),
291     $                           1, ONE, T( I+1, I ), 1 )
292                  ELSE
293*                    Skip any leading zeros.
294                     DO LASTV = 1, I-1
295                        IF( V( I, LASTV ).NE.ZERO ) EXIT
296                     END DO
297                     DO J = I+1, K
298                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
299                     END DO
300                     J = MAX( LASTV, PREVLASTV )
301*
302*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H
303*
304                     CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, -TAU( I ),
305     $                           V( I+1, J ), LDV, V( I, J ), LDV,
306     $                           ONE, T( I+1, I ), LDT )
307                  END IF
308*
309*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
310*
311                  CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
312     $                        T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
313                  IF( I.GT.1 ) THEN
314                     PREVLASTV = MIN( PREVLASTV, LASTV )
315                  ELSE
316                     PREVLASTV = LASTV
317                  END IF
318               END IF
319               T( I, I ) = TAU( I )
320            END IF
321         END DO
322      END IF
323      RETURN
324*
325*     End of ZLARFT
326*
327      END
328