MB04DL

Balancing a real pencil, optionally avoiding large norms for the scaled (sub)matrices

[Specification] [Arguments] [Method] [References] [Comments] [Example]

Purpose

  To balance a pair of N-by-N real matrices (A,B). This involves,
  first, permuting A and B by equivalence transformations to isolate
  eigenvalues in the first 1 to ILO-1 and last IHI+1 to N elements
  on the diagonal of A and B; and second, applying a diagonal
  equivalence transformation to rows and columns ILO to IHI to make
  the rows and columns as close in 1-norm as possible. Both steps
  are optional. Balancing may reduce the 1-norms of the matrices,
  and improve the accuracy of the computed eigenvalues and/or
  eigenvectors in the generalized eigenvalue problem
  A*x = lambda*B*x.

  This routine may optionally improve the conditioning of the
  scaling transformation compared to the LAPACK routine DGGBAL.

Specification
      SUBROUTINE MB04DL( JOB, N, THRESH, A, LDA, B, LDB, ILO, IHI,
     $                   LSCALE, RSCALE, DWORK, IWARN, INFO )
C     .. Scalar Arguments ..
      CHARACTER          JOB
      INTEGER            IHI, ILO, INFO, IWARN, LDA, LDB, N
      DOUBLE PRECISION   THRESH
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), B(LDB,*), DWORK(*), LSCALE(*),
     $                   RSCALE(*)

Arguments

Mode Parameters

  JOB     CHARACTER*1
          Specifies the operations to be performed on A and B:
          = 'N':  none:  simply set ILO = 1, LSCALE(I) = 1.0 and
                  RSCALE(I) = 1.0 for I = 1,...,N.
          = 'P':  permute only;
          = 'S':  scale only;
          = 'B':  both permute and scale.

Input/Output Parameters
  N       (input) INTEGER
          The order of matrices A and B.  N >= 0.

  THRESH  (input) DOUBLE PRECISION
          If JOB = 'S' or JOB = 'B', and THRESH >= 0, threshold
          value for magnitude of the elements to be considered in
          the scaling process: elements with magnitude less than or
          equal to THRESH*MXNORM are ignored for scaling, where
          MXNORM is the maximum of the 1-norms of the original
          submatrices A(s,s) and B(s,s), with s = ILO:IHI.
          If THRESH < 0, the subroutine finds the scaling factors
          for which some conditions, detailed below, are fulfilled.
          A sequence of increasing strictly positive threshold
          values is used.
          If THRESH = -1, the condition is that
             max( norm(A(s,s),1)/norm(B(s,s),1),
                  norm(B(s,s),1)/norm(S(s,s),1) )                (1)
          has the smallest value, for the threshold values used,
          where A(s,s) and B(s,s) are the scaled submatrices.
          If THRESH = -2, the norm ratio reduction (1) is tried, but
          the subroutine may return IWARN = 1 and reset the scaling
          factors to 1, if this seems suitable. See the description
          of the argument IWARN and FURTHER COMMENTS.
          If THRESH = -3, the condition is that
             norm(A(s,s),1)*norm(B(s,s),1)                       (2)
          has the smallest value for the scaled submatrices.
          If THRESH = -4, the norm reduction in (2) is tried, but
          the subroutine may return IWARN = 1 and reset the scaling
          factors to 1, as for THRESH = -2 above.
          If THRESH = -VALUE, with VALUE >= 10, the condition
          numbers of the left and right scaling transformations will
          be bounded by VALUE, i.e., the ratios between the largest
          and smallest entries in LSCALE(s) and RSCALE(s), will be
          at most VALUE. VALUE should be a power of 10.
          If JOB = 'N' or JOB = 'P', the value of THRESH is
          irrelevant.

  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the leading N-by-N part of this array must
          contain the matrix A.
          On exit, the leading N-by-N part of this array contains
          the balanced matrix A.
          In particular, the strictly lower triangular part of the
          first ILO-1 columns and the last N-IHI rows of A is zero.

  LDA     INTEGER
          The leading dimension of the array A.  LDA >= MAX(1,N).

  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
          On entry, the leading N-by-N part of this array must
          contain the matrix B.
          On exit, the leading N-by-N part of this array contains
          the balanced matrix B.
          In particular, the strictly lower triangular part of the
          first ILO-1 columns and the last N-IHI rows of B is zero.
          If JOB = 'N', the arrays A and B are not referenced.

  LDB    INTEGER
          The leading dimension of the array B.  LDB >= MAX(1, N).

  ILO     (output) INTEGER
  IHI     (output) 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.

  LSCALE  (output) 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.

  RSCALE  (output) 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.

