zgesvx.f.html

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

</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 N-by-NRHS 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 (2*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 has
</span><span class="comment">*</span><span class="comment">                       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">
</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
      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.261"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      DOUBLE PRECISION   <a name="DLAMCH.262"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.262"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>, <a name="ZLANTR.262"></a><a href="zlantr.f.html#ZLANTR.1">ZLANTR</a>
      EXTERNAL           <a name="LSAME.263"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, <a name="DLAMCH.263"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ZLANGE.263"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>, <a name="ZLANTR.263"></a><a href="zlantr.f.html#ZLANTR.1">ZLANTR</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="XERBLA.266"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>, <a name="ZGECON.266"></a><a href="zgecon.f.html#ZGECON.1">ZGECON</a>, <a name="ZGEEQU.266"></a><a href="zgeequ.f.html#ZGEEQU.1">ZGEEQU</a>, <a name="ZGERFS.266"></a><a href="zgerfs.f.html#ZGERFS.1">ZGERFS</a>, <a name="ZGETRF.266"></a><a href="zgetrf.f.html#ZGETRF.1">ZGETRF</a>, <a name="ZGETRS.266"></a><a href="zgetrs.f.html#ZGETRS.1">ZGETRS</a>,
     $                   <a name="ZLACPY.267"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>, <a name="ZLAQGE.267"></a><a href="zlaqge.f.html#ZLAQGE.1">ZLAQGE</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          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.275"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'N'</span> )
      EQUIL = <a name="LSAME.276"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( FACT, <span class="string">'E'</span> )
      NOTRAN = <a name="LSAME.277"></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.283"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'R'</span> ) .OR. <a name="LSAME.283"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</span> )
         COLEQU = <a name="LSAME.284"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'C'</span> ) .OR. <a name="LSAME.284"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'B'</span> )
         SMLNUM = <a name="DLAMCH.285"></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.291"></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.294"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( TRANS, <span class="string">'T'</span> ) .AND. .NOT.
     $         <a name="LSAME.295"></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( NRHS.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( <a name="LSAME.305"></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.306"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( EQUED, <span class="string">'N'</span> ) ) ) THEN
         INFO = -10
      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 = -11
            ELSE IF( N.GT.0 ) THEN
               ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )