Annotation of rpl/lapack/lapack/dlaln2.f, revision 1.9

1.9     ! bertrand    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: *  =====================================================================
1.1       bertrand  218:       SUBROUTINE DLALN2( LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B,
                    219:      $                   LDB, WR, WI, X, LDX, SCALE, XNORM, INFO )
                    220: *
1.9     ! bertrand  221: *  -- LAPACK auxiliary routine (version 3.4.0) --
1.1       bertrand  222: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    223: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9     ! bertrand  224: *     November 2011
1.1       bertrand  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: *
1.8       bertrand  363: *        Compute the real part of  C = ca A - w D  (or  ca A**T - w D )
1.1       bertrand  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>