001:       SUBROUTINE SLAHRD( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY )
002: *
003: *  -- LAPACK auxiliary routine (version 3.2) --
004: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
005: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            K, LDA, LDT, LDY, N, NB
010: *     ..
011: *     .. Array Arguments ..
012:       REAL               A( LDA, * ), T( LDT, NB ), TAU( NB ),
013:      $                   Y( LDY, NB )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  SLAHRD reduces the first NB columns of a real general n-by-(n-k+1)
020: *  matrix A so that elements below the k-th subdiagonal are zero. The
021: *  reduction is performed by an orthogonal similarity transformation
022: *  Q' * A * Q. The routine returns the matrices V and T which determine
023: *  Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T.
024: *
025: *  This is an OBSOLETE auxiliary routine. 
026: *  This routine will be 'deprecated' in a  future release.
027: *  Please use the new routine SLAHR2 instead.
028: *
029: *  Arguments
030: *  =========
031: *
032: *  N       (input) INTEGER
033: *          The order of the matrix A.
034: *
035: *  K       (input) INTEGER
036: *          The offset for the reduction. Elements below the k-th
037: *          subdiagonal in the first NB columns are reduced to zero.
038: *
039: *  NB      (input) INTEGER
040: *          The number of columns to be reduced.
041: *
042: *  A       (input/output) REAL array, dimension (LDA,N-K+1)
043: *          On entry, the n-by-(n-k+1) general matrix A.
044: *          On exit, the elements on and above the k-th subdiagonal in
045: *          the first NB columns are overwritten with the corresponding
046: *          elements of the reduced matrix; the elements below the k-th
047: *          subdiagonal, with the array TAU, represent the matrix Q as a
048: *          product of elementary reflectors. The other columns of A are
049: *          unchanged. See Further Details.
050: *
051: *  LDA     (input) INTEGER
052: *          The leading dimension of the array A.  LDA >= max(1,N).
053: *
054: *  TAU     (output) REAL array, dimension (NB)
055: *          The scalar factors of the elementary reflectors. See Further
056: *          Details.
057: *
058: *  T       (output) REAL array, dimension (LDT,NB)
059: *          The upper triangular matrix T.
060: *
061: *  LDT     (input) INTEGER
062: *          The leading dimension of the array T.  LDT >= NB.
063: *
064: *  Y       (output) REAL array, dimension (LDY,NB)
065: *          The n-by-nb matrix Y.
066: *
067: *  LDY     (input) INTEGER
068: *          The leading dimension of the array Y. LDY >= N.
069: *
070: *  Further Details
071: *  ===============
072: *
073: *  The matrix Q is represented as a product of nb elementary reflectors
074: *
075: *     Q = H(1) H(2) . . . H(nb).
076: *
077: *  Each H(i) has the form
078: *
079: *     H(i) = I - tau * v * v'
080: *
081: *  where tau is a real scalar, and v is a real vector with
082: *  v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in
083: *  A(i+k+1:n,i), and tau in TAU(i).
084: *
085: *  The elements of the vectors v together form the (n-k+1)-by-nb matrix
086: *  V which is needed, with T and Y, to apply the transformation to the
087: *  unreduced part of the matrix, using an update of the form:
088: *  A := (I - V*T*V') * (A - Y*V').
089: *
090: *  The contents of A on exit are illustrated by the following example
091: *  with n = 7, k = 3 and nb = 2:
092: *
093: *     ( a   h   a   a   a )
094: *     ( a   h   a   a   a )
095: *     ( a   h   a   a   a )
096: *     ( h   h   a   a   a )
097: *     ( v1  h   a   a   a )
098: *     ( v1  v2  a   a   a )
099: *     ( v1  v2  a   a   a )
100: *
101: *  where a denotes an element of the original matrix A, h denotes a
102: *  modified element of the upper Hessenberg matrix H, and vi denotes an
103: *  element of the vector defining H(i).
104: *
105: *  =====================================================================
106: *
107: *     .. Parameters ..
108:       REAL               ZERO, ONE
109:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
110: *     ..
111: *     .. Local Scalars ..
112:       INTEGER            I
113:       REAL               EI
114: *     ..
115: *     .. External Subroutines ..
116:       EXTERNAL           SAXPY, SCOPY, SGEMV, SLARFG, SSCAL, STRMV
117: *     ..
118: *     .. Intrinsic Functions ..
119:       INTRINSIC          MIN
120: *     ..
121: *     .. Executable Statements ..
122: *
123: *     Quick return if possible
124: *
125:       IF( N.LE.1 )
126:      $   RETURN
127: *
128:       DO 10 I = 1, NB
129:          IF( I.GT.1 ) THEN
130: *
131: *           Update A(1:n,i)
132: *
133: *           Compute i-th column of A - Y * V'
134: *
135:             CALL SGEMV( 'No transpose', N, I-1, -ONE, Y, LDY,
136:      $                  A( K+I-1, 1 ), LDA, ONE, A( 1, I ), 1 )
137: *
138: *           Apply I - V * T' * V' to this column (call it b) from the
139: *           left, using the last column of T as workspace
140: *
141: *           Let  V = ( V1 )   and   b = ( b1 )   (first I-1 rows)
142: *                    ( V2 )             ( b2 )
143: *
144: *           where V1 is unit lower triangular
145: *
146: *           w := V1' * b1
147: *
148:             CALL SCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 )
149:             CALL STRMV( 'Lower', 'Transpose', 'Unit', I-1, A( K+1, 1 ),
150:      $                  LDA, T( 1, NB ), 1 )
151: *
152: *           w := w + V2'*b2
153: *
154:             CALL SGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ),
155:      $                  LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 )
156: *
157: *           w := T'*w
158: *
159:             CALL STRMV( 'Upper', 'Transpose', 'Non-unit', I-1, T, LDT,
160:      $                  T( 1, NB ), 1 )
161: *
162: *           b2 := b2 - V2*w
163: *
164:             CALL SGEMV( 'No transpose', N-K-I+1, I-1, -ONE, A( K+I, 1 ),
165:      $                  LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 )
166: *
167: *           b1 := b1 - V1*w
168: *
169:             CALL STRMV( 'Lower', 'No transpose', 'Unit', I-1,
170:      $                  A( K+1, 1 ), LDA, T( 1, NB ), 1 )
171:             CALL SAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 )
172: *
173:             A( K+I-1, I-1 ) = EI
174:          END IF
175: *
176: *        Generate the elementary reflector H(i) to annihilate
177: *        A(k+i+1:n,i)
178: *
179:          CALL SLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1,
180:      $                TAU( I ) )
181:          EI = A( K+I, I )
182:          A( K+I, I ) = ONE
183: *
184: *        Compute  Y(1:n,i)
185: *
186:          CALL SGEMV( 'No transpose', N, N-K-I+1, ONE, A( 1, I+1 ), LDA,
187:      $               A( K+I, I ), 1, ZERO, Y( 1, I ), 1 )
188:          CALL SGEMV( 'Transpose', N-K-I+1, I-1, ONE, A( K+I, 1 ), LDA,
189:      $               A( K+I, I ), 1, ZERO, T( 1, I ), 1 )
190:          CALL SGEMV( 'No transpose', N, I-1, -ONE, Y, LDY, T( 1, I ), 1,
191:      $               ONE, Y( 1, I ), 1 )
192:          CALL SSCAL( N, TAU( I ), Y( 1, I ), 1 )
193: *
194: *        Compute T(1:i,i)
195: *
196:          CALL SSCAL( I-1, -TAU( I ), T( 1, I ), 1 )
197:          CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, T, LDT,
198:      $               T( 1, I ), 1 )
199:          T( I, I ) = TAU( I )
200: *
201:    10 CONTINUE
202:       A( K+NB, NB ) = EI
203: *
204:       RETURN
205: *
206: *     End of SLAHRD
207: *
208:       END
209: