zlatdf.f.html

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

</span>      INTEGER            I, INFO, J, K
      DOUBLE PRECISION   RTEMP, SCALE, SMINU, SPLUS
      COMPLEX*16         BM, BP, PMONE, TEMP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      DOUBLE PRECISION   RWORK( MAXDIM )
      COMPLEX*16         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           ZAXPY, ZCOPY, <a name="ZGECON.124"></a><a href="zgecon.f.html#ZGECON.1">ZGECON</a>, <a name="ZGESC2.124"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>, <a name="ZLASSQ.124"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>, <a name="ZLASWP.124"></a><a href="zlaswp.f.html#ZLASWP.1">ZLASWP</a>,
     $                   ZSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   DZASUM
      COMPLEX*16         ZDOTC
      EXTERNAL           DZASUM, ZDOTC
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, 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="ZLASWP.141"></a><a href="zlaswp.f.html#ZLASWP.1">ZLASWP</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 = -CONE
         DO 10 J = 1, N - 1
            BP = RHS( J ) + CONE
            BM = RHS( J ) - CONE
            SPLUS = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Lockahead for L- part RHS(1:N-1) = +-1
</span><span class="comment">*</span><span class="comment">           SPLUS and SMIN computed more efficiently than in BSOLVE[1].
</span><span class="comment">*</span><span class="comment">
</span>            SPLUS = SPLUS + DBLE( ZDOTC( N-J, Z( J+1, J ), 1, Z( J+1,
     $              J ), 1 ) )
            SMINU = DBLE( ZDOTC( N-J, Z( J+1, J ), 1, RHS( J+1 ), 1 ) )
            SPLUS = SPLUS*DBLE( 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 we
</span><span class="comment">*</span><span class="comment">              choose -1, thereafter +1. This is a simple way to get
</span><span class="comment">*</span><span class="comment">              good estimates of matrices like Byers well-known example
</span><span class="comment">*</span><span class="comment">              (see [1]). (Not done in BSOLVE.)
</span><span class="comment">*</span><span class="comment">
</span>               RHS( J ) = RHS( J ) + PMONE
               PMONE = CONE
            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 ZAXPY( N-J, TEMP, Z( J+1, J ), 1, RHS( J+1 ), 1 )
   10    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Solve for U- part, lockahead 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 ZCOPY( N-1, RHS, 1, WORK, 1 )
         WORK( N ) = RHS( N ) + CONE
         RHS( N ) = RHS( N ) - CONE
         SPLUS = ZERO
         SMINU = ZERO
         DO 30 I = N, 1, -1
            TEMP = CONE / Z( I, I )
            WORK( I ) = WORK( I )*TEMP
            RHS( I ) = RHS( I )*TEMP
            DO 20 K = I + 1, N
               WORK( I ) = WORK( I ) - WORK( K )*( Z( I, K )*TEMP )
               RHS( I ) = RHS( I ) - RHS( K )*( Z( I, K )*TEMP )
   20       CONTINUE
            SPLUS = SPLUS + ABS( WORK( I ) )
            SMINU = SMINU + ABS( RHS( I ) )
   30    CONTINUE
         IF( SPLUS.GT.SMINU )
     $      CALL ZCOPY( N, WORK, 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="ZLASWP.206"></a><a href="zlaswp.f.html#ZLASWP.1">ZLASWP</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="ZLASSQ.210"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>( N, RHS, 1, RDSCAL, RDSUM )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     ENTRY IJOB = 2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute approximate nullvector XM of Z
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZGECON.218"></a><a href="zgecon.f.html#ZGECON.1">ZGECON</a>( <span class="string">'I'</span>, N, Z, LDZ, ONE, RTEMP, WORK, RWORK, INFO )
      CALL ZCOPY( 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="ZLASWP.223"></a><a href="zlaswp.f.html#ZLASWP.1">ZLASWP</a>( 1, XM, LDZ, 1, N-1, IPIV, -1 )
      TEMP = CONE / SQRT( ZDOTC( N, XM, 1, XM, 1 ) )
      CALL ZSCAL( N, TEMP, XM, 1 )
      CALL ZCOPY( N, XM, 1, XP, 1 )
      CALL ZAXPY( N, CONE, RHS, 1, XP, 1 )
      CALL ZAXPY( N, -CONE, XM, 1, RHS, 1 )
      CALL <a name="ZGESC2.229"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>( N, Z, LDZ, RHS, IPIV, JPIV, SCALE )
      CALL <a name="ZGESC2.230"></a><a href="zgesc2.f.html#ZGESC2.1">ZGESC2</a>( N, Z, LDZ, XP, IPIV, JPIV, SCALE )
      IF( DZASUM( N, XP, 1 ).GT.DZASUM( N, RHS, 1 ) )
     $   CALL ZCOPY( 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="ZLASSQ.236"></a><a href="zlassq.f.html#ZLASSQ.1">ZLASSQ</a>( N, RHS, 1, RDSCAL, RDSUM )
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZLATDF.239"></a><a href="zlatdf.f.html#ZLATDF.1">ZLATDF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>