|
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 |
SUBROUTINE CTZRQF( M, N, A, LDA, TAU, INFO )
* * -- 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 INFO, LDA, M, N * .. * .. Array Arguments .. COMPLEX A( LDA, * ), TAU( * ) * .. * * Purpose * ======= * * This routine is deprecated and has been replaced by routine CTZRZF. * * CTZRQF reduces the M-by-N ( M<=N ) complex upper trapezoidal matrix A * to upper triangular form by means of unitary transformations. * * The upper trapezoidal matrix A is factored as * * A = ( R 0 ) * Z, * * where Z is an N-by-N unitary matrix and R is an M-by-M upper * triangular matrix. * * 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 >= M. * * A (input/output) COMPLEX 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 M+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 array, dimension (M) * The scalar factors of the elementary reflectors. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * * Further Details * =============== * * The factorization is obtained by Householder's method. The kth * transformation matrix, Z( k ), whose conjugate transpose 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 ( n - m ) element vector. * tau and z( k ) are chosen to annihilate the elements of the kth row * of X. * * The scalar tau is returned in the kth element of TAU and the vector * u( k ) in the kth row of A, such that the elements of z( k ) are * in a( k, m + 1 ), ..., a( k, n ). The elements of R are returned in * the upper triangular part of A. * * Z is given by * * Z = Z( 1 ) * Z( 2 ) * ... * Z( m ). * * ===================================================================== * * .. Parameters .. COMPLEX CONE, CZERO PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ), $ CZERO = ( 0.0E+0, 0.0E+0 ) ) * .. * .. Local Scalars .. INTEGER I, K, M1 COMPLEX ALPHA * .. * .. Intrinsic Functions .. INTRINSIC CONJG, MAX, MIN * .. * .. External Subroutines .. EXTERNAL CAXPY, CCOPY, CGEMV, CGERC, CLACGV, CLARFG, $ XERBLA * .. * .. Executable Statements .. * * Test the input parameters. * INFO = 0 IF( M.LT.0 ) THEN INFO = -1 ELSE IF( N.LT.M ) THEN INFO = -2 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN INFO = -4 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'CTZRQF', -INFO ) RETURN END IF * * Perform the factorization. * IF( M.EQ.0 ) $ RETURN IF( M.EQ.N ) THEN DO 10 I = 1, N TAU( I ) = CZERO 10 CONTINUE ELSE M1 = MIN( M+1, N ) DO 20 K = M, 1, -1 * * Use a Householder reflection to zero the kth row of A. * First set up the reflection. * A( K, K ) = CONJG( A( K, K ) ) CALL CLACGV( N-M, A( K, M1 ), LDA ) ALPHA = A( K, K ) CALL CLARFG( N-M+1, ALPHA, A( K, M1 ), LDA, TAU( K ) ) A( K, K ) = ALPHA TAU( K ) = CONJG( TAU( K ) ) * IF( TAU( K ).NE.CZERO .AND. K.GT.1 ) THEN * * We now perform the operation A := A*P( k )**H. * * Use the first ( k - 1 ) elements of TAU to store a( k ), * where a( k ) consists of the first ( k - 1 ) elements of * the kth column of A. Also let B denote the first * ( k - 1 ) rows of the last ( n - m ) columns of A. * CALL CCOPY( K-1, A( 1, K ), 1, TAU, 1 ) * * Form w = a( k ) + B*z( k ) in TAU. * CALL CGEMV( 'No transpose', K-1, N-M, CONE, A( 1, M1 ), $ LDA, A( K, M1 ), LDA, CONE, TAU, 1 ) * * Now form a( k ) := a( k ) - conjg(tau)*w * and B := B - conjg(tau)*w*z( k )**H. * CALL CAXPY( K-1, -CONJG( TAU( K ) ), TAU, 1, A( 1, K ), $ 1 ) CALL CGERC( K-1, N-M, -CONJG( TAU( K ) ), TAU, 1, $ A( K, M1 ), LDA, A( 1, M1 ), LDA ) END IF 20 CONTINUE END IF * RETURN * * End of CTZRQF * END |