slatdf.f.html

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

<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      INTEGER            IWORK( MAXDIM )
      REAL               WORK( 4*MAXDIM ), XM( MAXDIM ), XP( MAXDIM )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SAXPY, SCOPY, <a name="SGECON.119"></a><a href="sgecon.f.html#SGECON.1">SGECON</a>, <a name="SGESC2.119"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>, <a name="SLASSQ.119"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</a>, <a name="SLASWP.119"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a>,
     $                   SSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      REAL               SASUM, SDOT
      EXTERNAL           SASUM, SDOT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, SQRT
<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>      IF( IJOB.NE.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply permutations IPIV to RHS
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASWP.135"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a>( 1, RHS, LDZ, 1, N-1, IPIV, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve for L-part choosing RHS either to +1 or -1.
</span><span class="comment">*</span><span class="comment">
</span>         PMONE = -ONE
<span class="comment">*</span><span class="comment">
</span>         DO 10 J = 1, N - 1
            BP = RHS( J ) + ONE
            BM = RHS( J ) - ONE
            SPLUS = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Look-ahead for L-part RHS(1:N-1) = + or -1, SPLUS and
</span><span class="comment">*</span><span class="comment">           SMIN computed more efficiently than in BSOLVE [1].
</span><span class="comment">*</span><span class="comment">
</span>            SPLUS = SPLUS + SDOT( N-J, Z( J+1, J ), 1, Z( J+1, J ), 1 )
            SMINU = SDOT( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 )
            SPLUS = SPLUS*RHS( J )
            IF( SPLUS.GT.SMINU ) THEN
               RHS( J ) = BP
            ELSE IF( SMINU.GT.SPLUS ) THEN
               RHS( J ) = BM
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              In this case the updating sums are equal and we can
</span><span class="comment">*</span><span class="comment">              choose RHS(J) +1 or -1. The first time this happens
</span><span class="comment">*</span><span class="comment">              we choose -1, thereafter +1. This is a simple way to
</span><span class="comment">*</span><span class="comment">              get good estimates of matrices like Byers well-known
</span><span class="comment">*</span><span class="comment">              example (see [1]). (Not done in BSOLVE.)
</span><span class="comment">*</span><span class="comment">
</span>               RHS( J ) = RHS( J ) + PMONE
               PMONE = ONE
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the remaining r.h.s.
</span><span class="comment">*</span><span class="comment">
</span>            TEMP = -RHS( J )
            CALL SAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve for U-part, look-ahead for RHS(N) = +-1. This is not done
</span><span class="comment">*</span><span class="comment">        in BSOLVE and will hopefully give us a better estimate because
</span><span class="comment">*</span><span class="comment">        any ill-conditioning of the original matrix is transfered to U
</span><span class="comment">*</span><span class="comment">        and not to L. U(N, N) is an approximation to sigma_min(LU).
</span><span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( N-1, RHS, 1, XP, 1 )
         XP( N ) = RHS( N ) + ONE
         RHS( N ) = RHS( N ) - ONE
         SPLUS = ZERO
         SMINU = ZERO
         DO 30 I = N, 1, -1
            TEMP = ONE / Z( I, I )
            XP( I ) = XP( I )*TEMP
            RHS( I ) = RHS( I )*TEMP
            DO 20 K = I + 1, N
               XP( I ) = XP( I ) - XP( K )*( Z( I, K )*TEMP )
               RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
   20       CONTINUE
            SPLUS = SPLUS + ABS( XP( I ) )
            SMINU = SMINU + ABS( RHS( I ) )
   30    CONTINUE
         IF( SPLUS.GT.SMINU )
     $      CALL SCOPY( N, XP, 1, RHS, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply the permutations JPIV to the computed solution (RHS)
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASWP.201"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a>( 1, RHS, LDZ, 1, N-1, JPIV, -1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the sum of squares
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASSQ.205"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</a>( N, RHS, 1, RDSCAL, RDSUM )
<span class="comment">*</span><span class="comment">
</span>      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        IJOB = 2, Compute approximate nullvector XM of Z
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SGECON.211"></a><a href="sgecon.f.html#SGECON.1">SGECON</a>( <span class="string">'I'</span>, N, Z, LDZ, ONE, TEMP, WORK, IWORK, INFO )
         CALL SCOPY( N, WORK( N+1 ), 1, XM, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute RHS
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASWP.216"></a><a href="slaswp.f.html#SLASWP.1">SLASWP</a>( 1, XM, LDZ, 1, N-1, IPIV, -1 )
         TEMP = ONE / SQRT( SDOT( N, XM, 1, XM, 1 ) )
         CALL SSCAL( N, TEMP, XM, 1 )
         CALL SCOPY( N, XM, 1, XP, 1 )
         CALL SAXPY( N, ONE, RHS, 1, XP, 1 )
         CALL SAXPY( N, -ONE, XM, 1, RHS, 1 )
         CALL <a name="SGESC2.222"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( N, Z, LDZ, RHS, IPIV, JPIV, TEMP )
         CALL <a name="SGESC2.223"></a><a href="sgesc2.f.html#SGESC2.1">SGESC2</a>( N, Z, LDZ, XP, IPIV, JPIV, TEMP )
         IF( SASUM( N, XP, 1 ).GT.SASUM( N, RHS, 1 ) )
     $      CALL SCOPY( N, XP, 1, RHS, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the sum of squares
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="SLASSQ.229"></a><a href="slassq.f.html#SLASSQ.1">SLASSQ</a>( N, RHS, 1, RDSCAL, RDSUM )
<span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLATDF.235"></a><a href="slatdf.f.html#SLATDF.1">SLATDF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>