diff --git a/SRC/clarft.f b/SRC/clarft.f index dbd2e49a2..78c754db8 100644 --- a/SRC/clarft.f +++ b/SRC/clarft.f @@ -48,6 +48,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**H * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -60,6 +63,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -67,8 +72,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -139,7 +146,7 @@ *> the H(i) is best illustrated by the following example with n = 5 and *> k = 3. The elements equal to 1 are not stored. *> -*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and (STOREV = 'R' or STOREV = 'T'): *> *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) *> ( v1 1 ) ( 1 v2 v2 v2 ) @@ -147,13 +154,36 @@ *> ( v1 v2 v3 ) *> ( v1 v2 v3 ) *> -*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> DIRECT = 'B' and STOREV = 'C': (DIRECT = 'B' or DIRECT = 'T') and STOREV = 'R': *> *> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) *> ( v1 v2 v3 ) ( v2 v2 v2 1 ) *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -167,28 +197,30 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * .. Scalar Arguments * - CHARACTER DIRECT, STOREV - INTEGER K, LDT, LDV, N + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. * - COMPLEX T( LDT, * ), TAU( * ), V( LDV, * ) + COMPLEX T( LDT, * ), TAU( * ), V( LDV, * ) * .. * * .. Parameters .. * - COMPLEX ONE, NEG_ONE, ZERO - PARAMETER(ONE=(1.0E+0,0.0E+0), ZERO = (0.0E+0,0.0E+0), - $ NEG_ONE=(-1.0E+0,0.0E+0)) + COMPLEX ONE, NEG_ONE + PARAMETER(ONE=(1.0E+0,0.0E+0), + $ NEG_ONE=(-1.0E+0,0.0E+0)) * * .. Local Scalars .. * - INTEGER I,J,L,NX - LOGICAL QR,LQ,QL,DIRF,COLV + INTEGER I,J,L,NX,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET * * .. External Subroutines .. * - EXTERNAL CTRMM, CGEMM, CLACPY, CLARFT_LVL2 + EXTERNAL CLARFT_LVL2, CTRMM, CGEMM, CLACPY, XERBLA * * .. External Functions.. * @@ -199,7 +231,7 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * .. Intrinsic Functions.. * INTRINSIC CONJG -* +* * The general scheme used is inspired by the approach inside DGEQRT3 * which was (at the time of writing this code): * Based on the algorithm of Elmroth and Gustavson, @@ -207,57 +239,115 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * .. * .. Executable Statements .. * -* Quick return if possible +* Convert our character flags to logical flags for later * - IF(N.EQ.0.OR.K.EQ.0) THEN - RETURN - END IF + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') * -* Base case +* Error handling for our character flags * - IF(N.EQ.1.OR.K.EQ.1) THEN - T(1,1) = TAU(1) + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('CLARFT', INFO) RETURN END IF * -* Determine when to cross over into the level 2 based implementation +* Quick return if possible * - NX = ILAENV(3, "CLARFT", DIRECT // STOREV, N, K, -1, -1) - IF(K.LT.NX) THEN - CALL CLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) + IF(N.EQ.0.OR.K.EQ.0) THEN RETURN END IF * -* Beginning of executable statements -* - L = K / 2 +* Now we determine what factorization our flags are associated with * -* Determine what kind of Q we need to compute -* We assume that if the user doesn't provide 'F' for DIRECT, -* then they meant to provide 'B' and if they don't provide -* 'C' for STOREV, then they meant to provide 'R' +* QR happens when we have forward direction in column storage * - DIRF = LSAME(DIRECT,'F') - COLV = LSAME(STOREV,'C') + QR = DIRF.AND.STOREC * -* QR happens when we have forward direction in column storage +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute * - QR = DIRF.AND.COLV + LQT = DIRF.AND.STORET * -* LQ happens when we have forward direction in row storage +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute * - LQ = DIRF.AND.(.NOT.COLV) + LQ = DIRF.AND.STORER * * QL happens when we have backward direction in column storage * - QL = (.NOT.DIRF).AND.COLV + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute * -* The last case is RQ. Due to how we structured this, if the -* above 3 are false, then RQ must be true, so we never store -* this -* RQ happens when we have backward direction in row storage -* RQ = (.NOT.DIRF).AND.(.NOT.COLV) + RQT = DIRT.AND.STORER * +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER +* +* Base case +* + IF(N.EQ.1.OR.K.EQ.1) THEN + IF( LQT.OR.RQT ) THEN + T(1,1) = CONJG(TAU(1)) + ELSE + T(1,1) = TAU(1) + END IF + RETURN + END IF +* +* Determine crossover point from level 2 to level 3 BLAS implementation +* + NX = ILAENV(3, "CLARFT", DIRECT // STOREV, N, K, -1, -1) + IF(K.LT.NX) THEN +* +* Finish this component with a level 2 BLAS implementation +* + CALL CLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + L = K / 2 IF(QR) THEN * * Break V apart into 6 components @@ -271,17 +361,17 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\C^{l,l} unit lower triangular * V_{2,1}\in\C^{k-l,l} rectangular * V_{3,1}\in\C^{n-k,l} rectangular -* +* * V_{2,2}\in\C^{k-l,k-l} unit lower triangular * V_{3,2}\in\C^{n-k,k-l} rectangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -292,17 +382,17 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_1*T_{1,1}*V_1')*(I - V_2*T_{2,2}*V_2') * = I - V_1*T_{1,1}*V_1' - V_2*T_{2,2}*V_2' + V_1*T_{1,1}*V_1'*V_2*T_{2,2}*V_2' -* -* Define T{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2} -* -* Then, we can define the matrix V as +* +* Define T_{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2} +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -314,15 +404,15 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL CLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * -* Compute T_{1,2} +* Compute T_{1,2} * T_{1,2} = V_{2,1}' * DO J = 1, L DO I = 1, K-L - T(J, L+I) = CONJG(V(L+I, J)) + T(J,L+I) = CONJG(V(L+I,J)) END DO END DO * @@ -335,9 +425,8 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = V_{3,1}'*V_{3,2} + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL CGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, - $ V(K+1, 1), LDV, V(K+1, L+1), LDV, ONE, T(1, L+1), - $ LDT) + CALL CGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, + $ V(K+1, 1), LDV, V(K+1,L+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1'*V_2 * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -350,8 +439,8 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * T_{1,2} = T_{1,2}*T_{2,2} * - CALL CTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + CALL CTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) ELSE IF(LQ) THEN * @@ -365,19 +454,19 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\C^{l,l} unit upper triangular * V_{1,2}\in\C^{l,k-l} rectangular * V_{1,3}\in\C^{l,n-k} rectangular -* +* * V_{2,2}\in\C^{k-l,k-l} unit upper triangular * V_{2,3}\in\C^{k-l,n-k} rectangular * * Where l = floor(k/2) * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -386,20 +475,20 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2}\in\C^{l, k-l} rectangular * * Then, consider the product: -* -* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) -* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 -* +* +* (I - V_1'*T_{1,1}'*V_1)*(I - V_2'*T_{2,2}'*V_2) +* = I - V_1'*T_{1,1}'*V_1 - V_2'*T_{2,2}'*V_2 + V_1'*T_{1,1}'*V_1*V_2'*T_{2,2}'*V_2 +* * Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{1,2} * @@ -409,14 +498,14 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL CLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * * Compute T_{1,2} * T_{1,2} = V_{1,2} * - CALL CLACPY('All', L, K-L, V(1, L+1), LDV, T(1, L+1), LDT) + CALL CLACPY('All', L, K - L, V(1,L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*V_{2,2}' * @@ -444,6 +533,99 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * CALL CTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) + ELSE IF(LQT) THEN +* +* Break V apart into 6 components +* +* V = |----------------------| +* |V_{1,1} V_{1,2} V{1,3}| +* |0 V_{2,2} V{2,3}| +* |----------------------| +* +* V_{1,1}\in\C^{l,l} unit upper triangular +* V_{1,2}\in\C^{l,k-l} rectangular +* V_{1,3}\in\C^{l,n-k} rectangular +* +* V_{2,2}\in\C^{k-l,k-l} unit upper triangular +* V_{2,3}\in\C^{k-l,n-k} rectangular +* +* Where l = floor(k/2) +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} 0 | +* |T_{2,1} T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\C^{l, l} lower triangular +* T_{2,2}\in\C^{k-l, k-l} lower triangular +* T_{2,1}\in\C^{k-l, l} rectangular +* +* Then, consider the product: +* +* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) +* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 +* +* Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{2,1} +* +* Compute T_{1,1} recursively +* + CALL CLARFT(DIRECT, STOREV, N, L, V, LDV, TAU, T, LDT) +* +* Compute T_{2,2} recursively +* + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) + +* +* Compute T_{2,1} +* T_{2,1} = V_{1,2}' +* + DO I = 1, K-L + DO J = 1, L + T(L+I,J) = CONJG(V(J,L+I)) + END DO + END DO +* +* T_{2,1} = V_{2,2}*T_{2,1} +* + CALL CTRMM('Left', 'Upper', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(L+1,L+1), LDV, T(L+1,1), LDT) +* +* T_{2,1} = V_{2,3}*V_{1,3}' + T_{2,1} +* Note: We assume K <= N, and GEMM will do nothing if N=K +* + CALL CGEMM('No Transpose', 'Conjugate', K-L, L, N-K, ONE, + $ V(L+1,K+1), LDV, V(1, K+1), LDV, ONE, T(L+1,1), LDT) +* +* At this point, we have that T_{2,1} = V_2*V_1' +* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} +* respectively. +* +* T_{2,1} = -T_{2,2}*T_{2,1} +* + CALL CTRMM('Left', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, NEG_ONE, T(L+1,L+1), LDT, T(L+1,1), LDT) +* +* T_{2,1} = T_{2,1}*T_{1,1} +* + CALL CTRMM('Right', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, ONE, T, LDT, T(L+1,1), LDT) ELSE IF(QL) THEN * * Break V apart into 6 components @@ -456,18 +638,18 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * V_{1,1}\in\C^{n-k,k-l} rectangular * V_{2,1}\in\C^{k-l,k-l} unit upper triangular -* +* * V_{1,2}\in\C^{n-k,l} rectangular * V_{2,2}\in\C^{k-l,l} rectangular * V_{3,2}\in\C^{l,l} unit upper triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -478,17 +660,17 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_2*T_{2,2}*V_2')*(I - V_1*T_{1,1}*V_1') * = I - V_2*T_{2,2}*V_2' - V_1*T_{1,1}*V_1' + V_2*T_{2,2}*V_2'*V_1*T_{1,1}*V_1' -* +* * Define T_{2,1} = -T_{2,2}*V_2'*V_1*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -496,34 +678,34 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{1,1} recursively * - CALL CLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * CALL CLARFT(DIRECT, STOREV, N, L, V(1, K-L+1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2}' * DO J = 1, K-L DO I = 1, L - T(K-L+I, J) = CONJG(V(N-K+J, K-L+I)) + T(K-L+I,J) = CONJG(V(N-K+J, K-L+I)) END DO END DO * * T_{2,1} = T_{2,1}*V_{2,1} * CALL CTRMM('Right', 'Upper', 'No transpose', 'Unit', L, - $ K-L, ONE, V(N-K+1, 1), LDV, T(K-L+1, 1), LDT) + $ K-L, ONE, V(N-K+1,1), LDV, T(K-L+1,1), LDT) * * T_{2,1} = V_{2,2}'*V_{2,1} + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL CGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, - $ V(1, K-L+1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + $ V(1,K-L+1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2'*V_1 * All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} @@ -532,17 +714,13 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{2,1} = -T_{2,2}*T_{2,1} * CALL CTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * CALL CTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) - ELSE -* -* Else means RQ case -* + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQ) THEN * Break V apart into 6 components * * V = |-----------------------| @@ -557,13 +735,13 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{2,2}\in\C^{l,k-l} rectangular * V_{2,3}\in\C^{l,l} unit lower triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -574,51 +752,51 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* -* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) -* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 -* +* +* (I - V_2'*T_{2,2}'*V_2)*(I - V_1'*T_{1,1}'*V_1) +* = I - V_2'*T_{2,2}'*V_2 - V_1'*T_{1,1}'*V_1 + V_2'*T_{2,2}'*V_2*V_1'*T_{1,1}'*V_1 +* * Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{2,1} * * Compute T_{1,1} recursively * - CALL CLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * CALL CLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2} * - CALL CLACPY('All', L, K-L, V(K-L+1, N-K+1), LDV, - $ T(K-L+1, 1), LDT) + CALL CLACPY('All', L, K-L, V(K-L+1,N-K+1), LDV, T(K-L+1,1), + $ LDT) * * T_{2,1} = T_{2,1}*V_{1,2}' * CALL CTRMM('Right', 'Lower', 'Conjugate', 'Unit', L, K-L, - $ ONE, V(1, N-K+1), LDV, T(K-L+1,1), LDT) + $ ONE, V(1, N-K+1), LDV, T(K-L+1,1), LDT) * * T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL CGEMM('No transpose', 'Conjugate', L, K-L, N-K, ONE, - $ V(K-L+1, 1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + CALL CGEMM('No transpose', 'Conjugate', L, K-L, N-K, ONE, + $ V(K-L+1,1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2*V_1' @@ -627,14 +805,104 @@ RECURSIVE SUBROUTINE CLARFT( DIRECT, STOREV, N, K, V, LDV, * * T_{2,1} = -T_{2,2}*T_{2,1} * - CALL CTRMM('Left', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + CALL CTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * - CALL CTRMM('Right', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) + CALL CTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQT) THEN +* Break V apart into 6 components +* +* V = |-----------------------| +* |V_{1,1} V_{1,2} 0 | +* |V_{2,1} V_{2,2} V_{2,3}| +* |-----------------------| +* +* V_{1,1}\in\C^{k-l,n-k} rectangular +* V_{1,2}\in\C^{k-l,k-l} unit lower triangular +* +* V_{2,1}\in\C^{l,n-k} rectangular +* V_{2,2}\in\C^{l,k-l} rectangular +* V_{2,3}\in\C^{l,l} unit lower triangular +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} T_{1,2}| +* | 0 T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\C^{k-l, k-l} non-unit upper triangular +* T_{2,2}\in\C^{l, l} non-unit upper triangular +* T_{1,2}\in\C^{k-l, l} rectangular +* +* Where l = floor(k/2) +* +* Then, consider the product: +* +* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) +* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 +* +* Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{1,2} +* +* Compute T_{1,1} recursively +* + CALL CLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) +* +* Compute T_{2,2} recursively +* + CALL CLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) +* +* Compute T_{1,2} +* T_{1,2} = V_{2,2}' +* + DO I = 1, K-L + DO J = 1, L + T(I,K-L+J) = CONJG(V(K-L+J, N-K+I)) + END DO + END DO +* +* T_{1,2} = V_{1,2}T_{1,2} +* + CALL CTRMM('Left', 'Lower', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(1,N-K+1), LDV, T(1,K-L+1), LDT) +* +* T_{1,2} = V_{1,1}V_{2,1}' + T_{1,2} +* + CALL CGEMM('No Transpose', 'Conjugate', K-L, L, N-K, ONE, V, + $ LDV, V(K-L+1,1), LDV, ONE, T(1, K-L+1), LDT) +* +* At this point, we have that T_{1,2} = V_1*V_2' +* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} +* respectively. +* +* T_{1,2} = -T_{1,1}*T_{1,2} +* + CALL CTRMM('Left', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, NEG_ONE, T, LDT, T(1, K-L+1), LDT) +* +* T_{1,2} = T_{1,2}*T_{2,2} +* + CALL CTRMM('Right', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, ONE, T(K-L+1,K-L+1), LDT, T(1, K-L+1), LDT) END IF END SUBROUTINE diff --git a/SRC/clarft_lvl2.f b/SRC/clarft_lvl2.f index 3b0aea113..0dca6c0dc 100644 --- a/SRC/clarft_lvl2.f +++ b/SRC/clarft_lvl2.f @@ -1,23 +1,14 @@ -*> \brief \b CLARFT_LVL2: Level 2 BLAS version for terminating case of CLARFT +*> \brief \b CLARFT_LVL2 forms the triangular factor T of a block reflector H = I - vtvH * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> Download CLARFT_LVL2 + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -* * Definition: * =========== * -* SUBROUTINE CLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, -* T, LDT ) +* SUBROUTINE CLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT ) * * .. Scalar Arguments .. * CHARACTER DIRECT, STOREV @@ -49,6 +40,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**H * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -61,6 +55,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -68,8 +64,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -155,6 +153,29 @@ *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -170,159 +191,568 @@ SUBROUTINE CLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. +* COMPLEX T( LDT, * ), TAU( * ), V( LDV, * ) * .. * -* ===================================================================== -* * .. Parameters .. - COMPLEX ONE, ZERO - PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ), - $ ZERO = ( 0.0E+0, 0.0E+0 ) ) -* .. +* + COMPLEX ONE + PARAMETER(ONE=(1.0E+0,0.0E+0)) +* * .. Local Scalars .. - INTEGER I, J, PREVLASTV, LASTV -* .. +* + INTEGER I,J,KMI,NMI,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET +* * .. External Subroutines .. - EXTERNAL CGEMM, CGEMV, CTRMV -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME +* + EXTERNAL CTRMV,CGEMV,CGEMM,XERBLA +* +* .. External Functions.. +* + LOGICAL LSAME + EXTERNAL LSAME +* +* .. Intrinsic Functions.. +* + INTRINSIC CONJG * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('CLARFT_LVL2', INFO) + RETURN + END IF +* * Quick return if possible * - IF( N.EQ.0 ) - $ RETURN -* - IF( LSAME( DIRECT, 'F' ) ) THEN - PREVLASTV = N - DO I = 1, K - PREVLASTV = MAX( PREVLASTV, I ) - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = 1, I - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) -* - CALL CGEMV( 'Conjugate transpose', J-I, I-1, - $ -TAU( I ), V( I+1, 1 ), LDV, - $ V( I+1, I ), 1, - $ ONE, T( 1, I ), 1 ) - ELSE -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( J , I ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H -* - CALL CGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), - $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, - $ ONE, T( 1, I ), LDT ) - END IF -* -* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) -* - CALL CTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, - $ T, - $ LDT, T( 1, I ), 1 ) - T( I, I ) = TAU( I ) - IF( I.GT.1 ) THEN - PREVLASTV = MAX( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF + IF(N.EQ.0.OR.K.EQ.0) THEN + RETURN + END IF +* +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER + IF( N.EQ.1.OR.K.EQ.1) THEN + IF( LQT.OR.RQT ) THEN + T(1,1) = CONJG(TAU(1)) + ELSE + T(1,1) = TAU(1) + END IF + RETURN + END IF + IF (QR) THEN +* +* Break V into 9 components +* +* V = |-----------------------| +* |V_{1,1} 0 0 | i-1 +* |V_{2,1} V_{2,2} 0 | 1 +* |V_{3,1} V_{3,2} V_{3,3}| n-i +* |-----------------------| +* i-1 1 k-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1}\\V_{2,1}\\V_{3,1}]' +* [0\\V_{2,2}\\V_{3,2}]T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{2,1}' + V_{3,1}'V_{3,2})T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{2,1}'T_{2,2} = -\tau_{i}V_{2,1}' +* T_{1,2} = -\tau_{i}V_{3,2}' V_{3,1} + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + + DO I = 2, K +* +* T_{1,2} = -V_{2,1}'V_{2,2}T_{2,2} = -\tau_i V_{2,1}' +* We must do this at copy time as otherwise gemv will do nothing +* on the last column when n=k, but we neet to make sure we are +* scaled by this value +* + DO J = 1, I-1 + T(J,I) = -TAU(I)*CONJG(V(I,J)) + END DO + +* +* T_{1,2} = -V_{3,1}'V_{3,2}T_{2,2} + T_{1,2} +* = -\tau{i} V_{3,2}'V_{3,1} + T_{1,2} +* + CALL CGEMV('Conjugate Transpose', N-I, I-1, -TAU(I), + $ V(I+1,1), LDV, V(I+1,I), 1, ONE, T(1, I), 1) + + +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL CTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) + +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - ELSE - PREVLASTV = 1 - DO I = K, 1, -1 - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = I, K - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( I.LT.K ) THEN - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) -* - CALL CGEMV( 'Conjugate transpose', N-K+I-J, K-I, - $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), - $ 1, ONE, T( I+1, I ), 1 ) - ELSE -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( J, N-K+I ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H -* - CALL CGEMM( 'N', 'C', K-I, 1, N-K+I-J, - $ -TAU( I ), - $ V( I+1, J ), LDV, V( I, J ), LDV, - $ ONE, T( I+1, I ), LDT ) - END IF -* -* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) -* - CALL CTRMV( 'Lower', 'No transpose', 'Non-unit', - $ K-I, - $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) - IF( I.GT.1 ) THEN - PREVLASTV = MIN( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF - T( I, I ) = TAU( I ) - END IF + ELSE IF (LQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1} V_{1,2} V_{1,3}][ 0 V_{2,2} V_{2,3} ]'T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{1,2} + V_{1,3}V_{2,3}')T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{1,2}T_{2,2} = -\tau_{i}V_{1,2} +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + + DO I = 2, K +* +* T_{1,2} = -\tau_{i}V_{1,2} +* + DO J = 1, I-1 + T(J, I) = -TAU(I)*V(J, I) + END DO + +* +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* + CALL CGEMM('No Transpose', 'Conjugate Transpose', I-1, + $ 1, N-I, -TAU(I), V(1,I+1), LDV, V(I, I+1), LDV, ONE, + $ T(1, I), LDT) +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL CTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) + +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - END IF - RETURN + ELSE IF (LQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | i-1 +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{2,1} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{3,1:3}. On +* each iteration, T_{3,1:3} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{2,1} = -T_{2,2}[0 V_{2,2} V_{2,3}][V_{1,1} V_{1,2} V_{1,3}]'T_{1,1} +* +* T_{2,1} = -T_{2,2}(V_{1,2}' + V_{2,3}V_{1,3}')T_{1,1} +* +* This means we will do the following +* +* T_{2,1} = -T_{2,2}V_{1,2}' = -\tau_{i}V_{1,2}' +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* T_{2,1} = T_{1,1}'T_{2,1} +* T_{2,2} = \tau{i} +* + T(1,1) = CONJG(TAU(1)) + + DO I = 2, K +* +* T_{2,1} = -\tau_{i}V_{1,2}' +* + DO J = 1, I-1 + T(I,J) = -CONJG(TAU(I)*V(J,I)) + END DO +* +* T_{2,1} = -\tau_{i}V_{2,3}V_{1,3}' + T_{2,1} +* + CALL CGEMM('No Transpose', 'Conjugate Transpose', 1, + $ I-1, N-I, -CONJG(TAU(I)), V(I,I+1), LDV, V(1, I+1), + $ LDV, ONE, T(I, 1), LDT) +* +* T_{2,1} = T_{1,1}'T_{2,1} +* + CALL CTRMV('Lower', 'Transpose', 'Non-unit', + $ I-1, T, LDT, T(I,1), LDT) + + T(I,I) = CONJG(TAU(I)) + END DO + ELSE IF (QL) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | n-i +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{1,3}\\V_{2,3}\\V_{3,3}]' +* [V_{1,2}\\V_{2,2}\\0]T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{1,3}'V_{1,2} + V_{2,3}')T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -T_{3,3}V_{3,2}' = -\tau_{k-i+1}V_{3,2}' +* T_{3,2} = -\tau_{k-i+1}V_{1,3}'V_{1,2} + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau(k-i+1) +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau(k-i+1)V_{2,3}' +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*CONJG(V(NMI, KMI + J)) + END DO +* +* T_{3,2} = -\tau(k-i+1)V_{1,3}'V_{1,2} + T_{3,2} +* + CALL CGEMV('Conjugate Transpose', N-I, I-1, -TAU(KMI), + $ V(1, KMI + 1), LDV, V(1, KMI), 1, ONE, + $ T(KMI+1, KMI), 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL CTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI + 1, KMI + 1), LDT, T(KMI + 1, KMI), 1) + + END DO + ELSE IF (RQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{3,1} V_{3,2} V_{3,3}][V_{2,1} V_{2,2} 0]'T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{3,1}V_{2,1}' + V_{3,2})T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* + CALL CGEMM('No Transpose', 'Conjugate Transpose', I-1, + $ 1, N-I, -TAU(KMI), V(KMI+1, 1), LDV, V(KMI, 1), LDV, + $ ONE, T(KMI+1, KMI), LDT) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL CTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI+1, KMI), 1) + END DO + ELSE IF (RQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | k-i +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{2,3} at +* each iteration i = 2, \dots k, and then grow into T_{1,1:3}. On +* each iteration, T_{1,1:3} are not referenced. See clarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{2,3} = -T_{2,2}[V_{2,1} V_{2,2} 0][V_{3,1} V_{3,2} V_{3,3}]'T_{3,3} +* +* T_{3,2} = -T_{2,2}(V_{2,1}V_{3,1}' + V_{3,2}')T_{3,3} * -* End of CLARFT_LVL2 +* Thus, we will compute * - END +* T_{2,2} = \tau_{k-i+1} +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* T_{2,3} = -\tau_{k-i+1}V_{2,1}V_{3,1}' + T_{2,3} +* T_{2,3} = T_{3,3}'T_{2,3} +* + T(K,K) = CONJG(TAU(K)) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = CONJG(TAU(KMI)) +* +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* + DO J = 1, I-1 + T(KMI, KMI + J) = -CONJG(TAU(KMI)*V(KMI + J, NMI)) + END DO +* +* T_{2,3} = -\tau_{k-i+1}V_{2,1}V_{3,1}' + T_{2,3} +* + CALL CGEMM('No Transpose', 'Conjugate Transpose', 1, + $ I-1, N-I, -CONJG(TAU(KMI)), V(KMI, 1), LDV, + $ V(KMI+1,1), LDV, ONE, T(KMI, KMI+1), LDT) +* +* T_{2,3} = T_{3,3}'T_{2,3} +* + CALL CTRMV('Upper', 'Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI, KMI+1), LDT) + END DO + END IF + END SUBROUTINE diff --git a/SRC/dlarft.f b/SRC/dlarft.f index 567b1e52e..154f7226d 100644 --- a/SRC/dlarft.f +++ b/SRC/dlarft.f @@ -48,6 +48,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**T * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -60,6 +63,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -67,8 +72,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -139,7 +146,7 @@ *> the H(i) is best illustrated by the following example with n = 5 and *> k = 3. The elements equal to 1 are not stored. *> -*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and (STOREV = 'R' or STOREV = 'T'): *> *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) *> ( v1 1 ) ( 1 v2 v2 v2 ) @@ -147,13 +154,36 @@ *> ( v1 v2 v3 ) *> ( v1 v2 v3 ) *> -*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> DIRECT = 'B' and STOREV = 'C': (DIRECT = 'B' or DIRECT = 'T') and STOREV = 'R': *> *> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) *> ( v1 v2 v3 ) ( v2 v2 v2 1 ) *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -177,24 +207,26 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * .. Parameters .. * - DOUBLE PRECISION ONE, NEG_ONE, ZERO - PARAMETER(ONE=1.0D+0, ZERO = 0.0D+0, NEG_ONE=-1.0D+0) + DOUBLE PRECISION ONE, NEG_ONE + PARAMETER(ONE=1.0D+0, NEG_ONE=-1.0D+0) * * .. Local Scalars .. * - INTEGER I,J,L,NX - LOGICAL QR,LQ,QL,DIRF,COLV + INTEGER I,J,L,NX,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET * * .. External Subroutines .. * - EXTERNAL DTRMM, DGEMM, DLACPY, DLARFT_LVL2 + EXTERNAL DLARFT_LVL2, DTRMM, DGEMM, DLACPY, XERBLA * * .. External Functions.. * LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV -* +* * The general scheme used is inspired by the approach inside DGEQRT3 * which was (at the time of writing this code): * Based on the algorithm of Elmroth and Gustavson, @@ -202,12 +234,90 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('DLARFT', INFO) + RETURN + END IF +* * Quick return if possible * IF(N.EQ.0.OR.K.EQ.0) THEN RETURN END IF * +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER +* * Base case * IF(N.EQ.1.OR.K.EQ.1) THEN @@ -215,44 +325,20 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, RETURN END IF * -* Determine when to cross over into the level 2 based implementation +* Determine crossover point from level 2 to level 3 BLAS implementation * NX = ILAENV(3, "DLARFT", DIRECT // STOREV, N, K, -1, -1) IF(K.LT.NX) THEN +* +* Finish this component with a level 2 BLAS implementation +* CALL DLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) RETURN END IF * -* Beginning of executable statements +* Beginning of executable statements for the recursive case * L = K / 2 -* -* Determine what kind of Q we need to compute -* We assume that if the user doesn't provide 'F' for DIRECT, -* then they meant to provide 'B' and if they don't provide -* 'C' for STOREV, then they meant to provide 'R' -* - DIRF = LSAME(DIRECT,'F') - COLV = LSAME(STOREV,'C') -* -* QR happens when we have forward direction in column storage -* - QR = DIRF.AND.COLV -* -* LQ happens when we have forward direction in row storage -* - LQ = DIRF.AND.(.NOT.COLV) -* -* QL happens when we have backward direction in column storage -* - QL = (.NOT.DIRF).AND.COLV -* -* The last case is RQ. Due to how we structured this, if the -* above 3 are false, then RQ must be true, so we never store -* this -* RQ happens when we have backward direction in row storage -* RQ = (.NOT.DIRF).AND.(.NOT.COLV) -* IF(QR) THEN * * Break V apart into 6 components @@ -266,17 +352,17 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\R^{l,l} unit lower triangular * V_{2,1}\in\R^{k-l,l} rectangular * V_{3,1}\in\R^{n-k,l} rectangular -* +* * V_{2,2}\in\R^{k-l,k-l} unit lower triangular * V_{3,2}\in\R^{n-k,k-l} rectangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -287,17 +373,17 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_1*T_{1,1}*V_1')*(I - V_2*T_{2,2}*V_2') * = I - V_1*T_{1,1}*V_1' - V_2*T_{2,2}*V_2' + V_1*T_{1,1}*V_1'*V_2*T_{2,2}*V_2' -* +* * Define T_{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -309,30 +395,29 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL DLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * -* Compute T_{1,2} +* Compute T_{1,2} * T_{1,2} = V_{2,1}' * DO J = 1, L DO I = 1, K-L - T(J, L+I) = V(L+I, J) + T(J,L+I) = V(L+I,J) END DO END DO * * T_{1,2} = T_{1,2}*V_{2,2} * CALL DTRMM('Right', 'Lower', 'No transpose', 'Unit', L, - $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{3,1}'*V_{3,2} + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL DGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, - $ V(K+1, 1), LDV, V(K+1, L+1), LDV, ONE, - $ T(1, L+1), LDT) + CALL DGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, + $ V(K+1, 1), LDV, V(K+1,L+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1'*V_2 * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -341,12 +426,12 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL DTRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * - CALL DTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + CALL DTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) ELSE IF(LQ) THEN * @@ -360,19 +445,19 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\R^{l,l} unit upper triangular * V_{1,2}\in\R^{l,k-l} rectangular * V_{1,3}\in\R^{l,n-k} rectangular -* +* * V_{2,2}\in\R^{k-l,k-l} unit upper triangular * V_{2,3}\in\R^{k-l,n-k} rectangular * * Where l = floor(k/2) * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -381,20 +466,20 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2}\in\R^{l, k-l} rectangular * * Then, consider the product: -* -* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) -* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 -* +* +* (I - V_1'*T_{1,1}'*V_1)*(I - V_2'*T_{2,2}'*V_2) +* = I - V_1'*T_{1,1}'*V_1 - V_2'*T_{2,2}'*V_2 + V_1'*T_{1,1}'*V_1*V_2'*T_{2,2}'*V_2 +* * Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{1,2} * @@ -404,27 +489,26 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL DLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * * Compute T_{1,2} * T_{1,2} = V_{1,2} * - CALL DLACPY('All', L, K-L, V(1, L+1), LDV, T(1, L+1), LDT) + CALL DLACPY('All', L, K - L, V(1,L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*V_{2,2}' * CALL DTRMM('Right', 'Upper', 'Transpose', 'Unit', L, K-L, - $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{1,3}*V_{2,3}' + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL DGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, - $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, - $ T(1, L+1), LDT) + $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1*V_2' * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -433,13 +517,106 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL DTRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * CALL DTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) + ELSE IF(LQT) THEN +* +* Break V apart into 6 components +* +* V = |----------------------| +* |V_{1,1} V_{1,2} V{1,3}| +* |0 V_{2,2} V{2,3}| +* |----------------------| +* +* V_{1,1}\in\R^{l,l} unit upper triangular +* V_{1,2}\in\R^{l,k-l} rectangular +* V_{1,3}\in\R^{l,n-k} rectangular +* +* V_{2,2}\in\R^{k-l,k-l} unit upper triangular +* V_{2,3}\in\R^{k-l,n-k} rectangular +* +* Where l = floor(k/2) +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} 0 | +* |T_{2,1} T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\R^{l, l} lower triangular +* T_{2,2}\in\R^{k-l, k-l} lower triangular +* T_{2,1}\in\R^{k-l, l} rectangular +* +* Then, consider the product: +* +* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) +* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 +* +* Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{2,1} +* +* Compute T_{1,1} recursively +* + CALL DLARFT(DIRECT, STOREV, N, L, V, LDV, TAU, T, LDT) +* +* Compute T_{2,2} recursively +* + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) + +* +* Compute T_{2,1} +* T_{2,1} = V_{1,2}' +* + DO I = 1, K-L + DO J = 1, L + T(L+I,J) = V(J,L+I) + END DO + END DO +* +* T_{2,1} = V_{2,2}*T_{2,1} +* + CALL DTRMM('Left', 'Upper', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(L+1,L+1), LDV, T(L+1,1), LDT) +* +* T_{2,1} = V_{2,3}*V_{1,3}' + T_{2,1} +* Note: We assume K <= N, and GEMM will do nothing if N=K +* + CALL DGEMM('No Transpose', 'Transpose', K-L, L, N-K, ONE, + $ V(L+1,K+1), LDV, V(1, K+1), LDV, ONE, T(L+1,1), LDT) +* +* At this point, we have that T_{2,1} = V_2*V_1' +* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} +* respectively. +* +* T_{2,1} = -T_{2,2}*T_{2,1} +* + CALL DTRMM('Left', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, NEG_ONE, T(L+1,L+1), LDT, T(L+1,1), LDT) +* +* T_{2,1} = T_{2,1}*T_{1,1} +* + CALL DTRMM('Right', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, ONE, T, LDT, T(L+1,1), LDT) ELSE IF(QL) THEN * * Break V apart into 6 components @@ -452,18 +629,18 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * V_{1,1}\in\R^{n-k,k-l} rectangular * V_{2,1}\in\R^{k-l,k-l} unit upper triangular -* +* * V_{1,2}\in\R^{n-k,l} rectangular * V_{2,2}\in\R^{k-l,l} rectangular * V_{3,2}\in\R^{l,l} unit upper triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -474,17 +651,17 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_2*T_{2,2}*V_2')*(I - V_1*T_{1,1}*V_1') * = I - V_2*T_{2,2}*V_2' - V_1*T_{1,1}*V_1' + V_2*T_{2,2}*V_2'*V_1*T_{1,1}*V_1' -* +* * Define T_{2,1} = -T_{2,2}*V_2'*V_1*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -492,34 +669,34 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{1,1} recursively * - CALL DLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * CALL DLARFT(DIRECT, STOREV, N, L, V(1, K-L+1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2}' * DO J = 1, K-L DO I = 1, L - T(K-L+I, J) = V(N-K+J, K-L+I) + T(K-L+I,J) = V(N-K+J, K-L+I) END DO END DO * * T_{2,1} = T_{2,1}*V_{2,1} * CALL DTRMM('Right', 'Upper', 'No transpose', 'Unit', L, - $ K-L, ONE, V(N-K+1, 1), LDV, T(K-L+1, 1), LDT) + $ K-L, ONE, V(N-K+1,1), LDV, T(K-L+1,1), LDT) * * T_{2,1} = V_{2,2}'*V_{2,1} + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL DGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, - $ V(1, K-L+1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + $ V(1,K-L+1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2'*V_1 * All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} @@ -528,17 +705,13 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{2,1} = -T_{2,2}*T_{2,1} * CALL DTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * CALL DTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) - ELSE -* -* Else means RQ case -* + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQ) THEN * Break V apart into 6 components * * V = |-----------------------| @@ -553,13 +726,13 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{2,2}\in\R^{l,k-l} rectangular * V_{2,3}\in\R^{l,l} unit lower triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -570,51 +743,51 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* -* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) -* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 -* +* +* (I - V_2'*T_{2,2}'*V_2)*(I - V_1'*T_{1,1}'*V_1) +* = I - V_2'*T_{2,2}'*V_2 - V_1'*T_{1,1}'*V_1 + V_2'*T_{2,2}'*V_2*V_1'*T_{1,1}'*V_1 +* * Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{2,1} * * Compute T_{1,1} recursively * - CALL DLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * - CALL DLARFT(DIRECT, STOREV, N, L, V(K-L+1, 1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + CALL DLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2} * - CALL DLACPY('All', L, K-L, V(K-L+1, N-K+1), LDV, - $ T(K-L+1, 1), LDT) + CALL DLACPY('All', L, K-L, V(K-L+1,N-K+1), LDV, T(K-L+1,1), + $ LDT) * * T_{2,1} = T_{2,1}*V_{1,2}' * CALL DTRMM('Right', 'Lower', 'Transpose', 'Unit', L, K-L, - $ ONE, V(1, N-K+1), LDV, T(K-L+1, 1), LDT) + $ ONE, V(1, N-K+1), LDV, T(K-L+1,1), LDT) * -* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} +* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL DGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, - $ V(K-L+1, 1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + CALL DGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, + $ V(K-L+1,1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2*V_1' @@ -623,14 +796,104 @@ RECURSIVE SUBROUTINE DLARFT( DIRECT, STOREV, N, K, V, LDV, * * T_{2,1} = -T_{2,2}*T_{2,1} * - CALL DTRMM('Left', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + CALL DTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * - CALL DTRMM('Right', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) + CALL DTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQT) THEN +* Break V apart into 6 components +* +* V = |-----------------------| +* |V_{1,1} V_{1,2} 0 | +* |V_{2,1} V_{2,2} V_{2,3}| +* |-----------------------| +* +* V_{1,1}\in\R^{k-l,n-k} rectangular +* V_{1,2}\in\R^{k-l,k-l} unit lower triangular +* +* V_{2,1}\in\R^{l,n-k} rectangular +* V_{2,2}\in\R^{l,k-l} rectangular +* V_{2,3}\in\R^{l,l} unit lower triangular +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} T_{1,2}| +* | 0 T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\R^{k-l, k-l} non-unit upper triangular +* T_{2,2}\in\R^{l, l} non-unit upper triangular +* T_{1,2}\in\R^{k-l, l} rectangular +* +* Where l = floor(k/2) +* +* Then, consider the product: +* +* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) +* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 +* +* Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{1,2} +* +* Compute T_{1,1} recursively +* + CALL DLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) +* +* Compute T_{2,2} recursively +* + CALL DLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) +* +* Compute T_{1,2} +* T_{1,2} = V_{2,2}' +* + DO I = 1, K-L + DO J = 1, L + T(I,K-L+J) = V(K-L+J, N-K+I) + END DO + END DO +* +* T_{1,2} = V_{1,2}T_{1,2} +* + CALL DTRMM('Left', 'Lower', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(1,N-K+1), LDV, T(1,K-L+1), LDT) +* +* T_{1,2} = V_{1,1}V_{2,1}' + T_{1,2} +* + CALL DGEMM('No Transpose', 'Transpose', K-L, L, N-K, ONE, V, + $ LDV, V(K-L+1,1), LDV, ONE, T(1, K-L+1), LDT) +* +* At this point, we have that T_{1,2} = V_1*V_2' +* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} +* respectively. +* +* T_{1,2} = -T_{1,1}*T_{1,2} +* + CALL DTRMM('Left', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, NEG_ONE, T, LDT, T(1, K-L+1), LDT) +* +* T_{1,2} = T_{1,2}*T_{2,2} +* + CALL DTRMM('Right', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, ONE, T(K-L+1,K-L+1), LDT, T(1, K-L+1), LDT) END IF END SUBROUTINE diff --git a/SRC/dlarft_lvl2.f b/SRC/dlarft_lvl2.f index 9614df466..df8b358d3 100644 --- a/SRC/dlarft_lvl2.f +++ b/SRC/dlarft_lvl2.f @@ -1,23 +1,15 @@ -*> \brief \b DLARFT_LVL2: Level 2 BLAS version for terminating case of DLARFT. +*> \brief \b DLARFT_LVL2 forms the triangular factor T of a block reflector H = I - vtvH * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> Download DLARFT_LVL2 + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -* * Definition: * =========== * * SUBROUTINE DLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, -* T, LDT ) +* $ T, LDT ) * * .. Scalar Arguments .. * CHARACTER DIRECT, STOREV @@ -49,6 +41,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**T * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -61,6 +56,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -68,8 +65,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -155,6 +154,29 @@ *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -170,157 +192,549 @@ SUBROUTINE DLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. +* DOUBLE PRECISION T( LDT, * ), TAU( * ), V( LDV, * ) * .. * -* ===================================================================== -* * .. Parameters .. - DOUBLE PRECISION ONE, ZERO - PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) -* .. +* + DOUBLE PRECISION ONE + PARAMETER(ONE=1.0D+0) +* * .. Local Scalars .. - INTEGER I, J, PREVLASTV, LASTV -* .. +* + INTEGER I,J,KMI,NMI,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET +* * .. External Subroutines .. - EXTERNAL DGEMV, DTRMV -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME +* + EXTERNAL DTRMV,DGEMV,XERBLA +* +* .. External Functions.. +* + LOGICAL LSAME + EXTERNAL LSAME * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('DLARFT_LVL2', INFO) + RETURN + END IF +* * Quick return if possible * - IF( N.EQ.0 ) - $ RETURN -* - IF( LSAME( DIRECT, 'F' ) ) THEN - PREVLASTV = N - DO I = 1, K - PREVLASTV = MAX( I, PREVLASTV ) - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = 1, I - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( I , J ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) -* - CALL DGEMV( 'Transpose', J-I, I-1, -TAU( I ), - $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, - $ T( 1, I ), 1 ) - ELSE -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( J , I ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T -* - CALL DGEMV( 'No transpose', I-1, J-I, -TAU( I ), - $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, ONE, - $ T( 1, I ), 1 ) - END IF -* -* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) -* - CALL DTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, - $ T, - $ LDT, T( 1, I ), 1 ) - T( I, I ) = TAU( I ) - IF( I.GT.1 ) THEN - PREVLASTV = MAX( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF + IF(N.EQ.0.OR.K.EQ.0) THEN + RETURN + END IF +* +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER + IF( N.EQ.1.OR.K.EQ.1) THEN + T(1,1) = TAU(1) + RETURN + END IF + IF (QR) THEN +* +* Break V into 9 components +* +* V = |-----------------------| +* |V_{1,1} 0 0 | i-1 +* |V_{2,1} V_{2,2} 0 | 1 +* |V_{3,1} V_{3,2} V_{3,3}| n-i +* |-----------------------| +* i-1 1 k-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1}\\V_{2,1}\\V_{3,1}]' +* [0\\V_{2,2}\\V_{3,2}]T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{2,1}' + V_{3,1}'V_{3,2})T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{2,1}'T_{2,2} = -\tau_{i}V_{2,1}' +* T_{1,2} = -\tau_{i}V_{3,2}' V_{3,1} + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + DO I = 2, K +* +* T_{1,2} = -V_{2,1}'V_{2,2}T_{2,2} = -\tau_i V_{2,1}' +* We must do this at copy time as otherwise gemv will do nothing +* on the last column when n=k, but we neet to make sure we are +* scaled by this value +* + DO J = 1, I-1 + T(J,I) = -V(I,J)*TAU(I) + END DO +* +* T_{1,2} = -V_{3,1}'V_{3,2}T_{2,2} + T_{1,2} +* = -\tau{i} V_{3,2}'V_{3,1} + T_{1,2} +* + CALL DGEMV('Transpose', N-I, I-1, -TAU(I), V(I+1,1), + $ LDV, V(I+1,I), 1, ONE, T(1, I), 1) +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL DTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - ELSE - PREVLASTV = 1 - DO I = K, 1, -1 - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = I, K - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( I.LT.K ) THEN - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( N-K+I , J ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) -* - CALL DGEMV( 'Transpose', N-K+I-J, K-I, - $ -TAU( I ), - $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, - $ T( I+1, I ), 1 ) - ELSE -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( J, N-K+I ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T -* - CALL DGEMV( 'No transpose', K-I, N-K+I-J, - $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, - $ ONE, T( I+1, I ), 1 ) - END IF -* -* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) -* - CALL DTRMV( 'Lower', 'No transpose', 'Non-unit', - $ K-I, - $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) - IF( I.GT.1 ) THEN - PREVLASTV = MIN( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF - T( I, I ) = TAU( I ) - END IF + ELSE IF (LQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1} V_{1,2} V_{1,3}][ 0 V_{2,2} V_{2,3} ]'T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{1,2} + V_{1,3}V_{2,3}')T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{1,2}T_{2,2} = -\tau_{i}V_{1,2} +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + DO I = 2, K +* +* T_{1,2} = -\tau_{i}V_{1,2} +* + DO J = 1, I-1 + T(J, I) = -TAU(I)*V(J, I) + END DO +* +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* + CALL DGEMV('No Transpose', I-1, N-I, -TAU(I), V(1, I+1), + $ LDV, V(I, I+1), LDV, ONE, T(1,I), 1) +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL DTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - END IF - RETURN + ELSE IF (LQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | i-1 +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{2,1} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{3,1:3}. On +* each iteration, T_{3,1:3} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{2,1} = -T_{2,2}[0 V_{2,2} V_{2,3}][V_{1,1} V_{1,2} V_{1,3}]'T_{1,1} +* +* T_{2,1} = -T_{2,2}(V_{1,2}' + V_{2,3}V_{1,3}')T_{1,1} +* +* This means we will do the following +* +* T_{2,1} = -T_{2,2}V_{1,2}' = -\tau_{i}V_{1,2}' +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* T_{2,1} = T_{1,1}'T_{2,1} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + DO I = 2, K +* +* T_{2,1} = -\tau_{i}V_{1,2}' +* + DO J = 1, I-1 + T(I,J) = -TAU(I)*V(J,I) + END DO +* +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* + CALL DGEMV('No transpose', I-1, N-I, -TAU(I), V(1, I+1), + $ LDV, V(I, I+1), LDV, ONE, T(I, 1), LDT) +* +* T_{2,1} = T_{1,1}'T_{2,1}' +* + CALL DTRMV('Lower', 'Transpose', 'Non-unit', I-1, + $ T, LDT, T(I,1), LDT) + T(I,I) = TAU(I) + END DO + ELSE IF (QL) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | n-i +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{1,3}\\V_{2,3}\\V_{3,3}]' +* [V_{1,2}\\V_{2,2}\\0]T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{1,3}'V_{1,2} + V_{2,3}')T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -T_{3,3}V_{3,2}' = -\tau_{k-i+1}V_{3,2}' +* T_{3,2} = -\tau_{k-i+1}V_{1,3}'V_{1,2} + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau(k-i+1) +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau(k-i+1)V_{2,3}' +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(NMI, KMI + J) + END DO +* +* T_{3,2} = -\tau(k-i+1)V_{1,3}'V_{1,2} + T_{3,2} +* + CALL DGEMV('Transpose', N-I, I-1, -TAU(KMI), + $ V(1, KMI + 1), LDV, V(1, KMI), 1, ONE, + $ T(KMI+1, KMI), 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL DTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI + 1, KMI + 1), LDT, T(KMI + 1, KMI), 1) + END DO + ELSE IF (RQ) THEN * -* End of DLARFT_LVL2 +* Break V into 9 components * - END +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{3,1} V_{3,2} V_{3,3}][V_{2,1} V_{2,2} 0]'T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{3,1}V_{2,1}' + V_{3,2})T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* + CALL DGEMV('No Transpose', I-1, N-I, -TAU(KMI), + $ V(KMI+1, 1), LDV, V(KMI, 1), LDV, ONE, T(KMI+1, KMI), + $ 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL DTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI+1, KMI), 1) + END DO + ELSE IF (RQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | k-i +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{2,3} at +* each iteration i = 2, \dots k, and then grow into T_{1,1:3}. On +* each iteration, T_{1,1:3} are not referenced. See dlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{2,3} = -T_{2,2}[V_{2,1} V_{2,2} 0][V_{3,1} V_{3,2} V_{3,3}]'T_{3,3} +* +* After transposing when necessary to fit our blas routines' interface, +* we get +* T_{3,2} = -T_{2,2}(V_{3,1}V_{2,1} + V_{3,2})T_{3,3} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* T_{2,3} = -\tau_{k-i+1}V_{3,1}V_{2,1} + T_{2,3} +* T_{2,3} = T_{3,3}'T_{2,3} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* + DO J = 1, I-1 + T(KMI, KMI + J) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{2,3} = -\tau_{k-i+1}V_{3,1}V_{2,1} + T_{2,3} +* + CALL DGEMV('No Transpose', I-1, N-I, -TAU(KMI), + $ V(KMI+1, 1), LDV, V(KMI, 1), LDV, ONE, T(KMI, KMI+1), + $ LDT) +* +* T_{2,3} = T_{3,3}'T_{2,3} +* + CALL DTRMV('Upper', 'Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI, KMI+1), LDT) + END DO + END IF + END SUBROUTINE diff --git a/SRC/slarft.f b/SRC/slarft.f index ba13c068b..a994a9249 100644 --- a/SRC/slarft.f +++ b/SRC/slarft.f @@ -48,6 +48,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**T * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -60,6 +63,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -67,8 +72,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -125,7 +132,7 @@ * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley -*> \author Johnathan Rhyne, Univ. of Colorado Denver (original author, 2024) +*> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \ingroup larft @@ -139,7 +146,7 @@ *> the H(i) is best illustrated by the following example with n = 5 and *> k = 3. The elements equal to 1 are not stored. *> -*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and STOREV = 'R': +*> DIRECT = 'F' and STOREV = 'C': DIRECT = 'F' and (STOREV = 'R' or STOREV = 'T'): *> *> V = ( 1 ) V = ( 1 v1 v1 v1 v1 ) *> ( v1 1 ) ( 1 v2 v2 v2 ) @@ -147,13 +154,36 @@ *> ( v1 v2 v3 ) *> ( v1 v2 v3 ) *> -*> DIRECT = 'B' and STOREV = 'C': DIRECT = 'B' and STOREV = 'R': +*> DIRECT = 'B' and STOREV = 'C': (DIRECT = 'B' or DIRECT = 'T') and STOREV = 'R': *> *> V = ( v1 v2 v3 ) V = ( v1 v1 1 ) *> ( v1 v2 v3 ) ( v2 v2 v2 1 ) *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -177,24 +207,26 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * .. Parameters .. * - REAL ONE, NEG_ONE, ZERO - PARAMETER(ONE=1.0E+0, ZERO = 0.0E+0, NEG_ONE=-1.0E+0) + REAL ONE, NEG_ONE + PARAMETER(ONE=1.0E+0, NEG_ONE=-1.0E+0) * * .. Local Scalars .. * - INTEGER I,J,L,NX - LOGICAL QR,LQ,QL,DIRF,COLV + INTEGER I,J,L,NX,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET * * .. External Subroutines .. * - EXTERNAL STRMM, SGEMM, SLACPY, SLARFT_LVL2 + EXTERNAL SLARFT_LVL2, STRMM, SGEMM, SLACPY, XERBLA * * .. External Functions.. * LOGICAL LSAME INTEGER ILAENV EXTERNAL LSAME, ILAENV -* +* * The general scheme used is inspired by the approach inside DGEQRT3 * which was (at the time of writing this code): * Based on the algorithm of Elmroth and Gustavson, @@ -202,12 +234,90 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('SLARFT', INFO) + RETURN + END IF +* * Quick return if possible * IF(N.EQ.0.OR.K.EQ.0) THEN RETURN END IF * +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER +* * Base case * IF(N.EQ.1.OR.K.EQ.1) THEN @@ -215,44 +325,20 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, RETURN END IF * -* Determine when to cross over into the level 2 based implementation +* Determine crossover point from level 2 to level 3 BLAS implementation * NX = ILAENV(3, "SLARFT", DIRECT // STOREV, N, K, -1, -1) IF(K.LT.NX) THEN +* +* Finish this component with a level 2 BLAS implementation +* CALL SLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) RETURN END IF * -* Beginning of executable statements +* Beginning of executable statements for the recursive case * L = K / 2 -* -* Determine what kind of Q we need to compute -* We assume that if the user doesn't provide 'F' for DIRECT, -* then they meant to provide 'B' and if they don't provide -* 'C' for STOREV, then they meant to provide 'R' -* - DIRF = LSAME(DIRECT,'F') - COLV = LSAME(STOREV,'C') -* -* QR happens when we have forward direction in column storage -* - QR = DIRF.AND.COLV -* -* LQ happens when we have forward direction in row storage -* - LQ = DIRF.AND.(.NOT.COLV) -* -* QL happens when we have backward direction in column storage -* - QL = (.NOT.DIRF).AND.COLV -* -* The last case is RQ. Due to how we structured this, if the -* above 3 are false, then RQ must be true, so we never store -* this -* RQ happens when we have backward direction in row storage -* RQ = (.NOT.DIRF).AND.(.NOT.COLV) -* IF(QR) THEN * * Break V apart into 6 components @@ -266,17 +352,17 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\R^{l,l} unit lower triangular * V_{2,1}\in\R^{k-l,l} rectangular * V_{3,1}\in\R^{n-k,l} rectangular -* +* * V_{2,2}\in\R^{k-l,k-l} unit lower triangular * V_{3,2}\in\R^{n-k,k-l} rectangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -287,17 +373,17 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_1*T_{1,1}*V_1')*(I - V_2*T_{2,2}*V_2') * = I - V_1*T_{1,1}*V_1' - V_2*T_{2,2}*V_2' + V_1*T_{1,1}*V_1'*V_2*T_{2,2}*V_2' -* +* * Define T_{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -309,30 +395,29 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL SLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * -* Compute T_{1,2} +* Compute T_{1,2} * T_{1,2} = V_{2,1}' * DO J = 1, L DO I = 1, K-L - T(J, L+I) = V(L+I, J) + T(J,L+I) = V(L+I,J) END DO END DO * * T_{1,2} = T_{1,2}*V_{2,2} * CALL STRMM('Right', 'Lower', 'No transpose', 'Unit', L, - $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{3,1}'*V_{3,2} + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL SGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, - $ V(K+1, 1), LDV, V(K+1, L+1), LDV, ONE, - $ T(1, L+1), LDT) + CALL SGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, + $ V(K+1, 1), LDV, V(K+1,L+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1'*V_2 * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -341,12 +426,12 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL STRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * - CALL STRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + CALL STRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) ELSE IF(LQ) THEN * @@ -360,19 +445,19 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\R^{l,l} unit upper triangular * V_{1,2}\in\R^{l,k-l} rectangular * V_{1,3}\in\R^{l,n-k} rectangular -* +* * V_{2,2}\in\R^{k-l,k-l} unit upper triangular * V_{2,3}\in\R^{k-l,n-k} rectangular * * Where l = floor(k/2) * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -381,20 +466,20 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2}\in\R^{l, k-l} rectangular * * Then, consider the product: -* -* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) -* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 -* +* +* (I - V_1'*T_{1,1}'*V_1)*(I - V_2'*T_{2,2}'*V_2) +* = I - V_1'*T_{1,1}'*V_1 - V_2'*T_{2,2}'*V_2 + V_1'*T_{1,1}'*V_1*V_2'*T_{2,2}'*V_2 +* * Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{1,2} * @@ -404,27 +489,26 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL SLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * * Compute T_{1,2} * T_{1,2} = V_{1,2} * - CALL SLACPY('All', L, K-L, V(1, L+1), LDV, T(1, L+1), LDT) + CALL SLACPY('All', L, K - L, V(1,L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*V_{2,2}' * CALL STRMM('Right', 'Upper', 'Transpose', 'Unit', L, K-L, - $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{1,3}*V_{2,3}' + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL SGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, - $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, - $ T(1, L+1), LDT) + $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1*V_2' * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -433,13 +517,106 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL STRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * CALL STRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) + ELSE IF(LQT) THEN +* +* Break V apart into 6 components +* +* V = |----------------------| +* |V_{1,1} V_{1,2} V{1,3}| +* |0 V_{2,2} V{2,3}| +* |----------------------| +* +* V_{1,1}\in\R^{l,l} unit upper triangular +* V_{1,2}\in\R^{l,k-l} rectangular +* V_{1,3}\in\R^{l,n-k} rectangular +* +* V_{2,2}\in\R^{k-l,k-l} unit upper triangular +* V_{2,3}\in\R^{k-l,n-k} rectangular +* +* Where l = floor(k/2) +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} 0 | +* |T_{2,1} T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\R^{l, l} lower triangular +* T_{2,2}\in\R^{k-l, k-l} lower triangular +* T_{2,1}\in\R^{k-l, l} rectangular +* +* Then, consider the product: +* +* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) +* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 +* +* Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{2,1} +* +* Compute T_{1,1} recursively +* + CALL SLARFT(DIRECT, STOREV, N, L, V, LDV, TAU, T, LDT) +* +* Compute T_{2,2} recursively +* + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) + +* +* Compute T_{2,1} +* T_{2,1} = V_{1,2}' +* + DO I = 1, K-L + DO J = 1, L + T(L+I,J) = V(J,L+I) + END DO + END DO +* +* T_{2,1} = V_{2,2}*T_{2,1} +* + CALL STRMM('Left', 'Upper', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(L+1,L+1), LDV, T(L+1,1), LDT) +* +* T_{2,1} = V_{2,3}*V_{1,3}' + T_{2,1} +* Note: We assume K <= N, and GEMM will do nothing if N=K +* + CALL SGEMM('No Transpose', 'Transpose', K-L, L, N-K, ONE, + $ V(L+1,K+1), LDV, V(1, K+1), LDV, ONE, T(L+1,1), LDT) +* +* At this point, we have that T_{2,1} = V_2*V_1' +* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} +* respectively. +* +* T_{2,1} = -T_{2,2}*T_{2,1} +* + CALL STRMM('Left', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, NEG_ONE, T(L+1,L+1), LDT, T(L+1,1), LDT) +* +* T_{2,1} = T_{2,1}*T_{1,1} +* + CALL STRMM('Right', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, ONE, T, LDT, T(L+1,1), LDT) ELSE IF(QL) THEN * * Break V apart into 6 components @@ -452,18 +629,18 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * V_{1,1}\in\R^{n-k,k-l} rectangular * V_{2,1}\in\R^{k-l,k-l} unit upper triangular -* +* * V_{1,2}\in\R^{n-k,l} rectangular * V_{2,2}\in\R^{k-l,l} rectangular * V_{3,2}\in\R^{l,l} unit upper triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -474,17 +651,17 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_2*T_{2,2}*V_2')*(I - V_1*T_{1,1}*V_1') * = I - V_2*T_{2,2}*V_2' - V_1*T_{1,1}*V_1' + V_2*T_{2,2}*V_2'*V_1*T_{1,1}*V_1' -* +* * Define T_{2,1} = -T_{2,2}*V_2'*V_1*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -492,34 +669,34 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{1,1} recursively * - CALL SLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * CALL SLARFT(DIRECT, STOREV, N, L, V(1, K-L+1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2}' * DO J = 1, K-L DO I = 1, L - T(K-L+I, J) = V(N-K+J, K-L+I) + T(K-L+I,J) = V(N-K+J, K-L+I) END DO END DO * * T_{2,1} = T_{2,1}*V_{2,1} * CALL STRMM('Right', 'Upper', 'No transpose', 'Unit', L, - $ K-L, ONE, V(N-K+1, 1), LDV, T(K-L+1, 1), LDT) + $ K-L, ONE, V(N-K+1,1), LDV, T(K-L+1,1), LDT) * * T_{2,1} = V_{2,2}'*V_{2,1} + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL SGEMM('Transpose', 'No transpose', L, K-L, N-K, ONE, - $ V(1, K-L+1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + $ V(1,K-L+1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2'*V_1 * All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} @@ -528,17 +705,13 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{2,1} = -T_{2,2}*T_{2,1} * CALL STRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * CALL STRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) - ELSE -* -* Else means RQ case -* + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQ) THEN * Break V apart into 6 components * * V = |-----------------------| @@ -553,13 +726,13 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{2,2}\in\R^{l,k-l} rectangular * V_{2,3}\in\R^{l,l} unit lower triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -570,51 +743,51 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* -* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) -* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 -* +* +* (I - V_2'*T_{2,2}'*V_2)*(I - V_1'*T_{1,1}'*V_1) +* = I - V_2'*T_{2,2}'*V_2 - V_1'*T_{1,1}'*V_1 + V_2'*T_{2,2}'*V_2*V_1'*T_{1,1}'*V_1 +* * Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'TV +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{2,1} * * Compute T_{1,1} recursively * - CALL SLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * - CALL SLARFT(DIRECT, STOREV, N, L, V(K-L+1, 1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + CALL SLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2} * - CALL SLACPY('All', L, K-L, V(K-L+1, N-K+1), LDV, - $ T(K-L+1, 1), LDT) + CALL SLACPY('All', L, K-L, V(K-L+1,N-K+1), LDV, T(K-L+1,1), + $ LDT) * * T_{2,1} = T_{2,1}*V_{1,2}' * CALL STRMM('Right', 'Lower', 'Transpose', 'Unit', L, K-L, - $ ONE, V(1, N-K+1), LDV, T(K-L+1, 1), LDT) + $ ONE, V(1, N-K+1), LDV, T(K-L+1,1), LDT) * -* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} +* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL SGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, - $ V(K-L+1, 1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + CALL SGEMM('No transpose', 'Transpose', L, K-L, N-K, ONE, + $ V(K-L+1,1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2*V_1' @@ -623,14 +796,104 @@ RECURSIVE SUBROUTINE SLARFT( DIRECT, STOREV, N, K, V, LDV, * * T_{2,1} = -T_{2,2}*T_{2,1} * - CALL STRMM('Left', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + CALL STRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * - CALL STRMM('Right', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) + CALL STRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQT) THEN +* Break V apart into 6 components +* +* V = |-----------------------| +* |V_{1,1} V_{1,2} 0 | +* |V_{2,1} V_{2,2} V_{2,3}| +* |-----------------------| +* +* V_{1,1}\in\R^{k-l,n-k} rectangular +* V_{1,2}\in\R^{k-l,k-l} unit lower triangular +* +* V_{2,1}\in\R^{l,n-k} rectangular +* V_{2,2}\in\R^{l,k-l} rectangular +* V_{2,3}\in\R^{l,l} unit lower triangular +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} T_{1,2}| +* | 0 T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\R^{k-l, k-l} non-unit upper triangular +* T_{2,2}\in\R^{l, l} non-unit upper triangular +* T_{1,2}\in\R^{k-l, l} rectangular +* +* Where l = floor(k/2) +* +* Then, consider the product: +* +* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) +* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 +* +* Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{1,2} +* +* Compute T_{1,1} recursively +* + CALL SLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) +* +* Compute T_{2,2} recursively +* + CALL SLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) +* +* Compute T_{1,2} +* T_{1,2} = V_{2,2}' +* + DO I = 1, K-L + DO J = 1, L + T(I,K-L+J) = V(K-L+J, N-K+I) + END DO + END DO +* +* T_{1,2} = V_{1,2}T_{1,2} +* + CALL STRMM('Left', 'Lower', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(1,N-K+1), LDV, T(1,K-L+1), LDT) +* +* T_{1,2} = V_{1,1}V_{2,1}' + T_{1,2} +* + CALL SGEMM('No Transpose', 'Transpose', K-L, L, N-K, ONE, V, + $ LDV, V(K-L+1,1), LDV, ONE, T(1, K-L+1), LDT) +* +* At this point, we have that T_{1,2} = V_1*V_2' +* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} +* respectively. +* +* T_{1,2} = -T_{1,1}*T_{1,2} +* + CALL STRMM('Left', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, NEG_ONE, T, LDT, T(1, K-L+1), LDT) +* +* T_{1,2} = T_{1,2}*T_{2,2} +* + CALL STRMM('Right', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, ONE, T(K-L+1,K-L+1), LDT, T(1, K-L+1), LDT) END IF END SUBROUTINE diff --git a/SRC/slarft_lvl2.f b/SRC/slarft_lvl2.f index 7107a91d5..2a4f5bd8e 100644 --- a/SRC/slarft_lvl2.f +++ b/SRC/slarft_lvl2.f @@ -1,18 +1,10 @@ -*> \brief \b SLARFT_LVL2: Level 2 BLAS version for terminating case of SLARFT. +*> \brief \b SLARFT_LVL2 forms the triangular factor T of a block reflector H = I - vtvH * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> Download SLARFT_LVL2 + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -* * Definition: * =========== * @@ -49,6 +41,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**T * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -61,6 +56,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -68,8 +65,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -155,6 +154,29 @@ *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -170,157 +192,555 @@ SUBROUTINE SLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. +* REAL T( LDT, * ), TAU( * ), V( LDV, * ) * .. * -* ===================================================================== -* * .. Parameters .. - REAL ONE, ZERO - PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) -* .. +* + REAL ONE + PARAMETER(ONE=1.0E+0) +* * .. Local Scalars .. - INTEGER I, J, PREVLASTV, LASTV -* .. +* + INTEGER I,J,KMI,NMI,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET +* * .. External Subroutines .. - EXTERNAL SGEMV, STRMV -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME +* + EXTERNAL STRMV,SGEMV,XERBLA +* +* .. External Functions.. +* + LOGICAL LSAME + EXTERNAL LSAME * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('SLARFT_LVL2', INFO) + RETURN + END IF +* * Quick return if possible * - IF( N.EQ.0 ) - $ RETURN -* - IF( LSAME( DIRECT, 'F' ) ) THEN - PREVLASTV = N - DO I = 1, K - PREVLASTV = MAX( I, PREVLASTV ) - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = 1, I - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( I , J ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i) -* - CALL SGEMV( 'Transpose', J-I, I-1, -TAU( I ), - $ V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE, - $ T( 1, I ), 1 ) - ELSE -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( J , I ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T -* - CALL SGEMV( 'No transpose', I-1, J-I, -TAU( I ), - $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, - $ ONE, T( 1, I ), 1 ) - END IF -* -* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) -* - CALL STRMV( 'Upper', 'No transpose', 'Non-unit', I-1, - $ T, - $ LDT, T( 1, I ), 1 ) - T( I, I ) = TAU( I ) - IF( I.GT.1 ) THEN - PREVLASTV = MAX( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF + IF(N.EQ.0.OR.K.EQ.0) THEN + RETURN + END IF +* +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER + IF( N.EQ.1.OR.K.EQ.1) THEN + T(1,1) = TAU(1) + RETURN + END IF + IF (QR) THEN +* +* Break V into 9 components +* +* V = |-----------------------| +* |V_{1,1} 0 0 | i-1 +* |V_{2,1} V_{2,2} 0 | 1 +* |V_{3,1} V_{3,2} V_{3,3}| n-i +* |-----------------------| +* i-1 1 k-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1}\\V_{2,1}\\V_{3,1}]' +* [0\\V_{2,2}\\V_{3,2}]T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{2,1}' + V_{3,1}'V_{3,2})T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{2,1}'T_{2,2} = -\tau_{i}V_{2,1}' +* T_{1,2} = -\tau_{i}V_{3,2}' V_{3,1} + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + + DO I = 2, K +* +* T_{1,2} = -V_{2,1}'V_{2,2}T_{2,2} = -\tau_i V_{2,1}' +* We must do this at copy time as otherwise gemv will do nothing +* on the last column when n=k, but we neet to make sure we are +* scaled by this value +* + DO J = 1, I-1 + T(J,I) = -V(I,J)*TAU(I) + END DO + +* +* T_{1,2} = -V_{3,1}'V_{3,2}T_{2,2} + T_{1,2} +* = -\tau{i} V_{3,2}'V_{3,1} + T_{1,2} +* + CALL SGEMV('Transpose', N-I, I-1, -TAU(I), V(I+1,1), + $ LDV, V(I+1,I), 1, ONE, T(1, I), 1) + + +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL STRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) + +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - ELSE - PREVLASTV = 1 - DO I = K, 1, -1 - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = I, K - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( I.LT.K ) THEN - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( N-K+I , J ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i) -* - CALL SGEMV( 'Transpose', N-K+I-J, K-I, - $ -TAU( I ), - $ V( J, I+1 ), LDV, V( J, I ), 1, ONE, - $ T( I+1, I ), 1 ) - ELSE -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( J, N-K+I ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T -* - CALL SGEMV( 'No transpose', K-I, N-K+I-J, - $ -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV, - $ ONE, T( I+1, I ), 1 ) - END IF -* -* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) -* - CALL STRMV( 'Lower', 'No transpose', 'Non-unit', - $ K-I, - $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) - IF( I.GT.1 ) THEN - PREVLASTV = MIN( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF - T( I, I ) = TAU( I ) - END IF + ELSE IF (LQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1} V_{1,2} V_{1,3}][ 0 V_{2,2} V_{2,3} ]'T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{1,2} + V_{1,3}V_{2,3}')T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{1,2}T_{2,2} = -\tau_{i}V_{1,2} +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + DO I = 2, K +* +* T_{1,2} = -\tau_{i}V_{1,2} +* + DO J = 1, I-1 + T(J, I) = -TAU(I)*V(J, I) + END DO +* +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* + CALL SGEMV('No Transpose', I-1, N-I, -TAU(I), V(1, I+1), + $ LDV, V(I, I+1), LDV, ONE, T(1,I), 1) +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL STRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - END IF - RETURN + ELSE IF (LQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | i-1 +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{2,1} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{3,1:3}. On +* each iteration, T_{3,1:3} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{2,1} = -T_{2,2}[0 V_{2,2} V_{2,3}][V_{1,1} V_{1,2} V_{1,3}]'T_{1,1} +* +* T_{2,1} = -T_{2,2}(V_{1,2}' + V_{2,3}V_{1,3}')T_{1,1} +* +* This means we will do the following +* +* T_{2,1} = -T_{2,2}V_{1,2}' = -\tau_{i}V_{1,2}' +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* T_{2,1} = T_{1,1}'T_{2,1} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + DO I = 2, K +* +* T_{2,1} = -\tau_{i}V_{1,2}' +* + DO J = 1, I-1 + T(I,J) = -TAU(I)*V(J,I) + END DO +* +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* + CALL SGEMV('No transpose', I-1, N-I, -TAU(I), V(1, I+1), + $ LDV, V(I, I+1), LDV, ONE, T(I, 1), LDT) + +* +* T_{2,1}' = T_{1,1}'T_{2,1}' +* + CALL STRMV('Lower', 'Transpose', 'Non-unit', I-1, + $ T, LDT, T(I,1), LDT) + T(I,I) = TAU(I) + END DO + ELSE IF (QL) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | n-i +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{1,3}\\V_{2,3}\\V_{3,3}]' +* [V_{1,2}\\V_{2,2}\\0]T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{1,3}'V_{1,2} + V_{2,3}')T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -T_{3,3}V_{3,2}' = -\tau_{k-i+1}V_{3,2}' +* T_{3,2} = -\tau_{k-i+1}V_{1,3}'V_{1,2} + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau(k-i+1) +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau(k-i+1)V_{2,3}' +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(NMI, KMI + J) + END DO +* +* T_{3,2} = -\tau(k-i+1)V_{1,3}'V_{1,2} + T_{3,2} +* + CALL SGEMV('Transpose', N-I, I-1, -TAU(KMI), + $ V(1, KMI + 1), LDV, V(1, KMI), 1, ONE, + $ T(KMI+1, KMI), 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL STRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI + 1, KMI + 1), LDT, T(KMI + 1, KMI), 1) + END DO + ELSE IF (RQ) THEN * -* End of SLARFT_LVL2 +* Break V into 9 components * - END +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{3,1} V_{3,2} V_{3,3}][V_{2,1} V_{2,2} 0]'T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{3,1}V_{2,1}' + V_{3,2})T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* + CALL SGEMV('No Transpose', I-1, N-I, -TAU(KMI), + $ V(KMI+1, 1), LDV, V(KMI, 1), LDV, ONE, T(KMI+1, KMI), + $ 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL STRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI+1, KMI), 1) + END DO + ELSE IF (RQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | k-i +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{2,3} at +* each iteration i = 2, \dots k, and then grow into T_{1,1:3}. On +* each iteration, T_{1,1:3} are not referenced. See slarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{2,3} = -T_{2,2}[V_{2,1} V_{2,2} 0][V_{3,1} V_{3,2} V_{3,3}]'T_{3,3} +* +* After transposing when necessary to fit our blas routines' interface, +* we get +* T_{3,2} = -T_{2,2}(V_{3,1}V_{2,1} + V_{3,2})T_{3,3} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* T_{2,3} = -\tau_{k-i+1}V_{3,1}V_{2,1} + T_{2,3} +* T_{2,3} = T_{3,3}'T_{2,3} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* + DO J = 1, I-1 + T(KMI, KMI + J) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{2,3} = -\tau_{k-i+1}V_{3,1}V_{2,1} + T_{2,3} +* + CALL SGEMV('No Transpose', I-1, N-I, -TAU(KMI), + $ V(KMI+1, 1), LDV, V(KMI, 1), LDV, ONE, T(KMI, KMI+1), + $ LDT) +* +* T_{2,3} = T_{3,3}'T_{2,3} +* + CALL STRMV('Upper', 'Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI, KMI+1), LDT) + END DO + END IF + END SUBROUTINE diff --git a/SRC/zlarft.f b/SRC/zlarft.f index 626a3c4b9..8a0d3f0cf 100644 --- a/SRC/zlarft.f +++ b/SRC/zlarft.f @@ -48,6 +48,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**H * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -60,6 +63,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -67,8 +72,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -154,6 +161,29 @@ *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -167,28 +197,30 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * .. Scalar Arguments * - CHARACTER DIRECT, STOREV - INTEGER K, LDT, LDV, N + CHARACTER DIRECT, STOREV + INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. * - COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) + COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) * .. * * .. Parameters .. * - COMPLEX*16 ONE, NEG_ONE, ZERO - PARAMETER(ONE=(1.0D+0,0.0D+0), ZERO = (0.0D+0,0.0D+0), - $ NEG_ONE=(-1.0D+0,0.0D+0)) + COMPLEX*16 ONE, NEG_ONE + PARAMETER(ONE=(1.0D+0,0.0D+0), + $ NEG_ONE=(-1.0D+0,0.0D+0)) * * .. Local Scalars .. * - INTEGER I,J,L,NX - LOGICAL QR,LQ,QL,DIRF,COLV + INTEGER I,J,L,NX,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET * * .. External Subroutines .. * - EXTERNAL ZTRMM, ZGEMM, ZLACPY, ZLARFT_LVL2 + EXTERNAL ZLARFT_LVL2, ZTRMM, ZGEMM, ZLACPY, XERBLA * * .. External Functions.. * @@ -198,8 +230,8 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * .. Intrinsic Functions.. * - INTRINSIC CONJG -* + INTRINSIC DCONJG +* * The general scheme used is inspired by the approach inside DGEQRT3 * which was (at the time of writing this code): * Based on the algorithm of Elmroth and Gustavson, @@ -207,57 +239,115 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * .. * .. Executable Statements .. * -* Quick return if possible +* Convert our character flags to logical flags for later * - IF(N.EQ.0.OR.K.EQ.0) THEN - RETURN - END IF + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') * -* Base case +* Error handling for our character flags * - IF(N.EQ.1.OR.K.EQ.1) THEN - T(1,1) = TAU(1) + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('ZLARFT', INFO) RETURN END IF * -* Determine when to cross over into the level 2 based implementation +* Quick return if possible * - NX = ILAENV(3, "ZLARFT", DIRECT // STOREV, N, K, -1, -1) - IF(K.LT.NX) THEN - CALL ZLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) + IF(N.EQ.0.OR.K.EQ.0) THEN RETURN END IF * -* Beginning of executable statements -* - L = K / 2 +* Now we determine what factorization our flags are associated with * -* Determine what kind of Q we need to compute -* We assume that if the user doesn't provide 'F' for DIRECT, -* then they meant to provide 'B' and if they don't provide -* 'C' for STOREV, then they meant to provide 'R' +* QR happens when we have forward direction in column storage * - DIRF = LSAME(DIRECT,'F') - COLV = LSAME(STOREV,'C') + QR = DIRF.AND.STOREC * -* QR happens when we have forward direction in column storage +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute * - QR = DIRF.AND.COLV + LQT = DIRF.AND.STORET * -* LQ happens when we have forward direction in row storage +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute * - LQ = DIRF.AND.(.NOT.COLV) + LQ = DIRF.AND.STORER * * QL happens when we have backward direction in column storage * - QL = (.NOT.DIRF).AND.COLV + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER +* +* Base case +* + IF(N.EQ.1.OR.K.EQ.1) THEN + IF( LQT.OR.RQT ) THEN + T(1,1) = DCONJG(TAU(1)) + ELSE + T(1,1) = TAU(1) + END IF + RETURN + END IF +* +* Determine crossover point from level 2 to level 3 BLAS implementation +* + NX = ILAENV(3, "ZLARFT", DIRECT // STOREV, N, K, -1, -1) + IF(K.LT.NX) THEN * -* The last case is RQ. Due to how we structured this, if the -* above 3 are false, then RQ must be true, so we never store -* this -* RQ happens when we have backward direction in row storage -* RQ = (.NOT.DIRF).AND.(.NOT.COLV) +* Finish this component with a level 2 BLAS implementation * + CALL ZLARFT_LVL2(DIRECT, STOREV, N, K, V, LDV, TAU, T, LDT) + RETURN + END IF +* +* Beginning of executable statements for the recursive case +* + L = K / 2 IF(QR) THEN * * Break V apart into 6 components @@ -271,17 +361,17 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\C^{l,l} unit lower triangular * V_{2,1}\in\C^{k-l,l} rectangular * V_{3,1}\in\C^{n-k,l} rectangular -* +* * V_{2,2}\in\C^{k-l,k-l} unit lower triangular * V_{3,2}\in\C^{n-k,k-l} rectangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -292,17 +382,17 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_1*T_{1,1}*V_1')*(I - V_2*T_{2,2}*V_2') * = I - V_1*T_{1,1}*V_1' - V_2*T_{2,2}*V_2' + V_1*T_{1,1}*V_1'*V_2*T_{2,2}*V_2' -* +* * Define T_{1,2} = -T_{1,1}*V_1'*V_2*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -314,30 +404,29 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * -* Compute T_{1,2} +* Compute T_{1,2} * T_{1,2} = V_{2,1}' * DO J = 1, L DO I = 1, K-L - T(J, L+I) = CONJG(V(L+I, J)) + T(J,L+I) = CONJG(V(L+I,J)) END DO END DO * * T_{1,2} = T_{1,2}*V_{2,2} * CALL ZTRMM('Right', 'Lower', 'No transpose', 'Unit', L, - $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ K-L, ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{3,1}'*V_{3,2} + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL ZGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, - $ V(K+1, 1), LDV, V(K+1, L+1), LDV, ONE, - $ T(1, L+1), LDT) + CALL ZGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, + $ V(K+1, 1), LDV, V(K+1,L+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1'*V_2 * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -346,12 +435,12 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL ZTRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * - CALL ZTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + CALL ZTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) ELSE IF(LQ) THEN * @@ -365,41 +454,41 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{1,1}\in\C^{l,l} unit upper triangular * V_{1,2}\in\C^{l,k-l} rectangular * V_{1,3}\in\C^{l,n-k} rectangular -* +* * V_{2,2}\in\C^{k-l,k-l} unit upper triangular * V_{2,3}\in\C^{k-l,n-k} rectangular * * Where l = floor(k/2) * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} T_{1,2}| * |0 T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * -* T_{1,1}\in\C^{l, l} upper triangular -* T_{2,2}\in\C^{k-l, k-l} upper triangular -* T_{1,2}\in\C^{l, k-l} rectangular +* T_{1,1}\in\C^{l, l} upper triangular +* T_{2,2}\in\C^{k-l, k-l} upper triangular +* T_{1,2}\in\C^{l, k-l} rectangular * * Then, consider the product: -* -* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) -* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 -* +* +* (I - V_1'*T_{1,1}'*V_1)*(I - V_2'*T_{2,2}'*V_2) +* = I - V_1'*T_{1,1}'*V_1 - V_2'*T_{2,2}'*V_2 + V_1'*T_{1,1}'*V_1*V_2'*T_{2,2}'*V_2 +* * Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{1,2} * @@ -409,27 +498,26 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{2,2} recursively * - CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V(L+1, L+1), LDV, - $ TAU(L+1), T(L+1, L+1), LDT) + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) * * Compute T_{1,2} * T_{1,2} = V_{1,2} * - CALL ZLACPY('All', L, K-L, V(1, L+1), LDV, T(1, L+1), LDT) + CALL ZLACPY('All', L, K - L, V(1,L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*V_{2,2}' * CALL ZTRMM('Right', 'Upper', 'Conjugate', 'Unit', L, K-L, - $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) + $ ONE, V(L+1, L+1), LDV, T(1, L+1), LDT) * * T_{1,2} = V_{1,3}*V_{2,3}' + T_{1,2} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL ZGEMM('No transpose', 'Conjugate', L, K-L, N-K, ONE, - $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, - $ T(1, L+1), LDT) + $ V(1, K+1), LDV, V(L+1, K+1), LDV, ONE, T(1, L+1), LDT) * * At this point, we have that T_{1,2} = V_1*V_2' * All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} @@ -438,13 +526,106 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{1,2} = -T_{1,1}*T_{1,2} * CALL ZTRMM('Left', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) + $ K-L, NEG_ONE, T, LDT, T(1, L+1), LDT) * * T_{1,2} = T_{1,2}*T_{2,2} * CALL ZTRMM('Right', 'Upper', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T(L+1, L+1), LDT, T(1, L+1), LDT) + $ K-L, ONE, T(L+1,L+1), LDT, T(1, L+1), LDT) + ELSE IF(LQT) THEN +* +* Break V apart into 6 components +* +* V = |----------------------| +* |V_{1,1} V_{1,2} V{1,3}| +* |0 V_{2,2} V{2,3}| +* |----------------------| +* +* V_{1,1}\in\C^{l,l} unit upper triangular +* V_{1,2}\in\C^{l,k-l} rectangular +* V_{1,3}\in\C^{l,n-k} rectangular +* +* V_{2,2}\in\C^{k-l,k-l} unit upper triangular +* V_{2,3}\in\C^{k-l,n-k} rectangular +* +* Where l = floor(k/2) +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} 0 | +* |T_{2,1} T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\C^{l, l} lower triangular +* T_{2,2}\in\C^{k-l, k-l} lower triangular +* T_{2,1}\in\C^{k-l, l} rectangular +* +* Then, consider the product: +* +* (I - V_1'*T_{1,1}*V_1)*(I - V_2'*T_{2,2}*V_2) +* = I - V_1'*T_{1,1}*V_1 - V_2'*T_{2,2}*V_2 + V_1'*T_{1,1}*V_1*V_2'*T_{2,2}*V_2 +* +* Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{2,1} +* +* Compute T_{1,1} recursively +* + CALL ZLARFT(DIRECT, STOREV, N, L, V, LDV, TAU, T, LDT) +* +* Compute T_{2,2} recursively +* + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V(L+1,L+1), LDV, + $ TAU(L+1), T(L+1,L+1), LDT) + +* +* Compute T_{2,1} +* T_{2,1} = V_{1,2}' +* + DO I = 1, K-L + DO J = 1, L + T(L+I,J) = CONJG(V(J,L+I)) + END DO + END DO +* +* T_{2,1} = V_{2,2}*T_{2,1} +* + CALL ZTRMM('Left', 'Upper', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(L+1,L+1), LDV, T(L+1,1), LDT) +* +* T_{2,1} = V_{2,3}*V_{1,3}' + T_{2,1} +* Note: We assume K <= N, and GEMM will do nothing if N=K +* + CALL ZGEMM('No Transpose', 'Conjugate', K-L, L, N-K, ONE, + $ V(L+1,K+1), LDV, V(1, K+1), LDV, ONE, T(L+1,1), LDT) +* +* At this point, we have that T_{2,1} = V_2*V_1' +* All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} +* respectively. +* +* T_{2,1} = -T_{2,2}*T_{2,1} +* + CALL ZTRMM('Left', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, NEG_ONE, T(L+1,L+1), LDT, T(L+1,1), LDT) +* +* T_{2,1} = T_{2,1}*T_{1,1} +* + CALL ZTRMM('Right', 'Lower', 'No Transpose', 'Non-unit', + $ K-L, L, ONE, T, LDT, T(L+1,1), LDT) ELSE IF(QL) THEN * * Break V apart into 6 components @@ -457,18 +638,18 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * V_{1,1}\in\C^{n-k,k-l} rectangular * V_{2,1}\in\C^{k-l,k-l} unit upper triangular -* +* * V_{1,2}\in\C^{n-k,l} rectangular * V_{2,2}\in\C^{k-l,l} rectangular * V_{3,2}\in\C^{l,l} unit upper triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -479,17 +660,17 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* +* * (I - V_2*T_{2,2}*V_2')*(I - V_1*T_{1,1}*V_1') * = I - V_2*T_{2,2}*V_2' - V_1*T_{1,1}*V_1' + V_2*T_{2,2}*V_2'*V_1*T_{1,1}*V_1' -* +* * Define T_{2,1} = -T_{2,2}*V_2'*V_1*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |-------| * |V_1 V_2| * |-------| -* +* * So, our product is equivalent to the matrix product * I - V*T*V' * This means, we can compute T_{1,1} and T_{2,2}, then use this information @@ -497,34 +678,34 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * Compute T_{1,1} recursively * - CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * CALL ZLARFT(DIRECT, STOREV, N, L, V(1, K-L+1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2}' * DO J = 1, K-L DO I = 1, L - T(K-L+I, J) = CONJG(V(N-K+J, K-L+I)) + T(K-L+I,J) = CONJG(V(N-K+J, K-L+I)) END DO END DO * * T_{2,1} = T_{2,1}*V_{2,1} * CALL ZTRMM('Right', 'Upper', 'No transpose', 'Unit', L, - $ K-L, ONE, V(N-K+1, 1), LDV, T(K-L+1, 1), LDT) + $ K-L, ONE, V(N-K+1,1), LDV, T(K-L+1,1), LDT) * * T_{2,1} = V_{2,2}'*V_{2,1} + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * CALL ZGEMM('Conjugate', 'No transpose', L, K-L, N-K, ONE, - $ V(1, K-L+1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + $ V(1,K-L+1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2'*V_1 * All that is left is to pre and post multiply by -T_{2,2} and T_{1,1} @@ -533,17 +714,13 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * T_{2,1} = -T_{2,2}*T_{2,1} * CALL ZTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * CALL ZTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) - ELSE -* -* Else means RQ case -* + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQ) THEN * Break V apart into 6 components * * V = |-----------------------| @@ -558,13 +735,13 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * V_{2,2}\in\C^{l,k-l} rectangular * V_{2,3}\in\C^{l,l} unit lower triangular * -* We will construct the T matrix +* We will construct the T matrix * T = |---------------| * |T_{1,1} 0 | * |T_{2,1} T_{2,2}| * |---------------| * -* T is the triangular factor obtained from block reflectors. +* T is the triangular factor obtained from block reflectors. * To motivate the structure, assume we have already computed T_{1,1} * and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 * @@ -575,51 +752,51 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * Where l = floor(k/2) * * Then, consider the product: -* -* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) -* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 -* +* +* (I - V_2'*T_{2,2}'*V_2)*(I - V_1'*T_{1,1}'*V_1) +* = I - V_2'*T_{2,2}'*V_2 - V_1'*T_{1,1}'*V_1 + V_2'*T_{2,2}'*V_2*V_1'*T_{1,1}'*V_1 +* * Define T_{2,1} = -T_{2,2}*V_2*V_1'*T_{1,1} -* -* Then, we can define the matrix V as +* +* Then, we can define the matrix V as * V = |---| * |V_1| * |V_2| * |---| -* +* * So, our product is equivalent to the matrix product -* I - V'*T*V +* I - V'*T'*V * This means, we can compute T_{1,1} and T_{2,2}, then use this information * to compute T_{2,1} * * Compute T_{1,1} recursively * - CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, LDT) + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) * * Compute T_{2,2} recursively * - CALL ZLARFT(DIRECT, STOREV, N, L, V(K-L+1, 1), LDV, - $ TAU(K-L+1), T(K-L+1, K-L+1), LDT) + CALL ZLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) * * Compute T_{2,1} * T_{2,1} = V_{2,2} * - CALL ZLACPY('All', L, K-L, V(K-L+1, N-K+1), LDV, - $ T(K-L+1, 1), LDT) + CALL ZLACPY('All', L, K-L, V(K-L+1,N-K+1), LDV, T(K-L+1,1), + $ LDT) * * T_{2,1} = T_{2,1}*V_{1,2}' * CALL ZTRMM('Right', 'Lower', 'Conjugate', 'Unit', L, K-L, - $ ONE, V(1, N-K+1), LDV, T(K-L+1, 1), LDT) + $ ONE, V(1, N-K+1), LDV, T(K-L+1,1), LDT) * -* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} +* T_{2,1} = V_{2,1}*V_{1,1}' + T_{2,1} * Note: We assume K <= N, and GEMM will do nothing if N=K * - CALL ZGEMM('No transpose', 'Conjugate', L, K-L, N-K, ONE, - $ V(K-L+1, 1), LDV, V, LDV, ONE, T(K-L+1, 1), - $ LDT) + CALL ZGEMM('No transpose', 'Conjugate', L, K-L, N-K, ONE, + $ V(K-L+1,1), LDV, V, LDV, ONE, T(K-L+1,1), LDT) * * At this point, we have that T_{2,1} = V_2*V_1' @@ -628,14 +805,104 @@ RECURSIVE SUBROUTINE ZLARFT( DIRECT, STOREV, N, K, V, LDV, * * T_{2,1} = -T_{2,2}*T_{2,1} * - CALL ZTRMM('Left', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, NEG_ONE, T(K-L+1, K-L+1), LDT, - $ T(K-L+1, 1), LDT) + CALL ZTRMM('Left', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, NEG_ONE, T(K-L+1,K-L+1), LDT, T(K-L+1,1), LDT) * * T_{2,1} = T_{2,1}*T_{1,1} * - CALL ZTRMM('Right', 'Lower', 'No tranpose', 'Non-unit', L, - $ K-L, ONE, T, LDT, T(K-L+1, 1), LDT) + CALL ZTRMM('Right', 'Lower', 'No transpose', 'Non-unit', L, + $ K-L, ONE, T, LDT, T(K-L+1,1), LDT) + ELSE IF(RQT) THEN +* Break V apart into 6 components +* +* V = |-----------------------| +* |V_{1,1} V_{1,2} 0 | +* |V_{2,1} V_{2,2} V_{2,3}| +* |-----------------------| +* +* V_{1,1}\in\C^{k-l,n-k} rectangular +* V_{1,2}\in\C^{k-l,k-l} unit lower triangular +* +* V_{2,1}\in\C^{l,n-k} rectangular +* V_{2,2}\in\C^{l,k-l} rectangular +* V_{2,3}\in\C^{l,l} unit lower triangular +* +* We will construct the T matrix +* T = |---------------| +* |T_{1,1} T_{1,2}| +* | 0 T_{2,2}| +* |---------------| +* +* T is the triangular factor obtained from block reflectors. +* To motivate the structure, assume we have already computed T_{1,1} +* and T_{2,2}. Then collect the associated reflectors in V_1 and V_2 +* +* T_{1,1}\in\C^{k-l, k-l} non-unit upper triangular +* T_{2,2}\in\C^{l, l} non-unit upper triangular +* T_{1,2}\in\C^{k-l, l} rectangular +* +* Where l = floor(k/2) +* +* Then, consider the product: +* +* (I - V_2'*T_{2,2}*V_2)*(I - V_1'*T_{1,1}*V_1) +* = I - V_2'*T_{2,2}*V_2 - V_1'*T_{1,1}*V_1 + V_2'*T_{2,2}*V_2*V_1'*T_{1,1}*V_1 +* +* Define T_{1,2} = -T_{1,1}*V_1*V_2'*T_{2,2} +* +* Then, we can define the matrix V as +* V = |---| +* |V_1| +* |V_2| +* |---| +* +* So, our product is equivalent to the matrix product +* I - V'*T*V +* This means, we can compute T_{1,1} and T_{2,2}, then use this information +* to compute T_{1,2} +* +* Compute T_{1,1} recursively +* + CALL ZLARFT(DIRECT, STOREV, N-L, K-L, V, LDV, TAU, T, + $ LDT) +* +* Compute T_{2,2} recursively +* + CALL ZLARFT(DIRECT, STOREV, N, L, V(K-L+1,1), LDV, + $ TAU(K-L+1), T(K-L+1,K-L+1), LDT) +* +* Compute T_{1,2} +* T_{1,2} = V_{2,2}' +* + DO I = 1, K-L + DO J = 1, L + T(I,K-L+J) = CONJG(V(K-L+J, N-K+I)) + END DO + END DO +* +* T_{1,2} = V_{1,2}T_{1,2} +* + CALL ZTRMM('Left', 'Lower', 'No Transpose', 'Unit', K-L, L, + $ ONE, V(1,N-K+1), LDV, T(1,K-L+1), LDT) +* +* T_{1,2} = V_{1,1}V_{2,1}' + T_{1,2} +* + CALL ZGEMM('No Transpose', 'Conjugate', K-L, L, N-K, ONE, V, + $ LDV, V(K-L+1,1), LDV, ONE, T(1, K-L+1), LDT) +* +* At this point, we have that T_{1,2} = V_1*V_2' +* All that is left is to pre and post multiply by -T_{1,1} and T_{2,2} +* respectively. +* +* T_{1,2} = -T_{1,1}*T_{1,2} +* + CALL ZTRMM('Left', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, NEG_ONE, T, LDT, T(1, K-L+1), LDT) +* +* T_{1,2} = T_{1,2}*T_{2,2} +* + CALL ZTRMM('Right', 'Upper', 'No Transpose', 'Non-Unit', + $ K-L, L, ONE, T(K-L+1,K-L+1), LDT, T(1, K-L+1), LDT) END IF END SUBROUTINE diff --git a/SRC/zlarft_lvl2.f b/SRC/zlarft_lvl2.f index 808c7fdb2..b16350b61 100644 --- a/SRC/zlarft_lvl2.f +++ b/SRC/zlarft_lvl2.f @@ -1,18 +1,10 @@ -*> \brief \b ZLARFT_LVL2: Level 2 BLAS version for terminating case of ZLARFT. +*> \brief \b ZLARFT_LVL2 forms the triangular factor T of a block reflector H = I - vtvH * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> Download ZLARFT_LVL2 + dependencies -*> -*> [TGZ] -*> -*> [ZIP] -*> -*> [TXT] -* * Definition: * =========== * @@ -49,6 +41,9 @@ *> H(i) is stored in the i-th row of the array V, and *> *> H = I - V**H * T * V +*> +*> If DIRECT or STOREV = 'T', see Further Details for the shape of T +*> *> \endverbatim * * Arguments: @@ -61,6 +56,8 @@ *> multiplied to form the block reflector: *> = 'F': H = H(1) H(2) . . . H(k) (Forward) *> = 'B': H = H(k) . . . H(2) H(1) (Backward) +*> = 'T' (with STOREV='R'): H = H(k) . . . H(2) H(1) (backward) +*> but we return the T matrix that is already (conjugate) transposed *> \endverbatim *> *> \param[in] STOREV @@ -68,8 +65,10 @@ *> STOREV is CHARACTER*1 *> Specifies how the vectors which define the elementary *> reflectors are stored (see also Further Details): -*> = 'C': columnwise -*> = 'R': rowwise +*> = 'C': column-wise +*> = 'R': row-wise +*> = 'T': (With DIRECT='F') Row-wise, but we return the T +*> matrix that is already (conjugate) transposed. *> \endverbatim *> *> \param[in] N @@ -155,6 +154,29 @@ *> ( 1 v2 v3 ) ( v3 v3 v3 v3 1 ) *> ( 1 v3 ) *> ( 1 ) +*> +*> In addition, the shape of T is determined by these same flags as +*> in the below table. +*> 'U' denotes upper triangular +*> 'L' denotes lower triangular +*> 'X' denotes no current implementation +*> We also provide the logical variable that represents the case +*> in the code if it is implemented +*> +*> |-----------------------------------------------------------| +*> | | DIRECT = 'F' | DIRECT = 'B' | DIRECT = 'T' | +*> |--------------+--------------+--------------+--------------| +*> | STOREV = 'C' | U (QR) | L (QL) | X | +*> | STOREV = 'R' | U (LQ) | L (RQ) | U (RQT) | +*> | STOREV = 'T' | L (LQT) | X | X | +*> |-----------------------------------------------------------| +*> +*> Finally, the relationship between the (conjugate) transposed T matrices +*> are as follows: (Note that T_{FC} denotes the T associated with calling +*> this routine with DIRECT = 'F' and STOREV = 'C') +*> +*> T_{TR} = (T_{BR})**H +*> T_{FT} = (T_{FR})**H *> \endverbatim *> * ===================================================================== @@ -170,158 +192,568 @@ SUBROUTINE ZLARFT_LVL2( DIRECT, STOREV, N, K, V, LDV, TAU, INTEGER K, LDT, LDV, N * .. * .. Array Arguments .. +* COMPLEX*16 T( LDT, * ), TAU( * ), V( LDV, * ) * .. * -* ===================================================================== -* * .. Parameters .. - COMPLEX*16 ONE, ZERO - PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ), - $ ZERO = ( 0.0D+0, 0.0D+0 ) ) -* .. +* + COMPLEX*16 ONE + PARAMETER(ONE=(1.0D+0,0.0D+0)) +* * .. Local Scalars .. - INTEGER I, J, PREVLASTV, LASTV -* .. +* + INTEGER I,J,KMI,NMI,INFO + LOGICAL QR, LQ, QL, RQ, LQT, RQT, + $ DIRF, DIRB, DIRT, + $ STOREC, STORER, STORET +* * .. External Subroutines .. - EXTERNAL ZGEMV, ZTRMV, ZGEMM -* .. -* .. External Functions .. - LOGICAL LSAME - EXTERNAL LSAME +* + EXTERNAL ZTRMV,ZGEMV,ZGEMM,XERBLA +* +* .. External Functions.. +* + LOGICAL LSAME + EXTERNAL LSAME +* +* .. Intrinsic Functions.. +* + INTRINSIC CONJG * .. * .. Executable Statements .. * +* Convert our character flags to logical flags for later +* + DIRF = LSAME(DIRECT,'F') + DIRB = LSAME(DIRECT,'B') + DIRT = LSAME(DIRECT,'T') + STOREC = LSAME(STOREV,'C') + STORER = LSAME(STOREV,'R') + STORET = LSAME(STOREV,'T') +* +* Error handling for our character flags +* + INFO = 0 + IF( .NOT.(DIRF.OR.DIRB.OR.DIRT) ) THEN +* +* DIRECT holds an illegal value +* + INFO = 1 + ELSE IF( .NOT.(STOREC.OR.STORER.OR.STORET) ) THEN +* +* STOREV holds an illegal value +* + INFO = 2 + ELSE IF( DIRB.AND.STORET ) THEN +* +* This case is purposefully not implemented, but any other value for +* STOREV is valid, so we report STOREV as the invalid input +* + INFO = 2 + ELSE IF( DIRT.AND.STOREC ) THEN +* +* This case is purposefully not implemented, but any other value for +* DIRECT is valid, so we report DIRECT as the invalid input +* + INFO = 1 + ELSE IF( DIRT.AND.STORET ) THEN +* +* This case is purposefully not implemented, and is ambiguous what +* the user wants to do, so we arbitrarily say DIRECT is the incorrect +* character flag. +* + INFO = 1 + END IF + + IF( INFO.NE.0 ) THEN + CALL XERBLA('ZLARFT_LVL2', INFO) + RETURN + END IF +* * Quick return if possible * - IF( N.EQ.0 ) - $ RETURN -* - IF( LSAME( DIRECT, 'F' ) ) THEN - PREVLASTV = N - DO I = 1, K - PREVLASTV = MAX( PREVLASTV, I ) - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = 1, I - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * CONJG( V( I , J ) ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**H * V(i:j,i) -* - CALL ZGEMV( 'Conjugate transpose', J-I, I-1, - $ -TAU( I ), V( I+1, 1 ), LDV, - $ V( I+1, I ), 1, ONE, T( 1, I ), 1 ) - ELSE -* Skip any trailing zeros. - DO LASTV = N, I+1, -1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = 1, I-1 - T( J, I ) = -TAU( I ) * V( J , I ) - END DO - J = MIN( LASTV, PREVLASTV ) -* -* T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**H -* - CALL ZGEMM( 'N', 'C', I-1, 1, J-I, -TAU( I ), - $ V( 1, I+1 ), LDV, V( I, I+1 ), LDV, - $ ONE, T( 1, I ), LDT ) - END IF -* -* T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i) -* - CALL ZTRMV( 'Upper', 'No transpose', 'Non-unit', I-1, - $ T, - $ LDT, T( 1, I ), 1 ) - T( I, I ) = TAU( I ) - IF( I.GT.1 ) THEN - PREVLASTV = MAX( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF + IF(N.EQ.0.OR.K.EQ.0) THEN + RETURN + END IF +* +* Now we determine what factorization our flags are associated with +* +* QR happens when we have forward direction in column storage +* + QR = DIRF.AND.STOREC +* +* LQT happens when we have forward direction in row storage and want to compute the transpose of +* the T we would normally compute +* + LQT = DIRF.AND.STORET +* +* LQ happens when we have forward direction in row storage and want to compute the T we would +* normally compute +* + LQ = DIRF.AND.STORER +* +* QL happens when we have backward direction in column storage +* + QL = DIRB.AND.STOREC +* +* RQT happens when we have backward direction in row storage and want to compute the transpose +* of the T we would normally compute +* + RQT = DIRT.AND.STORER +* +* RQ happens when we have backward direction in row storage and want to compute the T that we +* would normally compute +* + RQ = DIRB.AND.STORER + IF( N.EQ.1.OR.K.EQ.1) THEN + IF( LQT.OR.RQT ) THEN + T(1,1) = CONJG(TAU(1)) + ELSE + T(1,1) = TAU(1) + END IF + RETURN + END IF + IF (QR) THEN +* +* Break V into 9 components +* +* V = |-----------------------| +* |V_{1,1} 0 0 | i-1 +* |V_{2,1} V_{2,2} 0 | 1 +* |V_{3,1} V_{3,2} V_{3,3}| n-i +* |-----------------------| +* i-1 1 k-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1}\\V_{2,1}\\V_{3,1}]' +* [0\\V_{2,2}\\V_{3,2}]T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{2,1}' + V_{3,1}'V_{3,2})T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{2,1}'T_{2,2} = -\tau_{i}V_{2,1}' +* T_{1,2} = -\tau_{i}V_{3,2}' V_{3,1} + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + + DO I = 2, K +* +* T_{1,2} = -V_{2,1}'V_{2,2}T_{2,2} = -\tau_i V_{2,1}' +* We must do this at copy time as otherwise gemv will do nothing +* on the last column when n=k, but we neet to make sure we are +* scaled by this value +* + DO J = 1, I-1 + T(J,I) = -TAU(I)*CONJG(V(I,J)) + END DO + +* +* T_{1,2} = -V_{3,1}'V_{3,2}T_{2,2} + T_{1,2} +* = -\tau{i} V_{3,2}'V_{3,1} + T_{1,2} +* + CALL ZGEMV('Conjugate Transpose', N-I, I-1, -TAU(I), + $ V(I+1,1), LDV, V(I+1,I), 1, ONE, T(1, I), 1) + + +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL ZTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) + +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - ELSE - PREVLASTV = 1 - DO I = K, 1, -1 - IF( TAU( I ).EQ.ZERO ) THEN -* -* H(i) = I -* - DO J = I, K - T( J, I ) = ZERO - END DO - ELSE -* -* general case -* - IF( I.LT.K ) THEN - IF( LSAME( STOREV, 'C' ) ) THEN -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( LASTV, I ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * CONJG( V( N-K+I , J ) ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**H * V(j:n-k+i,i) -* - CALL ZGEMV( 'Conjugate transpose', N-K+I-J, K-I, - $ -TAU( I ), V( J, I+1 ), LDV, V( J, I ), - $ 1, ONE, T( I+1, I ), 1 ) - ELSE -* Skip any leading zeros. - DO LASTV = 1, I-1 - IF( V( I, LASTV ).NE.ZERO ) EXIT - END DO - DO J = I+1, K - T( J, I ) = -TAU( I ) * V( J, N-K+I ) - END DO - J = MAX( LASTV, PREVLASTV ) -* -* T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**H -* - CALL ZGEMM( 'N', 'C', K-I, 1, N-K+I-J, - $ -TAU( I ), - $ V( I+1, J ), LDV, V( I, J ), LDV, - $ ONE, T( I+1, I ), LDT ) - END IF -* -* T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i) -* - CALL ZTRMV( 'Lower', 'No transpose', 'Non-unit', - $ K-I, - $ T( I+1, I+1 ), LDT, T( I+1, I ), 1 ) - IF( I.GT.1 ) THEN - PREVLASTV = MIN( PREVLASTV, LASTV ) - ELSE - PREVLASTV = LASTV - END IF - END IF - T( I, I ) = TAU( I ) - END IF + ELSE IF (LQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | i-1 +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{1,2} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,3}. On +* each iteration, T_{1:3,3} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{1,2} = -T_{1,1}[V_{1,1} V_{1,2} V_{1,3}][ 0 V_{2,2} V_{2,3} ]'T_{2,2} +* +* T_{1,2} = -T_{1,1}(V_{1,2} + V_{1,3}V_{2,3}')T_{2,2} +* +* This means we will do the following +* +* T_{1,2} = -V_{1,2}T_{2,2} = -\tau_{i}V_{1,2} +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* T_{1,2} = T_{1,1}T_{1,2} +* T_{2,2} = \tau{i} +* + T(1,1) = TAU(1) + + DO I = 2, K +* +* T_{1,2} = -\tau_{i}V_{1,2} +* + DO J = 1, I-1 + T(J, I) = -TAU(I)*V(J, I) + END DO + +* +* T_{1,2} = -\tau_{i}V_{1,3}V_{2,3}' + T_{1,2} +* + CALL ZGEMM('No Transpose', 'Conjugate Transpose', I-1, + $ 1, N-I, -TAU(I), V(1,I+1), LDV, V(I, I+1), LDV, ONE, + $ T(1, I), LDT) +* +* T_{1,2} = T_{1,1}T_{1,2} +* + CALL ZTRMV('Upper', 'No Transpose', 'Non-unit', I-1, + $ T, LDT, T(1,I), 1) + +* +* T_{2,2} = \tau{i} +* + T(I,I) = TAU(I) END DO - END IF - RETURN + ELSE IF (LQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | i-1 +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | k-i +* |-------------------------| +* i-1 1 n-i +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{3,3} and repeat until we +* reach the end. On each iteration V_{3,3} is not referenced +* +* We will construct T one column at a time from left to right +* after initializing T(1,1) = TAU(1) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | i-1 +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | k-i +* |-------------------------| +* i-1 1 k-i +* +* Similarly as above, we will construct T_{2,1} and T_{2,2} at +* each iteration i = 2, \dots k, and then grow into T_{3,1:3}. On +* each iteration, T_{3,1:3} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We now get +* +* T_{2,1} = -T_{2,2}[0 V_{2,2} V_{2,3}][V_{1,1} V_{1,2} V_{1,3}]'T_{1,1} +* +* T_{2,1} = -T_{2,2}(V_{1,2}' + V_{2,3}V_{1,3}')T_{1,1} +* +* This means we will do the following +* +* T_{2,1} = -T_{2,2}V_{1,2}' = -\tau_{i}V_{1,2}' +* T_{2,1} = -\tau_{i}V_{1,3}V_{2,3}' + T_{2,1} +* T_{2,1} = T_{1,1}'T_{2,1} +* T_{2,2} = \tau{i} +* + T(1,1) = CONJG(TAU(1)) + + DO I = 2, K +* +* T_{2,1} = -\tau_{i}V_{1,2}' +* + DO J = 1, I-1 + T(I,J) = -CONJG(TAU(I)*V(J,I)) + END DO +* +* T_{2,1} = -\tau_{i}V_{2,3}V_{1,3}' + T_{2,1} +* + CALL ZGEMM('No Transpose', 'Conjugate Transpose', 1, + $ I-1, N-I, -CONJG(TAU(I)), V(I,I+1), LDV, V(1, I+1), + $ LDV, ONE, T(I, 1), LDT) +* +* T_{2,1} = T_{1,1}'T_{2,1} +* + CALL ZTRMV('Lower', 'Transpose', 'Non-unit', + $ I-1, T, LDT, T(I,1), LDT) + + T(I,I) = CONJG(TAU(I)) + END DO + ELSE IF (QL) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} V_{1,2} V_{1,3} | n-i +* | 0 V_{2,2} V_{2,3} | 1 +* | 0 0 V_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit upper triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{1,3}\\V_{2,3}\\V_{3,3}]' +* [V_{1,2}\\V_{2,2}\\0]T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{1,3}'V_{1,2} + V_{2,3}')T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -T_{3,3}V_{3,2}' = -\tau_{k-i+1}V_{3,2}' +* T_{3,2} = -\tau_{k-i+1}V_{1,3}'V_{1,2} + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau(k-i+1) +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau(k-i+1)V_{2,3}' +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*CONJG(V(NMI, KMI + J)) + END DO +* +* T_{3,2} = -\tau(k-i+1)V_{1,3}'V_{1,2} + T_{3,2} +* + CALL ZGEMV('Conjugate Transpose', N-I, I-1, -TAU(KMI), + $ V(1, KMI + 1), LDV, V(1, KMI), 1, ONE, + $ T(KMI+1, KMI), 1) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL ZTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI + 1, KMI + 1), LDT, T(KMI + 1, KMI), 1) + + END DO + ELSE IF (RQ) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} 0 0 | k-i +* | T_{2,1} T_{2,2} 0 | 1 +* | T_{3,1} T_{3,2} T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{3,2} at +* each iteration i = 2, \dots k, and then grow into T_{1:3,1}. On +* each iteration, T_{1:3,1} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{3,2} = -T_{3,3}[V_{3,1} V_{3,2} V_{3,3}][V_{2,1} V_{2,2} 0]'T_{2,2} +* +* T_{3,2} = -T_{3,3}(V_{3,1}V_{2,1}' + V_{3,2})T_{2,2} +* +* Thus, we will compute +* +* T_{2,2} = \tau_{k-i+1} +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* T_{3,2} = T_{3,3}T_{3,2} +* + T(K,K) = TAU(K) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = TAU(KMI) +* +* T_{3,2} = -\tau_{k-i+1}V_{3,2} +* + DO J = 1, I-1 + T(KMI + J, KMI) = -TAU(KMI)*V(KMI + J, NMI) + END DO +* +* T_{3,2} = -\tau_{k-i+1}V_{3,1}V_{2,1}' + T_{3,2} +* + CALL ZGEMM('No Transpose', 'Conjugate Transpose', I-1, + $ 1, N-I, -TAU(KMI), V(KMI+1, 1), LDV, V(KMI, 1), LDV, + $ ONE, T(KMI+1, KMI), LDT) +* +* T_{3,2} = T_{3,3}T_{3,2} +* + CALL ZTRMV('Lower', 'No Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI+1, KMI), 1) + END DO + ELSE IF (RQT) THEN +* +* Break V into 9 components +* +* V = |-------------------------| +* | V_{1,1} 0 0 | k-i +* | V_{2,1} V_{2,2} 0 | 1 +* | V_{3,1} V_{3,2} V_{3,3} | i-1 +* |-------------------------| +* n-i 1 i-1 +* +* V_{1,1}, V_{2,2} and V_{3,3} are unit lower triangular +* +* This is how we are going to view the matrix V at each step +* i=2,\dots,k, then we grow into V_{1,1} and repeat until we +* reach the end. On each iteration V_{1,1} is not referenced +* +* We will construct T one column at a time from right to left +* after initializing T(K,K) = TAU(K) +* +* T = |-------------------------| +* | T_{1,1} T_{1,2} T_{1,3} | k-i +* | 0 T_{2,2} T_{2,3} | 1 +* | 0 0 T_{3,3} | i-1 +* |-------------------------| +* k-i 1 i-1 +* +* T_{1,1}, T_{2,2}, and T_{3,3} are non-unit lower triangular +* +* Similarly as above, we will construct T_{2,2} and T_{2,3} at +* each iteration i = 2, \dots k, and then grow into T_{1,1:3}. On +* each iteration, T_{1,1:3} are not referenced. See zlarft.f +* for details on how these formulae were constructed. +* +* We get that +* +* T_{2,3} = -T_{2,2}[V_{2,1} V_{2,2} 0][V_{3,1} V_{3,2} V_{3,3}]'T_{3,3} +* +* T_{3,2} = -T_{2,2}(V_{2,1}V_{3,1}' + V_{3,2}')T_{3,3} * -* End of ZLARFT_LVL2 +* Thus, we will compute * - END +* T_{2,2} = \tau_{k-i+1} +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* T_{2,3} = -\tau_{k-i+1}V_{2,1}V_{3,1}' + T_{2,3} +* T_{2,3} = T_{3,3}'T_{2,3} +* + T(K,K) = CONJG(TAU(K)) + DO I = 2, K + KMI = K-I+1 + NMI = N-I+1 +* +* T_{2,2} = \tau_{k-i+1} +* + T(KMI,KMI) = CONJG(TAU(KMI)) +* +* T_{2,3} = -\tau_{k-i+1}V_{3,2}' +* + DO J = 1, I-1 + T(KMI, KMI + J) = -CONJG(TAU(KMI)*V(KMI + J, NMI)) + END DO +* +* T_{2,3} = -\tau_{k-i+1}V_{2,1}V_{3,1}' + T_{2,3} +* + CALL ZGEMM('No Transpose', 'Conjugate Transpose', 1, + $ I-1, N-I, -CONJG(TAU(KMI)), V(KMI, 1), LDV, + $ V(KMI+1,1), LDV, ONE, T(KMI, KMI+1), LDT) +* +* T_{2,3} = T_{3,3}'T_{2,3} +* + CALL ZTRMV('Upper', 'Transpose', 'Non-unit', I-1, + $ T(KMI+1, KMI+1), LDT, T(KMI, KMI+1), LDT) + END DO + END IF + END SUBROUTINE