LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zgbbrd.f
Go to the documentation of this file.
1 *> \brief \b ZGBBRD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZGBBRD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgbbrd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgbbrd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgbbrd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGBBRD( 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 * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
30 * COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
31 * $ Q( LDQ, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> ZGBBRD 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 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*16 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*16 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*16 array, dimension (max(M,N))
166 *> \endverbatim
167 *>
168 *> \param[out] RWORK
169 *> \verbatim
170 *> RWORK is DOUBLE PRECISION 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 *> \date November 2011
189 *
190 *> \ingroup complex16GBcomputational
191 *
192 * =====================================================================
193  SUBROUTINE zgbbrd( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
194  $ ldq, pt, ldpt, c, ldc, work, rwork, info )
195 *
196 * -- LAPACK computational routine (version 3.4.0) --
197 * -- LAPACK is a software package provided by Univ. of Tennessee, --
198 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
199 * November 2011
200 *
201 * .. Scalar Arguments ..
202  CHARACTER vect
203  INTEGER info, kl, ku, ldab, ldc, ldpt, ldq, m, n, ncc
204 * ..
205 * .. Array Arguments ..
206  DOUBLE PRECISION d( * ), e( * ), rwork( * )
207  COMPLEX*16 ab( ldab, * ), c( ldc, * ), pt( ldpt, * ),
208  $ q( ldq, * ), work( * )
209 * ..
210 *
211 * =====================================================================
212 *
213 * .. Parameters ..
214  DOUBLE PRECISION zero
215  parameter( zero = 0.0d+0 )
216  COMPLEX*16 czero, cone
217  parameter( czero = ( 0.0d+0, 0.0d+0 ),
218  $ cone = ( 1.0d+0, 0.0d+0 ) )
219 * ..
220 * .. Local Scalars ..
221  LOGICAL wantb, wantc, wantpt, wantq
222  INTEGER i, inca, j, j1, j2, kb, kb1, kk, klm, klu1,
223  $ kun, l, minmn, ml, ml0, mu, mu0, nr, nrt
224  DOUBLE PRECISION abst, rc
225  COMPLEX*16 ra, rb, rs, t
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL xerbla, zlargv, zlartg, zlartv, zlaset, zrot,
229  $ zscal
230 * ..
231 * .. Intrinsic Functions ..
232  INTRINSIC abs, dconjg, max, min
233 * ..
234 * .. External Functions ..
235  LOGICAL lsame
236  EXTERNAL lsame
237 * ..
238 * .. Executable Statements ..
239 *
240 * Test the input parameters
241 *
242  wantb = lsame( vect, 'B' )
243  wantq = lsame( vect, 'Q' ) .OR. wantb
244  wantpt = lsame( vect, 'P' ) .OR. wantb
245  wantc = ncc.GT.0
246  klu1 = kl + ku + 1
247  info = 0
248  IF( .NOT.wantq .AND. .NOT.wantpt .AND. .NOT.lsame( vect, 'N' ) )
249  $ THEN
250  info = -1
251  ELSE IF( m.LT.0 ) THEN
252  info = -2
253  ELSE IF( n.LT.0 ) THEN
254  info = -3
255  ELSE IF( ncc.LT.0 ) THEN
256  info = -4
257  ELSE IF( kl.LT.0 ) THEN
258  info = -5
259  ELSE IF( ku.LT.0 ) THEN
260  info = -6
261  ELSE IF( ldab.LT.klu1 ) THEN
262  info = -8
263  ELSE IF( ldq.LT.1 .OR. wantq .AND. ldq.LT.max( 1, m ) ) THEN
264  info = -12
265  ELSE IF( ldpt.LT.1 .OR. wantpt .AND. ldpt.LT.max( 1, n ) ) THEN
266  info = -14
267  ELSE IF( ldc.LT.1 .OR. wantc .AND. ldc.LT.max( 1, m ) ) THEN
268  info = -16
269  END IF
270  IF( info.NE.0 ) THEN
271  CALL xerbla( 'ZGBBRD', -info )
272  return
273  END IF
274 *
275 * Initialize Q and P**H to the unit matrix, if needed
276 *
277  IF( wantq )
278  $ CALL zlaset( 'Full', m, m, czero, cone, q, ldq )
279  IF( wantpt )
280  $ CALL zlaset( 'Full', n, n, czero, cone, pt, ldpt )
281 *
282 * Quick return if possible.
283 *
284  IF( m.EQ.0 .OR. n.EQ.0 )
285  $ return
286 *
287  minmn = min( m, n )
288 *
289  IF( kl+ku.GT.1 ) THEN
290 *
291 * Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
292 * first to lower bidiagonal form and then transform to upper
293 * bidiagonal
294 *
295  IF( ku.GT.0 ) THEN
296  ml0 = 1
297  mu0 = 2
298  ELSE
299  ml0 = 2
300  mu0 = 1
301  END IF
302 *
303 * Wherever possible, plane rotations are generated and applied in
304 * vector operations of length NR over the index set J1:J2:KLU1.
305 *
306 * The complex sines of the plane rotations are stored in WORK,
307 * and the real cosines in RWORK.
308 *
309  klm = min( m-1, kl )
310  kun = min( n-1, ku )
311  kb = klm + kun
312  kb1 = kb + 1
313  inca = kb1*ldab
314  nr = 0
315  j1 = klm + 2
316  j2 = 1 - kun
317 *
318  DO 90 i = 1, minmn
319 *
320 * Reduce i-th column and i-th row of matrix to bidiagonal form
321 *
322  ml = klm + 1
323  mu = kun + 1
324  DO 80 kk = 1, kb
325  j1 = j1 + kb
326  j2 = j2 + kb
327 *
328 * generate plane rotations to annihilate nonzero elements
329 * which have been created below the band
330 *
331  IF( nr.GT.0 )
332  $ CALL zlargv( nr, ab( klu1, j1-klm-1 ), inca,
333  $ work( j1 ), kb1, rwork( j1 ), kb1 )
334 *
335 * apply plane rotations from the left
336 *
337  DO 10 l = 1, kb
338  IF( j2-klm+l-1.GT.n ) THEN
339  nrt = nr - 1
340  ELSE
341  nrt = nr
342  END IF
343  IF( nrt.GT.0 )
344  $ CALL zlartv( nrt, ab( klu1-l, j1-klm+l-1 ), inca,
345  $ ab( klu1-l+1, j1-klm+l-1 ), inca,
346  $ rwork( j1 ), work( j1 ), kb1 )
347  10 continue
348 *
349  IF( ml.GT.ml0 ) THEN
350  IF( ml.LE.m-i+1 ) THEN
351 *
352 * generate plane rotation to annihilate a(i+ml-1,i)
353 * within the band, and apply rotation from the left
354 *
355  CALL zlartg( ab( ku+ml-1, i ), ab( ku+ml, i ),
356  $ rwork( i+ml-1 ), work( i+ml-1 ), ra )
357  ab( ku+ml-1, i ) = ra
358  IF( i.LT.n )
359  $ CALL zrot( min( ku+ml-2, n-i ),
360  $ ab( ku+ml-2, i+1 ), ldab-1,
361  $ ab( ku+ml-1, i+1 ), ldab-1,
362  $ rwork( i+ml-1 ), work( i+ml-1 ) )
363  END IF
364  nr = nr + 1
365  j1 = j1 - kb1
366  END IF
367 *
368  IF( wantq ) THEN
369 *
370 * accumulate product of plane rotations in Q
371 *
372  DO 20 j = j1, j2, kb1
373  CALL zrot( m, q( 1, j-1 ), 1, q( 1, j ), 1,
374  $ rwork( j ), dconjg( work( j ) ) )
375  20 continue
376  END IF
377 *
378  IF( wantc ) THEN
379 *
380 * apply plane rotations to C
381 *
382  DO 30 j = j1, j2, kb1
383  CALL zrot( ncc, c( j-1, 1 ), ldc, c( j, 1 ), 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 zlargv( 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 zlartv( 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 zlartg( 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 zrot( 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 zrot( n, pt( j+kun-1, 1 ), ldpt,
452  $ pt( j+kun, 1 ), ldpt, rwork( j+kun ),
453  $ dconjg( 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 zlartg( 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 zrot( m, q( 1, i ), 1, q( 1, i+1 ), 1, rc,
500  $ dconjg( rs ) )
501  IF( wantc )
502  $ CALL zrot( 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 zlartg( ab( ku+1, i ), rb, rc, rs, ra )
518  ab( ku+1, i ) = ra
519  IF( i.GT.1 ) THEN
520  rb = -dconjg( rs )*ab( ku, i )
521  ab( ku, i ) = rc*ab( ku, i )
522  END IF
523  IF( wantpt )
524  $ CALL zrot( n, pt( i, 1 ), ldpt, pt( m+1, 1 ), ldpt,
525  $ rc, dconjg( 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 zscal( m, t, q( 1, i ), 1 )
544  IF( wantc )
545  $ CALL zscal( ncc, dconjg( 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 )*dconjg( t )
553  ELSE
554  t = ab( ku, i+1 )*dconjg( 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 zscal( n, t, pt( i+1, 1 ), ldpt )
565  t = ab( ku+1, i+1 )*dconjg( t )
566  END IF
567  END IF
568  120 continue
569  return
570 *
571 * End of ZGBBRD
572 *
573  END