LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
sgebrd.f
Go to the documentation of this file.
1*> \brief \b SGEBRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgebrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgebrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgebrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
22* INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LWORK, M, N
26* ..
27* .. Array Arguments ..
28* REAL A( LDA, * ), D( * ), E( * ), TAUP( * ),
29* \$ TAUQ( * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> SGEBRD reduces a general real M-by-N matrix A to upper or lower
39*> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
40*>
41*> If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
42*> \endverbatim
43*
44* Arguments:
45* ==========
46*
47*> \param[in] M
48*> \verbatim
49*> M is INTEGER
50*> The number of rows in the matrix A. M >= 0.
51*> \endverbatim
52*>
53*> \param[in] N
54*> \verbatim
55*> N is INTEGER
56*> The number of columns in the matrix A. N >= 0.
57*> \endverbatim
58*>
59*> \param[in,out] A
60*> \verbatim
61*> A is REAL array, dimension (LDA,N)
62*> On entry, the M-by-N general matrix to be reduced.
63*> On exit,
64*> if m >= n, the diagonal and the first superdiagonal are
65*> overwritten with the upper bidiagonal matrix B; the
66*> elements below the diagonal, with the array TAUQ, represent
67*> the orthogonal matrix Q as a product of elementary
68*> reflectors, and the elements above the first superdiagonal,
69*> with the array TAUP, represent the orthogonal matrix P as
70*> a product of elementary reflectors;
71*> if m < n, the diagonal and the first subdiagonal are
72*> overwritten with the lower bidiagonal matrix B; the
73*> elements below the first subdiagonal, with the array TAUQ,
74*> represent the orthogonal matrix Q as a product of
75*> elementary reflectors, and the elements above the diagonal,
76*> with the array TAUP, represent the orthogonal matrix P as
77*> a product of elementary reflectors.
78*> See Further Details.
79*> \endverbatim
80*>
81*> \param[in] LDA
82*> \verbatim
83*> LDA is INTEGER
84*> The leading dimension of the array A. LDA >= max(1,M).
85*> \endverbatim
86*>
87*> \param[out] D
88*> \verbatim
89*> D is REAL array, dimension (min(M,N))
90*> The diagonal elements of the bidiagonal matrix B:
91*> D(i) = A(i,i).
92*> \endverbatim
93*>
94*> \param[out] E
95*> \verbatim
96*> E is REAL array, dimension (min(M,N)-1)
97*> The off-diagonal elements of the bidiagonal matrix B:
98*> if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
99*> if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
100*> \endverbatim
101*>
102*> \param[out] TAUQ
103*> \verbatim
104*> TAUQ is REAL array, dimension (min(M,N))
105*> The scalar factors of the elementary reflectors which
106*> represent the orthogonal matrix Q. See Further Details.
107*> \endverbatim
108*>
109*> \param[out] TAUP
110*> \verbatim
111*> TAUP is REAL array, dimension (min(M,N))
112*> The scalar factors of the elementary reflectors which
113*> represent the orthogonal matrix P. See Further Details.
114*> \endverbatim
115*>
116*> \param[out] WORK
117*> \verbatim
118*> WORK is REAL array, dimension (MAX(1,LWORK))
119*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
120*> \endverbatim
121*>
122*> \param[in] LWORK
123*> \verbatim
124*> LWORK is INTEGER
125*> The length of the array WORK. LWORK >= max(1,M,N).
126*> For optimum performance LWORK >= (M+N)*NB, where NB
127*> is the optimal blocksize.
128*>
129*> If LWORK = -1, then a workspace query is assumed; the routine
130*> only calculates the optimal size of the WORK array, returns
131*> this value as the first entry of the WORK array, and no error
132*> message related to LWORK is issued by XERBLA.
133*> \endverbatim
134*>
135*> \param[out] INFO
136*> \verbatim
137*> INFO is INTEGER
138*> = 0: successful exit
139*> < 0: if INFO = -i, the i-th argument had an illegal value.
140*> \endverbatim
141*
142* Authors:
143* ========
144*
145*> \author Univ. of Tennessee
146*> \author Univ. of California Berkeley
147*> \author Univ. of Colorado Denver
148*> \author NAG Ltd.
149*
150*> \ingroup gebrd
151*
152*> \par Further Details:
153* =====================
154*>
155*> \verbatim
156*>
157*> The matrices Q and P are represented as products of elementary
158*> reflectors:
159*>
160*> If m >= n,
161*>
162*> Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1)
163*>
164*> Each H(i) and G(i) has the form:
165*>
166*> H(i) = I - tauq * v * v**T and G(i) = I - taup * u * u**T
167*>
168*> where tauq and taup are real scalars, and v and u are real vectors;
169*> v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
170*> u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
171*> tauq is stored in TAUQ(i) and taup in TAUP(i).
172*>
173*> If m < n,
174*>
175*> Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m)
176*>
177*> Each H(i) and G(i) has the form:
178*>
179*> H(i) = I - tauq * v * v**T and G(i) = I - taup * u * u**T
180*>
181*> where tauq and taup are real scalars, and v and u are real vectors;
182*> v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
183*> u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
184*> tauq is stored in TAUQ(i) and taup in TAUP(i).
185*>
186*> The contents of A on exit are illustrated by the following examples:
187*>
188*> m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
189*>
190*> ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
191*> ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
192*> ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
193*> ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
194*> ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
195*> ( v1 v2 v3 v4 v5 )
196*>
197*> where d and e denote diagonal and off-diagonal elements of B, vi
198*> denotes an element of the vector defining H(i), and ui an element of
199*> the vector defining G(i).
200*> \endverbatim
201*>
202* =====================================================================
203 SUBROUTINE sgebrd( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
204 \$ INFO )
205*
206* -- LAPACK computational routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 INTEGER INFO, LDA, LWORK, M, N
212* ..
213* .. Array Arguments ..
214 REAL A( LDA, * ), D( * ), E( * ), TAUP( * ),
215 \$ tauq( * ), work( * )
216* ..
217*
218* =====================================================================
219*
220* .. Parameters ..
221 REAL ONE
222 parameter( one = 1.0e+0 )
223* ..
224* .. Local Scalars ..
225 LOGICAL LQUERY
226 INTEGER I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
227 \$ nbmin, nx, ws
228* ..
229* .. External Subroutines ..
230 EXTERNAL sgebd2, sgemm, slabrd, xerbla
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC max, min
234* ..
235* .. External Functions ..
236 INTEGER ILAENV
237 REAL SROUNDUP_LWORK
238 EXTERNAL ilaenv, sroundup_lwork
239* ..
240* .. Executable Statements ..
241*
242* Test the input parameters
243*
244 info = 0
245 nb = max( 1, ilaenv( 1, 'SGEBRD', ' ', m, n, -1, -1 ) )
246 lwkopt = ( m+n )*nb
247 work( 1 ) = sroundup_lwork(lwkopt)
248 lquery = ( lwork.EQ.-1 )
249 IF( m.LT.0 ) THEN
250 info = -1
251 ELSE IF( n.LT.0 ) THEN
252 info = -2
253 ELSE IF( lda.LT.max( 1, m ) ) THEN
254 info = -4
255 ELSE IF( lwork.LT.max( 1, m, n ) .AND. .NOT.lquery ) THEN
256 info = -10
257 END IF
258 IF( info.LT.0 ) THEN
259 CALL xerbla( 'SGEBRD', -info )
260 RETURN
261 ELSE IF( lquery ) THEN
262 RETURN
263 END IF
264*
265* Quick return if possible
266*
267 minmn = min( m, n )
268 IF( minmn.EQ.0 ) THEN
269 work( 1 ) = 1
270 RETURN
271 END IF
272*
273 ws = max( m, n )
274 ldwrkx = m
275 ldwrky = n
276*
277 IF( nb.GT.1 .AND. nb.LT.minmn ) THEN
278*
279* Set the crossover point NX.
280*
281 nx = max( nb, ilaenv( 3, 'SGEBRD', ' ', m, n, -1, -1 ) )
282*
283* Determine when to switch from blocked to unblocked code.
284*
285 IF( nx.LT.minmn ) THEN
286 ws = ( m+n )*nb
287 IF( lwork.LT.ws ) THEN
288*
289* Not enough work space for the optimal NB, consider using
290* a smaller block size.
291*
292 nbmin = ilaenv( 2, 'SGEBRD', ' ', m, n, -1, -1 )
293 IF( lwork.GE.( m+n )*nbmin ) THEN
294 nb = lwork / ( m+n )
295 ELSE
296 nb = 1
297 nx = minmn
298 END IF
299 END IF
300 END IF
301 ELSE
302 nx = minmn
303 END IF
304*
305 DO 30 i = 1, minmn - nx, nb
306*
307* Reduce rows and columns i:i+nb-1 to bidiagonal form and return
308* the matrices X and Y which are needed to update the unreduced
309* part of the matrix
310*
311 CALL slabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),
312 \$ tauq( i ), taup( i ), work, ldwrkx,
313 \$ work( ldwrkx*nb+1 ), ldwrky )
314*
315* Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
316* of the form A := A - V*Y**T - X*U**T
317*
318 CALL sgemm( 'No transpose', 'Transpose', m-i-nb+1, n-i-nb+1,
319 \$ nb, -one, a( i+nb, i ), lda,
320 \$ work( ldwrkx*nb+nb+1 ), ldwrky, one,
321 \$ a( i+nb, i+nb ), lda )
322 CALL sgemm( 'No transpose', 'No transpose', m-i-nb+1, n-i-nb+1,
323 \$ nb, -one, work( nb+1 ), ldwrkx, a( i, i+nb ), lda,
324 \$ one, a( i+nb, i+nb ), lda )
325*
326* Copy diagonal and off-diagonal elements of B back into A
327*
328 IF( m.GE.n ) THEN
329 DO 10 j = i, i + nb - 1
330 a( j, j ) = d( j )
331 a( j, j+1 ) = e( j )
332 10 CONTINUE
333 ELSE
334 DO 20 j = i, i + nb - 1
335 a( j, j ) = d( j )
336 a( j+1, j ) = e( j )
337 20 CONTINUE
338 END IF
339 30 CONTINUE
340*
341* Use unblocked code to reduce the remainder of the matrix
342*
343 CALL sgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),
344 \$ tauq( i ), taup( i ), work, iinfo )
345 work( 1 ) = sroundup_lwork(ws)
346 RETURN
347*
348* End of SGEBRD
349*
350 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgebd2(m, n, a, lda, d, e, tauq, taup, work, info)
SGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
Definition sgebd2.f:189
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
Definition sgebrd.f:205
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slabrd(m, n, nb, a, lda, d, e, tauq, taup, x, ldx, y, ldy)
SLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
Definition slabrd.f:210