File:  [local] / rpl / lapack / lapack / dlaln2.f
Revision 1.10: download - view: text, annotated - select for diffs - revision graph
Mon Nov 21 22:19:32 2011 UTC (12 years, 6 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_8, rpl-4_1_7, rpl-4_1_6, rpl-4_1_5, rpl-4_1_4, HEAD
Cohérence

    1: *> \brief \b DLALN2
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLALN2 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaln2.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaln2.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaln2.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
   22: *                          LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
   23:    24: *       .. Scalar Arguments ..
   25: *       LOGICAL            LTRANS
   26: *       INTEGER            INFO, LDA, LDB, LDX, NA, NW
   27: *       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
   28: *       ..
   29: *       .. Array Arguments ..
   30: *       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DLALN2 solves a system of the form  (ca A - w D ) X = s B
   40: *> or (ca A**T - w D) X = s B   with possible scaling ("s") and
   41: *> perturbation of A.  (A**T means A-transpose.)
   42: *>
   43: *> A is an NA x NA real matrix, ca is a real scalar, D is an NA x NA
   44: *> real diagonal matrix, w is a real or complex value, and X and B are
   45: *> NA x 1 matrices -- real if w is real, complex if w is complex.  NA
   46: *> may be 1 or 2.
   47: *>
   48: *> If w is complex, X and B are represented as NA x 2 matrices,
   49: *> the first column of each being the real part and the second
   50: *> being the imaginary part.
   51: *>
   52: *> "s" is a scaling factor (.LE. 1), computed by DLALN2, which is
   53: *> so chosen that X can be computed without overflow.  X is further
   54: *> scaled if necessary to assure that norm(ca A - w D)*norm(X) is less
   55: *> than overflow.
   56: *>
   57: *> If both singular values of (ca A - w D) are less than SMIN,
   58: *> SMIN*identity will be used instead of (ca A - w D).  If only one
   59: *> singular value is less than SMIN, one element of (ca A - w D) will be
   60: *> perturbed enough to make the smallest singular value roughly SMIN.
   61: *> If both singular values are at least SMIN, (ca A - w D) will not be
   62: *> perturbed.  In any case, the perturbation will be at most some small
   63: *> multiple of max( SMIN, ulp*norm(ca A - w D) ).  The singular values
   64: *> are computed by infinity-norm approximations, and thus will only be
   65: *> correct to a factor of 2 or so.
   66: *>
   67: *> Note: all input quantities are assumed to be smaller than overflow
   68: *> by a reasonable factor.  (See BIGNUM.)
   69: *> \endverbatim
   70: *
   71: *  Arguments:
   72: *  ==========
   73: *
   74: *> \param[in] LTRANS
   75: *> \verbatim
   76: *>          LTRANS is LOGICAL
   77: *>          =.TRUE.:  A-transpose will be used.
   78: *>          =.FALSE.: A will be used (not transposed.)
   79: *> \endverbatim
   80: *>
   81: *> \param[in] NA
   82: *> \verbatim
   83: *>          NA is INTEGER
   84: *>          The size of the matrix A.  It may (only) be 1 or 2.
   85: *> \endverbatim
   86: *>
   87: *> \param[in] NW
   88: *> \verbatim
   89: *>          NW is INTEGER
   90: *>          1 if "w" is real, 2 if "w" is complex.  It may only be 1
   91: *>          or 2.
   92: *> \endverbatim
   93: *>
   94: *> \param[in] SMIN
   95: *> \verbatim
   96: *>          SMIN is DOUBLE PRECISION
   97: *>          The desired lower bound on the singular values of A.  This
   98: *>          should be a safe distance away from underflow or overflow,
   99: *>          say, between (underflow/machine precision) and  (machine
  100: *>          precision * overflow ).  (See BIGNUM and ULP.)
  101: *> \endverbatim
  102: *>
  103: *> \param[in] CA
  104: *> \verbatim
  105: *>          CA is DOUBLE PRECISION
  106: *>          The coefficient c, which A is multiplied by.
  107: *> \endverbatim
  108: *>
  109: *> \param[in] A
  110: *> \verbatim
  111: *>          A is DOUBLE PRECISION array, dimension (LDA,NA)
  112: *>          The NA x NA matrix A.
  113: *> \endverbatim
  114: *>
  115: *> \param[in] LDA
  116: *> \verbatim
  117: *>          LDA is INTEGER
  118: *>          The leading dimension of A.  It must be at least NA.
  119: *> \endverbatim
  120: *>
  121: *> \param[in] D1
  122: *> \verbatim
  123: *>          D1 is DOUBLE PRECISION
  124: *>          The 1,1 element in the diagonal matrix D.
  125: *> \endverbatim
  126: *>
  127: *> \param[in] D2
  128: *> \verbatim
  129: *>          D2 is DOUBLE PRECISION
  130: *>          The 2,2 element in the diagonal matrix D.  Not used if NW=1.
  131: *> \endverbatim
  132: *>
  133: *> \param[in] B
  134: *> \verbatim
  135: *>          B is DOUBLE PRECISION array, dimension (LDB,NW)
  136: *>          The NA x NW matrix B (right-hand side).  If NW=2 ("w" is
  137: *>          complex), column 1 contains the real part of B and column 2
  138: *>          contains the imaginary part.
  139: *> \endverbatim
  140: *>
  141: *> \param[in] LDB
  142: *> \verbatim
  143: *>          LDB is INTEGER
  144: *>          The leading dimension of B.  It must be at least NA.
  145: *> \endverbatim
  146: *>
  147: *> \param[in] WR
  148: *> \verbatim
  149: *>          WR is DOUBLE PRECISION
  150: *>          The real part of the scalar "w".
  151: *> \endverbatim
  152: *>
  153: *> \param[in] WI
  154: *> \verbatim
  155: *>          WI is DOUBLE PRECISION
  156: *>          The imaginary part of the scalar "w".  Not used if NW=1.
  157: *> \endverbatim
  158: *>
  159: *> \param[out] X
  160: *> \verbatim
  161: *>          X is DOUBLE PRECISION array, dimension (LDX,NW)
  162: *>          The NA x NW matrix X (unknowns), as computed by DLALN2.
  163: *>          If NW=2 ("w" is complex), on exit, column 1 will contain
  164: *>          the real part of X and column 2 will contain the imaginary
  165: *>          part.
  166: *> \endverbatim
  167: *>
  168: *> \param[in] LDX
  169: *> \verbatim
  170: *>          LDX is INTEGER
  171: *>          The leading dimension of X.  It must be at least NA.
  172: *> \endverbatim
  173: *>
  174: *> \param[out] SCALE
  175: *> \verbatim
  176: *>          SCALE is DOUBLE PRECISION
  177: *>          The scale factor that B must be multiplied by to insure
  178: *>          that overflow does not occur when computing X.  Thus,
  179: *>          (ca A - w D) X  will be SCALE*B, not B (ignoring
  180: *>          perturbations of A.)  It will be at most 1.
  181: *> \endverbatim
  182: *>
  183: *> \param[out] XNORM
  184: *> \verbatim
  185: *>          XNORM is DOUBLE PRECISION
  186: *>          The infinity-norm of X, when X is regarded as an NA x NW
  187: *>          real matrix.
  188: *> \endverbatim
  189: *>
  190: *> \param[out] INFO
  191: *> \verbatim
  192: *>          INFO is INTEGER
  193: *>          An error flag.  It will be set to zero if no error occurs,
  194: *>          a negative number if an argument is in error, or a positive
  195: *>          number if  ca A - w D  had to be perturbed.
  196: *>          The possible values are:
  197: *>          = 0: No error occurred, and (ca A - w D) did not have to be
  198: *>                 perturbed.
  199: *>          = 1: (ca A - w D) had to be perturbed to make its smallest
  200: *>               (or only) singular value greater than SMIN.
  201: *>          NOTE: In the interests of speed, this routine does not
  202: *>                check the inputs for errors.
  203: *> \endverbatim
  204: *
  205: *  Authors:
  206: *  ========
  207: *
  208: *> \author Univ. of Tennessee 
  209: *> \author Univ. of California Berkeley 
  210: *> \author Univ. of Colorado Denver 
  211: *> \author NAG Ltd. 
  212: *
  213: *> \date November 2011
  214: *
  215: *> \ingroup doubleOTHERauxiliary
  216: *
  217: *  =====================================================================
  218:       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
  219:      $                   LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
  220: *
  221: *  -- LAPACK auxiliary routine (version 3.4.0) --
  222: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  223: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  224: *     November 2011
  225: *
  226: *     .. Scalar Arguments ..
  227:       LOGICAL            LTRANS
  228:       INTEGER            INFO, LDA, LDB, LDX, NA, NW
  229:       DOUBLE PRECISION   CA, D1, D2, SCALE, SMIN, WI, WR, XNORM
  230: *     ..
  231: *     .. Array Arguments ..
  232:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), X( LDX, * )
  233: *     ..
  234: *
  235: * =====================================================================
  236: *
  237: *     .. Parameters ..
  238:       DOUBLE PRECISION   ZERO, ONE
  239:       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
  240:       DOUBLE PRECISION   TWO
  241:       PARAMETER          ( TWO = 2.0D0 )
  242: *     ..
  243: *     .. Local Scalars ..
  244:       INTEGER            ICMAX, J
  245:       DOUBLE PRECISION   BBND, BI1, BI2, BIGNUM, BNORM, BR1, BR2, CI21,
  246:      $                   CI22, CMAX, CNORM, CR21, CR22, CSI, CSR, LI21,
  247:      $                   LR21, SMINI, SMLNUM, TEMP, U22ABS, UI11, UI11R,
  248:      $                   UI12, UI12S, UI22, UR11, UR11R, UR12, UR12S,
  249:      $                   UR22, XI1, XI2, XR1, XR2
  250: *     ..
  251: *     .. Local Arrays ..
  252:       LOGICAL            RSWAP( 4 ), ZSWAP( 4 )
  253:       INTEGER            IPIVOT( 4, 4 )
  254:       DOUBLE PRECISION   CI( 2, 2 ), CIV( 4 ), CR( 2, 2 ), CRV( 4 )
  255: *     ..
  256: *     .. External Functions ..
  257:       DOUBLE PRECISION   DLAMCH
  258:       EXTERNAL           DLAMCH
  259: *     ..
  260: *     .. External Subroutines ..
  261:       EXTERNAL           DLADIV
  262: *     ..
  263: *     .. Intrinsic Functions ..
  264:       INTRINSIC          ABS, MAX
  265: *     ..
  266: *     .. Equivalences ..
  267:       EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),
  268:      $                   ( CR( 1, 1 ), CRV( 1 ) )
  269: *     ..
  270: *     .. Data statements ..
  271:       DATA               ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
  272:       DATA               RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
  273:       DATA               IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
  274:      $                   3, 2, 1 /
  275: *     ..
  276: *     .. Executable Statements ..
  277: *
  278: *     Compute BIGNUM
  279: *
  280:       SMLNUM = TWO*DLAMCH( 'Safe minimum' )
  281:       BIGNUM = ONE / SMLNUM
  282:       SMINI = MAX( SMIN, SMLNUM )
  283: *
  284: *     Don't check for input errors
  285: *
  286:       INFO = 0
  287: *
  288: *     Standard Initializations
  289: *
  290:       SCALE = ONE
  291: *
  292:       IF( NA.EQ.1 ) THEN
  293: *
  294: *        1 x 1  (i.e., scalar) system   C X = B
  295: *
  296:          IF( NW.EQ.1 ) THEN
  297: *
  298: *           Real 1x1 system.
  299: *
  300: *           C = ca A - w D
  301: *
  302:             CSR = CA*A( 1, 1 ) - WR*D1
  303:             CNORM = ABS( CSR )
  304: *
  305: *           If | C | < SMINI, use C = SMINI
  306: *
  307:             IF( CNORM.LT.SMINI ) THEN
  308:                CSR = SMINI
  309:                CNORM = SMINI
  310:                INFO = 1
  311:             END IF
  312: *
  313: *           Check scaling for  X = B / C
  314: *
  315:             BNORM = ABS( B( 1, 1 ) )
  316:             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
  317:                IF( BNORM.GT.BIGNUM*CNORM )
  318:      $            SCALE = ONE / BNORM
  319:             END IF
  320: *
  321: *           Compute X
  322: *
  323:             X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
  324:             XNORM = ABS( X( 1, 1 ) )
  325:          ELSE
  326: *
  327: *           Complex 1x1 system (w is complex)
  328: *
  329: *           C = ca A - w D
  330: *
  331:             CSR = CA*A( 1, 1 ) - WR*D1
  332:             CSI = -WI*D1
  333:             CNORM = ABS( CSR ) + ABS( CSI )
  334: *
  335: *           If | C | < SMINI, use C = SMINI
  336: *
  337:             IF( CNORM.LT.SMINI ) THEN
  338:                CSR = SMINI
  339:                CSI = ZERO
  340:                CNORM = SMINI
  341:                INFO = 1
  342:             END IF
  343: *
  344: *           Check scaling for  X = B / C
  345: *
  346:             BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
  347:             IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
  348:                IF( BNORM.GT.BIGNUM*CNORM )
  349:      $            SCALE = ONE / BNORM
  350:             END IF
  351: *
  352: *           Compute X
  353: *
  354:             CALL DLADIV( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
  355:      $                   X( 1, 1 ), X( 1, 2 ) )
  356:             XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
  357:          END IF
  358: *
  359:       ELSE
  360: *
  361: *        2x2 System
  362: *
  363: *        Compute the real part of  C = ca A - w D  (or  ca A**T - w D )
  364: *
  365:          CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
  366:          CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
  367:          IF( LTRANS ) THEN
  368:             CR( 1, 2 ) = CA*A( 2, 1 )
  369:             CR( 2, 1 ) = CA*A( 1, 2 )
  370:          ELSE
  371:             CR( 2, 1 ) = CA*A( 2, 1 )
  372:             CR( 1, 2 ) = CA*A( 1, 2 )
  373:          END IF
  374: *
  375:          IF( NW.EQ.1 ) THEN
  376: *
  377: *           Real 2x2 system  (w is real)
  378: *
  379: *           Find the largest element in C
  380: *
  381:             CMAX = ZERO
  382:             ICMAX = 0
  383: *
  384:             DO 10 J = 1, 4
  385:                IF( ABS( CRV( J ) ).GT.CMAX ) THEN
  386:                   CMAX = ABS( CRV( J ) )
  387:                   ICMAX = J
  388:                END IF
  389:    10       CONTINUE
  390: *
  391: *           If norm(C) < SMINI, use SMINI*identity.
  392: *
  393:             IF( CMAX.LT.SMINI ) THEN
  394:                BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
  395:                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
  396:                   IF( BNORM.GT.BIGNUM*SMINI )
  397:      $               SCALE = ONE / BNORM
  398:                END IF
  399:                TEMP = SCALE / SMINI
  400:                X( 1, 1 ) = TEMP*B( 1, 1 )
  401:                X( 2, 1 ) = TEMP*B( 2, 1 )
  402:                XNORM = TEMP*BNORM
  403:                INFO = 1
  404:                RETURN
  405:             END IF
  406: *
  407: *           Gaussian elimination with complete pivoting.
  408: *
  409:             UR11 = CRV( ICMAX )
  410:             CR21 = CRV( IPIVOT( 2, ICMAX ) )
  411:             UR12 = CRV( IPIVOT( 3, ICMAX ) )
  412:             CR22 = CRV( IPIVOT( 4, ICMAX ) )
  413:             UR11R = ONE / UR11
  414:             LR21 = UR11R*CR21
  415:             UR22 = CR22 - UR12*LR21
  416: *
  417: *           If smaller pivot < SMINI, use SMINI
  418: *
  419:             IF( ABS( UR22 ).LT.SMINI ) THEN
  420:                UR22 = SMINI
  421:                INFO = 1
  422:             END IF
  423:             IF( RSWAP( ICMAX ) ) THEN
  424:                BR1 = B( 2, 1 )
  425:                BR2 = B( 1, 1 )
  426:             ELSE
  427:                BR1 = B( 1, 1 )
  428:                BR2 = B( 2, 1 )
  429:             END IF
  430:             BR2 = BR2 - LR21*BR1
  431:             BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
  432:             IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
  433:                IF( BBND.GE.BIGNUM*ABS( UR22 ) )
  434:      $            SCALE = ONE / BBND
  435:             END IF
  436: *
  437:             XR2 = ( BR2*SCALE ) / UR22
  438:             XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
  439:             IF( ZSWAP( ICMAX ) ) THEN
  440:                X( 1, 1 ) = XR2
  441:                X( 2, 1 ) = XR1
  442:             ELSE
  443:                X( 1, 1 ) = XR1
  444:                X( 2, 1 ) = XR2
  445:             END IF
  446:             XNORM = MAX( ABS( XR1 ), ABS( XR2 ) )
  447: *
  448: *           Further scaling if  norm(A) norm(X) > overflow
  449: *
  450:             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
  451:                IF( XNORM.GT.BIGNUM / CMAX ) THEN
  452:                   TEMP = CMAX / BIGNUM
  453:                   X( 1, 1 ) = TEMP*X( 1, 1 )
  454:                   X( 2, 1 ) = TEMP*X( 2, 1 )
  455:                   XNORM = TEMP*XNORM
  456:                   SCALE = TEMP*SCALE
  457:                END IF
  458:             END IF
  459:          ELSE
  460: *
  461: *           Complex 2x2 system  (w is complex)
  462: *
  463: *           Find the largest element in C
  464: *
  465:             CI( 1, 1 ) = -WI*D1
  466:             CI( 2, 1 ) = ZERO
  467:             CI( 1, 2 ) = ZERO
  468:             CI( 2, 2 ) = -WI*D2
  469:             CMAX = ZERO
  470:             ICMAX = 0
  471: *
  472:             DO 20 J = 1, 4
  473:                IF( ABS( CRV( J ) )+ABS( CIV( J ) ).GT.CMAX ) THEN
  474:                   CMAX = ABS( CRV( J ) ) + ABS( CIV( J ) )
  475:                   ICMAX = J
  476:                END IF
  477:    20       CONTINUE
  478: *
  479: *           If norm(C) < SMINI, use SMINI*identity.
  480: *
  481:             IF( CMAX.LT.SMINI ) THEN
  482:                BNORM = MAX( ABS( B( 1, 1 ) )+ABS( B( 1, 2 ) ),
  483:      $                 ABS( B( 2, 1 ) )+ABS( B( 2, 2 ) ) )
  484:                IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
  485:                   IF( BNORM.GT.BIGNUM*SMINI )
  486:      $               SCALE = ONE / BNORM
  487:                END IF
  488:                TEMP = SCALE / SMINI
  489:                X( 1, 1 ) = TEMP*B( 1, 1 )
  490:                X( 2, 1 ) = TEMP*B( 2, 1 )
  491:                X( 1, 2 ) = TEMP*B( 1, 2 )
  492:                X( 2, 2 ) = TEMP*B( 2, 2 )
  493:                XNORM = TEMP*BNORM
  494:                INFO = 1
  495:                RETURN
  496:             END IF
  497: *
  498: *           Gaussian elimination with complete pivoting.
  499: *
  500:             UR11 = CRV( ICMAX )
  501:             UI11 = CIV( ICMAX )
  502:             CR21 = CRV( IPIVOT( 2, ICMAX ) )
  503:             CI21 = CIV( IPIVOT( 2, ICMAX ) )
  504:             UR12 = CRV( IPIVOT( 3, ICMAX ) )
  505:             UI12 = CIV( IPIVOT( 3, ICMAX ) )
  506:             CR22 = CRV( IPIVOT( 4, ICMAX ) )
  507:             CI22 = CIV( IPIVOT( 4, ICMAX ) )
  508:             IF( ICMAX.EQ.1 .OR. ICMAX.EQ.4 ) THEN
  509: *
  510: *              Code when off-diagonals of pivoted C are real
  511: *
  512:                IF( ABS( UR11 ).GT.ABS( UI11 ) ) THEN
  513:                   TEMP = UI11 / UR11
  514:                   UR11R = ONE / ( UR11*( ONE+TEMP**2 ) )
  515:                   UI11R = -TEMP*UR11R
  516:                ELSE
  517:                   TEMP = UR11 / UI11
  518:                   UI11R = -ONE / ( UI11*( ONE+TEMP**2 ) )
  519:                   UR11R = -TEMP*UI11R
  520:                END IF
  521:                LR21 = CR21*UR11R
  522:                LI21 = CR21*UI11R
  523:                UR12S = UR12*UR11R
  524:                UI12S = UR12*UI11R
  525:                UR22 = CR22 - UR12*LR21
  526:                UI22 = CI22 - UR12*LI21
  527:             ELSE
  528: *
  529: *              Code when diagonals of pivoted C are real
  530: *
  531:                UR11R = ONE / UR11
  532:                UI11R = ZERO
  533:                LR21 = CR21*UR11R
  534:                LI21 = CI21*UR11R
  535:                UR12S = UR12*UR11R
  536:                UI12S = UI12*UR11R
  537:                UR22 = CR22 - UR12*LR21 + UI12*LI21
  538:                UI22 = -UR12*LI21 - UI12*LR21
  539:             END IF
  540:             U22ABS = ABS( UR22 ) + ABS( UI22 )
  541: *
  542: *           If smaller pivot < SMINI, use SMINI
  543: *
  544:             IF( U22ABS.LT.SMINI ) THEN
  545:                UR22 = SMINI
  546:                UI22 = ZERO
  547:                INFO = 1
  548:             END IF
  549:             IF( RSWAP( ICMAX ) ) THEN
  550:                BR2 = B( 1, 1 )
  551:                BR1 = B( 2, 1 )
  552:                BI2 = B( 1, 2 )
  553:                BI1 = B( 2, 2 )
  554:             ELSE
  555:                BR1 = B( 1, 1 )
  556:                BR2 = B( 2, 1 )
  557:                BI1 = B( 1, 2 )
  558:                BI2 = B( 2, 2 )
  559:             END IF
  560:             BR2 = BR2 - LR21*BR1 + LI21*BI1
  561:             BI2 = BI2 - LI21*BR1 - LR21*BI1
  562:             BBND = MAX( ( ABS( BR1 )+ABS( BI1 ) )*
  563:      $             ( U22ABS*( ABS( UR11R )+ABS( UI11R ) ) ),
  564:      $             ABS( BR2 )+ABS( BI2 ) )
  565:             IF( BBND.GT.ONE .AND. U22ABS.LT.ONE ) THEN
  566:                IF( BBND.GE.BIGNUM*U22ABS ) THEN
  567:                   SCALE = ONE / BBND
  568:                   BR1 = SCALE*BR1
  569:                   BI1 = SCALE*BI1
  570:                   BR2 = SCALE*BR2
  571:                   BI2 = SCALE*BI2
  572:                END IF
  573:             END IF
  574: *
  575:             CALL DLADIV( BR2, BI2, UR22, UI22, XR2, XI2 )
  576:             XR1 = UR11R*BR1 - UI11R*BI1 - UR12S*XR2 + UI12S*XI2
  577:             XI1 = UI11R*BR1 + UR11R*BI1 - UI12S*XR2 - UR12S*XI2
  578:             IF( ZSWAP( ICMAX ) ) THEN
  579:                X( 1, 1 ) = XR2
  580:                X( 2, 1 ) = XR1
  581:                X( 1, 2 ) = XI2
  582:                X( 2, 2 ) = XI1
  583:             ELSE
  584:                X( 1, 1 ) = XR1
  585:                X( 2, 1 ) = XR2
  586:                X( 1, 2 ) = XI1
  587:                X( 2, 2 ) = XI2
  588:             END IF
  589:             XNORM = MAX( ABS( XR1 )+ABS( XI1 ), ABS( XR2 )+ABS( XI2 ) )
  590: *
  591: *           Further scaling if  norm(A) norm(X) > overflow
  592: *
  593:             IF( XNORM.GT.ONE .AND. CMAX.GT.ONE ) THEN
  594:                IF( XNORM.GT.BIGNUM / CMAX ) THEN
  595:                   TEMP = CMAX / BIGNUM
  596:                   X( 1, 1 ) = TEMP*X( 1, 1 )
  597:                   X( 2, 1 ) = TEMP*X( 2, 1 )
  598:                   X( 1, 2 ) = TEMP*X( 1, 2 )
  599:                   X( 2, 2 ) = TEMP*X( 2, 2 )
  600:                   XNORM = TEMP*XNORM
  601:                   SCALE = TEMP*SCALE
  602:                END IF
  603:             END IF
  604:          END IF
  605:       END IF
  606: *
  607:       RETURN
  608: *
  609: *     End of DLALN2
  610: *
  611:       END

CVSweb interface <joel.bertrand@systella.fr>