LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgbbrd()

subroutine dgbbrd ( character vect,
integer m,
integer n,
integer ncc,
integer kl,
integer ku,
double precision, dimension( ldab, * ) ab,
integer ldab,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( ldpt, * ) pt,
integer ldpt,
double precision, dimension( ldc, * ) c,
integer ldc,
double precision, dimension( * ) work,
integer info )

DGBBRD

Download DGBBRD + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGBBRD reduces a real general m-by-n band matrix A to upper
!> bidiagonal form B by an orthogonal transformation: Q**T * A * P = B.
!>
!> The routine computes B, and optionally forms Q or P**T, or computes
!> Q**T*C for a given matrix C.
!>
Parameters
[in]VECT
!>          VECT is CHARACTER*1
!>          Specifies whether or not the matrices Q and P**T are to be
!>          formed.
!>          = 'N': do not form Q or P**T;
!>          = 'Q': form Q only;
!>          = 'P': form P**T only;
!>          = 'B': form both.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the matrix A.  N >= 0.
!>
[in]NCC
!>          NCC is INTEGER
!>          The number of columns of the matrix C.  NCC >= 0.
!>
[in]KL
!>          KL is INTEGER
!>          The number of subdiagonals of the matrix A. KL >= 0.
!>
[in]KU
!>          KU is INTEGER
!>          The number of superdiagonals of the matrix A. KU >= 0.
!>
[in,out]AB
!>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
!>          On entry, the m-by-n band matrix A, stored in rows 1 to
!>          KL+KU+1. The j-th column of A is stored in the j-th column of
!>          the array AB as follows:
!>          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
!>          On exit, A is overwritten by values generated during the
!>          reduction.
!>
[in]LDAB
!>          LDAB is INTEGER
!>          The leading dimension of the array A. LDAB >= KL+KU+1.
!>
[out]D
!>          D is DOUBLE PRECISION array, dimension (min(M,N))
!>          The diagonal elements of the bidiagonal matrix B.
!>
[out]E
!>          E is DOUBLE PRECISION array, dimension (min(M,N)-1)
!>          The superdiagonal elements of the bidiagonal matrix B.
!>
[out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ,M)
!>          If VECT = 'Q' or 'B', the m-by-m orthogonal matrix Q.
!>          If VECT = 'N' or 'P', the array Q is not referenced.
!>
[in]LDQ
!>          LDQ is INTEGER
!>          The leading dimension of the array Q.
!>          LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
!>
[out]PT
!>          PT is DOUBLE PRECISION array, dimension (LDPT,N)
!>          If VECT = 'P' or 'B', the n-by-n orthogonal matrix P'.
!>          If VECT = 'N' or 'Q', the array PT is not referenced.
!>
[in]LDPT
!>          LDPT is INTEGER
!>          The leading dimension of the array PT.
!>          LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
!>
[in,out]C
!>          C is DOUBLE PRECISION array, dimension (LDC,NCC)
!>          On entry, an m-by-ncc matrix C.
!>          On exit, C is overwritten by Q**T*C.
!>          C is not referenced if NCC = 0.
!>
[in]LDC
!>          LDC is INTEGER
!>          The leading dimension of the array C.
!>          LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (2*max(M,N))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 183 of file dgbbrd.f.

