SRC/cgeqp3.f

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