Annotation of rpl/lapack/lapack/dtpqrt2.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DTPQRT2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DTPQRT2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtpqrt2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtpqrt2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtpqrt2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DTPQRT2( 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: *> DTPQRT2 computes a QR 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 upper 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 upper triangular N-by-N matrix A.
! 70: *> On exit, the elements on and above the diagonal of the array
! 71: *> contain the upper triangular matrix R.
! 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 M-L rows
! 84: *> are rectangular, and the last L rows are upper 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,N)
! 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,N)
! 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 April 2012
! 123: *
! 124: *> \ingroup doubleOTHERcomputational
! 125: *
! 126: *> \par Further Details:
! 127: * =====================
! 128: *>
! 129: *> \verbatim
! 130: *>
! 131: *> The input matrix C is a (N+M)-by-N matrix
! 132: *>
! 133: *> C = [ A ]
! 134: *> [ B ]
! 135: *>
! 136: *> where A is an upper triangular N-by-N matrix, and B is M-by-N pentagonal
! 137: *> matrix consisting of a (M-L)-by-N rectangular matrix B1 on top of a L-by-N
! 138: *> upper trapezoidal matrix B2:
! 139: *>
! 140: *> B = [ B1 ] <- (M-L)-by-N rectangular
! 141: *> [ B2 ] <- L-by-N upper trapezoidal.
! 142: *>
! 143: *> The upper trapezoidal matrix B2 consists of the first L rows of a
! 144: *> N-by-N upper triangular matrix, where 0 <= L <= MIN(M,N). If L=0,
! 145: *> B is rectangular M-by-N; if M=L=N, B is upper triangular.
! 146: *>
! 147: *> The matrix W stores the elementary reflectors H(i) in the i-th column
! 148: *> below the diagonal (of A) in the (N+M)-by-N input matrix C
! 149: *>
! 150: *> C = [ A ] <- upper triangular N-by-N
! 151: *> [ B ] <- M-by-N pentagonal
! 152: *>
! 153: *> so that W can be represented as
! 154: *>
! 155: *> W = [ I ] <- identity, N-by-N
! 156: *> [ V ] <- M-by-N, same form as B.
! 157: *>
! 158: *> Thus, all of information needed for W is contained on exit in B, which
! 159: *> we call V above. Note that V has the same form as B; that is,
! 160: *>
! 161: *> V = [ V1 ] <- (M-L)-by-N rectangular
! 162: *> [ V2 ] <- L-by-N upper trapezoidal.
! 163: *>
! 164: *> The columns of V represent the vectors which define the H(i)'s.
! 165: *> The (M+N)-by-(M+N) block reflector H is then given by
! 166: *>
! 167: *> H = I - W * T * W**T
! 168: *>
! 169: *> where W^H is the conjugate transpose of W and T is the upper triangular
! 170: *> factor of the block reflector.
! 171: *> \endverbatim
! 172: *>
! 173: * =====================================================================
! 174: SUBROUTINE DTPQRT2( M, N, L, A, LDA, B, LDB, T, LDT, INFO )
! 175: *
! 176: * -- LAPACK computational routine (version 3.4.1) --
! 177: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 178: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 179: * April 2012
! 180: *
! 181: * .. Scalar Arguments ..
! 182: INTEGER INFO, LDA, LDB, LDT, N, M, L
! 183: * ..
! 184: * .. Array Arguments ..
! 185: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * )
! 186: * ..
! 187: *
! 188: * =====================================================================
! 189: *
! 190: * .. Parameters ..
! 191: DOUBLE PRECISION ONE, ZERO
! 192: PARAMETER( ONE = 1.0, ZERO = 0.0 )
! 193: * ..
! 194: * .. Local Scalars ..
! 195: INTEGER I, J, P, MP, NP
! 196: DOUBLE PRECISION ALPHA
! 197: * ..
! 198: * .. External Subroutines ..
! 199: EXTERNAL DLARFG, DGEMV, DGER, DTRMV, XERBLA
! 200: * ..
! 201: * .. Intrinsic Functions ..
! 202: INTRINSIC MAX, MIN
! 203: * ..
! 204: * .. Executable Statements ..
! 205: *
! 206: * Test the input arguments
! 207: *
! 208: INFO = 0
! 209: IF( M.LT.0 ) THEN
! 210: INFO = -1
! 211: ELSE IF( N.LT.0 ) THEN
! 212: INFO = -2
! 213: ELSE IF( L.LT.0 .OR. L.GT.MIN(M,N) ) THEN
! 214: INFO = -3
! 215: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 216: INFO = -5
! 217: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
! 218: INFO = -7
! 219: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
! 220: INFO = -9
! 221: END IF
! 222: IF( INFO.NE.0 ) THEN
! 223: CALL XERBLA( 'DTPQRT2', -INFO )
! 224: RETURN
! 225: END IF
! 226: *
! 227: * Quick return if possible
! 228: *
! 229: IF( N.EQ.0 .OR. M.EQ.0 ) RETURN
! 230: *
! 231: DO I = 1, N
! 232: *
! 233: * Generate elementary reflector H(I) to annihilate B(:,I)
! 234: *
! 235: P = M-L+MIN( L, I )
! 236: CALL DLARFG( P+1, A( I, I ), B( 1, I ), 1, T( I, 1 ) )
! 237: IF( I.LT.N ) THEN
! 238: *
! 239: * W(1:N-I) := C(I:M,I+1:N)^H * C(I:M,I) [use W = T(:,N)]
! 240: *
! 241: DO J = 1, N-I
! 242: T( J, N ) = (A( I, I+J ))
! 243: END DO
! 244: CALL DGEMV( 'T', P, N-I, ONE, B( 1, I+1 ), LDB,
! 245: $ B( 1, I ), 1, ONE, T( 1, N ), 1 )
! 246: *
! 247: * C(I:M,I+1:N) = C(I:m,I+1:N) + alpha*C(I:M,I)*W(1:N-1)^H
! 248: *
! 249: ALPHA = -(T( I, 1 ))
! 250: DO J = 1, N-I
! 251: A( I, I+J ) = A( I, I+J ) + ALPHA*(T( J, N ))
! 252: END DO
! 253: CALL DGER( P, N-I, ALPHA, B( 1, I ), 1,
! 254: $ T( 1, N ), 1, B( 1, I+1 ), LDB )
! 255: END IF
! 256: END DO
! 257: *
! 258: DO I = 2, N
! 259: *
! 260: * T(1:I-1,I) := C(I:M,1:I-1)^H * (alpha * C(I:M,I))
! 261: *
! 262: ALPHA = -T( I, 1 )
! 263:
! 264: DO J = 1, I-1
! 265: T( J, I ) = ZERO
! 266: END DO
! 267: P = MIN( I-1, L )
! 268: MP = MIN( M-L+1, M )
! 269: NP = MIN( P+1, N )
! 270: *
! 271: * Triangular part of B2
! 272: *
! 273: DO J = 1, P
! 274: T( J, I ) = ALPHA*B( M-L+J, I )
! 275: END DO
! 276: CALL DTRMV( 'U', 'T', 'N', P, B( MP, 1 ), LDB,
! 277: $ T( 1, I ), 1 )
! 278: *
! 279: * Rectangular part of B2
! 280: *
! 281: CALL DGEMV( 'T', L, I-1-P, ALPHA, B( MP, NP ), LDB,
! 282: $ B( MP, I ), 1, ZERO, T( NP, I ), 1 )
! 283: *
! 284: * B1
! 285: *
! 286: CALL DGEMV( 'T', M-L, I-1, ALPHA, B, LDB, B( 1, I ), 1,
! 287: $ ONE, T( 1, I ), 1 )
! 288: *
! 289: * T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
! 290: *
! 291: CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
! 292: *
! 293: * T(I,I) = tau(I)
! 294: *
! 295: T( I, I ) = T( I, 1 )
! 296: T( I, 1 ) = ZERO
! 297: END DO
! 298:
! 299: *
! 300: * End of DTPQRT2
! 301: *
! 302: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtpqrt2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack