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