LAPACK  3.9.0
LAPACK: Linear Algebra PACKage

◆ zunhr_col()

subroutine zunhr_col ( integer  M,
integer  N,
integer  NB,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldt, * )  T,
integer  LDT,
complex*16, dimension( * )  D,
integer  INFO 
)

ZUNHR_COL

Download ZUNHR_COL + dependencies [TGZ] [ZIP] [TXT] \par Purpose: @verbatim ZUNHR_COL takes an M-by-N complex matrix Q_in with orthonormal columns as input, stored in A, and performs Householder Reconstruction (HR), i.e. reconstructs Householder vectors V(i) implicitly representing another M-by-N matrix Q_out, with the property that Q_in = Q_out*S, where S is an N-by-N diagonal matrix with diagonal entries equal to +1 or -1. The Householder vectors (columns V(i) of V) are stored in A on output, and the diagonal entries of S are stored in D. Block reflectors are also returned in T (same output format as ZGEQRT). \endverbatim @param [in] M @verbatim M is INTEGER The number of rows of the matrix A. M >= 0. \endverbatim @param [in] N @verbatim N is INTEGER The number of columns of the matrix A. M >= N >= 0. \endverbatim @param [in] NB @verbatim NB is INTEGER The column block size to be used in the reconstruction of Householder column vector blocks in the array A and corresponding block reflectors in the array T. NB >= 1. (Note that if NB > N, then N is used instead of NB as the column block size.) \endverbatim @param [in,out] A @verbatim A is COMPLEX*16 array, dimension (LDA,N) On entry: The array A contains an M-by-N orthonormal matrix Q_in, i.e the columns of A are orthogonal unit vectors. On exit: The elements below the diagonal of A represent the unit lower-trapezoidal matrix V of Householder column vectors V(i). The unit diagonal entries of V are not stored (same format as the output below the diagonal in A from ZGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out. The elements above the diagonal contain the factor U of the "modified" LU-decomposition: Q_in - ( S ) = V * U ( 0 ) where 0 is a (M-N)-by-(M-N) zero matrix. \endverbatim @param [in] LDA @verbatim LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M). \endverbatim @param [out] T @verbatim T is COMPLEX*16 array, dimension (LDT, N) Let NOCB = Number_of_output_col_blocks = CEIL(N/NB) On exit, T(1:NB, 1:N) contains NOCB upper-triangular block reflectors used to define Q_out stored in compact form as a sequence of upper-triangular NB-by-NB column blocks (same format as the output T in ZGEQRT). The matrix T and the matrix V stored on output in A implicitly define Q_out. NOTE: The lower triangles below the upper-triangular blcoks will be filled with zeros. See Further Details. \endverbatim @param [in] LDT @verbatim LDT is INTEGER The leading dimension of the array T. LDT >= max(1,min(NB,N)). \endverbatim @param [out] D @verbatim D is COMPLEX*16 array, dimension min(M,N). The elements can be only plus or minus one. D(i) is constructed as D(i) = -SIGN(Q_in_i(i,i)), where 1 <= i <= min(M,N), and Q_in_i is Q_in after performing i-1 steps of “modified” Gaussian elimination. See Further Details. \endverbatim @param [out] INFO @verbatim INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value \endverbatim \par Further Details: @verbatim The computed M-by-M unitary factor Q_out is defined implicitly as a product of unitary matrices Q_out(i). Each Q_out(i) is stored in the compact WY-representation format in the corresponding blocks of matrices V (stored in A) and T. The M-by-N unit lower-trapezoidal matrix V stored in the M-by-N matrix A contains the column vectors V(i) in NB-size column blocks VB(j). For example, VB(1) contains the columns V(1), V(2), ... V(NB). NOTE: The unit entries on the diagonal of Y are not stored in A. The number of column blocks is NOCB = Number_of_output_col_blocks = CEIL(N/NB) where each block is of order NB except for the last block, which is of order LAST_NB = N - (NOCB-1)*NB. For example, if M=6, N=5 and NB=2, the matrix V is V = ( VB(1), VB(2), VB(3) ) = = ( 1 ) ( v21 1 ) ( v31 v32 1 ) ( v41 v42 v43 1 ) ( v51 v52 v53 v54 1 ) ( v61 v62 v63 v54 v65 ) For each of the column blocks VB(i), an upper-triangular block reflector TB(i) is computed. These blocks are stored as a sequence of upper-triangular column blocks in the NB-by-N matrix T. The size of each TB(i) block is NB-by-NB, except for the last block, whose size is LAST_NB-by-LAST_NB. For example, if M=6, N=5 and NB=2, the matrix T is T = ( TB(1), TB(2), TB(3) ) = = ( t11 t12 t13 t14 t15 ) ( t22 t24 ) The M-by-M factor Q_out is given as a product of NOCB unitary M-by-M matrices Q_out(i). Q_out = Q_out(1) * Q_out(2) * ... * Q_out(NOCB), where each matrix Q_out(i) is given by the WY-representation using corresponding blocks from the matrices V and T: Q_out(i) = I - VB(i) * TB(i) * (VB(i))**T, where I is the identity matrix. Here is the formula with matrix dimensions: Q(i){M-by-M} = I{M-by-M} - VB(i){M-by-INB} * TB(i){INB-by-INB} * (VB(i))**T {INB-by-M}, where INB = NB, except for the last block NOCB for which INB=LAST_NB. ===== NOTE: ===== If Q_in is the result of doing a QR factorization B = Q_in * R_in, then: B = (Q_out*S) * R_in = Q_out * (S * R_in) = O_out * R_out. So if one wants to interpret Q_out as the result of the QR factorization of B, then corresponding R_out should be obtained by R_out = S * R_in, i.e. some rows of R_in should be multiplied by -1. For the details of the algorithm, see [1]. [1] "Reconstructing Householder vectors from tall-skinny QR", G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen, E. Solomonik, J. Parallel Distrib. Comput., vol. 85, pp. 3-31, 2015. \endverbatim \author Univ. of Tennessee \author Univ. of California Berkeley \author Univ. of Colorado Denver \author NAG Ltd. \date November 2019 \par Contributors: @verbatim November 2019, Igor Kozachenko, Computer Science Division, University of California, Berkeley \endverbatim

