Annotation of rpl/lapack/lapack/zgelqt3.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZGELQT3 recursively computes a LQ factorization of a general real or complex matrix using the compact WY representation of Q.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGEQRT3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelqt3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelqt3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelqt3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * RECURSIVE SUBROUTINE ZGELQT3( M, N, A, LDA, T, LDT, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER INFO, LDA, M, N, LDT
! 25: * ..
! 26: * .. Array Arguments ..
! 27: * COMPLEX*16 A( LDA, * ), T( LDT, * )
! 28: * ..
! 29: *
! 30: *
! 31: *> \par Purpose:
! 32: * =============
! 33: *>
! 34: *> \verbatim
! 35: *>
! 36: *> DGELQT3 recursively computes a LQ factorization of a complex M-by-N
! 37: *> matrix A, using the compact WY representation of Q.
! 38: *>
! 39: *> Based on the algorithm of Elmroth and Gustavson,
! 40: *> IBM J. Res. Develop. Vol 44 No. 4 July 2000.
! 41: *> \endverbatim
! 42: *
! 43: * Arguments:
! 44: * ==========
! 45: *
! 46: *> \param[in] M
! 47: *> \verbatim
! 48: *> M is INTEGER
! 49: *> The number of rows of the matrix A. M =< N.
! 50: *> \endverbatim
! 51: *>
! 52: *> \param[in] N
! 53: *> \verbatim
! 54: *> N is INTEGER
! 55: *> The number of columns of the matrix A. N >= 0.
! 56: *> \endverbatim
! 57: *>
! 58: *> \param[in,out] A
! 59: *> \verbatim
! 60: *> A is COMPLEX*16 array, dimension (LDA,N)
! 61: *> On entry, the real M-by-N matrix A. On exit, the elements on and
! 62: *> below the diagonal contain the N-by-N lower triangular matrix L; the
! 63: *> elements above the diagonal are the rows of V. See below for
! 64: *> further details.
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in] LDA
! 68: *> \verbatim
! 69: *> LDA is INTEGER
! 70: *> The leading dimension of the array A. LDA >= max(1,M).
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[out] T
! 74: *> \verbatim
! 75: *> T is COMPLEX*16 array, dimension (LDT,N)
! 76: *> The N-by-N upper triangular factor of the block reflector.
! 77: *> The elements on and above the diagonal contain the block
! 78: *> reflector T; the elements below the diagonal are not used.
! 79: *> See below for further details.
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] LDT
! 83: *> \verbatim
! 84: *> LDT is INTEGER
! 85: *> The leading dimension of the array T. LDT >= max(1,N).
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[out] INFO
! 89: *> \verbatim
! 90: *> INFO is INTEGER
! 91: *> = 0: successful exit
! 92: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 93: *> \endverbatim
! 94: *
! 95: * Authors:
! 96: * ========
! 97: *
! 98: *> \author Univ. of Tennessee
! 99: *> \author Univ. of California Berkeley
! 100: *> \author Univ. of Colorado Denver
! 101: *> \author NAG Ltd.
! 102: *
! 103: *> \date December 2016
! 104: *
! 105: *> \ingroup doubleGEcomputational
! 106: *
! 107: *> \par Further Details:
! 108: * =====================
! 109: *>
! 110: *> \verbatim
! 111: *>
! 112: *> The matrix V stores the elementary reflectors H(i) in the i-th column
! 113: *> below the diagonal. For example, if M=5 and N=3, the matrix V is
! 114: *>
! 115: *> V = ( 1 v1 v1 v1 v1 )
! 116: *> ( 1 v2 v2 v2 )
! 117: *> ( 1 v3 v3 v3 )
! 118: *>
! 119: *>
! 120: *> where the vi's represent the vectors which define H(i), which are returned
! 121: *> in the matrix A. The 1's along the diagonal of V are not stored in A. The
! 122: *> block reflector H is then given by
! 123: *>
! 124: *> H = I - V * T * V**T
! 125: *>
! 126: *> where V**T is the transpose of V.
! 127: *>
! 128: *> For details of the algorithm, see Elmroth and Gustavson (cited above).
! 129: *> \endverbatim
! 130: *>
! 131: * =====================================================================
! 132: RECURSIVE SUBROUTINE ZGELQT3( M, N, A, LDA, T, LDT, INFO )
! 133: *
! 134: * -- LAPACK computational routine (version 3.7.0) --
! 135: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 136: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 137: * December 2016
! 138: *
! 139: * .. Scalar Arguments ..
! 140: INTEGER INFO, LDA, M, N, LDT
! 141: * ..
! 142: * .. Array Arguments ..
! 143: COMPLEX*16 A( LDA, * ), T( LDT, * )
! 144: * ..
! 145: *
! 146: * =====================================================================
! 147: *
! 148: * .. Parameters ..
! 149: COMPLEX*16 ONE, ZERO
! 150: PARAMETER ( ONE = (1.0D+00,0.0D+00) )
! 151: PARAMETER ( ZERO = (0.0D+00,0.0D+00))
! 152: * ..
! 153: * .. Local Scalars ..
! 154: INTEGER I, I1, J, J1, M1, M2, N1, N2, IINFO
! 155: * ..
! 156: * .. External Subroutines ..
! 157: EXTERNAL ZLARFG, ZTRMM, ZGEMM, XERBLA
! 158: * ..
! 159: * .. Executable Statements ..
! 160: *
! 161: INFO = 0
! 162: IF( M .LT. 0 ) THEN
! 163: INFO = -1
! 164: ELSE IF( N .LT. M ) THEN
! 165: INFO = -2
! 166: ELSE IF( LDA .LT. MAX( 1, M ) ) THEN
! 167: INFO = -4
! 168: ELSE IF( LDT .LT. MAX( 1, M ) ) THEN
! 169: INFO = -6
! 170: END IF
! 171: IF( INFO.NE.0 ) THEN
! 172: CALL XERBLA( 'ZGELQT3', -INFO )
! 173: RETURN
! 174: END IF
! 175: *
! 176: IF( M.EQ.1 ) THEN
! 177: *
! 178: * Compute Householder transform when N=1
! 179: *
! 180: CALL ZLARFG( N, A, A( 1, MIN( 2, N ) ), LDA, T )
! 181: T(1,1)=CONJG(T(1,1))
! 182: *
! 183: ELSE
! 184: *
! 185: * Otherwise, split A into blocks...
! 186: *
! 187: M1 = M/2
! 188: M2 = M-M1
! 189: I1 = MIN( M1+1, M )
! 190: J1 = MIN( M+1, N )
! 191: *
! 192: * Compute A(1:M1,1:N) <- (Y1,R1,T1), where Q1 = I - Y1 T1 Y1^H
! 193: *
! 194: CALL ZGELQT3( M1, N, A, LDA, T, LDT, IINFO )
! 195: *
! 196: * Compute A(J1:M,1:N) = A(J1:M,1:N) Q1^H [workspace: T(1:N1,J1:N)]
! 197: *
! 198: DO I=1,M2
! 199: DO J=1,M1
! 200: T( I+M1, J ) = A( I+M1, J )
! 201: END DO
! 202: END DO
! 203: CALL ZTRMM( 'R', 'U', 'C', 'U', M2, M1, ONE,
! 204: & A, LDA, T( I1, 1 ), LDT )
! 205: *
! 206: CALL ZGEMM( 'N', 'C', M2, M1, N-M1, ONE, A( I1, I1 ), LDA,
! 207: & A( 1, I1 ), LDA, ONE, T( I1, 1 ), LDT)
! 208: *
! 209: CALL ZTRMM( 'R', 'U', 'N', 'N', M2, M1, ONE,
! 210: & T, LDT, T( I1, 1 ), LDT )
! 211: *
! 212: CALL ZGEMM( 'N', 'N', M2, N-M1, M1, -ONE, T( I1, 1 ), LDT,
! 213: & A( 1, I1 ), LDA, ONE, A( I1, I1 ), LDA )
! 214: *
! 215: CALL ZTRMM( 'R', 'U', 'N', 'U', M2, M1 , ONE,
! 216: & A, LDA, T( I1, 1 ), LDT )
! 217: *
! 218: DO I=1,M2
! 219: DO J=1,M1
! 220: A( I+M1, J ) = A( I+M1, J ) - T( I+M1, J )
! 221: T( I+M1, J )= ZERO
! 222: END DO
! 223: END DO
! 224: *
! 225: * Compute A(J1:M,J1:N) <- (Y2,R2,T2) where Q2 = I - Y2 T2 Y2^H
! 226: *
! 227: CALL ZGELQT3( M2, N-M1, A( I1, I1 ), LDA,
! 228: & T( I1, I1 ), LDT, IINFO )
! 229: *
! 230: * Compute T3 = T(J1:N1,1:N) = -T1 Y1^H Y2 T2
! 231: *
! 232: DO I=1,M2
! 233: DO J=1,M1
! 234: T( J, I+M1 ) = (A( J, I+M1 ))
! 235: END DO
! 236: END DO
! 237: *
! 238: CALL ZTRMM( 'R', 'U', 'C', 'U', M1, M2, ONE,
! 239: & A( I1, I1 ), LDA, T( 1, I1 ), LDT )
! 240: *
! 241: CALL ZGEMM( 'N', 'C', M1, M2, N-M, ONE, A( 1, J1 ), LDA,
! 242: & A( I1, J1 ), LDA, ONE, T( 1, I1 ), LDT )
! 243: *
! 244: CALL ZTRMM( 'L', 'U', 'N', 'N', M1, M2, -ONE, T, LDT,
! 245: & T( 1, I1 ), LDT )
! 246: *
! 247: CALL ZTRMM( 'R', 'U', 'N', 'N', M1, M2, ONE,
! 248: & T( I1, I1 ), LDT, T( 1, I1 ), LDT )
! 249: *
! 250: *
! 251: *
! 252: * Y = (Y1,Y2); L = [ L1 0 ]; T = [T1 T3]
! 253: * [ A(1:N1,J1:N) L2 ] [ 0 T2]
! 254: *
! 255: END IF
! 256: *
! 257: RETURN
! 258: *
! 259: * End of ZGELQT3
! 260: *
! 261: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelqt3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack