LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
slarft.f
Go to the documentation of this file.
1*> \brief \b SLARFT 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*> Download SLARFT + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarft.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarft.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarft.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
20*
21* .. Scalar Arguments ..
22* CHARACTER DIRECT, STOREV
23* INTEGER K, LDT, LDV, N
24* ..
25* .. Array Arguments ..
26* REAL T( LDT, * ), TAU( * ), V( LDV, * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> SLARFT forms the triangular factor T of a real block reflector H
36*> of order n, which is defined as a product of k elementary reflectors.
37*>
38*> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
39*>
40*> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
41*>
42*> If STOREV = 'C', the vector which defines the elementary reflector
43*> H(i) is stored in the i-th column of the array V, and
44*>
45*> H = I - V * T * V**T
46*>
47*> If STOREV = 'R', the vector which defines the elementary reflector
48*> H(i) is stored in the i-th row of the array V, and
49*>
50*> H = I - V**T * T * V
51*> \endverbatim
52*
53* Arguments:
54* ==========
55*
56*> \param[in] DIRECT
57*> \verbatim
58*> DIRECT is CHARACTER*1
59*> Specifies the order in which the elementary reflectors are
60*> multiplied to form the block reflector:
61*> = 'F': H = H(1) H(2) . . . H(k) (Forward)
62*> = 'B': H = H(k) . . . H(2) H(1) (Backward)
63*> \endverbatim
64*>
65*> \param[in] STOREV
66*> \verbatim
67*> STOREV is CHARACTER*1
68*> Specifies how the vectors which define the elementary
69*> reflectors are stored (see also Further Details):
70*> = 'C': columnwise
71*> = 'R': rowwise
72*> \endverbatim
73*>
74*> \param[in] N
75*> \verbatim
76*> N is INTEGER
77*> The order of the block reflector H. N >= 0.
78*> \endverbatim
79*>
80*> \param[in] K
81*> \verbatim
82*> K is INTEGER
83*> The order of the triangular factor T (= the number of
84*> elementary reflectors). K >= 1.
85*> \endverbatim
86*>
87*> \param[in] V
88*> \verbatim
89*> V is REAL array, dimension
90*> (LDV,K) if STOREV = 'C'
91*> (LDV,N) if STOREV = 'R'
92*> The matrix V. See further details.
93*> \endverbatim
94*>
95*> \param[in] LDV
96*> \verbatim
97*> LDV is INTEGER
98*> The leading dimension of the array V.
99*> If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
100*> \endverbatim
101*>
102*> \param[in] TAU
103*> \verbatim
104*> TAU is REAL array, dimension (K)
105*> TAU(i) must contain the scalar factor of the elementary
106*> reflector H(i).
107*> \endverbatim
108*>
109*> \param[out] T
110*> \verbatim
111*> T is REAL array, dimension (LDT,K)
112*> The k by k triangular factor T of the block reflector.
113*> If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
114*> lower triangular. The rest of the array is not used.
115*> \endverbatim
116*>
117*> \param[in] LDT
118*> \verbatim
119*> LDT is INTEGER
120*> The leading dimension of the array T. LDT >= K.
121*> \endverbatim
122*
123* Authors:
124* ========
125*
126*> \author Univ. of Tennessee
127*> \author Univ. of California Berkeley
128*> \author Johnathan Rhyne, Univ. of Colorado Denver (original author, 2024)
129*> \author NAG Ltd.
130*
131*> \ingroup larft
132*
133*> \par Further Details:
134* =====================
135*>
136*> \verbatim
137*>
138*> The shape of the matrix V and the storage of the vectors which define
139*> the H(i) is best illustrated by the following example with n = 5 and
140*> k = 3. The elements equal to 1 are not stored.
141*>
142*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R':
143*>
144*> V = ( 1 ) V = ( 1 v1 v1 v1 v1 )
145*> ( v1 1 ) ( 1 v2 v2 v2 )
146*> ( v1 v2 1 ) ( 1 v3 v3 )
147*> ( v1 v2 v3 )
148*> ( v1 v2 v3 )
149*>
150*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R':
151*>
152*> V = ( v1 v2 v3 ) V = ( v1 v1 1 )
153*> ( v1 v2 v3 ) ( v2 v2 v2 1 )
154*> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 )
155*> ( 1 v3 )
156*> ( 1 )
157*> \endverbatim
158*>
159* =====================================================================
160 RECURSIVE SUBROUTINE slarft( DIRECT, STOREV, N, K, V, LDV,
161 $ TAU, T, LDT )
162*
163* -- LAPACK auxiliary routine --
164* -- LAPACK is a software package provided by Univ. of Tennessee, --
165* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*
167* .. Scalar Arguments
168*
169 CHARACTER direct, storev
170 INTEGER k, ldt, ldv, n
171* ..
172* .. Array Arguments ..
173*
174 REAL t( ldt, * ), tau( * ), v( ldv, * )
175* ..
176*
177* .. Parameters ..
178*
179 REAL one, neg_one, zero
180 parameter(one=1.0e+0, zero = 0.0e+0, neg_one=-1.0e+0)
181*
182* .. Local Scalars ..
183*
184 INTEGER i,j,l
185 LOGICAL qr,lq,ql,dirf,colv
186*
187* .. External Subroutines ..
188*
189 EXTERNAL strmm,sgemm,slacpy
190*
191* .. External Functions..
192*
193 LOGICAL lsame
194 EXTERNAL lsame
195*
196* The general scheme used is inspired by the approach inside DGEQRT3
197* which was (at the time of writing this code):
198* Based on the algorithm of Elmroth and Gustavson,
199* IBM J. Res. Develop. Vol 44 No. 4 July 2000.
200* ..
201* .. Executable Statements ..
202*
203* Quick return if possible
204*
205 IF(n.EQ.0.OR.k.EQ.0) THEN
206 RETURN
207 END IF
208*
209* Base case
210*
211 IF(n.EQ.1.OR.k.EQ.1) THEN
212 t(1,1) = tau(1)
213 RETURN
214 END IF
215*
216* Beginning of executable statements
217*
218 l = k / 2
219*
220* Determine what kind of Q we need to compute
221* We assume that if the user doesn't provide 'F' for DIRECT,
222* then they meant to provide 'B' and if they don't provide
223* 'C' for STOREV, then they meant to provide 'R'
224*
225 dirf = lsame(direct,'F')
226 colv = lsame(storev,'C')
227*
228* QR happens when we have forward direction in column storage
229*
230 qr = dirf.AND.colv
231*
232* LQ happens when we have forward direction in row storage
233*
234 lq = dirf.AND.(.NOT.colv)
235*
236* QL happens when we have backward direction in column storage
237*
238 ql = (.NOT.dirf).AND.colv
239*
240* The last case is RQ. Due to how we structured this, if the
241* above 3 are false, then RQ must be true, so we never store
242* this
243* RQ happens when we have backward direction in row storage
244* RQ = (.NOT.DIRF).AND.(.NOT.COLV)
245*
246 IF(qr) THEN
247*
248* Break V apart into 6 components
249*
250* V = |---------------|
251* |V_{1,1} 0 |
252* |V_{2,1} V_{2,2}|
253* |V_{3,1} V_{3,2}|
254* |---------------|
255*
256* V_{1,1}\in\R^{l,l} unit lower triangular
257* V_{2,1}\in\R^{k-l,l} rectangular
258* V_{3,1}\in\R^{n-k,l} rectangular
259*
260* V_{2,2}\in\R^{k-l,k-l} unit lower triangular
261* V_{3,2}\in\R^{n-k,k-l} rectangular
262*
263* We will construct the T matrix
264* T = |---------------|
265* |T_{1,1} T_{1,2}|
266* |0 T_{2,2}|
267* |---------------|
268*
269* T is the triangular factor obtained from block reflectors.
270* To motivate the structure, assume we have already computed T_{1,1}
271* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2
272*
273* T_{1,1}\in\R^{l, l} upper triangular
274* T_{2,2}\in\R^{k-l, k-l} upper triangular
275* T_{1,2}\in\R^{l, k-l} rectangular
276*
277* Where l = floor(k/2)
278*
279* Then, consider the product:
280*
281* (I - V_1*T_{1,1}*V_1')*(I - V_2*T_{2,2}*V_2')
282* = I - V_1*T_{1,1}*V_1' - V_2*T_{2,2}*V_2' + V_1*T_{1,1}*V_1'*V_2*T_{2,2}*V_2'
283*
284* Define T_{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2}
285*
286* Then, we can define the matrix V as
287* V = |-------|
288* |V_1 V_2|
289* |-------|
290*
291* So, our product is equivalent to the matrix product
292* I - V*T*V'
293* This means, we can compute T_{1,1} and T_{2,2}, then use this information
294* to compute T_{1,2}
295*
296* Compute T_{1,1} recursively
297*
298 CALL slarft(direct, storev, n, l, v, ldv, tau, t, ldt)
299*
300* Compute T_{2,2} recursively
301*
302 CALL slarft(direct, storev, n-l, k-l, v(l+1, l+1), ldv,
303 $ tau(l+1), t(l+1, l+1), ldt)
304*
305* Compute T_{1,2}
306* T_{1,2} = V_{2,1}'
307*
308 DO j = 1, l
309 DO i = 1, k-l
310 t(j, l+i) = v(l+i, j)
311 END DO
312 END DO
313*
314* T_{1,2} = T_{1,2}*V_{2,2}
315*
316 CALL strmm('Right', 'Lower', 'No transpose', 'Unit', l,
317 $ k-l, one, v(l+1, l+1), ldv, t(1, l+1), ldt)
318
319*
320* T_{1,2} = V_{3,1}'*V_{3,2} + T_{1,2}
321* Note: We assume K <= N, and GEMM will do nothing if N=K
322*
323 CALL sgemm('Transpose', 'No transpose', l, k-l, n-k, one,
324 $ v(k+1, 1), ldv, v(k+1, l+1), ldv, one,
325 $ t(1, l+1), ldt)
326*
327* At this point, we have that T_{1,2} = V_1'*V_2
328* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2}
329* respectively.
330*
331* T_{1,2} = -T_{1,1}*T_{1,2}
332*
333 CALL strmm('Left', 'Upper', 'No transpose', 'Non-unit', l,
334 $ k-l, neg_one, t, ldt, t(1, l+1), ldt)
335*
336* T_{1,2} = T_{1,2}*T_{2,2}
337*
338 CALL strmm('Right', 'Upper', 'No transpose', 'Non-unit', l,
339 $ k-l, one, t(l+1, l+1), ldt, t(1, l+1), ldt)
340
341 ELSE IF(lq) THEN
342*
343* Break V apart into 6 components
344*
345* V = |----------------------|
346* |V_{1,1} V_{1,2} V{1,3}|
347* |0 V_{2,2} V{2,3}|
348* |----------------------|
349*
350* V_{1,1}\in\R^{l,l} unit upper triangular
351* V_{1,2}\in\R^{l,k-l} rectangular
352* V_{1,3}\in\R^{l,n-k} rectangular
353*
354* V_{2,2}\in\R^{k-l,k-l} unit upper triangular
355* V_{2,3}\in\R^{k-l,n-k} rectangular
356*
357* Where l = floor(k/2)
358*
359* We will construct the T matrix
360* T = |---------------|
361* |T_{1,1} T_{1,2}|
362* |0 T_{2,2}|
363* |---------------|
364*
365* T is the triangular factor obtained from block reflectors.
366* To motivate the structure, assume we have already computed T_{1,1}
367* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2
368*
369* T_{1,1}\in\R^{l, l} upper triangular
370* T_{2,2}\in\R^{k-l, k-l} upper triangular
371* T_{1,2}\in\R^{l, k-l} rectangular
372*
373* Then, consider the product:
374*
375* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2)
376* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2
377*
378* Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2}
379*
380* Then, we can define the matrix V as
381* V = |---|
382* |V_1|
383* |V_2|
384* |---|
385*
386* So, our product is equivalent to the matrix product
387* I - V'*T*V
388* This means, we can compute T_{1,1} and T_{2,2}, then use this information
389* to compute T_{1,2}
390*
391* Compute T_{1,1} recursively
392*
393 CALL slarft(direct, storev, n, l, v, ldv, tau, t, ldt)
394*
395* Compute T_{2,2} recursively
396*
397 CALL slarft(direct, storev, n-l, k-l, v(l+1, l+1), ldv,
398 $ tau(l+1), t(l+1, l+1), ldt)
399
400*
401* Compute T_{1,2}
402* T_{1,2} = V_{1,2}
403*
404 CALL slacpy('All', l, k-l, v(1, l+1), ldv, t(1, l+1), ldt)
405*
406* T_{1,2} = T_{1,2}*V_{2,2}'
407*
408 CALL strmm('Right', 'Upper', 'Transpose', 'Unit', l, k-l,
409 $ one, v(l+1, l+1), ldv, t(1, l+1), ldt)
410
411*
412* T_{1,2} = V_{1,3}*V_{2,3}' + T_{1,2}
413* Note: We assume K <= N, and GEMM will do nothing if N=K
414*
415 CALL sgemm('No transpose', 'Transpose', l, k-l, n-k, one,
416 $ v(1, k+1), ldv, v(l+1, k+1), ldv, one,
417 $ t(1, l+1), ldt)
418*
419* At this point, we have that T_{1,2} = V_1*V_2'
420* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2}
421* respectively.
422*
423* T_{1,2} = -T_{1,1}*T_{1,2}
424*
425 CALL strmm('Left', 'Upper', 'No transpose', 'Non-unit', l,
426 $ k-l, neg_one, t, ldt, t(1, l+1), ldt)
427
428*
429* T_{1,2} = T_{1,2}*T_{2,2}
430*
431 CALL strmm('Right', 'Upper', 'No transpose', 'Non-unit', l,
432 $ k-l, one, t(l+1, l+1), ldt, t(1, l+1), ldt)
433 ELSE IF(ql) THEN
434*
435* Break V apart into 6 components
436*
437* V = |---------------|
438* |V_{1,1} V_{1,2}|
439* |V_{2,1} V_{2,2}|
440* |0 V_{3,2}|
441* |---------------|
442*
443* V_{1,1}\in\R^{n-k,k-l} rectangular
444* V_{2,1}\in\R^{k-l,k-l} unit upper triangular
445*
446* V_{1,2}\in\R^{n-k,l} rectangular
447* V_{2,2}\in\R^{k-l,l} rectangular
448* V_{3,2}\in\R^{l,l} unit upper triangular
449*
450* We will construct the T matrix
451* T = |---------------|
452* |T_{1,1} 0 |
453* |T_{2,1} T_{2,2}|
454* |---------------|
455*
456* T is the triangular factor obtained from block reflectors.
457* To motivate the structure, assume we have already computed T_{1,1}
458* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2
459*
460* T_{1,1}\in\R^{k-l, k-l} non-unit lower triangular
461* T_{2,2}\in\R^{l, l} non-unit lower triangular
462* T_{2,1}\in\R^{k-l, l} rectangular
463*
464* Where l = floor(k/2)
465*
466* Then, consider the product:
467*
468* (I - V_2*T_{2,2}*V_2')*(I - V_1*T_{1,1}*V_1')
469* = I - V_2*T_{2,2}*V_2' - V_1*T_{1,1}*V_1' + V_2*T_{2,2}*V_2'*V_1*T_{1,1}*V_1'
470*
471* Define T_{2,1} = -T_{2,2}*V_2'*V_1*T_{1,1}
472*
473* Then, we can define the matrix V as
474* V = |-------|
475* |V_1 V_2|
476* |-------|
477*
478* So, our product is equivalent to the matrix product
479* I - V*T*V'
480* This means, we can compute T_{1,1} and T_{2,2}, then use this information
481* to compute T_{2,1}
482*
483* Compute T_{1,1} recursively
484*
485 CALL slarft(direct, storev, n-l, k-l, v, ldv, tau, t, ldt)
486*
487* Compute T_{2,2} recursively
488*
489 CALL slarft(direct, storev, n, l, v(1, k-l+1), ldv,
490 $ tau(k-l+1), t(k-l+1, k-l+1), ldt)
491*
492* Compute T_{2,1}
493* T_{2,1} = V_{2,2}'
494*
495 DO j = 1, k-l
496 DO i = 1, l
497 t(k-l+i, j) = v(n-k+j, k-l+i)
498 END DO
499 END DO
500*
501* T_{2,1} = T_{2,1}*V_{2,1}
502*
503 CALL strmm('Right', 'Upper', 'No transpose', 'Unit', l,
504 $ k-l, one, v(n-k+1, 1), ldv, t(k-l+1, 1), ldt)
505
506*
507* T_{2,1} = V_{2,2}'*V_{2,1} + T_{2,1}
508* Note: We assume K <= N, and GEMM will do nothing if N=K
509*
510 CALL sgemm('Transpose', 'No transpose', l, k-l, n-k, one,
511 $ v(1, k-l+1), ldv, v, ldv, one, t(k-l+1, 1),
512 $ ldt)
513*
514* At this point, we have that T_{2,1} = V_2'*V_1
515* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1}
516* respectively.
517*
518* T_{2,1} = -T_{2,2}*T_{2,1}
519*
520 CALL strmm('Left', 'Lower', 'No transpose', 'Non-unit', l,
521 $ k-l, neg_one, t(k-l+1, k-l+1), ldt,
522 $ t(k-l+1, 1), ldt)
523*
524* T_{2,1} = T_{2,1}*T_{1,1}
525*
526 CALL strmm('Right', 'Lower', 'No transpose', 'Non-unit', l,
527 $ k-l, one, t, ldt, t(k-l+1, 1), ldt)
528 ELSE
529*
530* Else means RQ case
531*
532* Break V apart into 6 components
533*
534* V = |-----------------------|
535* |V_{1,1} V_{1,2} 0 |
536* |V_{2,1} V_{2,2} V_{2,3}|
537* |-----------------------|
538*
539* V_{1,1}\in\R^{k-l,n-k} rectangular
540* V_{1,2}\in\R^{k-l,k-l} unit lower triangular
541*
542* V_{2,1}\in\R^{l,n-k} rectangular
543* V_{2,2}\in\R^{l,k-l} rectangular
544* V_{2,3}\in\R^{l,l} unit lower triangular
545*
546* We will construct the T matrix
547* T = |---------------|
548* |T_{1,1} 0 |
549* |T_{2,1} T_{2,2}|
550* |---------------|
551*
552* T is the triangular factor obtained from block reflectors.
553* To motivate the structure, assume we have already computed T_{1,1}
554* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2
555*
556* T_{1,1}\in\R^{k-l, k-l} non-unit lower triangular
557* T_{2,2}\in\R^{l, l} non-unit lower triangular
558* T_{2,1}\in\R^{k-l, l} rectangular
559*
560* Where l = floor(k/2)
561*
562* Then, consider the product:
563*
564* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1)
565* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1
566*
567* Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1}
568*
569* Then, we can define the matrix V as
570* V = |---|
571* |V_1|
572* |V_2|
573* |---|
574*
575* So, our product is equivalent to the matrix product
576* I - V'TV
577* This means, we can compute T_{1,1} and T_{2,2}, then use this information
578* to compute T_{2,1}
579*
580* Compute T_{1,1} recursively
581*
582 CALL slarft(direct, storev, n-l, k-l, v, ldv, tau, t, ldt)
583*
584* Compute T_{2,2} recursively
585*
586 CALL slarft(direct, storev, n, l, v(k-l+1, 1), ldv,
587 $ tau(k-l+1), t(k-l+1, k-l+1), ldt)
588*
589* Compute T_{2,1}
590* T_{2,1} = V_{2,2}
591*
592 CALL slacpy('All', l, k-l, v(k-l+1, n-k+1), ldv,
593 $ t(k-l+1, 1), ldt)
594
595*
596* T_{2,1} = T_{2,1}*V_{1,2}'
597*
598 CALL strmm('Right', 'Lower', 'Transpose', 'Unit', l, k-l,
599 $ one, v(1, n-k+1), ldv, t(k-l+1, 1), ldt)
600
601*
602* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1}
603* Note: We assume K <= N, and GEMM will do nothing if N=K
604*
605 CALL sgemm('No transpose', 'Transpose', l, k-l, n-k, one,
606 $ v(k-l+1, 1), ldv, v, ldv, one, t(k-l+1, 1),
607 $ ldt)
608
609*
610* At this point, we have that T_{2,1} = V_2*V_1'
611* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1}
612* respectively.
613*
614* T_{2,1} = -T_{2,2}*T_{2,1}
615*
616 CALL strmm('Left', 'Lower', 'No tranpose', 'Non-unit', l,
617 $ k-l, neg_one, t(k-l+1, k-l+1), ldt,
618 $ t(k-l+1, 1), ldt)
619
620*
621* T_{2,1} = T_{2,1}*T_{1,1}
622*
623 CALL strmm('Right', 'Lower', 'No tranpose', 'Non-unit', l,
624 $ k-l, one, t, ldt, t(k-l+1, 1), ldt)
625 END IF
626 END SUBROUTINE
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
slacpy
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:101
slarft
recursive subroutine slarft(direct, storev, n, k, v, ldv, tau, t, ldt)
SLARFT forms the triangular factor T of a block reflector H = I - vtvH
Definition slarft.f:162
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177