SUBROUTINE MVMRWK( N, M, X, LDX, Y, LDY ) * .. * .. Scalar Arguments .. INTEGER LDY, LDX, M, N * .. * .. Array Arguments .. DOUBLE PRECISION Y( LDY, * ), X( LDX, * ) * .. * * Purpose * ======= * * Compute * * Y(:,1:M) = A*X(:,1:M) * * The matrix is from a random walk example. The order of the random * walk test matrix is N = (K+1)*(K+2)/2, where K is the size of * triangular random walk grid. * * Arguments * ========= * * N (input) INTEGER * The order of the matrix A. * N should be equal to (K+1)*(K+2)/2 for an integer K > 0. * * M (input) INTEGERS * The number of columns of X to multiply. * * X (input) REAL array, dimension ( LDX, M ) * X contains the matrix X. * * LDX (input) INTEGER * The leading dimension of the array X, LDX >= max(1,N) * * Y (output) REAL array, dimension (LDY, M ) * contains the product of the matrix A with X. * * LDY (input) INTEGER * The leading dimension of the array y, LDY >= max(1,N) * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, HALF PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0 ) * .. * .. Local Scalars .. INTEGER I, J, K, KK, LK, NN DOUBLE PRECISION UPPER * .. * .. Intrinsic Functions .. INTRINSIC REAL, SQRT * .. * .. Statement Functions .. INTEGER IINDEX * .. * .. Statement Function definitions .. IINDEX( K, I, J ) = I*K - ( ( I-1 )*( I-2 ) ) / 2 + 1 + J + 1 * .. * .. Executable Statements .. * NN = ( -3+SQRT( REAL( 8*N+1 ) ) ) / 2 * * * Compute Y( :,1:M ) = A*X( :,1:M ) * DO 30 KK = 1, M * * Compute the vector Y(1:N,KK) = A*X(1:N,KK) * DO 20 I = 0, NN DO 10 J = 0, NN - I * K = IINDEX( NN, I, J ) UPPER = REAL( I+J ) / REAL( 2*NN ) Y( K, KK ) = ZERO * IF( I.EQ.0 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I, J+1 ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) LK = IINDEX( NN, I+1, J ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.LE.NN-1 ) THEN LK = IINDEX( NN, I, J-1 ) Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK ) LK = IINDEX( NN, I+1, J ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) LK = IINDEX( NN, I, J+1 ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) * ELSE IF( I.EQ.0 .AND. J.EQ.NN ) THEN LK = IINDEX( NN, I, J-1 ) Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK ) * ELSE IF( I.LE.NN-1 .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK ) LK = IINDEX( NN, I+1, J ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) LK = IINDEX( NN, I, J+1 ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) * ELSE IF( I.EQ.NN .AND. J.EQ.0 ) THEN LK = IINDEX( NN, I-1, J ) Y( K, KK ) = Y( K, KK ) + 2*UPPER*X( LK, KK ) * ELSE IF( ( I+J ).EQ.NN ) THEN LK = IINDEX( NN, I-1, J ) Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK ) LK = IINDEX( NN, I, J-1 ) Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK ) * ELSE LK = IINDEX( NN, I-1, J ) Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK ) LK = IINDEX( NN, I, J-1 ) Y( K, KK ) = Y( K, KK ) + UPPER*X( LK, KK ) LK = IINDEX( NN, I, J+1 ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) LK = IINDEX( NN, I+1, J ) Y( K, KK ) = Y( K, KK ) + \$ ( HALF-UPPER )*X( LK, KK ) * END IF 10 CONTINUE 20 CONTINUE * 30 CONTINUE * RETURN * * End of MVMRWK * END