zcgesv.f.html

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

</span><span class="comment">*</span><span class="comment">     Test the input parameters.
</span><span class="comment">*</span><span class="comment">
</span>      IF (N.LT.0) THEN
          INFO = -1
      ELSE IF (NRHS.LT.0) THEN
          INFO = -2
      ELSE IF (LDA.LT.MAX(1,N)) THEN
          INFO = -4
      ELSE IF (LDB.LT.MAX(1,N)) THEN
          INFO = -7
      ELSE IF (LDX.LT.MAX(1,N)) THEN
          INFO = -9
      END IF
      IF (INFO.NE.0) THEN
          CALL <a name="XERBLA.181"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>(<span class="string">'<a name="ZCGESV.181"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a>'</span>,-INFO)
          RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if (N.EQ.0).
</span><span class="comment">*</span><span class="comment">
</span>      IF (N.EQ.0) RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Skip single precision iterative refinement if a priori slower
</span><span class="comment">*</span><span class="comment">     than double precision factorization.
</span><span class="comment">*</span><span class="comment">
</span>      IF (.NOT.DOITREF) THEN
          ITER = -1
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute some constants.
</span><span class="comment">*</span><span class="comment">
</span>      ANRM = <a name="ZLANGE.199"></a><a href="zlange.f.html#ZLANGE.1">ZLANGE</a>(<span class="string">'I'</span>,N,N,A,LDA,WORK)
      EPS = <a name="DLAMCH.200"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>(<span class="string">'Epsilon'</span>)
      CTE = ANRM*EPS*SQRT(DBLE(N))*BWDMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set the pointers PTSA, PTSX for referencing SA and SX in SWORK.
</span><span class="comment">*</span><span class="comment">
</span>      PTSA = 1
      PTSX = PTSA + N*N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert B from double precision to single precision and store the
</span><span class="comment">*</span><span class="comment">     result in SX.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLAG2C.211"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,NRHS,B,LDB,SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -2
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert A from double precision to single precision and store the
</span><span class="comment">*</span><span class="comment">     result in SA.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLAG2C.221"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,N,A,LDA,SWORK(PTSA),N,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -2
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the LU factorization of SA.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGETRF.230"></a><a href="cgetrf.f.html#CGETRF.1">CGETRF</a>(N,N,SWORK(PTSA),N,IPIV,INFO)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.NE.0) THEN
          ITER = -3
          GO TO 40
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve the system SA*SX = SB.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CGETRS.239"></a><a href="cgetrs.f.html#CGETRS.1">CGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,SWORK(PTSA),N,IPIV,
     +            SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Convert SX back to double precision
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="CLAG2Z.244"></a><a href="clag2z.f.html#CLAG2Z.1">CLAG2Z</a>(N,NRHS,SWORK(PTSX),N,X,LDX,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute R = B - AX (R is WORK).
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLACPY.248"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>      CALL ZGEMM(<span class="string">'No Transpose'</span>,<span class="string">'No Transpose'</span>,N,NRHS,N,NEGONE,A,LDA,X,
     +           LDX,ONE,WORK,N)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check whether the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">     stopping criterion. If yes, set ITER=0 and return.
</span><span class="comment">*</span><span class="comment">
</span>      DO I = 1,NRHS
          XNRM = CABS1(X(IZAMAX(N,X(1,I),1),I))
          RNRM = CABS1(WORK(IZAMAX(N,WORK(1,I),1),I))
          IF (RNRM.GT.XNRM*CTE) GOTO 10
      END DO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If we are here, the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">     stopping criterion. We are good to exit.
</span><span class="comment">*</span><span class="comment">
</span>      ITER = 0
      RETURN
<span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      DO 30 IITER = 1,ITERMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Convert R (in WORK) from double precision to single precision
</span><span class="comment">*</span><span class="comment">         and store the result in SX.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="ZLAG2C.275"></a><a href="zlag2c.f.html#ZLAG2C.1">ZLAG2C</a>(N,NRHS,WORK,N,SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span>          IF (INFO.NE.0) THEN
              ITER = -2
              GO TO 40
          END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Solve the system SA*SX = SR.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="CGETRS.284"></a><a href="cgetrs.f.html#CGETRS.1">CGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,SWORK(PTSA),N,IPIV,
     +                SWORK(PTSX),N,INFO)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Convert SX back to double precision and update the current
</span><span class="comment">*</span><span class="comment">         iterate.
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="CLAG2Z.290"></a><a href="clag2z.f.html#CLAG2Z.1">CLAG2Z</a>(N,NRHS,SWORK(PTSX),N,WORK,N,INFO)
<span class="comment">*</span><span class="comment">
</span>          CALL ZAXPY(N*NRHS,ONE,WORK,1,X,1)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Compute R = B - AX (R is WORK).
</span><span class="comment">*</span><span class="comment">
</span>          CALL <a name="ZLACPY.296"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>          CALL ZGEMM(<span class="string">'No Transpose'</span>,<span class="string">'No Transpose'</span>,N,NRHS,N,NEGONE,A,
     +               LDA,X,LDX,ONE,WORK,N)
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         Check whether the NRHS normwised backward errors satisfy the
</span><span class="comment">*</span><span class="comment">         stopping criterion. If yes, set ITER=IITER&gt;0 and return.
</span><span class="comment">*</span><span class="comment">
</span>          DO I = 1,NRHS
              XNRM = CABS1(X(IZAMAX(N,X(1,I),1),I))
              RNRM = CABS1(WORK(IZAMAX(N,WORK(1,I),1),I))
              IF (RNRM.GT.XNRM*CTE) GOTO 20
          END DO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">         If we are here, the NRHS normwised backward errors satisfy the 
</span><span class="comment">*</span><span class="comment">         stopping criterion, we are good to exit.
</span><span class="comment">*</span><span class="comment">
</span>          ITER = IITER
<span class="comment">*</span><span class="comment">
</span>          RETURN
<span class="comment">*</span><span class="comment">
</span>   20     CONTINUE
<span class="comment">*</span><span class="comment">
</span>   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If we are at this place of the code, this is because we have
</span><span class="comment">*</span><span class="comment">     performed ITER=ITERMAX iterations and never satisified the stopping
</span><span class="comment">*</span><span class="comment">     criterion, set up the ITER flag accordingly and follow up on double
</span><span class="comment">*</span><span class="comment">     precision routine.
</span><span class="comment">*</span><span class="comment">
</span>      ITER = -ITERMAX - 1
<span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Single-precision iterative refinement failed to converge to a
</span><span class="comment">*</span><span class="comment">     satisfactory solution, so we resort to double precision.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZGETRF.333"></a><a href="zgetrf.f.html#ZGETRF.1">ZGETRF</a>(N,N,A,LDA,IPIV,INFO)
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZLACPY.335"></a><a href="zlacpy.f.html#ZLACPY.1">ZLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,X,LDX)
<span class="comment">*</span><span class="comment">
</span>      IF (INFO.EQ.0) THEN
          CALL <a name="ZGETRS.338"></a><a href="zgetrs.f.html#ZGETRS.1">ZGETRS</a>(<span class="string">'No transpose'</span>,N,NRHS,A,LDA,IPIV,X,LDX,INFO)
      END IF
<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="ZCGESV.343"></a><a href="zcgesv.f.html#ZCGESV.1">ZCGESV</a>.
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>