Annotation of rpl/lapack/lapack/dgeqp3.f, revision 1.14

1.9       bertrand    1: *> \brief \b DGEQP3
                      2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download DGEQP3 + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeqp3.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeqp3.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeqp3.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
                     22: * 
                     23: *       .. Scalar Arguments ..
                     24: *       INTEGER            INFO, LDA, LWORK, M, N
                     25: *       ..
                     26: *       .. Array Arguments ..
                     27: *       INTEGER            JPVT( * )
                     28: *       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
                     29: *       ..
                     30: *  
                     31: *
                     32: *> \par Purpose:
                     33: *  =============
                     34: *>
                     35: *> \verbatim
                     36: *>
                     37: *> DGEQP3 computes a QR factorization with column pivoting of a
                     38: *> matrix A:  A*P = Q*R  using Level 3 BLAS.
                     39: *> \endverbatim
                     40: *
                     41: *  Arguments:
                     42: *  ==========
                     43: *
                     44: *> \param[in] M
                     45: *> \verbatim
                     46: *>          M is INTEGER
                     47: *>          The number of rows of the matrix A. M >= 0.
                     48: *> \endverbatim
                     49: *>
                     50: *> \param[in] N
                     51: *> \verbatim
                     52: *>          N is INTEGER
                     53: *>          The number of columns of the matrix A.  N >= 0.
                     54: *> \endverbatim
                     55: *>
                     56: *> \param[in,out] A
                     57: *> \verbatim
                     58: *>          A is DOUBLE PRECISION array, dimension (LDA,N)
                     59: *>          On entry, the M-by-N matrix A.
                     60: *>          On exit, the upper triangle of the array contains the
                     61: *>          min(M,N)-by-N upper trapezoidal matrix R; the elements below
                     62: *>          the diagonal, together with the array TAU, represent the
                     63: *>          orthogonal matrix Q as a product of min(M,N) elementary
                     64: *>          reflectors.
                     65: *> \endverbatim
                     66: *>
                     67: *> \param[in] LDA
                     68: *> \verbatim
                     69: *>          LDA is INTEGER
                     70: *>          The leading dimension of the array A. LDA >= max(1,M).
                     71: *> \endverbatim
                     72: *>
                     73: *> \param[in,out] JPVT
                     74: *> \verbatim
                     75: *>          JPVT is INTEGER array, dimension (N)
                     76: *>          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
                     77: *>          to the front of A*P (a leading column); if JPVT(J)=0,
                     78: *>          the J-th column of A is a free column.
                     79: *>          On exit, if JPVT(J)=K, then the J-th column of A*P was the
                     80: *>          the K-th column of A.
                     81: *> \endverbatim
                     82: *>
                     83: *> \param[out] TAU
                     84: *> \verbatim
                     85: *>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
                     86: *>          The scalar factors of the elementary reflectors.
                     87: *> \endverbatim
                     88: *>
                     89: *> \param[out] WORK
                     90: *> \verbatim
                     91: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
                     92: *>          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
                     93: *> \endverbatim
                     94: *>
                     95: *> \param[in] LWORK
                     96: *> \verbatim
                     97: *>          LWORK is INTEGER
                     98: *>          The dimension of the array WORK. LWORK >= 3*N+1.
                     99: *>          For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
                    100: *>          is the optimal blocksize.
                    101: *>
                    102: *>          If LWORK = -1, then a workspace query is assumed; the routine
                    103: *>          only calculates the optimal size of the WORK array, returns
                    104: *>          this value as the first entry of the WORK array, and no error
                    105: *>          message related to LWORK is issued by XERBLA.
                    106: *> \endverbatim
                    107: *>
                    108: *> \param[out] INFO
                    109: *> \verbatim
                    110: *>          INFO is INTEGER
                    111: *>          = 0: successful exit.
                    112: *>          < 0: if INFO = -i, the i-th argument had an illegal value.
                    113: *> \endverbatim
                    114: *
                    115: *  Authors:
                    116: *  ========
                    117: *
                    118: *> \author Univ. of Tennessee 
                    119: *> \author Univ. of California Berkeley 
                    120: *> \author Univ. of Colorado Denver 
                    121: *> \author NAG Ltd. 
                    122: *