Definition at line 260 of file zunhr_col.f.

260  IMPLICIT NONE
261 *
262 * -- LAPACK computational routine (version 3.9.0) --
263 * -- LAPACK is a software package provided by Univ. of Tennessee, --
264 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265 * November 2019
266 *
267 * .. Scalar Arguments ..
268  INTEGER INFO, LDA, LDT, M, N, NB
269 * ..
270 * .. Array Arguments ..
271  COMPLEX*16 A( LDA, * ), D( * ), T( LDT, * )
272 * ..
273 *
274 * =====================================================================
275 *
276 * .. Parameters ..
277  COMPLEX*16 CONE, CZERO
278  parameter( cone = ( 1.0d+0, 0.0d+0 ),
279  $ czero = ( 0.0d+0, 0.0d+0 ) )
280 * ..
281 * .. Local Scalars ..
282  INTEGER I, IINFO, J, JB, JBTEMP1, JBTEMP2, JNB,
283  $ NPLUSONE
284 * ..
285 * .. External Subroutines ..
286  EXTERNAL zcopy, zlaunhr_col_getrfnp, zscal, ztrsm,
287  $ xerbla
288 * ..
289 * .. Intrinsic Functions ..
290  INTRINSIC max, min
291 * ..
292 * .. Executable Statements ..
293 *
294 * Test the input parameters
295 *
296  info = 0
297  IF( m.LT.0 ) THEN
298  info = -1
299  ELSE IF( n.LT.0 .OR. n.GT.m ) THEN
300  info = -2
301  ELSE IF( nb.LT.1 ) THEN
302  info = -3
303  ELSE IF( lda.LT.max( 1, m ) ) THEN
304  info = -5
305  ELSE IF( ldt.LT.max( 1, min( nb, n ) ) ) THEN
306  info = -7
307  END IF
308 *
309 * Handle error in the input parameters.
310 *
311  IF( info.NE.0 ) THEN
312  CALL xerbla( 'ZUNHR_COL', -info )
313  RETURN
314  END IF
315 *
316 * Quick return if possible
317 *
318  IF( min( m, n ).EQ.0 ) THEN
319  RETURN
320  END IF
321 *
322 * On input, the M-by-N matrix A contains the unitary
323 * M-by-N matrix Q_in.
324 *
325 * (1) Compute the unit lower-trapezoidal V (ones on the diagonal
326 * are not stored) by performing the "modified" LU-decomposition.
327 *
328 * Q_in - ( S ) = V * U = ( V1 ) * U,
329 * ( 0 ) ( V2 )
330 *
331 * where 0 is an (M-N)-by-N zero matrix.
332 *
333 * (1-1) Factor V1 and U.
334 
335  CALL zlaunhr_col_getrfnp( n, n, a, lda, d, iinfo )
336 *
337 * (1-2) Solve for V2.
338 *
339  IF( m.GT.n ) THEN
340  CALL ztrsm( 'R', 'U', 'N', 'N', m-n, n, cone, a, lda,
341  $ a( n+1, 1 ), lda )
342  END IF
343 *
344 * (2) Reconstruct the block reflector T stored in T(1:NB, 1:N)
345 * as a sequence of upper-triangular blocks with NB-size column
346 * blocking.
347 *
348 * Loop over the column blocks of size NB of the array A(1:M,1:N)
349 * and the array T(1:NB,1:N), JB is the column index of a column
350 * block, JNB is the column block size at each step JB.
351 *
352  nplusone = n + 1
353  DO jb = 1, n, nb
354 *
355 * (2-0) Determine the column block size JNB.
356 *
357  jnb = min( nplusone-jb, nb )
358 *
359 * (2-1) Copy the upper-triangular part of the current JNB-by-JNB
360 * diagonal block U(JB) (of the N-by-N matrix U) stored
361 * in A(JB:JB+JNB-1,JB:JB+JNB-1) into the upper-triangular part
362 * of the current JNB-by-JNB block T(1:JNB,JB:JB+JNB-1)
363 * column-by-column, total JNB*(JNB+1)/2 elements.
364 *
365  jbtemp1 = jb - 1
366  DO j = jb, jb+jnb-1
367  CALL zcopy( j-jbtemp1, a( jb, j ), 1, t( 1, j ), 1 )
368  END DO
369 *
370 * (2-2) Perform on the upper-triangular part of the current
371 * JNB-by-JNB diagonal block U(JB) (of the N-by-N matrix U) stored
372 * in T(1:JNB,JB:JB+JNB-1) the following operation in place:
373 * (-1)*U(JB)*S(JB), i.e the result will be stored in the upper-
374 * triangular part of T(1:JNB,JB:JB+JNB-1). This multiplication
375 * of the JNB-by-JNB diagonal block U(JB) by the JNB-by-JNB
376 * diagonal block S(JB) of the N-by-N sign matrix S from the
377 * right means changing the sign of each J-th column of the block
378 * U(JB) according to the sign of the diagonal element of the block
379 * S(JB), i.e. S(J,J) that is stored in the array element D(J).
380 *
381  DO j = jb, jb+jnb-1
382  IF( d( j ).EQ.cone ) THEN
383  CALL zscal( j-jbtemp1, -cone, t( 1, j ), 1 )
384  END IF
385  END DO
386 *
387 * (2-3) Perform the triangular solve for the current block
388 * matrix X(JB):
389 *
390 * X(JB) * (A(JB)**T) = B(JB), where:
391 *
392 * A(JB)**T is a JNB-by-JNB unit upper-triangular
393 * coefficient block, and A(JB)=V1(JB), which
394 * is a JNB-by-JNB unit lower-triangular block
395 * stored in A(JB:JB+JNB-1,JB:JB+JNB-1).
396 * The N-by-N matrix V1 is the upper part
397 * of the M-by-N lower-trapezoidal matrix V
398 * stored in A(1:M,1:N);
399 *
400 * B(JB) is a JNB-by-JNB upper-triangular right-hand
401 * side block, B(JB) = (-1)*U(JB)*S(JB), and
402 * B(JB) is stored in T(1:JNB,JB:JB+JNB-1);
403 *
404 * X(JB) is a JNB-by-JNB upper-triangular solution
405 * block, X(JB) is the upper-triangular block
406 * reflector T(JB), and X(JB) is stored
407 * in T(1:JNB,JB:JB+JNB-1).
408 *
409 * In other words, we perform the triangular solve for the
410 * upper-triangular block T(JB):
411 *
412 * T(JB) * (V1(JB)**T) = (-1)*U(JB)*S(JB).
413 *
414 * Even though the blocks X(JB) and B(JB) are upper-
415 * triangular, the routine ZTRSM will access all JNB**2
416 * elements of the square T(1:JNB,JB:JB+JNB-1). Therefore,
417 * we need to set to zero the elements of the block
418 * T(1:JNB,JB:JB+JNB-1) below the diagonal before the call
419 * to ZTRSM.
420 *
421 * (2-3a) Set the elements to zero.
422 *
423  jbtemp2 = jb - 2
424  DO j = jb, jb+jnb-2
425  DO i = j-jbtemp2, nb
426  t( i, j ) = czero
427  END DO
428  END DO
429 *
430 * (2-3b) Perform the triangular solve.
431 *
432  CALL ztrsm( 'R', 'L', 'C', 'U', jnb, jnb, cone,
433  $ a( jb, jb ), lda, t( 1, jb ), ldt )
434 *
435  END DO
436 *
437  RETURN
438 *
439 * End of ZUNHR_COL
440 *
Here is the call graph for this function:
Here is the caller graph for this function:
ztrsm
subroutine ztrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
ZTRSM
Definition: ztrsm.f:182
zcopy
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
Definition: zcopy.f:83
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zlaunhr_col_getrfnp
subroutine zlaunhr_col_getrfnp(M, N, A, LDA, D, INFO)
ZLAUNHR_COL_GETRFNP
Definition: zlaunhr_col_getrfnp.f:148
zscal
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80