1*> \brief \b CGGBAL
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CGGBAL + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggbal.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggbal.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggbal.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGGBAL( 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*       REAL               LSCALE( * ), RSCALE( * ), WORK( * )
30*       COMPLEX            A( LDA, * ), B( LDB, * )
31*       ..
32*
33*
34*> \par Purpose:
35*  =============
36*>
37*> \verbatim
38*>
39*> CGGBAL balances a pair of general complex 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 COMPLEX 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 COMPLEX 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 REAL 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) is the scaling factor
119*>          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 REAL 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) is the scaling
133*>          factor applied to column j, then
134*>            RSCALE(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 REAL 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 complexGBcomputational
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 CGGBAL( 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      REAL               LSCALE( * ), RSCALE( * ), WORK( * )
188      COMPLEX            A( LDA, * ), B( LDB, * )
189*     ..
190*
191*  =====================================================================
192*
193*     .. Parameters ..
194      REAL               ZERO, HALF, ONE
195      PARAMETER          ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 )
196      REAL               THREE, SCLFAC
197      PARAMETER          ( THREE = 3.0E+0, SCLFAC = 1.0E+1 )
198      COMPLEX            CZERO
199      PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
200*     ..
201*     .. Local Scalars ..
202      INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
203     $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
204     $                   M, NR, NRP2
205      REAL               ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
206     $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
207     $                   SFMIN, SUM, T, TA, TB, TC
208      COMPLEX            CDUM
209*     ..
210*     .. External Functions ..
211      LOGICAL            LSAME
212      INTEGER            ICAMAX
213      REAL               SDOT, SLAMCH
214      EXTERNAL           LSAME, ICAMAX, SDOT, SLAMCH
215*     ..
216*     .. External Subroutines ..
217      EXTERNAL           CSSCAL, CSWAP, SAXPY, SSCAL, XERBLA
218*     ..
219*     .. Intrinsic Functions ..
220      INTRINSIC          ABS, AIMAG, INT, LOG10, MAX, MIN, REAL, SIGN
221*     ..
222*     .. Statement Functions ..
223      REAL               CABS1
224*     ..
225*     .. Statement Function definitions ..
226      CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
227*     ..
228*     .. Executable Statements ..
229*
230*     Test the input parameters
231*
232      INFO = 0
233      IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
234     $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
235         INFO = -1
236      ELSE IF( N.LT.0 ) THEN
237         INFO = -2
238      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
239         INFO = -4
240      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
241         INFO = -6
242      END IF
243      IF( INFO.NE.0 ) THEN
244         CALL XERBLA( 'CGGBAL', -INFO )
245         RETURN
246      END IF
247*
248*     Quick return if possible
249*
250      IF( N.EQ.0 ) THEN
251         ILO = 1
252         IHI = N
253         RETURN
254      END IF
255*
256      IF( N.EQ.1 ) THEN
257         ILO = 1
258         IHI = N
259         LSCALE( 1 ) = ONE
260         RSCALE( 1 ) = ONE
261         RETURN
262      END IF
263*
264      IF( LSAME( JOB, 'N' ) ) THEN
265         ILO = 1
266         IHI = N
267         DO 10 I = 1, N
268            LSCALE( I ) = ONE
269            RSCALE( I ) = ONE
270   10    CONTINUE
271         RETURN
272      END IF
273*
274      K = 1
275      L = N
276      IF( LSAME( JOB, 'S' ) )
277     $   GO TO 190
278*
279      GO TO 30
280*
281*     Permute the matrices A and B to isolate the eigenvalues.
282*
283*     Find row with one nonzero in columns 1 through L
284*
285   20 CONTINUE
286      L = LM1
287      IF( L.NE.1 )
288     $   GO TO 30
289*
290      RSCALE( 1 ) = ONE
291      LSCALE( 1 ) = ONE
292      GO TO 190
293*
294   30 CONTINUE
295      LM1 = L - 1
296      DO 80 I = L, 1, -1
297         DO 40 J = 1, LM1
298            JP1 = J + 1
299            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
300     $         GO TO 50
301   40    CONTINUE
302         J = L
303         GO TO 70
304*
305   50    CONTINUE
306         DO 60 J = JP1, L
307            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
308     $         GO TO 80
309   60    CONTINUE
310         J = JP1 - 1
311*
312   70    CONTINUE
313         M = L
314         IFLOW = 1
315         GO TO 160
316   80 CONTINUE
317      GO TO 100
318*
319*     Find column with one nonzero in rows K through N
320*
321   90 CONTINUE
322      K = K + 1
323*
324  100 CONTINUE
325      DO 150 J = K, L
326         DO 110 I = K, LM1
327            IP1 = I + 1
328            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
329     $         GO TO 120
330  110    CONTINUE
331         I = L
332         GO TO 140
333  120    CONTINUE
334         DO 130 I = IP1, L
335            IF( A( I, J ).NE.CZERO .OR. B( I, J ).NE.CZERO )
336     $         GO TO 150
337  130    CONTINUE
338         I = IP1 - 1
339  140    CONTINUE
340         M = K
341         IFLOW = 2
342         GO TO 160
343  150 CONTINUE
344      GO TO 190
345*
346*     Permute rows M and I
347*
348  160 CONTINUE
349      LSCALE( M ) = I
350      IF( I.EQ.M )
351     $   GO TO 170
352      CALL CSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
353      CALL CSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
354*
355*     Permute columns M and J
356*
357  170 CONTINUE
358      RSCALE( M ) = J
359      IF( J.EQ.M )
360     $   GO TO 180
361      CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
362      CALL CSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
363*
364  180 CONTINUE
365      GO TO ( 20, 90 )IFLOW
366*
367  190 CONTINUE
368      ILO = K
369      IHI = L
370*
371      IF( LSAME( JOB, 'P' ) ) THEN
372         DO 195 I = ILO, IHI
373            LSCALE( I ) = ONE
374            RSCALE( I ) = ONE
375  195    CONTINUE
376         RETURN
377      END IF
378*
379      IF( ILO.EQ.IHI )
380     $   RETURN
381*
382*     Balance the submatrix in rows ILO to IHI.
383*
384      NR = IHI - ILO + 1
385      DO 200 I = ILO, IHI
386         RSCALE( I ) = ZERO
387         LSCALE( I ) = ZERO
388*
389         WORK( I ) = ZERO
390         WORK( I+N ) = ZERO
391         WORK( I+2*N ) = ZERO
392         WORK( I+3*N ) = ZERO
393         WORK( I+4*N ) = ZERO
394         WORK( I+5*N ) = ZERO
395  200 CONTINUE
396*
397*     Compute right side vector in resulting linear equations
398*
399      BASL = LOG10( SCLFAC )
400      DO 240 I = ILO, IHI
401         DO 230 J = ILO, IHI
402            IF( A( I, J ).EQ.CZERO ) THEN
403               TA = ZERO
404               GO TO 210
405            END IF
406            TA = LOG10( CABS1( A( I, J ) ) ) / BASL
407*
408  210       CONTINUE
409            IF( B( I, J ).EQ.CZERO ) THEN
410               TB = ZERO
411               GO TO 220
412            END IF
413            TB = LOG10( CABS1( B( I, J ) ) ) / BASL
414*
415  220       CONTINUE
416            WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
417            WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
418  230    CONTINUE
419  240 CONTINUE
420*
421      COEF = ONE / REAL( 2*NR )
422      COEF2 = COEF*COEF
423      COEF5 = HALF*COEF2
424      NRP2 = NR + 2
425      BETA = ZERO
426      IT = 1
427*
428*     Start generalized conjugate gradient iteration
429*
430  250 CONTINUE
431*
432      GAMMA = SDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
433     $        SDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
434*
435      EW = ZERO
436      EWC = ZERO
437      DO 260 I = ILO, IHI
438         EW = EW + WORK( I+4*N )
439         EWC = EWC + WORK( I+5*N )
440  260 CONTINUE
441*
442      GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
443      IF( GAMMA.EQ.ZERO )
444     $   GO TO 350
445      IF( IT.NE.1 )
446     $   BETA = GAMMA / PGAMMA
447      T = COEF5*( EWC-THREE*EW )
448      TC = COEF5*( EW-THREE*EWC )
449*
450      CALL SSCAL( NR, BETA, WORK( ILO ), 1 )
451      CALL SSCAL( NR, BETA, WORK( ILO+N ), 1 )
452*
453      CALL SAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
454      CALL SAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
455*
456      DO 270 I = ILO, IHI
457         WORK( I ) = WORK( I ) + TC
458         WORK( I+N ) = WORK( I+N ) + T
459  270 CONTINUE
460*
461*     Apply matrix to vector
462*
463      DO 300 I = ILO, IHI
464         KOUNT = 0
465         SUM = ZERO
466         DO 290 J = ILO, IHI
467            IF( A( I, J ).EQ.CZERO )
468     $         GO TO 280
469            KOUNT = KOUNT + 1
470            SUM = SUM + WORK( J )
471  280       CONTINUE
472            IF( B( I, J ).EQ.CZERO )
473     $         GO TO 290
474            KOUNT = KOUNT + 1
475            SUM = SUM + WORK( J )
476  290    CONTINUE
477         WORK( I+2*N ) = REAL( KOUNT )*WORK( I+N ) + SUM
478  300 CONTINUE
479*
480      DO 330 J = ILO, IHI
481         KOUNT = 0
482         SUM = ZERO
483         DO 320 I = ILO, IHI
484            IF( A( I, J ).EQ.CZERO )
485     $         GO TO 310
486            KOUNT = KOUNT + 1
487            SUM = SUM + WORK( I+N )
488  310       CONTINUE
489            IF( B( I, J ).EQ.CZERO )
490     $         GO TO 320
491            KOUNT = KOUNT + 1
492            SUM = SUM + WORK( I+N )
493  320    CONTINUE
494         WORK( J+3*N ) = REAL( KOUNT )*WORK( J ) + SUM
495  330 CONTINUE
496*
497      SUM = SDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
498     $      SDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
499      ALPHA = GAMMA / SUM
500*
501*     Determine correction to current iteration
502*
503      CMAX = ZERO
504      DO 340 I = ILO, IHI
505         COR = ALPHA*WORK( I+N )
506         IF( ABS( COR ).GT.CMAX )
507     $      CMAX = ABS( COR )
508         LSCALE( I ) = LSCALE( I ) + COR
509         COR = ALPHA*WORK( I )
510         IF( ABS( COR ).GT.CMAX )
511     $      CMAX = ABS( COR )
512         RSCALE( I ) = RSCALE( I ) + COR
513  340 CONTINUE
514      IF( CMAX.LT.HALF )
515     $   GO TO 350
516*
517      CALL SAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
518      CALL SAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
519*
520      PGAMMA = GAMMA
521      IT = IT + 1
522      IF( IT.LE.NRP2 )
523     $   GO TO 250
524*
525*     End generalized conjugate gradient iteration
526*
527  350 CONTINUE
528      SFMIN = SLAMCH( 'S' )
529      SFMAX = ONE / SFMIN
530      LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
531      LSFMAX = INT( LOG10( SFMAX ) / BASL )
532      DO 360 I = ILO, IHI
533         IRAB = ICAMAX( N-ILO+1, A( I, ILO ), LDA )
534         RAB = ABS( A( I, IRAB+ILO-1 ) )
535         IRAB = ICAMAX( N-ILO+1, B( I, ILO ), LDB )
536         RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
537         LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
538         IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
539         IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
540         LSCALE( I ) = SCLFAC**IR
541         ICAB = ICAMAX( IHI, A( 1, I ), 1 )
542         CAB = ABS( A( ICAB, I ) )
543         ICAB = ICAMAX( IHI, B( 1, I ), 1 )
544         CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
545         LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
546         JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
547         JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
548         RSCALE( I ) = SCLFAC**JC
549  360 CONTINUE
550*
551*     Row scaling of matrices A and B
552*
553      DO 370 I = ILO, IHI
554         CALL CSSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
555         CALL CSSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
556  370 CONTINUE
557*
558*     Column scaling of matrices A and B
559*
560      DO 380 J = ILO, IHI
561         CALL CSSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
562         CALL CSSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
563  380 CONTINUE
564*
565      RETURN
566*
567*     End of CGGBAL
568*
569      END
570