1.12      bertrand  123: *> \date September 2012
1.9       bertrand  124: *
                    125: *> \ingroup doubleGEcomputational
                    126: *
                    127: *> \par Further Details:
                    128: *  =====================
                    129: *>
                    130: *> \verbatim
                    131: *>
                    132: *>  The matrix Q is represented as a product of elementary reflectors
                    133: *>
                    134: *>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
                    135: *>
                    136: *>  Each H(i) has the form
                    137: *>
                    138: *>     H(i) = I - tau * v * v**T
                    139: *>
1.12      bertrand  140: *>  where tau is a real scalar, and v is a real/complex vector
1.9       bertrand  141: *>  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
                    142: *>  A(i+1:m,i), and tau in TAU(i).
                    143: *> \endverbatim
                    144: *
                    145: *> \par Contributors:
                    146: *  ==================
                    147: *>
                    148: *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
                    149: *>    X. Sun, Computer Science Dept., Duke University, USA
                    150: *>
                    151: *  =====================================================================
1.1       bertrand  152:       SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
                    153: *
1.12      bertrand  154: *  -- LAPACK computational routine (version 3.4.2) --
1.1       bertrand  155: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    156: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12      bertrand  157: *     September 2012
1.1       bertrand  158: *
                    159: *     .. Scalar Arguments ..
                    160:       INTEGER            INFO, LDA, LWORK, M, N
                    161: *     ..
                    162: *     .. Array Arguments ..
                    163:       INTEGER            JPVT( * )
                    164:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
                    165: *     ..
                    166: *
                    167: *  =====================================================================
                    168: *
                    169: *     .. Parameters ..
                    170:       INTEGER            INB, INBMIN, IXOVER
                    171:       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
                    172: *     ..
                    173: *     .. Local Scalars ..
                    174:       LOGICAL            LQUERY
                    175:       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
                    176:      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
                    177: *     ..
                    178: *     .. External Subroutines ..
                    179:       EXTERNAL           DGEQRF, DLAQP2, DLAQPS, DORMQR, DSWAP, XERBLA
                    180: *     ..
                    181: *     .. External Functions ..
                    182:       INTEGER            ILAENV
                    183:       DOUBLE PRECISION   DNRM2
                    184:       EXTERNAL           ILAENV, DNRM2
                    185: *     ..
                    186: *     .. Intrinsic Functions ..
                    187:       INTRINSIC          INT, MAX, MIN
                    188: *     ..
                    189: *     .. Executable Statements ..
                    190: *
                    191: *     Test input arguments
1.9       bertrand  192: *  ====================
1.1       bertrand  193: *
                    194:       INFO = 0
                    195:       LQUERY = ( LWORK.EQ.-1 )
                    196:       IF( M.LT.0 ) THEN
                    197:          INFO = -1
                    198:       ELSE IF( N.LT.0 ) THEN
                    199:          INFO = -2
                    200:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
                    201:          INFO = -4
                    202:       END IF
                    203: *
                    204:       IF( INFO.EQ.0 ) THEN
                    205:          MINMN = MIN( M, N )
                    206:          IF( MINMN.EQ.0 ) THEN
                    207:             IWS = 1
                    208:             LWKOPT = 1
                    209:          ELSE
                    210:             IWS = 3*N + 1
                    211:             NB = ILAENV( INB, 'DGEQRF', ' ', M, N, -1, -1 )
                    212:             LWKOPT = 2*N + ( N + 1 )*NB
                    213:          END IF
                    214:          WORK( 1 ) = LWKOPT
                    215: *
                    216:          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
                    217:             INFO = -8
                    218:          END IF
                    219:       END IF
                    220: *
                    221:       IF( INFO.NE.0 ) THEN
                    222:          CALL XERBLA( 'DGEQP3', -INFO )
                    223:          RETURN
                    224:       ELSE IF( LQUERY ) THEN
                    225:          RETURN
                    226:       END IF
                    227: *
                    228: *     Quick return if possible.
                    229: *
                    230:       IF( MINMN.EQ.0 ) THEN
                    231:          RETURN
                    232:       END IF
                    233: *
                    234: *     Move initial columns up front.
                    235: *
                    236:       NFXD = 1
                    237:       DO 10 J = 1, N
                    238:          IF( JPVT( J ).NE.0 ) THEN
                    239:             IF( J.NE.NFXD ) THEN
                    240:                CALL DSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
                    241:                JPVT( J ) = JPVT( NFXD )
                    242:                JPVT( NFXD ) = J
                    243:             ELSE
                    244:                JPVT( J ) = J
                    245:             END IF
                    246:             NFXD = NFXD + 1
                    247:          ELSE
                    248:             JPVT( J ) = J
                    249:          END IF
                    250:    10 CONTINUE
                    251:       NFXD = NFXD - 1
                    252: *
                    253: *     Factorize fixed columns
1.9       bertrand  254: *  =======================
1.1       bertrand  255: *
                    256: *     Compute the QR factorization of fixed columns and update
                    257: *     remaining columns.
                    258: *
                    259:       IF( NFXD.GT.0 ) THEN
                    260:          NA = MIN( M, NFXD )
                    261: *CC      CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
                    262:          CALL DGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
                    263:          IWS = MAX( IWS, INT( WORK( 1 ) ) )
                    264:          IF( NA.LT.N ) THEN
                    265: *CC         CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
                    266: *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
                    267:             CALL DORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
                    268:      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
                    269:             IWS = MAX( IWS, INT( WORK( 1 ) ) )
                    270:          END IF
                    271:       END IF
                    272: *
                    273: *     Factorize free columns
1.9       bertrand  274: *  ======================
1.1       bertrand  275: *
                    276:       IF( NFXD.LT.MINMN ) THEN
                    277: *
                    278:          SM = M - NFXD
                    279:          SN = N - NFXD
                    280:          SMINMN = MINMN - NFXD
                    281: *
                    282: *        Determine the block size.
                    283: *
                    284:          NB = ILAENV( INB, 'DGEQRF', ' ', SM, SN, -1, -1 )
                    285:          NBMIN = 2
                    286:          NX = 0
                    287: *
                    288:          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
                    289: *
                    290: *           Determine when to cross over from blocked to unblocked code.
                    291: *
                    292:             NX = MAX( 0, ILAENV( IXOVER, 'DGEQRF', ' ', SM, SN, -1,
                    293:      $           -1 ) )
                    294: *
                    295: *
                    296:             IF( NX.LT.SMINMN ) THEN
                    297: *
                    298: *              Determine if workspace is large enough for blocked code.
                    299: *
                    300:                MINWS = 2*SN + ( SN+1 )*NB
                    301:                IWS = MAX( IWS, MINWS )
                    302:                IF( LWORK.LT.MINWS ) THEN
                    303: *
                    304: *                 Not enough workspace to use optimal NB: Reduce NB and
                    305: *                 determine the minimum value of NB.
                    306: *
                    307:                   NB = ( LWORK-2*SN ) / ( SN+1 )
                    308:                   NBMIN = MAX( 2, ILAENV( INBMIN, 'DGEQRF', ' ', SM, SN,
                    309:      $                    -1, -1 ) )
                    310: *
                    311: *
                    312:                END IF
                    313:             END IF
                    314:          END IF
                    315: *
                    316: *        Initialize partial column norms. The first N elements of work
                    317: *        store the exact column norms.
                    318: *
                    319:          DO 20 J = NFXD + 1, N
                    320:             WORK( J ) = DNRM2( SM, A( NFXD+1, J ), 1 )
                    321:             WORK( N+J ) = WORK( J )
                    322:    20    CONTINUE
                    323: *
                    324:          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
                    325:      $       ( NX.LT.SMINMN ) ) THEN
                    326: *
                    327: *           Use blocked code initially.
                    328: *
                    329:             J = NFXD + 1
                    330: *
                    331: *           Compute factorization: while loop.
                    332: *
                    333: *
                    334:             TOPBMN = MINMN - NX
                    335:    30       CONTINUE
                    336:             IF( J.LE.TOPBMN ) THEN
                    337:                JB = MIN( NB, TOPBMN-J+1 )
                    338: *
                    339: *              Factorize JB columns among columns J:N.
                    340: *
                    341:                CALL DLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
                    342:      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
                    343:      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
                    344: *
                    345:                J = J + FJB
                    346:                GO TO 30
                    347:             END IF
                    348:          ELSE
                    349:             J = NFXD + 1
                    350:          END IF
                    351: *
                    352: *        Use unblocked code to factor the last or only block.
                    353: *
                    354: *
                    355:          IF( J.LE.MINMN )
                    356:      $      CALL DLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
                    357:      $                   TAU( J ), WORK( J ), WORK( N+J ),
                    358:      $                   WORK( 2*N+1 ) )
                    359: *
                    360:       END IF
                    361: *
                    362:       WORK( 1 ) = IWS
                    363:       RETURN
                    364: *
                    365: *     End of DGEQP3
                    366: *
                    367:       END

CVSweb interface <joel.bertrand@systella.fr>