dlaln2.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 532 行 · 第 1/3 页

</span>      INTRINSIC          ABS, MAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Equivalences ..
</span>      EQUIVALENCE        ( CI( 1, 1 ), CIV( 1 ) ),
     $                   ( CR( 1, 1 ), CRV( 1 ) )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Data statements ..
</span>      DATA               ZSWAP / .FALSE., .FALSE., .TRUE., .TRUE. /
      DATA               RSWAP / .FALSE., .TRUE., .FALSE., .TRUE. /
      DATA               IPIVOT / 1, 2, 3, 4, 2, 1, 4, 3, 3, 4, 1, 2, 4,
     $                   3, 2, 1 /
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute BIGNUM
</span><span class="comment">*</span><span class="comment">
</span>      SMLNUM = TWO*<a name="DLAMCH.176"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
      BIGNUM = ONE / SMLNUM
      SMINI = MAX( SMIN, SMLNUM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Don't check for input errors
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Standard Initializations
</span><span class="comment">*</span><span class="comment">
</span>      SCALE = ONE
<span class="comment">*</span><span class="comment">
</span>      IF( NA.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        1 x 1  (i.e., scalar) system   C X = B
</span><span class="comment">*</span><span class="comment">
</span>         IF( NW.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Real 1x1 system.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           C = ca A - w D
</span><span class="comment">*</span><span class="comment">
</span>            CSR = CA*A( 1, 1 ) - WR*D1
            CNORM = ABS( CSR )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If | C | &lt; SMINI, use C = SMINI
</span><span class="comment">*</span><span class="comment">
</span>            IF( CNORM.LT.SMINI ) THEN
               CSR = SMINI
               CNORM = SMINI
               INFO = 1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Check scaling for  X = B / C
</span><span class="comment">*</span><span class="comment">
</span>            BNORM = ABS( B( 1, 1 ) )
            IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
               IF( BNORM.GT.BIGNUM*CNORM )
     $            SCALE = ONE / BNORM
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute X
</span><span class="comment">*</span><span class="comment">
</span>            X( 1, 1 ) = ( B( 1, 1 )*SCALE ) / CSR
            XNORM = ABS( X( 1, 1 ) )
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Complex 1x1 system (w is complex)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           C = ca A - w D
</span><span class="comment">*</span><span class="comment">
</span>            CSR = CA*A( 1, 1 ) - WR*D1
            CSI = -WI*D1
            CNORM = ABS( CSR ) + ABS( CSI )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If | C | &lt; SMINI, use C = SMINI
</span><span class="comment">*</span><span class="comment">
</span>            IF( CNORM.LT.SMINI ) THEN
               CSR = SMINI
               CSI = ZERO
               CNORM = SMINI
               INFO = 1
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Check scaling for  X = B / C
</span><span class="comment">*</span><span class="comment">
</span>            BNORM = ABS( B( 1, 1 ) ) + ABS( B( 1, 2 ) )
            IF( CNORM.LT.ONE .AND. BNORM.GT.ONE ) THEN
               IF( BNORM.GT.BIGNUM*CNORM )
     $            SCALE = ONE / BNORM
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute X
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLADIV.250"></a><a href="dladiv.f.html#DLADIV.1">DLADIV</a>( SCALE*B( 1, 1 ), SCALE*B( 1, 2 ), CSR, CSI,
     $                   X( 1, 1 ), X( 1, 2 ) )
            XNORM = ABS( X( 1, 1 ) ) + ABS( X( 1, 2 ) )
         END IF
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        2x2 System
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the real part of  C = ca A - w D  (or  ca A' - w D )
</span><span class="comment">*</span><span class="comment">
</span>         CR( 1, 1 ) = CA*A( 1, 1 ) - WR*D1
         CR( 2, 2 ) = CA*A( 2, 2 ) - WR*D2
         IF( LTRANS ) THEN
            CR( 1, 2 ) = CA*A( 2, 1 )
            CR( 2, 1 ) = CA*A( 1, 2 )
         ELSE
            CR( 2, 1 ) = CA*A( 2, 1 )
            CR( 1, 2 ) = CA*A( 1, 2 )
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( NW.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Real 2x2 system  (w is real)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Find the largest element in C
</span><span class="comment">*</span><span class="comment">
</span>            CMAX = ZERO
            ICMAX = 0
<span class="comment">*</span><span class="comment">
</span>            DO 10 J = 1, 4
               IF( ABS( CRV( J ) ).GT.CMAX ) THEN
                  CMAX = ABS( CRV( J ) )
                  ICMAX = J
               END IF
   10       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If norm(C) &lt; SMINI, use SMINI*identity.
</span><span class="comment">*</span><span class="comment">
</span>            IF( CMAX.LT.SMINI ) THEN
               BNORM = MAX( ABS( B( 1, 1 ) ), ABS( B( 2, 1 ) ) )
               IF( SMINI.LT.ONE .AND. BNORM.GT.ONE ) THEN
                  IF( BNORM.GT.BIGNUM*SMINI )
     $               SCALE = ONE / BNORM
               END IF
               TEMP = SCALE / SMINI
               X( 1, 1 ) = TEMP*B( 1, 1 )
               X( 2, 1 ) = TEMP*B( 2, 1 )
               XNORM = TEMP*BNORM
               INFO = 1
               RETURN
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Gaussian elimination with complete pivoting.
</span><span class="comment">*</span><span class="comment">
</span>            UR11 = CRV( ICMAX )
            CR21 = CRV( IPIVOT( 2, ICMAX ) )
            UR12 = CRV( IPIVOT( 3, ICMAX ) )
            CR22 = CRV( IPIVOT( 4, ICMAX ) )
            UR11R = ONE / UR11
            LR21 = UR11R*CR21
            UR22 = CR22 - UR12*LR21
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           If smaller pivot &lt; SMINI, use SMINI
</span><span class="comment">*</span><span class="comment">
</span>            IF( ABS( UR22 ).LT.SMINI ) THEN
               UR22 = SMINI
               INFO = 1
            END IF
            IF( RSWAP( ICMAX ) ) THEN
               BR1 = B( 2, 1 )
               BR2 = B( 1, 1 )
            ELSE
               BR1 = B( 1, 1 )
               BR2 = B( 2, 1 )
            END IF
            BR2 = BR2 - LR21*BR1
            BBND = MAX( ABS( BR1*( UR22*UR11R ) ), ABS( BR2 ) )
            IF( BBND.GT.ONE .AND. ABS( UR22 ).LT.ONE ) THEN
               IF( BBND.GE.BIGNUM*ABS( UR22 ) )
     $            SCALE = ONE / BBND
            END IF
<span class="comment">*</span><span class="comment">
</span>            XR2 = ( BR2*SCALE ) / UR22
            XR1 = ( SCALE*BR1 )*UR11R - XR2*( UR11R*UR12 )
            IF( ZSWAP( ICMAX ) ) THEN
               X( 1, 1 ) = XR2
               X( 2, 1 ) = XR1