Actual source code: mvmisg.f
1: SUBROUTINE MVMISG( TRANS, N, M, X, LDX, Y, LDY )
2: * ..
3: * .. Scalar Arguments ..
4: INTEGER LDY, LDX, M, N, TRANS
5: * ..
6: * .. Array Arguments ..
7: DOUBLE PRECISION Y( LDY, * ), X( LDX, * )
8: * ..
9: *
10: * Purpose
11: * =======
12: *
13: * Compute
14: *
15: * Y(:,1:M) = op(A)*X(:,1:M)
16: *
17: * where op(A) is A or A' (the transpose of A). The A is the Ising
18: * matrix.
19: *
20: * Arguments
21: * =========
22: *
23: * TRANS (input) INTEGER
24: * If TRANS = 0, compute Y(:,1:M) = A*X(:,1:M)
25: * If TRANS = 1, compute Y(:,1:M) = A'*X(:,1:M)
26: *
27: * N (input) INTEGER
28: * The order of the matrix A. N has to be an even number.
29: *
30: * M (input) INTEGERS
31: * The number of columns of X to multiply.
32: *
33: * X (input) DOUBLE PRECISION array, dimension ( LDX, M )
34: * X contains the matrix (vectors) X.
35: *
36: * LDX (input) INTEGER
37: * The leading dimension of array X, LDX >= max( 1, N )
38: *
39: * Y (output) DOUBLE PRECISION array, dimension (LDX, M )
40: * contains the product of the matrix op(A) with X.
41: *
42: * LDY (input) INTEGER
43: * The leading dimension of array Y, LDY >= max( 1, N )
44: *
45: * ===================================================================
46: *
47: *
48: * .. PARAMETERS ..
49: DOUBLE PRECISION PI
50: PARAMETER ( PI = 3.141592653589793D+00 )
51: DOUBLE PRECISION ALPHA, BETA
52: PARAMETER ( ALPHA = PI/4, BETA = PI/4 )
53: *
54: * .. Local Variables ..
55: INTEGER I, K
56: DOUBLE PRECISION COSA, COSB, SINA, SINB, TEMP, TEMP1
57: *
58: * .. Intrinsic functions ..
59: INTRINSIC COS, SIN
60: *
61: COSA = COS( ALPHA )
62: SINA = SIN( ALPHA )
63: COSB = COS( BETA )
64: SINB = SIN( BETA )
65: *
66: IF ( TRANS.EQ.0 ) THEN
67: *
68: * Compute Y(:,1:M) = A*X(:,1:M)
70: DO 30 K = 1, M
71: *
72: Y( 1, K ) = COSB*X( 1, K ) - SINB*X( N, K )
73: DO 10 I = 2, N-1, 2
74: Y( I, K ) = COSB*X( I, K ) + SINB*X( I+1, K )
75: Y( I+1, K ) = -SINB*X( I, K ) + COSB*X( I+1, K )
76: 10 CONTINUE
77: Y( N, K ) = SINB*X( 1, K ) + COSB*X( N, K )
78: *
79: DO 20 I = 1, N, 2
80: TEMP = COSA*Y( I, K ) + SINA*Y( I+1, K )
81: Y( I+1, K ) = -SINA*Y( I, K ) + COSA*Y( I+1, K )
82: Y( I, K ) = TEMP
83: 20 CONTINUE
84: *
85: 30 CONTINUE
86: *
87: ELSE IF ( TRANS.EQ.1 ) THEN
88: *
89: * Compute Y(:1:M) = A'*X(:,1:M)
90: *
91: DO 60 K = 1, M
92: *
93: DO 40 I = 1, N, 2
94: Y( I, K ) = COSA*X( I, K ) - SINA*X( I+1, K )
95: Y( I+1, K ) = SINA*X( I, K ) + COSA*X( I+1, K )
96: 40 CONTINUE
97: TEMP = COSB*Y(1,K) + SINB*Y(N,K)
98: DO 50 I = 2, N-1, 2
99: TEMP1 = COSB*Y( I, K ) - SINB*Y( I+1, K )
100: Y( I+1, K ) = SINB*Y( I, K ) + COSB*Y( I+1, K )
101: Y( I, K ) = TEMP1
102: 50 CONTINUE
103: Y( N, K ) = -SINB*Y( 1, K ) + COSB*Y( N, K )
104: Y( 1, K ) = TEMP
105: *
106: 60 CONTINUE
107: *
108: END IF
109: *
110: RETURN
111: *
112: * END OF MVMISG
113: END