1      SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
2     $                   RSCALE, WORK, INFO )
3*
4*  -- LAPACK routine (version 3.0) --
5*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
6*     Courant Institute, Argonne National Lab, and Rice University
7*     April 12, 2001
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 (6*N)
92*
93*  INFO    (output) INTEGER
94*          = 0:  successful exit
95*          < 0:  if INFO = -i, the i-th argument had an illegal value.
96*
97*  Further Details
98*  ===============
99*
100*  See R.C. WARD, Balancing the generalized eigenvalue problem,
101*                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
102*
103*  =====================================================================
104*
105*     .. Parameters ..
106      DOUBLE PRECISION   ZERO, HALF, ONE
107      PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
108      DOUBLE PRECISION   THREE, SCLFAC
109      PARAMETER          ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
110*     ..
111*     .. Local Scalars ..
112      INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
113     $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
114     $                   M, NR, NRP2
115      DOUBLE PRECISION   ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
116     $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
117     $                   SFMIN, SUM, T, TA, TB, TC
118*     ..
119*     .. External Functions ..
120      LOGICAL            LSAME
121      INTEGER            IDAMAX
122      DOUBLE PRECISION   DDOT, DLAMCH
123      EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
124*     ..
125*     .. External Subroutines ..
126      EXTERNAL           DAXPY, DSCAL, DSWAP, XERBLA
127*     ..
128*     .. Intrinsic Functions ..
129      INTRINSIC          ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
130*     ..
131*     .. Executable Statements ..
132*
133*     Test the input parameters
134*
135      INFO = 0
136      IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
137     $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
138         INFO = -1
139      ELSE IF( N.LT.0 ) THEN
140         INFO = -2
141      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
142         INFO = -4
143      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
144         INFO = -6
145      END IF
146      IF( INFO.NE.0 ) THEN
147         CALL XERBLA( 'DGGBAL', -INFO )
148         RETURN
149      END IF
150*
151      K = 1
152      L = N
153*
154*     Quick return if possible
155*
156      IF( N.EQ.0 )
157     $   RETURN
158*
159      IF( LSAME( JOB, 'N' ) ) THEN
160         ILO = 1
161         IHI = N
162         DO 10 I = 1, N
163            LSCALE( I ) = ONE
164            RSCALE( I ) = ONE
165   10    CONTINUE
166         RETURN
167      END IF
168*
169      IF( K.EQ.L ) THEN
170         ILO = 1
171         IHI = 1
172         LSCALE( 1 ) = ONE
173         RSCALE( 1 ) = ONE
174         RETURN
175      END IF
176*
177      IF( LSAME( JOB, 'S' ) )
178     $   GO TO 190
179*
180      GO TO 30
181*
182*     Permute the matrices A and B to isolate the eigenvalues.
183*
184*     Find row with one nonzero in columns 1 through L
185*
186   20 CONTINUE
187      L = LM1
188      IF( L.NE.1 )
189     $   GO TO 30
190*
191      RSCALE( 1 ) = ONE
192      LSCALE( 1 ) = ONE
193      GO TO 190
194*
195   30 CONTINUE
196      LM1 = L - 1
197      DO 80 I = L, 1, -1
198         DO 40 J = 1, LM1
199            JP1 = J + 1
200            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
201     $         GO TO 50
202   40    CONTINUE
203         J = L
204         GO TO 70
205*
206   50    CONTINUE
207         DO 60 J = JP1, L
208            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
209     $         GO TO 80
210   60    CONTINUE
211         J = JP1 - 1
212*
213   70    CONTINUE
214         M = L
215         IFLOW = 1
216         GO TO 160
217   80 CONTINUE
218      GO TO 100
219*
220*     Find column with one nonzero in rows K through N
221*
222   90 CONTINUE
223      K = K + 1
224*
225  100 CONTINUE
226      DO 150 J = K, L
227         DO 110 I = K, LM1
228            IP1 = I + 1
229            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
230     $         GO TO 120
231  110    CONTINUE
232         I = L
233         GO TO 140
234  120    CONTINUE
235         DO 130 I = IP1, L
236            IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
237     $         GO TO 150
238  130    CONTINUE
239         I = IP1 - 1
240  140    CONTINUE
241         M = K
242         IFLOW = 2
243         GO TO 160
244  150 CONTINUE
245      GO TO 190
246*
247*     Permute rows M and I
248*
249  160 CONTINUE
250      LSCALE( M ) = DBLE( I )
251      IF( I.EQ.M )
252     $   GO TO 170
253      CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
254      CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
255*
256*     Permute columns M and J
257*
258  170 CONTINUE
259      RSCALE( M ) = DBLE( J )
260      IF( J.EQ.M )
261     $   GO TO 180
262      CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
263      CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
264*
265  180 CONTINUE
266      GO TO ( 20, 90 )IFLOW
267*
268  190 CONTINUE
269      ILO = K
270      IHI = L
271*
272      IF( ILO.EQ.IHI )
273     $   RETURN
274*
275      IF( LSAME( JOB, 'P' ) )
276     $   RETURN
277*
278*     Balance the submatrix in rows ILO to IHI.
279*
280      NR = IHI - ILO + 1
281      DO 200 I = ILO, IHI
282         RSCALE( I ) = ZERO
283         LSCALE( I ) = ZERO
284*
285         WORK( I ) = ZERO
286         WORK( I+N ) = ZERO
287         WORK( I+2*N ) = ZERO
288         WORK( I+3*N ) = ZERO
289         WORK( I+4*N ) = ZERO
290         WORK( I+5*N ) = ZERO
291  200 CONTINUE
292*
293*     Compute right side vector in resulting linear equations
294*
295      BASL = LOG10( SCLFAC )
296      DO 240 I = ILO, IHI
297         DO 230 J = ILO, IHI
298            TB = B( I, J )
299            TA = A( I, J )
300            IF( TA.EQ.ZERO )
301     $         GO TO 210
302            TA = LOG10( ABS( TA ) ) / BASL
303  210       CONTINUE
304            IF( TB.EQ.ZERO )
305     $         GO TO 220
306            TB = LOG10( ABS( TB ) ) / BASL
307  220       CONTINUE
308            WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
309            WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
310  230    CONTINUE
311  240 CONTINUE
312*
313      COEF = ONE / DBLE( 2*NR )
314      COEF2 = COEF*COEF
315      COEF5 = HALF*COEF2
316      NRP2 = NR + 2
317      BETA = ZERO
318      IT = 1
319*
320*     Start generalized conjugate gradient iteration
321*
322  250 CONTINUE
323*
324      GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
325     $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
326*
327      EW = ZERO
328      EWC = ZERO
329      DO 260 I = ILO, IHI
330         EW = EW + WORK( I+4*N )
331         EWC = EWC + WORK( I+5*N )
332  260 CONTINUE
333*
334      GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
335      IF( GAMMA.EQ.ZERO )
336     $   GO TO 350
337      IF( IT.NE.1 )
338     $   BETA = GAMMA / PGAMMA
339      T = COEF5*( EWC-THREE*EW )
340      TC = COEF5*( EW-THREE*EWC )
341*
342      CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
343      CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
344*
345      CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
346      CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
347*
348      DO 270 I = ILO, IHI
349         WORK( I ) = WORK( I ) + TC
350         WORK( I+N ) = WORK( I+N ) + T
351  270 CONTINUE
352*
353*     Apply matrix to vector
354*
355      DO 300 I = ILO, IHI
356         KOUNT = 0
357         SUM = ZERO
358         DO 290 J = ILO, IHI
359            IF( A( I, J ).EQ.ZERO )
360     $         GO TO 280
361            KOUNT = KOUNT + 1
362            SUM = SUM + WORK( J )
363  280       CONTINUE
364            IF( B( I, J ).EQ.ZERO )
365     $         GO TO 290
366            KOUNT = KOUNT + 1
367            SUM = SUM + WORK( J )
368  290    CONTINUE
369         WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
370  300 CONTINUE
371*
372      DO 330 J = ILO, IHI
373         KOUNT = 0
374         SUM = ZERO
375         DO 320 I = ILO, IHI
376            IF( A( I, J ).EQ.ZERO )
377     $         GO TO 310
378            KOUNT = KOUNT + 1
379            SUM = SUM + WORK( I+N )
380  310       CONTINUE
381            IF( B( I, J ).EQ.ZERO )
382     $         GO TO 320
383            KOUNT = KOUNT + 1
384            SUM = SUM + WORK( I+N )
385  320    CONTINUE
386         WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
387  330 CONTINUE
388*
389      SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
390     $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
391      ALPHA = GAMMA / SUM
392*
393*     Determine correction to current iteration
394*
395      CMAX = ZERO
396      DO 340 I = ILO, IHI
397         COR = ALPHA*WORK( I+N )
398         IF( ABS( COR ).GT.CMAX )
399     $      CMAX = ABS( COR )
400         LSCALE( I ) = LSCALE( I ) + COR
401         COR = ALPHA*WORK( I )
402         IF( ABS( COR ).GT.CMAX )
403     $      CMAX = ABS( COR )
404         RSCALE( I ) = RSCALE( I ) + COR
405  340 CONTINUE
406      IF( CMAX.LT.HALF )
407     $   GO TO 350
408*
409      CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
410      CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
411*
412      PGAMMA = GAMMA
413      IT = IT + 1
414      IF( IT.LE.NRP2 )
415     $   GO TO 250
416*
417*     End generalized conjugate gradient iteration
418*
419  350 CONTINUE
420      SFMIN = DLAMCH( 'S' )
421      SFMAX = ONE / SFMIN
422      LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
423      LSFMAX = INT( LOG10( SFMAX ) / BASL )
424      DO 360 I = ILO, IHI
425         IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
426         RAB = ABS( A( I, IRAB+ILO-1 ) )
427         IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
428         RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
429         LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
430         IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
431         IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
432         LSCALE( I ) = SCLFAC**IR
433         ICAB = IDAMAX( IHI, A( 1, I ), 1 )
434         CAB = ABS( A( ICAB, I ) )
435         ICAB = IDAMAX( IHI, B( 1, I ), 1 )
436         CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
437         LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
438         JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
439         JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
440         RSCALE( I ) = SCLFAC**JC
441  360 CONTINUE
442*
443*     Row scaling of matrices A and B
444*
445      DO 370 I = ILO, IHI
446         CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
447         CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
448  370 CONTINUE
449*
450*     Column scaling of matrices A and B
451*
452      DO 380 J = ILO, IHI
453         CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
454         CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
455  380 CONTINUE
456*
457      RETURN
458*
459*     End of DGGBAL
460*
461      END
462