zgeqpf.f.html

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

</span>      INFO = 0
      IF( M.LT.0 ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
         INFO = -4
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.121"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="ZGEQPF.121"></a><a href="zgeqpf.f.html#ZGEQPF.1">ZGEQPF</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      MN = MIN( M, N )
      TOL3Z = SQRT(<a name="DLAMCH.126"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>(<span class="string">'Epsilon'</span>))
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Move initial columns up front
</span><span class="comment">*</span><span class="comment">
</span>      ITEMP = 1
      DO 10 I = 1, N
         IF( JPVT( I ).NE.0 ) THEN
            IF( I.NE.ITEMP ) THEN
               CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
               JPVT( I ) = JPVT( ITEMP )
               JPVT( ITEMP ) = I
            ELSE
               JPVT( I ) = I
            END IF
            ITEMP = ITEMP + 1
         ELSE
            JPVT( I ) = I
         END IF
   10 CONTINUE
      ITEMP = ITEMP - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the QR factorization and update remaining columns
</span><span class="comment">*</span><span class="comment">
</span>      IF( ITEMP.GT.0 ) THEN
         MA = MIN( ITEMP, M )
         CALL <a name="ZGEQR2.151"></a><a href="zgeqr2.f.html#ZGEQR2.1">ZGEQR2</a>( M, MA, A, LDA, TAU, WORK, INFO )
         IF( MA.LT.N ) THEN
            CALL <a name="ZUNM2R.153"></a><a href="zunm2r.f.html#ZUNM2R.1">ZUNM2R</a>( <span class="string">'Left'</span>, <span class="string">'Conjugate transpose'</span>, M, N-MA, MA, A,
     $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( ITEMP.LT.MN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Initialize partial column norms. The first n elements of
</span><span class="comment">*</span><span class="comment">        work store the exact column norms.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 I = ITEMP + 1, N
            RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
            RWORK( N+I ) = RWORK( I )
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute factorization
</span><span class="comment">*</span><span class="comment">
</span>         DO 40 I = ITEMP + 1, MN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Determine ith pivot column and swap if necessary
</span><span class="comment">*</span><span class="comment">
</span>            PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
<span class="comment">*</span><span class="comment">
</span>            IF( PVT.NE.I ) THEN
               CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
               ITEMP = JPVT( PVT )
               JPVT( PVT ) = JPVT( I )
               JPVT( I ) = ITEMP
               RWORK( PVT ) = RWORK( I )
               RWORK( N+PVT ) = RWORK( N+I )
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Generate elementary reflector H(i)
</span><span class="comment">*</span><span class="comment">
</span>            AII = A( I, I )
            CALL <a name="ZLARFG.188"></a><a href="zlarfg.f.html#ZLARFG.1">ZLARFG</a>( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
     $                   TAU( I ) )
            A( I, I ) = AII
<span class="comment">*</span><span class="comment">
</span>            IF( I.LT.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Apply H(i) to A(i:m,i+1:n) from the left
</span><span class="comment">*</span><span class="comment">
</span>               AII = A( I, I )
               A( I, I ) = DCMPLX( ONE )
               CALL <a name="ZLARF.198"></a><a href="zlarf.f.html#ZLARF.1">ZLARF</a>( <span class="string">'Left'</span>, M-I+1, N-I, A( I, I ), 1,
     $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
               A( I, I ) = AII
            END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Update partial column norms
</span><span class="comment">*</span><span class="comment">
</span>            DO 30 J = I + 1, N
               IF( RWORK( J ).NE.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 NOTE: The following 4 lines follow from the analysis in
</span><span class="comment">*</span><span class="comment">                 Lapack Working Note 176.
</span><span class="comment">*</span><span class="comment">                 
</span>                  TEMP = ABS( A( I, J ) ) / RWORK( J )
                  TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
                  TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
                  IF( TEMP2 .LE. TOL3Z ) THEN 
                     IF( M-I.GT.0 ) THEN
                        RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
                        RWORK( N+J ) = RWORK( J )
                     ELSE
                        RWORK( J ) = ZERO
                        RWORK( N+J ) = ZERO
                     END IF
                  ELSE
                     RWORK( J ) = RWORK( J )*SQRT( TEMP )
                  END IF
               END IF
   30       CONTINUE
<span class="comment">*</span><span class="comment">
</span>   40    CONTINUE
      END IF
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="ZGEQPF.232"></a><a href="zgeqpf.f.html#ZGEQPF.1">ZGEQPF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>