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

1.3       bertrand    1: *> \brief \b DGEQRT2 computes a QR factorization of a general real or complex matrix using the compact WY representation of 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 DGEQRT2 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqrt2.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqrt2.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqrt2.f">
1.1       bertrand   15: *> [TXT]</a>
1.7       bertrand   16: *> \endhtmlonly
1.1       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGEQRT2( M, N, A, LDA, T, LDT, INFO )
1.7       bertrand   22: *
1.1       bertrand   23: *       .. Scalar Arguments ..
                     24: *       INTEGER   INFO, LDA, LDT, M, N
                     25: *       ..
                     26: *       .. Array Arguments ..
                     27: *       DOUBLE PRECISION   A( LDA, * ), T( LDT, * )
                     28: *       ..
1.7       bertrand   29: *
1.1       bertrand   30: *
                     31: *> \par Purpose:
                     32: *  =============
                     33: *>
                     34: *> \verbatim
                     35: *>
1.7       bertrand   36: *> DGEQRT2 computes a QR factorization of a real M-by-N matrix A,
                     37: *> using the compact WY representation of Q.
1.1       bertrand   38: *> \endverbatim
                     39: *
                     40: *  Arguments:
                     41: *  ==========
                     42: *
                     43: *> \param[in] M
                     44: *> \verbatim
                     45: *>          M is INTEGER
                     46: *>          The number of rows of the matrix A.  M >= N.
                     47: *> \endverbatim
                     48: *>
                     49: *> \param[in] N
                     50: *> \verbatim
                     51: *>          N is INTEGER
                     52: *>          The number of columns of the matrix A.  N >= 0.
                     53: *> \endverbatim
                     54: *>
                     55: *> \param[in,out] A
                     56: *> \verbatim
                     57: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
                     58: *>          On entry, the real M-by-N matrix A.  On exit, the elements on and
                     59: *>          above the diagonal contain the N-by-N upper triangular matrix R; the
                     60: *>          elements below the diagonal are the columns of V.  See below for
                     61: *>          further details.
                     62: *> \endverbatim
                     63: *>
                     64: *> \param[in] LDA
                     65: *> \verbatim
                     66: *>          LDA is INTEGER
                     67: *>          The leading dimension of the array A.  LDA >= max(1,M).
                     68: *> \endverbatim
                     69: *>
                     70: *> \param[out] T
                     71: *> \verbatim
                     72: *>          T is DOUBLE PRECISION array, dimension (LDT,N)
                     73: *>          The N-by-N upper triangular factor of the block reflector.
                     74: *>          The elements on and above the diagonal contain the block
                     75: *>          reflector T; the elements below the diagonal are not used.
                     76: *>          See below for further details.
                     77: *> \endverbatim
                     78: *>
                     79: *> \param[in] LDT
                     80: *> \verbatim
                     81: *>          LDT is INTEGER
                     82: *>          The leading dimension of the array T.  LDT >= max(1,N).
                     83: *> \endverbatim
                     84: *>
                     85: *> \param[out] INFO
                     86: *> \verbatim
                     87: *>          INFO is INTEGER
                     88: *>          = 0: successful exit
                     89: *>          < 0: if INFO = -i, the i-th argument had an illegal value
                     90: *> \endverbatim
                     91: *
                     92: *  Authors:
                     93: *  ========
                     94: *
1.7       bertrand   95: *> \author Univ. of Tennessee
                     96: *> \author Univ. of California Berkeley
                     97: *> \author Univ. of Colorado Denver
                     98: *> \author NAG Ltd.
1.1       bertrand   99: *
                    100: *> \ingroup doubleGEcomputational
                    101: *
                    102: *> \par Further Details:
                    103: *  =====================
                    104: *>
                    105: *> \verbatim
                    106: *>
                    107: *>  The matrix V stores the elementary reflectors H(i) in the i-th column
                    108: *>  below the diagonal. For example, if M=5 and N=3, the matrix V is
                    109: *>
                    110: *>               V = (  1       )
                    111: *>                   ( v1  1    )
                    112: *>                   ( v1 v2  1 )
                    113: *>                   ( v1 v2 v3 )
                    114: *>                   ( v1 v2 v3 )
                    115: *>
                    116: *>  where the vi's represent the vectors which define H(i), which are returned
                    117: *>  in the matrix A.  The 1's along the diagonal of V are not stored in A.  The
                    118: *>  block reflector H is then given by
                    119: *>
                    120: *>               H = I - V * T * V**T
                    121: *>
                    122: *>  where V**T is the transpose of V.
                    123: *> \endverbatim
                    124: *>
                    125: *  =====================================================================
                    126:       SUBROUTINE DGEQRT2( M, N, A, LDA, T, LDT, INFO )
                    127: *
