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