Annotation of rpl/lapack/lapack/dtplqt2.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DTPLQT2 computes a LQ factorization of a real or complex "triangular-pentagonal" matrix, which is composed of a triangular block and a pentagonal block, using the compact WY representation for Q.
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTPLQT2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtplqt2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtplqt2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtplqt2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTPLQT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
! 22: *
! 23: * .. Scalar Arguments ..
! 24: * INTEGER INFO, LDA, LDB, LDT, N, M, L
! 25: * ..
! 26: * .. Array Arguments ..
! 27: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * )
! 28: * ..
! 29: *
! 30: *
! 31: *> \par Purpose:
! 32: * =============
! 33: *>
! 34: *> \verbatim
! 35: *>
! 36: *> DTPLQT2 computes a LQ a factorization of a real "triangular-pentagonal"
! 37: *> matrix C, which is composed of a triangular block A and pentagonal block B,
! 38: *> using the compact WY representation for Q.
! 39: *> \endverbatim
! 40: *
! 41: * Arguments:
! 42: * ==========
! 43: *
! 44: *> \param[in] M
! 45: *> \verbatim
! 46: *> M is INTEGER
! 47: *> The total number of rows of the matrix B.
! 48: *> M >= 0.
! 49: *> \endverbatim
! 50: *>
! 51: *> \param[in] N
! 52: *> \verbatim
! 53: *> N is INTEGER
! 54: *> The number of columns of the matrix B, and the order of
! 55: *> the triangular matrix A.
! 56: *> N >= 0.
! 57: *> \endverbatim
! 58: *>
! 59: *> \param[in] L
! 60: *> \verbatim
! 61: *> L is INTEGER
! 62: *> The number of rows of the lower trapezoidal part of B.
! 63: *> MIN(M,N) >= L >= 0. See Further Details.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in,out] A
! 67: *> \verbatim
! 68: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 69: *> On entry, the lower triangular M-by-M matrix A.
! 70: *> On exit, the elements on and below the diagonal of the array
! 71: *> contain the lower triangular matrix L.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] LDA
! 75: *> \verbatim
! 76: *> LDA is INTEGER
! 77: *> The leading dimension of the array A. LDA >= max(1,N).
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in,out] B
! 81: *> \verbatim
! 82: *> B is DOUBLE PRECISION array, dimension (LDB,N)
! 83: *> On entry, the pentagonal M-by-N matrix B. The first N-L columns
! 84: *> are rectangular, and the last L columns are lower trapezoidal.
! 85: *> On exit, B contains the pentagonal matrix V. See Further Details.
! 86: *> \endverbatim
! 87: *>
! 88: *> \param[in] LDB
! 89: *> \verbatim
! 90: *> LDB is INTEGER
! 91: *> The leading dimension of the array B. LDB >= max(1,M).
! 92: *> \endverbatim
! 93: *>
! 94: *> \param[out] T
! 95: *> \verbatim
! 96: *> T is DOUBLE PRECISION array, dimension (LDT,M)
! 97: *> The N-by-N upper triangular factor T of the block reflector.
! 98: *> See Further Details.
! 99: *> \endverbatim
! 100: *>
! 101: *> \param[in] LDT
! 102: *> \verbatim
! 103: *> LDT is INTEGER
! 104: *> The leading dimension of the array T. LDT >= max(1,M)
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[out] INFO
! 108: *> \verbatim
! 109: *> INFO is INTEGER
! 110: *> = 0: successful exit
! 111: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 112: *> \endverbatim
! 113: *
! 114: * Authors:
! 115: * ========
! 116: *
! 117: *> \author Univ. of Tennessee
! 118: *> \author Univ. of California Berkeley
! 119: *> \author Univ. of Colorado Denver
! 120: *> \author NAG Ltd.
! 121: *
! 122: *> \date December 2016
! 123: *
! 124: *> \ingroup doubleOTHERcomputational
! 125: *
! 126: *> \par Further Details:
! 127: * =====================
! 128: *>
! 129: *> \verbatim
! 130: *>
! 131: *> The input matrix C is a M-by-(M+N) matrix
! 132: *>
! 133: *> C = [ A ][ B ]
! 134: *>
! 135: *>
! 136: *> where A is an lower triangular N-by-N matrix, and B is M-by-N pentagonal
! 137: *> matrix consisting of a M-by-(N-L) rectangular matrix B1 left of a M-by-L
! 138: *> upper trapezoidal matrix B2:
! 139: *>
! 140: *> B = [ B1 ][ B2 ]
! 141: *> [ B1 ] <- M-by-(N-L) rectangular
! 142: *> [ B2 ] <- M-by-L lower trapezoidal.
! 143: *>
! 144: *> The lower trapezoidal matrix B2 consists of the first L columns of a
! 145: *> N-by-N lower triangular matrix, where 0 <= L <= MIN(M,N). If L=0,
! 146: *> B is rectangular M-by-N; if M=L=N, B is lower triangular.
! 147: *>
! 148: *> The matrix W stores the elementary reflectors H(i) in the i-th row
! 149: *> above the diagonal (of A) in the M-by-(M+N) input matrix C
! 150: *>
! 151: *> C = [ A ][ B ]
! 152: *> [ A ] <- lower triangular N-by-N
! 153: *> [ B ] <- M-by-N pentagonal
! 154: *>
! 155: *> so that W can be represented as
! 156: *>
! 157: *> W = [ I ][ V ]
! 158: *> [ I ] <- identity, N-by-N
! 159: *> [ V ] <- M-by-N, same form as B.
! 160: *>
! 161: *> Thus, all of information needed for W is contained on exit in B, which
! 162: *> we call V above. Note that V has the same form as B; that is,
! 163: *>
! 164: *> W = [ V1 ][ V2 ]
! 165: *> [ V1 ] <- M-by-(N-L) rectangular
! 166: *> [ V2 ] <- M-by-L lower trapezoidal.
! 167: *>
! 168: *> The rows of V represent the vectors which define the H(i)'s.
! 169: *> The (M+N)-by-(M+N) block reflector H is then given by
! 170: *>
! 171: *> H = I - W**T * T * W
! 172: *>
! 173: *> where W^H is the conjugate transpose of W and T is the upper triangular
! 174: *> factor of the block reflector.
! 175: *> \endverbatim
! 176: *>
! 177: * =====================================================================
! 178: SUBROUTINE DTPLQT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
! 179: *
! 180: * -- LAPACK computational routine (version 3.7.0) --
! 181: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 182: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 183: * December 2016
! 184: *
! 185: * .. Scalar Arguments ..
! 186: INTEGER INFO, LDA, LDB, LDT, N, M, L
! 187: * ..
! 188: * .. Array Arguments ..
! 189: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * )
! 190: * ..
! 191: *
! 192: * =====================================================================
! 193: *
! 194: * .. Parameters ..
! 195: DOUBLE PRECISION ONE, ZERO
! 196: PARAMETER( ONE = 1.0, ZERO = 0.0 )
! 197: * ..
! 198: * .. Local Scalars ..
! 199: INTEGER I, J, P, MP, NP
! 200: DOUBLE PRECISION ALPHA
! 201: * ..
! 202: * .. External Subroutines ..
! 203: EXTERNAL DLARFG, DGEMV, DGER, DTRMV, XERBLA
! 204: * ..
! 205: * .. Intrinsic Functions ..
! 206: INTRINSIC MAX, MIN
! 207: * ..
! 208: * .. Executable Statements ..
! 209: *
! 210: * Test the input arguments
! 211: *
! 212: INFO = 0
! 213: IF( M.LT.0 ) THEN
! 214: INFO = -1
! 215: ELSE IF( N.LT.0 ) THEN
! 216: INFO = -2
! 217: ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
! 218: INFO = -3
! 219: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 220: INFO = -5
! 221: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
! 222: INFO = -7
! 223: ELSE IF( LDT.LT.MAX( 1, M ) ) THEN
! 224: INFO = -9
! 225: END IF
! 226: IF( INFO.NE.0 ) THEN
! 227: CALL XERBLA( 'DTPLQT2', -INFO )
! 228: RETURN
! 229: END IF
! 230: *
! 231: * Quick return if possible
! 232: *
! 233: IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
! 234: *
! 235: DO I = 1, M
! 236: *
! 237: * Generate elementary reflector H(I) to annihilate B(I,:)
! 238: *
! 239: P = N-L+MIN( L, I )
! 240: CALL DLARFG( P+1, A( I, I ), B( I, 1 ), LDB, T( 1, I ) )
! 241: IF( I.LT.M ) THEN
! 242: *
! 243: * W(M-I:1) := C(I+1:M,I:N) * C(I,I:N) [use W = T(M,:)]
! 244: *
! 245: DO J = 1, M-I
! 246: T( M, J ) = (A( I+J, I ))
! 247: END DO
! 248: CALL DGEMV( 'N', M-I, P, ONE, B( I+1, 1 ), LDB,
! 249: $ B( I, 1 ), LDB, ONE, T( M, 1 ), LDT )
! 250: *
! 251: * C(I+1:M,I:N) = C(I+1:M,I:N) + alpha * C(I,I:N)*W(M-1:1)^H
! 252: *
! 253: ALPHA = -(T( 1, I ))
! 254: DO J = 1, M-I
! 255: A( I+J, I ) = A( I+J, I ) + ALPHA*(T( M, J ))
! 256: END DO
! 257: CALL DGER( M-I, P, ALPHA, T( M, 1 ), LDT,
! 258: $ B( I, 1 ), LDB, B( I+1, 1 ), LDB )
! 259: END IF
! 260: END DO
! 261: *
! 262: DO I = 2, M
! 263: *
! 264: * T(I,1:I-1) := C(I:I-1,1:N) * (alpha * C(I,I:N)^H)
! 265: *
! 266: ALPHA = -T( 1, I )
! 267:
! 268: DO J = 1, I-1
! 269: T( I, J ) = ZERO
! 270: END DO
! 271: P = MIN( I-1, L )
! 272: NP = MIN( N-L+1, N )
! 273: MP = MIN( P+1, M )
! 274: *
! 275: * Triangular part of B2
! 276: *
! 277: DO J = 1, P
! 278: T( I, J ) = ALPHA*B( I, N-L+J )
! 279: END DO
! 280: CALL DTRMV( 'L', 'N', 'N', P, B( 1, NP ), LDB,
! 281: $ T( I, 1 ), LDT )
! 282: *
! 283: * Rectangular part of B2
! 284: *
! 285: CALL DGEMV( 'N', I-1-P, L, ALPHA, B( MP, NP ), LDB,
! 286: $ B( I, NP ), LDB, ZERO, T( I,MP ), LDT )
! 287: *
! 288: * B1
! 289: *
! 290: CALL DGEMV( 'N', I-1, N-L, ALPHA, B, LDB, B( I, 1 ), LDB,
! 291: $ ONE, T( I, 1 ), LDT )
! 292: *
! 293: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(I,1:I-1)
! 294: *
! 295: CALL DTRMV( 'L', 'T', 'N', I-1, T, LDT, T( I, 1 ), LDT )
! 296: *
! 297: * T(I,I) = tau(I)
! 298: *
! 299: T( I, I ) = T( 1, I )
! 300: T( 1, I ) = ZERO
! 301: END DO
! 302: DO I=1,M
! 303: DO J= I+1,M
! 304: T(I,J)=T(J,I)
! 305: T(J,I)= ZERO
! 306: END DO
! 307: END DO
! 308:
! 309: *
! 310: * End of DTPLQT2
! 311: *
! 312: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtplqt2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack