Return to zgeqp3.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, RWORK,
    2:      $                   INFO )
    3: *
    4: *  -- LAPACK routine (version 3.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: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       INTEGER            INFO, LDA, LWORK, M, N
   11: *     ..
   12: *     .. Array Arguments ..
   13:       INTEGER            JPVT( * )
   14:       DOUBLE PRECISION   RWORK( * )
   15:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   16: *     ..
   17: *
   18: *  Purpose
   19: *  =======
   20: *
   21: *  ZGEQP3 computes a QR factorization with column pivoting of a
   22: *  matrix A:  A*P = Q*R  using Level 3 BLAS.
   23: *
   24: *  Arguments
   25: *  =========
   26: *
   27: *  M       (input) INTEGER
   28: *          The number of rows of the matrix A. M >= 0.
   29: *
   30: *  N       (input) INTEGER
   31: *          The number of columns of the matrix A.  N >= 0.
   32: *
   33: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
   34: *          On entry, the M-by-N matrix A.
   35: *          On exit, the upper triangle of the array contains the
   36: *          min(M,N)-by-N upper trapezoidal matrix R; the elements below
   37: *          the diagonal, together with the array TAU, represent the
   38: *          unitary matrix Q as a product of min(M,N) elementary
   39: *          reflectors.
   40: *
   41: *  LDA     (input) INTEGER
   42: *          The leading dimension of the array A. LDA >= max(1,M).
   43: *
   44: *  JPVT    (input/output) INTEGER array, dimension (N)
   45: *          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
   46: *          to the front of A*P (a leading column); if JPVT(J)=0,
   47: *          the J-th column of A is a free column.
   48: *          On exit, if JPVT(J)=K, then the J-th column of A*P was the
   49: *          the K-th column of A.
   50: *
   51: *  TAU     (output) COMPLEX*16 array, dimension (min(M,N))
   52: *          The scalar factors of the elementary reflectors.
   53: *
   54: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
   55: *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
   56: *
   57: *  LWORK   (input) INTEGER
   58: *          The dimension of the array WORK. LWORK >= N+1.
   59: *          For optimal performance LWORK >= ( N+1 )*NB, where NB
   60: *          is the optimal blocksize.
   61: *
   62: *          If LWORK = -1, then a workspace query is assumed; the routine
   63: *          only calculates the optimal size of the WORK array, returns
   64: *          this value as the first entry of the WORK array, and no error
   65: *          message related to LWORK is issued by XERBLA.
   66: *
   67: *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
   68: *
   69: *  INFO    (output) INTEGER
   70: *          = 0: successful exit.
   71: *          < 0: if INFO = -i, the i-th argument had an illegal value.
   72: *
   73: *  Further Details
   74: *  ===============
   75: *
   76: *  The matrix Q is represented as a product of elementary reflectors
   77: *
   78: *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
   79: *
   80: *  Each H(i) has the form
   81: *
   82: *     H(i) = I - tau * v * v'
   83: *
   84: *  where tau is a real/complex scalar, and v is a real/complex vector
   85: *  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
   86: *  A(i+1:m,i), and tau in TAU(i).
   87: *
   88: *  Based on contributions by
   89: *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
   90: *    X. Sun, Computer Science Dept., Duke University, USA
   91: *
   92: *  =====================================================================
   93: *
   94: *     .. Parameters ..
   95:       INTEGER            INB, INBMIN, IXOVER
   96:       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
   97: *     ..
   98: *     .. Local Scalars ..
   99:       LOGICAL            LQUERY
  100:       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
  101:      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
  102: *     ..
  103: *     .. External Subroutines ..
  104:       EXTERNAL           XERBLA, ZGEQRF, ZLAQP2, ZLAQPS, ZSWAP, ZUNMQR
  105: *     ..
  106: *     .. External Functions ..
  107:       INTEGER            ILAENV
  108:       DOUBLE PRECISION   DZNRM2
  109:       EXTERNAL           ILAENV, DZNRM2
  110: *     ..
  111: *     .. Intrinsic Functions ..
  112:       INTRINSIC          INT, MAX, MIN
  113: *     ..
  114: *     .. Executable Statements ..
  115: *
  116: *     Test input arguments
  117: *     ====================
  118: *
  119:       INFO = 0
  120:       LQUERY = ( LWORK.EQ.-1 )
  121:       IF( M.LT.0 ) THEN
  122:          INFO = -1
  123:       ELSE IF( N.LT.0 ) THEN
  124:          INFO = -2
  125:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  126:          INFO = -4
  127:       END IF
  128: *
  129:       IF( INFO.EQ.0 ) THEN
  130:          MINMN = MIN( M, N )
  131:          IF( MINMN.EQ.0 ) THEN
  132:             IWS = 1
  133:             LWKOPT = 1
  134:          ELSE
  135:             IWS = N + 1
  136:             NB = ILAENV( INB, 'ZGEQRF', ' ', M, N, -1, -1 )
  137:             LWKOPT = ( N + 1 )*NB
  138:          END IF
  139:          WORK( 1 ) = LWKOPT
  140: *
  141:          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
  142:             INFO = -8
  143:          END IF
  144:       END IF
  145: *
  146:       IF( INFO.NE.0 ) THEN
  147:          CALL XERBLA( 'ZGEQP3', -INFO )
  148:          RETURN
  149:       ELSE IF( LQUERY ) THEN
  150:          RETURN
  151:       END IF
  152: *
  153: *     Quick return if possible.
  154: *
  155:       IF( MINMN.EQ.0 ) THEN
  156:          RETURN
  157:       END IF
  158: *
  159: *     Move initial columns up front.
  160: *
  161:       NFXD = 1
  162:       DO 10 J = 1, N
  163:          IF( JPVT( J ).NE.0 ) THEN
  164:             IF( J.NE.NFXD ) THEN
  165:                CALL ZSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
  166:                JPVT( J ) = JPVT( NFXD )
  167:                JPVT( NFXD ) = J
  168:             ELSE
  169:                JPVT( J ) = J
  170:             END IF
  171:             NFXD = NFXD + 1
  172:          ELSE
  173:             JPVT( J ) = J
  174:          END IF
  175:    10 CONTINUE
  176:       NFXD = NFXD - 1
  177: *
  178: *     Factorize fixed columns
  179: *     =======================
  180: *
  181: *     Compute the QR factorization of fixed columns and update
  182: *     remaining columns.
  183: *
  184:       IF( NFXD.GT.0 ) THEN
  185:          NA = MIN( M, NFXD )
  186: *CC      CALL ZGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
  187:          CALL ZGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
  188:          IWS = MAX( IWS, INT( WORK( 1 ) ) )
  189:          IF( NA.LT.N ) THEN
  190: *CC         CALL ZUNM2R( 'Left', 'Conjugate Transpose', M, N-NA,
  191: *CC  $                   NA, A, LDA, TAU, A( 1, NA+1 ), LDA, WORK,
  192: *CC  $                   INFO )
  193:             CALL ZUNMQR( 'Left', 'Conjugate Transpose', M, N-NA, NA, A,
  194:      $                   LDA, TAU, A( 1, NA+1 ), LDA, WORK, LWORK,
  195:      $                   INFO )
  196:             IWS = MAX( IWS, INT( WORK( 1 ) ) )
  197:          END IF
  198:       END IF
  199: *
  200: *     Factorize free columns
  201: *     ======================
  202: *
  203:       IF( NFXD.LT.MINMN ) THEN
  204: *
  205:          SM = M - NFXD
  206:          SN = N - NFXD
  207:          SMINMN = MINMN - NFXD
  208: *
  209: *        Determine the block size.
  210: *
  211:          NB = ILAENV( INB, 'ZGEQRF', ' ', SM, SN, -1, -1 )
  212:          NBMIN = 2
  213:          NX = 0
  214: *
  215:          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
  216: *
  217: *           Determine when to cross over from blocked to unblocked code.
  218: *
  219:             NX = MAX( 0, ILAENV( IXOVER, 'ZGEQRF', ' ', SM, SN, -1,
  220:      $           -1 ) )
  221: *
  222: *
  223:             IF( NX.LT.SMINMN ) THEN
  224: *
  225: *              Determine if workspace is large enough for blocked code.
  226: *
  227:                MINWS = ( SN+1 )*NB
  228:                IWS = MAX( IWS, MINWS )
  229:                IF( LWORK.LT.MINWS ) THEN
  230: *
  231: *                 Not enough workspace to use optimal NB: Reduce NB and
  232: *                 determine the minimum value of NB.
  233: *
  234:                   NB = LWORK / ( SN+1 )
  235:                   NBMIN = MAX( 2, ILAENV( INBMIN, 'ZGEQRF', ' ', SM, SN,
  236:      $                    -1, -1 ) )
  237: *
  238: *
  239:                END IF
  240:             END IF
  241:          END IF
  242: *
  243: *        Initialize partial column norms. The first N elements of work
  244: *        store the exact column norms.
  245: *
  246:          DO 20 J = NFXD + 1, N
  247:             RWORK( J ) = DZNRM2( SM, A( NFXD+1, J ), 1 )
  248:             RWORK( N+J ) = RWORK( J )
  249:    20    CONTINUE
  250: *
  251:          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
  252:      $       ( NX.LT.SMINMN ) ) THEN
  253: *
  254: *           Use blocked code initially.
  255: *
  256:             J = NFXD + 1
  257: *
  258: *           Compute factorization: while loop.
  259: *
  260: *
  261:             TOPBMN = MINMN - NX
  262:    30       CONTINUE
  263:             IF( J.LE.TOPBMN ) THEN
  264:                JB = MIN( NB, TOPBMN-J+1 )
  265: *
  266: *              Factorize JB columns among columns J:N.
  267: *
  268:                CALL ZLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
  269:      $                      JPVT( J ), TAU( J ), RWORK( J ),
  270:      $                      RWORK( N+J ), WORK( 1 ), WORK( JB+1 ),
  271:      $                      N-J+1 )
  272: *
  273:                J = J + FJB
  274:                GO TO 30
  275:             END IF
  276:          ELSE
  277:             J = NFXD + 1
  278:          END IF
  279: *
  280: *        Use unblocked code to factor the last or only block.
  281: *
  282: *
  283:          IF( J.LE.MINMN )
  284:      $      CALL ZLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
  285:      $                   TAU( J ), RWORK( J ), RWORK( N+J ), WORK( 1 ) )
  286: *
  287:       END IF
  288: *
  289:       WORK( 1 ) = IWS
  290:       RETURN
  291: *
  292: *     End of ZGEQP3
  293: *
  294:       END