LAPACK 3.3.1 Linear Algebra PACKage

# dlaqp2.f

Go to the documentation of this file.
```00001       SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
00002      \$                   WORK )
00003 *
00004 *  -- LAPACK auxiliary routine (version 3.3.1) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *  -- April 2011                                                      --
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            LDA, M, N, OFFSET
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            JPVT( * )
00014       DOUBLE PRECISION   A( LDA, * ), TAU( * ), VN1( * ), VN2( * ),
00015      \$                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DLAQP2 computes a QR factorization with column pivoting of
00022 *  the block A(OFFSET+1:M,1:N).
00023 *  The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized.
00024 *
00025 *  Arguments
00026 *  =========
00027 *
00028 *  M       (input) INTEGER
00029 *          The number of rows of the matrix A. M >= 0.
00030 *
00031 *  N       (input) INTEGER
00032 *          The number of columns of the matrix A. N >= 0.
00033 *
00034 *  OFFSET  (input) INTEGER
00035 *          The number of rows of the matrix A that must be pivoted
00036 *          but no factorized. OFFSET >= 0.
00037 *
00038 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00039 *          On entry, the M-by-N matrix A.
00040 *          On exit, the upper triangle of block A(OFFSET+1:M,1:N) is
00041 *          the triangular factor obtained; the elements in block
00042 *          A(OFFSET+1:M,1:N) below the diagonal, together with the
00043 *          array TAU, represent the orthogonal matrix Q as a product of
00044 *          elementary reflectors. Block A(1:OFFSET,1:N) has been
00045 *          accordingly pivoted, but no factorized.
00046 *
00047 *  LDA     (input) INTEGER
00048 *          The leading dimension of the array A. LDA >= max(1,M).
00049 *
00050 *  JPVT    (input/output) INTEGER array, dimension (N)
00051 *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
00052 *          to the front of A*P (a leading column); if JPVT(i) = 0,
00053 *          the i-th column of A is a free column.
00054 *          On exit, if JPVT(i) = k, then the i-th column of A*P
00055 *          was the k-th column of A.
00056 *
00057 *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
00058 *          The scalar factors of the elementary reflectors.
00059 *
00060 *  VN1     (input/output) DOUBLE PRECISION array, dimension (N)
00061 *          The vector with the partial column norms.
00062 *
00063 *  VN2     (input/output) DOUBLE PRECISION array, dimension (N)
00064 *          The vector with the exact column norms.
00065 *
00066 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
00067 *
00068 *  Further Details
00069 *  ===============
00070 *
00071 *  Based on contributions by
00072 *    G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
00073 *    X. Sun, Computer Science Dept., Duke University, USA
00074 *
00075 *  Partial column norm updating strategy modified by
00076 *    Z. Drmac and Z. Bujanovic, Dept. of Mathematics,
00077 *    University of Zagreb, Croatia.
00078 *  -- April 2011                                                      --
00079 *  For more details see LAPACK Working Note 176.
00080 *  =====================================================================
00081 *
00082 *     .. Parameters ..
00083       DOUBLE PRECISION   ZERO, ONE
00084       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00085 *     ..
00086 *     .. Local Scalars ..
00087       INTEGER            I, ITEMP, J, MN, OFFPI, PVT
00088       DOUBLE PRECISION   AII, TEMP, TEMP2, TOL3Z
00089 *     ..
00090 *     .. External Subroutines ..
00091       EXTERNAL           DLARF, DLARFG, DSWAP
00092 *     ..
00093 *     .. Intrinsic Functions ..
00094       INTRINSIC          ABS, MAX, MIN, SQRT
00095 *     ..
00096 *     .. External Functions ..
00097       INTEGER            IDAMAX
00098       DOUBLE PRECISION   DLAMCH, DNRM2
00099       EXTERNAL           IDAMAX, DLAMCH, DNRM2
00100 *     ..
00101 *     .. Executable Statements ..
00102 *
00103       MN = MIN( M-OFFSET, N )
00104       TOL3Z = SQRT(DLAMCH('Epsilon'))
00105 *
00106 *     Compute factorization.
00107 *
00108       DO 20 I = 1, MN
00109 *
00110          OFFPI = OFFSET + I
00111 *
00112 *        Determine ith pivot column and swap if necessary.
00113 *
00114          PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 )
00115 *
00116          IF( PVT.NE.I ) THEN
00117             CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
00118             ITEMP = JPVT( PVT )
00119             JPVT( PVT ) = JPVT( I )
00120             JPVT( I ) = ITEMP
00121             VN1( PVT ) = VN1( I )
00122             VN2( PVT ) = VN2( I )
00123          END IF
00124 *
00125 *        Generate elementary reflector H(i).
00126 *
00127          IF( OFFPI.LT.M ) THEN
00128             CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1,
00129      \$                   TAU( I ) )
00130          ELSE
00131             CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) )
00132          END IF
00133 *
00134          IF( I.LE.N ) THEN
00135 *
00136 *           Apply H(i)**T to A(offset+i:m,i+1:n) from the left.
00137 *
00138             AII = A( OFFPI, I )
00139             A( OFFPI, I ) = ONE
00140             CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1,
00141      \$                  TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) )
00142             A( OFFPI, I ) = AII
00143          END IF
00144 *
00145 *        Update partial column norms.
00146 *
00147          DO 10 J = I + 1, N
00148             IF( VN1( J ).NE.ZERO ) THEN
00149 *
00150 *              NOTE: The following 4 lines follow from the analysis in
00151 *              Lapack Working Note 176.
00152 *
00153                TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2
00154                TEMP = MAX( TEMP, ZERO )
00155                TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2
00156                IF( TEMP2 .LE. TOL3Z ) THEN
00157                   IF( OFFPI.LT.M ) THEN
00158                      VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 )
00159                      VN2( J ) = VN1( J )
00160                   ELSE
00161                      VN1( J ) = ZERO
00162                      VN2( J ) = ZERO
00163                   END IF
00164                ELSE
00165                   VN1( J ) = VN1( J )*SQRT( TEMP )
00166                END IF
00167             END IF
00168    10    CONTINUE
00169 *
00170    20 CONTINUE
00171 *
00172       RETURN
00173 *
00174 *     End of DLAQP2
00175 *
00176       END
```