Annotation of rpl/lapack/lapack/zgbbrd.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE ZGBBRD( VECT, M, N, NCC, KL, KU, AB, LDAB, D, E, Q,
! 2: $ LDQ, PT, LDPT, C, LDC, WORK, RWORK, INFO )
! 3: *
! 4: * -- LAPACK routine (version 3.2) --
! 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 7: * November 2006
! 8: *
! 9: * .. Scalar Arguments ..
! 10: CHARACTER VECT
! 11: INTEGER INFO, KL, KU, LDAB, LDC, LDPT, LDQ, M, N, NCC
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION D( * ), E( * ), RWORK( * )
! 15: COMPLEX*16 AB( LDAB, * ), C( LDC, * ), PT( LDPT, * ),
! 16: $ Q( LDQ, * ), WORK( * )
! 17: * ..
! 18: *
! 19: * Purpose
! 20: * =======
! 21: *
! 22: * ZGBBRD reduces a complex general m-by-n band matrix A to real upper
! 23: * bidiagonal form B by a unitary transformation: Q' * A * P = B.
! 24: *
! 25: * The routine computes B, and optionally forms Q or P', or computes
! 26: * Q'*C for a given matrix C.
! 27: *
! 28: * Arguments
! 29: * =========
! 30: *
! 31: * VECT (input) CHARACTER*1
! 32: * Specifies whether or not the matrices Q and P' are to be
! 33: * formed.
! 34: * = 'N': do not form Q or P';
! 35: * = 'Q': form Q only;
! 36: * = 'P': form P' only;
! 37: * = 'B': form both.
! 38: *
! 39: * M (input) INTEGER
! 40: * The number of rows of the matrix A. M >= 0.
! 41: *
! 42: * N (input) INTEGER
! 43: * The number of columns of the matrix A. N >= 0.
! 44: *
! 45: * NCC (input) INTEGER
! 46: * The number of columns of the matrix C. NCC >= 0.
! 47: *
! 48: * KL (input) INTEGER
! 49: * The number of subdiagonals of the matrix A. KL >= 0.
! 50: *
! 51: * KU (input) INTEGER
! 52: * The number of superdiagonals of the matrix A. KU >= 0.
! 53: *
! 54: * AB (input/output) COMPLEX*16 array, dimension (LDAB,N)
! 55: * On entry, the m-by-n band matrix A, stored in rows 1 to
! 56: * KL+KU+1. The j-th column of A is stored in the j-th column of
! 57: * the array AB as follows:
! 58: * AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
! 59: * On exit, A is overwritten by values generated during the
! 60: * reduction.
! 61: *
! 62: * LDAB (input) INTEGER
! 63: * The leading dimension of the array A. LDAB >= KL+KU+1.
! 64: *
! 65: * D (output) DOUBLE PRECISION array, dimension (min(M,N))
! 66: * The diagonal elements of the bidiagonal matrix B.
! 67: *
! 68: * E (output) DOUBLE PRECISION array, dimension (min(M,N)-1)
! 69: * The superdiagonal elements of the bidiagonal matrix B.
! 70: *
! 71: * Q (output) COMPLEX*16 array, dimension (LDQ,M)
! 72: * If VECT = 'Q' or 'B', the m-by-m unitary matrix Q.
! 73: * If VECT = 'N' or 'P', the array Q is not referenced.
! 74: *
! 75: * LDQ (input) INTEGER
! 76: * The leading dimension of the array Q.
! 77: * LDQ >= max(1,M) if VECT = 'Q' or 'B'; LDQ >= 1 otherwise.
! 78: *
! 79: * PT (output) COMPLEX*16 array, dimension (LDPT,N)
! 80: * If VECT = 'P' or 'B', the n-by-n unitary matrix P'.
! 81: * If VECT = 'N' or 'Q', the array PT is not referenced.
! 82: *
! 83: * LDPT (input) INTEGER
! 84: * The leading dimension of the array PT.
! 85: * LDPT >= max(1,N) if VECT = 'P' or 'B'; LDPT >= 1 otherwise.
! 86: *
! 87: * C (input/output) COMPLEX*16 array, dimension (LDC,NCC)
! 88: * On entry, an m-by-ncc matrix C.
! 89: * On exit, C is overwritten by Q'*C.
! 90: * C is not referenced if NCC = 0.
! 91: *
! 92: * LDC (input) INTEGER
! 93: * The leading dimension of the array C.
! 94: * LDC >= max(1,M) if NCC > 0; LDC >= 1 if NCC = 0.
! 95: *
! 96: * WORK (workspace) COMPLEX*16 array, dimension (max(M,N))
! 97: *
! 98: * RWORK (workspace) DOUBLE PRECISION array, dimension (max(M,N))
! 99: *
! 100: * INFO (output) INTEGER
! 101: * = 0: successful exit.
! 102: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 103: *
! 104: * =====================================================================
! 105: *
! 106: * .. Parameters ..
! 107: DOUBLE PRECISION ZERO
! 108: PARAMETER ( ZERO = 0.0D+0 )
! 109: COMPLEX*16 CZERO, CONE
! 110: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
! 111: $ CONE = ( 1.0D+0, 0.0D+0 ) )
! 112: * ..
! 113: * .. Local Scalars ..
! 114: LOGICAL WANTB, WANTC, WANTPT, WANTQ
! 115: INTEGER I, INCA, J, J1, J2, KB, KB1, KK, KLM, KLU1,
! 116: $ KUN, L, MINMN, ML, ML0, MU, MU0, NR, NRT
! 117: DOUBLE PRECISION ABST, RC
! 118: COMPLEX*16 RA, RB, RS, T
! 119: * ..
! 120: * .. External Subroutines ..
! 121: EXTERNAL XERBLA, ZLARGV, ZLARTG, ZLARTV, ZLASET, ZROT,
! 122: $ ZSCAL
! 123: * ..
! 124: * .. Intrinsic Functions ..
! 125: INTRINSIC ABS, DCONJG, MAX, MIN
! 126: * ..
! 127: * .. External Functions ..
! 128: LOGICAL LSAME
! 129: EXTERNAL LSAME
! 130: * ..
! 131: * .. Executable Statements ..
! 132: *
! 133: * Test the input parameters
! 134: *
! 135: WANTB = LSAME( VECT, 'B' )
! 136: WANTQ = LSAME( VECT, 'Q' ) .OR. WANTB
! 137: WANTPT = LSAME( VECT, 'P' ) .OR. WANTB
! 138: WANTC = NCC.GT.0
! 139: KLU1 = KL + KU + 1
! 140: INFO = 0
! 141: IF( .NOT.WANTQ .AND. .NOT.WANTPT .AND. .NOT.LSAME( VECT, 'N' ) )
! 142: $ THEN
! 143: INFO = -1
! 144: ELSE IF( M.LT.0 ) THEN
! 145: INFO = -2
! 146: ELSE IF( N.LT.0 ) THEN
! 147: INFO = -3
! 148: ELSE IF( NCC.LT.0 ) THEN
! 149: INFO = -4
! 150: ELSE IF( KL.LT.0 ) THEN
! 151: INFO = -5
! 152: ELSE IF( KU.LT.0 ) THEN
! 153: INFO = -6
! 154: ELSE IF( LDAB.LT.KLU1 ) THEN
! 155: INFO = -8
! 156: ELSE IF( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.MAX( 1, M ) ) THEN
! 157: INFO = -12
! 158: ELSE IF( LDPT.LT.1 .OR. WANTPT .AND. LDPT.LT.MAX( 1, N ) ) THEN
! 159: INFO = -14
! 160: ELSE IF( LDC.LT.1 .OR. WANTC .AND. LDC.LT.MAX( 1, M ) ) THEN
! 161: INFO = -16
! 162: END IF
! 163: IF( INFO.NE.0 ) THEN
! 164: CALL XERBLA( 'ZGBBRD', -INFO )
! 165: RETURN
! 166: END IF
! 167: *
! 168: * Initialize Q and P' to the unit matrix, if needed
! 169: *
! 170: IF( WANTQ )
! 171: $ CALL ZLASET( 'Full', M, M, CZERO, CONE, Q, LDQ )
! 172: IF( WANTPT )
! 173: $ CALL ZLASET( 'Full', N, N, CZERO, CONE, PT, LDPT )
! 174: *
! 175: * Quick return if possible.
! 176: *
! 177: IF( M.EQ.0 .OR. N.EQ.0 )
! 178: $ RETURN
! 179: *
! 180: MINMN = MIN( M, N )
! 181: *
! 182: IF( KL+KU.GT.1 ) THEN
! 183: *
! 184: * Reduce to upper bidiagonal form if KU > 0; if KU = 0, reduce
! 185: * first to lower bidiagonal form and then transform to upper
! 186: * bidiagonal
! 187: *
! 188: IF( KU.GT.0 ) THEN
! 189: ML0 = 1
! 190: MU0 = 2
! 191: ELSE
! 192: ML0 = 2
! 193: MU0 = 1
! 194: END IF
! 195: *
! 196: * Wherever possible, plane rotations are generated and applied in
! 197: * vector operations of length NR over the index set J1:J2:KLU1.
! 198: *
! 199: * The complex sines of the plane rotations are stored in WORK,
! 200: * and the real cosines in RWORK.
! 201: *
! 202: KLM = MIN( M-1, KL )
! 203: KUN = MIN( N-1, KU )
! 204: KB = KLM + KUN
! 205: KB1 = KB + 1
! 206: INCA = KB1*LDAB
! 207: NR = 0
! 208: J1 = KLM + 2
! 209: J2 = 1 - KUN
! 210: *
! 211: DO 90 I = 1, MINMN
! 212: *
! 213: * Reduce i-th column and i-th row of matrix to bidiagonal form
! 214: *
! 215: ML = KLM + 1
! 216: MU = KUN + 1
! 217: DO 80 KK = 1, KB
! 218: J1 = J1 + KB
! 219: J2 = J2 + KB
! 220: *
! 221: * generate plane rotations to annihilate nonzero elements
! 222: * which have been created below the band
! 223: *
! 224: IF( NR.GT.0 )
! 225: $ CALL ZLARGV( NR, AB( KLU1, J1-KLM-1 ), INCA,
! 226: $ WORK( J1 ), KB1, RWORK( J1 ), KB1 )
! 227: *
! 228: * apply plane rotations from the left
! 229: *
! 230: DO 10 L = 1, KB
! 231: IF( J2-KLM+L-1.GT.N ) THEN
! 232: NRT = NR - 1
! 233: ELSE
! 234: NRT = NR
! 235: END IF
! 236: IF( NRT.GT.0 )
! 237: $ CALL ZLARTV( NRT, AB( KLU1-L, J1-KLM+L-1 ), INCA,
! 238: $ AB( KLU1-L+1, J1-KLM+L-1 ), INCA,
! 239: $ RWORK( J1 ), WORK( J1 ), KB1 )
! 240: 10 CONTINUE
! 241: *
! 242: IF( ML.GT.ML0 ) THEN
! 243: IF( ML.LE.M-I+1 ) THEN
! 244: *
! 245: * generate plane rotation to annihilate a(i+ml-1,i)
! 246: * within the band, and apply rotation from the left
! 247: *
! 248: CALL ZLARTG( AB( KU+ML-1, I ), AB( KU+ML, I ),
! 249: $ RWORK( I+ML-1 ), WORK( I+ML-1 ), RA )
! 250: AB( KU+ML-1, I ) = RA
! 251: IF( I.LT.N )
! 252: $ CALL ZROT( MIN( KU+ML-2, N-I ),
! 253: $ AB( KU+ML-2, I+1 ), LDAB-1,
! 254: $ AB( KU+ML-1, I+1 ), LDAB-1,
! 255: $ RWORK( I+ML-1 ), WORK( I+ML-1 ) )
! 256: END IF
! 257: NR = NR + 1
! 258: J1 = J1 - KB1
! 259: END IF
! 260: *
! 261: IF( WANTQ ) THEN
! 262: *
! 263: * accumulate product of plane rotations in Q
! 264: *
! 265: DO 20 J = J1, J2, KB1
! 266: CALL ZROT( M, Q( 1, J-1 ), 1, Q( 1, J ), 1,
! 267: $ RWORK( J ), DCONJG( WORK( J ) ) )
! 268: 20 CONTINUE
! 269: END IF
! 270: *
! 271: IF( WANTC ) THEN
! 272: *
! 273: * apply plane rotations to C
! 274: *
! 275: DO 30 J = J1, J2, KB1
! 276: CALL ZROT( NCC, C( J-1, 1 ), LDC, C( J, 1 ), LDC,
! 277: $ RWORK( J ), WORK( J ) )
! 278: 30 CONTINUE
! 279: END IF
! 280: *
! 281: IF( J2+KUN.GT.N ) THEN
! 282: *
! 283: * adjust J2 to keep within the bounds of the matrix
! 284: *
! 285: NR = NR - 1
! 286: J2 = J2 - KB1
! 287: END IF
! 288: *
! 289: DO 40 J = J1, J2, KB1
! 290: *
! 291: * create nonzero element a(j-1,j+ku) above the band
! 292: * and store it in WORK(n+1:2*n)
! 293: *
! 294: WORK( J+KUN ) = WORK( J )*AB( 1, J+KUN )
! 295: AB( 1, J+KUN ) = RWORK( J )*AB( 1, J+KUN )
! 296: 40 CONTINUE
! 297: *
! 298: * generate plane rotations to annihilate nonzero elements
! 299: * which have been generated above the band
! 300: *
! 301: IF( NR.GT.0 )
! 302: $ CALL ZLARGV( NR, AB( 1, J1+KUN-1 ), INCA,
! 303: $ WORK( J1+KUN ), KB1, RWORK( J1+KUN ),
! 304: $ KB1 )
! 305: *
! 306: * apply plane rotations from the right
! 307: *
! 308: DO 50 L = 1, KB
! 309: IF( J2+L-1.GT.M ) THEN
! 310: NRT = NR - 1
! 311: ELSE
! 312: NRT = NR
! 313: END IF
! 314: IF( NRT.GT.0 )
! 315: $ CALL ZLARTV( NRT, AB( L+1, J1+KUN-1 ), INCA,
! 316: $ AB( L, J1+KUN ), INCA,
! 317: $ RWORK( J1+KUN ), WORK( J1+KUN ), KB1 )
! 318: 50 CONTINUE
! 319: *
! 320: IF( ML.EQ.ML0 .AND. MU.GT.MU0 ) THEN
! 321: IF( MU.LE.N-I+1 ) THEN
! 322: *
! 323: * generate plane rotation to annihilate a(i,i+mu-1)
! 324: * within the band, and apply rotation from the right
! 325: *
! 326: CALL ZLARTG( AB( KU-MU+3, I+MU-2 ),
! 327: $ AB( KU-MU+2, I+MU-1 ),
! 328: $ RWORK( I+MU-1 ), WORK( I+MU-1 ), RA )
! 329: AB( KU-MU+3, I+MU-2 ) = RA
! 330: CALL ZROT( MIN( KL+MU-2, M-I ),
! 331: $ AB( KU-MU+4, I+MU-2 ), 1,
! 332: $ AB( KU-MU+3, I+MU-1 ), 1,
! 333: $ RWORK( I+MU-1 ), WORK( I+MU-1 ) )
! 334: END IF
! 335: NR = NR + 1
! 336: J1 = J1 - KB1
! 337: END IF
! 338: *
! 339: IF( WANTPT ) THEN
! 340: *
! 341: * accumulate product of plane rotations in P'
! 342: *
! 343: DO 60 J = J1, J2, KB1
! 344: CALL ZROT( N, PT( J+KUN-1, 1 ), LDPT,
! 345: $ PT( J+KUN, 1 ), LDPT, RWORK( J+KUN ),
! 346: $ DCONJG( WORK( J+KUN ) ) )
! 347: 60 CONTINUE
! 348: END IF
! 349: *
! 350: IF( J2+KB.GT.M ) THEN
! 351: *
! 352: * adjust J2 to keep within the bounds of the matrix
! 353: *
! 354: NR = NR - 1
! 355: J2 = J2 - KB1
! 356: END IF
! 357: *
! 358: DO 70 J = J1, J2, KB1
! 359: *
! 360: * create nonzero element a(j+kl+ku,j+ku-1) below the
! 361: * band and store it in WORK(1:n)
! 362: *
! 363: WORK( J+KB ) = WORK( J+KUN )*AB( KLU1, J+KUN )
! 364: AB( KLU1, J+KUN ) = RWORK( J+KUN )*AB( KLU1, J+KUN )
! 365: 70 CONTINUE
! 366: *
! 367: IF( ML.GT.ML0 ) THEN
! 368: ML = ML - 1
! 369: ELSE
! 370: MU = MU - 1
! 371: END IF
! 372: 80 CONTINUE
! 373: 90 CONTINUE
! 374: END IF
! 375: *
! 376: IF( KU.EQ.0 .AND. KL.GT.0 ) THEN
! 377: *
! 378: * A has been reduced to complex lower bidiagonal form
! 379: *
! 380: * Transform lower bidiagonal form to upper bidiagonal by applying
! 381: * plane rotations from the left, overwriting superdiagonal
! 382: * elements on subdiagonal elements
! 383: *
! 384: DO 100 I = 1, MIN( M-1, N )
! 385: CALL ZLARTG( AB( 1, I ), AB( 2, I ), RC, RS, RA )
! 386: AB( 1, I ) = RA
! 387: IF( I.LT.N ) THEN
! 388: AB( 2, I ) = RS*AB( 1, I+1 )
! 389: AB( 1, I+1 ) = RC*AB( 1, I+1 )
! 390: END IF
! 391: IF( WANTQ )
! 392: $ CALL ZROT( M, Q( 1, I ), 1, Q( 1, I+1 ), 1, RC,
! 393: $ DCONJG( RS ) )
! 394: IF( WANTC )
! 395: $ CALL ZROT( NCC, C( I, 1 ), LDC, C( I+1, 1 ), LDC, RC,
! 396: $ RS )
! 397: 100 CONTINUE
! 398: ELSE
! 399: *
! 400: * A has been reduced to complex upper bidiagonal form or is
! 401: * diagonal
! 402: *
! 403: IF( KU.GT.0 .AND. M.LT.N ) THEN
! 404: *
! 405: * Annihilate a(m,m+1) by applying plane rotations from the
! 406: * right
! 407: *
! 408: RB = AB( KU, M+1 )
! 409: DO 110 I = M, 1, -1
! 410: CALL ZLARTG( AB( KU+1, I ), RB, RC, RS, RA )
! 411: AB( KU+1, I ) = RA
! 412: IF( I.GT.1 ) THEN
! 413: RB = -DCONJG( RS )*AB( KU, I )
! 414: AB( KU, I ) = RC*AB( KU, I )
! 415: END IF
! 416: IF( WANTPT )
! 417: $ CALL ZROT( N, PT( I, 1 ), LDPT, PT( M+1, 1 ), LDPT,
! 418: $ RC, DCONJG( RS ) )
! 419: 110 CONTINUE
! 420: END IF
! 421: END IF
! 422: *
! 423: * Make diagonal and superdiagonal elements real, storing them in D
! 424: * and E
! 425: *
! 426: T = AB( KU+1, 1 )
! 427: DO 120 I = 1, MINMN
! 428: ABST = ABS( T )
! 429: D( I ) = ABST
! 430: IF( ABST.NE.ZERO ) THEN
! 431: T = T / ABST
! 432: ELSE
! 433: T = CONE
! 434: END IF
! 435: IF( WANTQ )
! 436: $ CALL ZSCAL( M, T, Q( 1, I ), 1 )
! 437: IF( WANTC )
! 438: $ CALL ZSCAL( NCC, DCONJG( T ), C( I, 1 ), LDC )
! 439: IF( I.LT.MINMN ) THEN
! 440: IF( KU.EQ.0 .AND. KL.EQ.0 ) THEN
! 441: E( I ) = ZERO
! 442: T = AB( 1, I+1 )
! 443: ELSE
! 444: IF( KU.EQ.0 ) THEN
! 445: T = AB( 2, I )*DCONJG( T )
! 446: ELSE
! 447: T = AB( KU, I+1 )*DCONJG( T )
! 448: END IF
! 449: ABST = ABS( T )
! 450: E( I ) = ABST
! 451: IF( ABST.NE.ZERO ) THEN
! 452: T = T / ABST
! 453: ELSE
! 454: T = CONE
! 455: END IF
! 456: IF( WANTPT )
! 457: $ CALL ZSCAL( N, T, PT( I+1, 1 ), LDPT )
! 458: T = AB( KU+1, I+1 )*DCONJG( T )
! 459: END IF
! 460: END IF
! 461: 120 CONTINUE
! 462: RETURN
! 463: *
! 464: * End of ZGBBRD
! 465: *
! 466: 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