SUBROUTINE MVMBWM( TRANS, N, M, X, LDX, Y, LDY ) * .. * .. Scalar Arguments .. INTEGER LDY, LDX, M, N, TRANS * .. * .. Array Arguments .. DOUBLE PRECISION Y( LDY, * ), X( LDX, * ) * .. * * Purpose * ======= * * Compute * * Y(:,1:M) = op(A)*X(:,1:M) * * where op(A) is A or A' (the transpose of A). The matrix A is a 2 by 2 * block matrix resulted from the finite difference discretization * of the Brusselator wave model. * * Arguments * ========= * * TRANS (input) INTEGER * If TRANS = 0, compute Y(:,1:M) = A*Q(:,1:M) * If TRANS = 1, compute Y(:,1:M) = A'*Q(:,1:M) * * N (input) INTEGER * The order of the matrix A. N has to be an even number. * * M (input) INTEGERS * The number of Q to multiply. * * X (input) DOUBLE PRECISION array, dimension ( LDX, M ) * contains the vectors X. * * LDX (input) INTEGER * The leading dimension of the array X, LDX >= max( 1, N ) * * Y (output) DOUBLE PRECISION array, dimension (LDY, M ) * contains the product of the matrix op(A) with X. * * LDY (input) INTEGER * The leading dimension of the array Y, LDY >= max( 1, N ) * * =================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, TWO PARAMETER ( ONE = 1.0D+0, TWO = 2.0D+0 ) DOUBLE PRECISION ALPHA, BETA, DELT1, DELT2, L PARAMETER ( ALPHA = 2.0D+0, BETA = 5.45D+0, $ DELT1 = 8.0D-3, DELT2 = 4.0D-3, $ L = 0.51302D+0 ) * * .. Local Variables .. INTEGER I, K, NHALF DOUBLE PRECISION HSTEP, TAU1, TAU2, DIAG1, DIAG2, ALPHA2 * * set up middle parameters * NHALF = N / 2 HSTEP = ONE / ( DBLE( NHALF ) + 1 ) TAU1 = DELT1 / (HSTEP*L)**2 TAU2 = DELT2 / (HSTEP*L)**2 DIAG1 = -TWO*TAU1 + BETA - ONE ALPHA2 = ALPHA**2 DIAG2 = -TWO*TAU2 - ALPHA2 * IF ( TRANS.EQ.0 ) THEN * * Compute Y(:,1*M) = A*X(:,1:M) * DO 30 K = 1, M * * the 1st row * Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K ) + $ ALPHA2*X( NHALF+1, K ) * * the 2nd row to (NHALF-1)th row * DO 10 I = 2, NHALF - 1 Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) + $ TAU1*X( I+1, K ) + ALPHA2*X( NHALF+I, K ) 10 CONTINUE * * the NHALFth row * Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) $ + ALPHA2*X( N,K ) * * the (NHALF+1)th row * Y( NHALF+1, K ) = -BETA*X( 1, K ) + DIAG2*X( NHALF+1, K ) $ + TAU2*X( NHALF+2, K ) * * the (NHALF+2)th row to (N-1)th row * DO 20 I = NHALF+2, N - 1 Y( I, K ) = -BETA*X( I-NHALF, K ) + TAU2*X( I-1, K ) + $ DIAG2*X( I, K ) + TAU2*X( I+1, K ) 20 CONTINUE * * the last (Nth) row * Y( N, K ) = -BETA*X( N-NHALF, K ) + TAU2*X( N-1,K ) + $ DIAG2*X( N, K ) * 30 CONTINUE * ELSE IF ( TRANS.EQ. 1 ) THEN * * Compute Y(:,1:M) = A'*X(:,1:M) * DO 60 K = 1, M * * the 1st row * Y( 1, K ) = DIAG1*X( 1, K ) + TAU1*X( 2, K ) $ - BETA*X( NHALF+1, K ) * * the 2nd row to (NHALF-1)th row * DO 40 I = 2, NHALF - 1 Y( I, K ) = TAU1*X( I-1, K ) + DIAG1*X( I, K ) + $ TAU1*X( I+1, K ) - BETA*X( NHALF+I, K ) 40 CONTINUE * * the NHALFth row * Y( NHALF, K ) = TAU1*X( NHALF-1, K ) + DIAG1*X( NHALF, K ) $ - BETA*X( N,K ) * * the (NHALF+1)th row * Y( NHALF+1, K ) = ALPHA2*X( 1, K ) + DIAG2*X( NHALF+1, K ) $ + TAU2*X( NHALF+2, K ) * * the (NHALF+2)th row to (N-1)th row * DO 50 I = NHALF+2, N - 1 Y( I, K ) = ALPHA2*X( I-NHALF, K ) + TAU2*X( I-1, K ) + $ DIAG2*X( I, K ) + TAU2*X( I+1, K ) 50 CONTINUE * * the last (Nth) row * Y( N, K ) = ALPHA2*X( N-NHALF, K ) + TAU2*X( N-1,K ) + $ DIAG2*X( N, K ) * 60 CONTINUE * END IF * RETURN * * End of MVMBWM * END