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

1.9       bertrand    1: *> \brief \b DLARZT
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download DLARZT + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarzt.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarzt.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarzt.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
                     22: * 
                     23: *       .. Scalar Arguments ..
                     24: *       CHARACTER          DIRECT, STOREV
                     25: *       INTEGER            K, LDT, LDV, N
                     26: *       ..
                     27: *       .. Array Arguments ..
                     28: *       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
                     29: *       ..
                     30: *  
                     31: *
                     32: *> \par Purpose:
                     33: *  =============
                     34: *>
                     35: *> \verbatim
                     36: *>
                     37: *> DLARZT forms the triangular factor T of a real block reflector
                     38: *> H of order > n, which is defined as a product of k elementary
                     39: *> reflectors.
                     40: *>
                     41: *> If DIRECT = 'F', H = H(1) H(2) . . . H(k) and T is upper triangular;
                     42: *>
                     43: *> If DIRECT = 'B', H = H(k) . . . H(2) H(1) and T is lower triangular.
                     44: *>
                     45: *> If STOREV = 'C', the vector which defines the elementary reflector
                     46: *> H(i) is stored in the i-th column of the array V, and
                     47: *>
                     48: *>    H  =  I - V * T * V**T
                     49: *>
                     50: *> If STOREV = 'R', the vector which defines the elementary reflector
                     51: *> H(i) is stored in the i-th row of the array V, and
                     52: *>
                     53: *>    H  =  I - V**T * T * V
                     54: *>
                     55: *> Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
                     56: *> \endverbatim
                     57: *
                     58: *  Arguments:
                     59: *  ==========
                     60: *
                     61: *> \param[in] DIRECT
                     62: *> \verbatim
                     63: *>          DIRECT is CHARACTER*1
                     64: *>          Specifies the order in which the elementary reflectors are
                     65: *>          multiplied to form the block reflector:
                     66: *>          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
                     67: *>          = 'B': H = H(k) . . . H(2) H(1) (Backward)
                     68: *> \endverbatim
                     69: *>
                     70: *> \param[in] STOREV
                     71: *> \verbatim
                     72: *>          STOREV is CHARACTER*1
                     73: *>          Specifies how the vectors which define the elementary
                     74: *>          reflectors are stored (see also Further Details):
                     75: *>          = 'C': columnwise                        (not supported yet)
                     76: *>          = 'R': rowwise
                     77: *> \endverbatim
                     78: *>
                     79: *> \param[in] N
                     80: *> \verbatim
                     81: *>          N is INTEGER
                     82: *>          The order of the block reflector H. N >= 0.
                     83: *> \endverbatim
                     84: *>
                     85: *> \param[in] K
                     86: *> \verbatim
                     87: *>          K is INTEGER
                     88: *>          The order of the triangular factor T (= the number of
                     89: *>          elementary reflectors). K >= 1.
                     90: *> \endverbatim
                     91: *>
                     92: *> \param[in,out] V
                     93: *> \verbatim
                     94: *>          V is DOUBLE PRECISION array, dimension
                     95: *>                               (LDV,K) if STOREV = 'C'
                     96: *>                               (LDV,N) if STOREV = 'R'
                     97: *>          The matrix V. See further details.
                     98: *> \endverbatim
                     99: *>
                    100: *> \param[in] LDV
                    101: *> \verbatim
                    102: *>          LDV is INTEGER
                    103: *>          The leading dimension of the array V.
                    104: *>          If STOREV = 'C', LDV >= max(1,N); if STOREV = 'R', LDV >= K.
                    105: *> \endverbatim
                    106: *>
                    107: *> \param[in] TAU
                    108: *> \verbatim
                    109: *>          TAU is DOUBLE PRECISION array, dimension (K)
                    110: *>          TAU(i) must contain the scalar factor of the elementary
                    111: *>          reflector H(i).
                    112: *> \endverbatim
                    113: *>
                    114: *> \param[out] T
                    115: *> \verbatim
                    116: *>          T is DOUBLE PRECISION array, dimension (LDT,K)
                    117: *>          The k by k triangular factor T of the block reflector.
                    118: *>          If DIRECT = 'F', T is upper triangular; if DIRECT = 'B', T is
                    119: *>          lower triangular. The rest of the array is not used.
                    120: *> \endverbatim
                    121: *>
                    122: *> \param[in] LDT
                    123: *> \verbatim
                    124: *>          LDT is INTEGER
                    125: *>          The leading dimension of the array T. LDT >= K.
                    126: *> \endverbatim
                    127: *
                    128: *  Authors:
                    129: *  ========
                    130: *
                    131: *> \author Univ. of Tennessee 
                    132: *> \author Univ. of California Berkeley 
                    133: *> \author Univ. of Colorado Denver 
                    134: *> \author NAG Ltd. 
                    135: *
                    136: *> \date November 2011
                    137: *
                    138: *> \ingroup doubleOTHERcomputational
                    139: *
                    140: *> \par Contributors:
                    141: *  ==================
                    142: *>
                    143: *>    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
                    144: *
                    145: *> \par Further Details:
                    146: *  =====================
                    147: *>
                    148: *> \verbatim
                    149: *>
                    150: *>  The shape of the matrix V and the storage of the vectors which define
                    151: *>  the H(i) is best illustrated by the following example with n = 5 and
                    152: *>  k = 3. The elements equal to 1 are not stored; the corresponding
                    153: *>  array elements are modified but restored on exit. The rest of the
                    154: *>  array is not used.
                    155: *>
                    156: *>  DIRECT = 'F' and STOREV = 'C':         DIRECT = 'F' and STOREV = 'R':
                    157: *>
                    158: *>                                              ______V_____
                    159: *>         ( v1 v2 v3 )                        /            \
                    160: *>         ( v1 v2 v3 )                      ( v1 v1 v1 v1 v1 . . . . 1 )
                    161: *>     V = ( v1 v2 v3 )                      ( v2 v2 v2 v2 v2 . . . 1   )
                    162: *>         ( v1 v2 v3 )                      ( v3 v3 v3 v3 v3 . . 1     )
                    163: *>         ( v1 v2 v3 )
                    164: *>            .  .  .
                    165: *>            .  .  .
                    166: *>            1  .  .
                    167: *>               1  .
                    168: *>                  1
                    169: *>
                    170: *>  DIRECT = 'B' and STOREV = 'C':         DIRECT = 'B' and STOREV = 'R':
                    171: *>
                    172: *>                                                        ______V_____
                    173: *>            1                                          /            \
                    174: *>            .  1                           ( 1 . . . . v1 v1 v1 v1 v1 )
                    175: *>            .  .  1                        ( . 1 . . . v2 v2 v2 v2 v2 )
                    176: *>            .  .  .                        ( . . 1 . . v3 v3 v3 v3 v3 )
                    177: *>            .  .  .
                    178: *>         ( v1 v2 v3 )
                    179: *>         ( v1 v2 v3 )
                    180: *>     V = ( v1 v2 v3 )
                    181: *>         ( v1 v2 v3 )
                    182: *>         ( v1 v2 v3 )
                    183: *> \endverbatim
                    184: *>
                    185: *  =====================================================================
1.1       bertrand  186:       SUBROUTINE DLARZT( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT )
                    187: *
1.9       bertrand  188: *  -- LAPACK computational routine (version 3.4.0) --
1.1       bertrand  189: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    190: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9       bertrand  191: *     November 2011
1.1       bertrand  192: *
                    193: *     .. Scalar Arguments ..
                    194:       CHARACTER          DIRECT, STOREV
                    195:       INTEGER            K, LDT, LDV, N
                    196: *     ..
                    197: *     .. Array Arguments ..
                    198:       DOUBLE PRECISION   T( LDT, * ), TAU( * ), V( LDV, * )
                    199: *     ..
                    200: *
                    201: *  =====================================================================
                    202: *
                    203: *     .. Parameters ..
                    204:       DOUBLE PRECISION   ZERO
                    205:       PARAMETER          ( ZERO = 0.0D+0 )
                    206: *     ..
                    207: *     .. Local Scalars ..
                    208:       INTEGER            I, INFO, J
                    209: *     ..
                    210: *     .. External Subroutines ..
                    211:       EXTERNAL           DGEMV, DTRMV, XERBLA
                    212: *     ..
                    213: *     .. External Functions ..
                    214:       LOGICAL            LSAME
                    215:       EXTERNAL           LSAME
                    216: *     ..
                    217: *     .. Executable Statements ..
                    218: *
                    219: *     Check for currently supported options
                    220: *
                    221:       INFO = 0
                    222:       IF( .NOT.LSAME( DIRECT, 'B' ) ) THEN
                    223:          INFO = -1
                    224:       ELSE IF( .NOT.LSAME( STOREV, 'R' ) ) THEN
                    225:          INFO = -2
                    226:       END IF
                    227:       IF( INFO.NE.0 ) THEN
                    228:          CALL XERBLA( 'DLARZT', -INFO )
                    229:          RETURN
                    230:       END IF
                    231: *
                    232:       DO 20 I = K, 1, -1
                    233:          IF( TAU( I ).EQ.ZERO ) THEN
                    234: *
                    235: *           H(i)  =  I
                    236: *
                    237:             DO 10 J = I, K
                    238:                T( J, I ) = ZERO
                    239:    10       CONTINUE
                    240:          ELSE
                    241: *
                    242: *           general case
                    243: *
                    244:             IF( I.LT.K ) THEN
                    245: *
1.8       bertrand  246: *              T(i+1:k,i) = - tau(i) * V(i+1:k,1:n) * V(i,1:n)**T
1.1       bertrand  247: *
                    248:                CALL DGEMV( 'No transpose', K-I, N, -TAU( I ),
                    249:      $                     V( I+1, 1 ), LDV, V( I, 1 ), LDV, ZERO,
                    250:      $                     T( I+1, I ), 1 )
                    251: *
                    252: *              T(i+1:k,i) = T(i+1:k,i+1:k) * T(i+1:k,i)
                    253: *
                    254:                CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', K-I,
                    255:      $                     T( I+1, I+1 ), LDT, T( I+1, I ), 1 )
                    256:             END IF
                    257:             T( I, I ) = TAU( I )
                    258:          END IF
                    259:    20 CONTINUE
                    260:       RETURN
                    261: *
                    262: *     End of DLARZT
                    263: *
                    264:       END

CVSweb interface <joel.bertrand@systella.fr>