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

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

CVSweb interface <joel.bertrand@systella.fr>