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

◆ zgbbrd()

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

ZGBBRD

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

Purpose:
 ZGBBRD reduces a complex general m-by-n band matrix A to real upper
 bidiagonal form B by a unitary transformation: Q**H * A * P = B.

 The routine computes B, and optionally forms Q or P**H, or computes
 Q**H*C for a given matrix C.
Parameters
[in]VECT
          VECT is CHARACTER*1
          Specifies whether or not the matrices Q and P**H are to be
          formed.
          = 'N': do not form Q or P**H;
          = 'Q': form Q only;
          = 'P': form P**H 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDQ,M)
          If VECT = 'Q' or 'B', the m-by-m unitary 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 COMPLEX*16 array, dimension (LDPT,N)
          If VECT = 'P' or 'B', the n-by-n unitary 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 COMPLEX*16 array, dimension (LDC,NCC)
          On entry, an m-by-ncc matrix C.
          On exit, C is overwritten by Q**H*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 COMPLEX*16 array, dimension (max(M,N))
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (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 191 of file zgbbrd.f.

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