zgbsvx.f.html

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

</span><span class="comment">*</span><span class="comment">          = 'N':  No equilibration (always true if FACT = 'N').
</span><span class="comment">*</span><span class="comment">          = 'R':  Row equilibration, i.e., A has been premultiplied by
</span><span class="comment">*</span><span class="comment">                  diag(R).
</span><span class="comment">*</span><span class="comment">          = 'C':  Column equilibration, i.e., A has been postmultiplied
</span><span class="comment">*</span><span class="comment">                  by diag(C).
</span><span class="comment">*</span><span class="comment">          = 'B':  Both row and column equilibration, i.e., A has been
</span><span class="comment">*</span><span class="comment">                  replaced by diag(R) * A * diag(C).
</span><span class="comment">*</span><span class="comment">          EQUED is an input argument if FACT = 'F'; otherwise, it is an
</span><span class="comment">*</span><span class="comment">          output argument.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  R       (input or output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The row scale factors for A.  If EQUED = 'R' or 'B', A is
</span><span class="comment">*</span><span class="comment">          multiplied on the left by diag(R); if EQUED = 'N' or 'C', R
</span><span class="comment">*</span><span class="comment">          is not accessed.  R is an input argument if FACT = 'F';
</span><span class="comment">*</span><span class="comment">          otherwise, R is an output argument.  If FACT = 'F' and
</span><span class="comment">*</span><span class="comment">          EQUED = 'R' or 'B', each element of R must be positive.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input or output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The column scale factors for A.  If EQUED = 'C' or 'B', A is
</span><span class="comment">*</span><span class="comment">          multiplied on the right by diag(C); if EQUED = 'N' or 'R', C
</span><span class="comment">*</span><span class="comment">          is not accessed.  C is an input argument if FACT = 'F';
</span><span class="comment">*</span><span class="comment">          otherwise, C is an output argument.  If FACT = 'F' and
</span><span class="comment">*</span><span class="comment">          EQUED = 'C' or 'B', each element of C must be positive.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
</span><span class="comment">*</span><span class="comment">          On entry, the right hand side matrix B.
</span><span class="comment">*</span><span class="comment">          On exit,
</span><span class="comment">*</span><span class="comment">          if EQUED = 'N', B is not modified;
</span><span class="comment">*</span><span class="comment">          if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by
</span><span class="comment">*</span><span class="comment">          diag(R)*B;
</span><span class="comment">*</span><span class="comment">          if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is
</span><span class="comment">*</span><span class="comment">          overwritten by diag(C)*B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B.  LDB &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  X       (output) COMPLEX*16 array, dimension (LDX,NRHS)
</span><span class="comment">*</span><span class="comment">          If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X
</span><span class="comment">*</span><span class="comment">          to the original system of equations.  Note that A and B are
</span><span class="comment">*</span><span class="comment">          modified on exit if EQUED .ne. 'N', and the solution to the
</span><span class="comment">*</span><span class="comment">          equilibrated system is inv(diag(C))*X if TRANS = 'N' and
</span><span class="comment">*</span><span class="comment">          EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C'
</span><span class="comment">*</span><span class="comment">          and EQUED = 'R' or 'B'.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDX     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array X.  LDX &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RCOND   (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The estimate of the reciprocal condition number of the matrix
</span><span class="comment">*</span><span class="comment">          A after equilibration (if done).  If RCOND is less than the
</span><span class="comment">*</span><span class="comment">          machine precision (in particular, if RCOND = 0), the matrix
</span><span class="comment">*</span><span class="comment">          is singular to working precision.  This condition is
</span><span class="comment">*</span><span class="comment">          indicated by a return code of INFO &gt; 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  FERR    (output) DOUBLE PRECISION array, dimension (NRHS)
</span><span class="comment">*</span><span class="comment">          The estimated forward error bound for each solution vector
</span><span class="comment">*</span><span class="comment">          X(j) (the j-th column of the solution matrix X).
</span><span class="comment">*</span><span class="comment">          If XTRUE is the true solution corresponding to X(j), FERR(j)
</span><span class="comment">*</span><span class="comment">          is an estimated upper bound for the magnitude of the largest
</span><span class="comment">*</span><span class="comment">          element in (X(j) - XTRUE) divided by the magnitude of the
</span><span class="comment">*</span><span class="comment">          largest element in X(j).  The estimate is as reliable as
</span><span class="comment">*</span><span class="comment">          the estimate for RCOND, and is almost always a slight
</span><span class="comment">*</span><span class="comment">          overestimate of the true error.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  BERR    (output) DOUBLE PRECISION array, dimension (NRHS)
</span><span class="comment">*</span><span class="comment">          The componentwise relative backward error of each solution
</span><span class="comment">*</span><span class="comment">          vector X(j) (i.e., the smallest relative change in
</span><span class="comment">*</span><span class="comment">          any element of A or B that makes X(j) an exact solution).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) COMPLEX*16 array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On exit, RWORK(1) contains the reciprocal pivot growth
</span><span class="comment">*</span><span class="comment">          factor norm(A)/norm(U). The &quot;max absolute element&quot; norm is
</span><span class="comment">*</span><span class="comment">          used. If RWORK(1) is much less than 1, then the stability
</span><span class="comment">*</span><span class="comment">          of the LU factorization of the (equilibrated) matrix A
</span><span class="comment">*</span><span class="comment">          could be poor. This also means that the solution X, condition
</span><span class="comment">*</span><span class="comment">          estimator RCOND, and forward error bound FERR could be
</span><span class="comment">*</span><span class="comment">          unreliable. If factorization fails with 0&lt;INFO&lt;=N, then
</span><span class="comment">*</span><span class="comment">          RWORK(1) contains the reciprocal pivot growth factor for the
</span><span class="comment">*</span><span class="comment">          leading INFO columns of A.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0:  if INFO = i, and i is
</span><span class="comment">*</span><span class="comment">                &lt;= N:  U(i,i) is exactly zero.  The factorization
</span><span class="comment">*</span><span class="comment">                       has been completed, but the factor U is exactly
</span><span class="comment">*</span><span class="comment">                       singular, so the solution and error bounds
</span><span class="comment">*</span><span class="comment">                       could not be computed. RCOND = 0 is returned.
</span><span class="comment">*</span><span class="comment">                = N+1: U is nonsingular, but RCOND is less than machine
</span><span class="comment">*</span><span class="comment">                       precision, meaning that the matrix is singular
</span><span class="comment">*</span><span class="comment">                       to working precision.  Nevertheless, the
</span><span class="comment">*</span><span class="comment">                       solution and error bounds are computed because
</span><span class="comment">*</span><span class="comment">                       there are a number of situations where the
</span><span class="comment">*</span><span class="comment">                       computed solution can be more accurate than the
</span><span class="comment">*</span><span class="comment">                       value of RCOND would suggest.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">  Moved setting of INFO = N+1 so INFO does not subsequently get
</span><span class="comment">*</span><span class="comment">  overwritten.  Sven, 17 Mar 05. 
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            COLEQU, EQUIL, NOFACT, NOTRAN, ROWEQU
      CHARACTER          NORM
      INTEGER            I, INFEQU, J, J1, J2
      DOUBLE PRECISION   AMAX, ANORM, BIGNUM, COLCND, RCMAX, RCMIN,
     $                   ROWCND, RPVGRW, SMLNUM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.278"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.279"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGB.279"></a><a href="zlangb.f.html#ZLANGB.1">ZLANGB</a>, <a name="ZLANTB.279"></a><a href="zlantb.f.html#ZLANTB.1">ZLANTB</a>
      EXTERNAL           <a name="LSAME.280"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.280"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGB.280"></a><a href="zlangb.f.html#ZLANGB.1">ZLANGB</a>, <a name="ZLANTB.280"></a><a href="zlantb.f.html#ZLANTB.1">ZLANTB</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="XERBLA.283"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, ZCOPY, <a name="ZGBCON.283"></a><a href="zgbcon.f.html#ZGBCON.1">ZGBCON</a>, <a name="ZGBEQU.283"></a><a href="zgbequ.f.html#ZGBEQU.1">ZGBEQU</a>, <a name="ZGBRFS.283"></a><a href="zgbrfs.f.html#ZGBRFS.1">ZGBRFS</a>, <a name="ZGBTRF.283"></a><a href="zgbtrf.f.html#ZGBTRF.1">ZGBTRF</a>,
     $                   <a name="ZGBTRS.284"></a><a href="zgbtrs.f.html#ZGBTRS.1">ZGBTRS</a>, <a name="ZLACPY.284"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLAQGB.284"></a><a href="zlaqgb.f.html#ZLAQGB.1">ZLAQGB</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN
<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>      INFO = 0
      NOFACT = <a name="LSAME.292"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'N'</span> )
      EQUIL = <a name="LSAME.293"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'E'</span> )
      NOTRAN = <a name="LSAME.294"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'N'</span> )
      IF( NOFACT .OR. EQUIL ) THEN
         EQUED = <span class="string">'N'</span>
         ROWEQU = .FALSE.
         COLEQU = .FALSE.
      ELSE
         ROWEQU = <a name="LSAME.300"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'R'</span> ) .OR. <a name="LSAME.300"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</span> )
         COLEQU = <a name="LSAME.301"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'C'</span> ) .OR. <a name="LSAME.301"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</span> )
         SMLNUM = <a name="DLAMCH.302"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
         BIGNUM = ONE / SMLNUM
      END IF
