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

◆ cgbbrd()

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

CGBBRD

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

Purpose:
!>
!> CGBBRD 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 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 REAL array, dimension (min(M,N))
!>          The diagonal elements of the bidiagonal matrix B.
!>
[out]E
!>          E is REAL array, dimension (min(M,N)-1)
!>          The superdiagonal elements of the bidiagonal matrix B.
!>
[out]Q
!>          Q is COMPLEX 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 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 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 array, dimension (max(M,N))
!>
[out]RWORK
!>          RWORK is REAL 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 189 of file cgbbrd.f.

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