dppsvx.f.html

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

<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.<a name="LSAME.257"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.<a name="LSAME.260"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) .AND. .NOT.<a name="LSAME.260"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) )
     $          THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( <a name="LSAME.267"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) .AND. .NOT.
     $         ( RCEQU .OR. <a name="LSAME.268"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'N'</span> ) ) ) THEN
         INFO = -7
      ELSE
         IF( RCEQU ) THEN
            SMIN = BIGNUM
            SMAX = ZERO
            DO 10 J = 1, N
               SMIN = MIN( SMIN, S( J ) )
               SMAX = MAX( SMAX, S( J ) )
   10       CONTINUE
            IF( SMIN.LE.ZERO ) THEN
               INFO = -8
            ELSE IF( N.GT.0 ) THEN
               SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM )
            ELSE
               SCOND = ONE
            END IF
         END IF
         IF( INFO.EQ.0 ) THEN
            IF( LDB.LT.MAX( 1, N ) ) THEN
               INFO = -10
            ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
               INFO = -12
            END IF
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.296"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DPPSVX.296"></a><a href="dppsvx.f.html#DPPSVX.1">DPPSVX</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( EQUIL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute row and column scalings to equilibrate the matrix A.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DPPEQU.304"></a><a href="dppequ.f.html#DPPEQU.1">DPPEQU</a>( UPLO, N, AP, S, SCOND, AMAX, INFEQU )
         IF( INFEQU.EQ.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Equilibrate the matrix.
</span><span class="comment">*</span><span class="comment">
</span>            CALL <a name="DLAQSP.309"></a><a href="dlaqsp.f.html#DLAQSP.1">DLAQSP</a>( UPLO, N, AP, S, SCOND, AMAX, EQUED )
            RCEQU = <a name="LSAME.310"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'Y'</span> )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Scale the right-hand side.
</span><span class="comment">*</span><span class="comment">
</span>      IF( RCEQU ) THEN
         DO 30 J = 1, NRHS
            DO 20 I = 1, N
               B( I, J ) = S( I )*B( I, J )
   20       CONTINUE
   30    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( NOFACT .OR. EQUIL ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the Cholesky factorization A = U'*U or A = L*L'.
</span><span class="comment">*</span><span class="comment">
</span>         CALL DCOPY( N*( N+1 ) / 2, AP, 1, AFP, 1 )
         CALL <a name="DPPTRF.329"></a><a href="dpptrf.f.html#DPPTRF.1">DPPTRF</a>( UPLO, N, AFP, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Return if INFO is non-zero.
</span><span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 )THEN
            RCOND = ZERO
            RETURN
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the norm of the matrix A.
</span><span class="comment">*</span><span class="comment">
</span>      ANORM = <a name="DLANSP.341"></a><a href="dlansp.f.html#DLANSP.1">DLANSP</a>( <span class="string">'I'</span>, UPLO, N, AP, WORK )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the reciprocal of the condition number of A.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DPPCON.345"></a><a href="dppcon.f.html#DPPCON.1">DPPCON</a>( UPLO, N, AFP, ANORM, RCOND, WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the solution matrix X.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLACPY.349"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>( <span class="string">'Full'</span>, N, NRHS, B, LDB, X, LDX )
      CALL <a name="DPPTRS.350"></a><a href="dpptrs.f.html#DPPTRS.1">DPPTRS</a>( UPLO, N, NRHS, AFP, X, LDX, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Use iterative refinement to improve the computed solution and
</span><span class="comment">*</span><span class="comment">     compute error bounds and backward error estimates for it.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="DPPRFS.355"></a><a href="dpprfs.f.html#DPPRFS.1">DPPRFS</a>( UPLO, N, NRHS, AP, AFP, B, LDB, X, LDX, FERR, BERR,
     $             WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Transform the solution matrix X to a solution of the original
</span><span class="comment">*</span><span class="comment">     system.
</span><span class="comment">*</span><span class="comment">
</span>      IF( RCEQU ) THEN
         DO 50 J = 1, NRHS
            DO 40 I = 1, N
               X( I, J ) = S( I )*X( I, J )
   40       CONTINUE
   50    CONTINUE
         DO 60 J = 1, NRHS
            FERR( J ) = FERR( J ) / SCOND
   60    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set INFO = N+1 if the matrix is singular to working precision.
</span><span class="comment">*</span><span class="comment">
</span>      IF( RCOND.LT.<a name="DLAMCH.374"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> ) )
     $   INFO = N + 1
<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="DPPSVX.379"></a><a href="dppsvx.f.html#DPPSVX.1">DPPSVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>