dorhr_col()

subroutine dorhr_col	(	integer	M,
		integer	N,
		integer	NB,
		double precision, dimension( lda, * )	A,
		integer	LDA,
		double precision, dimension( ldt, * )	T,
		integer	LDT,
		double precision, dimension( * )	D,
		integer	INFO
	)

DORHR_COL

Download DORHR_COL + dependencies [TGZ] [ZIP] [TXT] \par Purpose: @verbatim DORHR_COL takes an M-by-N real 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 DGEQRT). \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 DOUBLE PRECISION 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 DGEQRT). 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 DOUBLE PRECISION 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 DGEQRT). 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 DOUBLE PRECISION 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 orthogonal factor Q_out is defined implicitly as a product of orthogonal 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 orthogonal 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 dorhr_col.f.

      IMPLICIT NONE
*
*  -- LAPACK computational routine (version 3.9.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2019
*
*     .. Scalar Arguments ..
      INTEGER           INFO, LDA, LDT, M, N, NB
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION  A( LDA, * ), D( * ), T( LDT, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      parameter( one = 1.0d+0, zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, IINFO, J, JB, JBTEMP1, JBTEMP2, JNB,
     $                   NPLUSONE
*     ..
*     .. External Subroutines ..
      EXTERNAL           dcopy, dlaorhr_col_getrfnp, dscal, dtrsm,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min
*     ..
*     .. Executable Statements ..
*
*     Test the input parameters
*
      info = 0
      IF( m.LT.0 ) THEN
         info = -1
      ELSE IF( n.LT.0 .OR. n.GT.m ) THEN
         info = -2
      ELSE IF( nb.LT.1 ) THEN
         info = -3
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -5
      ELSE IF( ldt.LT.max( 1, min( nb, n ) ) ) THEN
         info = -7
      END IF
*
*     Handle error in the input parameters.
*
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DORHR_COL', -info )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( min( m, n ).EQ.0 ) THEN
         RETURN
      END IF
*
*     On input, the M-by-N matrix A contains the orthogonal
*     M-by-N matrix Q_in.
*
*     (1) Compute the unit lower-trapezoidal V (ones on the diagonal
*     are not stored) by performing the "modified" LU-decomposition.
*
*     Q_in - ( S ) = V * U = ( V1 ) * U,
*            ( 0 )           ( V2 )
*
*     where 0 is an (M-N)-by-N zero matrix.
*
*     (1-1) Factor V1 and U.
 
      CALL dlaorhr_col_getrfnp( n, n, a, lda, d, iinfo )
*
*     (1-2) Solve for V2.
*
      IF( m.GT.n ) THEN
         CALL dtrsm( 'R', 'U', 'N', 'N', m-n, n, one, a, lda,
     $               a( n+1, 1 ), lda )
      END IF
*
*     (2) Reconstruct the block reflector T stored in T(1:NB, 1:N)
*     as a sequence of upper-triangular blocks with NB-size column
*     blocking.
*
*     Loop over the column blocks of size NB of the array A(1:M,1:N)
*     and the array T(1:NB,1:N), JB is the column index of a column
*     block, JNB is the column block size at each step JB.
*
      nplusone = n + 1
      DO jb = 1, n, nb
*
*        (2-0) Determine the column block size JNB.
*
         jnb = min( nplusone-jb, nb )
*
*        (2-1) Copy the upper-triangular part of the current JNB-by-JNB
*        diagonal block U(JB) (of the N-by-N matrix U) stored
*        in A(JB:JB+JNB-1,JB:JB+JNB-1) into the upper-triangular part
*        of the current JNB-by-JNB block T(1:JNB,JB:JB+JNB-1)
*        column-by-column, total JNB*(JNB+1)/2 elements.
*
         jbtemp1 = jb - 1
         DO j = jb, jb+jnb-1
            CALL dcopy( j-jbtemp1, a( jb, j ), 1, t( 1, j ), 1 )
         END DO
*
*        (2-2) Perform on the upper-triangular part of the current
*        JNB-by-JNB diagonal block U(JB) (of the N-by-N matrix U) stored
*        in T(1:JNB,JB:JB+JNB-1) the following operation in place:
*        (-1)*U(JB)*S(JB), i.e the result will be stored in the upper-
*        triangular part of T(1:JNB,JB:JB+JNB-1). This multiplication
*        of the JNB-by-JNB diagonal block U(JB) by the JNB-by-JNB
*        diagonal block S(JB) of the N-by-N sign matrix S from the
*        right means changing the sign of each J-th column of the block
*        U(JB) according to the sign of the diagonal element of the block
*        S(JB), i.e. S(J,J) that is stored in the array element D(J).
*
         DO j = jb, jb+jnb-1
            IF( d( j ).EQ.one ) THEN
               CALL dscal( j-jbtemp1, -one, t( 1, j ), 1 )
            END IF
         END DO
*
*        (2-3) Perform the triangular solve for the current block
*        matrix X(JB):
*
*               X(JB) * (A(JB)**T) = B(JB), where:
*
*               A(JB)**T  is a JNB-by-JNB unit upper-triangular
*                         coefficient block, and A(JB)=V1(JB), which
*                         is a JNB-by-JNB unit lower-triangular block
*                         stored in A(JB:JB+JNB-1,JB:JB+JNB-1).
*                         The N-by-N matrix V1 is the upper part
*                         of the M-by-N lower-trapezoidal matrix V
*                         stored in A(1:M,1:N);
*
*               B(JB)     is a JNB-by-JNB  upper-triangular right-hand
*                         side block, B(JB) = (-1)*U(JB)*S(JB), and
*                         B(JB) is stored in T(1:JNB,JB:JB+JNB-1);
*
*               X(JB)     is a JNB-by-JNB upper-triangular solution
*                         block, X(JB) is the upper-triangular block
*                         reflector T(JB), and X(JB) is stored
*                         in T(1:JNB,JB:JB+JNB-1).
*
*             In other words, we perform the triangular solve for the
*             upper-triangular block T(JB):
*
*               T(JB) * (V1(JB)**T) = (-1)*U(JB)*S(JB).
*
*             Even though the blocks X(JB) and B(JB) are upper-
*             triangular, the routine DTRSM will access all JNB**2
*             elements of the square T(1:JNB,JB:JB+JNB-1). Therefore,
*             we need to set to zero the elements of the block
*             T(1:JNB,JB:JB+JNB-1) below the diagonal before the call
*             to DTRSM.
*
*        (2-3a) Set the elements to zero.
*
         jbtemp2 = jb - 2
         DO j = jb, jb+jnb-2
            DO i = j-jbtemp2, nb
               t( i, j ) = zero
            END DO
         END DO
*
*        (2-3b) Perform the triangular solve.
*
         CALL dtrsm( 'R', 'L', 'T', 'U', jnb, jnb, one,
     $               a( jb, jb ), lda, t( 1, jb ), ldt )
*
      END DO
*
      RETURN
*
*     End of DORHR_COL
*

Here is the call graph for this function:

Here is the caller graph for this function: