SUBROUTINE MVMMCD( 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 block * tridiagonal matrix resulted from the finite difference discretization * of a 2-D model convection diffusion operator, * * L[u] = - u_xx - u_yy + 2p_1 u_x + 2p_2 u_y - p_3 u. * * The constant parameters p1, p2 and p3 in the convection-diffusion * operator may be changed. * * The convection-diffusion operator is discretized on a K x K square * grid, and resulted in the matrix A of order N = K^2. * * Arguments * ========= * * TRANS (input) INTEGER * If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M) * If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M) * * N (input) INTEGER * The order of the matrix A. N has to be K*K for an integer K. * * M (input) INTEGERS * The number of columns of X to multiply. * * X (input) DOUBLE PRECISION array, dimension ( LDX, M ) * X contains the matrix (vectors) X. * * LDX (input) INTEGER * The leading dimension of array X, LDX >= max( 1, N ) * * Y (output) DOUBLE PRECISION array, dimension (LDY, M ) * contains the product of the matrix and vectors. * * LDY (input) INTEGER * The leading dimension of array Y, LDY >= max( 1, N ) * * =================================================================== * * .. Parameters .. DOUBLE PRECISION ONE, FOUR PARAMETER ( ONE = 1.0D+0, FOUR = 4.0D+0 ) DOUBLE PRECISION P1, P2, P3 PARAMETER ( P1 = 1.0D+0, P2 = 1.0D+0, P3 = 1.0D+0 ) * .. * .. Local Scalars .. INTEGER I, II, K, NS DOUBLE PRECISION BETA, GAMMA, HSTEP, SIGMA * .. * .. Intrinsic Functions .. INTRINSIC REAL, SQRT * .. * .. Executable Statements .. * NS = SQRT( REAL( N ) ) HSTEP = ONE / ( REAL( NS )+ONE ) SIGMA = P3*HSTEP*HSTEP GAMMA = P2*HSTEP BETA = P1*HSTEP * IF ( TRANS.EQ.0 ) THEN DO 70 K = 1, M * * Compute the vector Y(1:N,K) = A*X(1:N,K) * DO 60 I = 1, NS * * diagonal elements * DO 10 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = ( FOUR-SIGMA )*X( II, K ) 10 CONTINUE * * sub diagonal elements * DO 20 II = ( I-1 )*NS + 1 + 1, I*NS Y( II, K ) = Y( II, K ) + \$ ( -GAMMA-ONE )*X( II-1, K ) 20 CONTINUE * * super diagonal elements * DO 30 II = ( I-1 )*NS + 1, I*NS - 1 Y( II, K ) = Y( II, K ) + \$ ( GAMMA-ONE )*X( II+1, K ) 30 CONTINUE * * (NS+1)th sub diagonal * IF( I.GT.1 ) THEN DO 40 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = Y( II, K ) + \$ ( BETA-ONE )*X( II-NS, K ) 40 CONTINUE END IF * * (NS+1)th super diagonal * IF( I.LT.NS ) THEN DO 50 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = Y( II, K ) + \$ ( -BETA-ONE )*X( II+NS, K ) 50 CONTINUE END IF * 60 CONTINUE * 70 CONTINUE * ELSE IF ( TRANS.EQ.1 ) THEN * DO 140 K = 1, M * * Compute the vector Y(1:N,K) = trans(A)*X(1:N,K) * DO 130 I = 1, NS * DO 80 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = ( FOUR-SIGMA )*X( II, K ) 80 CONTINUE * DO 90 II = ( I-1 )*NS + 1, I*NS - 1 Y( II, K ) = Y( II, K ) + \$ ( -GAMMA-ONE )*X( II+1, K ) 90 CONTINUE * DO 100 II = ( I-1 )*NS + 1 + 1, I*NS Y( II, K ) = Y( II, K ) + \$ ( GAMMA-ONE )*X( II-1, K ) 100 CONTINUE * * (NS+1)th sub diagonal * IF( I.GT.1 ) THEN DO 110 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = Y( II, K ) + \$ (-BETA-ONE )*X( II-NS, K ) 110 CONTINUE END IF * * (NS+1)th super diagonal * IF( I.LT.NS ) THEN DO 120 II = ( I-1 )*NS + 1, I*NS Y( II, K ) = Y( II, K ) + \$ ( BETA-ONE )*X( II+NS, K ) 120 CONTINUE END IF * 130 CONTINUE * 140 CONTINUE * END IF * RETURN * * End of MVMMCD END