dlatdf.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 )
      DOUBLE PRECISION   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           DAXPY, DCOPY, <a name="DGECON.119"></a><a href="dgecon.f.html#DGECON.1">DGECON</a>, <a name="DGESC2.119"></a><a href="dgesc2.f.html#DGESC2.1">DGESC2</a>, <a name="DLASSQ.119"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</a>, <a name="DLASWP.119"></a><a href="dlaswp.f.html#DLASWP.1">DLASWP</a>,
     $                   DSCAL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   DASUM, DDOT
      EXTERNAL           DASUM, DDOT
<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="DLASWP.135"></a><a href="dlaswp.f.html#DLASWP.1">DLASWP</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 + DDOT( N-J, Z( J+1, J ), 1, Z( J+1, J ), 1 )
            SMINU = DDOT( 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 DAXPY( 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 DCOPY( 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 DCOPY( 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="DLASWP.201"></a><a href="dlaswp.f.html#DLASWP.1">DLASWP</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="DLASSQ.205"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</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="DGECON.211"></a><a href="dgecon.f.html#DGECON.1">DGECON</a>( <span class="string">'I'</span>, N, Z, LDZ, ONE, TEMP, WORK, IWORK, INFO )
         CALL DCOPY( 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="DLASWP.216"></a><a href="dlaswp.f.html#DLASWP.1">DLASWP</a>( 1, XM, LDZ, 1, N-1, IPIV, -1 )
         TEMP = ONE / SQRT( DDOT( N, XM, 1, XM, 1 ) )
         CALL DSCAL( N, TEMP, XM, 1 )
         CALL DCOPY( N, XM, 1, XP, 1 )
         CALL DAXPY( N, ONE, RHS, 1, XP, 1 )
         CALL DAXPY( N, -ONE, XM, 1, RHS, 1 )
         CALL <a name="DGESC2.222"></a><a href="dgesc2.f.html#DGESC2.1">DGESC2</a>( N, Z, LDZ, RHS, IPIV, JPIV, TEMP )
         CALL <a name="DGESC2.223"></a><a href="dgesc2.f.html#DGESC2.1">DGESC2</a>( N, Z, LDZ, XP, IPIV, JPIV, TEMP )
         IF( DASUM( N, XP, 1 ).GT.DASUM( N, RHS, 1 ) )
     $      CALL DCOPY( 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="DLASSQ.229"></a><a href="dlassq.f.html#DLASSQ.1">DLASSQ</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="DLATDF.235"></a><a href="dlatdf.f.html#DLATDF.1">DLATDF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>