LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zggbal()

subroutine zggbal ( character  JOB,
integer  N,
complex*16, dimension( lda, * )  A,
integer  LDA,
complex*16, dimension( ldb, * )  B,
integer  LDB,
integer  ILO,
integer  IHI,
double precision, dimension( * )  LSCALE,
double precision, dimension( * )  RSCALE,
double precision, dimension( * )  WORK,
integer  INFO 
)

ZGGBAL

Download ZGGBAL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 ZGGBAL balances a pair of general complex matrices (A,B).  This
 involves, first, permuting A and B by similarity transformations to
 isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
 elements on the diagonal; and second, applying a diagonal similarity
 transformation to rows and columns ILO to IHI to make the rows
 and columns as close in norm as possible. Both steps are optional.

 Balancing may reduce the 1-norm of the matrices, and improve the
 accuracy of the computed eigenvalues and/or eigenvectors in the
 generalized eigenvalue problem A*x = lambda*B*x.
Parameters
[in]JOB
          JOB is CHARACTER*1
          Specifies the operations to be performed on A and B:
          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
                  and RSCALE(I) = 1.0 for i=1,...,N;
          = 'P':  permute only;
          = 'S':  scale only;
          = 'B':  both permute and scale.
[in]N
          N is INTEGER
          The order of the matrices A and B.  N >= 0.
[in,out]A
          A is COMPLEX*16 array, dimension (LDA,N)
          On entry, the input matrix A.
          On exit, A is overwritten by the balanced matrix.
          If JOB = 'N', A is not referenced.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,N).
[in,out]B
          B is COMPLEX*16 array, dimension (LDB,N)
          On entry, the input matrix B.
          On exit, B is overwritten by the balanced matrix.
          If JOB = 'N', B is not referenced.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[out]ILO
          ILO is INTEGER
[out]IHI
          IHI is INTEGER
          ILO and IHI are set to integers such that on exit
          A(i,j) = 0 and B(i,j) = 0 if i > j and
          j = 1,...,ILO-1 or i = IHI+1,...,N.
          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
[out]LSCALE
          LSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the left side of A and B.  If P(j) is the index of the
          row interchanged with row j, and D(j) is the scaling factor
          applied to row j, then
            LSCALE(j) = P(j)    for J = 1,...,ILO-1
                      = D(j)    for J = ILO,...,IHI
                      = P(j)    for J = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]RSCALE
          RSCALE is DOUBLE PRECISION array, dimension (N)
          Details of the permutations and scaling factors applied
          to the right side of A and B.  If P(j) is the index of the
          column interchanged with column j, and D(j) is the scaling
          factor applied to column j, then
            RSCALE(j) = P(j)    for J = 1,...,ILO-1
                      = D(j)    for J = ILO,...,IHI
                      = P(j)    for J = IHI+1,...,N.
          The order in which the interchanges are made is N to IHI+1,
          then 1 to ILO-1.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (lwork)
          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
          at least 1 when JOB = 'N' or 'P'.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
  See R.C. WARD, Balancing the generalized eigenvalue problem,
                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.

Definition at line 175 of file zggbal.f.

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 LSCALE( * ), RSCALE( * ), WORK( * )
188  COMPLEX*16 A( LDA, * ), B( LDB, * )
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  COMPLEX*16 CZERO
199  parameter( czero = ( 0.0d+0, 0.0d+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  DOUBLE PRECISION 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*16 CDUM
209 * ..
210 * .. External Functions ..
211  LOGICAL LSAME
212  INTEGER IZAMAX
213  DOUBLE PRECISION DDOT, DLAMCH
214  EXTERNAL lsame, izamax, ddot, dlamch
215 * ..
216 * .. External Subroutines ..
217  EXTERNAL daxpy, dscal, xerbla, zdscal, zswap
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC abs, dble, dimag, int, log10, max, min, sign
221 * ..
222 * .. Statement Functions ..
223  DOUBLE PRECISION CABS1
224 * ..
225 * .. Statement Function definitions ..
226  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( 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( 'ZGGBAL', -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 ) = 1
291  lscale( 1 ) = 1
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 zswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
353  CALL zswap( 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 zswap( l, a( 1, j ), 1, a( 1, m ), 1 )
362  CALL zswap( 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 / dble( 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 = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
433  $ ddot( 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 dscal( nr, beta, work( ilo ), 1 )
451  CALL dscal( nr, beta, work( ilo+n ), 1 )
452 *
453  CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
454  CALL daxpy( 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 ) = dble( 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 ) = dble( kount )*work( j ) + sum
495  330 CONTINUE
496 *
497  sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
498  $ ddot( 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 daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
518  CALL daxpy( 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 = dlamch( '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 = izamax( n-ilo+1, a( i, ilo ), lda )
534  rab = abs( a( i, irab+ilo-1 ) )
535  irab = izamax( 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 = int(lscale( i ) + sign( half, lscale( i ) ))
539  ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
540  lscale( i ) = sclfac**ir
541  icab = izamax( ihi, a( 1, i ), 1 )
542  cab = abs( a( icab, i ) )
543  icab = izamax( ihi, b( 1, i ), 1 )
544  cab = max( cab, abs( b( icab, i ) ) )
545  lcab = int( log10( cab+sfmin ) / basl+one )
546  jc = int(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 zdscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
555  CALL zdscal( 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 zdscal( ihi, rscale( j ), a( 1, j ), 1 )
562  CALL zdscal( ihi, rscale( j ), b( 1, j ), 1 )
563  380 CONTINUE
564 *
565  RETURN
566 *
567 * End of ZGGBAL
568 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:78
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
Here is the call graph for this function:
Here is the caller graph for this function: