LAPACK: SRC/dggbal.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
00002      $                   RSCALE, WORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     June 2010
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOB
00011       INTEGER            IHI, ILO, INFO, LDA, LDB, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), LSCALE( * ),
00015      $                   RSCALE( * ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DGGBAL balances a pair of general real matrices (A,B).  This
00022 *  involves, first, permuting A and B by similarity transformations to
00023 *  isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
00024 *  elements on the diagonal; and second, applying a diagonal similarity
00025 *  transformation to rows and columns ILO to IHI to make the rows
00026 *  and columns as close in norm as possible. Both steps are optional.
00027 *
00028 *  Balancing may reduce the 1-norm of the matrices, and improve the
00029 *  accuracy of the computed eigenvalues and/or eigenvectors in the
00030 *  generalized eigenvalue problem A*x = lambda*B*x.
00031 *
00032 *  Arguments
00033 *  =========
00034 *
00035 *  JOB     (input) CHARACTER*1
00036 *          Specifies the operations to be performed on A and B:
00037 *          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
00038 *                  and RSCALE(I) = 1.0 for i = 1,...,N.
00039 *          = 'P':  permute only;
00040 *          = 'S':  scale only;
00041 *          = 'B':  both permute and scale.
00042 *
00043 *  N       (input) INTEGER
00044 *          The order of the matrices A and B.  N >= 0.
00045 *
00046 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00047 *          On entry, the input matrix A.
00048 *          On exit,  A is overwritten by the balanced matrix.
00049 *          If JOB = 'N', A is not referenced.
00050 *
00051 *  LDA     (input) INTEGER
00052 *          The leading dimension of the array A. LDA >= max(1,N).
00053 *
00054 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
00055 *          On entry, the input matrix B.
00056 *          On exit,  B is overwritten by the balanced matrix.
00057 *          If JOB = 'N', B is not referenced.
00058 *
00059 *  LDB     (input) INTEGER
00060 *          The leading dimension of the array B. LDB >= max(1,N).
00061 *
00062 *  ILO     (output) INTEGER
00063 *  IHI     (output) INTEGER
00064 *          ILO and IHI are set to integers such that on exit
00065 *          A(i,j) = 0 and B(i,j) = 0 if i > j and
00066 *          j = 1,...,ILO-1 or i = IHI+1,...,N.
00067 *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
00068 *
00069 *  LSCALE  (output) DOUBLE PRECISION array, dimension (N)
00070 *          Details of the permutations and scaling factors applied
00071 *          to the left side of A and B.  If P(j) is the index of the
00072 *          row interchanged with row j, and D(j)
00073 *          is the scaling factor applied to row j, then
00074 *            LSCALE(j) = P(j)    for J = 1,...,ILO-1
00075 *                      = D(j)    for J = ILO,...,IHI
00076 *                      = P(j)    for J = IHI+1,...,N.
00077 *          The order in which the interchanges are made is N to IHI+1,
00078 *          then 1 to ILO-1.
00079 *
00080 *  RSCALE  (output) DOUBLE PRECISION array, dimension (N)
00081 *          Details of the permutations and scaling factors applied
00082 *          to the right side of A and B.  If P(j) is the index of the
00083 *          column interchanged with column j, and D(j)
00084 *          is the scaling factor applied to column j, then
00085 *            LSCALE(j) = P(j)    for J = 1,...,ILO-1
00086 *                      = D(j)    for J = ILO,...,IHI
00087 *                      = P(j)    for J = IHI+1,...,N.
00088 *          The order in which the interchanges are made is N to IHI+1,
00089 *          then 1 to ILO-1.
00090 *
00091 *  WORK    (workspace) DOUBLE PRECISION array, dimension (lwork)
00092 *          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
00093 *          at least 1 when JOB = 'N' or 'P'.
00094 *
00095 *  INFO    (output) INTEGER
00096 *          = 0:  successful exit
00097 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00098 *
00099 *  Further Details
00100 *  ===============
00101 *
00102 *  See R.C. WARD, Balancing the generalized eigenvalue problem,
00103 *                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
00104 *
00105 *  =====================================================================
00106 *
00107 *     .. Parameters ..
00108       DOUBLE PRECISION   ZERO, HALF, ONE
00109       PARAMETER          ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 )
00110       DOUBLE PRECISION   THREE, SCLFAC
00111       PARAMETER          ( THREE = 3.0D+0, SCLFAC = 1.0D+1 )
00112 *     ..
00113 *     .. Local Scalars ..
00114       INTEGER            I, ICAB, IFLOW, IP1, IR, IRAB, IT, J, JC, JP1,
00115      $                   K, KOUNT, L, LCAB, LM1, LRAB, LSFMAX, LSFMIN,
00116      $                   M, NR, NRP2
00117       DOUBLE PRECISION   ALPHA, BASL, BETA, CAB, CMAX, COEF, COEF2,
00118      $                   COEF5, COR, EW, EWC, GAMMA, PGAMMA, RAB, SFMAX,
00119      $                   SFMIN, SUM, T, TA, TB, TC
00120 *     ..
00121 *     .. External Functions ..
00122       LOGICAL            LSAME
00123       INTEGER            IDAMAX
00124       DOUBLE PRECISION   DDOT, DLAMCH
00125       EXTERNAL           LSAME, IDAMAX, DDOT, DLAMCH
00126 *     ..
00127 *     .. External Subroutines ..
00128       EXTERNAL           DAXPY, DSCAL, DSWAP, XERBLA
00129 *     ..
00130 *     .. Intrinsic Functions ..
00131       INTRINSIC          ABS, DBLE, INT, LOG10, MAX, MIN, SIGN
00132 *     ..
00133 *     .. Executable Statements ..
00134 *
00135 *     Test the input parameters
00136 *
00137       INFO = 0
00138       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
00139      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
00140          INFO = -1
00141       ELSE IF( N.LT.0 ) THEN
00142          INFO = -2
00143       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00144          INFO = -4
00145       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00146          INFO = -6
00147       END IF
00148       IF( INFO.NE.0 ) THEN
00149          CALL XERBLA( 'DGGBAL', -INFO )
00150          RETURN
00151       END IF
00152 *
00153 *     Quick return if possible
00154 *
00155       IF( N.EQ.0 ) THEN
00156          ILO = 1
00157          IHI = N
00158          RETURN
00159       END IF
00160 *
00161       IF( N.EQ.1 ) THEN
00162          ILO = 1
00163          IHI = N
00164          LSCALE( 1 ) = ONE
00165          RSCALE( 1 ) = ONE
00166          RETURN
00167       END IF
00168 *
00169       IF( LSAME( JOB, 'N' ) ) THEN
00170          ILO = 1
00171          IHI = N
00172          DO 10 I = 1, N
00173             LSCALE( I ) = ONE
00174             RSCALE( I ) = ONE
00175    10    CONTINUE
00176          RETURN
00177       END IF
00178 *
00179       K = 1
00180       L = N
00181       IF( LSAME( JOB, 'S' ) )
00182      $   GO TO 190
00183 *
00184       GO TO 30
00185 *
00186 *     Permute the matrices A and B to isolate the eigenvalues.
00187 *
00188 *     Find row with one nonzero in columns 1 through L
00189 *
00190    20 CONTINUE
00191       L = LM1
00192       IF( L.NE.1 )
00193      $   GO TO 30
00194 *
00195       RSCALE( 1 ) = ONE
00196       LSCALE( 1 ) = ONE
00197       GO TO 190
00198 *
00199    30 CONTINUE
00200       LM1 = L - 1
00201       DO 80 I = L, 1, -1
00202          DO 40 J = 1, LM1
00203             JP1 = J + 1
00204             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
00205      $         GO TO 50
00206    40    CONTINUE
00207          J = L
00208          GO TO 70
00209 *
00210    50    CONTINUE
00211          DO 60 J = JP1, L
00212             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
00213      $         GO TO 80
00214    60    CONTINUE
00215          J = JP1 - 1
00216 *
00217    70    CONTINUE
00218          M = L
00219          IFLOW = 1
00220          GO TO 160
00221    80 CONTINUE
00222       GO TO 100
00223 *
00224 *     Find column with one nonzero in rows K through N
00225 *
00226    90 CONTINUE
00227       K = K + 1
00228 *
00229   100 CONTINUE
00230       DO 150 J = K, L
00231          DO 110 I = K, LM1
00232             IP1 = I + 1
00233             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
00234      $         GO TO 120
00235   110    CONTINUE
00236          I = L
00237          GO TO 140
00238   120    CONTINUE
00239          DO 130 I = IP1, L
00240             IF( A( I, J ).NE.ZERO .OR. B( I, J ).NE.ZERO )
00241      $         GO TO 150
00242   130    CONTINUE
00243          I = IP1 - 1
00244   140    CONTINUE
00245          M = K
00246          IFLOW = 2
00247          GO TO 160
00248   150 CONTINUE
00249       GO TO 190
00250 *
00251 *     Permute rows M and I
00252 *
00253   160 CONTINUE
00254       LSCALE( M ) = I
00255       IF( I.EQ.M )
00256      $   GO TO 170
00257       CALL DSWAP( N-K+1, A( I, K ), LDA, A( M, K ), LDA )
00258       CALL DSWAP( N-K+1, B( I, K ), LDB, B( M, K ), LDB )
00259 *
00260 *     Permute columns M and J
00261 *
00262   170 CONTINUE
00263       RSCALE( M ) = J
00264       IF( J.EQ.M )
00265      $   GO TO 180
00266       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
00267       CALL DSWAP( L, B( 1, J ), 1, B( 1, M ), 1 )
00268 *
00269   180 CONTINUE
00270       GO TO ( 20, 90 )IFLOW
00271 *
00272   190 CONTINUE
00273       ILO = K
00274       IHI = L
00275 *
00276       IF( LSAME( JOB, 'P' ) ) THEN
00277          DO 195 I = ILO, IHI
00278             LSCALE( I ) = ONE
00279             RSCALE( I ) = ONE
00280   195    CONTINUE
00281          RETURN
00282       END IF
00283 *
00284       IF( ILO.EQ.IHI )
00285      $   RETURN
00286 *
00287 *     Balance the submatrix in rows ILO to IHI.
00288 *
00289       NR = IHI - ILO + 1
00290       DO 200 I = ILO, IHI
00291          RSCALE( I ) = ZERO
00292          LSCALE( I ) = ZERO
00293 *
00294          WORK( I ) = ZERO
00295          WORK( I+N ) = ZERO
00296          WORK( I+2*N ) = ZERO
00297          WORK( I+3*N ) = ZERO
00298          WORK( I+4*N ) = ZERO
00299          WORK( I+5*N ) = ZERO
00300   200 CONTINUE
00301 *
00302 *     Compute right side vector in resulting linear equations
00303 *
00304       BASL = LOG10( SCLFAC )
00305       DO 240 I = ILO, IHI
00306          DO 230 J = ILO, IHI
00307             TB = B( I, J )
00308             TA = A( I, J )
00309             IF( TA.EQ.ZERO )
00310      $         GO TO 210
00311             TA = LOG10( ABS( TA ) ) / BASL
00312   210       CONTINUE
00313             IF( TB.EQ.ZERO )
00314      $         GO TO 220
00315             TB = LOG10( ABS( TB ) ) / BASL
00316   220       CONTINUE
00317             WORK( I+4*N ) = WORK( I+4*N ) - TA - TB
00318             WORK( J+5*N ) = WORK( J+5*N ) - TA - TB
00319   230    CONTINUE
00320   240 CONTINUE
00321 *
00322       COEF = ONE / DBLE( 2*NR )
00323       COEF2 = COEF*COEF
00324       COEF5 = HALF*COEF2
00325       NRP2 = NR + 2
00326       BETA = ZERO
00327       IT = 1
00328 *
00329 *     Start generalized conjugate gradient iteration
00330 *
00331   250 CONTINUE
00332 *
00333       GAMMA = DDOT( NR, WORK( ILO+4*N ), 1, WORK( ILO+4*N ), 1 ) +
00334      $        DDOT( NR, WORK( ILO+5*N ), 1, WORK( ILO+5*N ), 1 )
00335 *
00336       EW = ZERO
00337       EWC = ZERO
00338       DO 260 I = ILO, IHI
00339          EW = EW + WORK( I+4*N )
00340          EWC = EWC + WORK( I+5*N )
00341   260 CONTINUE
00342 *
00343       GAMMA = COEF*GAMMA - COEF2*( EW**2+EWC**2 ) - COEF5*( EW-EWC )**2
00344       IF( GAMMA.EQ.ZERO )
00345      $   GO TO 350
00346       IF( IT.NE.1 )
00347      $   BETA = GAMMA / PGAMMA
00348       T = COEF5*( EWC-THREE*EW )
00349       TC = COEF5*( EW-THREE*EWC )
00350 *
00351       CALL DSCAL( NR, BETA, WORK( ILO ), 1 )
00352       CALL DSCAL( NR, BETA, WORK( ILO+N ), 1 )
00353 *
00354       CALL DAXPY( NR, COEF, WORK( ILO+4*N ), 1, WORK( ILO+N ), 1 )
00355       CALL DAXPY( NR, COEF, WORK( ILO+5*N ), 1, WORK( ILO ), 1 )
00356 *
00357       DO 270 I = ILO, IHI
00358          WORK( I ) = WORK( I ) + TC
00359          WORK( I+N ) = WORK( I+N ) + T
00360   270 CONTINUE
00361 *
00362 *     Apply matrix to vector
00363 *
00364       DO 300 I = ILO, IHI
00365          KOUNT = 0
00366          SUM = ZERO
00367          DO 290 J = ILO, IHI
00368             IF( A( I, J ).EQ.ZERO )
00369      $         GO TO 280
00370             KOUNT = KOUNT + 1
00371             SUM = SUM + WORK( J )
00372   280       CONTINUE
00373             IF( B( I, J ).EQ.ZERO )
00374      $         GO TO 290
00375             KOUNT = KOUNT + 1
00376             SUM = SUM + WORK( J )
00377   290    CONTINUE
00378          WORK( I+2*N ) = DBLE( KOUNT )*WORK( I+N ) + SUM
00379   300 CONTINUE
00380 *
00381       DO 330 J = ILO, IHI
00382          KOUNT = 0
00383          SUM = ZERO
00384          DO 320 I = ILO, IHI
00385             IF( A( I, J ).EQ.ZERO )
00386      $         GO TO 310
00387             KOUNT = KOUNT + 1
00388             SUM = SUM + WORK( I+N )
00389   310       CONTINUE
00390             IF( B( I, J ).EQ.ZERO )
00391      $         GO TO 320
00392             KOUNT = KOUNT + 1
00393             SUM = SUM + WORK( I+N )
00394   320    CONTINUE
00395          WORK( J+3*N ) = DBLE( KOUNT )*WORK( J ) + SUM
00396   330 CONTINUE
00397 *
00398       SUM = DDOT( NR, WORK( ILO+N ), 1, WORK( ILO+2*N ), 1 ) +
00399      $      DDOT( NR, WORK( ILO ), 1, WORK( ILO+3*N ), 1 )
00400       ALPHA = GAMMA / SUM
00401 *
00402 *     Determine correction to current iteration
00403 *
00404       CMAX = ZERO
00405       DO 340 I = ILO, IHI
00406          COR = ALPHA*WORK( I+N )
00407          IF( ABS( COR ).GT.CMAX )
00408      $      CMAX = ABS( COR )
00409          LSCALE( I ) = LSCALE( I ) + COR
00410          COR = ALPHA*WORK( I )
00411          IF( ABS( COR ).GT.CMAX )
00412      $      CMAX = ABS( COR )
00413          RSCALE( I ) = RSCALE( I ) + COR
00414   340 CONTINUE
00415       IF( CMAX.LT.HALF )
00416      $   GO TO 350
00417 *
00418       CALL DAXPY( NR, -ALPHA, WORK( ILO+2*N ), 1, WORK( ILO+4*N ), 1 )
00419       CALL DAXPY( NR, -ALPHA, WORK( ILO+3*N ), 1, WORK( ILO+5*N ), 1 )
00420 *
00421       PGAMMA = GAMMA
00422       IT = IT + 1
00423       IF( IT.LE.NRP2 )
00424      $   GO TO 250
00425 *
00426 *     End generalized conjugate gradient iteration
00427 *
00428   350 CONTINUE
00429       SFMIN = DLAMCH( 'S' )
00430       SFMAX = ONE / SFMIN
00431       LSFMIN = INT( LOG10( SFMIN ) / BASL+ONE )
00432       LSFMAX = INT( LOG10( SFMAX ) / BASL )
00433       DO 360 I = ILO, IHI
00434          IRAB = IDAMAX( N-ILO+1, A( I, ILO ), LDA )
00435          RAB = ABS( A( I, IRAB+ILO-1 ) )
00436          IRAB = IDAMAX( N-ILO+1, B( I, ILO ), LDB )
00437          RAB = MAX( RAB, ABS( B( I, IRAB+ILO-1 ) ) )
00438          LRAB = INT( LOG10( RAB+SFMIN ) / BASL+ONE )
00439          IR = LSCALE( I ) + SIGN( HALF, LSCALE( I ) )
00440          IR = MIN( MAX( IR, LSFMIN ), LSFMAX, LSFMAX-LRAB )
00441          LSCALE( I ) = SCLFAC**IR
00442          ICAB = IDAMAX( IHI, A( 1, I ), 1 )
00443          CAB = ABS( A( ICAB, I ) )
00444          ICAB = IDAMAX( IHI, B( 1, I ), 1 )
00445          CAB = MAX( CAB, ABS( B( ICAB, I ) ) )
00446          LCAB = INT( LOG10( CAB+SFMIN ) / BASL+ONE )
00447          JC = RSCALE( I ) + SIGN( HALF, RSCALE( I ) )
00448          JC = MIN( MAX( JC, LSFMIN ), LSFMAX, LSFMAX-LCAB )
00449          RSCALE( I ) = SCLFAC**JC
00450   360 CONTINUE
00451 *
00452 *     Row scaling of matrices A and B
00453 *
00454       DO 370 I = ILO, IHI
00455          CALL DSCAL( N-ILO+1, LSCALE( I ), A( I, ILO ), LDA )
00456          CALL DSCAL( N-ILO+1, LSCALE( I ), B( I, ILO ), LDB )
00457   370 CONTINUE
00458 *
00459 *     Column scaling of matrices A and B
00460 *
00461       DO 380 J = ILO, IHI
00462          CALL DSCAL( IHI, RSCALE( J ), A( 1, J ), 1 )
00463          CALL DSCAL( IHI, RSCALE( J ), B( 1, J ), 1 )
00464   380 CONTINUE
00465 *
00466       RETURN
00467 *
00468 *     End of DGGBAL
00469 *
00470       END