sgeqp3.f.html

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

<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.139"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGEQP3.139"></a><a href="sgeqp3.f.html#SGEQP3.1">SGEQP3</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( MINMN.EQ.0 ) THEN
         RETURN
      END IF
<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>      NFXD = 1
      DO 10 J = 1, N
         IF( JPVT( J ).NE.0 ) THEN
            IF( J.NE.NFXD ) THEN
               CALL SSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
               JPVT( J ) = JPVT( NFXD )
               JPVT( NFXD ) = J
            ELSE
               JPVT( J ) = J
            END IF
            NFXD = NFXD + 1
         ELSE
            JPVT( J ) = J
         END IF
   10 CONTINUE
      NFXD = NFXD - 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Factorize fixed columns
</span><span class="comment">*</span><span class="comment">     =======================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the QR factorization of fixed columns and update
</span><span class="comment">*</span><span class="comment">     remaining columns.
</span><span class="comment">*</span><span class="comment">
</span>      IF( NFXD.GT.0 ) THEN
         NA = MIN( M, NFXD )
<span class="comment">*</span><span class="comment">CC      CALL <a name="SGEQR2.178"></a><a href="sgeqr2.f.html#SGEQR2.1">SGEQR2</a>( M, NA, A, LDA, TAU, WORK, INFO )
</span>         CALL <a name="SGEQRF.179"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
         IWS = MAX( IWS, INT( WORK( 1 ) ) )
         IF( NA.LT.N ) THEN
<span class="comment">*</span><span class="comment">CC         CALL <a name="SORM2R.182"></a><a href="sorm2r.f.html#SORM2R.1">SORM2R</a>( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
</span><span class="comment">*</span><span class="comment">CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
</span>            CALL <a name="SORMQR.184"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, M, N-NA, NA, A, LDA, TAU,
     $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
            IWS = MAX( IWS, INT( WORK( 1 ) ) )
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Factorize free columns
</span><span class="comment">*</span><span class="comment">     ======================
</span><span class="comment">*</span><span class="comment">
</span>      IF( NFXD.LT.MINMN ) THEN
<span class="comment">*</span><span class="comment">
</span>         SM = M - NFXD
         SN = N - NFXD
         SMINMN = MINMN - NFXD
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine the block size.
</span><span class="comment">*</span><span class="comment">
</span>         NB = <a name="ILAENV.201"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( INB, <span class="string">'<a name="SGEQRF.201"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>'</span>, <span class="string">' '</span>, SM, SN, -1, -1 )
         NBMIN = 2
         NX = 0
<span class="comment">*</span><span class="comment">
</span>         IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Determine when to cross over from blocked to unblocked code.
</span><span class="comment">*</span><span class="comment">
</span>            NX = MAX( 0, <a name="ILAENV.209"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( IXOVER, <span class="string">'<a name="SGEQRF.209"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>'</span>, <span class="string">' '</span>, SM, SN, -1,
     $           -1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>            IF( NX.LT.SMINMN ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Determine if workspace is large enough for blocked code.
</span><span class="comment">*</span><span class="comment">
</span>               MINWS = 2*SN + ( SN+1 )*NB
               IWS = MAX( IWS, MINWS )
               IF( LWORK.LT.MINWS ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Not enough workspace to use optimal NB: Reduce NB and
</span><span class="comment">*</span><span class="comment">                 determine the minimum value of NB.
</span><span class="comment">*</span><span class="comment">
</span>                  NB = ( LWORK-2*SN ) / ( SN+1 )
                  NBMIN = MAX( 2, <a name="ILAENV.225"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( INBMIN, <span class="string">'<a name="SGEQRF.225"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>'</span>, <span class="string">' '</span>, SM, SN,
     $                    -1, -1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>               END IF
            END IF
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Initialize partial column norms. The first N elements of work
</span><span class="comment">*</span><span class="comment">        store the exact column norms.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = NFXD + 1, N
            WORK( J ) = SNRM2( SM, A( NFXD+1, J ), 1 )
            WORK( N+J ) = WORK( J )
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
     $       ( NX.LT.SMINMN ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Use blocked code initially.
</span><span class="comment">*</span><span class="comment">
</span>            J = NFXD + 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Compute factorization: while loop.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>            TOPBMN = MINMN - NX
   30       CONTINUE
            IF( J.LE.TOPBMN ) THEN
               JB = MIN( NB, TOPBMN-J+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Factorize JB columns among columns J:N.
</span><span class="comment">*</span><span class="comment">
</span>               CALL <a name="SLAQPS.258"></a><a href="slaqps.f.html#SLAQPS.1">SLAQPS</a>( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
     $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
     $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
<span class="comment">*</span><span class="comment">
</span>               J = J + FJB
               GO TO 30
            END IF
         ELSE
            J = NFXD + 1
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Use unblocked code to factor the last or only block.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">
</span>         IF( J.LE.MINMN )
     $      CALL <a name="SLAQP2.273"></a><a href="slaqp2.f.html#SLAQP2.1">SLAQP2</a>( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
     $                   TAU( J ), WORK( J ), WORK( N+J ),
     $                   WORK( 2*N+1 ) )
<span class="comment">*</span><span class="comment">
</span>      END IF
<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="SGEQP3.282"></a><a href="sgeqp3.f.html#SGEQP3.1">SGEQP3</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>