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: *> \ingroup complex16GBcomputational
189: *
190: * =====================================================================
191: SUBROUTINE ZGBBRD( 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: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
204: COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
205: $ Q( LDQ, * ), WORK( * )
206: * ..
207: *
208: * =====================================================================
209: *
210: * .. Parameters ..
211: DOUBLE PRECISION ZERO
212: PARAMETER ( ZERO = 0.0D+0 )
213: COMPLEX*16 CZERO, CONE
214: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
215: $ CONE = ( 1.0D+0, 0.0D+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: DOUBLE PRECISION ABST, RC
222: COMPLEX*16 RA, RB, RS, T
223: * ..
224: * .. External Subroutines ..
225: EXTERNAL XERBLA, ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT,
226: $ ZSCAL
227: * ..
228: * .. Intrinsic Functions ..
229: INTRINSIC ABS, DCONJG, 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( 'ZGBBRD', -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 ZLASET( 'Full', M, M, CZERO, CONE, Q, LDQ )
276: IF( WANTPT )
277: $ CALL ZLASET( '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 ZLARGV( 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 ZLARTV( 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 ZLARTG( 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 ZROT( 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 ZROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
371: $ RWORK( J ), DCONJG( 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 ZROT( 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 ZLARGV( 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 ZLARTV( 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 ZLARTG( 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 ZROT( 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 ZROT( N, PT( J+KUN-1, 1 ), LDPT,
449: $ PT( J+KUN, 1 ), LDPT, RWORK( J+KUN ),
450: $ DCONJG( 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 ZLARTG( 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 ZROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC,
497: $ DCONJG( RS ) )
498: IF( WANTC )
499: $ CALL ZROT( 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 ZLARTG( AB( KU+1, I ), RB, RC, RS, RA )
515: AB( KU+1, I ) = RA
516: IF( I.GT.1 ) THEN
517: RB = -DCONJG( RS )*AB( KU, I )
518: AB( KU, I ) = RC*AB( KU, I )
519: END IF
520: IF( WANTPT )
521: $ CALL ZROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
522: $ RC, DCONJG( 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 ZSCAL( M, T, Q( 1, I ), 1 )
541: IF( WANTC )
542: $ CALL ZSCAL( NCC, DCONJG( 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 )*DCONJG( T )
550: ELSE
551: T = AB( KU, I+1 )*DCONJG( 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 ZSCAL( N, T, PT( I+1, 1 ), LDPT )
562: T = AB( KU+1, I+1 )*DCONJG( T )
563: END IF
564: END IF
565: 120 CONTINUE
566: RETURN
567: *
568: * End of ZGBBRD
569: *
570: 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