|
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 177 178 179 180 |
SUBROUTINE CLAQP2( 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 VN1( * ), VN2( * ) COMPLEX A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * CLAQP2 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) COMPLEX 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) COMPLEX 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) COMPLEX 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 COMPLEX CONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, $ CONE = ( 1.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. INTEGER I, ITEMP, J, MN, OFFPI, PVT REAL TEMP, TEMP2, TOL3Z COMPLEX AII * .. * .. External Subroutines .. EXTERNAL CLARF, CLARFG, CSWAP * .. * .. Intrinsic Functions .. INTRINSIC ABS, CONJG, MAX, MIN, SQRT * .. * .. External Functions .. INTEGER ISAMAX REAL SCNRM2, SLAMCH EXTERNAL ISAMAX, SCNRM2, SLAMCH * .. * .. 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 CSWAP( 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 CLARFG( M-OFFPI+1, A( OFFPI, I ), A( OFFPI+1, I ), 1, $ TAU( I ) ) ELSE CALL CLARFG( 1, A( M, I ), A( M, I ), 1, TAU( I ) ) END IF * IF( I.LT.N ) THEN * * Apply H(i)**H to A(offset+i:m,i+1:n) from the left. * AII = A( OFFPI, I ) A( OFFPI, I ) = CONE CALL CLARF( 'Left', M-OFFPI+1, N-I, A( OFFPI, I ), 1, $ CONJG( 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 ) = SCNRM2( 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 CLAQP2 * END |