slaqps.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>   10 CONTINUE
      IF( ( K.LT.NB ) .AND. ( LSTICC.EQ.0 ) ) THEN
         K = K + 1
         RK = OFFSET + K
<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 = ( K-1 ) + ISAMAX( N-K+1, VN1( K ), 1 )
         IF( PVT.NE.K ) THEN
            CALL SSWAP( M, A( 1, PVT ), 1, A( 1, K ), 1 )
            CALL SSWAP( K-1, F( PVT, 1 ), LDF, F( K, 1 ), LDF )
            ITEMP = JPVT( PVT )
            JPVT( PVT ) = JPVT( K )
            JPVT( K ) = ITEMP
            VN1( PVT ) = VN1( K )
            VN2( PVT ) = VN2( K )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Apply previous Householder reflectors to column K:
</span><span class="comment">*</span><span class="comment">        A(RK:M,K) := A(RK:M,K) - A(RK:M,1:K-1)*F(K,1:K-1)'.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.1 ) THEN
            CALL SGEMV( <span class="string">'No transpose'</span>, M-RK+1, K-1, -ONE, A( RK, 1 ),
     $                  LDA, F( K, 1 ), LDF, ONE, A( RK, K ), 1 )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Generate elementary reflector H(k).
</span><span class="comment">*</span><span class="comment">
</span>         IF( RK.LT.M ) THEN
            CALL <a name="SLARFG.154"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( M-RK+1, A( RK, K ), A( RK+1, K ), 1, TAU( K ) )
         ELSE
            CALL <a name="SLARFG.156"></a><a href="slarfg.f.html#SLARFG.1">SLARFG</a>( 1, A( RK, K ), A( RK, K ), 1, TAU( K ) )
         END IF
<span class="comment">*</span><span class="comment">
</span>         AKK = A( RK, K )
         A( RK, K ) = ONE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute Kth column of F:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute  F(K+1:N,K) := tau(K)*A(RK:M,K+1:N)'*A(RK:M,K).
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.N ) THEN
            CALL SGEMV( <span class="string">'Transpose'</span>, M-RK+1, N-K, TAU( K ),
     $                  A( RK, K+1 ), LDA, A( RK, K ), 1, ZERO,
     $                  F( K+1, K ), 1 )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Padding F(1:K,K) with zeros.
</span><span class="comment">*</span><span class="comment">
</span>         DO 20 J = 1, K
            F( J, K ) = ZERO
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Incremental updating of F:
</span><span class="comment">*</span><span class="comment">        F(1:N,K) := F(1:N,K) - tau(K)*F(1:N,1:K-1)*A(RK:M,1:K-1)'
</span><span class="comment">*</span><span class="comment">                    *A(RK:M,K).
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.1 ) THEN
            CALL SGEMV( <span class="string">'Transpose'</span>, M-RK+1, K-1, -TAU( K ), A( RK, 1 ),
     $                  LDA, A( RK, K ), 1, ZERO, AUXV( 1 ), 1 )
<span class="comment">*</span><span class="comment">
</span>            CALL SGEMV( <span class="string">'No transpose'</span>, N, K-1, ONE, F( 1, 1 ), LDF,
     $                  AUXV( 1 ), 1, ONE, F( 1, K ), 1 )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Update the current row of A:
</span><span class="comment">*</span><span class="comment">        A(RK,K+1:N) := A(RK,K+1:N) - A(RK,1:K)*F(K+1:N,1:K)'.
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.N ) THEN
            CALL SGEMV( <span class="string">'No transpose'</span>, N-K, K, -ONE, F( K+1, 1 ), LDF,
     $                  A( RK, 1 ), LDA, ONE, A( RK, K+1 ), LDA )
         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>         IF( RK.LT.LASTRK ) THEN
            DO 30 J = K + 1, N
               IF( VN1( 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( RK, J ) ) / VN1( J )
                  TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
                  TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
                  IF( TEMP2 .LE. TOL3Z ) THEN
                     VN2( J ) = REAL( LSTICC )
                     LSTICC = J
                  ELSE
                     VN1( J ) = VN1( J )*SQRT( TEMP )
                  END IF
               END IF
   30       CONTINUE
         END IF
<span class="comment">*</span><span class="comment">
</span>         A( RK, K ) = AKK
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        End of while loop.
</span><span class="comment">*</span><span class="comment">
</span>         GO TO 10
      END IF
      KB = K
      RK = OFFSET + KB
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Apply the block reflector to the rest of the matrix:
</span><span class="comment">*</span><span class="comment">     A(OFFSET+KB+1:M,KB+1:N) := A(OFFSET+KB+1:M,KB+1:N) -
</span><span class="comment">*</span><span class="comment">                         A(OFFSET+KB+1:M,1:KB)*F(KB+1:N,1:KB)'.
</span><span class="comment">*</span><span class="comment">
</span>      IF( KB.LT.MIN( N, M-OFFSET ) ) THEN
         CALL SGEMM( <span class="string">'No transpose'</span>, <span class="string">'Transpose'</span>, M-RK, N-KB, KB, -ONE,
     $               A( RK+1, 1 ), LDA, F( KB+1, 1 ), LDF, ONE,
     $               A( RK+1, KB+1 ), LDA )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Recomputation of difficult columns.
</span><span class="comment">*</span><span class="comment">
</span>   40 CONTINUE
      IF( LSTICC.GT.0 ) THEN
         ITEMP = NINT( VN2( LSTICC ) )
         VN1( LSTICC ) = SNRM2( M-RK, A( RK+1, LSTICC ), 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        NOTE: The computation of VN1( LSTICC ) relies on the fact that 
</span><span class="comment">*</span><span class="comment">        SNRM2 does not fail on vectors with norm below the value of
</span><span class="comment">*</span><span class="comment">        SQRT(<a name="DLAMCH.248"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>('S')) 
</span><span class="comment">*</span><span class="comment">
</span>         VN2( LSTICC ) = VN1( LSTICC )
         LSTICC = ITEMP
         GO TO 40
      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="SLAQPS.257"></a><a href="slaqps.f.html#SLAQPS.1">SLAQPS</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>