1: SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
2: $ RSCALE, 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 JOB
11: INTEGER IHI, ILO, INFO, LDA, LDB, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), LSCALE( * ),
15: $ RSCALE( * ), WORK( * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DGGBAL balances a pair of general real matrices (A,B). This
22: * involves, first, permuting A and B by similarity transformations to
23: * isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
24: * elements on the diagonal; and second, applying a diagonal similarity
25: * transformation to rows and columns ILO to IHI to make the rows
26: * and columns as close in norm as possible. Both steps are optional.
27: *
28: * Balancing may reduce the 1-norm of the matrices, and improve the
29: * accuracy of the computed eigenvalues and/or eigenvectors in the
30: * generalized eigenvalue problem A*x = lambda*B*x.
31: *
32: * Arguments
33: * =========
34: *
35: * JOB (input) CHARACTER*1
36: * Specifies the operations to be performed on A and B:
37: * = 'N': none: simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
38: * and RSCALE(I) = 1.0 for i = 1,...,N.
39: * = 'P': permute only;
40: * = 'S': scale only;
41: * = 'B': both permute and scale.
42: *
43: * N (input) INTEGER
44: * The order of the matrices A and B. N >= 0.
45: *
46: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
47: * On entry, the input matrix A.
48: * On exit, A is overwritten by the balanced matrix.
49: * If JOB = 'N', A is not referenced.
50: *
51: * LDA (input) INTEGER
52: * The leading dimension of the array A. LDA >= max(1,N).
53: *
54: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
55: * On entry, the input matrix B.
56: * On exit, B is overwritten by the balanced matrix.
57: * If JOB = 'N', B is not referenced.
58: *
59: * LDB (input) INTEGER
60: * The leading dimension of the array B. LDB >= max(1,N).
61: *
62: * ILO (output) INTEGER
63: * IHI (output) INTEGER
64: * ILO and IHI are set to integers such that on exit
65: * A(i,j) = 0 and B(i,j) = 0 if i > j and
66: * j = 1,...,ILO-1 or i = IHI+1,...,N.
67: * If JOB = 'N' or 'S', ILO = 1 and IHI = N.
68: *
69: * LSCALE (output) DOUBLE PRECISION array, dimension (N)
70: * Details of the permutations and scaling factors applied
71: * to the left side of A and B. If P(j) is the index of the
72: * row interchanged with row j, and D(j)
73: * is the scaling factor applied to row j, then
74: * LSCALE(j) = P(j) for J = 1,...,ILO-1
75: * = D(j) for J = ILO,...,IHI
76: * = P(j) for J = IHI+1,...,N.
77: * The order in which the interchanges are made is N to IHI+1,
78: * then 1 to ILO-1.
79: *
80: * RSCALE (output) DOUBLE PRECISION array, dimension (N)
81: * Details of the permutations and scaling factors applied
82: * to the right side of A and B. If P(j) is the index of the
83: * column interchanged with column j, and D(j)
84: * is the scaling factor applied to column j, then
85: * LSCALE(j) = P(j) for J = 1,...,ILO-1
86: * = D(j) for J = ILO,...,IHI
87: * = P(j) for J = IHI+1,...,N.
88: * The order in which the interchanges are made is N to IHI+1,
89: * then 1 to ILO-1.
90: *
91: * WORK (workspace) REAL array, dimension (lwork)
92: * lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
93: * at least 1 when JOB = 'N' or 'P'.
94: *
95: * INFO (output) INTEGER
96: * = 0: successful exit
97: * < 0: if INFO = -i, the i-th argument had an illegal value.
98: *
99: * Further Details
100: * ===============
101: *
102: * See R.C. WARD, Balancing the generalized eigenvalue problem,
103: * SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
104: *
105: * =====================================================================
106: *
107: * .. Parameters ..
108: DOUBLE PRECISION ZERO, HALF, ONE
109: PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
110: DOUBLE PRECISION THREE, SCLFAC
111: PARAMETER ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
112: * ..
113: * .. Local Scalars ..
114: INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
115: $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
116: $ M, NR, NRP2
117: DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
118: $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
119: $ SFMIN, SUM, T, TA, TB, TC
120: * ..
121: * .. External Functions ..
122: LOGICAL LSAME
123: INTEGER IDAMAX
124: DOUBLE PRECISION DDOT, DLAMCH
125: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
126: * ..
127: * .. External Subroutines ..
128: EXTERNAL DAXPY, DSCAL, DSWAP, XERBLA
129: * ..
130: * .. Intrinsic Functions ..
131: INTRINSIC ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
132: * ..
133: * .. Executable Statements ..
134: *
135: * Test the input parameters
136: *
137: INFO = 0
138: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
139: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
140: INFO = -1
141: ELSE IF( N.LT.0 ) THEN
142: INFO = -2
143: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
144: INFO = -4
145: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
146: INFO = -6
147: END IF
148: IF( INFO.NE.0 ) THEN
149: CALL XERBLA( 'DGGBAL', -INFO )
150: RETURN
151: END IF
152: *
153: * Quick return if possible
154: *
155: IF( N.EQ.0 ) THEN
156: ILO = 1
157: IHI = N
158: RETURN
159: END IF
160: *
161: IF( N.EQ.1 ) THEN
162: ILO = 1
163: IHI = N
164: LSCALE( 1 ) = ONE
165: RSCALE( 1 ) = ONE
166: RETURN
167: END IF
168: *
169: IF( LSAME( JOB, 'N' ) ) THEN
170: ILO = 1
171: IHI = N
172: DO 10 I = 1, N
173: LSCALE( I ) = ONE
174: RSCALE( I ) = ONE
175: 10 CONTINUE
176: RETURN
177: END IF
178: *
179: K = 1
180: L = N
181: IF( LSAME( JOB, 'S' ) )
182: $ GO TO 190
183: *
184: GO TO 30
185: *
186: * Permute the matrices A and B to isolate the eigenvalues.
187: *
188: * Find row with one nonzero in columns 1 through L
189: *
190: 20 CONTINUE
191: L = LM1
192: IF( L.NE.1 )
193: $ GO TO 30
194: *
195: RSCALE( 1 ) = ONE
196: LSCALE( 1 ) = ONE
197: GO TO 190
198: *
199: 30 CONTINUE
200: LM1 = L - 1
201: DO 80 I = L, 1, -1
202: DO 40 J = 1, LM1
203: JP1 = J + 1
204: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
205: $ GO TO 50
206: 40 CONTINUE
207: J = L
208: GO TO 70
209: *
210: 50 CONTINUE
211: DO 60 J = JP1, L
212: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
213: $ GO TO 80
214: 60 CONTINUE
215: J = JP1 - 1
216: *
217: 70 CONTINUE
218: M = L
219: IFLOW = 1
220: GO TO 160
221: 80 CONTINUE
222: GO TO 100
223: *
224: * Find column with one nonzero in rows K through N
225: *
226: 90 CONTINUE
227: K = K + 1
228: *
229: 100 CONTINUE
230: DO 150 J = K, L
231: DO 110 I = K, LM1
232: IP1 = I + 1
233: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
234: $ GO TO 120
235: 110 CONTINUE
236: I = L
237: GO TO 140
238: 120 CONTINUE
239: DO 130 I = IP1, L
240: IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
241: $ GO TO 150
242: 130 CONTINUE
243: I = IP1 - 1
244: 140 CONTINUE
245: M = K
246: IFLOW = 2
247: GO TO 160
248: 150 CONTINUE
249: GO TO 190
250: *
251: * Permute rows M and I
252: *
253: 160 CONTINUE
254: LSCALE( M ) = I
255: IF( I.EQ.M )
256: $ GO TO 170
257: CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
258: CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
259: *
260: * Permute columns M and J
261: *
262: 170 CONTINUE
263: RSCALE( M ) = J
264: IF( J.EQ.M )
265: $ GO TO 180
266: CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
267: CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
268: *
269: 180 CONTINUE
270: GO TO ( 20, 90 )IFLOW
271: *
272: 190 CONTINUE
273: ILO = K
274: IHI = L
275: *
276: IF( LSAME( JOB, 'P' ) ) THEN
277: DO 195 I = ILO, IHI
278: LSCALE( I ) = ONE
279: RSCALE( I ) = ONE
280: 195 CONTINUE
281: RETURN
282: END IF
283: *
284: IF( ILO.EQ.IHI )
285: $ RETURN
286: *
287: * Balance the submatrix in rows ILO to IHI.
288: *
289: NR = IHI - ILO + 1
290: DO 200 I = ILO, IHI
291: RSCALE( I ) = ZERO
292: LSCALE( I ) = ZERO
293: *
294: WORK( I ) = ZERO
295: WORK( I+N ) = ZERO
296: WORK( I+2*N ) = ZERO
297: WORK( I+3*N ) = ZERO
298: WORK( I+4*N ) = ZERO
299: WORK( I+5*N ) = ZERO
300: 200 CONTINUE
301: *
302: * Compute right side vector in resulting linear equations
303: *
304: BASL = LOG10( SCLFAC )
305: DO 240 I = ILO, IHI
306: DO 230 J = ILO, IHI
307: TB = B( I, J )
308: TA = A( I, J )
309: IF( TA.EQ.ZERO )
310: $ GO TO 210
311: TA = LOG10( ABS( TA ) ) / BASL
312: 210 CONTINUE
313: IF( TB.EQ.ZERO )
314: $ GO TO 220
315: TB = LOG10( ABS( TB ) ) / BASL
316: 220 CONTINUE
317: WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
318: WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
319: 230 CONTINUE
320: 240 CONTINUE
321: *
322: COEF = ONE / DBLE( 2*NR )
323: COEF2 = COEF*COEF
324: COEF5 = HALF*COEF2
325: NRP2 = NR + 2
326: BETA = ZERO
327: IT = 1
328: *
329: * Start generalized conjugate gradient iteration
330: *
331: 250 CONTINUE
332: *
333: GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
334: $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
335: *
336: EW = ZERO
337: EWC = ZERO
338: DO 260 I = ILO, IHI
339: EW = EW + WORK( I+4*N )
340: EWC = EWC + WORK( I+5*N )
341: 260 CONTINUE
342: *
343: GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
344: IF( GAMMA.EQ.ZERO )
345: $ GO TO 350
346: IF( IT.NE.1 )
347: $ BETA = GAMMA / PGAMMA
348: T = COEF5*( EWC-THREE*EW )
349: TC = COEF5*( EW-THREE*EWC )
350: *
351: CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
352: CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
353: *
354: CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
355: CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
356: *
357: DO 270 I = ILO, IHI
358: WORK( I ) = WORK( I ) + TC
359: WORK( I+N ) = WORK( I+N ) + T
360: 270 CONTINUE
361: *
362: * Apply matrix to vector
363: *
364: DO 300 I = ILO, IHI
365: KOUNT = 0
366: SUM = ZERO
367: DO 290 J = ILO, IHI
368: IF( A( I, J ).EQ.ZERO )
369: $ GO TO 280
370: KOUNT = KOUNT + 1
371: SUM = SUM + WORK( J )
372: 280 CONTINUE
373: IF( B( I, J ).EQ.ZERO )
374: $ GO TO 290
375: KOUNT = KOUNT + 1
376: SUM = SUM + WORK( J )
377: 290 CONTINUE
378: WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
379: 300 CONTINUE
380: *
381: DO 330 J = ILO, IHI
382: KOUNT = 0
383: SUM = ZERO
384: DO 320 I = ILO, IHI
385: IF( A( I, J ).EQ.ZERO )
386: $ GO TO 310
387: KOUNT = KOUNT + 1
388: SUM = SUM + WORK( I+N )
389: 310 CONTINUE
390: IF( B( I, J ).EQ.ZERO )
391: $ GO TO 320
392: KOUNT = KOUNT + 1
393: SUM = SUM + WORK( I+N )
394: 320 CONTINUE
395: WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
396: 330 CONTINUE
397: *
398: SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
399: $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
400: ALPHA = GAMMA / SUM
401: *
402: * Determine correction to current iteration
403: *
404: CMAX = ZERO
405: DO 340 I = ILO, IHI
406: COR = ALPHA*WORK( I+N )
407: IF( ABS( COR ).GT.CMAX )
408: $ CMAX = ABS( COR )
409: LSCALE( I ) = LSCALE( I ) + COR
410: COR = ALPHA*WORK( I )
411: IF( ABS( COR ).GT.CMAX )
412: $ CMAX = ABS( COR )
413: RSCALE( I ) = RSCALE( I ) + COR
414: 340 CONTINUE
415: IF( CMAX.LT.HALF )
416: $ GO TO 350
417: *
418: CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
419: CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
420: *
421: PGAMMA = GAMMA
422: IT = IT + 1
423: IF( IT.LE.NRP2 )
424: $ GO TO 250
425: *
426: * End generalized conjugate gradient iteration
427: *
428: 350 CONTINUE
429: SFMIN = DLAMCH( 'S' )
430: SFMAX = ONE / SFMIN
431: LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
432: LSFMAX = INT( LOG10( SFMAX ) / BASL )
433: DO 360 I = ILO, IHI
434: IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
435: RAB = ABS( A( I, IRAB+ILO-1 ) )
436: IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
437: RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
438: LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
439: IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
440: IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
441: LSCALE( I ) = SCLFAC**IR
442: ICAB = IDAMAX( IHI, A( 1, I ), 1 )
443: CAB = ABS( A( ICAB, I ) )
444: ICAB = IDAMAX( IHI, B( 1, I ), 1 )
445: CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
446: LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
447: JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
448: JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
449: RSCALE( I ) = SCLFAC**JC
450: 360 CONTINUE
451: *
452: * Row scaling of matrices A and B
453: *
454: DO 370 I = ILO, IHI
455: CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
456: CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
457: 370 CONTINUE
458: *
459: * Column scaling of matrices A and B
460: *
461: DO 380 J = ILO, IHI
462: CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
463: CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
464: 380 CONTINUE
465: *
466: RETURN
467: *
468: * End of DGGBAL
469: *
470: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggbal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack