diff options
author | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
---|---|---|
committer | jason <jason@8a072113-8704-0410-8d35-dd094bca7971> | 2008-10-28 01:38:50 +0000 |
commit | baba851215b44ac3b60b9248eb02bcce7eb76247 (patch) | |
tree | 8c0f5c006875532a30d4409f5e94b0f310ff00a7 /SRC/sopmtr.f | |
download | lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.gz lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.tar.bz2 lapack-baba851215b44ac3b60b9248eb02bcce7eb76247.zip |
Move LAPACK trunk into position.
Diffstat (limited to 'SRC/sopmtr.f')
-rw-r--r-- | SRC/sopmtr.f | 257 |
1 files changed, 257 insertions, 0 deletions
diff --git a/SRC/sopmtr.f b/SRC/sopmtr.f new file mode 100644 index 00000000..f6779dc3 --- /dev/null +++ b/SRC/sopmtr.f @@ -0,0 +1,257 @@ + SUBROUTINE SOPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, + $ INFO ) +* +* -- LAPACK routine (version 3.1) -- +* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. +* November 2006 +* +* .. Scalar Arguments .. + CHARACTER SIDE, TRANS, UPLO + INTEGER INFO, LDC, M, N +* .. +* .. Array Arguments .. + REAL AP( * ), C( LDC, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* SOPMTR overwrites the general real M-by-N matrix C with +* +* SIDE = 'L' SIDE = 'R' +* TRANS = 'N': Q * C C * Q +* TRANS = 'T': Q**T * C C * Q**T +* +* where Q is a real orthogonal matrix of order nq, with nq = m if +* SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of +* nq-1 elementary reflectors, as returned by SSPTRD using packed +* storage: +* +* if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1); +* +* if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1). +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q**T from the Left; +* = 'R': apply Q or Q**T from the Right. +* +* UPLO (input) CHARACTER*1 +* = 'U': Upper triangular packed storage used in previous +* call to SSPTRD; +* = 'L': Lower triangular packed storage used in previous +* call to SSPTRD. +* +* TRANS (input) CHARACTER*1 +* = 'N': No transpose, apply Q; +* = 'T': Transpose, apply Q**T. +* +* 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. +* +* AP (input) REAL array, dimension +* (M*(M+1)/2) if SIDE = 'L' +* (N*(N+1)/2) if SIDE = 'R' +* The vectors which define the elementary reflectors, as +* returned by SSPTRD. AP is modified by the routine but +* restored on exit. +* +* TAU (input) REAL array, dimension (M-1) if SIDE = 'L' +* or (N-1) if SIDE = 'R' +* TAU(i) must contain the scalar factor of the elementary +* reflector H(i), as returned by SSPTRD. +* +* C (input/output) REAL array, dimension (LDC,N) +* On entry, the M-by-N matrix C. +* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. +* +* LDC (input) INTEGER +* The leading dimension of the array C. LDC >= max(1,M). +* +* WORK (workspace) REAL array, dimension +* (N) if SIDE = 'L' +* (M) if SIDE = 'R' +* +* INFO (output) INTEGER +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value +* +* ===================================================================== +* +* .. Parameters .. + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) +* .. +* .. Local Scalars .. + LOGICAL FORWRD, LEFT, NOTRAN, UPPER + INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ + REAL AII +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL SLARF, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + NOTRAN = LSAME( TRANS, 'N' ) + UPPER = LSAME( UPLO, 'U' ) +* +* NQ is the order of Q +* + IF( LEFT ) THEN + NQ = M + ELSE + NQ = N + END IF + IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN + INFO = -1 + ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN + INFO = -2 + ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN + INFO = -3 + ELSE IF( M.LT.0 ) THEN + INFO = -4 + ELSE IF( N.LT.0 ) THEN + INFO = -5 + ELSE IF( LDC.LT.MAX( 1, M ) ) THEN + INFO = -9 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'SOPMTR', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 ) + $ RETURN +* + IF( UPPER ) THEN +* +* Q was determined by a call to SSPTRD with UPLO = 'U' +* + FORWRD = ( LEFT .AND. NOTRAN ) .OR. + $ ( .NOT.LEFT .AND. .NOT.NOTRAN ) +* + IF( FORWRD ) THEN + I1 = 1 + I2 = NQ - 1 + I3 = 1 + II = 2 + ELSE + I1 = NQ - 1 + I2 = 1 + I3 = -1 + II = NQ*( NQ+1 ) / 2 - 1 + END IF +* + IF( LEFT ) THEN + NI = N + ELSE + MI = M + END IF +* + DO 10 I = I1, I2, I3 + IF( LEFT ) THEN +* +* H(i) is applied to C(1:i,1:n) +* + MI = I + ELSE +* +* H(i) is applied to C(1:m,1:i) +* + NI = I + END IF +* +* Apply H(i) +* + AII = AP( II ) + AP( II ) = ONE + CALL SLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, LDC, + $ WORK ) + AP( II ) = AII +* + IF( FORWRD ) THEN + II = II + I + 2 + ELSE + II = II - I - 1 + END IF + 10 CONTINUE + ELSE +* +* Q was determined by a call to SSPTRD with UPLO = 'L'. +* + FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR. + $ ( .NOT.LEFT .AND. NOTRAN ) +* + IF( FORWRD ) THEN + I1 = 1 + I2 = NQ - 1 + I3 = 1 + II = 2 + ELSE + I1 = NQ - 1 + I2 = 1 + I3 = -1 + II = NQ*( NQ+1 ) / 2 - 1 + END IF +* + IF( LEFT ) THEN + NI = N + JC = 1 + ELSE + MI = M + IC = 1 + END IF +* + DO 20 I = I1, I2, I3 + AII = AP( II ) + AP( II ) = ONE + IF( LEFT ) THEN +* +* H(i) is applied to C(i+1:m,1:n) +* + MI = M - I + IC = I + 1 + ELSE +* +* H(i) is applied to C(1:m,i+1:n) +* + NI = N - I + JC = I + 1 + END IF +* +* Apply H(i) +* + CALL SLARF( SIDE, MI, NI, AP( II ), 1, TAU( I ), + $ C( IC, JC ), LDC, WORK ) + AP( II ) = AII +* + IF( FORWRD ) THEN + II = II + NQ - I + 1 + ELSE + II = II - NQ + I - 2 + END IF + 20 CONTINUE + END IF + RETURN +* +* End of SOPMTR +* + END |