zgbsvx.f.html

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

         END IF
         IF( COLEQU .AND. INFO.EQ.0 ) THEN
            RCMIN = BIGNUM
            RCMAX = ZERO
            DO 20 J = 1, N
               RCMIN = MIN( RCMIN, C( J ) )
               RCMAX = MAX( RCMAX, C( J ) )
   20       CONTINUE
            IF( RCMIN.LE.ZERO ) THEN
               INFO = -14
            ELSE IF( N.GT.0 ) THEN
               COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
            ELSE
               COLCND = ONE
            END IF
         END IF
         IF( INFO.EQ.0 ) THEN
            IF( LDB.LT.MAX( 1, N ) ) THEN
               INFO = -16
            ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
               INFO = -18
            END IF
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.370"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGBSVX.370"></a><a href="zgbsvx.f.html#ZGBSVX.1">ZGBSVX</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="ZGBEQU.378"></a><a href="zgbequ.f.html#ZGBEQU.1">ZGBEQU</a>( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
     $                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="ZLAQGB.384"></a><a href="zlaqgb.f.html#ZLAQGB.1">ZLAQGB</a>( N, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
     $                   AMAX, EQUED )
            ROWEQU = <a name="LSAME.386"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'R'</span> ) .OR. <a name="LSAME.386"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</span> )
            COLEQU = <a name="LSAME.387"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'C'</span> ) .OR. <a name="LSAME.387"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</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( NOTRAN ) THEN
         IF( ROWEQU ) THEN
            DO 40 J = 1, NRHS
               DO 30 I = 1, N
                  B( I, J ) = R( I )*B( I, J )
   30          CONTINUE
   40       CONTINUE
         END IF
      ELSE IF( COLEQU ) THEN
         DO 60 J = 1, NRHS
            DO 50 I = 1, N
               B( I, J ) = C( I )*B( I, J )
   50       CONTINUE
   60    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 LU factorization of the band matrix A.
</span><span class="comment">*</span><span class="comment">
</span>         DO 70 J = 1, N
            J1 = MAX( J-KU, 1 )
            J2 = MIN( J+KL, N )
            CALL ZCOPY( J2-J1+1, AB( KU+1-J+J1, J ), 1,
     $                  AFB( KL+KU+1-J+J1, J ), 1 )
   70    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         CALL <a name="ZGBTRF.420"></a><a href="zgbtrf.f.html#ZGBTRF.1">ZGBTRF</a>( N, N, KL, KU, AFB, LDAFB, IPIV, 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
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute the reciprocal pivot growth factor of the
</span><span class="comment">*</span><span class="comment">           leading rank-deficient INFO columns of A.
</span><span class="comment">*</span><span class="comment">
</span>            ANORM = ZERO
            DO 90 J = 1, INFO
               DO 80 I = MAX( KU+2-J, 1 ), MIN( N+KU+1-J, KL+KU+1 )
                  ANORM = MAX( ANORM, ABS( AB( I, J ) ) )
   80          CONTINUE
   90       CONTINUE
            RPVGRW = <a name="ZLANTB.435"></a><a href="zlantb.f.html#ZLANTB.1">ZLANTB</a>( <span class="string">'M'</span>, <span class="string">'U'</span>, <span class="string">'N'</span>, INFO, MIN( INFO-1, KL+KU ),
     $                       AFB( MAX( 1, KL+KU+2-INFO ), 1 ), LDAFB,
     $                       RWORK )
            IF( RPVGRW.EQ.ZERO ) THEN
               RPVGRW = ONE
            ELSE
               RPVGRW = ANORM / RPVGRW
            END IF
            RWORK( 1 ) = RPVGRW
            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 and the
</span><span class="comment">*</span><span class="comment">     reciprocal pivot growth factor RPVGRW.
</span><span class="comment">*</span><span class="comment">
</span>      IF( NOTRAN ) THEN
         NORM = <span class="string">'1'</span>
      ELSE
         NORM = <span class="string">'I'</span>
      END IF
      ANORM = <a name="ZLANGB.457"></a><a href="zlangb.f.html#ZLANGB.1">ZLANGB</a>( NORM, N, KL, KU, AB, LDAB, RWORK )
      RPVGRW = <a name="ZLANTB.458"></a><a href="zlantb.f.html#ZLANTB.1">ZLANTB</a>( <span class="string">'M'</span>, <span class="string">'U'</span>, <span class="string">'N'</span>, N, KL+KU, AFB, LDAFB, RWORK )
      IF( RPVGRW.EQ.ZERO ) THEN
         RPVGRW = ONE
      ELSE
         RPVGRW = <a name="ZLANGB.462"></a><a href="zlangb.f.html#ZLANGB.1">ZLANGB</a>( <span class="string">'M'</span>, N, KL, KU, AB, LDAB, RWORK ) / RPVGRW
      END IF
<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="ZGBCON.467"></a><a href="zgbcon.f.html#ZGBCON.1">ZGBCON</a>( NORM, N, KL, KU, AFB, LDAFB, IPIV, 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="ZLACPY.472"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>( <span class="string">'Full'</span>, N, NRHS, B, LDB, X, LDX )
      CALL <a name="ZGBTRS.473"></a><a href="zgbtrs.f.html#ZGBTRS.1">ZGBTRS</a>( TRANS, N, KL, KU, NRHS, AFB, LDAFB, IPIV, 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="ZGBRFS.479"></a><a href="zgbrfs.f.html#ZGBRFS.1">ZGBRFS</a>( TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV,
     $             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( NOTRAN ) THEN
         IF( COLEQU ) THEN
            DO 110 J = 1, NRHS
               DO 100 I = 1, N
                  X( I, J ) = C( I )*X( I, J )
  100          CONTINUE
  110       CONTINUE
            DO 120 J = 1, NRHS
               FERR( J ) = FERR( J ) / COLCND
  120       CONTINUE
         END IF
      ELSE IF( ROWEQU ) THEN
         DO 140 J = 1, NRHS
            DO 130 I = 1, N
               X( I, J ) = R( I )*X( I, J )
  130       CONTINUE
  140    CONTINUE
         DO 150 J = 1, NRHS
            FERR( J ) = FERR( J ) / ROWCND
  150    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.509"></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>      RWORK( 1 ) = RPVGRW
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZGBSVX.515"></a><a href="zgbsvx.f.html#ZGBSVX.1">ZGBSVX</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>