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