diff options
Diffstat (limited to '2.3-1/src/fortran/lapack/dlahr2.f')
| -rw-r--r-- | 2.3-1/src/fortran/lapack/dlahr2.f | 238 |
1 files changed, 238 insertions, 0 deletions
diff --git a/2.3-1/src/fortran/lapack/dlahr2.f b/2.3-1/src/fortran/lapack/dlahr2.f new file mode 100644 index 00000000..6af74977 --- /dev/null +++ b/2.3-1/src/fortran/lapack/dlahr2.f @@ -0,0 +1,238 @@ + SUBROUTINE DLAHR2( N, K, NB, A, LDA, TAU, T, LDT, Y, LDY ) +* +* -- LAPACK auxiliary routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + INTEGER K, LDA, LDT, LDY, N, NB +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), T( LDT, NB ), TAU( NB ), + $ Y( LDY, NB ) +* .. +* +* Purpose +* ======= +* +* DLAHR2 reduces the first NB columns of A real general n-BY-(n-k+1) +* matrix A so that elements below the k-th subdiagonal are zero. The +* reduction is performed by an orthogonal similarity transformation +* Q' * A * Q. The routine returns the matrices V and T which determine +* Q as a block reflector I - V*T*V', and also the matrix Y = A * V * T. +* +* This is an auxiliary routine called by DGEHRD. +* +* Arguments +* ========= +* +* N (input) INTEGER +* The order of the matrix A. +* +* K (input) INTEGER +* The offset for the reduction. Elements below the k-th +* subdiagonal in the first NB columns are reduced to zero. +* K < N. +* +* NB (input) INTEGER +* The number of columns to be reduced. +* +* A (input/output) DOUBLE PRECISION array, dimension (LDA,N-K+1) +* On entry, the n-by-(n-k+1) general matrix A. +* On exit, the elements on and above the k-th subdiagonal in +* the first NB columns are overwritten with the corresponding +* elements of the reduced matrix; the elements below the k-th +* subdiagonal, with the array TAU, represent the matrix Q as a +* product of elementary reflectors. The other columns of A are +* unchanged. See Further Details. +* +* LDA (input) INTEGER +* The leading dimension of the array A. LDA >= max(1,N). +* +* TAU (output) DOUBLE PRECISION array, dimension (NB) +* The scalar factors of the elementary reflectors. See Further +* Details. +* +* T (output) DOUBLE PRECISION array, dimension (LDT,NB) +* The upper triangular matrix T. +* +* LDT (input) INTEGER +* The leading dimension of the array T. LDT >= NB. +* +* Y (output) DOUBLE PRECISION array, dimension (LDY,NB) +* The n-by-nb matrix Y. +* +* LDY (input) INTEGER +* The leading dimension of the array Y. LDY >= N. +* +* Further Details +* =============== +* +* The matrix Q is represented as a product of nb elementary reflectors +* +* Q = H(1) H(2) . . . H(nb). +* +* Each H(i) has the form +* +* H(i) = I - tau * v * v' +* +* where tau is a real scalar, and v is a real vector with +* v(1:i+k-1) = 0, v(i+k) = 1; v(i+k+1:n) is stored on exit in +* A(i+k+1:n,i), and tau in TAU(i). +* +* The elements of the vectors v together form the (n-k+1)-by-nb matrix +* V which is needed, with T and Y, to apply the transformation to the +* unreduced part of the matrix, using an update of the form: +* A := (I - V*T*V') * (A - Y*V'). +* +* The contents of A on exit are illustrated by the following example +* with n = 7, k = 3 and nb = 2: +* +* ( a a a a a ) +* ( a a a a a ) +* ( a a a a a ) +* ( h h a a a ) +* ( v1 h a a a ) +* ( v1 v2 a a a ) +* ( v1 v2 a a a ) +* +* where a denotes an element of the original matrix A, h denotes a +* modified element of the upper Hessenberg matrix H, and vi denotes an +* element of the vector defining H(i). +* +* This file is a slight modification of LAPACK-3.0's DLAHRD +* incorporating improvements proposed by Quintana-Orti and Van de +* Gejin. Note that the entries of A(1:K,2:NB) differ from those +* returned by the original LAPACK routine. This function is +* not backward compatible with LAPACK3.0. +* +* ===================================================================== +* +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, + $ ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + INTEGER I + DOUBLE PRECISION EI +* .. +* .. External Subroutines .. + EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DLACPY, + $ DLARFG, DSCAL, DTRMM, DTRMV +* .. +* .. Intrinsic Functions .. + INTRINSIC MIN +* .. +* .. Executable Statements .. +* +* Quick return if possible +* + IF( N.LE.1 ) + $ RETURN +* + DO 10 I = 1, NB + IF( I.GT.1 ) THEN +* +* Update A(K+1:N,I) +* +* Update I-th column of A - Y * V' +* + CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, Y(K+1,1), LDY, + $ A( K+I-1, 1 ), LDA, ONE, A( K+1, I ), 1 ) +* +* Apply I - V * T' * V' to this column (call it b) from the +* left, using the last column of T as workspace +* +* Let V = ( V1 ) and b = ( b1 ) (first I-1 rows) +* ( V2 ) ( b2 ) +* +* where V1 is unit lower triangular +* +* w := V1' * b1 +* + CALL DCOPY( I-1, A( K+1, I ), 1, T( 1, NB ), 1 ) + CALL DTRMV( 'Lower', 'Transpose', 'UNIT', + $ I-1, A( K+1, 1 ), + $ LDA, T( 1, NB ), 1 ) +* +* w := w + V2'*b2 +* + CALL DGEMV( 'Transpose', N-K-I+1, I-1, + $ ONE, A( K+I, 1 ), + $ LDA, A( K+I, I ), 1, ONE, T( 1, NB ), 1 ) +* +* w := T'*w +* + CALL DTRMV( 'Upper', 'Transpose', 'NON-UNIT', + $ I-1, T, LDT, + $ T( 1, NB ), 1 ) +* +* b2 := b2 - V2*w +* + CALL DGEMV( 'NO TRANSPOSE', N-K-I+1, I-1, -ONE, + $ A( K+I, 1 ), + $ LDA, T( 1, NB ), 1, ONE, A( K+I, I ), 1 ) +* +* b1 := b1 - V1*w +* + CALL DTRMV( 'Lower', 'NO TRANSPOSE', + $ 'UNIT', I-1, + $ A( K+1, 1 ), LDA, T( 1, NB ), 1 ) + CALL DAXPY( I-1, -ONE, T( 1, NB ), 1, A( K+1, I ), 1 ) +* + A( K+I-1, I-1 ) = EI + END IF +* +* Generate the elementary reflector H(I) to annihilate +* A(K+I+1:N,I) +* + CALL DLARFG( N-K-I+1, A( K+I, I ), A( MIN( K+I+1, N ), I ), 1, + $ TAU( I ) ) + EI = A( K+I, I ) + A( K+I, I ) = ONE +* +* Compute Y(K+1:N,I) +* + CALL DGEMV( 'NO TRANSPOSE', N-K, N-K-I+1, + $ ONE, A( K+1, I+1 ), + $ LDA, A( K+I, I ), 1, ZERO, Y( K+1, I ), 1 ) + CALL DGEMV( 'Transpose', N-K-I+1, I-1, + $ ONE, A( K+I, 1 ), LDA, + $ A( K+I, I ), 1, ZERO, T( 1, I ), 1 ) + CALL DGEMV( 'NO TRANSPOSE', N-K, I-1, -ONE, + $ Y( K+1, 1 ), LDY, + $ T( 1, I ), 1, ONE, Y( K+1, I ), 1 ) + CALL DSCAL( N-K, TAU( I ), Y( K+1, I ), 1 ) +* +* Compute T(1:I,I) +* + CALL DSCAL( I-1, -TAU( I ), T( 1, I ), 1 ) + CALL DTRMV( 'Upper', 'No Transpose', 'NON-UNIT', + $ I-1, T, LDT, + $ T( 1, I ), 1 ) + T( I, I ) = TAU( I ) +* + 10 CONTINUE + A( K+NB, NB ) = EI +* +* Compute Y(1:K,1:NB) +* + CALL DLACPY( 'ALL', K, NB, A( 1, 2 ), LDA, Y, LDY ) + CALL DTRMM( 'RIGHT', 'Lower', 'NO TRANSPOSE', + $ 'UNIT', K, NB, + $ ONE, A( K+1, 1 ), LDA, Y, LDY ) + IF( N.GT.K+NB ) + $ CALL DGEMM( 'NO TRANSPOSE', 'NO TRANSPOSE', K, + $ NB, N-K-NB, ONE, + $ A( 1, 2+NB ), LDA, A( K+1+NB, 1 ), LDA, ONE, Y, + $ LDY ) + CALL DTRMM( 'RIGHT', 'Upper', 'NO TRANSPOSE', + $ 'NON-UNIT', K, NB, + $ ONE, T, LDT, Y, LDY ) +* + RETURN +* +* End of DLAHR2 +* + END |