SUBROUTINE ZBQRU_APPQ( SIDE, TRANS, M, N, K, V, LDV, TAU, C, LDC, \$ D, LDD, WORK, LWORK, INFO ) IMPLICIT NONE * * .. Scalar Arguments .. CHARACTER SIDE, TRANS INTEGER INFO, LDC, LDD, LDV, LWORK, K, M, N * .. * .. Array Arguments .. COMPLEX*16 C( LDC, * ), D( LDD, * ), V( LDV, * ), \$ TAU( * ), WORK( * ) * .. * * Purpose * ======= * * Applies the unitary factor Q (or its transpose) generated by the * routine ZBQRU. * If SIDE = 'L' then Q or Q' is applied from the left to the * (k+m) by n matrix [C]. * [D] * If SIDE = 'R' then Q or Q' is applied from the right to the * m by k+n matrix [C D]. * * * Arguments * ========= * * SIDE (input) CHARACTER*1 * = 'L': apply H or H' from the Left * = 'R': apply H or H' from the Right * * TRANS (input) CHARACTER*1 * = 'N': apply H (No transpose) * = 'C': apply H' (Conjugate transpose) * * M (input) INTEGER * The number of rows of the matrix D. * * N (input) INTEGER * The number of columns of the matrix D. * * K (input) INTEGER * The number of elementary reflectors whose product defines Q. * (Normally the order of the triangular matrix A in the call to * ZBQRU). * * V (input) COMPLEX*16 array, dimension (LDV,K) * The matrix V containing the Householder reflectors. * * LDV (input) INTEGER * The leading dimension of the array V. * If SIDE = 'L', LDV >= max(1,M); * if SIDE = 'R', LDV >= max(1,N). * * C (input/output) COMPLEX*16 array, dimension * (LDC,N) if SIDE = 'L' or (LDC,K) if SIDE = 'R'. * On entry, the k by n (if SIDE = 'L') or m by k * (if SIDE = 'R') matrix C. * On exit, C is updated by the matrix Q. * * LDC (input) INTEGER * The leading dimension of the array C. * If SIDE = 'L', LDC >= max(1,K); * if SIDE = 'R', LDC >= max(1,M). * * D (input/output) COMPLEX*16 array, dimension (LDD,N) * On entry, the m by n matrix D. * On exit, D is updated by the matrix Q. * * LDD (input) INTEGER * The leading dimension of the array D. LDD >= max(1,M); * * WORK (workspace/output) COMPLEX*16 array, dimension * (MAX(1,LWORK)) * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. * * LWORK (input) INTEGER * The dimension of the array WORK. * If SIDE = 'L', LDWORK >= max(1,N); * if SIDE = 'R', LDWORK >= max(1,M). * For optimum performance LWORK >= P*NB, where NB is * the optimal blocksize and P = M if SIDE = 'L' or P = N if * SIDE = 'R'. * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * INFO (output) INTEGER * = 0: successful exit * < 0: if INFO = -i, the i-th argument had an illegal value * If LWORK = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the WORK array, returns * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * * Further Details * =============== * * The matrix Q is represented as a product of elementary reflectors * * Q = H(1) H(2) . . . H(k). * * Each H(i) has the form * * H(i) = I - tau * v * v' * * where tau is a real scalar, and the nontrivial part of v is stored in * V(:,i), and tau in TAU(i). See also ZBQRU. * * ===================================================================== * * .. Parameters .. INTEGER NBMAX, LDT PARAMETER ( NBMAX = 64, LDT = NBMAX+1 ) * .. Local Scalars .. LOGICAL LBLOCK, LEFT, LQUERY, NOTRAN INTEGER I, J, JB, LWKOPT, MBLOCK, \$ NB, ROWV, UPD * .. Local Arrays .. COMPLEX*16 T( LDT, NBMAX ) * .. * .. External Subroutines .. EXTERNAL ZBQRU_LARF, ZBQRU_LARFB, ZBQRU_LARFT, \$ XERBLA * .. * .. External Functions .. INTEGER ILAENV LOGICAL LSAME EXTERNAL ILAENV, LSAME * .. * .. Intrinsic Functions .. INTRINSIC MAX, MIN * .. * .. Executable Statements .. * * * Test the input arguments * INFO = 0 LEFT = LSAME( SIDE, 'L' ) NOTRAN = LSAME( TRANS, 'N' ) LQUERY = ( LWORK.EQ.-1 ) * * ROWV is the number of rows in V. * UPD is the minimal length of WORK. * IF( LEFT ) THEN ROWV = M UPD = N ELSE ROWV = N UPD = M END IF * * Block parameters. * Tune these parameters to get optimal performance. * The default settings you find below are just a * heuristics and should be fairly OK. * NB = panel block size NB = ILAENV( 1, 'ZGEQRF', ' ', ROWV, UPD, -1, -1 ) NB = MIN( NB, NBMAX ) C PRINT*, 'NB = ', NB * MBLOCK = minimal ROWV for which the block algorithm is used MBLOCK = NB/2 LWKOPT = UPD*NB WORK( 1 ) = LWKOPT LQUERY = ( LWORK.EQ.-1 ) IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN INFO = -1 ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN INFO = -2 ELSE IF( M.LT.0 ) THEN INFO = -3 ELSE IF( N.LT.0 ) THEN INFO = -4 ELSE IF( K.LT.0 ) THEN INFO = -5 ELSE IF( LDV.LT.MAX( 1, ROWV ) ) THEN INFO = -7 ELSE IF( ( LEFT.AND.LDC.LT.MAX( 1, K ) ).OR. \$ ( ( .NOT.LEFT ).AND.LDC.LT.MAX( 1, M ) ) ) THEN INFO = -10 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN INFO = -12 ELSE IF( LWORK.LT.MAX( 1, UPD ) .AND. .NOT.LQUERY ) THEN INFO = -14 END IF IF( INFO.NE.0 ) THEN CALL XERBLA( 'ZBQRUA', -INFO ) RETURN ELSE IF( LQUERY ) THEN RETURN END IF * * Quick return if possible * IF( K.EQ.0 .OR. M.EQ.0 .OR. N.EQ.0 ) THEN WORK( 1 ) = 1 RETURN END IF LBLOCK = ( ROWV.GE.MBLOCK ).AND.( LWORK.GE.NB*UPD ) C PRINT*, 'LBLOCK = ', LBLOCK IF ( LEFT ) THEN IF ( NOTRAN ) THEN * * Apply Q from left. * DO 20 J = 1, K, NB JB = MIN( K-J+1, NB ) IF ( LBLOCK ) THEN CALL ZBQRU_LARFT( ROWV, JB, V(1,J), LDV, TAU(J), \$ T, LDT ) CALL ZBQRU_LARFB( 'Left', 'No Transpose', M, N, JB, \$ V(1,J), LDV, T, LDT, C(J,1), \$ LDC, D, LDD, WORK, UPD ) ELSE DO 10 I = J, J+JB-1 CALL ZBQRU_LARF( 'Left', ROWV, N, V( 1, I ), 1, \$ TAU( I ), C( I, 1 ), LDC, D, LDD, \$ WORK ) 10 CONTINUE END IF 20 CONTINUE ELSE * * Apply Q' from left. * DO 40 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB JB = MIN( K-J+1, NB ) IF ( LBLOCK ) THEN CALL ZBQRU_LARFT( ROWV, JB, V(1,J), LDV, TAU(J), \$ T, LDT ) CALL ZBQRU_LARFB( 'Left', 'Conjugate transpose', M, N, \$ JB, V(1,J), LDV, T, LDT, C(J,1), \$ LDC, D, LDD, WORK, UPD ) ELSE DO 30 I = J+JB-1, J, -1 CALL ZBQRU_LARF( 'Left', ROWV, N, V( 1, I ), 1, \$ DCONJG( TAU( I ) ), C( I, 1 ), \$ LDC, D, LDD, WORK ) 30 CONTINUE END IF 40 CONTINUE END IF ELSE IF ( NOTRAN ) THEN * * Apply Q from right. * DO 100 J = 1, K, NB JB = MIN( K-J+1, NB ) IF ( LBLOCK ) THEN CALL ZBQRU_LARFT( ROWV, JB, V(1,J), LDV, TAU(J), \$ T, LDT ) CALL ZBQRU_LARFB( 'Right', 'No Transpose', M, N, JB, \$ V(1,J), LDV, T, LDT, C(1,J), \$ LDC, D, LDD, WORK, UPD ) ELSE DO 90 I = J, J+JB-1 CALL ZBQRU_LARF( 'Right', M, ROWV, V( 1, I ), 1, \$ TAU( I ), C( 1, I ), 1, \$ D, LDD, WORK ) 90 CONTINUE END IF 100 CONTINUE ELSE DO 120 J = ( ( K-1 ) / NB )*NB + 1, 1, -NB JB = MIN( K-J+1, NB ) IF ( LBLOCK ) THEN CALL ZBQRU_LARFT( ROWV, JB, V(1,J), LDV, TAU(J), \$ T, LDT ) CALL ZBQRU_LARFB( 'Right', 'Conjugate Transpose', M, \$ N, JB, V(1,J), LDV, T, LDT, C(1,J), \$ LDC, D, LDD, WORK, UPD ) ELSE DO 110 I = J+JB-1, J, -1 CALL ZBQRU_LARF( 'Right', M, ROWV, V( 1, I ), 1, \$ DCONJG( TAU( I ) ), C( 1, I ), 1, \$ D, LDD, WORK ) 110 CONTINUE END IF 120 CONTINUE END IF END IF RETURN * * End of ZBQRU_APPQ * END