version 1.2, 2012/08/22 09:48:27
|
version 1.8, 2017/06/17 11:06:37
|
Line 1
|
Line 1
|
*> \brief \b DTPRFB |
*> \brief \b DTPRFB applies a real or complex "triangular-pentagonal" blocked reflector to a real or complex matrix, which is composed of two blocks. |
* |
* |
* =========== DOCUMENTATION =========== |
* =========== DOCUMENTATION =========== |
* |
* |
* Online html documentation available at |
* Online html documentation available at |
* http://www.netlib.org/lapack/explore-html/ |
* http://www.netlib.org/lapack/explore-html/ |
* |
* |
*> \htmlonly |
*> \htmlonly |
*> Download DTPRFB + dependencies |
*> Download DTPRFB + dependencies |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfb.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtprfb.f"> |
*> [TGZ]</a> |
*> [TGZ]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfb.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtprfb.f"> |
*> [ZIP]</a> |
*> [ZIP]</a> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfb.f"> |
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtprfb.f"> |
*> [TXT]</a> |
*> [TXT]</a> |
*> \endhtmlonly |
*> \endhtmlonly |
* |
* |
* Definition: |
* Definition: |
* =========== |
* =========== |
* |
* |
* SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, |
* SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, |
* V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) |
* V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
* CHARACTER DIRECT, SIDE, STOREV, TRANS |
* CHARACTER DIRECT, SIDE, STOREV, TRANS |
* INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N |
* INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), |
* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), |
* $ V( LDV, * ), WORK( LDWORK, * ) |
* $ V( LDV, * ), WORK( LDWORK, * ) |
* .. |
* .. |
* |
* |
* |
* |
*> \par Purpose: |
*> \par Purpose: |
* ============= |
* ============= |
*> |
*> |
*> \verbatim |
*> \verbatim |
*> |
*> |
*> DTPRFB applies a real "triangular-pentagonal" block reflector H or its |
*> DTPRFB applies a real "triangular-pentagonal" block reflector H or its |
*> transpose H**T to a real matrix C, which is composed of two |
*> transpose H**T to a real matrix C, which is composed of two |
*> blocks A and B, either from the left or right. |
*> blocks A and B, either from the left or right. |
*> |
*> |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Arguments: |
* Arguments: |
Line 80
|
Line 80
|
*> \param[in] M |
*> \param[in] M |
*> \verbatim |
*> \verbatim |
*> M is INTEGER |
*> M is INTEGER |
*> The number of rows of the matrix B. |
*> The number of rows of the matrix B. |
*> M >= 0. |
*> M >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] N |
*> \param[in] N |
*> \verbatim |
*> \verbatim |
*> N is INTEGER |
*> N is INTEGER |
*> The number of columns of the matrix B. |
*> The number of columns of the matrix B. |
*> N >= 0. |
*> N >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 95
|
Line 95
|
*> \verbatim |
*> \verbatim |
*> K is INTEGER |
*> K is INTEGER |
*> The order of the matrix T, i.e. the number of elementary |
*> The order of the matrix T, i.e. the number of elementary |
*> reflectors whose product defines the block reflector. |
*> reflectors whose product defines the block reflector. |
*> K >= 0. |
*> K >= 0. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] L |
*> \param[in] L |
*> \verbatim |
*> \verbatim |
*> L is INTEGER |
*> L is INTEGER |
*> The order of the trapezoidal part of V. |
*> The order of the trapezoidal part of V. |
*> K >= L >= 0. See Further Details. |
*> K >= L >= 0. See Further Details. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 129
|
Line 129
|
*> \verbatim |
*> \verbatim |
*> T is DOUBLE PRECISION array, dimension (LDT,K) |
*> T is DOUBLE PRECISION array, dimension (LDT,K) |
*> The triangular K-by-K matrix T in the representation of the |
*> The triangular K-by-K matrix T in the representation of the |
*> block reflector. |
*> block reflector. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDT |
*> \param[in] LDT |
*> \verbatim |
*> \verbatim |
*> LDT is INTEGER |
*> LDT is INTEGER |
*> The leading dimension of the array T. |
*> The leading dimension of the array T. |
*> LDT >= K. |
*> LDT >= K. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 144
|
Line 144
|
*> A is DOUBLE PRECISION array, dimension |
*> A is DOUBLE PRECISION array, dimension |
*> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' |
*> (LDA,N) if SIDE = 'L' or (LDA,K) if SIDE = 'R' |
*> On entry, the K-by-N or M-by-K matrix A. |
*> On entry, the K-by-N or M-by-K matrix A. |
*> On exit, A is overwritten by the corresponding block of |
*> On exit, A is overwritten by the corresponding block of |
*> H*C or H**T*C or C*H or C*H**T. See Futher Details. |
*> H*C or H**T*C or C*H or C*H**T. See Further Details. |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in] LDA |
*> \param[in] LDA |
*> \verbatim |
*> \verbatim |
*> LDA is INTEGER |
*> LDA is INTEGER |
*> The leading dimension of the array A. |
*> The leading dimension of the array A. |
*> If SIDE = 'L', LDC >= max(1,K); |
*> If SIDE = 'L', LDC >= max(1,K); |
*> If SIDE = 'R', LDC >= max(1,M). |
*> If SIDE = 'R', LDC >= max(1,M). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
*> \param[in,out] B |
*> \param[in,out] B |
Line 167
|
Line 167
|
*> \param[in] LDB |
*> \param[in] LDB |
*> \verbatim |
*> \verbatim |
*> LDB is INTEGER |
*> LDB is INTEGER |
*> The leading dimension of the array B. |
*> The leading dimension of the array B. |
*> LDB >= max(1,M). |
*> LDB >= max(1,M). |
*> \endverbatim |
*> \endverbatim |
*> |
*> |
Line 182
|
Line 182
|
*> \verbatim |
*> \verbatim |
*> LDWORK is INTEGER |
*> LDWORK is INTEGER |
*> The leading dimension of the array WORK. |
*> The leading dimension of the array WORK. |
*> If SIDE = 'L', LDWORK >= K; |
*> If SIDE = 'L', LDWORK >= K; |
*> if SIDE = 'R', LDWORK >= M. |
*> if SIDE = 'R', LDWORK >= M. |
*> \endverbatim |
*> \endverbatim |
* |
* |
* Authors: |
* Authors: |
* ======== |
* ======== |
* |
* |
*> \author Univ. of Tennessee |
*> \author Univ. of Tennessee |
*> \author Univ. of California Berkeley |
*> \author Univ. of California Berkeley |
*> \author Univ. of Colorado Denver |
*> \author Univ. of Colorado Denver |
*> \author NAG Ltd. |
*> \author NAG Ltd. |
* |
* |
*> \date November 2011 |
*> \date December 2016 |
* |
* |
*> \ingroup doubleOTHERauxiliary |
*> \ingroup doubleOTHERauxiliary |
* |
* |
Line 204
|
Line 204
|
*> \verbatim |
*> \verbatim |
*> |
*> |
*> The matrix C is a composite matrix formed from blocks A and B. |
*> The matrix C is a composite matrix formed from blocks A and B. |
*> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, |
*> The block B is of size M-by-N; if SIDE = 'R', A is of size M-by-K, |
*> and if SIDE = 'L', A is of size K-by-N. |
*> and if SIDE = 'L', A is of size K-by-N. |
*> |
*> |
*> If SIDE = 'R' and DIRECT = 'F', C = [A B]. |
*> If SIDE = 'R' and DIRECT = 'F', C = [A B]. |
*> |
*> |
*> If SIDE = 'L' and DIRECT = 'F', C = [A] |
*> If SIDE = 'L' and DIRECT = 'F', C = [A] |
*> [B]. |
*> [B]. |
*> |
*> |
*> If SIDE = 'R' and DIRECT = 'B', C = [B A]. |
*> If SIDE = 'R' and DIRECT = 'B', C = [B A]. |
*> |
*> |
*> If SIDE = 'L' and DIRECT = 'B', C = [B] |
*> If SIDE = 'L' and DIRECT = 'B', C = [B] |
*> [A]. |
*> [A]. |
*> |
*> |
*> The pentagonal matrix V is composed of a rectangular block V1 and a |
*> The pentagonal matrix V is composed of a rectangular block V1 and a |
*> trapezoidal block V2. The size of the trapezoidal block is determined by |
*> trapezoidal block V2. The size of the trapezoidal block is determined by |
*> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; |
*> the parameter L, where 0<=L<=K. If L=K, the V2 block of V is triangular; |
*> if L=0, there is no trapezoidal block, thus V = V1 is rectangular. |
*> if L=0, there is no trapezoidal block, thus V = V1 is rectangular. |
*> |
*> |
Line 235
|
Line 235
|
*> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) |
*> - V2 is lower trapezoidal (last L rows of K-by-K lower triangular) |
*> |
*> |
*> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] |
*> If DIRECT = 'B' and STOREV = 'R': V = [V2 V1] |
*> |
*> |
*> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) |
*> - V2 is upper trapezoidal (last L columns of K-by-K upper triangular) |
*> |
*> |
*> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. |
*> If STOREV = 'C' and SIDE = 'L', V is M-by-K with V2 L-by-K. |
Line 248
|
Line 248
|
*> \endverbatim |
*> \endverbatim |
*> |
*> |
* ===================================================================== |
* ===================================================================== |
SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, |
SUBROUTINE DTPRFB( SIDE, TRANS, DIRECT, STOREV, M, N, K, L, |
$ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) |
$ V, LDV, T, LDT, A, LDA, B, LDB, WORK, LDWORK ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.4.0) -- |
* -- LAPACK auxiliary routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2011 |
* December 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER DIRECT, SIDE, STOREV, TRANS |
CHARACTER DIRECT, SIDE, STOREV, TRANS |
INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N |
INTEGER K, L, LDA, LDB, LDT, LDV, LDWORK, M, N |
* .. |
* .. |
* .. Array Arguments .. |
* .. Array Arguments .. |
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), |
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), T( LDT, * ), |
$ V( LDV, * ), WORK( LDWORK, * ) |
$ V( LDV, * ), WORK( LDWORK, * ) |
* .. |
* .. |
* |
* |
Line 322
|
Line 322
|
END IF |
END IF |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN |
IF( COLUMN .AND. FORWARD .AND. LEFT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 336
|
Line 336
|
* H = I - W T W**T or H**T = I - W T**T W**T |
* H = I - W T W**T or H**T = I - W T**T W**T |
* |
* |
* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) |
* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) |
* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) |
* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
MP = MIN( M-L+1, M ) |
MP = MIN( M-L+1, M ) |
KP = MIN( L+1, K ) |
KP = MIN( L+1, K ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
DO I = 1, L |
DO I = 1, L |
WORK( I, J ) = B( M-L+I, J ) |
WORK( I, J ) = B( M-L+I, J ) |
END DO |
END DO |
END DO |
END DO |
CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV, |
CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( MP, 1 ), LDV, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, |
CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V, LDV, B, LDB, |
$ ONE, WORK, LDWORK ) |
$ ONE, WORK, LDWORK ) |
CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, |
CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V( 1, KP ), LDV, |
$ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) |
$ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
DO I = 1, K |
DO I = 1, K |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, |
CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
DO I = 1, K |
DO I = 1, K |
A( I, J ) = A( I, J ) - WORK( I, J ) |
A( I, J ) = A( I, J ) - WORK( I, J ) |
Line 373
|
Line 373
|
CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, |
CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, |
$ ONE, B, LDB ) |
$ ONE, B, LDB ) |
CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, |
CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V( MP, KP ), LDV, |
$ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) |
$ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) |
CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, |
CALL DTRMM( 'L', 'U', 'N', 'N', L, N, ONE, V( MP, 1 ), LDV, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
DO J = 1, N |
DO J = 1, N |
Line 383
|
Line 383
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN |
ELSE IF( COLUMN .AND. FORWARD .AND. RIGHT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 402
|
Line 402
|
* |
* |
NP = MIN( N-L+1, N ) |
NP = MIN( N-L+1, N ) |
KP = MIN( L+1, K ) |
KP = MIN( L+1, K ) |
* |
* |
DO J = 1, L |
DO J = 1, L |
DO I = 1, M |
DO I = 1, M |
WORK( I, J ) = B( I, N-L+J ) |
WORK( I, J ) = B( I, N-L+J ) |
Line 410
|
Line 410
|
END DO |
END DO |
CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, |
CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( NP, 1 ), LDV, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, |
CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B, LDB, |
$ V, LDV, ONE, WORK, LDWORK ) |
$ V, LDV, ONE, WORK, LDWORK ) |
CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, |
CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, |
$ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) |
$ V( 1, KP ), LDV, ZERO, WORK( 1, KP ), LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
DO I = 1, M |
DO I = 1, M |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, |
CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
DO I = 1, M |
DO I = 1, M |
A( I, J ) = A( I, J ) - WORK( I, J ) |
A( I, J ) = A( I, J ) - WORK( I, J ) |
Line 443
|
Line 443
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN |
ELSE IF( COLUMN .AND. BACKWARD .AND. LEFT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 457
|
Line 457
|
* H = I - W T W**T or H**T = I - W T**T W**T |
* H = I - W T W**T or H**T = I - W T**T W**T |
* |
* |
* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) |
* A = A - T (A + V**T B) or A = A - T**T (A + V**T B) |
* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) |
* B = B - V T (A + V**T B) or B = B - V T**T (A + V**T B) |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
Line 472
|
Line 472
|
* |
* |
CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV, |
CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, KP ), LDV, |
$ WORK( KP, 1 ), LDWORK ) |
$ WORK( KP, 1 ), LDWORK ) |
CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, |
CALL DGEMM( 'T', 'N', L, N, M-L, ONE, V( MP, KP ), LDV, |
$ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) |
$ B( MP, 1 ), LDB, ONE, WORK( KP, 1 ), LDWORK ) |
CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV, |
CALL DGEMM( 'T', 'N', K-L, N, M, ONE, V, LDV, |
$ B, LDB, ZERO, WORK, LDWORK ) |
$ B, LDB, ZERO, WORK, LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
DO I = 1, K |
DO I = 1, K |
Line 483
|
Line 483
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, |
CALL DTRMM( 'L', 'L', TRANS, 'N', K, N, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
DO I = 1, K |
DO I = 1, K |
A( I, J ) = A( I, J ) - WORK( I, J ) |
A( I, J ) = A( I, J ) - WORK( I, J ) |
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, |
CALL DGEMM( 'N', 'N', M-L, N, K, -ONE, V( MP, 1 ), LDV, |
$ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) |
$ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) |
CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, |
CALL DGEMM( 'N', 'N', L, N, K-L, -ONE, V, LDV, |
$ WORK, LDWORK, ONE, B, LDB ) |
$ WORK, LDWORK, ONE, B, LDB ) |
Line 505
|
Line 505
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN |
ELSE IF( COLUMN .AND. BACKWARD .AND. RIGHT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 524
|
Line 524
|
* |
* |
NP = MIN( L+1, N ) |
NP = MIN( L+1, N ) |
KP = MIN( K-L+1, K ) |
KP = MIN( K-L+1, K ) |
* |
* |
DO J = 1, L |
DO J = 1, L |
DO I = 1, M |
DO I = 1, M |
WORK( I, K-L+J ) = B( I, J ) |
WORK( I, K-L+J ) = B( I, J ) |
Line 532
|
Line 532
|
END DO |
END DO |
CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, |
CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, KP ), LDV, |
$ WORK( 1, KP ), LDWORK ) |
$ WORK( 1, KP ), LDWORK ) |
CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, |
CALL DGEMM( 'N', 'N', M, L, N-L, ONE, B( 1, NP ), LDB, |
$ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) |
$ V( NP, KP ), LDV, ONE, WORK( 1, KP ), LDWORK ) |
CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, |
CALL DGEMM( 'N', 'N', M, K-L, N, ONE, B, LDB, |
$ V, LDV, ZERO, WORK, LDWORK ) |
$ V, LDV, ZERO, WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
DO I = 1, M |
DO I = 1, M |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
WORK( I, J ) = WORK( I, J ) + A( I, J ) |
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, |
CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
DO I = 1, M |
DO I = 1, M |
A( I, J ) = A( I, J ) - WORK( I, J ) |
A( I, J ) = A( I, J ) - WORK( I, J ) |
Line 565
|
Line 565
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN |
ELSE IF( ROW .AND. FORWARD .AND. LEFT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 578
|
Line 578
|
* H = I - W**T T W or H**T = I - W**T T**T W |
* H = I - W**T T W or H**T = I - W**T T**T W |
* |
* |
* A = A - T (A + V B) or A = A - T**T (A + V B) |
* A = A - T (A + V B) or A = A - T**T (A + V B) |
* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) |
* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
Line 589
|
Line 589
|
DO I = 1, L |
DO I = 1, L |
WORK( I, J ) = B( M-L+I, J ) |
WORK( I, J ) = B( M-L+I, J ) |
END DO |
END DO |
END DO |
END DO |
CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, |
CALL DTRMM( 'L', 'L', 'N', 'N', L, N, ONE, V( 1, MP ), LDV, |
$ WORK, LDB ) |
$ WORK, LDB ) |
CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, |
CALL DGEMM( 'N', 'N', L, N, M-L, ONE, V, LDV,B, LDB, |
$ ONE, WORK, LDWORK ) |
$ ONE, WORK, LDWORK ) |
CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, |
CALL DGEMM( 'N', 'N', K-L, N, M, ONE, V( KP, 1 ), LDV, |
$ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) |
$ B, LDB, ZERO, WORK( KP, 1 ), LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
Line 603
|
Line 603
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, |
CALL DTRMM( 'L', 'U', TRANS, 'N', K, N, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, N |
DO J = 1, N |
Line 614
|
Line 614
|
* |
* |
CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, |
CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V, LDV, WORK, LDWORK, |
$ ONE, B, LDB ) |
$ ONE, B, LDB ) |
CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, |
CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V( KP, MP ), LDV, |
$ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) |
$ WORK( KP, 1 ), LDWORK, ONE, B( MP, 1 ), LDB ) |
CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV, |
CALL DTRMM( 'L', 'L', 'T', 'N', L, N, ONE, V( 1, MP ), LDV, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
Line 625
|
Line 625
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN |
ELSE IF( ROW .AND. FORWARD .AND. RIGHT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 653
|
Line 653
|
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV, |
CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B, LDB, V, LDV, |
$ ONE, WORK, LDWORK ) |
$ ONE, WORK, LDWORK ) |
CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, |
CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, |
$ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) |
$ V( KP, 1 ), LDV, ZERO, WORK( 1, KP ), LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
Line 662
|
Line 662
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, |
CALL DTRMM( 'R', 'U', TRANS, 'N', M, K, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
Line 671
|
Line 671
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, |
CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, |
$ V, LDV, ONE, B, LDB ) |
$ V, LDV, ONE, B, LDB ) |
CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, |
CALL DGEMM( 'N', 'N', M, L, K-L, -ONE, WORK( 1, KP ), LDWORK, |
$ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) |
$ V( KP, NP ), LDV, ONE, B( 1, NP ), LDB ) |
CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, |
CALL DTRMM( 'R', 'L', 'N', 'N', M, L, ONE, V( 1, NP ), LDV, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
DO J = 1, L |
DO J = 1, L |
Line 684
|
Line 684
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN |
ELSE IF( ROW .AND. BACKWARD .AND. LEFT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 697
|
Line 697
|
* H = I - W**T T W or H**T = I - W**T T**T W |
* H = I - W**T T W or H**T = I - W**T T**T W |
* |
* |
* A = A - T (A + V B) or A = A - T**T (A + V B) |
* A = A - T (A + V B) or A = A - T**T (A + V B) |
* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) |
* B = B - V**T T (A + V B) or B = B - V**T T**T (A + V B) |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
Line 733
|
Line 733
|
* |
* |
CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, |
CALL DGEMM( 'T', 'N', M-L, N, K, -ONE, V( 1, MP ), LDV, |
$ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) |
$ WORK, LDWORK, ONE, B( MP, 1 ), LDB ) |
CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, |
CALL DGEMM( 'T', 'N', L, N, K-L, -ONE, V, LDV, |
$ WORK, LDWORK, ONE, B, LDB ) |
$ WORK, LDWORK, ONE, B, LDB ) |
CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV, |
CALL DTRMM( 'L', 'U', 'T', 'N', L, N, ONE, V( KP, 1 ), LDV, |
$ WORK( KP, 1 ), LDWORK ) |
$ WORK( KP, 1 ), LDWORK ) |
DO J = 1, N |
DO J = 1, N |
DO I = 1, L |
DO I = 1, L |
B( I, J ) = B( I, J ) - WORK( K-L+I, J ) |
B( I, J ) = B( I, J ) - WORK( K-L+I, J ) |
Line 744
|
Line 744
|
END DO |
END DO |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
* |
* |
ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN |
ELSE IF( ROW .AND. BACKWARD .AND. RIGHT ) THEN |
* |
* |
* --------------------------------------------------------------------------- |
* --------------------------------------------------------------------------- |
Line 773
|
Line 773
|
CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB, |
CALL DGEMM( 'N', 'T', M, L, N-L, ONE, B( 1, NP ), LDB, |
$ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) |
$ V( KP, NP ), LDV, ONE, WORK( 1, KP ), LDWORK ) |
CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV, |
CALL DGEMM( 'N', 'T', M, K-L, N, ONE, B, LDB, V, LDV, |
$ ZERO, WORK, LDWORK ) |
$ ZERO, WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
DO I = 1, M |
DO I = 1, M |
Line 781
|
Line 781
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, |
CALL DTRMM( 'R', 'L', TRANS, 'N', M, K, ONE, T, LDT, |
$ WORK, LDWORK ) |
$ WORK, LDWORK ) |
* |
* |
DO J = 1, K |
DO J = 1, K |
Line 790
|
Line 790
|
END DO |
END DO |
END DO |
END DO |
* |
* |
CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, |
CALL DGEMM( 'N', 'N', M, N-L, K, -ONE, WORK, LDWORK, |
$ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) |
$ V( 1, NP ), LDV, ONE, B( 1, NP ), LDB ) |
CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, |
CALL DGEMM( 'N', 'N', M, L, K-L , -ONE, WORK, LDWORK, |
$ V, LDV, ONE, B, LDB ) |
$ V, LDV, ONE, B, LDB ) |
CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, |
CALL DTRMM( 'R', 'U', 'N', 'N', M, L, ONE, V( KP, 1 ), LDV, |
$ WORK( 1, KP ), LDWORK ) |
$ WORK( 1, KP ), LDWORK ) |