LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgeqp3.f
Go to the documentation of this file.
00001 *> \brief \b SGEQP3
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGEQP3 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqp3.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqp3.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqp3.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, LWORK, M, N
00025 *       ..
00026 *       .. Array Arguments ..
00027 *       INTEGER            JPVT( * )
00028 *       REAL               A( LDA, * ), TAU( * ), WORK( * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> SGEQP3 computes a QR factorization with column pivoting of a
00038 *> matrix A:  A*P = Q*R  using Level 3 BLAS.
00039 *> \endverbatim
00040 *
00041 *  Arguments:
00042 *  ==========
00043 *
00044 *> \param[in] M
00045 *> \verbatim
00046 *>          M is INTEGER
00047 *>          The number of rows of the matrix A. M >= 0.
00048 *> \endverbatim
00049 *>
00050 *> \param[in] N
00051 *> \verbatim
00052 *>          N is INTEGER
00053 *>          The number of columns of the matrix A.  N >= 0.
00054 *> \endverbatim
00055 *>
00056 *> \param[in,out] A
00057 *> \verbatim
00058 *>          A is REAL array, dimension (LDA,N)
00059 *>          On entry, the M-by-N matrix A.
00060 *>          On exit, the upper triangle of the array contains the
00061 *>          min(M,N)-by-N upper trapezoidal matrix R; the elements below
00062 *>          the diagonal, together with the array TAU, represent the
00063 *>          orthogonal matrix Q as a product of min(M,N) elementary
00064 *>          reflectors.
00065 *> \endverbatim
00066 *>
00067 *> \param[in] LDA
00068 *> \verbatim
00069 *>          LDA is INTEGER
00070 *>          The leading dimension of the array A. LDA >= max(1,M).
00071 *> \endverbatim
00072 *>
00073 *> \param[in,out] JPVT
00074 *> \verbatim
00075 *>          JPVT is INTEGER array, dimension (N)
00076 *>          On entry, if JPVT(J).ne.0, the J-th column of A is permuted
00077 *>          to the front of A*P (a leading column); if JPVT(J)=0,
00078 *>          the J-th column of A is a free column.
00079 *>          On exit, if JPVT(J)=K, then the J-th column of A*P was the
00080 *>          the K-th column of A.
00081 *> \endverbatim
00082 *>
00083 *> \param[out] TAU
00084 *> \verbatim
00085 *>          TAU is REAL array, dimension (min(M,N))
00086 *>          The scalar factors of the elementary reflectors.
00087 *> \endverbatim
00088 *>
00089 *> \param[out] WORK
00090 *> \verbatim
00091 *>          WORK is REAL array, dimension (MAX(1,LWORK))
00092 *>          On exit, if INFO=0, WORK(1) returns the optimal LWORK.
00093 *> \endverbatim
00094 *>
00095 *> \param[in] LWORK
00096 *> \verbatim
00097 *>          LWORK is INTEGER
00098 *>          The dimension of the array WORK. LWORK >= 3*N+1.
00099 *>          For optimal performance LWORK >= 2*N+( N+1 )*NB, where NB
00100 *>          is the optimal blocksize.
00101 *>
00102 *>          If LWORK = -1, then a workspace query is assumed; the routine
00103 *>          only calculates the optimal size of the WORK array, returns
00104 *>          this value as the first entry of the WORK array, and no error
00105 *>          message related to LWORK is issued by XERBLA.
00106 *> \endverbatim
00107 *>
00108 *> \param[out] INFO
00109 *> \verbatim
00110 *>          INFO is INTEGER
00111 *>          = 0: successful exit.
00112 *>          < 0: if INFO = -i, the i-th argument had an illegal value.
00113 *> \endverbatim
00114 *
00115 *  Authors:
00116 *  ========
00117 *
00118 *> \author Univ. of Tennessee 
00119 *> \author Univ. of California Berkeley 
00120 *> \author Univ. of Colorado Denver 
00121 *> \author NAG Ltd. 
00122 *
00123 *> \date November 2011
00124 *
00125 *> \ingroup realGEcomputational
00126 *
00127 *> \par Further Details:
00128 *  =====================
00129 *>
00130 *> \verbatim
00131 *>
00132 *>  The matrix Q is represented as a product of elementary reflectors
00133 *>
00134 *>     Q = H(1) H(2) . . . H(k), where k = min(m,n).
00135 *>
00136 *>  Each H(i) has the form
00137 *>
00138 *>     H(i) = I - tau * v * v**T
00139 *>
00140 *>  where tau is a real/complex scalar, and v is a real/complex vector
00141 *>  with v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in
00142 *>  A(i+1:m,i), and tau in TAU(i).
00143 *> \endverbatim
00144 *
00145 *> \par Contributors:
00146 *  ==================
00147 *>
00148 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00149 *>    X. Sun, Computer Science Dept., Duke University, USA
00150 *>
00151 *  =====================================================================
00152       SUBROUTINE SGEQP3( M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO )
00153 *
00154 *  -- LAPACK computational routine (version 3.4.0) --
00155 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00156 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00157 *     November 2011
00158 *
00159 *     .. Scalar Arguments ..
00160       INTEGER            INFO, LDA, LWORK, M, N
00161 *     ..
00162 *     .. Array Arguments ..
00163       INTEGER            JPVT( * )
00164       REAL               A( LDA, * ), TAU( * ), WORK( * )
00165 *     ..
00166 *
00167 *  =====================================================================
00168 *
00169 *     .. Parameters ..
00170       INTEGER            INB, INBMIN, IXOVER
00171       PARAMETER          ( INB = 1, INBMIN = 2, IXOVER = 3 )
00172 *     ..
00173 *     .. Local Scalars ..
00174       LOGICAL            LQUERY
00175       INTEGER            FJB, IWS, J, JB, LWKOPT, MINMN, MINWS, NA, NB,
00176      $                   NBMIN, NFXD, NX, SM, SMINMN, SN, TOPBMN
00177 *     ..
00178 *     .. External Subroutines ..
00179       EXTERNAL           SGEQRF, SLAQP2, SLAQPS, SORMQR, SSWAP, XERBLA
00180 *     ..
00181 *     .. External Functions ..
00182       INTEGER            ILAENV
00183       REAL               SNRM2
00184       EXTERNAL           ILAENV, SNRM2
00185 *     ..
00186 *     .. Intrinsic Functions ..
00187       INTRINSIC          INT, MAX, MIN
00188 *     ..
00189 *     .. Executable Statements ..
00190 *
00191       INFO = 0
00192       LQUERY = ( LWORK.EQ.-1 )
00193       IF( M.LT.0 ) THEN
00194          INFO = -1
00195       ELSE IF( N.LT.0 ) THEN
00196          INFO = -2
00197       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00198          INFO = -4
00199       END IF
00200 *
00201       IF( INFO.EQ.0 ) THEN
00202          MINMN = MIN( M, N )
00203          IF( MINMN.EQ.0 ) THEN
00204             IWS = 1
00205             LWKOPT = 1
00206          ELSE
00207             IWS = 3*N + 1
00208             NB = ILAENV( INB, 'SGEQRF', ' ', M, N, -1, -1 )
00209             LWKOPT = 2*N + ( N + 1 )*NB
00210          END IF
00211          WORK( 1 ) = LWKOPT
00212 *
00213          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
00214             INFO = -8
00215          END IF
00216       END IF
00217 *
00218       IF( INFO.NE.0 ) THEN
00219          CALL XERBLA( 'SGEQP3', -INFO )
00220          RETURN
00221       ELSE IF( LQUERY ) THEN
00222          RETURN
00223       END IF
00224 *
00225 *     Quick return if possible.
00226 *
00227       IF( MINMN.EQ.0 ) THEN
00228          RETURN
00229       END IF
00230 *
00231 *     Move initial columns up front.
00232 *
00233       NFXD = 1
00234       DO 10 J = 1, N
00235          IF( JPVT( J ).NE.0 ) THEN
00236             IF( J.NE.NFXD ) THEN
00237                CALL SSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
00238                JPVT( J ) = JPVT( NFXD )
00239                JPVT( NFXD ) = J
00240             ELSE
00241                JPVT( J ) = J
00242             END IF
00243             NFXD = NFXD + 1
00244          ELSE
00245             JPVT( J ) = J
00246          END IF
00247    10 CONTINUE
00248       NFXD = NFXD - 1
00249 *
00250 *     Factorize fixed columns
00251 *  =======================
00252 *
00253 *     Compute the QR factorization of fixed columns and update
00254 *     remaining columns.
00255 *
00256       IF( NFXD.GT.0 ) THEN
00257          NA = MIN( M, NFXD )
00258 *CC      CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
00259          CALL SGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
00260          IWS = MAX( IWS, INT( WORK( 1 ) ) )
00261          IF( NA.LT.N ) THEN
00262 *CC         CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
00263 *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
00264             CALL SORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
00265      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
00266             IWS = MAX( IWS, INT( WORK( 1 ) ) )
00267          END IF
00268       END IF
00269 *
00270 *     Factorize free columns
00271 *  ======================
00272 *
00273       IF( NFXD.LT.MINMN ) THEN
00274 *
00275          SM = M - NFXD
00276          SN = N - NFXD
00277          SMINMN = MINMN - NFXD
00278 *
00279 *        Determine the block size.
00280 *
00281          NB = ILAENV( INB, 'SGEQRF', ' ', SM, SN, -1, -1 )
00282          NBMIN = 2
00283          NX = 0
00284 *
00285          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
00286 *
00287 *           Determine when to cross over from blocked to unblocked code.
00288 *
00289             NX = MAX( 0, ILAENV( IXOVER, 'SGEQRF', ' ', SM, SN, -1,
00290      $           -1 ) )
00291 *
00292 *
00293             IF( NX.LT.SMINMN ) THEN
00294 *
00295 *              Determine if workspace is large enough for blocked code.
00296 *
00297                MINWS = 2*SN + ( SN+1 )*NB
00298                IWS = MAX( IWS, MINWS )
00299                IF( LWORK.LT.MINWS ) THEN
00300 *
00301 *                 Not enough workspace to use optimal NB: Reduce NB and
00302 *                 determine the minimum value of NB.
00303 *
00304                   NB = ( LWORK-2*SN ) / ( SN+1 )
00305                   NBMIN = MAX( 2, ILAENV( INBMIN, 'SGEQRF', ' ', SM, SN,
00306      $                    -1, -1 ) )
00307 *
00308 *
00309                END IF
00310             END IF
00311          END IF
00312 *
00313 *        Initialize partial column norms. The first N elements of work
00314 *        store the exact column norms.
00315 *
00316          DO 20 J = NFXD + 1, N
00317             WORK( J ) = SNRM2( SM, A( NFXD+1, J ), 1 )
00318             WORK( N+J ) = WORK( J )
00319    20    CONTINUE
00320 *
00321          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
00322      $       ( NX.LT.SMINMN ) ) THEN
00323 *
00324 *           Use blocked code initially.
00325 *
00326             J = NFXD + 1
00327 *
00328 *           Compute factorization: while loop.
00329 *
00330 *
00331             TOPBMN = MINMN - NX
00332    30       CONTINUE
00333             IF( J.LE.TOPBMN ) THEN
00334                JB = MIN( NB, TOPBMN-J+1 )
00335 *
00336 *              Factorize JB columns among columns J:N.
00337 *
00338                CALL SLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
00339      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
00340      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
00341 *
00342                J = J + FJB
00343                GO TO 30
00344             END IF
00345          ELSE
00346             J = NFXD + 1
00347          END IF
00348 *
00349 *        Use unblocked code to factor the last or only block.
00350 *
00351 *
00352          IF( J.LE.MINMN )
00353      $      CALL SLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
00354      $                   TAU( J ), WORK( J ), WORK( N+J ),
00355      $                   WORK( 2*N+1 ) )
00356 *
00357       END IF
00358 *
00359       WORK( 1 ) = IWS
00360       RETURN
00361 *
00362 *     End of SGEQP3
00363 *
00364       END
 All Files Functions