Return to zgeqpf.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour de lapack.

    1: *> \brief \b ZGEQPF
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download ZGEQPF + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeqpf.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeqpf.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeqpf.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
   22: *
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INFO, LDA, M, N
   25: *       ..
   26: *       .. Array Arguments ..
   27: *       INTEGER            JPVT( * )
   28: *       DOUBLE PRECISION   RWORK( * )
   29: *       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   30: *       ..
   31: *
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> This routine is deprecated and has been replaced by routine ZGEQP3.
   39: *>
   40: *> ZGEQPF computes a QR factorization with column pivoting of a
   41: *> complex M-by-N matrix A: A*P = Q*R.
   42: *> \endverbatim
   43: *
   44: *  Arguments:
   45: *  ==========
   46: *
   47: *> \param[in] M
   48: *> \verbatim
   49: *>          M is INTEGER
   50: *>          The number of rows of the matrix A. M >= 0.
   51: *> \endverbatim
   52: *>
   53: *> \param[in] N
   54: *> \verbatim
   55: *>          N is INTEGER
   56: *>          The number of columns of the matrix A. N >= 0
   57: *> \endverbatim
   58: *>
   59: *> \param[in,out] A
   60: *> \verbatim
   61: *>          A is COMPLEX*16 array, dimension (LDA,N)
   62: *>          On entry, the M-by-N matrix A.
   63: *>          On exit, the upper triangle of the array contains the
   64: *>          min(M,N)-by-N upper triangular matrix R; the elements
   65: *>          below the diagonal, together with the array TAU,
   66: *>          represent the unitary matrix Q as a product of
   67: *>          min(m,n) elementary reflectors.
   68: *> \endverbatim
   69: *>
   70: *> \param[in] LDA
   71: *> \verbatim
   72: *>          LDA is INTEGER
   73: *>          The leading dimension of the array A. LDA >= max(1,M).
   74: *> \endverbatim
   75: *>
   76: *> \param[in,out] JPVT
   77: *> \verbatim
   78: *>          JPVT is INTEGER array, dimension (N)
   79: *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
   80: *>          to the front of A*P (a leading column); if JPVT(i) = 0,
   81: *>          the i-th column of A is a free column.
   82: *>          On exit, if JPVT(i) = k, then the i-th column of A*P
   83: *>          was the k-th column of A.
   84: *> \endverbatim
   85: *>
   86: *> \param[out] TAU
   87: *> \verbatim
   88: *>          TAU is COMPLEX*16 array, dimension (min(M,N))
   89: *>          The scalar factors of the elementary reflectors.
   90: *> \endverbatim
   91: *>
   92: *> \param[out] WORK
   93: *> \verbatim
   94: *>          WORK is COMPLEX*16 array, dimension (N)
   95: *> \endverbatim
   96: *>
   97: *> \param[out] RWORK
   98: *> \verbatim
   99: *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
  100: *> \endverbatim
  101: *>
  102: *> \param[out] INFO
  103: *> \verbatim
  104: *>          INFO is INTEGER
  105: *>          = 0:  successful exit
  106: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  107: *> \endverbatim
  108: *
  109: *  Authors:
  110: *  ========
  111: *
  112: *> \author Univ. of Tennessee
  113: *> \author Univ. of California Berkeley
  114: *> \author Univ. of Colorado Denver
  115: *> \author NAG Ltd.
  116: *
  117: *> \date December 2016
  118: *
  119: *> \ingroup complex16GEcomputational
  120: *
  121: *> \par Further Details:
  122: *  =====================
  123: *>
  124: *> \verbatim
  125: *>
  126: *>  The matrix Q is represented as a product of elementary reflectors
  127: *>
  128: *>     Q = H(1) H(2) . . . H(n)
  129: *>
  130: *>  Each H(i) has the form
  131: *>
  132: *>     H = I - tau * v * v**H
  133: *>
  134: *>  where tau is a complex scalar, and v is a complex vector with
  135: *>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
  136: *>
  137: *>  The matrix P is represented in jpvt as follows: If
  138: *>     jpvt(j) = i
  139: *>  then the jth column of P is the ith canonical unit vector.
  140: *>
  141: *>  Partial column norm updating strategy modified by
  142: *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
  143: *>    University of Zagreb, Croatia.
  144: *>  -- April 2011                                                      --
  145: *>  For more details see LAPACK Working Note 176.
  146: *> \endverbatim
  147: *>
  148: *  =====================================================================
  149:       SUBROUTINE ZGEQPF( M, N, A, LDA, JPVT, TAU, WORK, RWORK, INFO )
  150: *
  151: *  -- LAPACK computational routine (version 3.7.0) --
  152: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  153: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  154: *     December 2016
  155: *
  156: *     .. Scalar Arguments ..
  157:       INTEGER            INFO, LDA, M, N
  158: *     ..
  159: *     .. Array Arguments ..
  160:       INTEGER            JPVT( * )
  161:       DOUBLE PRECISION   RWORK( * )
  162:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
  163: *     ..
  164: *
  165: *  =====================================================================
  166: *
  167: *     .. Parameters ..
  168:       DOUBLE PRECISION   ZERO, ONE
  169:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  170: *     ..
  171: *     .. Local Scalars ..
  172:       INTEGER            I, ITEMP, J, MA, MN, PVT
  173:       DOUBLE PRECISION   TEMP, TEMP2, TOL3Z
  174:       COMPLEX*16         AII
  175: *     ..
  176: *     .. External Subroutines ..
  177:       EXTERNAL           XERBLA, ZGEQR2, ZLARF, ZLARFG, ZSWAP, ZUNM2R
  178: *     ..
  179: *     .. Intrinsic Functions ..
  180:       INTRINSIC          ABS, DCMPLX, DCONJG, MAX, MIN, SQRT
  181: *     ..
  182: *     .. External Functions ..
  183:       INTEGER            IDAMAX
  184:       DOUBLE PRECISION   DLAMCH, DZNRM2
  185:       EXTERNAL           IDAMAX, DLAMCH, DZNRM2
  186: *     ..
  187: *     .. Executable Statements ..
  188: *
  189: *     Test the input arguments
  190: *
  191:       INFO = 0
  192:       IF( M.LT.0 ) THEN
  193:          INFO = -1
  194:       ELSE IF( N.LT.0 ) THEN
  195:          INFO = -2
  196:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  197:          INFO = -4
  198:       END IF
  199:       IF( INFO.NE.0 ) THEN
  200:          CALL XERBLA( 'ZGEQPF', -INFO )
  201:          RETURN
  202:       END IF
  203: *
  204:       MN = MIN( M, N )
  205:       TOL3Z = SQRT(DLAMCH('Epsilon'))
  206: *
  207: *     Move initial columns up front
  208: *
  209:       ITEMP = 1
  210:       DO 10 I = 1, N
  211:          IF( JPVT( I ).NE.0 ) THEN
  212:             IF( I.NE.ITEMP ) THEN
  213:                CALL ZSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
  214:                JPVT( I ) = JPVT( ITEMP )
  215:                JPVT( ITEMP ) = I
  216:             ELSE
  217:                JPVT( I ) = I
  218:             END IF
  219:             ITEMP = ITEMP + 1
  220:          ELSE
  221:             JPVT( I ) = I
  222:          END IF
  223:    10 CONTINUE
  224:       ITEMP = ITEMP - 1
  225: *
  226: *     Compute the QR factorization and update remaining columns
  227: *
  228:       IF( ITEMP.GT.0 ) THEN
  229:          MA = MIN( ITEMP, M )
  230:          CALL ZGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
  231:          IF( MA.LT.N ) THEN
  232:             CALL ZUNM2R( 'Left', 'Conjugate transpose', M, N-MA, MA, A,
  233:      $                   LDA, TAU, A( 1, MA+1 ), LDA, WORK, INFO )
  234:          END IF
  235:       END IF
  236: *
  237:       IF( ITEMP.LT.MN ) THEN
  238: *
  239: *        Initialize partial column norms. The first n elements of
  240: *        work store the exact column norms.
  241: *
  242:          DO 20 I = ITEMP + 1, N
  243:             RWORK( I ) = DZNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
  244:             RWORK( N+I ) = RWORK( I )
  245:    20    CONTINUE
  246: *
  247: *        Compute factorization
  248: *
  249:          DO 40 I = ITEMP + 1, MN
  250: *
  251: *           Determine ith pivot column and swap if necessary
  252: *
  253:             PVT = ( I-1 ) + IDAMAX( N-I+1, RWORK( I ), 1 )
  254: *
  255:             IF( PVT.NE.I ) THEN
  256:                CALL ZSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  257:                ITEMP = JPVT( PVT )
  258:                JPVT( PVT ) = JPVT( I )
  259:                JPVT( I ) = ITEMP
  260:                RWORK( PVT ) = RWORK( I )
  261:                RWORK( N+PVT ) = RWORK( N+I )
  262:             END IF
  263: *
  264: *           Generate elementary reflector H(i)
  265: *
  266:             AII = A( I, I )
  267:             CALL ZLARFG( M-I+1, AII, A( MIN( I+1, M ), I ), 1,
  268:      $                   TAU( I ) )
  269:             A( I, I ) = AII
  270: *
  271:             IF( I.LT.N ) THEN
  272: *
  273: *              Apply H(i) to A(i:m,i+1:n) from the left
  274: *
  275:                AII = A( I, I )
  276:                A( I, I ) = DCMPLX( ONE )
  277:                CALL ZLARF( 'Left', M-I+1, N-I, A( I, I ), 1,
  278:      $                     DCONJG( TAU( I ) ), A( I, I+1 ), LDA, WORK )
  279:                A( I, I ) = AII
  280:             END IF
  281: *
  282: *           Update partial column norms
  283: *
  284:             DO 30 J = I + 1, N
  285:                IF( RWORK( J ).NE.ZERO ) THEN
  286: *
  287: *                 NOTE: The following 4 lines follow from the analysis in
  288: *                 Lapack Working Note 176.
  289: *
  290:                   TEMP = ABS( A( I, J ) ) / RWORK( J )
  291:                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
  292:                   TEMP2 = TEMP*( RWORK( J ) / RWORK( N+J ) )**2
  293:                   IF( TEMP2 .LE. TOL3Z ) THEN
  294:                      IF( M-I.GT.0 ) THEN
  295:                         RWORK( J ) = DZNRM2( M-I, A( I+1, J ), 1 )
  296:                         RWORK( N+J ) = RWORK( J )
  297:                      ELSE
  298:                         RWORK( J ) = ZERO
  299:                         RWORK( N+J ) = ZERO
  300:                      END IF
  301:                   ELSE
  302:                      RWORK( J ) = RWORK( J )*SQRT( TEMP )
  303:                   END IF
  304:                END IF
  305:    30       CONTINUE
  306: *
  307:    40    CONTINUE
  308:       END IF
  309:       RETURN
  310: *
  311: *     End of ZGEQPF
  312: *
  313:       END