LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
slarfb_gett.f
Go to the documentation of this file.
1 *> \brief \b SLARFB_GETT
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SLARFB_GETT + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slarfb_gett.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slarfb_gett.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slarfb_gett.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARFB_GETT( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
22 * $ WORK, LDWORK )
23 * IMPLICIT NONE
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER IDENT
27 * INTEGER K, LDA, LDB, LDT, LDWORK, M, N
28 * ..
29 * .. Array Arguments ..
30 * REAL A( LDA, * ), B( LDB, * ), T( LDT, * ),
31 * $ WORK( LDWORK, * )
32 * ..
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> SLARFB_GETT applies a real Householder block reflector H from the
40 *> left to a real (K+M)-by-N "triangular-pentagonal" matrix
41 *> composed of two block matrices: an upper trapezoidal K-by-N matrix A
42 *> stored in the array A, and a rectangular M-by-(N-K) matrix B, stored
43 *> in the array B. The block reflector H is stored in a compact
44 *> WY-representation, where the elementary reflectors are in the
45 *> arrays A, B and T. See Further Details section.
46 *> \endverbatim
47 *
48 * Arguments:
49 * ==========
50 *
51 *> \param[in] IDENT
52 *> \verbatim
53 *> IDENT is CHARACTER*1
54 *> If IDENT = not 'I', or not 'i', then V1 is unit
55 *> lower-triangular and stored in the left K-by-K block of
56 *> the input matrix A,
57 *> If IDENT = 'I' or 'i', then V1 is an identity matrix and
58 *> not stored.
59 *> See Further Details section.
60 *> \endverbatim
61 *>
62 *> \param[in] M
63 *> \verbatim
64 *> M is INTEGER
65 *> The number of rows of the matrix B.
66 *> M >= 0.
67 *> \endverbatim
68 *>
69 *> \param[in] N
70 *> \verbatim
71 *> N is INTEGER
72 *> The number of columns of the matrices A and B.
73 *> N >= 0.
74 *> \endverbatim
75 *>
76 *> \param[in] K
77 *> \verbatim
78 *> K is INTEGER
79 *> The number or rows of the matrix A.
80 *> K is also order of the matrix T, i.e. the number of
81 *> elementary reflectors whose product defines the block
82 *> reflector. 0 <= K <= N.
83 *> \endverbatim
84 *>
85 *> \param[in] T
86 *> \verbatim
87 *> T is REAL array, dimension (LDT,K)
88 *> The upper-triangular K-by-K matrix T in the representation
89 *> of the block reflector.
90 *> \endverbatim
91 *>
92 *> \param[in] LDT
93 *> \verbatim
94 *> LDT is INTEGER
95 *> The leading dimension of the array T. LDT >= K.
96 *> \endverbatim
97 *>
98 *> \param[in,out] A
99 *> \verbatim
100 *> A is REAL array, dimension (LDA,N)
101 *>
102 *> On entry:
103 *> a) In the K-by-N upper-trapezoidal part A: input matrix A.
104 *> b) In the columns below the diagonal: columns of V1
105 *> (ones are not stored on the diagonal).
106 *>
107 *> On exit:
108 *> A is overwritten by rectangular K-by-N product H*A.
109 *>
110 *> See Further Details section.
111 *> \endverbatim
112 *>
113 *> \param[in] LDA
114 *> \verbatim
115 *> LDB is INTEGER
116 *> The leading dimension of the array A. LDA >= max(1,K).
117 *> \endverbatim
118 *>
119 *> \param[in,out] B
120 *> \verbatim
121 *> B is REAL array, dimension (LDB,N)
122 *>
123 *> On entry:
124 *> a) In the M-by-(N-K) right block: input matrix B.
125 *> b) In the M-by-N left block: columns of V2.
126 *>
127 *> On exit:
128 *> B is overwritten by rectangular M-by-N product H*B.
129 *>
130 *> See Further Details section.
131 *> \endverbatim
132 *>
133 *> \param[in] LDB
134 *> \verbatim
135 *> LDB is INTEGER
136 *> The leading dimension of the array B. LDB >= max(1,M).
137 *> \endverbatim
138 *>
139 *> \param[out] WORK
140 *> \verbatim
141 *> WORK is REAL array,
142 *> dimension (LDWORK,max(K,N-K))
143 *> \endverbatim
144 *>
145 *> \param[in] LDWORK
146 *> \verbatim
147 *> LDWORK is INTEGER
148 *> The leading dimension of the array WORK. LDWORK>=max(1,K).
149 *>
150 *> \endverbatim
151 *
152 * Authors:
153 * ========
154 *
155 *> \author Univ. of Tennessee
156 *> \author Univ. of California Berkeley
157 *> \author Univ. of Colorado Denver
158 *> \author NAG Ltd.
159 *
160 *> \ingroup singleOTHERauxiliary
161 *
162 *> \par Contributors:
163 * ==================
164 *>
165 *> \verbatim
166 *>
167 *> November 2020, Igor Kozachenko,
168 *> Computer Science Division,
169 *> University of California, Berkeley
170 *>
171 *> \endverbatim
172 *
173 *> \par Further Details:
174 * =====================
175 *>
176 *> \verbatim
177 *>
178 *> (1) Description of the Algebraic Operation.
179 *>
180 *> The matrix A is a K-by-N matrix composed of two column block
181 *> matrices, A1, which is K-by-K, and A2, which is K-by-(N-K):
182 *> A = ( A1, A2 ).
183 *> The matrix B is an M-by-N matrix composed of two column block
184 *> matrices, B1, which is M-by-K, and B2, which is M-by-(N-K):
185 *> B = ( B1, B2 ).
186 *>
187 *> Perform the operation:
188 *>
189 *> ( A_out ) := H * ( A_in ) = ( I - V * T * V**T ) * ( A_in ) =
190 *> ( B_out ) ( B_in ) ( B_in )
191 *> = ( I - ( V1 ) * T * ( V1**T, V2**T ) ) * ( A_in )
192 *> ( V2 ) ( B_in )
193 *> On input:
194 *>
195 *> a) ( A_in ) consists of two block columns:
196 *> ( B_in )
197 *>
198 *> ( A_in ) = (( A1_in ) ( A2_in )) = (( A1_in ) ( A2_in ))
199 *> ( B_in ) (( B1_in ) ( B2_in )) (( 0 ) ( B2_in )),
200 *>
201 *> where the column blocks are:
202 *>
203 *> ( A1_in ) is a K-by-K upper-triangular matrix stored in the
204 *> upper triangular part of the array A(1:K,1:K).
205 *> ( B1_in ) is an M-by-K rectangular ZERO matrix and not stored.
206 *>
207 *> ( A2_in ) is a K-by-(N-K) rectangular matrix stored
208 *> in the array A(1:K,K+1:N).
209 *> ( B2_in ) is an M-by-(N-K) rectangular matrix stored
210 *> in the array B(1:M,K+1:N).
211 *>
212 *> b) V = ( V1 )
213 *> ( V2 )
214 *>
215 *> where:
216 *> 1) if IDENT == 'I',V1 is a K-by-K identity matrix, not stored;
217 *> 2) if IDENT != 'I',V1 is a K-by-K unit lower-triangular matrix,
218 *> stored in the lower-triangular part of the array
219 *> A(1:K,1:K) (ones are not stored),
220 *> and V2 is an M-by-K rectangular stored the array B(1:M,1:K),
221 *> (because on input B1_in is a rectangular zero
222 *> matrix that is not stored and the space is
223 *> used to store V2).
224 *>
225 *> c) T is a K-by-K upper-triangular matrix stored
226 *> in the array T(1:K,1:K).
227 *>
228 *> On output:
229 *>
230 *> a) ( A_out ) consists of two block columns:
231 *> ( B_out )
232 *>
233 *> ( A_out ) = (( A1_out ) ( A2_out ))
234 *> ( B_out ) (( B1_out ) ( B2_out )),
235 *>
236 *> where the column blocks are:
237 *>
238 *> ( A1_out ) is a K-by-K square matrix, or a K-by-K
239 *> upper-triangular matrix, if V1 is an
240 *> identity matrix. AiOut is stored in
241 *> the array A(1:K,1:K).
242 *> ( B1_out ) is an M-by-K rectangular matrix stored
243 *> in the array B(1:M,K:N).
244 *>
245 *> ( A2_out ) is a K-by-(N-K) rectangular matrix stored
246 *> in the array A(1:K,K+1:N).
247 *> ( B2_out ) is an M-by-(N-K) rectangular matrix stored
248 *> in the array B(1:M,K+1:N).
249 *>
250 *>
251 *> The operation above can be represented as the same operation
252 *> on each block column:
253 *>
254 *> ( A1_out ) := H * ( A1_in ) = ( I - V * T * V**T ) * ( A1_in )
255 *> ( B1_out ) ( 0 ) ( 0 )
256 *>
257 *> ( A2_out ) := H * ( A2_in ) = ( I - V * T * V**T ) * ( A2_in )
258 *> ( B2_out ) ( B2_in ) ( B2_in )
259 *>
260 *> If IDENT != 'I':
261 *>
262 *> The computation for column block 1:
263 *>
264 *> A1_out: = A1_in - V1*T*(V1**T)*A1_in
265 *>
266 *> B1_out: = - V2*T*(V1**T)*A1_in
267 *>
268 *> The computation for column block 2, which exists if N > K:
269 *>
270 *> A2_out: = A2_in - V1*T*( (V1**T)*A2_in + (V2**T)*B2_in )
271 *>
272 *> B2_out: = B2_in - V2*T*( (V1**T)*A2_in + (V2**T)*B2_in )
273 *>
274 *> If IDENT == 'I':
275 *>
276 *> The operation for column block 1:
277 *>
278 *> A1_out: = A1_in - V1*T**A1_in
279 *>
280 *> B1_out: = - V2*T**A1_in
281 *>
282 *> The computation for column block 2, which exists if N > K:
283 *>
284 *> A2_out: = A2_in - T*( A2_in + (V2**T)*B2_in )
285 *>
286 *> B2_out: = B2_in - V2*T*( A2_in + (V2**T)*B2_in )
287 *>
288 *> (2) Description of the Algorithmic Computation.
289 *>
290 *> In the first step, we compute column block 2, i.e. A2 and B2.
291 *> Here, we need to use the K-by-(N-K) rectangular workspace
292 *> matrix W2 that is of the same size as the matrix A2.
293 *> W2 is stored in the array WORK(1:K,1:(N-K)).
294 *>
295 *> In the second step, we compute column block 1, i.e. A1 and B1.
296 *> Here, we need to use the K-by-K square workspace matrix W1
297 *> that is of the same size as the as the matrix A1.
298 *> W1 is stored in the array WORK(1:K,1:K).
299 *>
300 *> NOTE: Hence, in this routine, we need the workspace array WORK
301 *> only of size WORK(1:K,1:max(K,N-K)) so it can hold both W2 from
302 *> the first step and W1 from the second step.
303 *>
304 *> Case (A), when V1 is unit lower-triangular, i.e. IDENT != 'I',
305 *> more computations than in the Case (B).
306 *>
307 *> if( IDENT != 'I' ) then
308 *> if ( N > K ) then
309 *> (First Step - column block 2)
310 *> col2_(1) W2: = A2
311 *> col2_(2) W2: = (V1**T) * W2 = (unit_lower_tr_of_(A1)**T) * W2
312 *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
313 *> col2_(4) W2: = T * W2
314 *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
315 *> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
316 *> col2_(7) A2: = A2 - W2
317 *> else
318 *> (Second Step - column block 1)
319 *> col1_(1) W1: = A1
320 *> col1_(2) W1: = (V1**T) * W1 = (unit_lower_tr_of_(A1)**T) * W1
321 *> col1_(3) W1: = T * W1
322 *> col1_(4) B1: = - V2 * W1 = - B1 * W1
323 *> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
324 *> col1_(6) square A1: = A1 - W1
325 *> end if
326 *> end if
327 *>
328 *> Case (B), when V1 is an identity matrix, i.e. IDENT == 'I',
329 *> less computations than in the Case (A)
330 *>
331 *> if( IDENT == 'I' ) then
332 *> if ( N > K ) then
333 *> (First Step - column block 2)
334 *> col2_(1) W2: = A2
335 *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
336 *> col2_(4) W2: = T * W2
337 *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
338 *> col2_(7) A2: = A2 - W2
339 *> else
340 *> (Second Step - column block 1)
341 *> col1_(1) W1: = A1
342 *> col1_(3) W1: = T * W1
343 *> col1_(4) B1: = - V2 * W1 = - B1 * W1
344 *> col1_(6) upper-triangular_of_(A1): = A1 - W1
345 *> end if
346 *> end if
347 *>
348 *> Combine these cases (A) and (B) together, this is the resulting
349 *> algorithm:
350 *>
351 *> if ( N > K ) then
352 *>
353 *> (First Step - column block 2)
354 *>
355 *> col2_(1) W2: = A2
356 *> if( IDENT != 'I' ) then
357 *> col2_(2) W2: = (V1**T) * W2
358 *> = (unit_lower_tr_of_(A1)**T) * W2
359 *> end if
360 *> col2_(3) W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2]
361 *> col2_(4) W2: = T * W2
362 *> col2_(5) B2: = B2 - V2 * W2 = B2 - B1 * W2
363 *> if( IDENT != 'I' ) then
364 *> col2_(6) W2: = V1 * W2 = unit_lower_tr_of_(A1) * W2
365 *> end if
366 *> col2_(7) A2: = A2 - W2
367 *>
368 *> else
369 *>
370 *> (Second Step - column block 1)
371 *>
372 *> col1_(1) W1: = A1
373 *> if( IDENT != 'I' ) then
374 *> col1_(2) W1: = (V1**T) * W1
375 *> = (unit_lower_tr_of_(A1)**T) * W1
376 *> end if
377 *> col1_(3) W1: = T * W1
378 *> col1_(4) B1: = - V2 * W1 = - B1 * W1
379 *> if( IDENT != 'I' ) then
380 *> col1_(5) square W1: = V1 * W1 = unit_lower_tr_of_(A1) * W1
381 *> col1_(6_a) below_diag_of_(A1): = - below_diag_of_(W1)
382 *> end if
383 *> col1_(6_b) up_tr_of_(A1): = up_tr_of_(A1) - up_tr_of_(W1)
384 *>
385 *> end if
386 *>
387 *> \endverbatim
388 *>
389 * =====================================================================
390  SUBROUTINE slarfb_gett( IDENT, M, N, K, T, LDT, A, LDA, B, LDB,
391  $ WORK, LDWORK )
392  IMPLICIT NONE
393 *
394 * -- LAPACK auxiliary routine --
395 * -- LAPACK is a software package provided by Univ. of Tennessee, --
396 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397 *
398 * .. Scalar Arguments ..
399  CHARACTER IDENT
400  INTEGER K, LDA, LDB, LDT, LDWORK, M, N
401 * ..
402 * .. Array Arguments ..
403  REAL A( LDA, * ), B( LDB, * ), T( LDT, * ),
404  $ work( ldwork, * )
405 * ..
406 *
407 * =====================================================================
408 *
409 * .. Parameters ..
410  REAL ONE, ZERO
411  parameter( one = 1.0e+0, zero = 0.0e+0 )
412 * ..
413 * .. Local Scalars ..
414  LOGICAL LNOTIDENT
415  INTEGER I, J
416 * ..
417 * .. EXTERNAL FUNCTIONS ..
418  LOGICAL LSAME
419  EXTERNAL lsame
420 * ..
421 * .. External Subroutines ..
422  EXTERNAL scopy, sgemm, strmm
423 * ..
424 * .. Executable Statements ..
425 *
426 * Quick return if possible
427 *
428  IF( m.LT.0 .OR. n.LE.0 .OR. k.EQ.0 .OR. k.GT.n )
429  $ RETURN
430 *
431  lnotident = .NOT.lsame( ident, 'I' )
432 *
433 * ------------------------------------------------------------------
434 *
435 * First Step. Computation of the Column Block 2:
436 *
437 * ( A2 ) := H * ( A2 )
438 * ( B2 ) ( B2 )
439 *
440 * ------------------------------------------------------------------
441 *
442  IF( n.GT.k ) THEN
443 *
444 * col2_(1) Compute W2: = A2. Therefore, copy A2 = A(1:K, K+1:N)
445 * into W2=WORK(1:K, 1:N-K) column-by-column.
446 *
447  DO j = 1, n-k
448  CALL scopy( k, a( 1, k+j ), 1, work( 1, j ), 1 )
449  END DO
450 
451  IF( lnotident ) THEN
452 *
453 * col2_(2) Compute W2: = (V1**T) * W2 = (A1**T) * W2,
454 * V1 is not an identy matrix, but unit lower-triangular
455 * V1 stored in A1 (diagonal ones are not stored).
456 *
457 *
458  CALL strmm( 'L', 'L', 'T', 'U', k, n-k, one, a, lda,
459  $ work, ldwork )
460  END IF
461 *
462 * col2_(3) Compute W2: = W2 + (V2**T) * B2 = W2 + (B1**T) * B2
463 * V2 stored in B1.
464 *
465  IF( m.GT.0 ) THEN
466  CALL sgemm( 'T', 'N', k, n-k, m, one, b, ldb,
467  $ b( 1, k+1 ), ldb, one, work, ldwork )
468  END IF
469 *
470 * col2_(4) Compute W2: = T * W2,
471 * T is upper-triangular.
472 *
473  CALL strmm( 'L', 'U', 'N', 'N', k, n-k, one, t, ldt,
474  $ work, ldwork )
475 *
476 * col2_(5) Compute B2: = B2 - V2 * W2 = B2 - B1 * W2,
477 * V2 stored in B1.
478 *
479  IF( m.GT.0 ) THEN
480  CALL sgemm( 'N', 'N', m, n-k, k, -one, b, ldb,
481  $ work, ldwork, one, b( 1, k+1 ), ldb )
482  END IF
483 *
484  IF( lnotident ) THEN
485 *
486 * col2_(6) Compute W2: = V1 * W2 = A1 * W2,
487 * V1 is not an identity matrix, but unit lower-triangular,
488 * V1 stored in A1 (diagonal ones are not stored).
489 *
490  CALL strmm( 'L', 'L', 'N', 'U', k, n-k, one, a, lda,
491  $ work, ldwork )
492  END IF
493 *
494 * col2_(7) Compute A2: = A2 - W2 =
495 * = A(1:K, K+1:N-K) - WORK(1:K, 1:N-K),
496 * column-by-column.
497 *
498  DO j = 1, n-k
499  DO i = 1, k
500  a( i, k+j ) = a( i, k+j ) - work( i, j )
501  END DO
502  END DO
503 *
504  END IF
505 *
506 * ------------------------------------------------------------------
507 *
508 * Second Step. Computation of the Column Block 1:
509 *
510 * ( A1 ) := H * ( A1 )
511 * ( B1 ) ( 0 )
512 *
513 * ------------------------------------------------------------------
514 *
515 * col1_(1) Compute W1: = A1. Copy the upper-triangular
516 * A1 = A(1:K, 1:K) into the upper-triangular
517 * W1 = WORK(1:K, 1:K) column-by-column.
518 *
519  DO j = 1, k
520  CALL scopy( j, a( 1, j ), 1, work( 1, j ), 1 )
521  END DO
522 *
523 * Set the subdiagonal elements of W1 to zero column-by-column.
524 *
525  DO j = 1, k - 1
526  DO i = j + 1, k
527  work( i, j ) = zero
528  END DO
529  END DO
530 *
531  IF( lnotident ) THEN
532 *
533 * col1_(2) Compute W1: = (V1**T) * W1 = (A1**T) * W1,
534 * V1 is not an identity matrix, but unit lower-triangular
535 * V1 stored in A1 (diagonal ones are not stored),
536 * W1 is upper-triangular with zeroes below the diagonal.
537 *
538  CALL strmm( 'L', 'L', 'T', 'U', k, k, one, a, lda,
539  $ work, ldwork )
540  END IF
541 *
542 * col1_(3) Compute W1: = T * W1,
543 * T is upper-triangular,
544 * W1 is upper-triangular with zeroes below the diagonal.
545 *
546  CALL strmm( 'L', 'U', 'N', 'N', k, k, one, t, ldt,
547  $ work, ldwork )
548 *
549 * col1_(4) Compute B1: = - V2 * W1 = - B1 * W1,
550 * V2 = B1, W1 is upper-triangular with zeroes below the diagonal.
551 *
552  IF( m.GT.0 ) THEN
553  CALL strmm( 'R', 'U', 'N', 'N', m, k, -one, work, ldwork,
554  $ b, ldb )
555  END IF
556 *
557  IF( lnotident ) THEN
558 *
559 * col1_(5) Compute W1: = V1 * W1 = A1 * W1,
560 * V1 is not an identity matrix, but unit lower-triangular
561 * V1 stored in A1 (diagonal ones are not stored),
562 * W1 is upper-triangular on input with zeroes below the diagonal,
563 * and square on output.
564 *
565  CALL strmm( 'L', 'L', 'N', 'U', k, k, one, a, lda,
566  $ work, ldwork )
567 *
568 * col1_(6) Compute A1: = A1 - W1 = A(1:K, 1:K) - WORK(1:K, 1:K)
569 * column-by-column. A1 is upper-triangular on input.
570 * If IDENT, A1 is square on output, and W1 is square,
571 * if NOT IDENT, A1 is upper-triangular on output,
572 * W1 is upper-triangular.
573 *
574 * col1_(6)_a Compute elements of A1 below the diagonal.
575 *
576  DO j = 1, k - 1
577  DO i = j + 1, k
578  a( i, j ) = - work( i, j )
579  END DO
580  END DO
581 *
582  END IF
583 *
584 * col1_(6)_b Compute elements of A1 on and above the diagonal.
585 *
586  DO j = 1, k
587  DO i = 1, j
588  a( i, j ) = a( i, j ) - work( i, j )
589  END DO
590  END DO
591 *
592  RETURN
593 *
594 * End of SLARFB_GETT
595 *
596  END
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:177
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine slarfb_gett(IDENT, M, N, K, T, LDT, A, LDA, B, LDB, WORK, LDWORK)
SLARFB_GETT
Definition: slarfb_gett.f:392