1 SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
2 $ RSCALE, WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2.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 * June 2010
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) DOUBLE PRECISION 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
2 $ RSCALE, WORK, INFO )
3 *
4 * -- LAPACK routine (version 3.2.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 * June 2010
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) DOUBLE PRECISION 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