Annotation of rpl/lapack/lapack/dsbgst.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DSBGST( VECT, UPLO, N, KA, KB, AB, LDAB, BB, LDBB, X,
! 2: $ LDX, 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 UPLO, VECT
! 11: INTEGER INFO, KA, KB, LDAB, LDBB, LDX, N
! 12: * ..
! 13: * .. Array Arguments ..
! 14: DOUBLE PRECISION AB( LDAB, * ), BB( LDBB, * ), WORK( * ),
! 15: $ X( LDX, * )
! 16: * ..
! 17: *
! 18: * Purpose
! 19: * =======
! 20: *
! 21: * DSBGST reduces a real symmetric-definite banded generalized
! 22: * eigenproblem A*x = lambda*B*x to standard form C*y = lambda*y,
! 23: * such that C has the same bandwidth as A.
! 24: *
! 25: * B must have been previously factorized as S**T*S by DPBSTF, using a
! 26: * split Cholesky factorization. A is overwritten by C = X**T*A*X, where
! 27: * X = S**(-1)*Q and Q is an orthogonal matrix chosen to preserve the
! 28: * bandwidth of A.
! 29: *
! 30: * Arguments
! 31: * =========
! 32: *
! 33: * VECT (input) CHARACTER*1
! 34: * = 'N': do not form the transformation matrix X;
! 35: * = 'V': form X.
! 36: *
! 37: * UPLO (input) CHARACTER*1
! 38: * = 'U': Upper triangle of A is stored;
! 39: * = 'L': Lower triangle of A is stored.
! 40: *
! 41: * N (input) INTEGER
! 42: * The order of the matrices A and B. N >= 0.
! 43: *
! 44: * KA (input) INTEGER
! 45: * The number of superdiagonals of the matrix A if UPLO = 'U',
! 46: * or the number of subdiagonals if UPLO = 'L'. KA >= 0.
! 47: *
! 48: * KB (input) INTEGER
! 49: * The number of superdiagonals of the matrix B if UPLO = 'U',
! 50: * or the number of subdiagonals if UPLO = 'L'. KA >= KB >= 0.
! 51: *
! 52: * AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
! 53: * On entry, the upper or lower triangle of the symmetric band
! 54: * matrix A, stored in the first ka+1 rows of the array. The
! 55: * j-th column of A is stored in the j-th column of the array AB
! 56: * as follows:
! 57: * if UPLO = 'U', AB(ka+1+i-j,j) = A(i,j) for max(1,j-ka)<=i<=j;
! 58: * if UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+ka).
! 59: *
! 60: * On exit, the transformed matrix X**T*A*X, stored in the same
! 61: * format as A.
! 62: *
! 63: * LDAB (input) INTEGER
! 64: * The leading dimension of the array AB. LDAB >= KA+1.
! 65: *
! 66: * BB (input) DOUBLE PRECISION array, dimension (LDBB,N)
! 67: * The banded factor S from the split Cholesky factorization of
! 68: * B, as returned by DPBSTF, stored in the first KB+1 rows of
! 69: * the array.
! 70: *
! 71: * LDBB (input) INTEGER
! 72: * The leading dimension of the array BB. LDBB >= KB+1.
! 73: *
! 74: * X (output) DOUBLE PRECISION array, dimension (LDX,N)
! 75: * If VECT = 'V', the n-by-n matrix X.
! 76: * If VECT = 'N', the array X is not referenced.
! 77: *
! 78: * LDX (input) INTEGER
! 79: * The leading dimension of the array X.
! 80: * LDX >= max(1,N) if VECT = 'V'; LDX >= 1 otherwise.
! 81: *
! 82: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
! 83: *
! 84: * INFO (output) INTEGER
! 85: * = 0: successful exit
! 86: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 87: *
! 88: * =====================================================================
! 89: *
! 90: * .. Parameters ..
! 91: DOUBLE PRECISION ZERO, ONE
! 92: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 93: * ..
! 94: * .. Local Scalars ..
! 95: LOGICAL UPDATE, UPPER, WANTX
! 96: INTEGER I, I0, I1, I2, INCA, J, J1, J1T, J2, J2T, K,
! 97: $ KA1, KB1, KBT, L, M, NR, NRT, NX
! 98: DOUBLE PRECISION BII, RA, RA1, T
! 99: * ..
! 100: * .. External Functions ..
! 101: LOGICAL LSAME
! 102: EXTERNAL LSAME
! 103: * ..
! 104: * .. External Subroutines ..
! 105: EXTERNAL DGER, DLAR2V, DLARGV, DLARTG, DLARTV, DLASET,
! 106: $ DROT, DSCAL, XERBLA
! 107: * ..
! 108: * .. Intrinsic Functions ..
! 109: INTRINSIC MAX, MIN
! 110: * ..
! 111: * .. Executable Statements ..
! 112: *
! 113: * Test the input parameters
! 114: *
! 115: WANTX = LSAME( VECT, 'V' )
! 116: UPPER = LSAME( UPLO, 'U' )
! 117: KA1 = KA + 1
! 118: KB1 = KB + 1
! 119: INFO = 0
! 120: IF( .NOT.WANTX .AND. .NOT.LSAME( VECT, 'N' ) ) THEN
! 121: INFO = -1
! 122: ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
! 123: INFO = -2
! 124: ELSE IF( N.LT.0 ) THEN
! 125: INFO = -3
! 126: ELSE IF( KA.LT.0 ) THEN
! 127: INFO = -4
! 128: ELSE IF( KB.LT.0 .OR. KB.GT.KA ) THEN
! 129: INFO = -5
! 130: ELSE IF( LDAB.LT.KA+1 ) THEN
! 131: INFO = -7
! 132: ELSE IF( LDBB.LT.KB+1 ) THEN
! 133: INFO = -9
! 134: ELSE IF( LDX.LT.1 .OR. WANTX .AND. LDX.LT.MAX( 1, N ) ) THEN
! 135: INFO = -11
! 136: END IF
! 137: IF( INFO.NE.0 ) THEN
! 138: CALL XERBLA( 'DSBGST', -INFO )
! 139: RETURN
! 140: END IF
! 141: *
! 142: * Quick return if possible
! 143: *
! 144: IF( N.EQ.0 )
! 145: $ RETURN
! 146: *
! 147: INCA = LDAB*KA1
! 148: *
! 149: * Initialize X to the unit matrix, if needed
! 150: *
! 151: IF( WANTX )
! 152: $ CALL DLASET( 'Full', N, N, ZERO, ONE, X, LDX )
! 153: *
! 154: * Set M to the splitting point m. It must be the same value as is
! 155: * used in DPBSTF. The chosen value allows the arrays WORK and RWORK
! 156: * to be of dimension (N).
! 157: *
! 158: M = ( N+KB ) / 2
! 159: *
! 160: * The routine works in two phases, corresponding to the two halves
! 161: * of the split Cholesky factorization of B as S**T*S where
! 162: *
! 163: * S = ( U )
! 164: * ( M L )
! 165: *
! 166: * with U upper triangular of order m, and L lower triangular of
! 167: * order n-m. S has the same bandwidth as B.
! 168: *
! 169: * S is treated as a product of elementary matrices:
! 170: *
! 171: * S = S(m)*S(m-1)*...*S(2)*S(1)*S(m+1)*S(m+2)*...*S(n-1)*S(n)
! 172: *
! 173: * where S(i) is determined by the i-th row of S.
! 174: *
! 175: * In phase 1, the index i takes the values n, n-1, ... , m+1;
! 176: * in phase 2, it takes the values 1, 2, ... , m.
! 177: *
! 178: * For each value of i, the current matrix A is updated by forming
! 179: * inv(S(i))**T*A*inv(S(i)). This creates a triangular bulge outside
! 180: * the band of A. The bulge is then pushed down toward the bottom of
! 181: * A in phase 1, and up toward the top of A in phase 2, by applying
! 182: * plane rotations.
! 183: *
! 184: * There are kb*(kb+1)/2 elements in the bulge, but at most 2*kb-1
! 185: * of them are linearly independent, so annihilating a bulge requires
! 186: * only 2*kb-1 plane rotations. The rotations are divided into a 1st
! 187: * set of kb-1 rotations, and a 2nd set of kb rotations.
! 188: *
! 189: * Wherever possible, rotations are generated and applied in vector
! 190: * operations of length NR between the indices J1 and J2 (sometimes
! 191: * replaced by modified values NRT, J1T or J2T).
! 192: *
! 193: * The cosines and sines of the rotations are stored in the array
! 194: * WORK. The cosines of the 1st set of rotations are stored in
! 195: * elements n+2:n+m-kb-1 and the sines of the 1st set in elements
! 196: * 2:m-kb-1; the cosines of the 2nd set are stored in elements
! 197: * n+m-kb+1:2*n and the sines of the second set in elements m-kb+1:n.
! 198: *
! 199: * The bulges are not formed explicitly; nonzero elements outside the
! 200: * band are created only when they are required for generating new
! 201: * rotations; they are stored in the array WORK, in positions where
! 202: * they are later overwritten by the sines of the rotations which
! 203: * annihilate them.
! 204: *
! 205: * **************************** Phase 1 *****************************
! 206: *
! 207: * The logical structure of this phase is:
! 208: *
! 209: * UPDATE = .TRUE.
! 210: * DO I = N, M + 1, -1
! 211: * use S(i) to update A and create a new bulge
! 212: * apply rotations to push all bulges KA positions downward
! 213: * END DO
! 214: * UPDATE = .FALSE.
! 215: * DO I = M + KA + 1, N - 1
! 216: * apply rotations to push all bulges KA positions downward
! 217: * END DO
! 218: *
! 219: * To avoid duplicating code, the two loops are merged.
! 220: *
! 221: UPDATE = .TRUE.
! 222: I = N + 1
! 223: 10 CONTINUE
! 224: IF( UPDATE ) THEN
! 225: I = I - 1
! 226: KBT = MIN( KB, I-1 )
! 227: I0 = I - 1
! 228: I1 = MIN( N, I+KA )
! 229: I2 = I - KBT + KA1
! 230: IF( I.LT.M+1 ) THEN
! 231: UPDATE = .FALSE.
! 232: I = I + 1
! 233: I0 = M
! 234: IF( KA.EQ.0 )
! 235: $ GO TO 480
! 236: GO TO 10
! 237: END IF
! 238: ELSE
! 239: I = I + KA
! 240: IF( I.GT.N-1 )
! 241: $ GO TO 480
! 242: END IF
! 243: *
! 244: IF( UPPER ) THEN
! 245: *
! 246: * Transform A, working with the upper triangle
! 247: *
! 248: IF( UPDATE ) THEN
! 249: *
! 250: * Form inv(S(i))**T * A * inv(S(i))
! 251: *
! 252: BII = BB( KB1, I )
! 253: DO 20 J = I, I1
! 254: AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
! 255: 20 CONTINUE
! 256: DO 30 J = MAX( 1, I-KA ), I
! 257: AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
! 258: 30 CONTINUE
! 259: DO 60 K = I - KBT, I - 1
! 260: DO 40 J = I - KBT, K
! 261: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
! 262: $ BB( J-I+KB1, I )*AB( K-I+KA1, I ) -
! 263: $ BB( K-I+KB1, I )*AB( J-I+KA1, I ) +
! 264: $ AB( KA1, I )*BB( J-I+KB1, I )*
! 265: $ BB( K-I+KB1, I )
! 266: 40 CONTINUE
! 267: DO 50 J = MAX( 1, I-KA ), I - KBT - 1
! 268: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
! 269: $ BB( K-I+KB1, I )*AB( J-I+KA1, I )
! 270: 50 CONTINUE
! 271: 60 CONTINUE
! 272: DO 80 J = I, I1
! 273: DO 70 K = MAX( J-KA, I-KBT ), I - 1
! 274: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
! 275: $ BB( K-I+KB1, I )*AB( I-J+KA1, J )
! 276: 70 CONTINUE
! 277: 80 CONTINUE
! 278: *
! 279: IF( WANTX ) THEN
! 280: *
! 281: * post-multiply X by inv(S(i))
! 282: *
! 283: CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
! 284: IF( KBT.GT.0 )
! 285: $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
! 286: $ BB( KB1-KBT, I ), 1, X( M+1, I-KBT ), LDX )
! 287: END IF
! 288: *
! 289: * store a(i,i1) in RA1 for use in next loop over K
! 290: *
! 291: RA1 = AB( I-I1+KA1, I1 )
! 292: END IF
! 293: *
! 294: * Generate and apply vectors of rotations to chase all the
! 295: * existing bulges KA positions down toward the bottom of the
! 296: * band
! 297: *
! 298: DO 130 K = 1, KB - 1
! 299: IF( UPDATE ) THEN
! 300: *
! 301: * Determine the rotations which would annihilate the bulge
! 302: * which has in theory just been created
! 303: *
! 304: IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
! 305: *
! 306: * generate rotation to annihilate a(i,i-k+ka+1)
! 307: *
! 308: CALL DLARTG( AB( K+1, I-K+KA ), RA1,
! 309: $ WORK( N+I-K+KA-M ), WORK( I-K+KA-M ),
! 310: $ RA )
! 311: *
! 312: * create nonzero element a(i-k,i-k+ka+1) outside the
! 313: * band and store it in WORK(i-k)
! 314: *
! 315: T = -BB( KB1-K, I )*RA1
! 316: WORK( I-K ) = WORK( N+I-K+KA-M )*T -
! 317: $ WORK( I-K+KA-M )*AB( 1, I-K+KA )
! 318: AB( 1, I-K+KA ) = WORK( I-K+KA-M )*T +
! 319: $ WORK( N+I-K+KA-M )*AB( 1, I-K+KA )
! 320: RA1 = RA
! 321: END IF
! 322: END IF
! 323: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
! 324: NR = ( N-J2+KA ) / KA1
! 325: J1 = J2 + ( NR-1 )*KA1
! 326: IF( UPDATE ) THEN
! 327: J2T = MAX( J2, I+2*KA-K+1 )
! 328: ELSE
! 329: J2T = J2
! 330: END IF
! 331: NRT = ( N-J2T+KA ) / KA1
! 332: DO 90 J = J2T, J1, KA1
! 333: *
! 334: * create nonzero element a(j-ka,j+1) outside the band
! 335: * and store it in WORK(j-m)
! 336: *
! 337: WORK( J-M ) = WORK( J-M )*AB( 1, J+1 )
! 338: AB( 1, J+1 ) = WORK( N+J-M )*AB( 1, J+1 )
! 339: 90 CONTINUE
! 340: *
! 341: * generate rotations in 1st set to annihilate elements which
! 342: * have been created outside the band
! 343: *
! 344: IF( NRT.GT.0 )
! 345: $ CALL DLARGV( NRT, AB( 1, J2T ), INCA, WORK( J2T-M ), KA1,
! 346: $ WORK( N+J2T-M ), KA1 )
! 347: IF( NR.GT.0 ) THEN
! 348: *
! 349: * apply rotations in 1st set from the right
! 350: *
! 351: DO 100 L = 1, KA - 1
! 352: CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
! 353: $ AB( KA-L, J2+1 ), INCA, WORK( N+J2-M ),
! 354: $ WORK( J2-M ), KA1 )
! 355: 100 CONTINUE
! 356: *
! 357: * apply rotations in 1st set from both sides to diagonal
! 358: * blocks
! 359: *
! 360: CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
! 361: $ AB( KA, J2+1 ), INCA, WORK( N+J2-M ),
! 362: $ WORK( J2-M ), KA1 )
! 363: *
! 364: END IF
! 365: *
! 366: * start applying rotations in 1st set from the left
! 367: *
! 368: DO 110 L = KA - 1, KB - K + 1, -1
! 369: NRT = ( N-J2+L ) / KA1
! 370: IF( NRT.GT.0 )
! 371: $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
! 372: $ AB( L+1, J2+KA1-L ), INCA,
! 373: $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
! 374: 110 CONTINUE
! 375: *
! 376: IF( WANTX ) THEN
! 377: *
! 378: * post-multiply X by product of rotations in 1st set
! 379: *
! 380: DO 120 J = J2, J1, KA1
! 381: CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
! 382: $ WORK( N+J-M ), WORK( J-M ) )
! 383: 120 CONTINUE
! 384: END IF
! 385: 130 CONTINUE
! 386: *
! 387: IF( UPDATE ) THEN
! 388: IF( I2.LE.N .AND. KBT.GT.0 ) THEN
! 389: *
! 390: * create nonzero element a(i-kbt,i-kbt+ka+1) outside the
! 391: * band and store it in WORK(i-kbt)
! 392: *
! 393: WORK( I-KBT ) = -BB( KB1-KBT, I )*RA1
! 394: END IF
! 395: END IF
! 396: *
! 397: DO 170 K = KB, 1, -1
! 398: IF( UPDATE ) THEN
! 399: J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
! 400: ELSE
! 401: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
! 402: END IF
! 403: *
! 404: * finish applying rotations in 2nd set from the left
! 405: *
! 406: DO 140 L = KB - K, 1, -1
! 407: NRT = ( N-J2+KA+L ) / KA1
! 408: IF( NRT.GT.0 )
! 409: $ CALL DLARTV( NRT, AB( L, J2-L+1 ), INCA,
! 410: $ AB( L+1, J2-L+1 ), INCA, WORK( N+J2-KA ),
! 411: $ WORK( J2-KA ), KA1 )
! 412: 140 CONTINUE
! 413: NR = ( N-J2+KA ) / KA1
! 414: J1 = J2 + ( NR-1 )*KA1
! 415: DO 150 J = J1, J2, -KA1
! 416: WORK( J ) = WORK( J-KA )
! 417: WORK( N+J ) = WORK( N+J-KA )
! 418: 150 CONTINUE
! 419: DO 160 J = J2, J1, KA1
! 420: *
! 421: * create nonzero element a(j-ka,j+1) outside the band
! 422: * and store it in WORK(j)
! 423: *
! 424: WORK( J ) = WORK( J )*AB( 1, J+1 )
! 425: AB( 1, J+1 ) = WORK( N+J )*AB( 1, J+1 )
! 426: 160 CONTINUE
! 427: IF( UPDATE ) THEN
! 428: IF( I-K.LT.N-KA .AND. K.LE.KBT )
! 429: $ WORK( I-K+KA ) = WORK( I-K )
! 430: END IF
! 431: 170 CONTINUE
! 432: *
! 433: DO 210 K = KB, 1, -1
! 434: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
! 435: NR = ( N-J2+KA ) / KA1
! 436: J1 = J2 + ( NR-1 )*KA1
! 437: IF( NR.GT.0 ) THEN
! 438: *
! 439: * generate rotations in 2nd set to annihilate elements
! 440: * which have been created outside the band
! 441: *
! 442: CALL DLARGV( NR, AB( 1, J2 ), INCA, WORK( J2 ), KA1,
! 443: $ WORK( N+J2 ), KA1 )
! 444: *
! 445: * apply rotations in 2nd set from the right
! 446: *
! 447: DO 180 L = 1, KA - 1
! 448: CALL DLARTV( NR, AB( KA1-L, J2 ), INCA,
! 449: $ AB( KA-L, J2+1 ), INCA, WORK( N+J2 ),
! 450: $ WORK( J2 ), KA1 )
! 451: 180 CONTINUE
! 452: *
! 453: * apply rotations in 2nd set from both sides to diagonal
! 454: * blocks
! 455: *
! 456: CALL DLAR2V( NR, AB( KA1, J2 ), AB( KA1, J2+1 ),
! 457: $ AB( KA, J2+1 ), INCA, WORK( N+J2 ),
! 458: $ WORK( J2 ), KA1 )
! 459: *
! 460: END IF
! 461: *
! 462: * start applying rotations in 2nd set from the left
! 463: *
! 464: DO 190 L = KA - 1, KB - K + 1, -1
! 465: NRT = ( N-J2+L ) / KA1
! 466: IF( NRT.GT.0 )
! 467: $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
! 468: $ AB( L+1, J2+KA1-L ), INCA, WORK( N+J2 ),
! 469: $ WORK( J2 ), KA1 )
! 470: 190 CONTINUE
! 471: *
! 472: IF( WANTX ) THEN
! 473: *
! 474: * post-multiply X by product of rotations in 2nd set
! 475: *
! 476: DO 200 J = J2, J1, KA1
! 477: CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
! 478: $ WORK( N+J ), WORK( J ) )
! 479: 200 CONTINUE
! 480: END IF
! 481: 210 CONTINUE
! 482: *
! 483: DO 230 K = 1, KB - 1
! 484: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
! 485: *
! 486: * finish applying rotations in 1st set from the left
! 487: *
! 488: DO 220 L = KB - K, 1, -1
! 489: NRT = ( N-J2+L ) / KA1
! 490: IF( NRT.GT.0 )
! 491: $ CALL DLARTV( NRT, AB( L, J2+KA1-L ), INCA,
! 492: $ AB( L+1, J2+KA1-L ), INCA,
! 493: $ WORK( N+J2-M ), WORK( J2-M ), KA1 )
! 494: 220 CONTINUE
! 495: 230 CONTINUE
! 496: *
! 497: IF( KB.GT.1 ) THEN
! 498: DO 240 J = N - 1, I - KB + 2*KA + 1, -1
! 499: WORK( N+J-M ) = WORK( N+J-KA-M )
! 500: WORK( J-M ) = WORK( J-KA-M )
! 501: 240 CONTINUE
! 502: END IF
! 503: *
! 504: ELSE
! 505: *
! 506: * Transform A, working with the lower triangle
! 507: *
! 508: IF( UPDATE ) THEN
! 509: *
! 510: * Form inv(S(i))**T * A * inv(S(i))
! 511: *
! 512: BII = BB( 1, I )
! 513: DO 250 J = I, I1
! 514: AB( J-I+1, I ) = AB( J-I+1, I ) / BII
! 515: 250 CONTINUE
! 516: DO 260 J = MAX( 1, I-KA ), I
! 517: AB( I-J+1, J ) = AB( I-J+1, J ) / BII
! 518: 260 CONTINUE
! 519: DO 290 K = I - KBT, I - 1
! 520: DO 270 J = I - KBT, K
! 521: AB( K-J+1, J ) = AB( K-J+1, J ) -
! 522: $ BB( I-J+1, J )*AB( I-K+1, K ) -
! 523: $ BB( I-K+1, K )*AB( I-J+1, J ) +
! 524: $ AB( 1, I )*BB( I-J+1, J )*
! 525: $ BB( I-K+1, K )
! 526: 270 CONTINUE
! 527: DO 280 J = MAX( 1, I-KA ), I - KBT - 1
! 528: AB( K-J+1, J ) = AB( K-J+1, J ) -
! 529: $ BB( I-K+1, K )*AB( I-J+1, J )
! 530: 280 CONTINUE
! 531: 290 CONTINUE
! 532: DO 310 J = I, I1
! 533: DO 300 K = MAX( J-KA, I-KBT ), I - 1
! 534: AB( J-K+1, K ) = AB( J-K+1, K ) -
! 535: $ BB( I-K+1, K )*AB( J-I+1, I )
! 536: 300 CONTINUE
! 537: 310 CONTINUE
! 538: *
! 539: IF( WANTX ) THEN
! 540: *
! 541: * post-multiply X by inv(S(i))
! 542: *
! 543: CALL DSCAL( N-M, ONE / BII, X( M+1, I ), 1 )
! 544: IF( KBT.GT.0 )
! 545: $ CALL DGER( N-M, KBT, -ONE, X( M+1, I ), 1,
! 546: $ BB( KBT+1, I-KBT ), LDBB-1,
! 547: $ X( M+1, I-KBT ), LDX )
! 548: END IF
! 549: *
! 550: * store a(i1,i) in RA1 for use in next loop over K
! 551: *
! 552: RA1 = AB( I1-I+1, I )
! 553: END IF
! 554: *
! 555: * Generate and apply vectors of rotations to chase all the
! 556: * existing bulges KA positions down toward the bottom of the
! 557: * band
! 558: *
! 559: DO 360 K = 1, KB - 1
! 560: IF( UPDATE ) THEN
! 561: *
! 562: * Determine the rotations which would annihilate the bulge
! 563: * which has in theory just been created
! 564: *
! 565: IF( I-K+KA.LT.N .AND. I-K.GT.1 ) THEN
! 566: *
! 567: * generate rotation to annihilate a(i-k+ka+1,i)
! 568: *
! 569: CALL DLARTG( AB( KA1-K, I ), RA1, WORK( N+I-K+KA-M ),
! 570: $ WORK( I-K+KA-M ), RA )
! 571: *
! 572: * create nonzero element a(i-k+ka+1,i-k) outside the
! 573: * band and store it in WORK(i-k)
! 574: *
! 575: T = -BB( K+1, I-K )*RA1
! 576: WORK( I-K ) = WORK( N+I-K+KA-M )*T -
! 577: $ WORK( I-K+KA-M )*AB( KA1, I-K )
! 578: AB( KA1, I-K ) = WORK( I-K+KA-M )*T +
! 579: $ WORK( N+I-K+KA-M )*AB( KA1, I-K )
! 580: RA1 = RA
! 581: END IF
! 582: END IF
! 583: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
! 584: NR = ( N-J2+KA ) / KA1
! 585: J1 = J2 + ( NR-1 )*KA1
! 586: IF( UPDATE ) THEN
! 587: J2T = MAX( J2, I+2*KA-K+1 )
! 588: ELSE
! 589: J2T = J2
! 590: END IF
! 591: NRT = ( N-J2T+KA ) / KA1
! 592: DO 320 J = J2T, J1, KA1
! 593: *
! 594: * create nonzero element a(j+1,j-ka) outside the band
! 595: * and store it in WORK(j-m)
! 596: *
! 597: WORK( J-M ) = WORK( J-M )*AB( KA1, J-KA+1 )
! 598: AB( KA1, J-KA+1 ) = WORK( N+J-M )*AB( KA1, J-KA+1 )
! 599: 320 CONTINUE
! 600: *
! 601: * generate rotations in 1st set to annihilate elements which
! 602: * have been created outside the band
! 603: *
! 604: IF( NRT.GT.0 )
! 605: $ CALL DLARGV( NRT, AB( KA1, J2T-KA ), INCA, WORK( J2T-M ),
! 606: $ KA1, WORK( N+J2T-M ), KA1 )
! 607: IF( NR.GT.0 ) THEN
! 608: *
! 609: * apply rotations in 1st set from the left
! 610: *
! 611: DO 330 L = 1, KA - 1
! 612: CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
! 613: $ AB( L+2, J2-L ), INCA, WORK( N+J2-M ),
! 614: $ WORK( J2-M ), KA1 )
! 615: 330 CONTINUE
! 616: *
! 617: * apply rotations in 1st set from both sides to diagonal
! 618: * blocks
! 619: *
! 620: CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
! 621: $ INCA, WORK( N+J2-M ), WORK( J2-M ), KA1 )
! 622: *
! 623: END IF
! 624: *
! 625: * start applying rotations in 1st set from the right
! 626: *
! 627: DO 340 L = KA - 1, KB - K + 1, -1
! 628: NRT = ( N-J2+L ) / KA1
! 629: IF( NRT.GT.0 )
! 630: $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
! 631: $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
! 632: $ WORK( J2-M ), KA1 )
! 633: 340 CONTINUE
! 634: *
! 635: IF( WANTX ) THEN
! 636: *
! 637: * post-multiply X by product of rotations in 1st set
! 638: *
! 639: DO 350 J = J2, J1, KA1
! 640: CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
! 641: $ WORK( N+J-M ), WORK( J-M ) )
! 642: 350 CONTINUE
! 643: END IF
! 644: 360 CONTINUE
! 645: *
! 646: IF( UPDATE ) THEN
! 647: IF( I2.LE.N .AND. KBT.GT.0 ) THEN
! 648: *
! 649: * create nonzero element a(i-kbt+ka+1,i-kbt) outside the
! 650: * band and store it in WORK(i-kbt)
! 651: *
! 652: WORK( I-KBT ) = -BB( KBT+1, I-KBT )*RA1
! 653: END IF
! 654: END IF
! 655: *
! 656: DO 400 K = KB, 1, -1
! 657: IF( UPDATE ) THEN
! 658: J2 = I - K - 1 + MAX( 2, K-I0+1 )*KA1
! 659: ELSE
! 660: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
! 661: END IF
! 662: *
! 663: * finish applying rotations in 2nd set from the right
! 664: *
! 665: DO 370 L = KB - K, 1, -1
! 666: NRT = ( N-J2+KA+L ) / KA1
! 667: IF( NRT.GT.0 )
! 668: $ CALL DLARTV( NRT, AB( KA1-L+1, J2-KA ), INCA,
! 669: $ AB( KA1-L, J2-KA+1 ), INCA,
! 670: $ WORK( N+J2-KA ), WORK( J2-KA ), KA1 )
! 671: 370 CONTINUE
! 672: NR = ( N-J2+KA ) / KA1
! 673: J1 = J2 + ( NR-1 )*KA1
! 674: DO 380 J = J1, J2, -KA1
! 675: WORK( J ) = WORK( J-KA )
! 676: WORK( N+J ) = WORK( N+J-KA )
! 677: 380 CONTINUE
! 678: DO 390 J = J2, J1, KA1
! 679: *
! 680: * create nonzero element a(j+1,j-ka) outside the band
! 681: * and store it in WORK(j)
! 682: *
! 683: WORK( J ) = WORK( J )*AB( KA1, J-KA+1 )
! 684: AB( KA1, J-KA+1 ) = WORK( N+J )*AB( KA1, J-KA+1 )
! 685: 390 CONTINUE
! 686: IF( UPDATE ) THEN
! 687: IF( I-K.LT.N-KA .AND. K.LE.KBT )
! 688: $ WORK( I-K+KA ) = WORK( I-K )
! 689: END IF
! 690: 400 CONTINUE
! 691: *
! 692: DO 440 K = KB, 1, -1
! 693: J2 = I - K - 1 + MAX( 1, K-I0+1 )*KA1
! 694: NR = ( N-J2+KA ) / KA1
! 695: J1 = J2 + ( NR-1 )*KA1
! 696: IF( NR.GT.0 ) THEN
! 697: *
! 698: * generate rotations in 2nd set to annihilate elements
! 699: * which have been created outside the band
! 700: *
! 701: CALL DLARGV( NR, AB( KA1, J2-KA ), INCA, WORK( J2 ), KA1,
! 702: $ WORK( N+J2 ), KA1 )
! 703: *
! 704: * apply rotations in 2nd set from the left
! 705: *
! 706: DO 410 L = 1, KA - 1
! 707: CALL DLARTV( NR, AB( L+1, J2-L ), INCA,
! 708: $ AB( L+2, J2-L ), INCA, WORK( N+J2 ),
! 709: $ WORK( J2 ), KA1 )
! 710: 410 CONTINUE
! 711: *
! 712: * apply rotations in 2nd set from both sides to diagonal
! 713: * blocks
! 714: *
! 715: CALL DLAR2V( NR, AB( 1, J2 ), AB( 1, J2+1 ), AB( 2, J2 ),
! 716: $ INCA, WORK( N+J2 ), WORK( J2 ), KA1 )
! 717: *
! 718: END IF
! 719: *
! 720: * start applying rotations in 2nd set from the right
! 721: *
! 722: DO 420 L = KA - 1, KB - K + 1, -1
! 723: NRT = ( N-J2+L ) / KA1
! 724: IF( NRT.GT.0 )
! 725: $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
! 726: $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2 ),
! 727: $ WORK( J2 ), KA1 )
! 728: 420 CONTINUE
! 729: *
! 730: IF( WANTX ) THEN
! 731: *
! 732: * post-multiply X by product of rotations in 2nd set
! 733: *
! 734: DO 430 J = J2, J1, KA1
! 735: CALL DROT( N-M, X( M+1, J ), 1, X( M+1, J+1 ), 1,
! 736: $ WORK( N+J ), WORK( J ) )
! 737: 430 CONTINUE
! 738: END IF
! 739: 440 CONTINUE
! 740: *
! 741: DO 460 K = 1, KB - 1
! 742: J2 = I - K - 1 + MAX( 1, K-I0+2 )*KA1
! 743: *
! 744: * finish applying rotations in 1st set from the right
! 745: *
! 746: DO 450 L = KB - K, 1, -1
! 747: NRT = ( N-J2+L ) / KA1
! 748: IF( NRT.GT.0 )
! 749: $ CALL DLARTV( NRT, AB( KA1-L+1, J2 ), INCA,
! 750: $ AB( KA1-L, J2+1 ), INCA, WORK( N+J2-M ),
! 751: $ WORK( J2-M ), KA1 )
! 752: 450 CONTINUE
! 753: 460 CONTINUE
! 754: *
! 755: IF( KB.GT.1 ) THEN
! 756: DO 470 J = N - 1, I - KB + 2*KA + 1, -1
! 757: WORK( N+J-M ) = WORK( N+J-KA-M )
! 758: WORK( J-M ) = WORK( J-KA-M )
! 759: 470 CONTINUE
! 760: END IF
! 761: *
! 762: END IF
! 763: *
! 764: GO TO 10
! 765: *
! 766: 480 CONTINUE
! 767: *
! 768: * **************************** Phase 2 *****************************
! 769: *
! 770: * The logical structure of this phase is:
! 771: *
! 772: * UPDATE = .TRUE.
! 773: * DO I = 1, M
! 774: * use S(i) to update A and create a new bulge
! 775: * apply rotations to push all bulges KA positions upward
! 776: * END DO
! 777: * UPDATE = .FALSE.
! 778: * DO I = M - KA - 1, 2, -1
! 779: * apply rotations to push all bulges KA positions upward
! 780: * END DO
! 781: *
! 782: * To avoid duplicating code, the two loops are merged.
! 783: *
! 784: UPDATE = .TRUE.
! 785: I = 0
! 786: 490 CONTINUE
! 787: IF( UPDATE ) THEN
! 788: I = I + 1
! 789: KBT = MIN( KB, M-I )
! 790: I0 = I + 1
! 791: I1 = MAX( 1, I-KA )
! 792: I2 = I + KBT - KA1
! 793: IF( I.GT.M ) THEN
! 794: UPDATE = .FALSE.
! 795: I = I - 1
! 796: I0 = M + 1
! 797: IF( KA.EQ.0 )
! 798: $ RETURN
! 799: GO TO 490
! 800: END IF
! 801: ELSE
! 802: I = I - KA
! 803: IF( I.LT.2 )
! 804: $ RETURN
! 805: END IF
! 806: *
! 807: IF( I.LT.M-KBT ) THEN
! 808: NX = M
! 809: ELSE
! 810: NX = N
! 811: END IF
! 812: *
! 813: IF( UPPER ) THEN
! 814: *
! 815: * Transform A, working with the upper triangle
! 816: *
! 817: IF( UPDATE ) THEN
! 818: *
! 819: * Form inv(S(i))**T * A * inv(S(i))
! 820: *
! 821: BII = BB( KB1, I )
! 822: DO 500 J = I1, I
! 823: AB( J-I+KA1, I ) = AB( J-I+KA1, I ) / BII
! 824: 500 CONTINUE
! 825: DO 510 J = I, MIN( N, I+KA )
! 826: AB( I-J+KA1, J ) = AB( I-J+KA1, J ) / BII
! 827: 510 CONTINUE
! 828: DO 540 K = I + 1, I + KBT
! 829: DO 520 J = K, I + KBT
! 830: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
! 831: $ BB( I-J+KB1, J )*AB( I-K+KA1, K ) -
! 832: $ BB( I-K+KB1, K )*AB( I-J+KA1, J ) +
! 833: $ AB( KA1, I )*BB( I-J+KB1, J )*
! 834: $ BB( I-K+KB1, K )
! 835: 520 CONTINUE
! 836: DO 530 J = I + KBT + 1, MIN( N, I+KA )
! 837: AB( K-J+KA1, J ) = AB( K-J+KA1, J ) -
! 838: $ BB( I-K+KB1, K )*AB( I-J+KA1, J )
! 839: 530 CONTINUE
! 840: 540 CONTINUE
! 841: DO 560 J = I1, I
! 842: DO 550 K = I + 1, MIN( J+KA, I+KBT )
! 843: AB( J-K+KA1, K ) = AB( J-K+KA1, K ) -
! 844: $ BB( I-K+KB1, K )*AB( J-I+KA1, I )
! 845: 550 CONTINUE
! 846: 560 CONTINUE
! 847: *
! 848: IF( WANTX ) THEN
! 849: *
! 850: * post-multiply X by inv(S(i))
! 851: *
! 852: CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
! 853: IF( KBT.GT.0 )
! 854: $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( KB, I+1 ),
! 855: $ LDBB-1, X( 1, I+1 ), LDX )
! 856: END IF
! 857: *
! 858: * store a(i1,i) in RA1 for use in next loop over K
! 859: *
! 860: RA1 = AB( I1-I+KA1, I )
! 861: END IF
! 862: *
! 863: * Generate and apply vectors of rotations to chase all the
! 864: * existing bulges KA positions up toward the top of the band
! 865: *
! 866: DO 610 K = 1, KB - 1
! 867: IF( UPDATE ) THEN
! 868: *
! 869: * Determine the rotations which would annihilate the bulge
! 870: * which has in theory just been created
! 871: *
! 872: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
! 873: *
! 874: * generate rotation to annihilate a(i+k-ka-1,i)
! 875: *
! 876: CALL DLARTG( AB( K+1, I ), RA1, WORK( N+I+K-KA ),
! 877: $ WORK( I+K-KA ), RA )
! 878: *
! 879: * create nonzero element a(i+k-ka-1,i+k) outside the
! 880: * band and store it in WORK(m-kb+i+k)
! 881: *
! 882: T = -BB( KB1-K, I+K )*RA1
! 883: WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
! 884: $ WORK( I+K-KA )*AB( 1, I+K )
! 885: AB( 1, I+K ) = WORK( I+K-KA )*T +
! 886: $ WORK( N+I+K-KA )*AB( 1, I+K )
! 887: RA1 = RA
! 888: END IF
! 889: END IF
! 890: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
! 891: NR = ( J2+KA-1 ) / KA1
! 892: J1 = J2 - ( NR-1 )*KA1
! 893: IF( UPDATE ) THEN
! 894: J2T = MIN( J2, I-2*KA+K-1 )
! 895: ELSE
! 896: J2T = J2
! 897: END IF
! 898: NRT = ( J2T+KA-1 ) / KA1
! 899: DO 570 J = J1, J2T, KA1
! 900: *
! 901: * create nonzero element a(j-1,j+ka) outside the band
! 902: * and store it in WORK(j)
! 903: *
! 904: WORK( J ) = WORK( J )*AB( 1, J+KA-1 )
! 905: AB( 1, J+KA-1 ) = WORK( N+J )*AB( 1, J+KA-1 )
! 906: 570 CONTINUE
! 907: *
! 908: * generate rotations in 1st set to annihilate elements which
! 909: * have been created outside the band
! 910: *
! 911: IF( NRT.GT.0 )
! 912: $ CALL DLARGV( NRT, AB( 1, J1+KA ), INCA, WORK( J1 ), KA1,
! 913: $ WORK( N+J1 ), KA1 )
! 914: IF( NR.GT.0 ) THEN
! 915: *
! 916: * apply rotations in 1st set from the left
! 917: *
! 918: DO 580 L = 1, KA - 1
! 919: CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
! 920: $ AB( KA-L, J1+L ), INCA, WORK( N+J1 ),
! 921: $ WORK( J1 ), KA1 )
! 922: 580 CONTINUE
! 923: *
! 924: * apply rotations in 1st set from both sides to diagonal
! 925: * blocks
! 926: *
! 927: CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
! 928: $ AB( KA, J1 ), INCA, WORK( N+J1 ),
! 929: $ WORK( J1 ), KA1 )
! 930: *
! 931: END IF
! 932: *
! 933: * start applying rotations in 1st set from the right
! 934: *
! 935: DO 590 L = KA - 1, KB - K + 1, -1
! 936: NRT = ( J2+L-1 ) / KA1
! 937: J1T = J2 - ( NRT-1 )*KA1
! 938: IF( NRT.GT.0 )
! 939: $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
! 940: $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
! 941: $ WORK( J1T ), KA1 )
! 942: 590 CONTINUE
! 943: *
! 944: IF( WANTX ) THEN
! 945: *
! 946: * post-multiply X by product of rotations in 1st set
! 947: *
! 948: DO 600 J = J1, J2, KA1
! 949: CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
! 950: $ WORK( N+J ), WORK( J ) )
! 951: 600 CONTINUE
! 952: END IF
! 953: 610 CONTINUE
! 954: *
! 955: IF( UPDATE ) THEN
! 956: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
! 957: *
! 958: * create nonzero element a(i+kbt-ka-1,i+kbt) outside the
! 959: * band and store it in WORK(m-kb+i+kbt)
! 960: *
! 961: WORK( M-KB+I+KBT ) = -BB( KB1-KBT, I+KBT )*RA1
! 962: END IF
! 963: END IF
! 964: *
! 965: DO 650 K = KB, 1, -1
! 966: IF( UPDATE ) THEN
! 967: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
! 968: ELSE
! 969: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
! 970: END IF
! 971: *
! 972: * finish applying rotations in 2nd set from the right
! 973: *
! 974: DO 620 L = KB - K, 1, -1
! 975: NRT = ( J2+KA+L-1 ) / KA1
! 976: J1T = J2 - ( NRT-1 )*KA1
! 977: IF( NRT.GT.0 )
! 978: $ CALL DLARTV( NRT, AB( L, J1T+KA ), INCA,
! 979: $ AB( L+1, J1T+KA-1 ), INCA,
! 980: $ WORK( N+M-KB+J1T+KA ),
! 981: $ WORK( M-KB+J1T+KA ), KA1 )
! 982: 620 CONTINUE
! 983: NR = ( J2+KA-1 ) / KA1
! 984: J1 = J2 - ( NR-1 )*KA1
! 985: DO 630 J = J1, J2, KA1
! 986: WORK( M-KB+J ) = WORK( M-KB+J+KA )
! 987: WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
! 988: 630 CONTINUE
! 989: DO 640 J = J1, J2, KA1
! 990: *
! 991: * create nonzero element a(j-1,j+ka) outside the band
! 992: * and store it in WORK(m-kb+j)
! 993: *
! 994: WORK( M-KB+J ) = WORK( M-KB+J )*AB( 1, J+KA-1 )
! 995: AB( 1, J+KA-1 ) = WORK( N+M-KB+J )*AB( 1, J+KA-1 )
! 996: 640 CONTINUE
! 997: IF( UPDATE ) THEN
! 998: IF( I+K.GT.KA1 .AND. K.LE.KBT )
! 999: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
! 1000: END IF
! 1001: 650 CONTINUE
! 1002: *
! 1003: DO 690 K = KB, 1, -1
! 1004: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
! 1005: NR = ( J2+KA-1 ) / KA1
! 1006: J1 = J2 - ( NR-1 )*KA1
! 1007: IF( NR.GT.0 ) THEN
! 1008: *
! 1009: * generate rotations in 2nd set to annihilate elements
! 1010: * which have been created outside the band
! 1011: *
! 1012: CALL DLARGV( NR, AB( 1, J1+KA ), INCA, WORK( M-KB+J1 ),
! 1013: $ KA1, WORK( N+M-KB+J1 ), KA1 )
! 1014: *
! 1015: * apply rotations in 2nd set from the left
! 1016: *
! 1017: DO 660 L = 1, KA - 1
! 1018: CALL DLARTV( NR, AB( KA1-L, J1+L ), INCA,
! 1019: $ AB( KA-L, J1+L ), INCA,
! 1020: $ WORK( N+M-KB+J1 ), WORK( M-KB+J1 ), KA1 )
! 1021: 660 CONTINUE
! 1022: *
! 1023: * apply rotations in 2nd set from both sides to diagonal
! 1024: * blocks
! 1025: *
! 1026: CALL DLAR2V( NR, AB( KA1, J1 ), AB( KA1, J1-1 ),
! 1027: $ AB( KA, J1 ), INCA, WORK( N+M-KB+J1 ),
! 1028: $ WORK( M-KB+J1 ), KA1 )
! 1029: *
! 1030: END IF
! 1031: *
! 1032: * start applying rotations in 2nd set from the right
! 1033: *
! 1034: DO 670 L = KA - 1, KB - K + 1, -1
! 1035: NRT = ( J2+L-1 ) / KA1
! 1036: J1T = J2 - ( NRT-1 )*KA1
! 1037: IF( NRT.GT.0 )
! 1038: $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
! 1039: $ AB( L+1, J1T-1 ), INCA,
! 1040: $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
! 1041: $ KA1 )
! 1042: 670 CONTINUE
! 1043: *
! 1044: IF( WANTX ) THEN
! 1045: *
! 1046: * post-multiply X by product of rotations in 2nd set
! 1047: *
! 1048: DO 680 J = J1, J2, KA1
! 1049: CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
! 1050: $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
! 1051: 680 CONTINUE
! 1052: END IF
! 1053: 690 CONTINUE
! 1054: *
! 1055: DO 710 K = 1, KB - 1
! 1056: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
! 1057: *
! 1058: * finish applying rotations in 1st set from the right
! 1059: *
! 1060: DO 700 L = KB - K, 1, -1
! 1061: NRT = ( J2+L-1 ) / KA1
! 1062: J1T = J2 - ( NRT-1 )*KA1
! 1063: IF( NRT.GT.0 )
! 1064: $ CALL DLARTV( NRT, AB( L, J1T ), INCA,
! 1065: $ AB( L+1, J1T-1 ), INCA, WORK( N+J1T ),
! 1066: $ WORK( J1T ), KA1 )
! 1067: 700 CONTINUE
! 1068: 710 CONTINUE
! 1069: *
! 1070: IF( KB.GT.1 ) THEN
! 1071: DO 720 J = 2, MIN( I+KB, M ) - 2*KA - 1
! 1072: WORK( N+J ) = WORK( N+J+KA )
! 1073: WORK( J ) = WORK( J+KA )
! 1074: 720 CONTINUE
! 1075: END IF
! 1076: *
! 1077: ELSE
! 1078: *
! 1079: * Transform A, working with the lower triangle
! 1080: *
! 1081: IF( UPDATE ) THEN
! 1082: *
! 1083: * Form inv(S(i))**T * A * inv(S(i))
! 1084: *
! 1085: BII = BB( 1, I )
! 1086: DO 730 J = I1, I
! 1087: AB( I-J+1, J ) = AB( I-J+1, J ) / BII
! 1088: 730 CONTINUE
! 1089: DO 740 J = I, MIN( N, I+KA )
! 1090: AB( J-I+1, I ) = AB( J-I+1, I ) / BII
! 1091: 740 CONTINUE
! 1092: DO 770 K = I + 1, I + KBT
! 1093: DO 750 J = K, I + KBT
! 1094: AB( J-K+1, K ) = AB( J-K+1, K ) -
! 1095: $ BB( J-I+1, I )*AB( K-I+1, I ) -
! 1096: $ BB( K-I+1, I )*AB( J-I+1, I ) +
! 1097: $ AB( 1, I )*BB( J-I+1, I )*
! 1098: $ BB( K-I+1, I )
! 1099: 750 CONTINUE
! 1100: DO 760 J = I + KBT + 1, MIN( N, I+KA )
! 1101: AB( J-K+1, K ) = AB( J-K+1, K ) -
! 1102: $ BB( K-I+1, I )*AB( J-I+1, I )
! 1103: 760 CONTINUE
! 1104: 770 CONTINUE
! 1105: DO 790 J = I1, I
! 1106: DO 780 K = I + 1, MIN( J+KA, I+KBT )
! 1107: AB( K-J+1, J ) = AB( K-J+1, J ) -
! 1108: $ BB( K-I+1, I )*AB( I-J+1, J )
! 1109: 780 CONTINUE
! 1110: 790 CONTINUE
! 1111: *
! 1112: IF( WANTX ) THEN
! 1113: *
! 1114: * post-multiply X by inv(S(i))
! 1115: *
! 1116: CALL DSCAL( NX, ONE / BII, X( 1, I ), 1 )
! 1117: IF( KBT.GT.0 )
! 1118: $ CALL DGER( NX, KBT, -ONE, X( 1, I ), 1, BB( 2, I ), 1,
! 1119: $ X( 1, I+1 ), LDX )
! 1120: END IF
! 1121: *
! 1122: * store a(i,i1) in RA1 for use in next loop over K
! 1123: *
! 1124: RA1 = AB( I-I1+1, I1 )
! 1125: END IF
! 1126: *
! 1127: * Generate and apply vectors of rotations to chase all the
! 1128: * existing bulges KA positions up toward the top of the band
! 1129: *
! 1130: DO 840 K = 1, KB - 1
! 1131: IF( UPDATE ) THEN
! 1132: *
! 1133: * Determine the rotations which would annihilate the bulge
! 1134: * which has in theory just been created
! 1135: *
! 1136: IF( I+K-KA1.GT.0 .AND. I+K.LT.M ) THEN
! 1137: *
! 1138: * generate rotation to annihilate a(i,i+k-ka-1)
! 1139: *
! 1140: CALL DLARTG( AB( KA1-K, I+K-KA ), RA1,
! 1141: $ WORK( N+I+K-KA ), WORK( I+K-KA ), RA )
! 1142: *
! 1143: * create nonzero element a(i+k,i+k-ka-1) outside the
! 1144: * band and store it in WORK(m-kb+i+k)
! 1145: *
! 1146: T = -BB( K+1, I )*RA1
! 1147: WORK( M-KB+I+K ) = WORK( N+I+K-KA )*T -
! 1148: $ WORK( I+K-KA )*AB( KA1, I+K-KA )
! 1149: AB( KA1, I+K-KA ) = WORK( I+K-KA )*T +
! 1150: $ WORK( N+I+K-KA )*AB( KA1, I+K-KA )
! 1151: RA1 = RA
! 1152: END IF
! 1153: END IF
! 1154: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
! 1155: NR = ( J2+KA-1 ) / KA1
! 1156: J1 = J2 - ( NR-1 )*KA1
! 1157: IF( UPDATE ) THEN
! 1158: J2T = MIN( J2, I-2*KA+K-1 )
! 1159: ELSE
! 1160: J2T = J2
! 1161: END IF
! 1162: NRT = ( J2T+KA-1 ) / KA1
! 1163: DO 800 J = J1, J2T, KA1
! 1164: *
! 1165: * create nonzero element a(j+ka,j-1) outside the band
! 1166: * and store it in WORK(j)
! 1167: *
! 1168: WORK( J ) = WORK( J )*AB( KA1, J-1 )
! 1169: AB( KA1, J-1 ) = WORK( N+J )*AB( KA1, J-1 )
! 1170: 800 CONTINUE
! 1171: *
! 1172: * generate rotations in 1st set to annihilate elements which
! 1173: * have been created outside the band
! 1174: *
! 1175: IF( NRT.GT.0 )
! 1176: $ CALL DLARGV( NRT, AB( KA1, J1 ), INCA, WORK( J1 ), KA1,
! 1177: $ WORK( N+J1 ), KA1 )
! 1178: IF( NR.GT.0 ) THEN
! 1179: *
! 1180: * apply rotations in 1st set from the right
! 1181: *
! 1182: DO 810 L = 1, KA - 1
! 1183: CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
! 1184: $ INCA, WORK( N+J1 ), WORK( J1 ), KA1 )
! 1185: 810 CONTINUE
! 1186: *
! 1187: * apply rotations in 1st set from both sides to diagonal
! 1188: * blocks
! 1189: *
! 1190: CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
! 1191: $ AB( 2, J1-1 ), INCA, WORK( N+J1 ),
! 1192: $ WORK( J1 ), KA1 )
! 1193: *
! 1194: END IF
! 1195: *
! 1196: * start applying rotations in 1st set from the left
! 1197: *
! 1198: DO 820 L = KA - 1, KB - K + 1, -1
! 1199: NRT = ( J2+L-1 ) / KA1
! 1200: J1T = J2 - ( NRT-1 )*KA1
! 1201: IF( NRT.GT.0 )
! 1202: $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
! 1203: $ AB( KA1-L, J1T-KA1+L ), INCA,
! 1204: $ WORK( N+J1T ), WORK( J1T ), KA1 )
! 1205: 820 CONTINUE
! 1206: *
! 1207: IF( WANTX ) THEN
! 1208: *
! 1209: * post-multiply X by product of rotations in 1st set
! 1210: *
! 1211: DO 830 J = J1, J2, KA1
! 1212: CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
! 1213: $ WORK( N+J ), WORK( J ) )
! 1214: 830 CONTINUE
! 1215: END IF
! 1216: 840 CONTINUE
! 1217: *
! 1218: IF( UPDATE ) THEN
! 1219: IF( I2.GT.0 .AND. KBT.GT.0 ) THEN
! 1220: *
! 1221: * create nonzero element a(i+kbt,i+kbt-ka-1) outside the
! 1222: * band and store it in WORK(m-kb+i+kbt)
! 1223: *
! 1224: WORK( M-KB+I+KBT ) = -BB( KBT+1, I )*RA1
! 1225: END IF
! 1226: END IF
! 1227: *
! 1228: DO 880 K = KB, 1, -1
! 1229: IF( UPDATE ) THEN
! 1230: J2 = I + K + 1 - MAX( 2, K+I0-M )*KA1
! 1231: ELSE
! 1232: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
! 1233: END IF
! 1234: *
! 1235: * finish applying rotations in 2nd set from the left
! 1236: *
! 1237: DO 850 L = KB - K, 1, -1
! 1238: NRT = ( J2+KA+L-1 ) / KA1
! 1239: J1T = J2 - ( NRT-1 )*KA1
! 1240: IF( NRT.GT.0 )
! 1241: $ CALL DLARTV( NRT, AB( KA1-L+1, J1T+L-1 ), INCA,
! 1242: $ AB( KA1-L, J1T+L-1 ), INCA,
! 1243: $ WORK( N+M-KB+J1T+KA ),
! 1244: $ WORK( M-KB+J1T+KA ), KA1 )
! 1245: 850 CONTINUE
! 1246: NR = ( J2+KA-1 ) / KA1
! 1247: J1 = J2 - ( NR-1 )*KA1
! 1248: DO 860 J = J1, J2, KA1
! 1249: WORK( M-KB+J ) = WORK( M-KB+J+KA )
! 1250: WORK( N+M-KB+J ) = WORK( N+M-KB+J+KA )
! 1251: 860 CONTINUE
! 1252: DO 870 J = J1, J2, KA1
! 1253: *
! 1254: * create nonzero element a(j+ka,j-1) outside the band
! 1255: * and store it in WORK(m-kb+j)
! 1256: *
! 1257: WORK( M-KB+J ) = WORK( M-KB+J )*AB( KA1, J-1 )
! 1258: AB( KA1, J-1 ) = WORK( N+M-KB+J )*AB( KA1, J-1 )
! 1259: 870 CONTINUE
! 1260: IF( UPDATE ) THEN
! 1261: IF( I+K.GT.KA1 .AND. K.LE.KBT )
! 1262: $ WORK( M-KB+I+K-KA ) = WORK( M-KB+I+K )
! 1263: END IF
! 1264: 880 CONTINUE
! 1265: *
! 1266: DO 920 K = KB, 1, -1
! 1267: J2 = I + K + 1 - MAX( 1, K+I0-M )*KA1
! 1268: NR = ( J2+KA-1 ) / KA1
! 1269: J1 = J2 - ( NR-1 )*KA1
! 1270: IF( NR.GT.0 ) THEN
! 1271: *
! 1272: * generate rotations in 2nd set to annihilate elements
! 1273: * which have been created outside the band
! 1274: *
! 1275: CALL DLARGV( NR, AB( KA1, J1 ), INCA, WORK( M-KB+J1 ),
! 1276: $ KA1, WORK( N+M-KB+J1 ), KA1 )
! 1277: *
! 1278: * apply rotations in 2nd set from the right
! 1279: *
! 1280: DO 890 L = 1, KA - 1
! 1281: CALL DLARTV( NR, AB( L+1, J1 ), INCA, AB( L+2, J1-1 ),
! 1282: $ INCA, WORK( N+M-KB+J1 ), WORK( M-KB+J1 ),
! 1283: $ KA1 )
! 1284: 890 CONTINUE
! 1285: *
! 1286: * apply rotations in 2nd set from both sides to diagonal
! 1287: * blocks
! 1288: *
! 1289: CALL DLAR2V( NR, AB( 1, J1 ), AB( 1, J1-1 ),
! 1290: $ AB( 2, J1-1 ), INCA, WORK( N+M-KB+J1 ),
! 1291: $ WORK( M-KB+J1 ), KA1 )
! 1292: *
! 1293: END IF
! 1294: *
! 1295: * start applying rotations in 2nd set from the left
! 1296: *
! 1297: DO 900 L = KA - 1, KB - K + 1, -1
! 1298: NRT = ( J2+L-1 ) / KA1
! 1299: J1T = J2 - ( NRT-1 )*KA1
! 1300: IF( NRT.GT.0 )
! 1301: $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
! 1302: $ AB( KA1-L, J1T-KA1+L ), INCA,
! 1303: $ WORK( N+M-KB+J1T ), WORK( M-KB+J1T ),
! 1304: $ KA1 )
! 1305: 900 CONTINUE
! 1306: *
! 1307: IF( WANTX ) THEN
! 1308: *
! 1309: * post-multiply X by product of rotations in 2nd set
! 1310: *
! 1311: DO 910 J = J1, J2, KA1
! 1312: CALL DROT( NX, X( 1, J ), 1, X( 1, J-1 ), 1,
! 1313: $ WORK( N+M-KB+J ), WORK( M-KB+J ) )
! 1314: 910 CONTINUE
! 1315: END IF
! 1316: 920 CONTINUE
! 1317: *
! 1318: DO 940 K = 1, KB - 1
! 1319: J2 = I + K + 1 - MAX( 1, K+I0-M+1 )*KA1
! 1320: *
! 1321: * finish applying rotations in 1st set from the left
! 1322: *
! 1323: DO 930 L = KB - K, 1, -1
! 1324: NRT = ( J2+L-1 ) / KA1
! 1325: J1T = J2 - ( NRT-1 )*KA1
! 1326: IF( NRT.GT.0 )
! 1327: $ CALL DLARTV( NRT, AB( KA1-L+1, J1T-KA1+L ), INCA,
! 1328: $ AB( KA1-L, J1T-KA1+L ), INCA,
! 1329: $ WORK( N+J1T ), WORK( J1T ), KA1 )
! 1330: 930 CONTINUE
! 1331: 940 CONTINUE
! 1332: *
! 1333: IF( KB.GT.1 ) THEN
! 1334: DO 950 J = 2, MIN( I+KB, M ) - 2*KA - 1
! 1335: WORK( N+J ) = WORK( N+J+KA )
! 1336: WORK( J ) = WORK( J+KA )
! 1337: 950 CONTINUE
! 1338: END IF
! 1339: *
! 1340: END IF
! 1341: *
! 1342: GO TO 490
! 1343: *
! 1344: * End of DSBGST
! 1345: *
! 1346: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsbgst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack