KTH framework for Nek5000 toolboxes; testing version  0.0.1
dgebrd.f
Go to the documentation of this file.
1  SUBROUTINE dgebrd( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
2  $ INFO )
3 *
4 * -- LAPACK routine (version 3.0) --
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6 * Courant Institute, Argonne National Lab, and Rice University
7 * June 30, 1999
8 *
9 * .. Scalar Arguments ..
10  INTEGER INFO, LDA, LWORK, M, N
11 * ..
12 * .. Array Arguments ..
13  DOUBLE PRECISION A( LDA, * ), D( * ), E( * ), TAUP( * ),
14  $ tauq( * ), work( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DGEBRD reduces a general real M-by-N matrix A to upper or lower
21 * bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
22 *
23 * If m >= n, B is upper bidiagonal; if m < n, B is lower bidiagonal.
24 *
25 * Arguments
26 * =========
27 *
28 * M (input) INTEGER
29 * The number of rows in the matrix A. M >= 0.
30 *
31 * N (input) INTEGER
32 * The number of columns in the matrix A. N >= 0.
33 *
34 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
35 * On entry, the M-by-N general matrix to be reduced.
36 * On exit,
37 * if m >= n, the diagonal and the first superdiagonal are
38 * overwritten with the upper bidiagonal matrix B; the
39 * elements below the diagonal, with the array TAUQ, represent
40 * the orthogonal matrix Q as a product of elementary
41 * reflectors, and the elements above the first superdiagonal,
42 * with the array TAUP, represent the orthogonal matrix P as
43 * a product of elementary reflectors;
44 * if m < n, the diagonal and the first subdiagonal are
45 * overwritten with the lower bidiagonal matrix B; the
46 * elements below the first subdiagonal, with the array TAUQ,
47 * represent the orthogonal matrix Q as a product of
48 * elementary reflectors, and the elements above the diagonal,
49 * with the array TAUP, represent the orthogonal matrix P as
50 * a product of elementary reflectors.
51 * See Further Details.
52 *
53 * LDA (input) INTEGER
54 * The leading dimension of the array A. LDA >= max(1,M).
55 *
56 * D (output) DOUBLE PRECISION array, dimension (min(M,N))
57 * The diagonal elements of the bidiagonal matrix B:
58 * D(i) = A(i,i).
59 *
60 * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
61 * The off-diagonal elements of the bidiagonal matrix B:
62 * if m >= n, E(i) = A(i,i+1) for i = 1,2,...,n-1;
63 * if m < n, E(i) = A(i+1,i) for i = 1,2,...,m-1.
64 *
65 * TAUQ (output) DOUBLE PRECISION array dimension (min(M,N))
66 * The scalar factors of the elementary reflectors which
67 * represent the orthogonal matrix Q. See Further Details.
68 *
69 * TAUP (output) DOUBLE PRECISION array, dimension (min(M,N))
70 * The scalar factors of the elementary reflectors which
71 * represent the orthogonal matrix P. See Further Details.
72 *
73 * WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
74 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
75 *
76 * LWORK (input) INTEGER
77 * The length of the array WORK. LWORK >= max(1,M,N).
78 * For optimum performance LWORK >= (M+N)*NB, where NB
79 * is the optimal blocksize.
80 *
81 * If LWORK = -1, then a workspace query is assumed; the routine
82 * only calculates the optimal size of the WORK array, returns
83 * this value as the first entry of the WORK array, and no error
84 * message related to LWORK is issued by XERBLA.
85 *
86 * INFO (output) INTEGER
87 * = 0: successful exit
88 * < 0: if INFO = -i, the i-th argument had an illegal value.
89 *
90 * Further Details
91 * ===============
92 *
93 * The matrices Q and P are represented as products of elementary
94 * reflectors:
95 *
96 * If m >= n,
97 *
98 * Q = H(1) H(2) . . . H(n) and P = G(1) G(2) . . . G(n-1)
99 *
100 * Each H(i) and G(i) has the form:
101 *
102 * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
103 *
104 * where tauq and taup are real scalars, and v and u are real vectors;
105 * v(1:i-1) = 0, v(i) = 1, and v(i+1:m) is stored on exit in A(i+1:m,i);
106 * u(1:i) = 0, u(i+1) = 1, and u(i+2:n) is stored on exit in A(i,i+2:n);
107 * tauq is stored in TAUQ(i) and taup in TAUP(i).
108 *
109 * If m < n,
110 *
111 * Q = H(1) H(2) . . . H(m-1) and P = G(1) G(2) . . . G(m)
112 *
113 * Each H(i) and G(i) has the form:
114 *
115 * H(i) = I - tauq * v * v' and G(i) = I - taup * u * u'
116 *
117 * where tauq and taup are real scalars, and v and u are real vectors;
118 * v(1:i) = 0, v(i+1) = 1, and v(i+2:m) is stored on exit in A(i+2:m,i);
119 * u(1:i-1) = 0, u(i) = 1, and u(i+1:n) is stored on exit in A(i,i+1:n);
120 * tauq is stored in TAUQ(i) and taup in TAUP(i).
121 *
122 * The contents of A on exit are illustrated by the following examples:
123 *
124 * m = 6 and n = 5 (m > n): m = 5 and n = 6 (m < n):
125 *
126 * ( d e u1 u1 u1 ) ( d u1 u1 u1 u1 u1 )
127 * ( v1 d e u2 u2 ) ( e d u2 u2 u2 u2 )
128 * ( v1 v2 d e u3 ) ( v1 e d u3 u3 u3 )
129 * ( v1 v2 v3 d e ) ( v1 v2 e d u4 u4 )
130 * ( v1 v2 v3 v4 d ) ( v1 v2 v3 e d u5 )
131 * ( v1 v2 v3 v4 v5 )
132 *
133 * where d and e denote diagonal and off-diagonal elements of B, vi
134 * denotes an element of the vector defining H(i), and ui an element of
135 * the vector defining G(i).
136 *
137 * =====================================================================
138 *
139 * .. Parameters ..
140  DOUBLE PRECISION ONE
141  parameter( one = 1.0d+0 )
142 * ..
143 * .. Local Scalars ..
144  LOGICAL LQUERY
145  INTEGER I, IINFO, J, LDWRKX, LDWRKY, LWKOPT, MINMN, NB,
146  $ nbmin, nx
147  DOUBLE PRECISION WS
148 * ..
149 * .. External Subroutines ..
150  EXTERNAL dgebd2, dgemm, dlabrd, xerbla
151 * ..
152 * .. Intrinsic Functions ..
153  INTRINSIC dble, max, min
154 * ..
155 * .. External Functions ..
156  INTEGER ILAENV
157  EXTERNAL ilaenv
158 * ..
159 * .. Executable Statements ..
160 *
161 * Test the input parameters
162 *
163  info = 0
164  nb = max( 1, ilaenv( 1, 'DGEBRD', ' ', m, n, -1, -1 ) )
165  lwkopt = ( m+n )*nb
166  work( 1 ) = dble( lwkopt )
167  lquery = ( lwork.EQ.-1 )
168  IF( m.LT.0 ) THEN
169  info = -1
170  ELSE IF( n.LT.0 ) THEN
171  info = -2
172  ELSE IF( lda.LT.max( 1, m ) ) THEN
173  info = -4
174  ELSE IF( lwork.LT.max( 1, m, n ) .AND. .NOT.lquery ) THEN
175  info = -10
176  END IF
177  IF( info.LT.0 ) THEN
178  CALL xerbla( 'DGEBRD', -info )
179  RETURN
180  ELSE IF( lquery ) THEN
181  RETURN
182  END IF
183 *
184 * Quick return if possible
185 *
186  minmn = min( m, n )
187  IF( minmn.EQ.0 ) THEN
188  work( 1 ) = 1
189  RETURN
190  END IF
191 *
192  ws = max( m, n )
193  ldwrkx = m
194  ldwrky = n
195 *
196  IF( nb.GT.1 .AND. nb.LT.minmn ) THEN
197 *
198 * Set the crossover point NX.
199 *
200  nx = max( nb, ilaenv( 3, 'DGEBRD', ' ', m, n, -1, -1 ) )
201 *
202 * Determine when to switch from blocked to unblocked code.
203 *
204  IF( nx.LT.minmn ) THEN
205  ws = ( m+n )*nb
206  IF( lwork.LT.ws ) THEN
207 *
208 * Not enough work space for the optimal NB, consider using
209 * a smaller block size.
210 *
211  nbmin = ilaenv( 2, 'DGEBRD', ' ', m, n, -1, -1 )
212  IF( lwork.GE.( m+n )*nbmin ) THEN
213  nb = lwork / ( m+n )
214  ELSE
215  nb = 1
216  nx = minmn
217  END IF
218  END IF
219  END IF
220  ELSE
221  nx = minmn
222  END IF
223 *
224  DO 30 i = 1, minmn - nx, nb
225 *
226 * Reduce rows and columns i:i+nb-1 to bidiagonal form and return
227 * the matrices X and Y which are needed to update the unreduced
228 * part of the matrix
229 *
230  CALL dlabrd( m-i+1, n-i+1, nb, a( i, i ), lda, d( i ), e( i ),
231  $ tauq( i ), taup( i ), work, ldwrkx,
232  $ work( ldwrkx*nb+1 ), ldwrky )
233 *
234 * Update the trailing submatrix A(i+nb:m,i+nb:n), using an update
235 * of the form A := A - V*Y' - X*U'
236 *
237  CALL dgemm( 'No transpose', 'Transpose', m-i-nb+1, n-i-nb+1,
238  $ nb, -one, a( i+nb, i ), lda,
239  $ work( ldwrkx*nb+nb+1 ), ldwrky, one,
240  $ a( i+nb, i+nb ), lda )
241  CALL dgemm( 'No transpose', 'No transpose', m-i-nb+1, n-i-nb+1,
242  $ nb, -one, work( nb+1 ), ldwrkx, a( i, i+nb ), lda,
243  $ one, a( i+nb, i+nb ), lda )
244 *
245 * Copy diagonal and off-diagonal elements of B back into A
246 *
247  IF( m.GE.n ) THEN
248  DO 10 j = i, i + nb - 1
249  a( j, j ) = d( j )
250  a( j, j+1 ) = e( j )
251  10 CONTINUE
252  ELSE
253  DO 20 j = i, i + nb - 1
254  a( j, j ) = d( j )
255  a( j+1, j ) = e( j )
256  20 CONTINUE
257  END IF
258  30 CONTINUE
259 *
260 * Use unblocked code to reduce the remainder of the matrix
261 *
262  CALL dgebd2( m-i+1, n-i+1, a( i, i ), lda, d( i ), e( i ),
263  $ tauq( i ), taup( i ), work, iinfo )
264  work( 1 ) = ws
265  RETURN
266 *
267 * End of DGEBRD
268 *
269  END
subroutine dgebd2(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, INFO)
Definition: dgebd2.f:2
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
Definition: dgebrd.f:3
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgemm.f:3
subroutine dlabrd(M, N, NB, A, LDA, D, E, TAUQ, TAUP, X, LDX, Y, LDY)
Definition: dlabrd.f:3
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2