zgetri.f.html

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

      LQUERY = ( LWORK.EQ.-1 )
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -3
      ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
         INFO = -6
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.100"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGETRI.100"></a><a href="zgetri.f.html#ZGETRI.1">ZGETRI</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</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">     Form inv(U).  If INFO &gt; 0 from <a name="ZTRTRI.111"></a><a href="ztrtri.f.html#ZTRTRI.1">ZTRTRI</a>, then U is singular,
</span><span class="comment">*</span><span class="comment">     and the inverse is not computed.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="ZTRTRI.114"></a><a href="ztrtri.f.html#ZTRTRI.1">ZTRTRI</a>( <span class="string">'Upper'</span>, <span class="string">'Non-unit'</span>, N, A, LDA, INFO )
      IF( INFO.GT.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span>      NBMIN = 2
      LDWORK = N
      IF( NB.GT.1 .AND. NB.LT.N ) THEN
         IWS = MAX( LDWORK*NB, 1 )
         IF( LWORK.LT.IWS ) THEN
            NB = LWORK / LDWORK
            NBMIN = MAX( 2, <a name="ILAENV.124"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 2, <span class="string">'<a name="ZGETRI.124"></a><a href="zgetri.f.html#ZGETRI.1">ZGETRI</a>'</span>, <span class="string">' '</span>, N, -1, -1, -1 ) )
         END IF
      ELSE
         IWS = N
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve the equation inv(A)*L = inv(U) for inv(A).
</span><span class="comment">*</span><span class="comment">
</span>      IF( NB.LT.NBMIN .OR. NB.GE.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = N, 1, -1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy current column of L to WORK and replace with zeros.
</span><span class="comment">*</span><span class="comment">
</span>            DO 10 I = J + 1, N
               WORK( I ) = A( I, J )
               A( I, J ) = ZERO
   10       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute current column of inv(A).
</span><span class="comment">*</span><span class="comment">
</span>            IF( J.LT.N )
     $         CALL ZGEMV( <span class="string">'No transpose'</span>, N, N-J, -ONE, A( 1, J+1 ),
     $                     LDA, WORK( J+1 ), 1, ONE, A( 1, J ), 1 )
   20    CONTINUE
      ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use blocked code.
</span><span class="comment">*</span><span class="comment">
</span>         NN = ( ( N-1 ) / NB )*NB + 1
         DO 50 J = NN, 1, -NB
            JB = MIN( NB, N-J+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Copy current block column of L to WORK and replace with
</span><span class="comment">*</span><span class="comment">           zeros.
</span><span class="comment">*</span><span class="comment">
</span>            DO 40 JJ = J, J + JB - 1
               DO 30 I = JJ + 1, N
                  WORK( I+( JJ-J )*LDWORK ) = A( I, JJ )
                  A( I, JJ ) = ZERO
   30          CONTINUE
   40       CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute current block column of inv(A).
</span><span class="comment">*</span><span class="comment">
</span>            IF( J+JB.LE.N )
     $         CALL ZGEMM( <span class="string">'No transpose'</span>, <span class="string">'No transpose'</span>, N, JB,
     $                     N-J-JB+1, -ONE, A( 1, J+JB ), LDA,
     $                     WORK( J+JB ), LDWORK, ONE, A( 1, J ), LDA )
            CALL ZTRSM( <span class="string">'Right'</span>, <span class="string">'Lower'</span>, <span class="string">'No transpose'</span>, <span class="string">'Unit'</span>, N, JB,
     $                  ONE, WORK( J ), LDWORK, A( 1, J ), LDA )
   50    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply column interchanges.
</span><span class="comment">*</span><span class="comment">
</span>      DO 60 J = N - 1, 1, -1
         JP = IPIV( J )
         IF( JP.NE.J )
     $      CALL ZSWAP( N, A( 1, J ), 1, A( 1, JP ), 1 )
   60 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      WORK( 1 ) = IWS
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZGETRI.191"></a><a href="zgetri.f.html#ZGETRI.1">ZGETRI</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>