db/da4/slahrd_8f_source.html

00001       SUBROUTINE SLAHRD( 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       REAL               A( LDA, * ), T( LDT, NB ), TAU( NB ),
00013      $                   Y( LDY, NB )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SLAHRD reduces the first NB columns of a real 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 an orthogonal 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 SLAHR2 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) REAL 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) REAL array, dimension (NB)
00055 *          The scalar factors of the elementary reflectors. See Further
00056 *          Details.
00057 *
00058 *  T       (output) REAL 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) REAL 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 >= 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 real scalar, and v is a real 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       REAL               ZERO, ONE
00109       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00110 *     ..
00111 *     .. Local Scalars ..
00112       INTEGER            I
00113       REAL               EI
00114 *     ..
00115 *     .. External Subroutines ..
00116       EXTERNAL           SAXPY, SCOPY, SGEMV, SLARFG, SSCAL, STRMV
00117 *     ..
00118 *     .. Intrinsic Functions ..
00119       INTRINSIC          MIN
00120 *     ..
00121 *     .. Executable Statements ..
00122 *
00123 *     Quick return if possible
00124 *
00125       IF( N.LE.1 )
00126      $   RETURN
00127 *
00128       DO 10 I = 1, NB
00129          IF( I.GT.1 ) THEN
00130 *
00131 *           Update A(1:n,i)
00132 *
00133 *           Compute i-th column of A - Y * V'
00134 *
00135             CALL SGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
00136      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
00137 *
00138 *           Apply I - V * T' * V' to this column (call it b) from the
00139 *           left, using the last column of T as workspace
00140 *
00141 *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
00142 *                    ( V2 )             ( b2 )
00143 *
00144 *           where V1 is unit lower triangular
00145 *
00146 *           w := V1' * b1
00147 *
00148             CALL SCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
00149             CALL STRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
00150      $                  LDA, T( 1, NB ), 1 )
00151 *
00152 *           w := w + V2'*b2
00153 *
00154             CALL SGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
00155      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
00156 *
00157 *           w := T'*w
00158 *
00159             CALL STRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
00160      $                  T( 1, NB ), 1 )
00161 *
00162 *           b2 := b2 - V2*w
00163 *
00164             CALL SGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
00165      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
00166 *
00167 *           b1 := b1 - V1*w
00168 *
00169             CALL STRMV( 'Lower', 'No transpose', 'Unit', I-1,
00170      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
00171             CALL SAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
00172 *
00173             A( K+I-1, I-1 ) = EI
00174          END IF
00175 *
00176 *        Generate the elementary reflector H(i) to annihilate
00177 *        A(k+i+1:n,i)
00178 *
00179          CALL SLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
00180      $                TAU( I ) )
00181          EI = A( K+I, I )
00182          A( K+I, I ) = ONE
00183 *
00184 *        Compute  Y(1:n,i)
00185 *
00186          CALL SGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
00187      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
00188          CALL SGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
00189      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
00190          CALL SGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
00191      $               ONE, Y( 1, I ), 1 )
00192          CALL SSCAL( N, TAU( I ), Y( 1, I ), 1 )
00193 *
00194 *        Compute T(1:i,i)
00195 *
00196          CALL SSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
00197          CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
00198      $               T( 1, I ), 1 )
00199          T( I, I ) = TAU( I )
00200 *
00201    10 CONTINUE
00202       A( K+NB, NB ) = EI
00203 *
00204       RETURN
00205 *
00206 *     End of SLAHRD
00207 *
00208       END