<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.308"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) )
     $     THEN
         INFO = -1
      ELSE IF( .NOT.NOTRAN .AND. .NOT.<a name="LSAME.311"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) .AND. .NOT.
     $         <a name="LSAME.312"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'C'</span> ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 ) THEN
         INFO = -3
      ELSE IF( KL.LT.0 ) THEN
         INFO = -4
      ELSE IF( KU.LT.0 ) THEN
         INFO = -5
      ELSE IF( NRHS.LT.0 ) THEN
         INFO = -6
      ELSE IF( LDAB.LT.KL+KU+1 ) THEN
         INFO = -8
      ELSE IF( LDAFB.LT.2*KL+KU+1 ) THEN
         INFO = -10
      ELSE IF( <a name="LSAME.326"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'F'</span> ) .AND. .NOT.
     $         ( ROWEQU .OR. COLEQU .OR. <a name="LSAME.327"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'N'</span> ) ) ) THEN
         INFO = -12
      ELSE
         IF( ROWEQU ) THEN
            RCMIN = BIGNUM
            RCMAX = ZERO
            DO 10 J = 1, N
               RCMIN = MIN( RCMIN, R( J ) )
               RCMAX = MAX( RCMAX, R( J ) )
   10       CONTINUE
            IF( RCMIN.LE.ZERO ) THEN
               INFO = -13
            ELSE IF( N.GT.0 ) THEN
               ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
            ELSE
               ROWCND = ONE
            END IF