sgeqpf.f.html

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

         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.117"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGEQPF.117"></a><a href="sgeqpf.f.html#SGEQPF.1">SGEQPF</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span>      MN = MIN( M, N )
      TOL3Z = SQRT(<a name="SLAMCH.122"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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 SSWAP( 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="SGEQR2.147"></a><a href="sgeqr2.f.html#SGEQR2.1">SGEQR2</a>( M, MA, A, LDA, TAU, WORK, INFO )
         IF( MA.LT.N ) THEN
            CALL <a name="SORM2R.149"></a><a href="sorm2r.f.html#SORM2R.1">SORM2R</a>( <span class="string">'Left'</span>, <span class="string">'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
            WORK( I ) = SNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
            WORK( N+I ) = WORK( 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 ) + ISAMAX( N-I+1, WORK( I ), 1 )
<span class="comment">*</span><span class="comment">
</span>            IF( PVT.NE.I ) THEN
               CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
               ITEMP = JPVT( PVT )
               JPVT( PVT ) = JPVT( I )
               JPVT( I ) = ITEMP
               WORK( PVT ) = WORK( I )
               WORK( N+PVT ) = WORK( 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>            IF( I.LT.M ) THEN
               CALL <a name="SLARFG.184"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
            ELSE
               CALL <a name="SLARFG.186"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
            END IF
<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 ) = ONE
               CALL <a name="SLARF.195"></a><a href="slarf.f.html#SLARF.1">SLARF</a>( <span class="string">'LEFT'</span>, M-I+1, N-I, A( I, I ), 1, TAU( I ),
     $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
               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( WORK( 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 ) ) / WORK( J )
                  TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
                  TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2
                  IF( TEMP2 .LE. TOL3Z ) THEN 
                     IF( M-I.GT.0 ) THEN
                        WORK( J ) = SNRM2( M-I, A( I+1, J ), 1 )
                        WORK( N+J ) = WORK( J )
                     ELSE
                        WORK( J ) = ZERO
                        WORK( N+J ) = ZERO
                     END IF
                  ELSE
                     WORK( J ) = WORK( 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="SGEQPF.229"></a><a href="sgeqpf.f.html#SGEQPF.1">SGEQPF</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>