LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
cgbbrd.f
1 *> \brief \b CGBBRD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGBBRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbbrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbbrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbbrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
22 * LDQ, PT, LDPT, C, LDC, WORK, RWORK, 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 * REAL D( * ), E( * ), RWORK( * )
30 * COMPLEX AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
31 * $ Q( LDQ, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> CGBBRD reduces a complex general m-by-n band matrix A to real upper
41 *> bidiagonal form B by a unitary transformation: Q**H * A * P = B.
42 *>
43 *> The routine computes B, and optionally forms Q or P**H, or computes
44 *> Q**H*C for a given matrix C.
45 *> \endverbatim
46 *
47 * Arguments:
48 * ==========
49 *
50 *> \param[in] VECT
51 *> \verbatim
52 *> VECT is CHARACTER*1
53 *> Specifies whether or not the matrices Q and P**H are to be
54 *> formed.
55 *> = 'N': do not form Q or P**H;
56 *> = 'Q': form Q only;
57 *> = 'P': form P**H only;
58 *> = 'B': form both.
59 *> \endverbatim
60 *>
61 *> \param[in] M
62 *> \verbatim
63 *> M is INTEGER
64 *> The number of rows of the matrix A. M >= 0.
65 *> \endverbatim
66 *>
67 *> \param[in] N
68 *> \verbatim
69 *> N is INTEGER
70 *> The number of columns of the matrix A. N >= 0.
71 *> \endverbatim
72 *>
73 *> \param[in] NCC
74 *> \verbatim
75 *> NCC is INTEGER
76 *> The number of columns of the matrix C. NCC >= 0.
77 *> \endverbatim
78 *>
79 *> \param[in] KL
80 *> \verbatim
81 *> KL is INTEGER
82 *> The number of subdiagonals of the matrix A. KL >= 0.
83 *> \endverbatim
84 *>
85 *> \param[in] KU
86 *> \verbatim
87 *> KU is INTEGER
88 *> The number of superdiagonals of the matrix A. KU >= 0.
89 *> \endverbatim
90 *>
91 *> \param[in,out] AB
92 *> \verbatim
93 *> AB is COMPLEX array, dimension (LDAB,N)
94 *> On entry, the m-by-n band matrix A, stored in rows 1 to
95 *> KL+KU+1. The j-th column of A is stored in the j-th column of
96 *> the array AB as follows:
97 *> AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
98 *> On exit, A is overwritten by values generated during the
99 *> reduction.
100 *> \endverbatim
101 *>
102 *> \param[in] LDAB
103 *> \verbatim
104 *> LDAB is INTEGER
105 *> The leading dimension of the array A. LDAB >= KL+KU+1.
106 *> \endverbatim
107 *>
108 *> \param[out] D
109 *> \verbatim
110 *> D is REAL array, dimension (min(M,N))
111 *> The diagonal elements of the bidiagonal matrix B.
112 *> \endverbatim
113 *>
114 *> \param[out] E
115 *> \verbatim
116 *> E is REAL array, dimension (min(M,N)-1)
117 *> The superdiagonal elements of the bidiagonal matrix B.
118 *> \endverbatim
119 *>
120 *> \param[out] Q
121 *> \verbatim
122 *> Q is COMPLEX array, dimension (LDQ,M)
123 *> If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
124 *> If VECT = 'N' or 'P', the array Q is not referenced.
125 *> \endverbatim
126 *>
127 *> \param[in] LDQ
128 *> \verbatim
129 *> LDQ is INTEGER
130 *> The leading dimension of the array Q.
131 *> LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
132 *> \endverbatim
133 *>
134 *> \param[out] PT
135 *> \verbatim
136 *> PT is COMPLEX array, dimension (LDPT,N)
137 *> If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
138 *> If VECT = 'N' or 'Q', the array PT is not referenced.
139 *> \endverbatim
140 *>
141 *> \param[in] LDPT
142 *> \verbatim
143 *> LDPT is INTEGER
144 *> The leading dimension of the array PT.
145 *> LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
146 *> \endverbatim
147 *>
148 *> \param[in,out] C
149 *> \verbatim
150 *> C is COMPLEX array, dimension (LDC,NCC)
151 *> On entry, an m-by-ncc matrix C.
152 *> On exit, C is overwritten by Q**H*C.
153 *> C is not referenced if NCC = 0.
154 *> \endverbatim
155 *>
156 *> \param[in] LDC
157 *> \verbatim
158 *> LDC is INTEGER
159 *> The leading dimension of the array C.
160 *> LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
161 *> \endverbatim
162 *>
163 *> \param[out] WORK
164 *> \verbatim
165 *> WORK is COMPLEX array, dimension (max(M,N))
166 *> \endverbatim
167 *>
168 *> \param[out] RWORK
169 *> \verbatim
170 *> RWORK is REAL array, dimension (max(M,N))
171 *> \endverbatim
172 *>
173 *> \param[out] INFO
174 *> \verbatim
175 *> INFO is INTEGER
176 *> = 0: successful exit.
177 *> < 0: if INFO = -i, the i-th argument had an illegal value.
178 *> \endverbatim
179 *
180 * Authors:
181 * ========
182 *
183 *> \author Univ. of Tennessee
184 *> \author Univ. of California Berkeley
185 *> \author Univ. of Colorado Denver
186 *> \author NAG Ltd.
187 *
188 *> \ingroup gbbrd
189 *
190 * =====================================================================
191  SUBROUTINE cgbbrd( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
192  $ LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
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 *
570  END
cgbbrd
subroutine cgbbrd(VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q, LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO)
CGBBRD
Definition: cgbbrd.f:193
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
cscal
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
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:122
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:106
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:103
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:107