LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaqp2.f
Go to the documentation of this file.
00001 *> \brief \b DLAQP2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAQP2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqp2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqp2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqp2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
00022 *                          WORK )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            LDA, M, N, OFFSET
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            JPVT( * )
00029 *       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
00030 *      $                   WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLAQP2 computes a QR factorization with column pivoting of
00040 *> the block A(OFFSET+1:M,1:N).
00041 *> The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
00042 *> \endverbatim
00043 *
00044 *  Arguments:
00045 *  ==========
00046 *
00047 *> \param[in] M
00048 *> \verbatim
00049 *>          M is INTEGER
00050 *>          The number of rows of the matrix A. M >= 0.
00051 *> \endverbatim
00052 *>
00053 *> \param[in] N
00054 *> \verbatim
00055 *>          N is INTEGER
00056 *>          The number of columns of the matrix A. N >= 0.
00057 *> \endverbatim
00058 *>
00059 *> \param[in] OFFSET
00060 *> \verbatim
00061 *>          OFFSET is INTEGER
00062 *>          The number of rows of the matrix A that must be pivoted
00063 *>          but no factorized. OFFSET >= 0.
00064 *> \endverbatim
00065 *>
00066 *> \param[in,out] A
00067 *> \verbatim
00068 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00069 *>          On entry, the M-by-N matrix A.
00070 *>          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is 
00071 *>          the triangular factor obtained; the elements in block
00072 *>          A(OFFSET+1:M,1:N) below the diagonal, together with the
00073 *>          array TAU, represent the orthogonal matrix Q as a product of
00074 *>          elementary reflectors. Block A(1:OFFSET,1:N) has been
00075 *>          accordingly pivoted, but no factorized.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDA
00079 *> \verbatim
00080 *>          LDA is INTEGER
00081 *>          The leading dimension of the array A. LDA >= max(1,M).
00082 *> \endverbatim
00083 *>
00084 *> \param[in,out] JPVT
00085 *> \verbatim
00086 *>          JPVT is INTEGER array, dimension (N)
00087 *>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00088 *>          to the front of A*P (a leading column); if JPVT(i) = 0,
00089 *>          the i-th column of A is a free column.
00090 *>          On exit, if JPVT(i) = k, then the i-th column of A*P
00091 *>          was the k-th column of A.
00092 *> \endverbatim
00093 *>
00094 *> \param[out] TAU
00095 *> \verbatim
00096 *>          TAU is DOUBLE PRECISION array, dimension (min(M,N))
00097 *>          The scalar factors of the elementary reflectors.
00098 *> \endverbatim
00099 *>
00100 *> \param[in,out] VN1
00101 *> \verbatim
00102 *>          VN1 is DOUBLE PRECISION array, dimension (N)
00103 *>          The vector with the partial column norms.
00104 *> \endverbatim
00105 *>
00106 *> \param[in,out] VN2
00107 *> \verbatim
00108 *>          VN2 is DOUBLE PRECISION array, dimension (N)
00109 *>          The vector with the exact column norms.
00110 *> \endverbatim
00111 *>
00112 *> \param[out] WORK
00113 *> \verbatim
00114 *>          WORK is DOUBLE PRECISION array, dimension (N)
00115 *> \endverbatim
00116 *
00117 *  Authors:
00118 *  ========
00119 *
00120 *> \author Univ. of Tennessee 
00121 *> \author Univ. of California Berkeley 
00122 *> \author Univ. of Colorado Denver 
00123 *> \author NAG Ltd. 
00124 *
00125 *> \date November 2011
00126 *
00127 *> \ingroup doubleOTHERauxiliary
00128 *
00129 *> \par Contributors:
00130 *  ==================
00131 *>
00132 *>    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00133 *>    X. Sun, Computer Science Dept., Duke University, USA
00134 *> \n
00135 *>  Partial column norm updating strategy modified on April 2011
00136 *>    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00137 *>    University of Zagreb, Croatia.
00138 *
00139 *> \par References:
00140 *  ================
00141 *>
00142 *> LAPACK Working Note 176
00143 *
00144 *> \htmlonly
00145 *> <a href="http://www.netlib.org/lapack/lawnspdf/lawn176.pdf">[PDF]</a> 
00146 *> \endhtmlonly 
00147 *
00148 *  =====================================================================
00149       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
00150      $                   WORK )
00151 *
00152 *  -- LAPACK auxiliary routine (version 3.4.0) --
00153 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00154 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00155 *     November 2011
00156 *
00157 *     .. Scalar Arguments ..
00158       INTEGER            LDA, M, N, OFFSET
00159 *     ..
00160 *     .. Array Arguments ..
00161       INTEGER            JPVT( * )
00162       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
00163      $                   WORK( * )
00164 *     ..
00165 *
00166 *  =====================================================================
00167 *
00168 *     .. Parameters ..
00169       DOUBLE PRECISION   ZERO, ONE
00170       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00171 *     ..
00172 *     .. Local Scalars ..
00173       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
00174       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
00175 *     ..
00176 *     .. External Subroutines ..
00177       EXTERNAL           DLARF, DLARFG, DSWAP
00178 *     ..
00179 *     .. Intrinsic Functions ..
00180       INTRINSIC          ABS, MAX, MIN, SQRT
00181 *     ..
00182 *     .. External Functions ..
00183       INTEGER            IDAMAX
00184       DOUBLE PRECISION   DLAMCH, DNRM2
00185       EXTERNAL           IDAMAX, DLAMCH, DNRM2
00186 *     ..
00187 *     .. Executable Statements ..
00188 *
00189       MN = MIN( M-OFFSET, N )
00190       TOL3Z = SQRT(DLAMCH('Epsilon'))
00191 *
00192 *     Compute factorization.
00193 *
00194       DO 20 I = 1, MN
00195 *
00196          OFFPI = OFFSET + I
00197 *
00198 *        Determine ith pivot column and swap if necessary.
00199 *
00200          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
00201 *
00202          IF( PVT.NE.I ) THEN
00203             CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00204             ITEMP = JPVT( PVT )
00205             JPVT( PVT ) = JPVT( I )
00206             JPVT( I ) = ITEMP
00207             VN1( PVT ) = VN1( I )
00208             VN2( PVT ) = VN2( I )
00209          END IF
00210 *
00211 *        Generate elementary reflector H(i).
00212 *
00213          IF( OFFPI.LT.M ) THEN
00214             CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
00215      $                   TAU( I ) )
00216          ELSE
00217             CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
00218          END IF
00219 *
00220          IF( I.LE.N ) THEN
00221 *
00222 *           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
00223 *
00224             AII = A( OFFPI, I )
00225             A( OFFPI, I ) = ONE
00226             CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
00227      $                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
00228             A( OFFPI, I ) = AII
00229          END IF
00230 *
00231 *        Update partial column norms.
00232 *
00233          DO 10 J = I + 1, N
00234             IF( VN1( J ).NE.ZERO ) THEN
00235 *
00236 *              NOTE: The following 4 lines follow from the analysis in
00237 *              Lapack Working Note 176.
00238 *
00239                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
00240                TEMP = MAX( TEMP, ZERO )
00241                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
00242                IF( TEMP2 .LE. TOL3Z ) THEN
00243                   IF( OFFPI.LT.M ) THEN
00244                      VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
00245                      VN2( J ) = VN1( J )
00246                   ELSE
00247                      VN1( J ) = ZERO
00248                      VN2( J ) = ZERO
00249                   END IF
00250                ELSE
00251                   VN1( J ) = VN1( J )*SQRT( TEMP )
00252                END IF
00253             END IF
00254    10    CONTINUE
00255 *
00256    20 CONTINUE
00257 *
00258       RETURN
00259 *
00260 *     End of DLAQP2
00261 *
00262       END
 All Files Functions