LAPACK  3.8.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.
Date
June 2016
Further Details:
  See R.C. WARD, Balancing the generalized eigenvalue problem,
                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.

Definition at line 179 of file zggbal.f.

179 *
180 * -- LAPACK computational routine (version 3.7.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 * June 2016
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:84
subroutine zswap(N, ZX, INCX, ZY, INCY)
ZSWAP
Definition: zswap.f:83
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91
integer function izamax(N, ZX, INCX)
IZAMAX
Definition: izamax.f:73
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:81
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
Definition: zdscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: