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