LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgebrd.f
Go to the documentation of this file.
1*> \brief \b DGEBRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGEBRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgebrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgebrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgebrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGEBRD( 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* DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ),
29* $ TAUQ( * ), WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DGEBRD 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dgebrd( 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 DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ),
215 $ tauq( * ), work( * )
216* ..
217*
218* =====================================================================
219*
220* .. Parameters ..
221 DOUBLE PRECISION ONE
222 parameter( one = 1.0d+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 dgebd2, dgemm, dlabrd, xerbla
231* ..
232* .. Intrinsic Functions ..
233 INTRINSIC dble, max, min
234* ..
235* .. External Functions ..
236 INTEGER ILAENV
237 EXTERNAL ilaenv
238* ..
239* .. Executable Statements ..
240*
241* Test the input parameters
242*
243 info = 0
244 nb = max( 1, ilaenv( 1, 'DGEBRD', ' ', m, n, -1, -1 ) )
245 lwkopt = ( m+n )*nb
246 work( 1 ) = dble( lwkopt )
247 lquery = ( lwork.EQ.-1 )
248 IF( m.LT.0 ) THEN
249 info = -1
250 ELSE IF( n.LT.0 ) THEN
251 info = -2
252 ELSE IF( lda.LT.max( 1, m ) ) THEN
253 info = -4
254 ELSE IF( lwork.LT.max( 1, m, n ) .AND. .NOT.lquery ) THEN
255 info = -10
256 END IF
257 IF( info.LT.0 ) THEN
258 CALL xerbla( 'DGEBRD', -info )
259 RETURN
260 ELSE IF( lquery ) THEN
261 RETURN
262 END IF
263*
264* Quick return if possible
265*
266 minmn = min( m, n )
267 IF( minmn.EQ.0 ) THEN
268 work( 1 ) = 1
269 RETURN
270 END IF
271*
272 ws = max( m, n )
273 ldwrkx = m
274 ldwrky = n
275*
276 IF( nb.GT.1 .AND. nb.LT.minmn ) THEN
277*
278* Set the crossover point NX.
279*
280 nx = max( nb, ilaenv( 3, 'DGEBRD', ' ', m, n, -1, -1 ) )
281*
282* Determine when to switch from blocked to unblocked code.
283*
284 IF( nx.LT.minmn ) THEN
285 ws = ( m+n )*nb
286 IF( lwork.LT.ws ) THEN
287*
288* Not enough work space for the optimal NB, consider using
289* a smaller block size.
290*
291 nbmin = ilaenv( 2, 'DGEBRD', ' ', m, n, -1, -1 )
292 IF( lwork.GE.( m+n )*nbmin ) THEN
293 nb = lwork / ( m+n )
294 ELSE
295 nb = 1
296 nx = minmn
297 END IF
298 END IF
299 END IF
300 ELSE
301 nx = minmn
302 END IF
303*
304 DO 30 i = 1, minmn - nx, nb
305*
306* Reduce rows and columns i:i+nb-1 to bidiagonal form and return
307* the matrices X and Y which are needed to update the unreduced
308* part of the matrix
309*
310 CALL dlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),
311 $ tauq( i ), taup( i ), work, ldwrkx,
312 $ work( ldwrkx*nb+1 ), ldwrky )
313*
314* Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
315* of the form A := A - V*Y**T - X*U**T
316*
317 CALL dgemm( 'No transpose', 'Transpose', m-i-nb+1, n-i-nb+1,
318 $ nb, -one, a( i+nb, i ), lda,
319 $ work( ldwrkx*nb+nb+1 ), ldwrky, one,
320 $ a( i+nb, i+nb ), lda )
321 CALL dgemm( 'No transpose', 'No transpose', m-i-nb+1, n-i-nb+1,
322 $ nb, -one, work( nb+1 ), ldwrkx, a( i, i+nb ), lda,
323 $ one, a( i+nb, i+nb ), lda )
324*
325* Copy diagonal and off-diagonal elements of B back into A
326*
327 IF( m.GE.n ) THEN
328 DO 10 j = i, i + nb - 1
329 a( j, j ) = d( j )
330 a( j, j+1 ) = e( j )
331 10 CONTINUE
332 ELSE
333 DO 20 j = i, i + nb - 1
334 a( j, j ) = d( j )
335 a( j+1, j ) = e( j )
336 20 CONTINUE
337 END IF
338 30 CONTINUE
339*
340* Use unblocked code to reduce the remainder of the matrix
341*
342 CALL dgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),
343 $ tauq( i ), taup( i ), work, iinfo )
344 work( 1 ) = ws
345 RETURN
346*
347* End of DGEBRD
348*
349 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgebd2(m, n, a, lda, d, e, tauq, taup, work, info)
DGEBD2 reduces a general matrix to bidiagonal form using an unblocked algorithm.
Definition dgebd2.f:189
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:205
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlabrd(m, n, nb, a, lda, d, e, tauq, taup, x, ldx, y, ldy)
DLABRD reduces the first nb rows and columns of a general matrix to a bidiagonal form.
Definition dlabrd.f:210