Annotation of rpl/lapack/lapack/zgbbrd.f, revision 1.9
1.9 ! bertrand 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: * =====================================================================
1.1 bertrand 193: SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
194: $ LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
195: *
1.9 ! bertrand 196: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 197: * -- LAPACK is a software package provided by Univ. of Tennessee, --
198: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 199: * November 2011
1.1 bertrand 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: *
1.8 bertrand 275: * Initialize Q and P**H to the unit matrix, if needed
1.1 bertrand 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: *
1.8 bertrand 448: * accumulate product of plane rotations in P**H
1.1 bertrand 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>