LAPACK: SRC/sgeqp3.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGEQP3( 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       REAL               A( LDA, * ), TAU( * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SGEQP3 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) REAL 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) REAL array, dimension (min(M,N))
00050 *          The scalar factors of the elementary reflectors.
00051 *
00052 *  WORK    (workspace/output) REAL 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           SGEQRF, SLAQP2, SLAQPS, SORMQR, SSWAP, XERBLA
00101 *     ..
00102 *     .. External Functions ..
00103       INTEGER            ILAENV
00104       REAL               SNRM2
00105       EXTERNAL           ILAENV, SNRM2
00106 *     ..
00107 *     .. Intrinsic Functions ..
00108       INTRINSIC          INT, MAX, MIN
00109 *     ..
00110 *     .. Executable Statements ..
00111 *
00112       INFO = 0
00113       LQUERY = ( LWORK.EQ.-1 )
00114       IF( M.LT.0 ) THEN
00115          INFO = -1
00116       ELSE IF( N.LT.0 ) THEN
00117          INFO = -2
00118       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00119          INFO = -4
00120       END IF
00121 *
00122       IF( INFO.EQ.0 ) THEN
00123          MINMN = MIN( M, N )
00124          IF( MINMN.EQ.0 ) THEN
00125             IWS = 1
00126             LWKOPT = 1
00127          ELSE
00128             IWS = 3*N + 1
00129             NB = ILAENV( INB, 'SGEQRF', ' ', M, N, -1, -1 )
00130             LWKOPT = 2*N + ( N + 1 )*NB
00131          END IF
00132          WORK( 1 ) = LWKOPT
00133 *
00134          IF( ( LWORK.LT.IWS ) .AND. .NOT.LQUERY ) THEN
00135             INFO = -8
00136          END IF
00137       END IF
00138 *
00139       IF( INFO.NE.0 ) THEN
00140          CALL XERBLA( 'SGEQP3', -INFO )
00141          RETURN
00142       ELSE IF( LQUERY ) THEN
00143          RETURN
00144       END IF
00145 *
00146 *     Quick return if possible.
00147 *
00148       IF( MINMN.EQ.0 ) THEN
00149          RETURN
00150       END IF
00151 *
00152 *     Move initial columns up front.
00153 *
00154       NFXD = 1
00155       DO 10 J = 1, N
00156          IF( JPVT( J ).NE.0 ) THEN
00157             IF( J.NE.NFXD ) THEN
00158                CALL SSWAP( M, A( 1, J ), 1, A( 1, NFXD ), 1 )
00159                JPVT( J ) = JPVT( NFXD )
00160                JPVT( NFXD ) = J
00161             ELSE
00162                JPVT( J ) = J
00163             END IF
00164             NFXD = NFXD + 1
00165          ELSE
00166             JPVT( J ) = J
00167          END IF
00168    10 CONTINUE
00169       NFXD = NFXD - 1
00170 *
00171 *     Factorize fixed columns
00172 *     =======================
00173 *
00174 *     Compute the QR factorization of fixed columns and update
00175 *     remaining columns.
00176 *
00177       IF( NFXD.GT.0 ) THEN
00178          NA = MIN( M, NFXD )
00179 *CC      CALL SGEQR2( M, NA, A, LDA, TAU, WORK, INFO )
00180          CALL SGEQRF( M, NA, A, LDA, TAU, WORK, LWORK, INFO )
00181          IWS = MAX( IWS, INT( WORK( 1 ) ) )
00182          IF( NA.LT.N ) THEN
00183 *CC         CALL SORM2R( 'Left', 'Transpose', M, N-NA, NA, A, LDA,
00184 *CC  $                   TAU, A( 1, NA+1 ), LDA, WORK, INFO )
00185             CALL SORMQR( 'Left', 'Transpose', M, N-NA, NA, A, LDA, TAU,
00186      $                   A( 1, NA+1 ), LDA, WORK, LWORK, INFO )
00187             IWS = MAX( IWS, INT( WORK( 1 ) ) )
00188          END IF
00189       END IF
00190 *
00191 *     Factorize free columns
00192 *     ======================
00193 *
00194       IF( NFXD.LT.MINMN ) THEN
00195 *
00196          SM = M - NFXD
00197          SN = N - NFXD
00198          SMINMN = MINMN - NFXD
00199 *
00200 *        Determine the block size.
00201 *
00202          NB = ILAENV( INB, 'SGEQRF', ' ', SM, SN, -1, -1 )
00203          NBMIN = 2
00204          NX = 0
00205 *
00206          IF( ( NB.GT.1 ) .AND. ( NB.LT.SMINMN ) ) THEN
00207 *
00208 *           Determine when to cross over from blocked to unblocked code.
00209 *
00210             NX = MAX( 0, ILAENV( IXOVER, 'SGEQRF', ' ', SM, SN, -1,
00211      $           -1 ) )
00212 *
00213 *
00214             IF( NX.LT.SMINMN ) THEN
00215 *
00216 *              Determine if workspace is large enough for blocked code.
00217 *
00218                MINWS = 2*SN + ( SN+1 )*NB
00219                IWS = MAX( IWS, MINWS )
00220                IF( LWORK.LT.MINWS ) THEN
00221 *
00222 *                 Not enough workspace to use optimal NB: Reduce NB and
00223 *                 determine the minimum value of NB.
00224 *
00225                   NB = ( LWORK-2*SN ) / ( SN+1 )
00226                   NBMIN = MAX( 2, ILAENV( INBMIN, 'SGEQRF', ' ', SM, SN,
00227      $                    -1, -1 ) )
00228 *
00229 *
00230                END IF
00231             END IF
00232          END IF
00233 *
00234 *        Initialize partial column norms. The first N elements of work
00235 *        store the exact column norms.
00236 *
00237          DO 20 J = NFXD + 1, N
00238             WORK( J ) = SNRM2( SM, A( NFXD+1, J ), 1 )
00239             WORK( N+J ) = WORK( J )
00240    20    CONTINUE
00241 *
00242          IF( ( NB.GE.NBMIN ) .AND. ( NB.LT.SMINMN ) .AND.
00243      $       ( NX.LT.SMINMN ) ) THEN
00244 *
00245 *           Use blocked code initially.
00246 *
00247             J = NFXD + 1
00248 *
00249 *           Compute factorization: while loop.
00250 *
00251 *
00252             TOPBMN = MINMN - NX
00253    30       CONTINUE
00254             IF( J.LE.TOPBMN ) THEN
00255                JB = MIN( NB, TOPBMN-J+1 )
00256 *
00257 *              Factorize JB columns among columns J:N.
00258 *
00259                CALL SLAQPS( M, N-J+1, J-1, JB, FJB, A( 1, J ), LDA,
00260      $                      JPVT( J ), TAU( J ), WORK( J ), WORK( N+J ),
00261      $                      WORK( 2*N+1 ), WORK( 2*N+JB+1 ), N-J+1 )
00262 *
00263                J = J + FJB
00264                GO TO 30
00265             END IF
00266          ELSE
00267             J = NFXD + 1
00268          END IF
00269 *
00270 *        Use unblocked code to factor the last or only block.
00271 *
00272 *
00273          IF( J.LE.MINMN )
00274      $      CALL SLAQP2( M, N-J+1, J-1, A( 1, J ), LDA, JPVT( J ),
00275      $                   TAU( J ), WORK( J ), WORK( N+J ),
00276      $                   WORK( 2*N+1 ) )
00277 *
00278       END IF
00279 *
00280       WORK( 1 ) = IWS
00281       RETURN
00282 *
00283 *     End of SGEQP3
00284 *
00285       END