00001 SUBROUTINE DLAQP2( M, N, OFFSET, A, LDA, JPVT, TAU, VN1, VN2, 00002 $ WORK ) 00003 * 00004 * -- LAPACK auxiliary routine (version 3.2.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * June 2010 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 * June 2010 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)' 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
1.7.2