Workspace
  DWORK   DOUBLE PRECISION array, dimension (LDWORK) where
          LDWORK = 0,   if  JOB = 'N' or JOB = 'P', or N = 0;
          LDWORK = 6*N, if (JOB = 'S' or JOB = 'B') and THRESH >= 0;
          LDWORK = 8*N, if (JOB = 'S' or JOB = 'B') and THRESH <  0.
          On exit, if JOB = 'S' or JOB = 'B', DWORK(1) and DWORK(2)
          contain the initial 1-norms of A(s,s) and B(s,s), and
          DWORK(3) and DWORK(4) contain their final 1-norms,
          respectively. Moreover, DWORK(5) contains the THRESH value
          used (irrelevant if IWARN = 1 or ILO = IHI).

Warning Indicator
  IWARN   INTEGER
          = 0:  no warning;
          = 1:  scaling has been requested, for THRESH = -2 or
                THRESH = -4, but it most probably would not improve
                the accuracy of the computed solution for a related
                eigenproblem (since maximum norm increased
                significantly compared to the original pencil
                matrices and (very) high and/or small scaling
                factors occurred). The returned scaling factors have
                been reset to 1, but information about permutations,
                if requested, has been preserved.

Error Indicator
  INFO    INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal
                value.

Method
  Balancing consists of applying an equivalence transformation
  to isolate eigenvalues and/or to make the 1-norms of the rows
  and columns ILO,...,IHI of A and B nearly equal. If THRESH < 0,
  a search is performed to find those scaling factors giving the
  smallest norm ratio or product defined above (see the description
  of the parameter THRESH).

  Assuming JOB = 'S', let Dl and Dr be diagonal matrices containing
  the vectors LSCALE and RSCALE, respectively. The returned matrices
  are obtained using the equivalence transformation

     Dl*A*Dr and Dl*B*Dr.

  For THRESH = 0, the routine returns essentially the same results
  as the LAPACK subroutine DGGBAL [1]. Setting THRESH < 0, usually
  gives better results than DGGBAL for badly scaled matrix pencils.

References
  [1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
      Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
      Ostrouchov, S., and Sorensen, D.
      LAPACK Users' Guide: Second Edition.
      SIAM, Philadelphia, 1995.

Numerical Aspects
  No rounding errors appear if JOB = 'P'.

Further Comments
  If THRESH = -2, the increase of the maximum norm of the scaled
  submatrices, compared to the maximum norm of the initial
  submatrices, is bounded by MXGAIN = 100.
  If THRESH = -2, or THRESH = -4, the maximum condition number of
  the scaling transformations is bounded by MXCOND = 1/SQRT(EPS),
  where EPS is the machine precision (see LAPACK Library routine
  DLAMCH).

Example

Program Text

*     MB04DL EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2017 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 10 )
      INTEGER          LDA, LDB
      PARAMETER        ( LDA = NMAX, LDB = NMAX )
