Return to zgeqpf.f CVS log

Up to [local] / rpl / lapack / lapack

Initial revision

    1:       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
    2: *
    3: *  -- LAPACK deprecated driver routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, LDA, M, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       INTEGER            JPVT( * )
   13:       DOUBLE PRECISION   RWORK( * )
   14:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  This routine is deprecated and has been replaced by routine ZGEQP3.
   21: *
   22: *  ZGEQPF computes a QR factorization with column pivoting of a
   23: *  complex M-by-N matrix A: A*P = Q*R.
   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: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
   35: *          On entry, the M-by-N matrix A.
   36: *          On exit, the upper triangle of the array contains the
   37: *          min(M,N)-by-N upper triangular matrix R; the elements
   38: *          below the diagonal, together with the array TAU,
   39: *          represent the unitary matrix Q as a product of
   40: *          min(m,n) elementary reflectors.
   41: *
   42: *  LDA     (input) INTEGER
   43: *          The leading dimension of the array A. LDA >= max(1,M).
   44: *
   45: *  JPVT    (input/output) INTEGER array, dimension (N)
   46: *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
   47: *          to the front of A*P (a leading column); if JPVT(i) = 0,
   48: *          the i-th column of A is a free column.
   49: *          On exit, if JPVT(i) = k, then the i-th column of A*P
   50: *          was the k-th column of A.
   51: *
   52: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
   53: *          The scalar factors of the elementary reflectors.
   54: *
   55: *  WORK    (workspace) COMPLEX*16 array, dimension (N)
   56: *
   57: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
   58: *
   59: *  INFO    (output) INTEGER
   60: *          = 0:  successful exit
   61: *          < 0:  if INFO = -i, the i-th argument had an illegal value
   62: *
   63: *  Further Details
   64: *  ===============
   65: *
   66: *  The matrix Q is represented as a product of elementary reflectors
   67: *
   68: *     Q = H(1) H(2) . . . H(n)
   69: *
   70: *  Each H(i) has the form
   71: *
   72: *     H = I - tau * v * v'
   73: *
   74: *  where tau is a complex scalar, and v is a complex vector with
   75: *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
   76: *
   77: *  The matrix P is represented in jpvt as follows: If
   78: *     jpvt(j) = i
   79: *  then the jth column of P is the ith canonical unit vector.
   80: *
   81: *  Partial column norm updating strategy modified by
   82: *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
   83: *    University of Zagreb, Croatia.
   84: *    June 2006.
   85: *  For more details see LAPACK Working Note 176.
   86: *
   87: *  =====================================================================
   88: *
   89: *     .. Parameters ..
   90:       DOUBLE PRECISION   ZERO, ONE
   91:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
   92: *     ..
   93: *     .. Local Scalars ..
   94:       INTEGER            I, ITEMP, J, MA, MN, PVT
   95:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
   96:       COMPLEX*16         AII
   97: *     ..
   98: *     .. External Subroutines ..
   99:       EXTERNAL           XERBLA, ZGEQR2, ZLARF, ZLARFP, ZSWAP, ZUNM2R
  100: *     ..
  101: *     .. Intrinsic Functions ..
  102:       INTRINSIC          ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
  103: *     ..
  104: *     .. External Functions ..
  105:       INTEGER            IDAMAX
  106:       DOUBLE PRECISION   DLAMCH, DZNRM2
  107:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
  108: *     ..
  109: *     .. Executable Statements ..
  110: *
  111: *     Test the input arguments
  112: *
  113:       INFO = 0
  114:       IF( M.LT.0 ) THEN
  115:          INFO = -1
  116:       ELSE IF( N.LT.0 ) THEN
  117:          INFO = -2
  118:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  119:          INFO = -4
  120:       END IF
  121:       IF( INFO.NE.0 ) THEN
  122:          CALL XERBLA( 'ZGEQPF', -INFO )
  123:          RETURN
  124:       END IF
  125: *
  126:       MN = MIN( M, N )
  127:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  128: *
  129: *     Move initial columns up front
  130: *
  131:       ITEMP = 1
  132:       DO 10 I = 1, N
  133:          IF( JPVT( I ).NE.0 ) THEN
  134:             IF( I.NE.ITEMP ) THEN
  135:                CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
  136:                JPVT( I ) = JPVT( ITEMP )
  137:                JPVT( ITEMP ) = I
  138:             ELSE
  139:                JPVT( I ) = I
  140:             END IF
  141:             ITEMP = ITEMP + 1
  142:          ELSE
  143:             JPVT( I ) = I
  144:          END IF
  145:    10 CONTINUE
  146:       ITEMP = ITEMP - 1
  147: *
  148: *     Compute the QR factorization and update remaining columns
  149: *
  150:       IF( ITEMP.GT.0 ) THEN
  151:          MA = MIN( ITEMP, M )
  152:          CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
  153:          IF( MA.LT.N ) THEN
  154:             CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
  155:      $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
  156:          END IF
  157:       END IF
  158: *
  159:       IF( ITEMP.LT.MN ) THEN
  160: *
  161: *        Initialize partial column norms. The first n elements of
  162: *        work store the exact column norms.
  163: *
  164:          DO 20 I = ITEMP + 1, N
  165:             RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
  166:             RWORK( N+I ) = RWORK( I )
  167:    20    CONTINUE
  168: *
  169: *        Compute factorization
  170: *
  171:          DO 40 I = ITEMP + 1, MN
  172: *
  173: *           Determine ith pivot column and swap if necessary
  174: *
  175:             PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
  176: *
  177:             IF( PVT.NE.I ) THEN
  178:                CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  179:                ITEMP = JPVT( PVT )
  180:                JPVT( PVT ) = JPVT( I )
  181:                JPVT( I ) = ITEMP
  182:                RWORK( PVT ) = RWORK( I )
  183:                RWORK( N+PVT ) = RWORK( N+I )
  184:             END IF
  185: *
  186: *           Generate elementary reflector H(i)
  187: *
  188:             AII = A( I, I )
  189:             CALL ZLARFP( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
  190:      $                   TAU( I ) )
  191:             A( I, I ) = AII
  192: *
  193:             IF( I.LT.N ) THEN
  194: *
  195: *              Apply H(i) to A(i:m,i+1:n) from the left
  196: *
  197:                AII = A( I, I )
  198:                A( I, I ) = DCMPLX( ONE )
  199:                CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
  200:      $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
  201:                A( I, I ) = AII
  202:             END IF
  203: *
  204: *           Update partial column norms
  205: *
  206:             DO 30 J = I + 1, N
  207:                IF( RWORK( J ).NE.ZERO ) THEN
  208: *
  209: *                 NOTE: The following 4 lines follow from the analysis in
  210: *                 Lapack Working Note 176.
  211: *                 
  212:                   TEMP = ABS( A( I, J ) ) / RWORK( J )
  213:                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
  214:                   TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
  215:                   IF( TEMP2 .LE. TOL3Z ) THEN 
  216:                      IF( M-I.GT.0 ) THEN
  217:                         RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
  218:                         RWORK( N+J ) = RWORK( J )
  219:                      ELSE
  220:                         RWORK( J ) = ZERO
  221:                         RWORK( N+J ) = ZERO
  222:                      END IF
  223:                   ELSE
  224:                      RWORK( J ) = RWORK( J )*SQRT( TEMP )
  225:                   END IF
  226:                END IF
  227:    30       CONTINUE
  228: *
  229:    40    CONTINUE
  230:       END IF
  231:       RETURN
  232: *
  233: *     End of ZGEQPF
  234: *
  235:       END