clargv.f.html

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

         SAFMN2 = <a name="SLAMCH.108"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'B'</span> )**INT( LOG( SAFMIN / EPS ) /
     $            LOG( <a name="SLAMCH.109"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'B'</span> ) ) / TWO )
         SAFMX2 = ONE / SAFMN2
<span class="comment">*</span><span class="comment">     END IF
</span>      IX = 1
      IY = 1
      IC = 1
      DO 60 I = 1, N
         F = X( IX )
         G = Y( IY )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use identical algorithm as in <a name="CLARTG.119"></a><a href="clartg.f.html#CLARTG.1">CLARTG</a>
</span><span class="comment">*</span><span class="comment">
</span>         SCALE = MAX( ABS1( F ), ABS1( G ) )
         FS = F
         GS = G
         COUNT = 0
         IF( SCALE.GE.SAFMX2 ) THEN
   10       CONTINUE
            COUNT = COUNT + 1
            FS = FS*SAFMN2
            GS = GS*SAFMN2
            SCALE = SCALE*SAFMN2
            IF( SCALE.GE.SAFMX2 )
     $         GO TO 10
         ELSE IF( SCALE.LE.SAFMN2 ) THEN
            IF( G.EQ.CZERO ) THEN
               CS = ONE
               SN = CZERO
               R = F
               GO TO 50
            END IF
   20       CONTINUE
            COUNT = COUNT - 1
            FS = FS*SAFMX2
            GS = GS*SAFMX2
            SCALE = SCALE*SAFMX2
            IF( SCALE.LE.SAFMN2 )
     $         GO TO 20
         END IF
         F2 = ABSSQ( FS )
         G2 = ABSSQ( GS )
         IF( F2.LE.MAX( G2, ONE )*SAFMIN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           This is a rare case: F is very small.
</span><span class="comment">*</span><span class="comment">
</span>            IF( F.EQ.CZERO ) THEN
               CS = ZERO
               R = <a name="SLAPY2.156"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( REAL( G ), AIMAG( G ) )
<span class="comment">*</span><span class="comment">              Do complex/real division explicitly with two real
</span><span class="comment">*</span><span class="comment">              divisions
</span>               D = <a name="SLAPY2.159"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( REAL( GS ), AIMAG( GS ) )
               SN = CMPLX( REAL( GS ) / D, -AIMAG( GS ) / D )
               GO TO 50
            END IF
            F2S = <a name="SLAPY2.163"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( REAL( FS ), AIMAG( FS ) )
<span class="comment">*</span><span class="comment">           G2 and G2S are accurate
</span><span class="comment">*</span><span class="comment">           G2 is at least SAFMIN, and G2S is at least SAFMN2
</span>            G2S = SQRT( G2 )
<span class="comment">*</span><span class="comment">           Error in CS from underflow in F2S is at most
</span><span class="comment">*</span><span class="comment">           UNFL / SAFMN2 .lt. sqrt(UNFL*EPS) .lt. EPS
</span><span class="comment">*</span><span class="comment">           If MAX(G2,ONE)=G2, then F2 .lt. G2*SAFMIN,
</span><span class="comment">*</span><span class="comment">           and so CS .lt. sqrt(SAFMIN)
</span><span class="comment">*</span><span class="comment">           If MAX(G2,ONE)=ONE, then F2 .lt. SAFMIN
</span><span class="comment">*</span><span class="comment">           and so CS .lt. sqrt(SAFMIN)/SAFMN2 = sqrt(EPS)
</span><span class="comment">*</span><span class="comment">           Therefore, CS = F2S/G2S / sqrt( 1 + (F2S/G2S)**2 ) = F2S/G2S
</span>            CS = F2S / G2S
<span class="comment">*</span><span class="comment">           Make sure abs(FF) = 1
</span><span class="comment">*</span><span class="comment">           Do complex/real division explicitly with 2 real divisions
</span>            IF( ABS1( F ).GT.ONE ) THEN
               D = <a name="SLAPY2.178"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( REAL( F ), AIMAG( F ) )
               FF = CMPLX( REAL( F ) / D, AIMAG( F ) / D )
            ELSE
               DR = SAFMX2*REAL( F )
               DI = SAFMX2*AIMAG( F )
               D = <a name="SLAPY2.183"></a><a href="slapy2.f.html#SLAPY2.1">SLAPY2</a>( DR, DI )
               FF = CMPLX( DR / D, DI / D )
            END IF
            SN = FF*CMPLX( REAL( GS ) / G2S, -AIMAG( GS ) / G2S )
            R = CS*F + SN*G
         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           This is the most common case.
</span><span class="comment">*</span><span class="comment">           Neither F2 nor F2/G2 are less than SAFMIN
</span><span class="comment">*</span><span class="comment">           F2S cannot overflow, and it is accurate
</span><span class="comment">*</span><span class="comment">
</span>            F2S = SQRT( ONE+G2 / F2 )
<span class="comment">*</span><span class="comment">           Do the F2S(real)*FS(complex) multiply with two real
</span><span class="comment">*</span><span class="comment">           multiplies
</span>            R = CMPLX( F2S*REAL( FS ), F2S*AIMAG( FS ) )
            CS = ONE / F2S
            D = F2 + G2
<span class="comment">*</span><span class="comment">           Do complex/real division explicitly with two real divisions
</span>            SN = CMPLX( REAL( R ) / D, AIMAG( R ) / D )
            SN = SN*CONJG( GS )
            IF( COUNT.NE.0 ) THEN
               IF( COUNT.GT.0 ) THEN
                  DO 30 J = 1, COUNT
                     R = R*SAFMX2
   30             CONTINUE
               ELSE
                  DO 40 J = 1, -COUNT
                     R = R*SAFMN2
   40             CONTINUE
               END IF
            END IF
         END IF
   50    CONTINUE
         C( IC ) = CS
         Y( IY ) = SN
         X( IX ) = R
         IC = IC + INCC
         IY = IY + INCY
         IX = IX + INCX
   60 CONTINUE
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="CLARGV.225"></a><a href="clargv.f.html#CLARGV.1">CLARGV</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>