dsgesv.f.html

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

      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.174"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>(<span class="string">'<a name="DSGESV.174"></a><a href="dsgesv.f.html#DSGESV.1">DSGESV</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="DLANGE.192"></a><a href="dlange.f.html#DLANGE.1">DLANGE</a>(<span class="string">'I'</span>,N,N,A,LDA,WORK)
      EPS = <a name="DLAMCH.193"></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="DLAG2S.204"></a><a href="dlag2s.f.html#DLAG2S.1">DLAG2S</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="DLAG2S.214"></a><a href="dlag2s.f.html#DLAG2S.1">DLAG2S</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="SGETRF.223"></a><a href="sgetrf.f.html#SGETRF.1">SGETRF</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="SGETRS.232"></a><a href="sgetrs.f.html#SGETRS.1">SGETRS</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="SLAG2D.237"></a><a href="slag2d.f.html#SLAG2D.1">SLAG2D</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="DLACPY.241"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>      CALL DGEMM(<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 = ABS(X(IDAMAX(N,X(1,I),1),I))
          RNRM = ABS(WORK(IDAMAX(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="DLAG2S.268"></a><a href="dlag2s.f.html#DLAG2S.1">DLAG2S</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="SGETRS.277"></a><a href="sgetrs.f.html#SGETRS.1">SGETRS</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="SLAG2D.283"></a><a href="slag2d.f.html#SLAG2D.1">SLAG2D</a>(N,NRHS,SWORK(PTSX),N,WORK,N,INFO)
<span class="comment">*</span><span class="comment">
</span>          CALL DAXPY(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="DLACPY.289"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</a>(<span class="string">'All'</span>,N,NRHS,B,LDB,WORK,N)
<span class="comment">*</span><span class="comment">
</span>          CALL DGEMM(<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 = ABS(X(IDAMAX(N,X(1,I),1),I))
              RNRM = ABS(WORK(IDAMAX(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="DGETRF.326"></a><a href="dgetrf.f.html#DGETRF.1">DGETRF</a>(N,N,A,LDA,IPIV,INFO)
<span class="comment">*</span><span class="comment">
</span>      CALL <a name="DLACPY.328"></a><a href="dlacpy.f.html#DLACPY.1">DLACPY</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="DGETRS.331"></a><a href="dgetrs.f.html#DGETRS.1">DGETRS</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="DSGESV.336"></a><a href="dsgesv.f.html#DSGESV.1">DSGESV</a>.
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>