--- a +++ b/src/MB04OD.f @@ -0,0 +1,257 @@ + SUBROUTINE MB04OD( UPLO, N, M, P, R, LDR, A, LDA, B, LDB, C, LDC, + $ TAU, DWORK ) +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To calculate a QR factorization of the first block column and +C apply the orthogonal transformations (from the left) also to the +C second block column of a structured matrix, as follows +C _ _ +C [ R B ] [ R B ] +C Q' * [ ] = [ _ ] +C [ A C ] [ 0 C ] +C _ +C where R and R are upper triangular. The matrix A can be full or +C upper trapezoidal/triangular. The problem structure is exploited. +C +C ARGUMENTS +C +C Mode Parameters +C +C UPLO CHARACTER*1 +C Indicates if the matrix A is or not triangular as follows: +C = 'U': Matrix A is upper trapezoidal/triangular; +C = 'F': Matrix A is full. +C +C Input/Output Parameters +C +C N (input) INTEGER _ +C The order of the matrices R and R. N >= 0. +C +C M (input) INTEGER +C The number of columns of the matrices B and C. M >= 0. +C +C P (input) INTEGER +C The number of rows of the matrices A and C. P >= 0. +C +C R (input/output) DOUBLE PRECISION array, dimension (LDR,N) +C On entry, the leading N-by-N upper triangular part of this +C array must contain the upper triangular matrix R. +C On exit, the leading N-by-N upper triangular part of this +C _ +C array contains the upper triangular matrix R. +C The strict lower triangular part of this array is not +C referenced. +C +C LDR INTEGER +C The leading dimension of array R. LDR >= MAX(1,N). +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,N) +C On entry, if UPLO = 'F', the leading P-by-N part of this +C array must contain the matrix A. If UPLO = 'U', the +C leading MIN(P,N)-by-N part of this array must contain the +C upper trapezoidal (upper triangular if P >= N) matrix A, +C and the elements below the diagonal are not referenced. +C On exit, the leading P-by-N part (upper trapezoidal or +C triangular, if UPLO = 'U') of this array contains the +C trailing components (the vectors v, see Method) of the +C elementary reflectors used in the factorization. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= MAX(1,P). +C +C B (input/output) DOUBLE PRECISION array, dimension (LDB,M) +C On entry, the leading N-by-M part of this array must +C contain the matrix B. +C On exit, the leading N-by-M part of this array contains +C _ +C the computed matrix B. +C +C LDB INTEGER +C The leading dimension of array B. LDB >= MAX(1,N). +C +C C (input/output) DOUBLE PRECISION array, dimension (LDC,M) +C On entry, the leading P-by-M part of this array must +C contain the matrix C. +C On exit, the leading P-by-M part of this array contains +C _ +C the computed matrix C. +C +C LDC INTEGER +C The leading dimension of array C. LDC >= MAX(1,P). +C +C TAU (output) DOUBLE PRECISION array, dimension (N) +C The scalar factors of the elementary reflectors used. +C +C Workspace +C +C DWORK DOUBLE PRECISION array, dimension (MAX(N-1,M)) +C +C METHOD +C +C The routine uses N Householder transformations exploiting the zero +C pattern of the block matrix. A Householder matrix has the form +C +C ( 1 ) +C H = I - tau *u *u', u = ( v ), +C i i i i i ( i) +C +C where v is a P-vector, if UPLO = 'F', or a min(i,P)-vector, if +C i +C UPLO = 'U'. The components of v are stored in the i-th column +C i +C of A, and tau is stored in TAU(i). +C i +C In-line code for applying Householder transformations is used +C whenever possible (see MB04OY routine). +C +C NUMERICAL ASPECTS +C +C The algorithm is backward stable. +C +C CONTRIBUTORS +C +C V. Sima, Katholieke Univ. Leuven, Belgium, Feb. 1997. +C +C REVISIONS +C +C Dec. 1997. +C +C KEYWORDS +C +C Elementary reflector, QR factorization, orthogonal transformation. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) +C .. Scalar Arguments .. + CHARACTER UPLO + INTEGER LDA, LDB, LDC, LDR, M, N, P +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), + $ R(LDR,*), TAU(*) +C .. Local Scalars .. + LOGICAL LUPLO + INTEGER I, IM +C .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +C .. External Subroutines .. + EXTERNAL DLARFG, MB04OY +C .. Intrinsic Functions .. + INTRINSIC MIN +C .. Executable Statements .. +C +C For efficiency reasons, the parameters are not checked. +C + IF( MIN( N, P ).EQ.0 ) + $ RETURN +C + LUPLO = LSAME( UPLO, 'U' ) + IF ( LUPLO ) THEN +C + DO 10 I = 1, N +C +C Annihilate the I-th column of A and apply the +C transformations to the entire block matrix, exploiting +C its structure. +C + IM = MIN( I, P ) + CALL DLARFG( IM+1, R(I,I), A(1,I), 1, TAU(I) ) +C +C Compute +C [ R(I,I+1:N) ] +C w := [ 1 v' ] * [ ], +C [ A(1:IM,I+1:N) ] +C +C [ R(I,I+1:N) ] [ R(I,I+1:N) ] [ 1 ] +C [ ] := [ ] - tau * [ ] * w . +C [ A(1:IM,I+1:N) ] [ A(1:IM,I+1:N) ] [ v ] +C + IF ( N-I.GT.0 ) + $ CALL MB04OY( IM, N-I, A(1,I), TAU(I), R(I,I+1), LDR, + $ A(1,I+1), LDA, DWORK ) +C +C Compute +C [ B(I,:) ] +C w := [ 1 v' ] * [ ], +C [ C(1:IM,:) ] +C +C [ B(I,:) ] [ B(I,:) ] [ 1 ] +C [ ] := [ ] - tau * [ ] * w. +C [ C(1:IM,:) ] [ C(1:IM,:) ] [ v ] +C +C + IF ( M.GT.0 ) + $ CALL MB04OY( IM, M, A(1,I), TAU(I), B(I,1), LDB, C, LDC, + $ DWORK ) + 10 CONTINUE +C + ELSE +C + DO 20 I = 1, N - 1 +C +C Annihilate the I-th column of A and apply the +C transformations to the first block column, exploiting its +C structure. +C + CALL DLARFG( P+1, R(I,I), A(1,I), 1, TAU(I) ) +C +C Compute +C [ R(I,I+1:N) ] +C w := [ 1 v' ] * [ ], +C [ A(:,I+1:N) ] +C +C [ R(I,I+1:N) ] [ R(I,I+1:N) ] [ 1 ] +C [ ] := [ ] - tau * [ ] * w . +C [ A(:,I+1:N) ] [ A(:,I+1:N) ] [ v ] +C + CALL MB04OY( P, N-I, A(1,I), TAU(I), R(I,I+1), LDR, + $ A(1,I+1), LDA, DWORK ) + 20 CONTINUE +C + CALL DLARFG( P+1, R(N,N), A(1,N), 1, TAU(N) ) + IF ( M.GT.0 ) THEN +C +C Apply the transformations to the second block column. +C + DO 30 I = 1, N +C +C Compute +C [ B(I,:) ] +C w := [ 1 v' ] * [ ], +C [ C ] +C +C [ B(I,:) ] [ B(I,:) ] [ 1 ] +C [ ] := [ ] - tau * [ ] * w. +C [ C ] [ C ] [ v ] +C + CALL MB04OY( P, M, A(1,I), TAU(I), B(I,1), LDB, C, LDC, + $ DWORK ) + 30 CONTINUE +C + END IF + END IF + RETURN +C *** Last line of MB04OD *** + END