Return to dggbal.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b DGGBAL
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DGGBAL + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggbal.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggbal.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggbal.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
   22: *                          RSCALE, WORK, INFO )
   23: * 
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          JOB
   26: *       INTEGER            IHI, ILO, INFO, LDA, LDB, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), LSCALE( * ),
   30: *      $                   RSCALE( * ), WORK( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DGGBAL balances a pair of general real matrices (A,B).  This
   40: *> involves, first, permuting A and B by similarity transformations to
   41: *> isolate eigenvalues in the first 1 to ILO$-$1 and last IHI+1 to N
   42: *> elements on the diagonal; and second, applying a diagonal similarity
   43: *> transformation to rows and columns ILO to IHI to make the rows
   44: *> and columns as close in norm as possible. Both steps are optional.
   45: *>
   46: *> Balancing may reduce the 1-norm of the matrices, and improve the
   47: *> accuracy of the computed eigenvalues and/or eigenvectors in the
   48: *> generalized eigenvalue problem A*x = lambda*B*x.
   49: *> \endverbatim
   50: *
   51: *  Arguments:
   52: *  ==========
   53: *
   54: *> \param[in] JOB
   55: *> \verbatim
   56: *>          JOB is CHARACTER*1
   57: *>          Specifies the operations to be performed on A and B:
   58: *>          = 'N':  none:  simply set ILO = 1, IHI = N, LSCALE(I) = 1.0
   59: *>                  and RSCALE(I) = 1.0 for i = 1,...,N.
   60: *>          = 'P':  permute only;
   61: *>          = 'S':  scale only;
   62: *>          = 'B':  both permute and scale.
   63: *> \endverbatim
   64: *>
   65: *> \param[in] N
   66: *> \verbatim
   67: *>          N is INTEGER
   68: *>          The order of the matrices A and B.  N >= 0.
   69: *> \endverbatim
   70: *>
   71: *> \param[in,out] A
   72: *> \verbatim
   73: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
   74: *>          On entry, the input matrix A.
   75: *>          On exit,  A is overwritten by the balanced matrix.
   76: *>          If JOB = 'N', A is not referenced.
   77: *> \endverbatim
   78: *>
   79: *> \param[in] LDA
   80: *> \verbatim
   81: *>          LDA is INTEGER
   82: *>          The leading dimension of the array A. LDA >= max(1,N).
   83: *> \endverbatim
   84: *>
   85: *> \param[in,out] B
   86: *> \verbatim
   87: *>          B is DOUBLE PRECISION array, dimension (LDB,N)
   88: *>          On entry, the input matrix B.
   89: *>          On exit,  B is overwritten by the balanced matrix.
   90: *>          If JOB = 'N', B is not referenced.
   91: *> \endverbatim
   92: *>
   93: *> \param[in] LDB
   94: *> \verbatim
   95: *>          LDB is INTEGER
   96: *>          The leading dimension of the array B. LDB >= max(1,N).
   97: *> \endverbatim
   98: *>
   99: *> \param[out] ILO
  100: *> \verbatim
  101: *>          ILO is INTEGER
  102: *> \endverbatim
  103: *>
  104: *> \param[out] IHI
  105: *> \verbatim
  106: *>          IHI is INTEGER
  107: *>          ILO and IHI are set to integers such that on exit
  108: *>          A(i,j) = 0 and B(i,j) = 0 if i > j and
  109: *>          j = 1,...,ILO-1 or i = IHI+1,...,N.
  110: *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  111: *> \endverbatim
  112: *>
  113: *> \param[out] LSCALE
  114: *> \verbatim
  115: *>          LSCALE is DOUBLE PRECISION array, dimension (N)
  116: *>          Details of the permutations and scaling factors applied
  117: *>          to the left side of A and B.  If P(j) is the index of the
  118: *>          row interchanged with row j, and D(j)
  119: *>          is the scaling factor applied to row j, then
  120: *>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
  121: *>                      = D(j)    for J = ILO,...,IHI
  122: *>                      = P(j)    for J = IHI+1,...,N.
  123: *>          The order in which the interchanges are made is N to IHI+1,
  124: *>          then 1 to ILO-1.
  125: *> \endverbatim
  126: *>
  127: *> \param[out] RSCALE
  128: *> \verbatim
  129: *>          RSCALE is DOUBLE PRECISION array, dimension (N)
  130: *>          Details of the permutations and scaling factors applied
  131: *>          to the right side of A and B.  If P(j) is the index of the
  132: *>          column interchanged with column j, and D(j)
  133: *>          is the scaling factor applied to column j, then
  134: *>            LSCALE(j) = P(j)    for J = 1,...,ILO-1
  135: *>                      = D(j)    for J = ILO,...,IHI
  136: *>                      = P(j)    for J = IHI+1,...,N.
  137: *>          The order in which the interchanges are made is N to IHI+1,
  138: *>          then 1 to ILO-1.
  139: *> \endverbatim
  140: *>
  141: *> \param[out] WORK
  142: *> \verbatim
  143: *>          WORK is DOUBLE PRECISION array, dimension (lwork)
  144: *>          lwork must be at least max(1,6*N) when JOB = 'S' or 'B', and
  145: *>          at least 1 when JOB = 'N' or 'P'.
  146: *> \endverbatim
  147: *>
  148: *> \param[out] INFO
  149: *> \verbatim
  150: *>          INFO is INTEGER
  151: *>          = 0:  successful exit
  152: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
  153: *> \endverbatim
  154: *
  155: *  Authors:
  156: *  ========
  157: *
  158: *> \author Univ. of Tennessee 
  159: *> \author Univ. of California Berkeley 
  160: *> \author Univ. of Colorado Denver 
  161: *> \author NAG Ltd. 
  162: *
  163: *> \date November 2011
  164: *
  165: *> \ingroup doubleGBcomputational
  166: *
  167: *> \par Further Details:
  168: *  =====================
  169: *>
  170: *> \verbatim
  171: *>
  172: *>  See R.C. WARD, Balancing the generalized eigenvalue problem,
  173: *>                 SIAM J. Sci. Stat. Comp. 2 (1981), 141-152.
  174: *> \endverbatim
  175: *>
  176: *  =====================================================================
  177:       SUBROUTINE DGGBAL( JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE,
  178:      $                   RSCALE, WORK, INFO )
  179: *
  180: *  -- LAPACK computational routine (version 3.4.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 2011
  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 = 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 = 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: *
  559:       END