|
1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 |
SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2,
$ WORK ) * * -- LAPACK auxiliary routine (version 3.3.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * -- April 2011 -- * * .. Scalar Arguments .. INTEGER LDA, M, N, OFFSET * .. * .. Array Arguments .. INTEGER JPVT( * ) DOUBLE PRECISION A( LDA, * ), TAU( * ), VN1( * ), VN2( * ), $ WORK( * ) * .. * * Purpose * ======= * * DLAQP2 computes a QR factorization with column pivoting of * the block A(OFFSET+1:M,1:N). * The block A(1:OFFSET,1:N) is accordingly pivoted, but not factorized. * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix A. M >= 0. * * N (input) INTEGER * The number of columns of the matrix A. N >= 0. * * OFFSET (input) INTEGER * The number of rows of the matrix A that must be pivoted * but no factorized. OFFSET >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the M-by-N matrix A. * On exit, the upper triangle of block A(OFFSET+1:M,1:N) is * the triangular factor obtained; the elements in block * A(OFFSET+1:M,1:N) below the diagonal, together with the * array TAU, represent the orthogonal matrix Q as a product of * elementary reflectors. Block A(1:OFFSET,1:N) has been * accordingly pivoted, but no factorized. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * JPVT (input/output) INTEGER array, dimension (N) * On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted * to the front of A*P (a leading column); if JPVT(i) = 0, * the i-th column of A is a free column. * On exit, if JPVT(i) = k, then the i-th column of A*P * was the k-th column of A. * * TAU (output) DOUBLE PRECISION array, dimension (min(M,N)) * The scalar factors of the elementary reflectors. * * VN1 (input/output) DOUBLE PRECISION array, dimension (N) * The vector with the partial column norms. * * VN2 (input/output) DOUBLE PRECISION array, dimension (N) * The vector with the exact column norms. * * WORK (workspace) DOUBLE PRECISION array, dimension (N) * * Further Details * =============== * * Based on contributions by * G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain * X. Sun, Computer Science Dept., Duke University, USA * * Partial column norm updating strategy modified by * Z. Drmac and Z. Bujanovic, Dept. of Mathematics, * University of Zagreb, Croatia. * -- April 2011 -- * For more details see LAPACK Working Note 176. * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MN, OFFPI, PVT DOUBLE PRECISION AII, TEMP, TEMP2, TOL3Z * .. * .. External Subroutines .. EXTERNAL DLARF, DLARFG, DSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. External Functions .. INTEGER IDAMAX DOUBLE PRECISION DLAMCH, DNRM2 EXTERNAL IDAMAX, DLAMCH, DNRM2 * .. * .. Executable Statements .. * MN = MIN( M-OFFSET, N ) TOL3Z = SQRT(DLAMCH('Epsilon')) * * Compute factorization. * DO 20 I = 1, MN * OFFPI = OFFSET + I * * Determine ith pivot column and swap if necessary. * PVT = ( I-1 ) + IDAMAX( N-I+1, VN1( I ), 1 ) * IF( PVT.NE.I ) THEN CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 ) ITEMP = JPVT( PVT ) JPVT( PVT ) = JPVT( I ) JPVT( I ) = ITEMP VN1( PVT ) = VN1( I ) VN2( PVT ) = VN2( I ) END IF * * Generate elementary reflector H(i). * IF( OFFPI.LT.M ) THEN CALL DLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1, $ TAU( I ) ) ELSE CALL DLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) ) END IF * IF( I.LE.N ) THEN * * Apply H(i)**T to A(offset+i:m,i+1:n) from the left. * AII = A( OFFPI, I ) A( OFFPI, I ) = ONE CALL DLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, $ TAU( I ), A( OFFPI, I+1 ), LDA, WORK( 1 ) ) A( OFFPI, I ) = AII END IF * * Update partial column norms. * DO 10 J = I + 1, N IF( VN1( J ).NE.ZERO ) THEN * * NOTE: The following 4 lines follow from the analysis in * Lapack Working Note 176. * TEMP = ONE - ( ABS( A( OFFPI, J ) ) / VN1( J ) )**2 TEMP = MAX( TEMP, ZERO ) TEMP2 = TEMP*( VN1( J ) / VN2( J ) )**2 IF( TEMP2 .LE. TOL3Z ) THEN IF( OFFPI.LT.M ) THEN VN1( J ) = DNRM2( M-OFFPI, A( OFFPI+1, J ), 1 ) VN2( J ) = VN1( J ) ELSE VN1( J ) = ZERO VN2( J ) = ZERO END IF ELSE VN1( J ) = VN1( J )*SQRT( TEMP ) END IF END IF 10 CONTINUE * 20 CONTINUE * RETURN * * End of DLAQP2 * END |