spbsvx.f.html

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

</span>      IF( .NOT.NOFACT .AND. .NOT.EQUIL .AND. .NOT.<a name="LSAME.280"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.283"></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( KD.LT.0 ) THEN
         INFO = -4
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -5
      ELSE IF( LDAB.LT.KD+1 ) THEN
         INFO = -7
      ELSE IF( LDAFB.LT.KD+1 ) THEN
         INFO = -9
      ELSE IF( <a name="LSAME.295"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) .AND. .NOT.
     $         ( RCEQU .OR. <a name="LSAME.296"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'N'</span> ) ) ) THEN
         INFO = -10
      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 = -11
            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 = -13
            ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
               INFO = -15
            END IF
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.324"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SPBSVX.324"></a><a href="spbsvx.f.html#SPBSVX.1">SPBSVX</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="SPBEQU.332"></a><a href="spbequ.f.html#SPBEQU.1">SPBEQU</a>( UPLO, N, KD, AB, LDAB, 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="SLAQSB.337"></a><a href="slaqsb.f.html#SLAQSB.1">SLAQSB</a>( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED )
            RCEQU = <a name="LSAME.338"></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>         IF( UPPER ) THEN
            DO 40 J = 1, N
               J1 = MAX( J-KD, 1 )
               CALL SCOPY( J-J1+1, AB( KD+1-J+J1, J ), 1,
     $                     AFB( KD+1-J+J1, J ), 1 )
   40       CONTINUE
         ELSE
            DO 50 J = 1, N
               J2 = MIN( J+KD, N )
               CALL SCOPY( J2-J+1, AB( 1, J ), 1, AFB( 1, J ), 1 )
   50       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="SPBTRF.369"></a><a href="spbtrf.f.html#SPBTRF.1">SPBTRF</a>( UPLO, N, KD, AFB, LDAFB, 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="SLANSB.381"></a><a href="slansb.f.html#SLANSB.1">SLANSB</a>( <span class="string">'1'</span>, UPLO, N, KD, AB, LDAB, 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="SPBCON.385"></a><a href="spbcon.f.html#SPBCON.1">SPBCON</a>( UPLO, N, KD, AFB, LDAFB, 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="SLACPY.390"></a><a href="slacpy.f.html#SLACPY.1">SLACPY</a>( <span class="string">'Full'</span>, N, NRHS, B, LDB, X, LDX )
      CALL <a name="SPBTRS.391"></a><a href="spbtrs.f.html#SPBTRS.1">SPBTRS</a>( UPLO, N, KD, NRHS, AFB, LDAFB, 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="SPBRFS.396"></a><a href="spbrfs.f.html#SPBRFS.1">SPBRFS</a>( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, 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 70 J = 1, NRHS
            DO 60 I = 1, N
               X( I, J ) = S( I )*X( I, J )
   60       CONTINUE
   70    CONTINUE
         DO 80 J = 1, NRHS
            FERR( J ) = FERR( J ) / SCOND
   80    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="SLAMCH.415"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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="SPBSVX.420"></a><a href="spbsvx.f.html#SPBSVX.1">SPBSVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>