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