Mercurial > octave-nkf
diff libcruft/lapack/cungl2.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libcruft/lapack/cungl2.f Sun Apr 27 22:34:17 2008 +0200 @@ -0,0 +1,136 @@ + SUBROUTINE CUNGL2( M, N, K, A, LDA, TAU, WORK, INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER INFO, K, LDA, M, N +* .. +* .. Array Arguments .. + COMPLEX A( LDA, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* CUNGL2 generates an m-by-n complex matrix Q with orthonormal rows, +* which is defined as the first m rows of a product of k elementary +* reflectors of order n +* +* Q = H(k)' . . . H(2)' H(1)' +* +* as returned by CGELQF. +* +* Arguments +* ========= +* +* M (input) INTEGER +* The number of rows of the matrix Q. M >= 0. +* +* N (input) INTEGER +* The number of columns of the matrix Q. N >= M. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines the +* matrix Q. M >= K >= 0. +* +* A (input/output) COMPLEX array, dimension (LDA,N) +* On entry, the i-th row must contain the vector which defines +* the elementary reflector H(i), for i = 1,2,...,k, as returned +* by CGELQF in the first k rows of its array argument A. +* On exit, the m by n matrix Q. +* +* LDA (input) INTEGER +* The first dimension of the array A. LDA >= max(1,M). +* +* TAU (input) COMPLEX array, dimension (K) +* TAU(i) must contain the scalar factor of the elementary +* reflector H(i), as returned by CGELQF. +* +* WORK (workspace) COMPLEX array, dimension (M) +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument has an illegal value +* +* ===================================================================== +* +* .. Parameters .. + COMPLEX ONE, ZERO + PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), + $ ZERO = ( 0.0E+0, 0.0E+0 ) ) +* .. +* .. Local Scalars .. + INTEGER I, J, L +* .. +* .. External Subroutines .. + EXTERNAL CLACGV, CLARF, CSCAL, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC CONJG, MAX +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + IF( M.LT.0 ) THEN + INFO = -1 + ELSE IF( N.LT.M ) THEN + INFO = -2 + ELSE IF( K.LT.0 .OR. K.GT.M ) THEN + INFO = -3 + ELSE IF( LDA.LT.MAX( 1, M ) ) THEN + INFO = -5 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'CUNGL2', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.LE.0 ) + $ RETURN +* + IF( K.LT.M ) THEN +* +* Initialise rows k+1:m to rows of the unit matrix +* + DO 20 J = 1, N + DO 10 L = K + 1, M + A( L, J ) = ZERO + 10 CONTINUE + IF( J.GT.K .AND. J.LE.M ) + $ A( J, J ) = ONE + 20 CONTINUE + END IF +* + DO 40 I = K, 1, -1 +* +* Apply H(i)' to A(i:m,i:n) from the right +* + IF( I.LT.N ) THEN + CALL CLACGV( N-I, A( I, I+1 ), LDA ) + IF( I.LT.M ) THEN + A( I, I ) = ONE + CALL CLARF( 'Right', M-I, N-I+1, A( I, I ), LDA, + $ CONJG( TAU( I ) ), A( I+1, I ), LDA, WORK ) + END IF + CALL CSCAL( N-I, -TAU( I ), A( I, I+1 ), LDA ) + CALL CLACGV( N-I, A( I, I+1 ), LDA ) + END IF + A( I, I ) = ONE - CONJG( TAU( I ) ) +* +* Set A(i,1:i-1,i) to zero +* + DO 30 L = 1, I - 1 + A( I, L ) = ZERO + 30 CONTINUE + 40 CONTINUE + RETURN +* +* End of CUNGL2 +* + END