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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgbbrd.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack