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
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