LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
clahr2.f
Go to the documentation of this file.
00001 *> \brief \b CLAHR2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CLAHR2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clahr2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clahr2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clahr2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            K, LDA, LDT, LDY, N, NB
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       COMPLEX            A( LDA, * ), T( LDT, NB ), TAU( NB ),
00028 *      $                   Y( LDY, NB )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> CLAHR2 reduces the first NB columns of A complex general n-BY-(n-k+1)
00038 *> matrix A so that elements below the k-th subdiagonal are zero. The
00039 *> reduction is performed by an unitary similarity transformation
00040 *> Q**H * A * Q. The routine returns the matrices V and T which determine
00041 *> Q as a block reflector I - V*T*v**H, and also the matrix Y = A * V * T.
00042 *>
00043 *> This is an auxiliary routine called by CGEHRD.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] N
00050 *> \verbatim
00051 *>          N is INTEGER
00052 *>          The order of the matrix A.
00053 *> \endverbatim
00054 *>
00055 *> \param[in] K
00056 *> \verbatim
00057 *>          K is INTEGER
00058 *>          The offset for the reduction. Elements below the k-th
00059 *>          subdiagonal in the first NB columns are reduced to zero.
00060 *>          K < N.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] NB
00064 *> \verbatim
00065 *>          NB is INTEGER
00066 *>          The number of columns to be reduced.
00067 *> \endverbatim
00068 *>
00069 *> \param[in,out] A
00070 *> \verbatim
00071 *>          A is COMPLEX array, dimension (LDA,N-K+1)
00072 *>          On entry, the n-by-(n-k+1) general matrix A.
00073 *>          On exit, the elements on and above the k-th subdiagonal in
00074 *>          the first NB columns are overwritten with the corresponding
00075 *>          elements of the reduced matrix; the elements below the k-th
00076 *>          subdiagonal, with the array TAU, represent the matrix Q as a
00077 *>          product of elementary reflectors. The other columns of A are
00078 *>          unchanged. See Further Details.
00079 *> \endverbatim
00080 *>
00081 *> \param[in] LDA
00082 *> \verbatim
00083 *>          LDA is INTEGER
00084 *>          The leading dimension of the array A.  LDA >= max(1,N).
00085 *> \endverbatim
00086 *>
00087 *> \param[out] TAU
00088 *> \verbatim
00089 *>          TAU is COMPLEX array, dimension (NB)
00090 *>          The scalar factors of the elementary reflectors. See Further
00091 *>          Details.
00092 *> \endverbatim
00093 *>
00094 *> \param[out] T
00095 *> \verbatim
00096 *>          T is COMPLEX array, dimension (LDT,NB)
00097 *>          The upper triangular matrix T.
00098 *> \endverbatim
00099 *>
00100 *> \param[in] LDT
00101 *> \verbatim
00102 *>          LDT is INTEGER
00103 *>          The leading dimension of the array T.  LDT >= NB.
00104 *> \endverbatim
00105 *>
00106 *> \param[out] Y
00107 *> \verbatim
00108 *>          Y is COMPLEX array, dimension (LDY,NB)
00109 *>          The n-by-nb matrix Y.
00110 *> \endverbatim
00111 *>
00112 *> \param[in] LDY
00113 *> \verbatim
00114 *>          LDY is INTEGER
00115 *>          The leading dimension of the array Y. LDY >= N.
00116 *> \endverbatim
00117 *
00118 *  Authors:
00119 *  ========
00120 *
00121 *> \author Univ. of Tennessee 
00122 *> \author Univ. of California Berkeley 
00123 *> \author Univ. of Colorado Denver 
00124 *> \author NAG Ltd. 
00125 *
00126 *> \date November 2011
00127 *
00128 *> \ingroup complexOTHERauxiliary
00129 *
00130 *> \par Further Details:
00131 *  =====================
00132 *>
00133 *> \verbatim
00134 *>
00135 *>  The matrix Q is represented as a product of nb elementary reflectors
00136 *>
00137 *>     Q = H(1) H(2) . . . H(nb).
00138 *>
00139 *>  Each H(i) has the form
00140 *>
00141 *>     H(i) = I - tau * v * v**H
00142 *>
00143 *>  where tau is a complex scalar, and v is a complex vector with
00144 *>  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
00145 *>  A(i+k+1:n,i), and tau in TAU(i).
00146 *>
00147 *>  The elements of the vectors v together form the (n-k+1)-by-nb matrix
00148 *>  V which is needed, with T and Y, to apply the transformation to the
00149 *>  unreduced part of the matrix, using an update of the form:
00150 *>  A := (I - V*T*V**H) * (A - Y*V**H).
00151 *>
00152 *>  The contents of A on exit are illustrated by the following example
00153 *>  with n = 7, k = 3 and nb = 2:
00154 *>
00155 *>     ( a   a   a   a   a )
00156 *>     ( a   a   a   a   a )
00157 *>     ( a   a   a   a   a )
00158 *>     ( h   h   a   a   a )
00159 *>     ( v1  h   a   a   a )
00160 *>     ( v1  v2  a   a   a )
00161 *>     ( v1  v2  a   a   a )
00162 *>
00163 *>  where a denotes an element of the original matrix A, h denotes a
00164 *>  modified element of the upper Hessenberg matrix H, and vi denotes an
00165 *>  element of the vector defining H(i).
00166 *>
00167 *>  This subroutine is a slight modification of LAPACK-3.0's DLAHRD
00168 *>  incorporating improvements proposed by Quintana-Orti and Van de
00169 *>  Gejin. Note that the entries of A(1:K,2:NB) differ from those
00170 *>  returned by the original LAPACK-3.0's DLAHRD routine. (This
00171 *>  subroutine is not backward compatible with LAPACK-3.0's DLAHRD.)
00172 *> \endverbatim
00173 *
00174 *> \par References:
00175 *  ================
00176 *>
00177 *>  Gregorio Quintana-Orti and Robert van de Geijn, "Improving the
00178 *>  performance of reduction to Hessenberg form," ACM Transactions on
00179 *>  Mathematical Software, 32(2):180-194, June 2006.
00180 *>
00181 *  =====================================================================
00182       SUBROUTINE CLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
00183 *
00184 *  -- LAPACK auxiliary routine (version 3.4.0) --
00185 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00186 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00187 *     November 2011
00188 *
00189 *     .. Scalar Arguments ..
00190       INTEGER            K, LDA, LDT, LDY, N, NB
00191 *     ..
00192 *     .. Array Arguments ..
00193       COMPLEX            A( LDA, * ), T( LDT, NB ), TAU( NB ),
00194      $                   Y( LDY, NB )
00195 *     ..
00196 *
00197 *  =====================================================================
00198 *
00199 *     .. Parameters ..
00200       COMPLEX            ZERO, ONE
00201       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ), 
00202      $                     ONE = ( 1.0E+0, 0.0E+0 ) )
00203 *     ..
00204 *     .. Local Scalars ..
00205       INTEGER            I
00206       COMPLEX            EI
00207 *     ..
00208 *     .. External Subroutines ..
00209       EXTERNAL           CAXPY, CCOPY, CGEMM, CGEMV, CLACPY,
00210      $                   CLARFG, CSCAL, CTRMM, CTRMV, CLACGV
00211 *     ..
00212 *     .. Intrinsic Functions ..
00213       INTRINSIC          MIN
00214 *     ..
00215 *     .. Executable Statements ..
00216 *
00217 *     Quick return if possible
00218 *
00219       IF( N.LE.1 )
00220      $   RETURN
00221 *
00222       DO 10 I = 1, NB
00223          IF( I.GT.1 ) THEN
00224 *
00225 *           Update A(K+1:N,I)
00226 *
00227 *           Update I-th column of A - Y * V**H
00228 *
00229             CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) 
00230             CALL CGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY,
00231      $                  A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 )
00232             CALL CLACGV( I-1, A( K+I-1, 1 ), LDA ) 
00233 *
00234 *           Apply I - V * T**H * V**H to this column (call it b) from the
00235 *           left, using the last column of T as workspace
00236 *
00237 *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
00238 *                    ( V2 )             ( b2 )
00239 *
00240 *           where V1 is unit lower triangular
00241 *
00242 *           w := V1**H * b1
00243 *
00244             CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
00245             CALL CTRMV( 'Lower', 'Conjugate transpose', 'UNIT', 
00246      $                  I-1, A( K+1, 1 ),
00247      $                  LDA, T( 1, NB ), 1 )
00248 *
00249 *           w := w + V2**H * b2
00250 *
00251             CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, 
00252      $                  ONE, A( K+I, 1 ),
00253      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
00254 *
00255 *           w := T**H * w
00256 *
00257             CALL CTRMV( 'Upper', 'Conjugate transpose', 'NON-UNIT', 
00258      $                  I-1, T, LDT,
00259      $                  T( 1, NB ), 1 )
00260 *
00261 *           b2 := b2 - V2*w
00262 *
00263             CALL CGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, 
00264      $                  A( K+I, 1 ),
00265      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
00266 *
00267 *           b1 := b1 - V1*w
00268 *
00269             CALL CTRMV( 'Lower', 'NO TRANSPOSE', 
00270      $                  'UNIT', I-1,
00271      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
00272             CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
00273 *
00274             A( K+I-1, I-1 ) = EI
00275          END IF
00276 *
00277 *        Generate the elementary reflector H(I) to annihilate
00278 *        A(K+I+1:N,I)
00279 *
00280          CALL CLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
00281      $                TAU( I ) )
00282          EI = A( K+I, I )
00283          A( K+I, I ) = ONE
00284 *
00285 *        Compute  Y(K+1:N,I)
00286 *
00287          CALL CGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, 
00288      $               ONE, A( K+1, I+1 ),
00289      $               LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 )
00290          CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, 
00291      $               ONE, A( K+I, 1 ), LDA,
00292      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
00293          CALL CGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, 
00294      $               Y( K+1, 1 ), LDY,
00295      $               T( 1, I ), 1, ONE, Y( K+1, I ), 1 )
00296          CALL CSCAL( N-K, TAU( I ), Y( K+1, I ), 1 )
00297 *
00298 *        Compute T(1:I,I)
00299 *
00300          CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
00301          CALL CTRMV( 'Upper', 'No Transpose', 'NON-UNIT', 
00302      $               I-1, T, LDT,
00303      $               T( 1, I ), 1 )
00304          T( I, I ) = TAU( I )
00305 *
00306    10 CONTINUE
00307       A( K+NB, NB ) = EI
00308 *
00309 *     Compute Y(1:K,1:NB)
00310 *
00311       CALL CLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY )
00312       CALL CTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', 
00313      $            'UNIT', K, NB,
00314      $            ONE, A( K+1, 1 ), LDA, Y, LDY )
00315       IF( N.GT.K+NB )
00316      $   CALL CGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, 
00317      $               NB, N-K-NB, ONE,
00318      $               A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y,
00319      $               LDY )
00320       CALL CTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', 
00321      $            'NON-UNIT', K, NB,
00322      $            ONE, T, LDT, Y, LDY )
00323 *
00324       RETURN
00325 *
00326 *     End of CLAHR2
00327 *
00328       END
 All Files Functions