slasv2.f.html

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

         HT = TEMP
         TEMP = FA
         FA = HA
         HA = TEMP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now FA .ge. HA
</span><span class="comment">*</span><span class="comment">
</span>      END IF
      GT = G
      GA = ABS( GT )
      IF( GA.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Diagonal matrix
</span><span class="comment">*</span><span class="comment">
</span>         SSMIN = HA
         SSMAX = FA
         CLT = ONE
         CRT = ONE
         SLT = ZERO
         SRT = ZERO
      ELSE
         GASMAL = .TRUE.
         IF( GA.GT.FA ) THEN
            PMAX = 2
            IF( ( FA / GA ).LT.<a name="SLAMCH.143"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>( <span class="string">'EPS'</span> ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Case of very large GA
</span><span class="comment">*</span><span class="comment">
</span>               GASMAL = .FALSE.
               SSMAX = GA
               IF( HA.GT.ONE ) THEN
                  SSMIN = FA / ( GA / HA )
               ELSE
                  SSMIN = ( FA / GA )*HA
               END IF
               CLT = ONE
               SLT = HT / GT
               SRT = ONE
               CRT = FT / GT
            END IF
         END IF
         IF( GASMAL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Normal case
</span><span class="comment">*</span><span class="comment">
</span>            D = FA - HA
            IF( D.EQ.FA ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copes with infinite F or H
</span><span class="comment">*</span><span class="comment">
</span>               L = ONE
            ELSE
               L = D / FA
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that 0 .le. L .le. 1
</span><span class="comment">*</span><span class="comment">
</span>            M = GT / FT
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that abs(M) .le. 1/macheps
</span><span class="comment">*</span><span class="comment">
</span>            T = TWO - L
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that T .ge. 1
</span><span class="comment">*</span><span class="comment">
</span>            MM = M*M
            TT = T*T
            S = SQRT( TT+MM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that 1 .le. S .le. 1 + 1/macheps
</span><span class="comment">*</span><span class="comment">
</span>            IF( L.EQ.ZERO ) THEN
               R = ABS( M )
            ELSE
               R = SQRT( L*L+MM )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that 0 .le. R .le. 1 + 1/macheps
</span><span class="comment">*</span><span class="comment">
</span>            A = HALF*( S+R )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Note that 1 .le. A .le. 1 + abs(M)
</span><span class="comment">*</span><span class="comment">
</span>            SSMIN = HA / A
            SSMAX = FA*A
            IF( MM.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Note that M is very tiny
</span><span class="comment">*</span><span class="comment">
</span>               IF( L.EQ.ZERO ) THEN
                  T = SIGN( TWO, FT )*SIGN( ONE, GT )
               ELSE
                  T = GT / SIGN( D, FT ) + M / T
               END IF
            ELSE
               T = ( M / ( S+T )+M / ( R+L ) )*( ONE+A )
            END IF
            L = SQRT( T*T+FOUR )
            CRT = TWO / L
            SRT = T / L
            CLT = ( CRT+SRT*M ) / A
            SLT = ( HT / FT )*SRT / A
         END IF
      END IF
      IF( SWAP ) THEN
         CSL = SRT
         SNL = CRT
         CSR = SLT
         SNR = CLT
      ELSE
         CSL = CLT
         SNL = SLT
         CSR = CRT
         SNR = SRT
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Correct signs of SSMAX and SSMIN
</span><span class="comment">*</span><span class="comment">
</span>      IF( PMAX.EQ.1 )
     $   TSIGN = SIGN( ONE, CSR )*SIGN( ONE, CSL )*SIGN( ONE, F )
      IF( PMAX.EQ.2 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, CSL )*SIGN( ONE, G )
      IF( PMAX.EQ.3 )
     $   TSIGN = SIGN( ONE, SNR )*SIGN( ONE, SNL )*SIGN( ONE, H )
      SSMAX = SIGN( SSMAX, TSIGN )
      SSMIN = SIGN( SSMIN, TSIGN*SIGN( ONE, F )*SIGN( ONE, H ) )
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLASV2.247"></a><a href="slasv2.f.html#SLASV2.1">SLASV2</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>