|
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 |
SUBROUTINE ZLATRZ( M, N, L, A, LDA, TAU, WORK )
* * -- LAPACK 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 L, LDA, M, N * .. * .. Array Arguments .. COMPLEX*16 A( LDA, * ), TAU( * ), WORK( * ) * .. * * Purpose * ======= * * ZLATRZ factors the M-by-(M+L) complex upper trapezoidal matrix * [ A1 A2 ] = [ A(1:M,1:M) A(1:M,N-L+1:N) ] as ( R 0 ) * Z by means * of unitary transformations, where Z is an (M+L)-by-(M+L) unitary * matrix and, R and A1 are M-by-M upper triangular matrices. * * 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. * * L (input) INTEGER * The number of columns of the matrix A containing the * meaningful part of the Householder vectors. N-M >= L >= 0. * * A (input/output) COMPLEX*16 array, dimension (LDA,N) * On entry, the leading M-by-N upper trapezoidal part of the * array A must contain the matrix to be factorized. * On exit, the leading M-by-M upper triangular part of A * contains the upper triangular matrix R, and elements N-L+1 to * N of the first M rows of A, with the array TAU, represent the * unitary matrix Z as a product of M elementary reflectors. * * LDA (input) INTEGER * The leading dimension of the array A. LDA >= max(1,M). * * TAU (output) COMPLEX*16 array, dimension (M) * The scalar factors of the elementary reflectors. * * WORK (workspace) COMPLEX*16 array, dimension (M) * * Further Details * =============== * * Based on contributions by * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA * * The factorization is obtained by Householder's method. The kth * transformation matrix, Z( k ), which is used to introduce zeros into * the ( m - k + 1 )th row of A, is given in the form * * Z( k ) = ( I 0 ), * ( 0 T( k ) ) * * where * * T( k ) = I - tau*u( k )*u( k )**H, u( k ) = ( 1 ), * ( 0 ) * ( z( k ) ) * * tau is a scalar and z( k ) is an l element vector. tau and z( k ) * are chosen to annihilate the elements of the kth row of A2. * * The scalar tau is returned in the kth element of TAU and the vector * u( k ) in the kth row of A2, such that the elements of z( k ) are * in a( k, l + 1 ), ..., a( k, n ). The elements of R are returned in * the upper triangular part of A1. * * Z is given by * * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). * * ===================================================================== * * .. Parameters .. COMPLEX*16 ZERO PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) ) * .. * .. Local Scalars .. INTEGER I COMPLEX*16 ALPHA * .. * .. External Subroutines .. EXTERNAL ZLACGV, ZLARFG, ZLARZ * .. * .. Intrinsic Functions .. INTRINSIC DCONJG * .. * .. Executable Statements .. * * Quick return if possible * IF( M.EQ.0 ) THEN RETURN ELSE IF( M.EQ.N ) THEN DO 10 I = 1, N TAU( I ) = ZERO 10 CONTINUE RETURN END IF * DO 20 I = M, 1, -1 * * Generate elementary reflector H(i) to annihilate * [ A(i,i) A(i,n-l+1:n) ] * CALL ZLACGV( L, A( I, N-L+1 ), LDA ) ALPHA = DCONJG( A( I, I ) ) CALL ZLARFG( L+1, ALPHA, A( I, N-L+1 ), LDA, TAU( I ) ) TAU( I ) = DCONJG( TAU( I ) ) * * Apply H(i) to A(1:i-1,i:n) from the right * CALL ZLARZ( 'Right', I-1, N-I+1, L, A( I, N-L+1 ), LDA, $ DCONJG( TAU( I ) ), A( 1, I ), LDA, WORK ) A( I, I ) = DCONJG( ALPHA ) * 20 CONTINUE * RETURN * * End of ZLATRZ * END |