1.10    ! bertrand  128: *  -- LAPACK computational routine --
1.1       bertrand  129: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    130: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    131: *
                    132: *     .. Scalar Arguments ..
                    133:       INTEGER   INFO, LDA, LDT, M, N
                    134: *     ..
                    135: *     .. Array Arguments ..
                    136:       DOUBLE PRECISION   A( LDA, * ), T( LDT, * )
                    137: *     ..
                    138: *
                    139: *  =====================================================================
                    140: *
                    141: *     .. Parameters ..
                    142:       DOUBLE PRECISION  ONE, ZERO
                    143:       PARAMETER( ONE = 1.0D+00, ZERO = 0.0D+00 )
                    144: *     ..
                    145: *     .. Local Scalars ..
                    146:       INTEGER   I, K
                    147:       DOUBLE PRECISION   AII, ALPHA
                    148: *     ..
                    149: *     .. External Subroutines ..
                    150:       EXTERNAL  DLARFG, DGEMV, DGER, DTRMV, XERBLA
                    151: *     ..
                    152: *     .. Executable Statements ..
                    153: *
                    154: *     Test the input arguments
                    155: *
                    156:       INFO = 0
1.10    ! bertrand  157:       IF( N.LT.0 ) THEN
        !           158:          INFO = -2
        !           159:       ELSE IF( M.LT.N ) THEN
1.1       bertrand  160:          INFO = -1
                    161:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    162:          INFO = -4
                    163:       ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
                    164:          INFO = -6
                    165:       END IF
                    166:       IF( INFO.NE.0 ) THEN
                    167:          CALL XERBLA( 'DGEQRT2', -INFO )
                    168:          RETURN
                    169:       END IF
1.7       bertrand  170: *
1.1       bertrand  171:       K = MIN( M, N )
                    172: *
                    173:       DO I = 1, K
                    174: *
                    175: *        Generate elem. refl. H(i) to annihilate A(i+1:m,i), tau(I) -> T(I,1)
                    176: *
                    177:          CALL DLARFG( M-I+1, A( I, I ), A( MIN( I+1, M ), I ), 1,
                    178:      $                T( I, 1 ) )
                    179:          IF( I.LT.N ) THEN
                    180: *
                    181: *           Apply H(i) to A(I:M,I+1:N) from the left
                    182: *
                    183:             AII = A( I, I )
                    184:             A( I, I ) = ONE
                    185: *
                    186: *           W(1:N-I) := A(I:M,I+1:N)^H * A(I:M,I) [W = T(:,N)]
                    187: *
1.7       bertrand  188:             CALL DGEMV( 'T',M-I+1, N-I, ONE, A( I, I+1 ), LDA,
1.1       bertrand  189:      $                  A( I, I ), 1, ZERO, T( 1, N ), 1 )
                    190: *
                    191: *           A(I:M,I+1:N) = A(I:m,I+1:N) + alpha*A(I:M,I)*W(1:N-1)^H
                    192: *
                    193:             ALPHA = -(T( I, 1 ))
1.7       bertrand  194:             CALL DGER( M-I+1, N-I, ALPHA, A( I, I ), 1,
1.1       bertrand  195:      $           T( 1, N ), 1, A( I, I+1 ), LDA )
                    196:             A( I, I ) = AII
                    197:          END IF
                    198:       END DO
                    199: *
                    200:       DO I = 2, N
                    201:          AII = A( I, I )
                    202:          A( I, I ) = ONE
                    203: *
                    204: *        T(1:I-1,I) := alpha * A(I:M,1:I-1)**T * A(I:M,I)
                    205: *
                    206:          ALPHA = -T( I, 1 )
1.7       bertrand  207:          CALL DGEMV( 'T', M-I+1, I-1, ALPHA, A( I, 1 ), LDA,
1.1       bertrand  208:      $               A( I, I ), 1, ZERO, T( 1, I ), 1 )
                    209:          A( I, I ) = AII
                    210: *
                    211: *        T(1:I-1,I) := T(1:I-1,1:I-1) * T(1:I-1,I)
                    212: *
                    213:          CALL DTRMV( 'U', 'N', 'N', I-1, T, LDT, T( 1, I ), 1 )
                    214: *
                    215: *           T(I,I) = tau(I)
                    216: *
                    217:             T( I, I ) = T( I, 1 )
                    218:             T( I, 1) = ZERO
                    219:       END DO
1.7       bertrand  220: 
1.1       bertrand  221: *
                    222: *     End of DGEQRT2
                    223: *
                    224:       END

CVSweb interface <joel.bertrand@systella.fr>