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