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