Annotation of rpl/lapack/lapack/zggbal.f, revision 1.7
1.1 bertrand 1: SUBROUTINE ZGGBAL( 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 LSCALE( * ), RSCALE( * ), WORK( * )
15: COMPLEX*16 A( LDA, * ), B( LDB, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * ZGGBAL balances a pair of general complex 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) COMPLEX*16 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) COMPLEX*16 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) is the scaling factor
73: * 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) is the scaling
84: * factor applied to column j, then
85: * RSCALE(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: COMPLEX*16 CZERO
113: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
114: * ..
115: * .. Local Scalars ..
116: INTEGER I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
117: $ K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
118: $ M, NR, NRP2
119: DOUBLE PRECISION ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
120: $ COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
121: $ SFMIN, SUM, T, TA, TB, TC
122: COMPLEX*16 CDUM
123: * ..
124: * .. External Functions ..
125: LOGICAL LSAME
126: INTEGER IZAMAX
127: DOUBLE PRECISION DDOT, DLAMCH
128: EXTERNAL LSAME, IZAMAX, DDOT, DLAMCH
129: * ..
130: * .. External Subroutines ..
131: EXTERNAL DAXPY, DSCAL, XERBLA, ZDSCAL, ZSWAP
132: * ..
133: * .. Intrinsic Functions ..
134: INTRINSIC ABS, DBLE, DIMAG, INT, LOG10, MAX, MIN, SIGN
135: * ..
136: * .. Statement Functions ..
137: DOUBLE PRECISION CABS1
138: * ..
139: * .. Statement Function definitions ..
140: CABS1( CDUM ) = ABS( DBLE( CDUM ) ) + ABS( DIMAG( CDUM ) )
141: * ..
142: * .. Executable Statements ..
143: *
144: * Test the input parameters
145: *
146: INFO = 0
147: IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
148: $ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
149: INFO = -1
150: ELSE IF( N.LT.0 ) THEN
151: INFO = -2
152: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
153: INFO = -4
154: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
155: INFO = -6
156: END IF
157: IF( INFO.NE.0 ) THEN
158: CALL XERBLA( 'ZGGBAL', -INFO )
159: RETURN
160: END IF
161: *
162: * Quick return if possible
163: *
164: IF( N.EQ.0 ) THEN
165: ILO = 1
166: IHI = N
167: RETURN
168: END IF
169: *
170: IF( N.EQ.1 ) THEN
171: ILO = 1
172: IHI = N
173: LSCALE( 1 ) = ONE
174: RSCALE( 1 ) = ONE
175: RETURN
176: END IF
177: *
178: IF( LSAME( JOB, 'N' ) ) THEN
179: ILO = 1
180: IHI = N
181: DO 10 I = 1, N
182: LSCALE( I ) = ONE
183: RSCALE( I ) = ONE
184: 10 CONTINUE
185: RETURN
186: END IF
187: *
188: K = 1
189: L = N
190: IF( LSAME( JOB, 'S' ) )
191: $ GO TO 190
192: *
193: GO TO 30
194: *
195: * Permute the matrices A and B to isolate the eigenvalues.
196: *
197: * Find row with one nonzero in columns 1 through L
198: *
199: 20 CONTINUE
200: L = LM1
201: IF( L.NE.1 )
202: $ GO TO 30
203: *
204: RSCALE( 1 ) = 1
205: LSCALE( 1 ) = 1
206: GO TO 190
207: *
208: 30 CONTINUE
209: LM1 = L - 1
210: DO 80 I = L, 1, -1
211: DO 40 J = 1, LM1
212: JP1 = J + 1
213: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
214: $ GO TO 50
215: 40 CONTINUE
216: J = L
217: GO TO 70
218: *
219: 50 CONTINUE
220: DO 60 J = JP1, L
221: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
222: $ GO TO 80
223: 60 CONTINUE
224: J = JP1 - 1
225: *
226: 70 CONTINUE
227: M = L
228: IFLOW = 1
229: GO TO 160
230: 80 CONTINUE
231: GO TO 100
232: *
233: * Find column with one nonzero in rows K through N
234: *
235: 90 CONTINUE
236: K = K + 1
237: *
238: 100 CONTINUE
239: DO 150 J = K, L
240: DO 110 I = K, LM1
241: IP1 = I + 1
242: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
243: $ GO TO 120
244: 110 CONTINUE
245: I = L
246: GO TO 140
247: 120 CONTINUE
248: DO 130 I = IP1, L
249: IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
250: $ GO TO 150
251: 130 CONTINUE
252: I = IP1 - 1
253: 140 CONTINUE
254: M = K
255: IFLOW = 2
256: GO TO 160
257: 150 CONTINUE
258: GO TO 190
259: *
260: * Permute rows M and I
261: *
262: 160 CONTINUE
263: LSCALE( M ) = I
264: IF( I.EQ.M )
265: $ GO TO 170
266: CALL ZSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
267: CALL ZSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
268: *
269: * Permute columns M and J
270: *
271: 170 CONTINUE
272: RSCALE( M ) = J
273: IF( J.EQ.M )
274: $ GO TO 180
275: CALL ZSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
276: CALL ZSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
277: *
278: 180 CONTINUE
279: GO TO ( 20, 90 )IFLOW
280: *
281: 190 CONTINUE
282: ILO = K
283: IHI = L
284: *
285: IF( LSAME( JOB, 'P' ) ) THEN
286: DO 195 I = ILO, IHI
287: LSCALE( I ) = ONE
288: RSCALE( I ) = ONE
289: 195 CONTINUE
290: RETURN
291: END IF
292: *
293: IF( ILO.EQ.IHI )
294: $ RETURN
295: *
296: * Balance the submatrix in rows ILO to IHI.
297: *
298: NR = IHI - ILO + 1
299: DO 200 I = ILO, IHI
300: RSCALE( I ) = ZERO
301: LSCALE( I ) = ZERO
302: *
303: WORK( I ) = ZERO
304: WORK( I+N ) = ZERO
305: WORK( I+2*N ) = ZERO
306: WORK( I+3*N ) = ZERO
307: WORK( I+4*N ) = ZERO
308: WORK( I+5*N ) = ZERO
309: 200 CONTINUE
310: *
311: * Compute right side vector in resulting linear equations
312: *
313: BASL = LOG10( SCLFAC )
314: DO 240 I = ILO, IHI
315: DO 230 J = ILO, IHI
316: IF( A( I, J ).EQ.CZERO ) THEN
317: TA = ZERO
318: GO TO 210
319: END IF
320: TA = LOG10( CABS1( A( I, J ) ) ) / BASL
321: *
322: 210 CONTINUE
323: IF( B( I, J ).EQ.CZERO ) THEN
324: TB = ZERO
325: GO TO 220
326: END IF
327: TB = LOG10( CABS1( B( I, J ) ) ) / BASL
328: *
329: 220 CONTINUE
330: WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
331: WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
332: 230 CONTINUE
333: 240 CONTINUE
334: *
335: COEF = ONE / DBLE( 2*NR )
336: COEF2 = COEF*COEF
337: COEF5 = HALF*COEF2
338: NRP2 = NR + 2
339: BETA = ZERO
340: IT = 1
341: *
342: * Start generalized conjugate gradient iteration
343: *
344: 250 CONTINUE
345: *
346: GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
347: $ DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
348: *
349: EW = ZERO
350: EWC = ZERO
351: DO 260 I = ILO, IHI
352: EW = EW + WORK( I+4*N )
353: EWC = EWC + WORK( I+5*N )
354: 260 CONTINUE
355: *
356: GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
357: IF( GAMMA.EQ.ZERO )
358: $ GO TO 350
359: IF( IT.NE.1 )
360: $ BETA = GAMMA / PGAMMA
361: T = COEF5*( EWC-THREE*EW )
362: TC = COEF5*( EW-THREE*EWC )
363: *
364: CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
365: CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
366: *
367: CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
368: CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
369: *
370: DO 270 I = ILO, IHI
371: WORK( I ) = WORK( I ) + TC
372: WORK( I+N ) = WORK( I+N ) + T
373: 270 CONTINUE
374: *
375: * Apply matrix to vector
376: *
377: DO 300 I = ILO, IHI
378: KOUNT = 0
379: SUM = ZERO
380: DO 290 J = ILO, IHI
381: IF( A( I, J ).EQ.CZERO )
382: $ GO TO 280
383: KOUNT = KOUNT + 1
384: SUM = SUM + WORK( J )
385: 280 CONTINUE
386: IF( B( I, J ).EQ.CZERO )
387: $ GO TO 290
388: KOUNT = KOUNT + 1
389: SUM = SUM + WORK( J )
390: 290 CONTINUE
391: WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
392: 300 CONTINUE
393: *
394: DO 330 J = ILO, IHI
395: KOUNT = 0
396: SUM = ZERO
397: DO 320 I = ILO, IHI
398: IF( A( I, J ).EQ.CZERO )
399: $ GO TO 310
400: KOUNT = KOUNT + 1
401: SUM = SUM + WORK( I+N )
402: 310 CONTINUE
403: IF( B( I, J ).EQ.CZERO )
404: $ GO TO 320
405: KOUNT = KOUNT + 1
406: SUM = SUM + WORK( I+N )
407: 320 CONTINUE
408: WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
409: 330 CONTINUE
410: *
411: SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
412: $ DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
413: ALPHA = GAMMA / SUM
414: *
415: * Determine correction to current iteration
416: *
417: CMAX = ZERO
418: DO 340 I = ILO, IHI
419: COR = ALPHA*WORK( I+N )
420: IF( ABS( COR ).GT.CMAX )
421: $ CMAX = ABS( COR )
422: LSCALE( I ) = LSCALE( I ) + COR
423: COR = ALPHA*WORK( I )
424: IF( ABS( COR ).GT.CMAX )
425: $ CMAX = ABS( COR )
426: RSCALE( I ) = RSCALE( I ) + COR
427: 340 CONTINUE
428: IF( CMAX.LT.HALF )
429: $ GO TO 350
430: *
431: CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
432: CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
433: *
434: PGAMMA = GAMMA
435: IT = IT + 1
436: IF( IT.LE.NRP2 )
437: $ GO TO 250
438: *
439: * End generalized conjugate gradient iteration
440: *
441: 350 CONTINUE
442: SFMIN = DLAMCH( 'S' )
443: SFMAX = ONE / SFMIN
444: LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
445: LSFMAX = INT( LOG10( SFMAX ) / BASL )
446: DO 360 I = ILO, IHI
447: IRAB = IZAMAX( N-ILO+1, A( I, ILO ), LDA )
448: RAB = ABS( A( I, IRAB+ILO-1 ) )
449: IRAB = IZAMAX( N-ILO+1, B( I, ILO ), LDB )
450: RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
451: LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
452: IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
453: IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
454: LSCALE( I ) = SCLFAC**IR
455: ICAB = IZAMAX( IHI, A( 1, I ), 1 )
456: CAB = ABS( A( ICAB, I ) )
457: ICAB = IZAMAX( IHI, B( 1, I ), 1 )
458: CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
459: LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
460: JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
461: JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
462: RSCALE( I ) = SCLFAC**JC
463: 360 CONTINUE
464: *
465: * Row scaling of matrices A and B
466: *
467: DO 370 I = ILO, IHI
468: CALL ZDSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
469: CALL ZDSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
470: 370 CONTINUE
471: *
472: * Column scaling of matrices A and B
473: *
474: DO 380 J = ILO, IHI
475: CALL ZDSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
476: CALL ZDSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
477: 380 CONTINUE
478: *
479: RETURN
480: *
481: * End of ZGGBAL
482: *
483: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zggbal.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack