SUBROUTINE MB04IY( SIDE, TRANS, N, M, K, P, A, LDA, TAU, C, LDC, $ DWORK, LDWORK, INFO ) C C RELEASE 4.0, WGS COPYRIGHT 1999. C C PURPOSE C C To overwrite the real n-by-m matrix C with Q' * C, Q * C, C C * Q', or C * Q, according to the following table C C SIDE = 'L' SIDE = 'R' C TRANS = 'N': Q * C C * Q C TRANS = 'T': Q'* C C * Q' C C where Q is a real orthogonal matrix defined as the product of C k elementary reflectors C C Q = H(1) H(2) . . . H(k) C C as returned by SLICOT Library routine MB04ID. Q is of order n C if SIDE = 'L' and of order m if SIDE = 'R'. C C ARGUMENTS C C Mode Parameters C C SIDE CHARACTER*1 C Specify if Q or Q' is applied from the left or right, C as follows: C = 'L': apply Q or Q' from the left; C = 'R': apply Q or Q' from the right. C C TRANS CHARACTER*1 C Specify if Q or Q' is to be applied, as follows: C = 'N': apply Q (No transpose); C = 'T': apply Q' (Transpose). C C Input/Output Parameters C C N (input) INTEGER C The number of rows of the matrix C. N >= 0. C C M (input) INTEGER C The number of columns of the matrix C. M >= 0. C C K (input) INTEGER C The number of elementary reflectors whose product defines C the matrix Q. C N >= K >= 0, if SIDE = 'L'; C M >= K >= 0, if SIDE = 'R'. C C P (input) INTEGER C The order of the zero triagle (or the number of rows of C the zero trapezoid) in the matrix triangularized by SLICOT C Library routine MB04ID. P >= 0. C C A (input) DOUBLE PRECISION array, dimension (LDA,K) C On input, the elements in the rows i+1:min(n,n-p-1+i) of C the i-th column, and TAU(i), represent the orthogonal C reflector H(i), so that matrix Q is the product of C elementary reflectors: Q = H(1) H(2) . . . H(k). C A is modified by the routine but restored on exit. C C LDA INTEGER C The leading dimension of the array A. C LDA >= max(1,N), if SIDE = 'L'; C LDA >= max(1,M), if SIDE = 'R'. C C TAU (input) DOUBLE PRECISION array, dimension (K) C The scalar factors of the elementary reflectors. C C C (input/output) DOUBLE PRECISION array, dimension (LDC,M) C On entry, the leading N-by-M part of this array must C contain the matrix C. C On exit, the leading N-by-M part of this array contains C the updated matrix C. C C LDC INTEGER C The leading dimension of the array C. LDC >= max(1,N). C C Workspace C C DWORK DOUBLE PRECISION array, dimension (LDWORK) C On exit, if INFO = 0, DWORK(1) returns the optimal value C of LDWORK. C C LDWORK INTEGER C The length of the array DWORK. C LDWORK >= MAX(1,M), if SIDE = 'L'; C LDWORK >= MAX(1,N), if SIDE = 'R'. C For optimum performance LDWORK >= M*NB if SIDE = 'L', C or LDWORK >= N*NB if SIDE = 'R', where NB is the optimal C block size. C C Error Indicator C C INFO INTEGER C = 0: successful exit; C < 0: if INFO = -i, the i-th argument had an illegal C value. C C METHOD C C If SIDE = 'L', each elementary reflector H(i) modifies C n-p elements of each column of C, for i = 1:p+1, and C n-i+1 elements, for i = p+2:k. C If SIDE = 'R', each elementary reflector H(i) modifies C m-p elements of each row of C, for i = 1:p+1, and C m-i+1 elements, for i = p+2:k. C C NUMERICAL ASPECTS C C The implemented method is numerically stable. C C CONTRIBUTOR C C V. Sima, Research Institute for Informatics, Bucharest, Aug. 1999. C C REVISIONS C C - C C KEYWORDS C C Matrix operations, QR decomposition. C C ****************************************************************** C C .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) C .. Scalar Arguments .. INTEGER INFO, K, LDA, LDC, LDWORK, M, N, P CHARACTER SIDE, TRANS C .. Array Arguments .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ), DWORK( * ), TAU( * ) C .. Local Scalars .. LOGICAL LEFT, TRAN INTEGER I DOUBLE PRECISION AII, WRKOPT C .. External Functions .. LOGICAL LSAME EXTERNAL LSAME C .. External Subroutines .. EXTERNAL DLARF, DORMQR, XERBLA C .. Intrinsic Functions .. INTRINSIC DBLE, MAX, MIN C .. Executable Statements .. C C Check the scalar input arguments. C INFO = 0 LEFT = LSAME( SIDE, 'L' ) TRAN = LSAME( TRANS, 'T' ) C IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.TRAN .AND. .NOT.LSAME( TRANS, 'N' ) ) THEN INFO = -2 ELSE IF( N.LT.0 ) THEN INFO = -3 ELSE IF( M.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 .OR. ( LEFT .AND. K.GT.N ) .OR. $ ( .NOT.LEFT .AND. K.GT.M ) ) THEN INFO = -5 ELSE IF( P.LT.0 ) THEN INFO = -6 ELSE IF( ( LEFT .AND. LDA.LT.MAX( 1, N ) ) .OR. $ ( .NOT.LEFT .AND. LDA.LT.MAX( 1, M ) ) ) THEN INFO = -8 ELSE IF( LDC.LT.MAX( 1, N ) ) THEN INFO = -11 ELSE IF( ( LEFT .AND. LDWORK.LT.MAX( 1, M ) ) .OR. $ ( .NOT.LEFT .AND. LDWORK.LT.MAX( 1, N ) ) ) THEN INFO = -13 END IF C IF( INFO.NE.0 ) THEN CALL XERBLA( 'MB04IY', -INFO ) RETURN END IF C C Quick return if possible. C IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 .OR. ( LEFT .AND. N.LT.P ) $ .OR. ( .NOT.LEFT .AND. M.LT.P ) ) THEN DWORK(1) = ONE RETURN END IF C C (Note: Comments in the code beginning "Workspace:" describe the C minimal amount of real workspace needed at that point in the C code, as well as the preferred amount for good performance. C NB refers to the optimal block size for the immediately C following subroutine, as returned by ILAENV.) C IF( LEFT ) THEN WRKOPT = DBLE( M ) IF( TRAN ) THEN C DO 10 I = 1, MIN( K, P ) C C Apply H(i) to C(i:i+n-p-1,1:m), from the left. C Workspace: need M. C AII = A( I, I ) A( I, I ) = ONE CALL DLARF( SIDE, N-P, M, A( I, I ), 1, TAU( I ), $ C( I, 1 ), LDC, DWORK ) A( I, I ) = AII 10 CONTINUE C IF ( P.LE.MIN( N, K ) ) THEN C C Apply H(i) to C, i = p+1:k, from the left. C Workspace: need M; prefer M*NB. C CALL DORMQR( SIDE, TRANS, N-P, M, K-P, A( P+1, P+1 ), $ LDA, TAU( P+1 ), C( P+1, 1 ), LDC, DWORK, $ LDWORK, I ) WRKOPT = MAX( WRKOPT, DWORK( 1 ) ) END IF C ELSE C IF ( P.LE.MIN( N, K ) ) THEN C C Apply H(i) to C, i = k:p+1:-1, from the left. C Workspace: need M; prefer M*NB. C CALL DORMQR( SIDE, TRANS, N-P, M, K-P, A( P+1, P+1 ), $ LDA, TAU( P+1 ), C( P+1, 1 ), LDC, DWORK, $ LDWORK, I ) WRKOPT = MAX( WRKOPT, DWORK( 1 ) ) END IF C DO 20 I = MIN( K, P ), 1, -1 C C Apply H(i) to C(i:i+n-p-1,1:m), from the left. C Workspace: need M. C AII = A( I, I ) A( I, I ) = ONE CALL DLARF( SIDE, N-P, M, A( I, I ), 1, TAU( I ), $ C( I, 1 ), LDC, DWORK ) A( I, I ) = AII 20 CONTINUE END IF C ELSE C WRKOPT = DBLE( N ) IF( TRAN ) THEN C IF ( P.LE.MIN( M, K ) ) THEN C C Apply H(i) to C, i = k:p+1:-1, from the right. C Workspace: need N; prefer N*NB. C CALL DORMQR( SIDE, TRANS, N, M-P, K-P, A( P+1, P+1 ), $ LDA, TAU( P+1 ), C( 1, P+1 ), LDC, DWORK, $ LDWORK, I ) WRKOPT = MAX( WRKOPT, DWORK( 1 ) ) END IF C DO 30 I = MIN( K, P ), 1, -1 C C Apply H(i) to C(1:n,i:i+m-p-1), from the right. C Workspace: need N. C AII = A( I, I ) A( I, I ) = ONE CALL DLARF( SIDE, N, M-P, A( I, I ), 1, TAU( I ), $ C( 1, I ), LDC, DWORK ) A( I, I ) = AII 30 CONTINUE C ELSE C DO 40 I = 1, MIN( K, P ) C C Apply H(i) to C(1:n,i:i+m-p-1), from the right. C Workspace: need N. C AII = A( I, I ) A( I, I ) = ONE CALL DLARF( SIDE, N, M-P, A( I, I ), 1, TAU( I ), $ C( 1, I ), LDC, DWORK ) A( I, I ) = AII 40 CONTINUE C IF ( P.LE.MIN( M, K ) ) THEN C C Apply H(i) to C, i = p+1:k, from the right. C Workspace: need N; prefer N*NB. C CALL DORMQR( SIDE, TRANS, N, M-P, K-P, A( P+1, P+1 ), $ LDA, TAU( P+1 ), C( 1, P+1 ), LDC, DWORK, $ LDWORK, I ) WRKOPT = MAX( WRKOPT, DWORK( 1 ) ) END IF C END IF END IF C DWORK( 1 ) = WRKOPT RETURN C C *** Last line of MB04IY *** END