LAPACK  3.9.1
LAPACK: Linear Algebra PACKage

◆ 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 191 of file cgbbrd.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  REAL D( * ), E( * ), RWORK( * )
204  COMPLEX AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
205  $ Q( LDQ, * ), WORK( * )
206 * ..
207 *
208 * =====================================================================
209 *
210 * .. Parameters ..
211  REAL ZERO
212  parameter( zero = 0.0e+0 )
213  COMPLEX CZERO, CONE
214  parameter( czero = ( 0.0e+0, 0.0e+0 ),
215  $ cone = ( 1.0e+0, 0.0e+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  REAL ABST, RC
222  COMPLEX RA, RB, RS, T
223 * ..
224 * .. External Subroutines ..
225  EXTERNAL clargv, clartg, clartv, claset, crot, cscal,
226  $ xerbla
227 * ..
228 * .. Intrinsic Functions ..
229  INTRINSIC abs, conjg, 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( 'CGBBRD', -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 claset( 'Full', m, m, czero, cone, q, ldq )
276  IF( wantpt )
277  $ CALL claset( '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 clargv( 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 clartv( 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 clartg( 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 crot( 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 crot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
371  $ rwork( j ), conjg( 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 crot( 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 clargv( 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 clartv( 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 clartg( 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 crot( 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 crot( n, pt( j+kun-1, 1 ), ldpt,
449  $ pt( j+kun, 1 ), ldpt, rwork( j+kun ),
450  $ conjg( 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 clartg( 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 crot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc,
497  $ conjg( rs ) )
498  IF( wantc )
499  $ CALL crot( 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 clartg( ab( ku+1, i ), rb, rc, rs, ra )
515  ab( ku+1, i ) = ra
516  IF( i.GT.1 ) THEN
517  rb = -conjg( rs )*ab( ku, i )
518  ab( ku, i ) = rc*ab( ku, i )
519  END IF
520  IF( wantpt )
521  $ CALL crot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
522  $ rc, conjg( 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 cscal( m, t, q( 1, i ), 1 )
541  IF( wantc )
542  $ CALL cscal( ncc, conjg( 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 )*conjg( t )
550  ELSE
551  t = ab( ku, i+1 )*conjg( 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 cscal( n, t, pt( i+1, 1 ), ldpt )
562  t = ab( ku+1, i+1 )*conjg( t )
563  END IF
564  END IF
565  120 CONTINUE
566  RETURN
567 *
568 * End of CGBBRD
569 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine clartg(F, G, CS, SN, R)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f:103
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:107
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:106
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:103
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:122
Here is the call graph for this function:
Here is the caller graph for this function: