Return to zlaqp2.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack vers la version 3.2.2.

    1:       SUBROUTINE ZLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
    2:      $                   WORK )
    3: *
    4: *  -- LAPACK auxiliary routine (version 3.2.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     June 2010
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            LDA, M, N, OFFSET
   11: *     ..
   12: *     .. Array Arguments ..
   13:       INTEGER            JPVT( * )
   14:       DOUBLE PRECISION   VN1( * ), VN2( * )
   15:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  ZLAQP2 computes a QR factorization with column pivoting of
   22: *  the block A(OFFSET+1:M,1:N).
   23: *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
   24: *
   25: *  Arguments
   26: *  =========
   27: *
   28: *  M       (input) INTEGER
   29: *          The number of rows of the matrix A. M >= 0.
   30: *
   31: *  N       (input) INTEGER
   32: *          The number of columns of the matrix A. N >= 0.
   33: *
   34: *  OFFSET  (input) INTEGER
   35: *          The number of rows of the matrix A that must be pivoted
   36: *          but no factorized. OFFSET >= 0.
   37: *
   38: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
   39: *          On entry, the M-by-N matrix A.
   40: *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
   41: *          the triangular factor obtained; the elements in block
   42: *          A(OFFSET+1:M,1:N) below the diagonal, together with the
   43: *          array TAU, represent the orthogonal matrix Q as a product of
   44: *          elementary reflectors. Block A(1:OFFSET,1:N) has been
   45: *          accordingly pivoted, but no factorized.
   46: *
   47: *  LDA     (input) INTEGER
   48: *          The leading dimension of the array A. LDA >= max(1,M).
   49: *
   50: *  JPVT    (input/output) INTEGER array, dimension (N)
   51: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
   52: *          to the front of A*P (a leading column); if JPVT(i) = 0,
   53: *          the i-th column of A is a free column.
   54: *          On exit, if JPVT(i) = k, then the i-th column of A*P
   55: *          was the k-th column of A.
   56: *
   57: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
   58: *          The scalar factors of the elementary reflectors.
   59: *
   60: *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)
   61: *          The vector with the partial column norms.
   62: *
   63: *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)
   64: *          The vector with the exact column norms.
   65: *
   66: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
   67: *
   68: *  Further Details
   69: *  ===============
   70: *
   71: *  Based on contributions by
   72: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
   73: *    X. Sun, Computer Science Dept., Duke University, USA
   74: *
   75: *  Partial column norm updating strategy modified by
   76: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
   77: *    University of Zagreb, Croatia.
   78: *     June 2010
   79: *  For more details see LAPACK Working Note 176.
   80: *  =====================================================================
   81: *
   82: *     .. Parameters ..
   83:       DOUBLE PRECISION   ZERO, ONE
   84:       COMPLEX*16         CONE
   85:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0,
   86:      $                   CONE = ( 1.0D+0, 0.0D+0 ) )
   87: *     ..
   88: *     .. Local Scalars ..
   89:       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
   90:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
   91:       COMPLEX*16         AII
   92: *     ..
   93: *     .. External Subroutines ..
   94:       EXTERNAL           ZLARF, ZLARFG, ZSWAP
   95: *     ..
   96: *     .. Intrinsic Functions ..
   97:       INTRINSIC          ABS, DCONJG, MAX, MIN, SQRT
   98: *     ..
   99: *     .. External Functions ..
  100:       INTEGER            IDAMAX
  101:       DOUBLE PRECISION   DLAMCH, DZNRM2
  102:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
  103: *     ..
  104: *     .. Executable Statements ..
  105: *
  106:       MN = MIN( M-OFFSET, N )
  107:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  108: *
  109: *     Compute factorization.
  110: *
  111:       DO 20 I = 1, MN
  112: *
  113:          OFFPI = OFFSET + I
  114: *
  115: *        Determine ith pivot column and swap if necessary.
  116: *
  117:          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
  118: *
  119:          IF( PVT.NE.I ) THEN
  120:             CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  121:             ITEMP = JPVT( PVT )
  122:             JPVT( PVT ) = JPVT( I )
  123:             JPVT( I ) = ITEMP
  124:             VN1( PVT ) = VN1( I )
  125:             VN2( PVT ) = VN2( I )
  126:          END IF
  127: *
  128: *        Generate elementary reflector H(i).
  129: *
  130:          IF( OFFPI.LT.M ) THEN
  131:             CALL ZLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
  132:      $                   TAU( I ) )
  133:          ELSE
  134:             CALL ZLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
  135:          END IF
  136: *
  137:          IF( I.LT.N ) THEN
  138: *
  139: *           Apply H(i)' to A(offset+i:m,i+1:n) from the left.
  140: *
  141:             AII = A( OFFPI, I )
  142:             A( OFFPI, I ) = CONE
  143:             CALL ZLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
  144:      $                  DCONJG( TAU( I ) ), A( OFFPI, I+1 ), LDA,
  145:      $                  WORK( 1 ) )
  146:             A( OFFPI, I ) = AII
  147:          END IF
  148: *
  149: *        Update partial column norms.
  150: *
  151:          DO 10 J = I + 1, N
  152:             IF( VN1( J ).NE.ZERO ) THEN
  153: *
  154: *              NOTE: The following 4 lines follow from the analysis in
  155: *              Lapack Working Note 176.
  156: *
  157:                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
  158:                TEMP = MAX( TEMP, ZERO )
  159:                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
  160:                IF( TEMP2 .LE. TOL3Z ) THEN
  161:                   IF( OFFPI.LT.M ) THEN
  162:                      VN1( J ) = DZNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
  163:                      VN2( J ) = VN1( J )
  164:                   ELSE
  165:                      VN1( J ) = ZERO
  166:                      VN2( J ) = ZERO
  167:                   END IF
  168:                ELSE
  169:                   VN1( J ) = VN1( J )*SQRT( TEMP )
  170:                END IF
  171:             END IF
  172:    10    CONTINUE
  173: *
  174:    20 CONTINUE
  175: *
  176:       RETURN
  177: *
  178: *     End of ZLAQP2
  179: *
  180:       END