diff options
Diffstat (limited to 'src/fortran/lapack/dorm2r.f')
| -rw-r--r-- | src/fortran/lapack/dorm2r.f | 197 |
1 files changed, 197 insertions, 0 deletions
diff --git a/src/fortran/lapack/dorm2r.f b/src/fortran/lapack/dorm2r.f new file mode 100644 index 0000000..79c9ef3 --- /dev/null +++ b/src/fortran/lapack/dorm2r.f @@ -0,0 +1,197 @@ + SUBROUTINE DORM2R( SIDE, TRANS, M, N, K, A, LDA, 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 + INTEGER INFO, K, LDA, LDC, M, N +* .. +* .. Array Arguments .. + DOUBLE PRECISION A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * ) +* .. +* +* Purpose +* ======= +* +* DORM2R overwrites the general real m by n matrix C with +* +* Q * C if SIDE = 'L' and TRANS = 'N', or +* +* Q'* C if SIDE = 'L' and TRANS = 'T', or +* +* C * Q if SIDE = 'R' and TRANS = 'N', or +* +* C * Q' if SIDE = 'R' and TRANS = 'T', +* +* where Q is a real orthogonal matrix defined as the product of k +* elementary reflectors +* +* Q = H(1) H(2) . . . H(k) +* +* as returned by DGEQRF. Q is of order m if SIDE = 'L' and of order n +* if SIDE = 'R'. +* +* Arguments +* ========= +* +* SIDE (input) CHARACTER*1 +* = 'L': apply Q or Q' from the Left +* = 'R': apply Q or Q' from the Right +* +* TRANS (input) CHARACTER*1 +* = 'N': apply Q (No transpose) +* = 'T': apply Q' (Transpose) +* +* 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. +* +* K (input) INTEGER +* The number of elementary reflectors whose product defines +* the matrix Q. +* If SIDE = 'L', M >= K >= 0; +* if SIDE = 'R', N >= K >= 0. +* +* A (input) DOUBLE PRECISION array, dimension (LDA,K) +* The i-th column must contain the vector which defines the +* elementary reflector H(i), for i = 1,2,...,k, as returned by +* DGEQRF in the first k columns of its array argument A. +* A is modified by the routine but restored on exit. +* +* LDA (input) INTEGER +* The leading dimension of the array A. +* If SIDE = 'L', LDA >= max(1,M); +* if SIDE = 'R', LDA >= max(1,N). +* +* TAU (input) DOUBLE PRECISION array, dimension (K) +* TAU(i) must contain the scalar factor of the elementary +* reflector H(i), as returned by DGEQRF. +* +* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) +* On entry, the m by n matrix C. +* On exit, C is overwritten by Q*C or Q'*C or C*Q' or C*Q. +* +* LDC (input) INTEGER +* The leading dimension of the array C. LDC >= max(1,M). +* +* WORK (workspace) DOUBLE PRECISION 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 .. + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) +* .. +* .. Local Scalars .. + LOGICAL LEFT, NOTRAN + INTEGER I, I1, I2, I3, IC, JC, MI, NI, NQ + DOUBLE PRECISION AII +* .. +* .. External Functions .. + LOGICAL LSAME + EXTERNAL LSAME +* .. +* .. External Subroutines .. + EXTERNAL DLARF, XERBLA +* .. +* .. Intrinsic Functions .. + INTRINSIC MAX +* .. +* .. Executable Statements .. +* +* Test the input arguments +* + INFO = 0 + LEFT = LSAME( SIDE, 'L' ) + NOTRAN = LSAME( TRANS, 'N' ) +* +* 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.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN + INFO = -2 + ELSE IF( M.LT.0 ) THEN + INFO = -3 + ELSE IF( N.LT.0 ) THEN + INFO = -4 + ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN + INFO = -5 + ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN + INFO = -7 + ELSE IF( LDC.LT.MAX( 1, M ) ) THEN + INFO = -10 + END IF + IF( INFO.NE.0 ) THEN + CALL XERBLA( 'DORM2R', -INFO ) + RETURN + END IF +* +* Quick return if possible +* + IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) + $ RETURN +* + IF( ( LEFT .AND. .NOT.NOTRAN ) .OR. ( .NOT.LEFT .AND. NOTRAN ) ) + $ THEN + I1 = 1 + I2 = K + I3 = 1 + ELSE + I1 = K + I2 = 1 + I3 = -1 + END IF +* + IF( LEFT ) THEN + NI = N + JC = 1 + ELSE + MI = M + IC = 1 + END IF +* + DO 10 I = I1, I2, I3 + IF( LEFT ) THEN +* +* H(i) is applied to C(i:m,1:n) +* + MI = M - I + 1 + IC = I + ELSE +* +* H(i) is applied to C(1:m,i:n) +* + NI = N - I + 1 + JC = I + END IF +* +* Apply H(i) +* + AII = A( I, I ) + A( I, I ) = ONE + CALL DLARF( SIDE, MI, NI, A( I, I ), 1, TAU( I ), C( IC, JC ), + $ LDC, WORK ) + A( I, I ) = AII + 10 CONTINUE + RETURN +* +* End of DORM2R +* + END |