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