LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dggbal ( character  JOB,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
integer  ILO,
integer  IHI,
double precision, dimension( * )  LSCALE,
double precision, dimension( * )  RSCALE,
double precision, dimension( * )  WORK,
integer  INFO 
)

DGGBAL

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

Purpose:
 DGGBAL balances a pair of general real 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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
            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]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
November 2015
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 dggbal.f.

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 a( lda, * ), b( ldb, * ), lscale( * ),
191  $ rscale( * ), work( * )
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 * ..
202 * .. Local Scalars ..
203  INTEGER i, icab, iflow, ip1, ir, irab, it, j, jc, jp1,
204  $ k, kount, l, lcab, lm1, lrab, lsfmax, lsfmin,
205  $ m, nr, nrp2
206  DOUBLE PRECISION alpha, basl, beta, cab, cmax, coef, coef2,
207  $ coef5, cor, ew, ewc, gamma, pgamma, rab, sfmax,
208  $ sfmin, sum, t, ta, tb, tc
209 * ..
210 * .. External Functions ..
211  LOGICAL lsame
212  INTEGER idamax
213  DOUBLE PRECISION ddot, dlamch
214  EXTERNAL lsame, idamax, ddot, dlamch
215 * ..
216 * .. External Subroutines ..
217  EXTERNAL daxpy, dscal, dswap, xerbla
218 * ..
219 * .. Intrinsic Functions ..
220  INTRINSIC abs, dble, int, log10, max, min, sign
221 * ..
222 * .. Executable Statements ..
223 *
224 * Test the input parameters
225 *
226  info = 0
227  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.lsame( job, 'P' ) .AND.
228  $ .NOT.lsame( job, 'S' ) .AND. .NOT.lsame( job, 'B' ) ) THEN
229  info = -1
230  ELSE IF( n.LT.0 ) THEN
231  info = -2
232  ELSE IF( lda.LT.max( 1, n ) ) THEN
233  info = -4
234  ELSE IF( ldb.LT.max( 1, n ) ) THEN
235  info = -6
236  END IF
237  IF( info.NE.0 ) THEN
238  CALL xerbla( 'DGGBAL', -info )
239  RETURN
240  END IF
241 *
242 * Quick return if possible
243 *
244  IF( n.EQ.0 ) THEN
245  ilo = 1
246  ihi = n
247  RETURN
248  END IF
249 *
250  IF( n.EQ.1 ) THEN
251  ilo = 1
252  ihi = n
253  lscale( 1 ) = one
254  rscale( 1 ) = one
255  RETURN
256  END IF
257 *
258  IF( lsame( job, 'N' ) ) THEN
259  ilo = 1
260  ihi = n
261  DO 10 i = 1, n
262  lscale( i ) = one
263  rscale( i ) = one
264  10 CONTINUE
265  RETURN
266  END IF
267 *
268  k = 1
269  l = n
270  IF( lsame( job, 'S' ) )
271  $ GO TO 190
272 *
273  GO TO 30
274 *
275 * Permute the matrices A and B to isolate the eigenvalues.
276 *
277 * Find row with one nonzero in columns 1 through L
278 *
279  20 CONTINUE
280  l = lm1
281  IF( l.NE.1 )
282  $ GO TO 30
283 *
284  rscale( 1 ) = one
285  lscale( 1 ) = one
286  GO TO 190
287 *
288  30 CONTINUE
289  lm1 = l - 1
290  DO 80 i = l, 1, -1
291  DO 40 j = 1, lm1
292  jp1 = j + 1
293  IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
294  $ GO TO 50
295  40 CONTINUE
296  j = l
297  GO TO 70
298 *
299  50 CONTINUE
300  DO 60 j = jp1, l
301  IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
302  $ GO TO 80
303  60 CONTINUE
304  j = jp1 - 1
305 *
306  70 CONTINUE
307  m = l
308  iflow = 1
309  GO TO 160
310  80 CONTINUE
311  GO TO 100
312 *
313 * Find column with one nonzero in rows K through N
314 *
315  90 CONTINUE
316  k = k + 1
317 *
318  100 CONTINUE
319  DO 150 j = k, l
320  DO 110 i = k, lm1
321  ip1 = i + 1
322  IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
323  $ GO TO 120
324  110 CONTINUE
325  i = l
326  GO TO 140
327  120 CONTINUE
328  DO 130 i = ip1, l
329  IF( a( i, j ).NE.zero .OR. b( i, j ).NE.zero )
330  $ GO TO 150
331  130 CONTINUE
332  i = ip1 - 1
333  140 CONTINUE
334  m = k
335  iflow = 2
336  GO TO 160
337  150 CONTINUE
338  GO TO 190
339 *
340 * Permute rows M and I
341 *
342  160 CONTINUE
343  lscale( m ) = i
344  IF( i.EQ.m )
345  $ GO TO 170
346  CALL dswap( n-k+1, a( i, k ), lda, a( m, k ), lda )
347  CALL dswap( n-k+1, b( i, k ), ldb, b( m, k ), ldb )
348 *
349 * Permute columns M and J
350 *
351  170 CONTINUE
352  rscale( m ) = j
353  IF( j.EQ.m )
354  $ GO TO 180
355  CALL dswap( l, a( 1, j ), 1, a( 1, m ), 1 )
356  CALL dswap( l, b( 1, j ), 1, b( 1, m ), 1 )
357 *
358  180 CONTINUE
359  GO TO ( 20, 90 )iflow
360 *
361  190 CONTINUE
362  ilo = k
363  ihi = l
364 *
365  IF( lsame( job, 'P' ) ) THEN
366  DO 195 i = ilo, ihi
367  lscale( i ) = one
368  rscale( i ) = one
369  195 CONTINUE
370  RETURN
371  END IF
372 *
373  IF( ilo.EQ.ihi )
374  $ RETURN
375 *
376 * Balance the submatrix in rows ILO to IHI.
377 *
378  nr = ihi - ilo + 1
379  DO 200 i = ilo, ihi
380  rscale( i ) = zero
381  lscale( i ) = zero
382 *
383  work( i ) = zero
384  work( i+n ) = zero
385  work( i+2*n ) = zero
386  work( i+3*n ) = zero
387  work( i+4*n ) = zero
388  work( i+5*n ) = zero
389  200 CONTINUE
390 *
391 * Compute right side vector in resulting linear equations
392 *
393  basl = log10( sclfac )
394  DO 240 i = ilo, ihi
395  DO 230 j = ilo, ihi
396  tb = b( i, j )
397  ta = a( i, j )
398  IF( ta.EQ.zero )
399  $ GO TO 210
400  ta = log10( abs( ta ) ) / basl
401  210 CONTINUE
402  IF( tb.EQ.zero )
403  $ GO TO 220
404  tb = log10( abs( tb ) ) / basl
405  220 CONTINUE
406  work( i+4*n ) = work( i+4*n ) - ta - tb
407  work( j+5*n ) = work( j+5*n ) - ta - tb
408  230 CONTINUE
409  240 CONTINUE
410 *
411  coef = one / dble( 2*nr )
412  coef2 = coef*coef
413  coef5 = half*coef2
414  nrp2 = nr + 2
415  beta = zero
416  it = 1
417 *
418 * Start generalized conjugate gradient iteration
419 *
420  250 CONTINUE
421 *
422  gamma = ddot( nr, work( ilo+4*n ), 1, work( ilo+4*n ), 1 ) +
423  $ ddot( nr, work( ilo+5*n ), 1, work( ilo+5*n ), 1 )
424 *
425  ew = zero
426  ewc = zero
427  DO 260 i = ilo, ihi
428  ew = ew + work( i+4*n )
429  ewc = ewc + work( i+5*n )
430  260 CONTINUE
431 *
432  gamma = coef*gamma - coef2*( ew**2+ewc**2 ) - coef5*( ew-ewc )**2
433  IF( gamma.EQ.zero )
434  $ GO TO 350
435  IF( it.NE.1 )
436  $ beta = gamma / pgamma
437  t = coef5*( ewc-three*ew )
438  tc = coef5*( ew-three*ewc )
439 *
440  CALL dscal( nr, beta, work( ilo ), 1 )
441  CALL dscal( nr, beta, work( ilo+n ), 1 )
442 *
443  CALL daxpy( nr, coef, work( ilo+4*n ), 1, work( ilo+n ), 1 )
444  CALL daxpy( nr, coef, work( ilo+5*n ), 1, work( ilo ), 1 )
445 *
446  DO 270 i = ilo, ihi
447  work( i ) = work( i ) + tc
448  work( i+n ) = work( i+n ) + t
449  270 CONTINUE
450 *
451 * Apply matrix to vector
452 *
453  DO 300 i = ilo, ihi
454  kount = 0
455  sum = zero
456  DO 290 j = ilo, ihi
457  IF( a( i, j ).EQ.zero )
458  $ GO TO 280
459  kount = kount + 1
460  sum = sum + work( j )
461  280 CONTINUE
462  IF( b( i, j ).EQ.zero )
463  $ GO TO 290
464  kount = kount + 1
465  sum = sum + work( j )
466  290 CONTINUE
467  work( i+2*n ) = dble( kount )*work( i+n ) + sum
468  300 CONTINUE
469 *
470  DO 330 j = ilo, ihi
471  kount = 0
472  sum = zero
473  DO 320 i = ilo, ihi
474  IF( a( i, j ).EQ.zero )
475  $ GO TO 310
476  kount = kount + 1
477  sum = sum + work( i+n )
478  310 CONTINUE
479  IF( b( i, j ).EQ.zero )
480  $ GO TO 320
481  kount = kount + 1
482  sum = sum + work( i+n )
483  320 CONTINUE
484  work( j+3*n ) = dble( kount )*work( j ) + sum
485  330 CONTINUE
486 *
487  sum = ddot( nr, work( ilo+n ), 1, work( ilo+2*n ), 1 ) +
488  $ ddot( nr, work( ilo ), 1, work( ilo+3*n ), 1 )
489  alpha = gamma / sum
490 *
491 * Determine correction to current iteration
492 *
493  cmax = zero
494  DO 340 i = ilo, ihi
495  cor = alpha*work( i+n )
496  IF( abs( cor ).GT.cmax )
497  $ cmax = abs( cor )
498  lscale( i ) = lscale( i ) + cor
499  cor = alpha*work( i )
500  IF( abs( cor ).GT.cmax )
501  $ cmax = abs( cor )
502  rscale( i ) = rscale( i ) + cor
503  340 CONTINUE
504  IF( cmax.LT.half )
505  $ GO TO 350
506 *
507  CALL daxpy( nr, -alpha, work( ilo+2*n ), 1, work( ilo+4*n ), 1 )
508  CALL daxpy( nr, -alpha, work( ilo+3*n ), 1, work( ilo+5*n ), 1 )
509 *
510  pgamma = gamma
511  it = it + 1
512  IF( it.LE.nrp2 )
513  $ GO TO 250
514 *
515 * End generalized conjugate gradient iteration
516 *
517  350 CONTINUE
518  sfmin = dlamch( 'S' )
519  sfmax = one / sfmin
520  lsfmin = int( log10( sfmin ) / basl+one )
521  lsfmax = int( log10( sfmax ) / basl )
522  DO 360 i = ilo, ihi
523  irab = idamax( n-ilo+1, a( i, ilo ), lda )
524  rab = abs( a( i, irab+ilo-1 ) )
525  irab = idamax( n-ilo+1, b( i, ilo ), ldb )
526  rab = max( rab, abs( b( i, irab+ilo-1 ) ) )
527  lrab = int( log10( rab+sfmin ) / basl+one )
528  ir = int(lscale( i ) + sign( half, lscale( i ) ))
529  ir = min( max( ir, lsfmin ), lsfmax, lsfmax-lrab )
530  lscale( i ) = sclfac**ir
531  icab = idamax( ihi, a( 1, i ), 1 )
532  cab = abs( a( icab, i ) )
533  icab = idamax( ihi, b( 1, i ), 1 )
534  cab = max( cab, abs( b( icab, i ) ) )
535  lcab = int( log10( cab+sfmin ) / basl+one )
536  jc = int(rscale( i ) + sign( half, rscale( i ) ))
537  jc = min( max( jc, lsfmin ), lsfmax, lsfmax-lcab )
538  rscale( i ) = sclfac**jc
539  360 CONTINUE
540 *
541 * Row scaling of matrices A and B
542 *
543  DO 370 i = ilo, ihi
544  CALL dscal( n-ilo+1, lscale( i ), a( i, ilo ), lda )
545  CALL dscal( n-ilo+1, lscale( i ), b( i, ilo ), ldb )
546  370 CONTINUE
547 *
548 * Column scaling of matrices A and B
549 *
550  DO 380 j = ilo, ihi
551  CALL dscal( ihi, rscale( j ), a( 1, j ), 1 )
552  CALL dscal( ihi, rscale( j ), b( 1, j ), 1 )
553  380 CONTINUE
554 *
555  RETURN
556 *
557 * End of DGGBAL
558 *
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:53
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:53
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:54
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:55
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: