001:       SUBROUTINE CGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
002:      $                   INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       INTEGER            INFO, LDA, LWORK, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       REAL               D( * ), E( * )
014:       COMPLEX            A( LDA, * ), TAUP( * ), TAUQ( * ),
015:      $                   WORK( * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CGEBRD reduces a general complex M-by-N matrix A to upper or lower
022: *  bidiagonal form B by a unitary transformation: Q**H * A * P = B.
023: *
024: *  If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
025: *
026: *  Arguments
027: *  =========
028: *
029: *  M       (input) INTEGER
030: *          The number of rows in the matrix A.  M >= 0.
031: *
032: *  N       (input) INTEGER
033: *          The number of columns in the matrix A.  N >= 0.
034: *
035: *  A       (input/output) COMPLEX array, dimension (LDA,N)
036: *          On entry, the M-by-N general matrix to be reduced.
037: *          On exit,
038: *          if m >= n, the diagonal and the first superdiagonal are
039: *            overwritten with the upper bidiagonal matrix B; the
040: *            elements below the diagonal, with the array TAUQ, represent
041: *            the unitary matrix Q as a product of elementary
042: *            reflectors, and the elements above the first superdiagonal,
043: *            with the array TAUP, represent the unitary matrix P as
044: *            a product of elementary reflectors;
045: *          if m < n, the diagonal and the first subdiagonal are
046: *            overwritten with the lower bidiagonal matrix B; the
047: *            elements below the first subdiagonal, with the array TAUQ,
048: *            represent the unitary matrix Q as a product of
049: *            elementary reflectors, and the elements above the diagonal,
050: *            with the array TAUP, represent the unitary matrix P as
051: *            a product of elementary reflectors.
052: *          See Further Details.
053: *
054: *  LDA     (input) INTEGER
055: *          The leading dimension of the array A.  LDA >= max(1,M).
056: *
057: *  D       (output) REAL array, dimension (min(M,N))
058: *          The diagonal elements of the bidiagonal matrix B:
059: *          D(i) = A(i,i).
060: *
061: *  E       (output) REAL array, dimension (min(M,N)-1)
062: *          The off-diagonal elements of the bidiagonal matrix B:
063: *          if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
064: *          if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
065: *
066: *  TAUQ    (output) COMPLEX array dimension (min(M,N))
067: *          The scalar factors of the elementary reflectors which
068: *          represent the unitary matrix Q. See Further Details.
069: *
070: *  TAUP    (output) COMPLEX array, dimension (min(M,N))
071: *          The scalar factors of the elementary reflectors which
072: *          represent the unitary matrix P. See Further Details.
073: *
074: *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
075: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
076: *
077: *  LWORK   (input) INTEGER
078: *          The length of the array WORK.  LWORK >= max(1,M,N).
079: *          For optimum performance LWORK >= (M+N)*NB, where NB
080: *          is the optimal blocksize.
081: *
082: *          If LWORK = -1, then a workspace query is assumed; the routine
083: *          only calculates the optimal size of the WORK array, returns
084: *          this value as the first entry of the WORK array, and no error
085: *          message related to LWORK is issued by XERBLA.
086: *
087: *  INFO    (output) INTEGER
088: *          = 0:  successful exit.
089: *          < 0:  if INFO = -i, the i-th argument had an illegal value.
090: *
091: *  Further Details
092: *  ===============
093: *
094: *  The matrices Q and P are represented as products of elementary
095: *  reflectors:
096: *
097: *  If m >= n,
098: *
099: *     Q = H(1) H(2) . . . H(n)  and  P = G(1) G(2) . . . G(n-1)
100: *
101: *  Each H(i) and G(i) has the form:
102: *
103: *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
104: *
105: *  where tauq and taup are complex scalars, and v and u are complex
106: *  vectors; v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in
107: *  A(i+1:m,i); u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in
108: *  A(i,i+2:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
109: *
110: *  If m < n,
111: *
112: *     Q = H(1) H(2) . . . H(m-1)  and  P = G(1) G(2) . . . G(m)
113: *
114: *  Each H(i) and G(i) has the form:
115: *
116: *     H(i) = I - tauq * v * v'  and G(i) = I - taup * u * u'
117: *
118: *  where tauq and taup are complex scalars, and v and u are complex
119: *  vectors; v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in
120: *  A(i+2:m,i); u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in
121: *  A(i,i+1:n); tauq is stored in TAUQ(i) and taup in TAUP(i).
122: *
123: *  The contents of A on exit are illustrated by the following examples:
124: *
125: *  m = 6 and n = 5 (m > n):          m = 5 and n = 6 (m < n):
126: *
127: *    (  d   e   u1  u1  u1 )           (  d   u1  u1  u1  u1  u1 )
128: *    (  v1  d   e   u2  u2 )           (  e   d   u2  u2  u2  u2 )
129: *    (  v1  v2  d   e   u3 )           (  v1  e   d   u3  u3  u3 )
130: *    (  v1  v2  v3  d   e  )           (  v1  v2  e   d   u4  u4 )
131: *    (  v1  v2  v3  v4  d  )           (  v1  v2  v3  e   d   u5 )
132: *    (  v1  v2  v3  v4  v5 )
133: *
134: *  where d and e denote diagonal and off-diagonal elements of B, vi
135: *  denotes an element of the vector defining H(i), and ui an element of
136: *  the vector defining G(i).
137: *
138: *  =====================================================================
139: *
140: *     .. Parameters ..
141:       COMPLEX            ONE
142:       PARAMETER          ( ONE = ( 1.0E+0, 0.0E+0 ) )
143: *     ..
144: *     .. Local Scalars ..
145:       LOGICAL            LQUERY
146:       INTEGER            I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
147:      $                   NBMIN, NX
148:       REAL               WS
149: *     ..
150: *     .. External Subroutines ..
151:       EXTERNAL           CGEBD2, CGEMM, CLABRD, XERBLA
152: *     ..
153: *     .. Intrinsic Functions ..
154:       INTRINSIC          MAX, MIN, REAL
155: *     ..
156: *     .. External Functions ..
157:       INTEGER            ILAENV
158:       EXTERNAL           ILAENV
159: *     ..
160: *     .. Executable Statements ..
161: *
162: *     Test the input parameters
163: *
164:       INFO = 0
165:       NB = MAX( 1, ILAENV( 1, 'CGEBRD', ' ', M, N, -1, -1 ) )
166:       LWKOPT = ( M+N )*NB
167:       WORK( 1 ) = REAL( LWKOPT )
168:       LQUERY = ( LWORK.EQ.-1 )
169:       IF( M.LT.0 ) THEN
170:          INFO = -1
171:       ELSE IF( N.LT.0 ) THEN
172:          INFO = -2
173:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
174:          INFO = -4
175:       ELSE IF( LWORK.LT.MAX( 1, M, N ) .AND. .NOT.LQUERY ) THEN
176:          INFO = -10
177:       END IF
178:       IF( INFO.LT.0 ) THEN
179:          CALL XERBLA( 'CGEBRD', -INFO )
180:          RETURN
181:       ELSE IF( LQUERY ) THEN
182:          RETURN
183:       END IF
184: *
185: *     Quick return if possible
186: *
187:       MINMN = MIN( M, N )
188:       IF( MINMN.EQ.0 ) THEN
189:          WORK( 1 ) = 1
190:          RETURN
191:       END IF
192: *
193:       WS = MAX( M, N )
194:       LDWRKX = M
195:       LDWRKY = N
196: *
197:       IF( NB.GT.1 .AND. NB.LT.MINMN ) THEN
198: *
199: *        Set the crossover point NX.
200: *
201:          NX = MAX( NB, ILAENV( 3, 'CGEBRD', ' ', M, N, -1, -1 ) )
202: *
203: *        Determine when to switch from blocked to unblocked code.
204: *
205:          IF( NX.LT.MINMN ) THEN
206:             WS = ( M+N )*NB
207:             IF( LWORK.LT.WS ) THEN
208: *
209: *              Not enough work space for the optimal NB, consider using
210: *              a smaller block size.
211: *
212:                NBMIN = ILAENV( 2, 'CGEBRD', ' ', M, N, -1, -1 )
213:                IF( LWORK.GE.( M+N )*NBMIN ) THEN
214:                   NB = LWORK / ( M+N )
215:                ELSE
216:                   NB = 1
217:                   NX = MINMN
218:                END IF
219:             END IF
220:          END IF
221:       ELSE
222:          NX = MINMN
223:       END IF
224: *
225:       DO 30 I = 1, MINMN - NX, NB
226: *
227: *        Reduce rows and columns i:i+ib-1 to bidiagonal form and return
228: *        the matrices X and Y which are needed to update the unreduced
229: *        part of the matrix
230: *
231:          CALL CLABRD( M-I+1, N-I+1, NB, A( I, I ), LDA, D( I ), E( I ),
232:      $                TAUQ( I ), TAUP( I ), WORK, LDWRKX,
233:      $                WORK( LDWRKX*NB+1 ), LDWRKY )
234: *
235: *        Update the trailing submatrix A(i+ib:m,i+ib:n), using
236: *        an update of the form  A := A - V*Y' - X*U'
237: *
238:          CALL CGEMM( 'No transpose', 'Conjugate transpose', M-I-NB+1,
239:      $               N-I-NB+1, NB, -ONE, A( I+NB, I ), LDA,
240:      $               WORK( LDWRKX*NB+NB+1 ), LDWRKY, ONE,
241:      $               A( I+NB, I+NB ), LDA )
242:          CALL CGEMM( 'No transpose', 'No transpose', M-I-NB+1, N-I-NB+1,
243:      $               NB, -ONE, WORK( NB+1 ), LDWRKX, A( I, I+NB ), LDA,
244:      $               ONE, A( I+NB, I+NB ), LDA )
245: *
246: *        Copy diagonal and off-diagonal elements of B back into A
247: *
248:          IF( M.GE.N ) THEN
249:             DO 10 J = I, I + NB - 1
250:                A( J, J ) = D( J )
251:                A( J, J+1 ) = E( J )
252:    10       CONTINUE
253:          ELSE
254:             DO 20 J = I, I + NB - 1
255:                A( J, J ) = D( J )
256:                A( J+1, J ) = E( J )
257:    20       CONTINUE
258:          END IF
259:    30 CONTINUE
260: *
261: *     Use unblocked code to reduce the remainder of the matrix
262: *
263:       CALL CGEBD2( M-I+1, N-I+1, A( I, I ), LDA, D( I ), E( I ),
264:      $             TAUQ( I ), TAUP( I ), WORK, IINFO )
265:       WORK( 1 ) = WS
266:       RETURN
267: *
268: *     End of CGEBRD
269: *
270:       END
271: