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