LAPACK: SRC/dgeqp3.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, LDA, LWORK, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            JPVT( * )
00013       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DGEQP3 computes a QR factorization with column pivoting of a
00020 *  matrix A:  A*P = Q*R  using Level 3 BLAS.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  M       (input) INTEGER
00026 *          The number of rows of the matrix A. M >= 0.
00027 *
00028 *  N       (input) INTEGER
00029 *          The number of columns of the matrix A.  N >= 0.
00030 *
00031 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00032 *          On entry, the M-by-N matrix A.
00033 *          On exit, the upper triangle of the array contains the
00034 *          min(M,N)-by-N upper trapezoidal matrix R; the elements below
00035 *          the diagonal, together with the array TAU, represent the
00036 *          orthogonal matrix Q as a product of min(M,N) elementary
00037 *          reflectors.
00038 *
00039 *  LDA     (input) INTEGER
00040 *          The leading dimension of the array A. LDA >= max(1,M).
00041 *
00042 *  JPVT    (input/output) INTEGER array, dimension (N)
00043 *          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
00044 *          to the front of A*P (a leading column); if JPVT(J)=0,
00045 *          the J-th column of A is a free column.
00046 *          On exit, if JPVT(J)=K, then the J-th column of A*P was the
00047 *          the K-th column of A.
00048 *
00049 *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
00050 *          The scalar factors of the elementary reflectors.
00051 *
00052 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00053 *          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
00054 *
00055 *  LWORK   (input) INTEGER
00056 *          The dimension of the array WORK. LWORK >= 3*N+1.
00057 *          For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
00058 *          is the optimal blocksize.
00059 *
00060 *          If LWORK = -1, then a workspace query is assumed; the routine
00061 *          only calculates the optimal size of the WORK array, returns
00062 *          this value as the first entry of the WORK array, and no error
00063 *          message related to LWORK is issued by XERBLA.
00064 *
00065 *  INFO    (output) INTEGER
00066 *          = 0: successful exit.
00067 *          < 0: if INFO = -i, the i-th argument had an illegal value.
00068 *
00069 *  Further Details
00070 *  ===============
00071 *
00072 *  The matrix Q is represented as a product of elementary reflectors
00073 *
00074 *     Q = H(1) H(2) . . . H(k), where k = min(m,n).
00075 *
00076 *  Each H(i) has the form
00077 *
00078 *     H(i) = I - tau * v * v'
00079 *
00080 *  where tau is a real/complex scalar, and v is a real/complex vector
00081 *  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
00082 *  A(i+1:m,i), and tau in TAU(i).
00083 *
00084 *  Based on contributions by
00085 *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00086 *    X. Sun, Computer Science Dept., Duke University, USA
00087 *
00088 *  =====================================================================
00089 *
00090 *     .. Parameters ..
00091       INTEGER            INB, INBMIN, IXOVER
00092       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
00093 *     ..
00094 *     .. Local Scalars ..
00095       LOGICAL            LQUERY
00096       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
00097      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
00098 *     ..
00099 *     .. External Subroutines ..
00100       EXTERNAL           DGEQRF, DLAQP2, DLAQPS, DORMQR, DSWAP, XERBLA
00101 *     ..
00102 *     .. External Functions ..
00103       INTEGER            ILAENV
00104       DOUBLE PRECISION   DNRM2
00105       EXTERNAL           ILAENV, DNRM2
00106 *     ..
00107 *     .. Intrinsic Functions ..
00108       INTRINSIC          INT, MAX, MIN
00109 *     ..
00110 *     .. Executable Statements ..
00111 *
00112 *     Test input arguments
00113 *     ====================
00114 *
00115       INFO = 0
00116       LQUERY = ( LWORK.EQ.-1 )
00117       IF( M.LT.0 ) THEN
00118          INFO = -1
00119       ELSE IF( N.LT.0 ) THEN
00120          INFO = -2
00121       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00122          INFO = -4
00123       END IF
00124 *
00125       IF( INFO.EQ.0 ) THEN
00126          MINMN = MIN( M, N )
00127          IF( MINMN.EQ.0 ) THEN
00128             IWS = 1
00129             LWKOPT = 1
00130          ELSE
00131             IWS = 3*N + 1
00132             NB = ILAENV( INB, 'DGEQRF', ' ', M, N, -1, -1 )
00133             LWKOPT = 2*N + ( N + 1 )*NB
00134          END IF
00135          WORK( 1 ) = LWKOPT
00136 *
00137          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
00138             INFO = -8
00139          END IF
00140       END IF
00141 *
00142       IF( INFO.NE.0 ) THEN
00143          CALL XERBLA( 'DGEQP3', -INFO )
00144          RETURN
00145       ELSE IF( LQUERY ) THEN
00146          RETURN
00147       END IF
00148 *
00149 *     Quick return if possible.
00150 *
00151       IF( MINMN.EQ.0 ) THEN
00152          RETURN
00153       END IF
00154 *
00155 *     Move initial columns up front.
00156 *
00157       NFXD = 1
00158       DO 10 J = 1, N
00159          IF( JPVT( J ).NE.0 ) THEN
00160             IF( J.NE.NFXD ) THEN
00161                CALL DSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
00162                JPVT( J ) = JPVT( NFXD )
00163                JPVT( NFXD ) = J
00164             ELSE
00165                JPVT( J ) = J
00166             END IF
00167             NFXD = NFXD + 1
00168          ELSE
00169             JPVT( J ) = J
00170          END IF
00171    10 CONTINUE
00172       NFXD = NFXD - 1
00173 *
00174 *     Factorize fixed columns
00175 *     =======================
00176 *
00177 *     Compute the QR factorization of fixed columns and update
00178 *     remaining columns.
00179 *
00180       IF( NFXD.GT.0 ) THEN
00181          NA = MIN( M, NFXD )
00182 *CC      CALL DGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
00183          CALL DGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
00184          IWS = MAX( IWS, INT( WORK( 1 ) ) )
00185          IF( NA.LT.N ) THEN
00186 *CC         CALL DORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
00187 *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
00188             CALL DORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
00189      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
00190             IWS = MAX( IWS, INT( WORK( 1 ) ) )
00191          END IF
00192       END IF
00193 *
00194 *     Factorize free columns
00195 *     ======================
00196 *
00197       IF( NFXD.LT.MINMN ) THEN
00198 *
00199          SM = M - NFXD
00200          SN = N - NFXD
00201          SMINMN = MINMN - NFXD
00202 *
00203 *        Determine the block size.
00204 *
00205          NB = ILAENV( INB, 'DGEQRF', ' ', SM, SN, -1, -1 )
00206          NBMIN = 2
00207          NX = 0
00208 *
00209          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
00210 *
00211 *           Determine when to cross over from blocked to unblocked code.
00212 *
00213             NX = MAX( 0, ILAENV( IXOVER, 'DGEQRF', ' ', SM, SN, -1,
00214      $           -1 ) )
00215 *
00216 *
00217             IF( NX.LT.SMINMN ) THEN
00218 *
00219 *              Determine if workspace is large enough for blocked code.
00220 *
00221                MINWS = 2*SN + ( SN+1 )*NB
00222                IWS = MAX( IWS, MINWS )
00223                IF( LWORK.LT.MINWS ) THEN
00224 *
00225 *                 Not enough workspace to use optimal NB: Reduce NB and
00226 *                 determine the minimum value of NB.
00227 *
00228                   NB = ( LWORK-2*SN ) / ( SN+1 )
00229                   NBMIN = MAX( 2, ILAENV( INBMIN, 'DGEQRF', ' ', SM, SN,
00230      $                    -1, -1 ) )
00231 *
00232 *
00233                END IF
00234             END IF
00235          END IF
00236 *
00237 *        Initialize partial column norms. The first N elements of work
00238 *        store the exact column norms.
00239 *
00240          DO 20 J = NFXD + 1, N
00241             WORK( J ) = DNRM2( SM, A( NFXD+1, J ), 1 )
00242             WORK( N+J ) = WORK( J )
00243    20    CONTINUE
00244 *
00245          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
00246      $       ( NX.LT.SMINMN ) ) THEN
00247 *
00248 *           Use blocked code initially.
00249 *
00250             J = NFXD + 1
00251 *
00252 *           Compute factorization: while loop.
00253 *
00254 *
00255             TOPBMN = MINMN - NX
00256    30       CONTINUE
00257             IF( J.LE.TOPBMN ) THEN
00258                JB = MIN( NB, TOPBMN-J+1 )
00259 *
00260 *              Factorize JB columns among columns J:N.
00261 *
00262                CALL DLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
00263      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
00264      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
00265 *
00266                J = J + FJB
00267                GO TO 30
00268             END IF
00269          ELSE
00270             J = NFXD + 1
00271          END IF
00272 *
00273 *        Use unblocked code to factor the last or only block.
00274 *
00275 *
00276          IF( J.LE.MINMN )
00277      $      CALL DLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
00278      $                   TAU( J ), WORK( J ), WORK( N+J ),
00279      $                   WORK( 2*N+1 ) )
00280 *
00281       END IF
00282 *
00283       WORK( 1 ) = IWS
00284       RETURN
00285 *
00286 *     End of DGEQP3
00287 *
00288       END