185*
186* -- LAPACK computational routine --
187* -- LAPACK is a software package provided by Univ. of Tennessee, --
188* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
189*
190* .. Scalar Arguments ..
191 CHARACTER VECT
192 INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
193* ..
194* .. Array Arguments ..
195 DOUBLE PRECISION AB( LDAB, * ), C( LDC, * ), D( * ), E( * ),
196 $ PT( LDPT, * ), Q( LDQ, * ), WORK( * )
197* ..
198*
199* =====================================================================
200*
201* .. Parameters ..
202 DOUBLE PRECISION ZERO, ONE
203 parameter( zero = 0.0d+0, one = 1.0d+0 )
204* ..
205* .. Local Scalars ..
206 LOGICAL WANTB, WANTC, WANTPT, WANTQ
207 INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
208 $ KUN, L, MINMN, ML, ML0, MN, MU, MU0, NR, NRT
209 DOUBLE PRECISION RA, RB, RC, RS
210* ..
211* .. External Subroutines ..
212 EXTERNAL dlargv, dlartg, dlartv, dlaset, drot,
213 $ xerbla
214* ..
215* .. Intrinsic Functions ..
216 INTRINSIC max, min
217* ..
218* .. External Functions ..
219 LOGICAL LSAME
220 EXTERNAL lsame
221* ..
222* .. Executable Statements ..
223*
224* Test the input parameters
225*
226 wantb = lsame( vect, 'B' )
227 wantq = lsame( vect, 'Q' ) .OR. wantb
228 wantpt = lsame( vect, 'P' ) .OR. wantb
229 wantc = ncc.GT.0
230 klu1 = kl + ku + 1
231 info = 0
232 IF( .NOT.wantq .AND.
233 $ .NOT.wantpt .AND.
234 $ .NOT.lsame( vect, 'N' ) )
235 $ THEN
236 info = -1
237 ELSE IF( m.LT.0 ) THEN
238 info = -2
239 ELSE IF( n.LT.0 ) THEN
240 info = -3
241 ELSE IF( ncc.LT.0 ) THEN
242 info = -4
243 ELSE IF( kl.LT.0 ) THEN
244 info = -5
245 ELSE IF( ku.LT.0 ) THEN
246 info = -6
247 ELSE IF( ldab.LT.klu1 ) THEN
248 info = -8
249 ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
250 info = -12
251 ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
252 info = -14
253 ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
254 info = -16
255 END IF
256 IF( info.NE.0 ) THEN
257 CALL xerbla( 'DGBBRD', -info )
258 RETURN
259 END IF
260*
261* Initialize Q and P**T to the unit matrix, if needed
262*
263 IF( wantq )
264 $ CALL dlaset( 'Full', m, m, zero, one, q, ldq )
265 IF( wantpt )
266 $ CALL dlaset( 'Full', n, n, zero, one, pt, ldpt )
267*
268* Quick return if possible.
269*
270 IF( m.EQ.0 .OR. n.EQ.0 )
271 $ RETURN
272*
273 minmn = min( m, n )
274*
275 IF( kl+ku.GT.1 ) THEN
276*
277* Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
278* first to lower bidiagonal form and then transform to upper
279* bidiagonal
280*
281 IF( ku.GT.0 ) THEN
282 ml0 = 1
283 mu0 = 2
284 ELSE
285 ml0 = 2
286 mu0 = 1
287 END IF
288*
289* Wherever possible, plane rotations are generated and applied in
290* vector operations of length NR over the index set J1:J2:KLU1.
291*
292* The sines of the plane rotations are stored in WORK(1:max(m,n))
293* and the cosines in WORK(max(m,n)+1:2*max(m,n)).
294*
295 mn = max( m, n )
296 klm = min( m-1, kl )
297 kun = min( n-1, ku )
298 kb = klm + kun
299 kb1 = kb + 1
300 inca = kb1*ldab
301 nr = 0
302 j1 = klm + 2
303 j2 = 1 - kun
304*
305 DO 90 i = 1, minmn
306*
307* Reduce i-th column and i-th row of matrix to bidiagonal form
308*
309 ml = klm + 1
310 mu = kun + 1
311 DO 80 kk = 1, kb
312 j1 = j1 + kb
313 j2 = j2 + kb
314*
315* generate plane rotations to annihilate nonzero elements
316* which have been created below the band
317*
318 IF( nr.GT.0 )
319 $ CALL dlargv( nr, ab( klu1, j1-klm-1 ), inca,
320 $ work( j1 ), kb1, work( mn+j1 ), kb1 )
321*
322* apply plane rotations from the left
323*
324 DO 10 l = 1, kb
325 IF( j2-klm+l-1.GT.n ) THEN
326 nrt = nr - 1
327 ELSE
328 nrt = nr
329 END IF
330 IF( nrt.GT.0 )
331 $ CALL dlartv( nrt, ab( klu1-l, j1-klm+l-1 ),
332 $ inca,
333 $ ab( klu1-l+1, j1-klm+l-1 ), inca,
334 $ work( mn+j1 ), work( j1 ), kb1 )
335 10 CONTINUE
336*
337 IF( ml.GT.ml0 ) THEN
338 IF( ml.LE.m-i+1 ) THEN
339*
340* generate plane rotation to annihilate a(i+ml-1,i)
341* within the band, and apply rotation from the left
342*
343 CALL dlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
344 $ work( mn+i+ml-1 ), work( i+ml-1 ),
345 $ ra )
346 ab( ku+ml-1, i ) = ra
347 IF( i.LT.n )
348 $ CALL drot( min( ku+ml-2, n-i ),
349 $ ab( ku+ml-2, i+1 ), ldab-1,
350 $ ab( ku+ml-1, i+1 ), ldab-1,
351 $ work( mn+i+ml-1 ), work( i+ml-1 ) )
352 END IF
353 nr = nr + 1
354 j1 = j1 - kb1
355 END IF
356*
357 IF( wantq ) THEN
358*
359* accumulate product of plane rotations in Q
360*
361 DO 20 j = j1, j2, kb1
362 CALL drot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
363 $ work( mn+j ), work( j ) )
364 20 CONTINUE
365 END IF
366*
367 IF( wantc ) THEN
368*
369* apply plane rotations to C
370*
371 DO 30 j = j1, j2, kb1
372 CALL drot( ncc, c( j-1, 1 ), ldc, c( j, 1 ),
373 $ ldc,
374 $ work( mn+j ), work( j ) )
375 30 CONTINUE
376 END IF
377*
378 IF( j2+kun.GT.n ) THEN
379*
380* adjust J2 to keep within the bounds of the matrix
381*
382 nr = nr - 1
383 j2 = j2 - kb1
384 END IF
385*
386 DO 40 j = j1, j2, kb1
387*
388* create nonzero element a(j-1,j+ku) above the band
389* and store it in WORK(n+1:2*n)
390*
391 work( j+kun ) = work( j )*ab( 1, j+kun )
392 ab( 1, j+kun ) = work( mn+j )*ab( 1, j+kun )
393 40 CONTINUE
394*
395* generate plane rotations to annihilate nonzero elements
396* which have been generated above the band
397*
398 IF( nr.GT.0 )
399 $ CALL dlargv( nr, ab( 1, j1+kun-1 ), inca,
400 $ work( j1+kun ), kb1, work( mn+j1+kun ),
401 $ kb1 )
402*
403* apply plane rotations from the right
404*
405 DO 50 l = 1, kb
406 IF( j2+l-1.GT.m ) THEN
407 nrt = nr - 1
408 ELSE
409 nrt = nr
410 END IF
411 IF( nrt.GT.0 )
412 $ CALL dlartv( nrt, ab( l+1, j1+kun-1 ), inca,
413 $ ab( l, j1+kun ), inca,
414 $ work( mn+j1+kun ), work( j1+kun ),
415 $ kb1 )
416 50 CONTINUE
417*
418 IF( ml.EQ.ml0 .AND. mu.GT.mu0 ) THEN
419 IF( mu.LE.n-i+1 ) THEN
420*
421* generate plane rotation to annihilate a(i,i+mu-1)
422* within the band, and apply rotation from the right
423*
424 CALL dlartg( ab( ku-mu+3, i+mu-2 ),
425 $ ab( ku-mu+2, i+mu-1 ),
426 $ work( mn+i+mu-1 ), work( i+mu-1 ),
427 $ ra )
428 ab( ku-mu+3, i+mu-2 ) = ra
429 CALL drot( min( kl+mu-2, m-i ),
430 $ ab( ku-mu+4, i+mu-2 ), 1,
431 $ ab( ku-mu+3, i+mu-1 ), 1,
432 $ work( mn+i+mu-1 ), work( i+mu-1 ) )
433 END IF
434 nr = nr + 1
435 j1 = j1 - kb1
436 END IF
437*
438 IF( wantpt ) THEN
439*
440* accumulate product of plane rotations in P**T
441*
442 DO 60 j = j1, j2, kb1
443 CALL drot( n, pt( j+kun-1, 1 ), ldpt,
444 $ pt( j+kun, 1 ), ldpt, work( mn+j+kun ),
445 $ work( j+kun ) )
446 60 CONTINUE
447 END IF
448*
449 IF( j2+kb.GT.m ) THEN
450*
451* adjust J2 to keep within the bounds of the matrix
452*
453 nr = nr - 1
454 j2 = j2 - kb1
455 END IF
456*
457 DO 70 j = j1, j2, kb1
458*
459* create nonzero element a(j+kl+ku,j+ku-1) below the
460* band and store it in WORK(1:n)
461*
462 work( j+kb ) = work( j+kun )*ab( klu1, j+kun )
463 ab( klu1, j+kun ) = work( mn+j+kun )*ab( klu1, j+kun )
464 70 CONTINUE
465*
466 IF( ml.GT.ml0 ) THEN
467 ml = ml - 1
468 ELSE
469 mu = mu - 1
470 END IF
471 80 CONTINUE
472 90 CONTINUE
473 END IF
474*
475 IF( ku.EQ.0 .AND. kl.GT.0 ) THEN
476*
477* A has been reduced to lower bidiagonal form
478*
479* Transform lower bidiagonal form to upper bidiagonal by applying
480* plane rotations from the left, storing diagonal elements in D
481* and off-diagonal elements in E
482*
483 DO 100 i = 1, min( m-1, n )
484 CALL dlartg( ab( 1, i ), ab( 2, i ), rc, rs, ra )
485 d( i ) = ra
486 IF( i.LT.n ) THEN
487 e( i ) = rs*ab( 1, i+1 )
488 ab( 1, i+1 ) = rc*ab( 1, i+1 )
489 END IF
490 IF( wantq )
491 $ CALL drot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc, rs )
492 IF( wantc )
493 $ CALL drot( ncc, c( i, 1 ), ldc, c( i+1, 1 ), ldc, rc,
494 $ rs )
495 100 CONTINUE
496 IF( m.LE.n )
497 $ d( m ) = ab( 1, m )
498 ELSE IF( ku.GT.0 ) THEN
499*
500* A has been reduced to upper bidiagonal form
501*
502 IF( m.LT.n ) THEN
503*
504* Annihilate a(m,m+1) by applying plane rotations from the
505* right, storing diagonal elements in D and off-diagonal
506* elements in E
507*
508 rb = ab( ku, m+1 )
509 DO 110 i = m, 1, -1
510 CALL dlartg( ab( ku+1, i ), rb, rc, rs, ra )
511 d( i ) = ra
512 IF( i.GT.1 ) THEN
513 rb = -rs*ab( ku, i )
514 e( i-1 ) = rc*ab( ku, i )
515 END IF
516 IF( wantpt )
517 $ CALL drot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
518 $ rc, rs )
519 110 CONTINUE
520 ELSE
521*
522* Copy off-diagonal elements to E and diagonal elements to D
523*
524 DO 120 i = 1, minmn - 1
525 e( i ) = ab( ku, i+1 )
526 120 CONTINUE
527 DO 130 i = 1, minmn
528 d( i ) = ab( ku+1, i )
529 130 CONTINUE
530 END IF
531 ELSE
532*
533* A is diagonal. Set elements of E to zero and copy diagonal
534* elements to D.
535*
536 DO 140 i = 1, minmn - 1
537 e( i ) = zero
538 140 CONTINUE
539 DO 150 i = 1, minmn
540 d( i ) = ab( 1, i )
541 150 CONTINUE
542 END IF
543 RETURN
544*
545* End of DGBBRD
546*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dlargv
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:102
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlartv
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:106
dlaset
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:108
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
Here is the call graph for this function:
Here is the caller graph for this function: