SRC/sgeqp3.f

   1 *> \brief \b SGEQP3
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download SGEQP3 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqp3.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqp3.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqp3.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SGEQP3( 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 *       REAL               A( LDA, * ), TAU( * ), WORK( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> SGEQP3 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 REAL 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 REAL 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 REAL 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 *
 123 *> \date November 2015
 124 *
 125 *> \ingroup realGEcomputational
 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 *>
 140 *>  where tau is a real scalar, and v is a real/complex vector
 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 *  =====================================================================
 152       SUBROUTINE SGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
 153 *
 154 *  -- LAPACK computational routine (version 3.6.0) --
 155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 157 *     November 2015
 158 *
 159 *     .. Scalar Arguments ..
 160       INTEGER            INFO, LDA, LWORK, M, N
 161 *     ..
 162 *     .. Array Arguments ..
 163       INTEGER            JPVT( * )
 164       REAL               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           SGEQRF, SLAQP2, SLAQPS, SORMQR, SSWAP, XERBLA
 180 *     ..
 181 *     .. External Functions ..
 182       INTEGER            ILAENV
 183       REAL               SNRM2
 184       EXTERNAL           ILAENV, SNRM2
 185 *     ..
 186 *     .. Intrinsic Functions ..
 187       INTRINSIC          INT, MAX, MIN
 188 *     Test input arguments
 189 *  ====================
 190 *
 191       INFO = 0
 192       LQUERY = ( LWORK.EQ.-1 )
 193       IF( M.LT.0 ) THEN
 194          INFO = -1
 195       ELSE IF( N.LT.0 ) THEN
 196          INFO = -2
 197       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 198          INFO = -4
 199       END IF
 200 *
 201       IF( INFO.EQ.0 ) THEN
 202          MINMN = MIN( M, N )
 203          IF( MINMN.EQ.0 ) THEN
 204             IWS = 1
 205             LWKOPT = 1
 206          ELSE
 207             IWS = 3*N + 1
 208             NB = ILAENV( INB, 'SGEQRF', ' ', M, N, -1, -1 )
 209             LWKOPT = 2*N + ( N + 1 )*NB
 210          END IF
 211          WORK( 1 ) = LWKOPT
 212 *
 213          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
 214             INFO = -8
 215          END IF
 216       END IF
 217 *
 218       IF( INFO.NE.0 ) THEN
 219          CALL XERBLA( 'SGEQP3', -INFO )
 220          RETURN
 221       ELSE IF( LQUERY ) THEN
 222          RETURN
 223       END IF
 224 *
 225 *     Move initial columns up front.
 226 *
 227       NFXD = 1
 228       DO 10 J = 1, N
 229          IF( JPVT( J ).NE.0 ) THEN
 230             IF( J.NE.NFXD ) THEN
 231                CALL SSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
 232                JPVT( J ) = JPVT( NFXD )
 233                JPVT( NFXD ) = J
 234             ELSE
 235                JPVT( J ) = J
 236             END IF
 237             NFXD = NFXD + 1
 238          ELSE
 239             JPVT( J ) = J
 240          END IF
 241    10 CONTINUE
 242       NFXD = NFXD - 1
 243 *
 244 *     Factorize fixed columns
 245 *  =======================
 246 *
 247 *     Compute the QR factorization of fixed columns and update
 248 *     remaining columns.
 249 *
 250       IF( NFXD.GT.0 ) THEN
 251          NA = MIN( M, NFXD )
 252 *CC      CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
 253          CALL SGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
 254          IWS = MAX( IWS, INT( WORK( 1 ) ) )
 255          IF( NA.LT.N ) THEN
 256 *CC         CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
 257 *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
 258             CALL SORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
 259      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
 260             IWS = MAX( IWS, INT( WORK( 1 ) ) )
 261          END IF
 262       END IF
 263 *
 264 *     Factorize free columns
 265 *  ======================
 266 *
 267       IF( NFXD.LT.MINMN ) THEN
 268 *
 269          SM = M - NFXD
 270          SN = N - NFXD
 271          SMINMN = MINMN - NFXD
 272 *
 273 *        Determine the block size.
 274 *
 275          NB = ILAENV( INB, 'SGEQRF', ' ', SM, SN, -1, -1 )
 276          NBMIN = 2
 277          NX = 0
 278 *
 279          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
 280 *
 281 *           Determine when to cross over from blocked to unblocked code.
 282 *
 283             NX = MAX( 0, ILAENV( IXOVER, 'SGEQRF', ' ', SM, SN, -1,
 284      $           -1 ) )
 285 *
 286 *
 287             IF( NX.LT.SMINMN ) THEN
 288 *
 289 *              Determine if workspace is large enough for blocked code.
 290 *
 291                MINWS = 2*SN + ( SN+1 )*NB
 292                IWS = MAX( IWS, MINWS )
 293                IF( LWORK.LT.MINWS ) THEN
 294 *
 295 *                 Not enough workspace to use optimal NB: Reduce NB and
 296 *                 determine the minimum value of NB.
 297 *
 298                   NB = ( LWORK-2*SN ) / ( SN+1 )
 299                   NBMIN = MAX( 2, ILAENV( INBMIN, 'SGEQRF', ' ', SM, SN,
 300      $                    -1, -1 ) )
 301 *
 302 *
 303                END IF
 304             END IF
 305          END IF
 306 *
 307 *        Initialize partial column norms. The first N elements of work
 308 *        store the exact column norms.
 309 *
 310          DO 20 J = NFXD + 1, N
 311             WORK( J ) = SNRM2( SM, A( NFXD+1, J ), 1 )
 312             WORK( N+J ) = WORK( J )
 313    20    CONTINUE
 314 *
 315          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
 316      $       ( NX.LT.SMINMN ) ) THEN
 317 *
 318 *           Use blocked code initially.
 319 *
 320             J = NFXD + 1
 321 *
 322 *           Compute factorization: while loop.
 323 *
 324 *
 325             TOPBMN = MINMN - NX
 326    30       CONTINUE
 327             IF( J.LE.TOPBMN ) THEN
 328                JB = MIN( NB, TOPBMN-J+1 )
 329 *
 330 *              Factorize JB columns among columns J:N.
 331 *
 332                CALL SLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
 333      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
 334      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
 335 *
 336                J = J + FJB
 337                GO TO 30
 338             END IF
 339          ELSE
 340             J = NFXD + 1
 341          END IF
 342 *
 343 *        Use unblocked code to factor the last or only block.
 344 *
 345 *
 346          IF( J.LE.MINMN )
 347      $      CALL SLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
 348      $                   TAU( J ), WORK( J ), WORK( N+J ),
 349      $                   WORK( 2*N+1 ) )
 350 *
 351       END IF
 352 *
 353       WORK( 1 ) = IWS
 354       RETURN
 355 *
 356 *     End of SGEQP3
 357 *
 358       END