Annotation of rpl/lapack/lapack/dtpqrt2.f, revision 1.10

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

CVSweb interface <joel.bertrand@systella.fr>