|
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 SLAQP2( 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( * ) REAL A( LDA, * ), TAU( * ), VN1( * ), VN2( * ), $ WORK( * ) * .. * * Purpose * ======= * * SLAQP2 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) REAL 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) REAL array, dimension (min(M,N)) * The scalar factors of the elementary reflectors. * * VN1 (input/output) REAL array, dimension (N) * The vector with the partial column norms. * * VN2 (input/output) REAL array, dimension (N) * The vector with the exact column norms. * * WORK (workspace) REAL 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 .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MN, OFFPI, PVT REAL AII, TEMP, TEMP2, TOL3Z * .. * .. External Subroutines .. EXTERNAL SLARF, SLARFG, SSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN, SQRT * .. * .. External Functions .. INTEGER ISAMAX REAL SLAMCH, SNRM2 EXTERNAL ISAMAX, SLAMCH, SNRM2 * .. * .. Executable Statements .. * MN = MIN( M-OFFSET, N ) TOL3Z = SQRT(SLAMCH('Epsilon')) * * Compute factorization. * DO 20 I = 1, MN * OFFPI = OFFSET + I * * Determine ith pivot column and swap if necessary. * PVT = ( I-1 ) + ISAMAX( N-I+1, VN1( I ), 1 ) * IF( PVT.NE.I ) THEN CALL SSWAP( 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 SLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1, $ TAU( I ) ) ELSE CALL SLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) ) END IF * IF( I.LT.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 SLARF( '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 ) = SNRM2( 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 SLAQP2 * END |