LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgeqpf.f
Go to the documentation of this file.
00001 *> \brief \b SGEQPF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGEQPF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeqpf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeqpf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeqpf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       INTEGER            INFO, LDA, 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 *> This routine is deprecated and has been replaced by routine SGEQP3.
00038 *>
00039 *> SGEQPF computes a QR factorization with column pivoting of a
00040 *> real M-by-N matrix A: A*P = Q*R.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] M
00047 *> \verbatim
00048 *>          M is INTEGER
00049 *>          The number of rows of the matrix A. M >= 0.
00050 *> \endverbatim
00051 *>
00052 *> \param[in] N
00053 *> \verbatim
00054 *>          N is INTEGER
00055 *>          The number of columns of the matrix A. N >= 0
00056 *> \endverbatim
00057 *>
00058 *> \param[in,out] A
00059 *> \verbatim
00060 *>          A is REAL array, dimension (LDA,N)
00061 *>          On entry, the M-by-N matrix A.
00062 *>          On exit, the upper triangle of the array contains the
00063 *>          min(M,N)-by-N upper triangular matrix R; the elements
00064 *>          below the diagonal, together with the array TAU,
00065 *>          represent the orthogonal matrix Q as a product of
00066 *>          min(m,n) elementary reflectors.
00067 *> \endverbatim
00068 *>
00069 *> \param[in] LDA
00070 *> \verbatim
00071 *>          LDA is INTEGER
00072 *>          The leading dimension of the array A. LDA >= max(1,M).
00073 *> \endverbatim
00074 *>
00075 *> \param[in,out] JPVT
00076 *> \verbatim
00077 *>          JPVT is INTEGER array, dimension (N)
00078 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00079 *>          to the front of A*P (a leading column); if JPVT(i) = 0,
00080 *>          the i-th column of A is a free column.
00081 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
00082 *>          was the k-th column of A.
00083 *> \endverbatim
00084 *>
00085 *> \param[out] TAU
00086 *> \verbatim
00087 *>          TAU is REAL array, dimension (min(M,N))
00088 *>          The scalar factors of the elementary reflectors.
00089 *> \endverbatim
00090 *>
00091 *> \param[out] WORK
00092 *> \verbatim
00093 *>          WORK is REAL array, dimension (3*N)
00094 *> \endverbatim
00095 *>
00096 *> \param[out] INFO
00097 *> \verbatim
00098 *>          INFO is INTEGER
00099 *>          = 0:  successful exit
00100 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00101 *> \endverbatim
00102 *
00103 *  Authors:
00104 *  ========
00105 *
00106 *> \author Univ. of Tennessee 
00107 *> \author Univ. of California Berkeley 
00108 *> \author Univ. of Colorado Denver 
00109 *> \author NAG Ltd. 
00110 *
00111 *> \date November 2011
00112 *
00113 *> \ingroup realGEcomputational
00114 *
00115 *> \par Further Details:
00116 *  =====================
00117 *>
00118 *> \verbatim
00119 *>
00120 *>  The matrix Q is represented as a product of elementary reflectors
00121 *>
00122 *>     Q = H(1) H(2) . . . H(n)
00123 *>
00124 *>  Each H(i) has the form
00125 *>
00126 *>     H = I - tau * v * v**T
00127 *>
00128 *>  where tau is a real scalar, and v is a real vector with
00129 *>  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
00130 *>
00131 *>  The matrix P is represented in jpvt as follows: If
00132 *>     jpvt(j) = i
00133 *>  then the jth column of P is the ith canonical unit vector.
00134 *>
00135 *>  Partial column norm updating strategy modified by
00136 *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00137 *>    University of Zagreb, Croatia.
00138 *>  -- April 2011                                                      --
00139 *>  For more details see LAPACK Working Note 176.
00140 *> \endverbatim
00141 *>
00142 *  =====================================================================
00143       SUBROUTINE SGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
00144 *
00145 *  -- LAPACK computational routine (version 3.4.0) --
00146 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00147 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00148 *     November 2011
00149 *
00150 *     .. Scalar Arguments ..
00151       INTEGER            INFO, LDA, M, N
00152 *     ..
00153 *     .. Array Arguments ..
00154       INTEGER            JPVT( * )
00155       REAL               A( LDA, * ), TAU( * ), WORK( * )
00156 *     ..
00157 *
00158 *  =====================================================================
00159 *
00160 *     .. Parameters ..
00161       REAL               ZERO, ONE
00162       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       INTEGER            I, ITEMP, J, MA, MN, PVT
00166       REAL               AII, TEMP, TEMP2, TOL3Z
00167 *     ..
00168 *     .. External Subroutines ..
00169       EXTERNAL           SGEQR2, SLARF, SLARFG, SORM2R, SSWAP, XERBLA
00170 *     ..
00171 *     .. Intrinsic Functions ..
00172       INTRINSIC          ABS, MAX, MIN, SQRT
00173 *     ..
00174 *     .. External Functions ..
00175       INTEGER            ISAMAX
00176       REAL               SLAMCH, SNRM2
00177       EXTERNAL           ISAMAX, SLAMCH, SNRM2
00178 *     ..
00179 *     .. Executable Statements ..
00180 *
00181 *     Test the input arguments
00182 *
00183       INFO = 0
00184       IF( M.LT.0 ) THEN
00185          INFO = -1
00186       ELSE IF( N.LT.0 ) THEN
00187          INFO = -2
00188       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00189          INFO = -4
00190       END IF
00191       IF( INFO.NE.0 ) THEN
00192          CALL XERBLA( 'SGEQPF', -INFO )
00193          RETURN
00194       END IF
00195 *
00196       MN = MIN( M, N )
00197       TOL3Z = SQRT(SLAMCH('Epsilon'))
00198 *
00199 *     Move initial columns up front
00200 *
00201       ITEMP = 1
00202       DO 10 I = 1, N
00203          IF( JPVT( I ).NE.0 ) THEN
00204             IF( I.NE.ITEMP ) THEN
00205                CALL SSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
00206                JPVT( I ) = JPVT( ITEMP )
00207                JPVT( ITEMP ) = I
00208             ELSE
00209                JPVT( I ) = I
00210             END IF
00211             ITEMP = ITEMP + 1
00212          ELSE
00213             JPVT( I ) = I
00214          END IF
00215    10 CONTINUE
00216       ITEMP = ITEMP - 1
00217 *
00218 *     Compute the QR factorization and update remaining columns
00219 *
00220       IF( ITEMP.GT.0 ) THEN
00221          MA = MIN( ITEMP, M )
00222          CALL SGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
00223          IF( MA.LT.N ) THEN
00224             CALL SORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU,
00225      $                   A( 1, MA+1 ), LDA, WORK, INFO )
00226          END IF
00227       END IF
00228 *
00229       IF( ITEMP.LT.MN ) THEN
00230 *
00231 *        Initialize partial column norms. The first n elements of
00232 *        work store the exact column norms.
00233 *
00234          DO 20 I = ITEMP + 1, N
00235             WORK( I ) = SNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
00236             WORK( N+I ) = WORK( I )
00237    20    CONTINUE
00238 *
00239 *        Compute factorization
00240 *
00241          DO 40 I = ITEMP + 1, MN
00242 *
00243 *           Determine ith pivot column and swap if necessary
00244 *
00245             PVT = ( I-1 ) + ISAMAX( N-I+1, WORK( I ), 1 )
00246 *
00247             IF( PVT.NE.I ) THEN
00248                CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00249                ITEMP = JPVT( PVT )
00250                JPVT( PVT ) = JPVT( I )
00251                JPVT( I ) = ITEMP
00252                WORK( PVT ) = WORK( I )
00253                WORK( N+PVT ) = WORK( N+I )
00254             END IF
00255 *
00256 *           Generate elementary reflector H(i)
00257 *
00258             IF( I.LT.M ) THEN
00259                CALL SLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
00260             ELSE
00261                CALL SLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
00262             END IF
00263 *
00264             IF( I.LT.N ) THEN
00265 *
00266 *              Apply H(i) to A(i:m,i+1:n) from the left
00267 *
00268                AII = A( I, I )
00269                A( I, I ) = ONE
00270                CALL SLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ),
00271      $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
00272                A( I, I ) = AII
00273             END IF
00274 *
00275 *           Update partial column norms
00276 *
00277             DO 30 J = I + 1, N
00278                IF( WORK( J ).NE.ZERO ) THEN
00279 *
00280 *                 NOTE: The following 4 lines follow from the analysis in
00281 *                 Lapack Working Note 176.
00282 *                 
00283                   TEMP = ABS( A( I, J ) ) / WORK( J )
00284                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
00285                   TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2
00286                   IF( TEMP2 .LE. TOL3Z ) THEN 
00287                      IF( M-I.GT.0 ) THEN
00288                         WORK( J ) = SNRM2( M-I, A( I+1, J ), 1 )
00289                         WORK( N+J ) = WORK( J )
00290                      ELSE
00291                         WORK( J ) = ZERO
00292                         WORK( N+J ) = ZERO
00293                      END IF
00294                   ELSE
00295                      WORK( J ) = WORK( J )*SQRT( TEMP )
00296                   END IF
00297                END IF
00298    30       CONTINUE
00299 *
00300    40    CONTINUE
00301       END IF
00302       RETURN
00303 *
00304 *     End of SGEQPF
00305 *
00306       END
 All Files Functions