LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgbbrd.f
Go to the documentation of this file.
1*> \brief \b DGBBRD
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGBBRD + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbbrd.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbbrd.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbbrd.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
22* LDQ, PT, LDPT, C, LDC, WORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER VECT
26* INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
30* $ PT( LDPT, * ), Q( LDQ, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DGBBRD reduces a real general m-by-n band matrix A to upper
40*> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
41*>
42*> The routine computes B, and optionally forms Q or P**T, or computes
43*> Q**T*C for a given matrix C.
44*> \endverbatim
45*
46* Arguments:
47* ==========
48*
49*> \param[in] VECT
50*> \verbatim
51*> VECT is CHARACTER*1
52*> Specifies whether or not the matrices Q and P**T are to be
53*> formed.
54*> = 'N': do not form Q or P**T;
55*> = 'Q': form Q only;
56*> = 'P': form P**T only;
57*> = 'B': form both.
58*> \endverbatim
59*>
60*> \param[in] M
61*> \verbatim
62*> M is INTEGER
63*> The number of rows of the matrix A. M >= 0.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*> N is INTEGER
69*> The number of columns of the matrix A. N >= 0.
70*> \endverbatim
71*>
72*> \param[in] NCC
73*> \verbatim
74*> NCC is INTEGER
75*> The number of columns of the matrix C. NCC >= 0.
76*> \endverbatim
77*>
78*> \param[in] KL
79*> \verbatim
80*> KL is INTEGER
81*> The number of subdiagonals of the matrix A. KL >= 0.
82*> \endverbatim
83*>
84*> \param[in] KU
85*> \verbatim
86*> KU is INTEGER
87*> The number of superdiagonals of the matrix A. KU >= 0.
88*> \endverbatim
89*>
90*> \param[in,out] AB
91*> \verbatim
92*> AB is DOUBLE PRECISION array, dimension (LDAB,N)
93*> On entry, the m-by-n band matrix A, stored in rows 1 to
94*> KL+KU+1. The j-th column of A is stored in the j-th column of
95*> the array AB as follows:
96*> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
97*> On exit, A is overwritten by values generated during the
98*> reduction.
99*> \endverbatim
100*>
101*> \param[in] LDAB
102*> \verbatim
103*> LDAB is INTEGER
104*> The leading dimension of the array A. LDAB >= KL+KU+1.
105*> \endverbatim
106*>
107*> \param[out] D
108*> \verbatim
109*> D is DOUBLE PRECISION array, dimension (min(M,N))
110*> The diagonal elements of the bidiagonal matrix B.
111*> \endverbatim
112*>
113*> \param[out] E
114*> \verbatim
115*> E is DOUBLE PRECISION array, dimension (min(M,N)-1)
116*> The superdiagonal elements of the bidiagonal matrix B.
117*> \endverbatim
118*>
119*> \param[out] Q
120*> \verbatim
121*> Q is DOUBLE PRECISION array, dimension (LDQ,M)
122*> If VECT = 'Q' or 'B', the m-by-m orthogonal matrix Q.
123*> If VECT = 'N' or 'P', the array Q is not referenced.
124*> \endverbatim
125*>
126*> \param[in] LDQ
127*> \verbatim
128*> LDQ is INTEGER
129*> The leading dimension of the array Q.
130*> LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
131*> \endverbatim
132*>
133*> \param[out] PT
134*> \verbatim
135*> PT is DOUBLE PRECISION array, dimension (LDPT,N)
136*> If VECT = 'P' or 'B', the n-by-n orthogonal matrix P'.
137*> If VECT = 'N' or 'Q', the array PT is not referenced.
138*> \endverbatim
139*>
140*> \param[in] LDPT
141*> \verbatim
142*> LDPT is INTEGER
143*> The leading dimension of the array PT.
144*> LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
145*> \endverbatim
146*>
147*> \param[in,out] C
148*> \verbatim
149*> C is DOUBLE PRECISION array, dimension (LDC,NCC)
150*> On entry, an m-by-ncc matrix C.
151*> On exit, C is overwritten by Q**T*C.
152*> C is not referenced if NCC = 0.
153*> \endverbatim
154*>
155*> \param[in] LDC
156*> \verbatim
157*> LDC is INTEGER
158*> The leading dimension of the array C.
159*> LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
160*> \endverbatim
161*>
162*> \param[out] WORK
163*> \verbatim
164*> WORK is DOUBLE PRECISION array, dimension (2*max(M,N))
165*> \endverbatim
166*>
167*> \param[out] INFO
168*> \verbatim
169*> INFO is INTEGER
170*> = 0: successful exit.
171*> < 0: if INFO = -i, the i-th argument had an illegal value.
172*> \endverbatim
173*
174* Authors:
175* ========
176*
177*> \author Univ. of Tennessee
178*> \author Univ. of California Berkeley
179*> \author Univ. of Colorado Denver
180*> \author NAG Ltd.
181*
182*> \ingroup gbbrd
183*
184* =====================================================================
185 SUBROUTINE dgbbrd( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
186 $ LDQ, PT, LDPT, C, LDC, WORK, INFO )
187*
188* -- LAPACK computational routine --
189* -- LAPACK is a software package provided by Univ. of Tennessee, --
190* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
191*
192* .. Scalar Arguments ..
193 CHARACTER VECT
194 INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
195* ..
196* .. Array Arguments ..
197 DOUBLE PRECISION AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
198 $ pt( ldpt, * ), q( ldq, * ), work( * )
199* ..
200*
201* =====================================================================
202*
203* .. Parameters ..
204 DOUBLE PRECISION ZERO, ONE
205 parameter( zero = 0.0d+0, one = 1.0d+0 )
206* ..
207* .. Local Scalars ..
208 LOGICAL WANTB, WANTC, WANTPT, WANTQ
209 INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
210 $ kun, l, minmn, ml, ml0, mn, mu, mu0, nr, nrt
211 DOUBLE PRECISION RA, RB, RC, RS
212* ..
213* .. External Subroutines ..
214 EXTERNAL dlargv, dlartg, dlartv, dlaset, drot, xerbla
215* ..
216* .. Intrinsic Functions ..
217 INTRINSIC max, min
218* ..
219* .. External Functions ..
220 LOGICAL LSAME
221 EXTERNAL lsame
222* ..
223* .. Executable Statements ..
224*
225* Test the input parameters
226*
227 wantb = lsame( vect, 'B' )
228 wantq = lsame( vect, 'Q' ) .OR. wantb
229 wantpt = lsame( vect, 'P' ) .OR. wantb
230 wantc = ncc.GT.0
231 klu1 = kl + ku + 1
232 info = 0
233 IF( .NOT.wantq .AND. .NOT.wantpt .AND. .NOT.lsame( vect, 'N' ) )
234 $ THEN
235 info = -1
236 ELSE IF( m.LT.0 ) THEN
237 info = -2
238 ELSE IF( n.LT.0 ) THEN
239 info = -3
240 ELSE IF( ncc.LT.0 ) THEN
241 info = -4
242 ELSE IF( kl.LT.0 ) THEN
243 info = -5
244 ELSE IF( ku.LT.0 ) THEN
245 info = -6
246 ELSE IF( ldab.LT.klu1 ) THEN
247 info = -8
248 ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
249 info = -12
250 ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
251 info = -14
252 ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
253 info = -16
254 END IF
255 IF( info.NE.0 ) THEN
256 CALL xerbla( 'DGBBRD', -info )
257 RETURN
258 END IF
259*
260* Initialize Q and P**T to the unit matrix, if needed
261*
262 IF( wantq )
263 $ CALL dlaset( 'Full', m, m, zero, one, q, ldq )
264 IF( wantpt )
265 $ CALL dlaset( 'Full', n, n, zero, one, pt, ldpt )
266*
267* Quick return if possible.
268*
269 IF( m.EQ.0 .OR. n.EQ.0 )
270 $ RETURN
271*
272 minmn = min( m, n )
273*
274 IF( kl+ku.GT.1 ) THEN
275*
276* Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
277* first to lower bidiagonal form and then transform to upper
278* bidiagonal
279*
280 IF( ku.GT.0 ) THEN
281 ml0 = 1
282 mu0 = 2
283 ELSE
284 ml0 = 2
285 mu0 = 1
286 END IF
287*
288* Wherever possible, plane rotations are generated and applied in
289* vector operations of length NR over the index set J1:J2:KLU1.
290*
291* The sines of the plane rotations are stored in WORK(1:max(m,n))
292* and the cosines in WORK(max(m,n)+1:2*max(m,n)).
293*
294 mn = max( m, n )
295 klm = min( m-1, kl )
296 kun = min( n-1, ku )
297 kb = klm + kun
298 kb1 = kb + 1
299 inca = kb1*ldab
300 nr = 0
301 j1 = klm + 2
302 j2 = 1 - kun
303*
304 DO 90 i = 1, minmn
305*
306* Reduce i-th column and i-th row of matrix to bidiagonal form
307*
308 ml = klm + 1
309 mu = kun + 1
310 DO 80 kk = 1, kb
311 j1 = j1 + kb
312 j2 = j2 + kb
313*
314* generate plane rotations to annihilate nonzero elements
315* which have been created below the band
316*
317 IF( nr.GT.0 )
318 $ CALL dlargv( nr, ab( klu1, j1-klm-1 ), inca,
319 $ work( j1 ), kb1, work( mn+j1 ), kb1 )
320*
321* apply plane rotations from the left
322*
323 DO 10 l = 1, kb
324 IF( j2-klm+l-1.GT.n ) THEN
325 nrt = nr - 1
326 ELSE
327 nrt = nr
328 END IF
329 IF( nrt.GT.0 )
330 $ CALL dlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,
331 $ ab( klu1-l+1, j1-klm+l-1 ), inca,
332 $ work( mn+j1 ), work( j1 ), kb1 )
333 10 CONTINUE
334*
335 IF( ml.GT.ml0 ) THEN
336 IF( ml.LE.m-i+1 ) THEN
337*
338* generate plane rotation to annihilate a(i+ml-1,i)
339* within the band, and apply rotation from the left
340*
341 CALL dlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
342 $ work( mn+i+ml-1 ), work( i+ml-1 ),
343 $ ra )
344 ab( ku+ml-1, i ) = ra
345 IF( i.LT.n )
346 $ CALL drot( min( ku+ml-2, n-i ),
347 $ ab( ku+ml-2, i+1 ), ldab-1,
348 $ ab( ku+ml-1, i+1 ), ldab-1,
349 $ work( mn+i+ml-1 ), work( i+ml-1 ) )
350 END IF
351 nr = nr + 1
352 j1 = j1 - kb1
353 END IF
354*
355 IF( wantq ) THEN
356*
357* accumulate product of plane rotations in Q
358*
359 DO 20 j = j1, j2, kb1
360 CALL drot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
361 $ work( mn+j ), work( j ) )
362 20 CONTINUE
363 END IF
364*
365 IF( wantc ) THEN
366*
367* apply plane rotations to C
368*
369 DO 30 j = j1, j2, kb1
370 CALL drot( ncc, c( j-1, 1 ), ldc, c( j, 1 ), ldc,
371 $ work( mn+j ), work( j ) )
372 30 CONTINUE
373 END IF
374*
375 IF( j2+kun.GT.n ) THEN
376*
377* adjust J2 to keep within the bounds of the matrix
378*
379 nr = nr - 1
380 j2 = j2 - kb1
381 END IF
382*
383 DO 40 j = j1, j2, kb1
384*
385* create nonzero element a(j-1,j+ku) above the band
386* and store it in WORK(n+1:2*n)
387*
388 work( j+kun ) = work( j )*ab( 1, j+kun )
389 ab( 1, j+kun ) = work( mn+j )*ab( 1, j+kun )
390 40 CONTINUE
391*
392* generate plane rotations to annihilate nonzero elements
393* which have been generated above the band
394*
395 IF( nr.GT.0 )
396 $ CALL dlargv( nr, ab( 1, j1+kun-1 ), inca,
397 $ work( j1+kun ), kb1, work( mn+j1+kun ),
398 $ kb1 )
399*
400* apply plane rotations from the right
401*
402 DO 50 l = 1, kb
403 IF( j2+l-1.GT.m ) THEN
404 nrt = nr - 1
405 ELSE
406 nrt = nr
407 END IF
408 IF( nrt.GT.0 )
409 $ CALL dlartv( nrt, ab( l+1, j1+kun-1 ), inca,
410 $ ab( l, j1+kun ), inca,
411 $ work( mn+j1+kun ), work( j1+kun ),
412 $ kb1 )
413 50 CONTINUE
414*
415 IF( ml.EQ.ml0 .AND. mu.GT.mu0 ) THEN
416 IF( mu.LE.n-i+1 ) THEN
417*
418* generate plane rotation to annihilate a(i,i+mu-1)
419* within the band, and apply rotation from the right
420*
421 CALL dlartg( ab( ku-mu+3, i+mu-2 ),
422 $ ab( ku-mu+2, i+mu-1 ),
423 $ work( mn+i+mu-1 ), work( i+mu-1 ),
424 $ ra )
425 ab( ku-mu+3, i+mu-2 ) = ra
426 CALL drot( min( kl+mu-2, m-i ),
427 $ ab( ku-mu+4, i+mu-2 ), 1,
428 $ ab( ku-mu+3, i+mu-1 ), 1,
429 $ work( mn+i+mu-1 ), work( i+mu-1 ) )
430 END IF
431 nr = nr + 1
432 j1 = j1 - kb1
433 END IF
434*
435 IF( wantpt ) THEN
436*
437* accumulate product of plane rotations in P**T
438*
439 DO 60 j = j1, j2, kb1
440 CALL drot( n, pt( j+kun-1, 1 ), ldpt,
441 $ pt( j+kun, 1 ), ldpt, work( mn+j+kun ),
442 $ work( j+kun ) )
443 60 CONTINUE
444 END IF
445*
446 IF( j2+kb.GT.m ) THEN
447*
448* adjust J2 to keep within the bounds of the matrix
449*
450 nr = nr - 1
451 j2 = j2 - kb1
452 END IF
453*
454 DO 70 j = j1, j2, kb1
455*
456* create nonzero element a(j+kl+ku,j+ku-1) below the
457* band and store it in WORK(1:n)
458*
459 work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
460 ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
461 70 CONTINUE
462*
463 IF( ml.GT.ml0 ) THEN
464 ml = ml - 1
465 ELSE
466 mu = mu - 1
467 END IF
468 80 CONTINUE
469 90 CONTINUE
470 END IF
471*
472 IF( ku.EQ.0 .AND. kl.GT.0 ) THEN
473*
474* A has been reduced to lower bidiagonal form
475*
476* Transform lower bidiagonal form to upper bidiagonal by applying
477* plane rotations from the left, storing diagonal elements in D
478* and off-diagonal elements in E
479*
480 DO 100 i = 1, min( m-1, n )
481 CALL dlartg( ab( 1, i ), ab( 2, i ), rc, rs, ra )
482 d( i ) = ra
483 IF( i.LT.n ) THEN
484 e( i ) = rs*ab( 1, i+1 )
485 ab( 1, i+1 ) = rc*ab( 1, i+1 )
486 END IF
487 IF( wantq )
488 $ CALL drot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc, rs )
489 IF( wantc )
490 $ CALL drot( ncc, c( i, 1 ), ldc, c( i+1, 1 ), ldc, rc,
491 $ rs )
492 100 CONTINUE
493 IF( m.LE.n )
494 $ d( m ) = ab( 1, m )
495 ELSE IF( ku.GT.0 ) THEN
496*
497* A has been reduced to upper bidiagonal form
498*
499 IF( m.LT.n ) THEN
500*
501* Annihilate a(m,m+1) by applying plane rotations from the
502* right, storing diagonal elements in D and off-diagonal
503* elements in E
504*
505 rb = ab( ku, m+1 )
506 DO 110 i = m, 1, -1
507 CALL dlartg( ab( ku+1, i ), rb, rc, rs, ra )
508 d( i ) = ra
509 IF( i.GT.1 ) THEN
510 rb = -rs*ab( ku, i )
511 e( i-1 ) = rc*ab( ku, i )
512 END IF
513 IF( wantpt )
514 $ CALL drot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
515 $ rc, rs )
516 110 CONTINUE
517 ELSE
518*
519* Copy off-diagonal elements to E and diagonal elements to D
520*
521 DO 120 i = 1, minmn - 1
522 e( i ) = ab( ku, i+1 )
523 120 CONTINUE
524 DO 130 i = 1, minmn
525 d( i ) = ab( ku+1, i )
526 130 CONTINUE
527 END IF
528 ELSE
529*
530* A is diagonal. Set elements of E to zero and copy diagonal
531* elements to D.
532*
533 DO 140 i = 1, minmn - 1
534 e( i ) = zero
535 140 CONTINUE
536 DO 150 i = 1, minmn
537 d( i ) = ab( 1, i )
538 150 CONTINUE
539 END IF
540 RETURN
541*
542* End of DGBBRD
543*
544 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgbbrd(vect, m, n, ncc, kl, ku, ab, ldab, d, e, q, ldq, pt, ldpt, c, ldc, work, info)
DGBBRD
Definition dgbbrd.f:187
subroutine dlargv(n, x, incx, y, incy, c, incc)
DLARGV generates a vector of plane rotations with real cosines and real sines.
Definition dlargv.f:104
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
subroutine dlartv(n, x, incx, y, incy, c, s, incc)
DLARTV applies a vector of plane rotations with real cosines and real sines to the elements of a pair...
Definition dlartv.f:108
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92