![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
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