*     .. Local Scalars ..
      CHARACTER*1      JOB
      INTEGER          I, ILO, INFO, IWARN, J, N
      DOUBLE PRECISION THRESH
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA, NMAX), B(LDB, NMAX), DWORK(8*NMAX),
     $                 LSCALE(NMAX), RSCALE(NMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         MB04DL
*     .. Executable Statements ..
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * )  N, JOB, THRESH
      IF( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99985 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
         CALL MB04DL( JOB, N, THRESH, A, LDA, B, LDB, ILO, IHI, LSCALE,
     $                RSCALE, DWORK, IWARN, INFO )
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            WRITE ( NOUT, FMT = 99997 )
            DO 10  I = 1, N
               WRITE ( NOUT, FMT = 99995 ) ( A(I,J), J = 1,N )
10          CONTINUE
            WRITE ( NOUT, FMT = 99996 )
            DO 40  I = 1, N
               WRITE ( NOUT, FMT = 99995 ) ( B(I,J), J = 1,N )
40          CONTINUE
            WRITE ( NOUT, FMT = 99994 )  ILO
            WRITE ( NOUT, FMT = 99993 )  IHI
            WRITE ( NOUT, FMT = 99991 )
            WRITE ( NOUT, FMT = 99995 ) ( LSCALE(I), I = 1,N )
            WRITE ( NOUT, FMT = 99990 )
            WRITE ( NOUT, FMT = 99995 ) ( RSCALE(I), I = 1,N )
            IF ( LSAME( JOB, 'S' ) .OR. LSAME( JOB, 'B' ) ) THEN
               IF ( .NOT.( THRESH.EQ.-2 .OR. THRESH.EQ.-4 ) ) THEN
                  WRITE ( NOUT, FMT = 99989 )
                  WRITE ( NOUT, FMT = 99995 ) ( DWORK(I), I = 1,2 )
                  WRITE ( NOUT, FMT = 99988 )
                  WRITE ( NOUT, FMT = 99995 ) ( DWORK(I), I = 3,4 )
                  WRITE ( NOUT, FMT = 99987 )
                  WRITE ( NOUT, FMT = 99995 ) ( DWORK(5) )
               ELSE
                  WRITE ( NOUT, FMT = 99986 ) IWARN
               END IF
            END IF
         END IF
      END IF
*
99999 FORMAT (' MB04DL EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04DL = ',I2)
99997 FORMAT (' The balanced matrix A is ')
99996 FORMAT (/' The balanced matrix B is ')
99995 FORMAT (20(1X,G12.4))
99994 FORMAT (/' ILO = ',I4)
99993 FORMAT (/' IHI = ',I4)
99991 FORMAT (/' The permutations and left scaling factors are ')
99990 FORMAT (/' The permutations and right scaling factors are ')
99989 FORMAT (/' The initial 1-norms of the (sub)matrices are ')
99988 FORMAT (/' The final 1-norms of the (sub)matrices are ')
99987 FORMAT (/' The threshold value finally used is ')
99986 FORMAT (/' IWARN on exit from MB04DL = ',I2)
99985 FORMAT (/' N is out of range.',/' N = ',I5)
      END
Program Data
MB04DL EXAMPLE PROGRAM DATA
       4       B      -3
         1         0    -1e-12         0
         0        -2         0         0
        -1        -1        -1         0
        -1        -1         0         2
         1         0         0         0
         0         1         0         0
         0         0         1         0
         0         0         0         1

Program Results
 MB04DL EXAMPLE PROGRAM RESULTS

 The balanced matrix A is 
    2.000       -1.000        0.000       -1.000    
    0.000        1.000      -0.1000E-11    0.000    
    0.000       -1.000       -1.000       -1.000    
    0.000        0.000        0.000       -2.000    

 The balanced matrix B is 
    1.000        0.000        0.000        0.000    
    0.000        1.000        0.000        0.000    
    0.000        0.000        1.000        0.000    
    0.000        0.000        0.000        1.000    

 ILO =    2

 IHI =    3

 The permutations and left scaling factors are 
    2.000        1.000        1.000        2.000    

 The permutations and right scaling factors are 
    2.000        1.000        1.000        2.000    

 The initial 1-norms of the (sub)matrices are 
    2.000        1.000    

 The final 1-norms of the (sub)matrices are 
    2.000        1.000    

 The threshold value finally used is 
   0.2500E-12

Return to index