SUBROUTINE DELCOLS( M, N, A, LDA, K, P, TAU, WORK, INFO ) * * Craig Lucas, University of Manchester * March, 2004 * * .. Scalar Arguments .. INTEGER INFO, K, LDA, M, N, P * .. * .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * Given a real m by (n+p) matrix, B, and the QR factorization * B = Q_B * R_B, DELCOLS computes the QR factorization * C = Q * R where C is the matrix B with p columns deleted * from the kth column onwards. * * The input to this routine is Q_B' * C * * Arguments * ========= * * M (input) INTEGER * The number of rows of the matrix C. M >= 0. * * N (input) INTEGER * The number of columns of the matrix C. N >= 0. * * A (input/output) DOUBLE PRECISION array, dimension (LDA,N) * On entry, the matrix Q_B' * C. The elements in columns * 1:K-1 are not referenced. * * On exit, the elements on and above the diagonal contain * the n by n upper triangular part of the matrix R. The * elements below the diagonal in columns k:n, together with * TAU represent the orthogonal matrix Q as a product of * elementary reflectors (see Further Details). * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * K (input) INTEGER * The position of the first column deleted from B. * 0 < K <= N+P. * * P (input) INTEGER * The number of columns deleted from B. P > 0. * * TAU (output) DOUBLE PRECISION array, dimension(N-K+1) * The scalar factors of the elementary reflectors * (see Further Details). * * WORK DOUBLE PRECISION array, dimension (P+1) * Work space. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -I, the I-th argument had an illegal value. * * Further Details * =============== * * The matrix Q is represented as a product of Q_B and elementary * reflectors * * Q = Q_B * H(k) * H(k+1) *...* H(last), last = min( m-1, n ). * * Each H(j) has the form * * H(j) = I - tau*v*v' * * where tau is a real scalar, and v is a real vector with * v(1:j-1) = 0, v(j) = 1, v(j+1:j+lenh-1), lenh = min( p+1, m-j+1 ), * stored on exit in A(j+1:j+lenh-1,j) and v(j+lenh:m) = 0, tau is * stored in TAU(j). * * The matrix Q can be formed with DELCOLSQ * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) * .. * .. Local Scalars .. DOUBLE PRECISION AJJ INTEGER J, LAST, LENH * .. * .. External Subroutines .. EXTERNAL DLARF, DLARFG, XERBLA * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.0 ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 ELSE IF( K.GT.N+P .OR. K.LE.0 ) THEN INFO = -5 ELSE IF( P.LE.0 ) THEN INFO = -6 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'DELCOLS', -INFO ) RETURN END IF * LAST = MIN( M-1, N ) * DO 10 J = K, LAST * * Generate elementary reflector H(J) to annihilate the nonzero * entries below A(J,J) * LENH = MIN( P+1, M-J+1 ) CALL DLARFG( LENH, A( J, J ), A( J+1, J ), 1, TAU( J-K+1 ) ) * IF( J.LT.N ) THEN * * Apply H(J) to trailing matrix from left * AJJ = A( J, J ) A( J, J ) = ONE CALL DLARF( 'L', LENH, N-J, A( J, J ), 1, TAU( J-K+1 ), \$ A( J, J+1 ), LDA, WORK ) A( J, J ) = AJJ * END IF * 10 CONTINUE * RETURN * * End of DELCOLS * END