cpbsvx.f.html

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

     $     THEN
         INFO = -1
      ELSE IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.282"></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.294"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) .AND. .NOT.
     $         ( RCEQU .OR. <a name="LSAME.295"></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.323"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="CPBSVX.323"></a><a href="cpbsvx.f.html#CPBSVX.1">CPBSVX</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="CPBEQU.331"></a><a href="cpbequ.f.html#CPBEQU.1">CPBEQU</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="CLAQHB.336"></a><a href="claqhb.f.html#CLAQHB.1">CLAQHB</a>( UPLO, N, KD, AB, LDAB, S, SCOND, AMAX, EQUED )
            RCEQU = <a name="LSAME.337"></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 CCOPY( 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 CCOPY( 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="CPBTRF.368"></a><a href="cpbtrf.f.html#CPBTRF.1">CPBTRF</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="CLANHB.380"></a><a href="clanhb.f.html#CLANHB.1">CLANHB</a>( <span class="string">'1'</span>, UPLO, N, KD, AB, LDAB, RWORK )
<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="CPBCON.384"></a><a href="cpbcon.f.html#CPBCON.1">CPBCON</a>( UPLO, N, KD, AFB, LDAFB, ANORM, RCOND, WORK, RWORK,
     $             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="CLACPY.389"></a><a href="clacpy.f.html#CLACPY.1">CLACPY</a>( <span class="string">'Full'</span>, N, NRHS, B, LDB, X, LDX )
      CALL <a name="CPBTRS.390"></a><a href="cpbtrs.f.html#CPBTRS.1">CPBTRS</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="CPBRFS.395"></a><a href="cpbrfs.f.html#CPBRFS.1">CPBRFS</a>( UPLO, N, KD, NRHS, AB, LDAB, AFB, LDAFB, B, LDB, X,
     $             LDX, FERR, BERR, WORK, RWORK, 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.414"></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="CPBSVX.419"></a><a href="cpbsvx.f.html#CPBSVX.1">CPBSVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>