LAPACK: SRC/sgeqpf.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
00002 *
00003 *  -- LAPACK deprecated computational routine (version 3.2.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 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            INFO, LDA, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            JPVT( * )
00013       REAL               A( LDA, * ), TAU( * ), WORK( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  This routine is deprecated and has been replaced by routine SGEQP3.
00020 *
00021 *  SGEQPF computes a QR factorization with column pivoting of a
00022 *  real M-by-N matrix A: A*P = Q*R.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  M       (input) INTEGER
00028 *          The number of rows of the matrix A. M >= 0.
00029 *
00030 *  N       (input) INTEGER
00031 *          The number of columns of the matrix A. N >= 0
00032 *
00033 *  A       (input/output) REAL array, dimension (LDA,N)
00034 *          On entry, the M-by-N matrix A.
00035 *          On exit, the upper triangle of the array contains the
00036 *          min(M,N)-by-N upper triangular matrix R; the elements
00037 *          below the diagonal, together with the array TAU,
00038 *          represent the orthogonal matrix Q as a product of
00039 *          min(m,n) elementary reflectors.
00040 *
00041 *  LDA     (input) INTEGER
00042 *          The leading dimension of the array A. LDA >= max(1,M).
00043 *
00044 *  JPVT    (input/output) INTEGER array, dimension (N)
00045 *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00046 *          to the front of A*P (a leading column); if JPVT(i) = 0,
00047 *          the i-th column of A is a free column.
00048 *          On exit, if JPVT(i) = k, then the i-th column of A*P
00049 *          was the k-th column of A.
00050 *
00051 *  TAU     (output) REAL array, dimension (min(M,N))
00052 *          The scalar factors of the elementary reflectors.
00053 *
00054 *  WORK    (workspace) REAL array, dimension (3*N)
00055 *
00056 *  INFO    (output) INTEGER
00057 *          = 0:  successful exit
00058 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00059 *
00060 *  Further Details
00061 *  ===============
00062 *
00063 *  The matrix Q is represented as a product of elementary reflectors
00064 *
00065 *     Q = H(1) H(2) . . . H(n)
00066 *
00067 *  Each H(i) has the form
00068 *
00069 *     H = I - tau * v * v'
00070 *
00071 *  where tau is a real scalar, and v is a real vector with
00072 *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
00073 *
00074 *  The matrix P is represented in jpvt as follows: If
00075 *     jpvt(j) = i
00076 *  then the jth column of P is the ith canonical unit vector.
00077 *
00078 *  Partial column norm updating strategy modified by
00079 *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00080 *    University of Zagreb, Croatia.
00081 *     June 2010
00082 *  For more details see LAPACK Working Note 176.
00083 *
00084 *  =====================================================================
00085 *
00086 *     .. Parameters ..
00087       REAL               ZERO, ONE
00088       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00089 *     ..
00090 *     .. Local Scalars ..
00091       INTEGER            I, ITEMP, J, MA, MN, PVT
00092       REAL               AII, TEMP, TEMP2, TOL3Z
00093 *     ..
00094 *     .. External Subroutines ..
00095       EXTERNAL           SGEQR2, SLARF, SLARFG, SORM2R, SSWAP, XERBLA
00096 *     ..
00097 *     .. Intrinsic Functions ..
00098       INTRINSIC          ABS, MAX, MIN, SQRT
00099 *     ..
00100 *     .. External Functions ..
00101       INTEGER            ISAMAX
00102       REAL               SLAMCH, SNRM2
00103       EXTERNAL           ISAMAX, SLAMCH, SNRM2
00104 *     ..
00105 *     .. Executable Statements ..
00106 *
00107 *     Test the input arguments
00108 *
00109       INFO = 0
00110       IF( M.LT.0 ) THEN
00111          INFO = -1
00112       ELSE IF( N.LT.0 ) THEN
00113          INFO = -2
00114       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00115          INFO = -4
00116       END IF
00117       IF( INFO.NE.0 ) THEN
00118          CALL XERBLA( 'SGEQPF', -INFO )
00119          RETURN
00120       END IF
00121 *
00122       MN = MIN( M, N )
00123       TOL3Z = SQRT(SLAMCH('Epsilon'))
00124 *
00125 *     Move initial columns up front
00126 *
00127       ITEMP = 1
00128       DO 10 I = 1, N
00129          IF( JPVT( I ).NE.0 ) THEN
00130             IF( I.NE.ITEMP ) THEN
00131                CALL SSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
00132                JPVT( I ) = JPVT( ITEMP )
00133                JPVT( ITEMP ) = I
00134             ELSE
00135                JPVT( I ) = I
00136             END IF
00137             ITEMP = ITEMP + 1
00138          ELSE
00139             JPVT( I ) = I
00140          END IF
00141    10 CONTINUE
00142       ITEMP = ITEMP - 1
00143 *
00144 *     Compute the QR factorization and update remaining columns
00145 *
00146       IF( ITEMP.GT.0 ) THEN
00147          MA = MIN( ITEMP, M )
00148          CALL SGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
00149          IF( MA.LT.N ) THEN
00150             CALL SORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU,
00151      $                   A( 1, MA+1 ), LDA, WORK, INFO )
00152          END IF
00153       END IF
00154 *
00155       IF( ITEMP.LT.MN ) THEN
00156 *
00157 *        Initialize partial column norms. The first n elements of
00158 *        work store the exact column norms.
00159 *
00160          DO 20 I = ITEMP + 1, N
00161             WORK( I ) = SNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
00162             WORK( N+I ) = WORK( I )
00163    20    CONTINUE
00164 *
00165 *        Compute factorization
00166 *
00167          DO 40 I = ITEMP + 1, MN
00168 *
00169 *           Determine ith pivot column and swap if necessary
00170 *
00171             PVT = ( I-1 ) + ISAMAX( N-I+1, WORK( I ), 1 )
00172 *
00173             IF( PVT.NE.I ) THEN
00174                CALL SSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00175                ITEMP = JPVT( PVT )
00176                JPVT( PVT ) = JPVT( I )
00177                JPVT( I ) = ITEMP
00178                WORK( PVT ) = WORK( I )
00179                WORK( N+PVT ) = WORK( N+I )
00180             END IF
00181 *
00182 *           Generate elementary reflector H(i)
00183 *
00184             IF( I.LT.M ) THEN
00185                CALL SLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
00186             ELSE
00187                CALL SLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
00188             END IF
00189 *
00190             IF( I.LT.N ) THEN
00191 *
00192 *              Apply H(i) to A(i:m,i+1:n) from the left
00193 *
00194                AII = A( I, I )
00195                A( I, I ) = ONE
00196                CALL SLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ),
00197      $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
00198                A( I, I ) = AII
00199             END IF
00200 *
00201 *           Update partial column norms
00202 *
00203             DO 30 J = I + 1, N
00204                IF( WORK( J ).NE.ZERO ) THEN
00205 *
00206 *                 NOTE: The following 4 lines follow from the analysis in
00207 *                 Lapack Working Note 176.
00208 *                 
00209                   TEMP = ABS( A( I, J ) ) / WORK( J )
00210                   TEMP = MAX( ZERO, ( ONE+TEMP )*( ONE-TEMP ) )
00211                   TEMP2 = TEMP*( WORK( J ) / WORK( N+J ) )**2
00212                   IF( TEMP2 .LE. TOL3Z ) THEN 
00213                      IF( M-I.GT.0 ) THEN
00214                         WORK( J ) = SNRM2( M-I, A( I+1, J ), 1 )
00215                         WORK( N+J ) = WORK( J )
00216                      ELSE
00217                         WORK( J ) = ZERO
00218                         WORK( N+J ) = ZERO
00219                      END IF
00220                   ELSE
00221                      WORK( J ) = WORK( J )*SQRT( TEMP )
00222                   END IF
00223                END IF
00224    30       CONTINUE
00225 *
00226    40    CONTINUE
00227       END IF
00228       RETURN
00229 *
00230 *     End of SGEQPF
00231 *
00232       END