LAPACK 3.3.0

clahrd.f

Go to the documentation of this file.
00001       SUBROUTINE CLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            K, LDA, LDT, LDY, N, NB
00010 *     ..
00011 *     .. Array Arguments ..
00012       COMPLEX            A( LDA, * ), T( LDT, NB ), TAU( NB ),
00013      $                   Y( LDY, NB )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  CLAHRD reduces the first NB columns of a complex general n-by-(n-k+1)
00020 *  matrix A so that elements below the k-th subdiagonal are zero. The
00021 *  reduction is performed by a unitary similarity transformation
00022 *  Q' * A * Q. The routine returns the matrices V and T which determine
00023 *  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
00024 *
00025 *  This is an OBSOLETE auxiliary routine. 
00026 *  This routine will be 'deprecated' in a  future release.
00027 *  Please use the new routine CLAHR2 instead.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  N       (input) INTEGER
00033 *          The order of the matrix A.
00034 *
00035 *  K       (input) INTEGER
00036 *          The offset for the reduction. Elements below the k-th
00037 *          subdiagonal in the first NB columns are reduced to zero.
00038 *
00039 *  NB      (input) INTEGER
00040 *          The number of columns to be reduced.
00041 *
00042 *  A       (input/output) COMPLEX array, dimension (LDA,N-K+1)
00043 *          On entry, the n-by-(n-k+1) general matrix A.
00044 *          On exit, the elements on and above the k-th subdiagonal in
00045 *          the first NB columns are overwritten with the corresponding
00046 *          elements of the reduced matrix; the elements below the k-th
00047 *          subdiagonal, with the array TAU, represent the matrix Q as a
00048 *          product of elementary reflectors. The other columns of A are
00049 *          unchanged. See Further Details.
00050 *
00051 *  LDA     (input) INTEGER
00052 *          The leading dimension of the array A.  LDA >= max(1,N).
00053 *
00054 *  TAU     (output) COMPLEX array, dimension (NB)
00055 *          The scalar factors of the elementary reflectors. See Further
00056 *          Details.
00057 *
00058 *  T       (output) COMPLEX array, dimension (LDT,NB)
00059 *          The upper triangular matrix T.
00060 *
00061 *  LDT     (input) INTEGER
00062 *          The leading dimension of the array T.  LDT >= NB.
00063 *
00064 *  Y       (output) COMPLEX array, dimension (LDY,NB)
00065 *          The n-by-nb matrix Y.
00066 *
00067 *  LDY     (input) INTEGER
00068 *          The leading dimension of the array Y. LDY >= max(1,N).
00069 *
00070 *  Further Details
00071 *  ===============
00072 *
00073 *  The matrix Q is represented as a product of nb elementary reflectors
00074 *
00075 *     Q = H(1) H(2) . . . H(nb).
00076 *
00077 *  Each H(i) has the form
00078 *
00079 *     H(i) = I - tau * v * v'
00080 *
00081 *  where tau is a complex scalar, and v is a complex vector with
00082 *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
00083 *  A(i+k+1:n,i), and tau in TAU(i).
00084 *
00085 *  The elements of the vectors v together form the (n-k+1)-by-nb matrix
00086 *  V which is needed, with T and Y, to apply the transformation to the
00087 *  unreduced part of the matrix, using an update of the form:
00088 *  A := (I - V*T*V') * (A - Y*V').
00089 *
00090 *  The contents of A on exit are illustrated by the following example
00091 *  with n = 7, k = 3 and nb = 2:
00092 *
00093 *     ( a   h   a   a   a )
00094 *     ( a   h   a   a   a )
00095 *     ( a   h   a   a   a )
00096 *     ( h   h   a   a   a )
00097 *     ( v1  h   a   a   a )
00098 *     ( v1  v2  a   a   a )
00099 *     ( v1  v2  a   a   a )
00100 *
00101 *  where a denotes an element of the original matrix A, h denotes a
00102 *  modified element of the upper Hessenberg matrix H, and vi denotes an
00103 *  element of the vector defining H(i).
00104 *
00105 *  =====================================================================
00106 *
00107 *     .. Parameters ..
00108       COMPLEX            ZERO, ONE
00109       PARAMETER          ( ZERO = ( 0.0E+0, 0.0E+0 ),
00110      $                   ONE = ( 1.0E+0, 0.0E+0 ) )
00111 *     ..
00112 *     .. Local Scalars ..
00113       INTEGER            I
00114       COMPLEX            EI
00115 *     ..
00116 *     .. External Subroutines ..
00117       EXTERNAL           CAXPY, CCOPY, CGEMV, CLACGV, CLARFG, CSCAL,
00118      $                   CTRMV
00119 *     ..
00120 *     .. Intrinsic Functions ..
00121       INTRINSIC          MIN
00122 *     ..
00123 *     .. Executable Statements ..
00124 *
00125 *     Quick return if possible
00126 *
00127       IF( N.LE.1 )
00128      $   RETURN
00129 *
00130       DO 10 I = 1, NB
00131          IF( I.GT.1 ) THEN
00132 *
00133 *           Update A(1:n,i)
00134 *
00135 *           Compute i-th column of A - Y * V'
00136 *
00137             CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
00138             CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
00139      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
00140             CALL CLACGV( I-1, A( K+I-1, 1 ), LDA )
00141 *
00142 *           Apply I - V * T' * V' to this column (call it b) from the
00143 *           left, using the last column of T as workspace
00144 *
00145 *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
00146 *                    ( V2 )             ( b2 )
00147 *
00148 *           where V1 is unit lower triangular
00149 *
00150 *           w := V1' * b1
00151 *
00152             CALL CCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
00153             CALL CTRMV( 'Lower', 'Conjugate transpose', 'Unit', I-1,
00154      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
00155 *
00156 *           w := w + V2'*b2
00157 *
00158             CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
00159      $                  A( K+I, 1 ), LDA, A( K+I, I ), 1, ONE,
00160      $                  T( 1, NB ), 1 )
00161 *
00162 *           w := T'*w
00163 *
00164             CALL CTRMV( 'Upper', 'Conjugate transpose', 'Non-unit', I-1,
00165      $                  T, LDT, T( 1, NB ), 1 )
00166 *
00167 *           b2 := b2 - V2*w
00168 *
00169             CALL CGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
00170      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
00171 *
00172 *           b1 := b1 - V1*w
00173 *
00174             CALL CTRMV( 'Lower', 'No transpose', 'Unit', I-1,
00175      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
00176             CALL CAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
00177 *
00178             A( K+I-1, I-1 ) = EI
00179          END IF
00180 *
00181 *        Generate the elementary reflector H(i) to annihilate
00182 *        A(k+i+1:n,i)
00183 *
00184          EI = A( K+I, I )
00185          CALL CLARFG( N-K-I+1, EI, A( MIN( K+I+1, N ), I ), 1,
00186      $                TAU( I ) )
00187          A( K+I, I ) = ONE
00188 *
00189 *        Compute  Y(1:n,i)
00190 *
00191          CALL CGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
00192      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
00193          CALL CGEMV( 'Conjugate transpose', N-K-I+1, I-1, ONE,
00194      $               A( K+I, 1 ), LDA, A( K+I, I ), 1, ZERO, T( 1, I ),
00195      $               1 )
00196          CALL CGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
00197      $               ONE, Y( 1, I ), 1 )
00198          CALL CSCAL( N, TAU( I ), Y( 1, I ), 1 )
00199 *
00200 *        Compute T(1:i,i)
00201 *
00202          CALL CSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
00203          CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
00204      $               T( 1, I ), 1 )
00205          T( I, I ) = TAU( I )
00206 *
00207    10 CONTINUE
00208       A( K+NB, NB ) = EI
00209 *
00210       RETURN
00211 *
00212 *     End of CLAHRD
00213 *
00214       END
 All Files Functions