1: SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA,BETA
4: INTEGER INCX,INCY,LDA,M,N
5: CHARACTER TRANS
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),X(*),Y(*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DGEMV performs one of the matrix-vector operations
15: *
16: * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y,
17: *
18: * where alpha and beta are scalars, x and y are vectors and A is an
19: * m by n matrix.
20: *
21: * Arguments
22: * ==========
23: *
24: * TRANS - CHARACTER*1.
25: * On entry, TRANS specifies the operation to be performed as
26: * follows:
27: *
28: * TRANS = 'N' or 'n' y := alpha*A*x + beta*y.
29: *
30: * TRANS = 'T' or 't' y := alpha*A'*x + beta*y.
31: *
32: * TRANS = 'C' or 'c' y := alpha*A'*x + beta*y.
33: *
34: * Unchanged on exit.
35: *
36: * M - INTEGER.
37: * On entry, M specifies the number of rows of the matrix A.
38: * M must be at least zero.
39: * Unchanged on exit.
40: *
41: * N - INTEGER.
42: * On entry, N specifies the number of columns of the matrix A.
43: * N must be at least zero.
44: * Unchanged on exit.
45: *
46: * ALPHA - DOUBLE PRECISION.
47: * On entry, ALPHA specifies the scalar alpha.
48: * Unchanged on exit.
49: *
50: * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
51: * Before entry, the leading m by n part of the array A must
52: * contain the matrix of coefficients.
53: * Unchanged on exit.
54: *
55: * LDA - INTEGER.
56: * On entry, LDA specifies the first dimension of A as declared
57: * in the calling (sub) program. LDA must be at least
58: * max( 1, m ).
59: * Unchanged on exit.
60: *
61: * X - DOUBLE PRECISION array of DIMENSION at least
62: * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
63: * and at least
64: * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
65: * Before entry, the incremented array X must contain the
66: * vector x.
67: * Unchanged on exit.
68: *
69: * INCX - INTEGER.
70: * On entry, INCX specifies the increment for the elements of
71: * X. INCX must not be zero.
72: * Unchanged on exit.
73: *
74: * BETA - DOUBLE PRECISION.
75: * On entry, BETA specifies the scalar beta. When BETA is
76: * supplied as zero then Y need not be set on input.
77: * Unchanged on exit.
78: *
79: * Y - DOUBLE PRECISION array of DIMENSION at least
80: * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
81: * and at least
82: * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
83: * Before entry with BETA non-zero, the incremented array Y
84: * must contain the vector y. On exit, Y is overwritten by the
85: * updated vector y.
86: *
87: * INCY - INTEGER.
88: * On entry, INCY specifies the increment for the elements of
89: * Y. INCY must not be zero.
90: * Unchanged on exit.
91: *
92: * Further Details
93: * ===============
94: *
95: * Level 2 Blas routine.
96: *
97: * -- Written on 22-October-1986.
98: * Jack Dongarra, Argonne National Lab.
99: * Jeremy Du Croz, Nag Central Office.
100: * Sven Hammarling, Nag Central Office.
101: * Richard Hanson, Sandia National Labs.
102: *
103: * =====================================================================
104: *
105: * .. Parameters ..
106: DOUBLE PRECISION ONE,ZERO
107: PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
108: * ..
109: * .. Local Scalars ..
110: DOUBLE PRECISION TEMP
111: INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
112: * ..
113: * .. External Functions ..
114: LOGICAL LSAME
115: EXTERNAL LSAME
116: * ..
117: * .. External Subroutines ..
118: EXTERNAL XERBLA
119: * ..
120: * .. Intrinsic Functions ..
121: INTRINSIC MAX
122: * ..
123: *
124: * Test the input parameters.
125: *
126: INFO = 0
127: IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
128: + .NOT.LSAME(TRANS,'C')) THEN
129: INFO = 1
130: ELSE IF (M.LT.0) THEN
131: INFO = 2
132: ELSE IF (N.LT.0) THEN
133: INFO = 3
134: ELSE IF (LDA.LT.MAX(1,M)) THEN
135: INFO = 6
136: ELSE IF (INCX.EQ.0) THEN
137: INFO = 8
138: ELSE IF (INCY.EQ.0) THEN
139: INFO = 11
140: END IF
141: IF (INFO.NE.0) THEN
142: CALL XERBLA('DGEMV ',INFO)
143: RETURN
144: END IF
145: *
146: * Quick return if possible.
147: *
148: IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
149: + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
150: *
151: * Set LENX and LENY, the lengths of the vectors x and y, and set
152: * up the start points in X and Y.
153: *
154: IF (LSAME(TRANS,'N')) THEN
155: LENX = N
156: LENY = M
157: ELSE
158: LENX = M
159: LENY = N
160: END IF
161: IF (INCX.GT.0) THEN
162: KX = 1
163: ELSE
164: KX = 1 - (LENX-1)*INCX
165: END IF
166: IF (INCY.GT.0) THEN
167: KY = 1
168: ELSE
169: KY = 1 - (LENY-1)*INCY
170: END IF
171: *
172: * Start the operations. In this version the elements of A are
173: * accessed sequentially with one pass through A.
174: *
175: * First form y := beta*y.
176: *
177: IF (BETA.NE.ONE) THEN
178: IF (INCY.EQ.1) THEN
179: IF (BETA.EQ.ZERO) THEN
180: DO 10 I = 1,LENY
181: Y(I) = ZERO
182: 10 CONTINUE
183: ELSE
184: DO 20 I = 1,LENY
185: Y(I) = BETA*Y(I)
186: 20 CONTINUE
187: END IF
188: ELSE
189: IY = KY
190: IF (BETA.EQ.ZERO) THEN
191: DO 30 I = 1,LENY
192: Y(IY) = ZERO
193: IY = IY + INCY
194: 30 CONTINUE
195: ELSE
196: DO 40 I = 1,LENY
197: Y(IY) = BETA*Y(IY)
198: IY = IY + INCY
199: 40 CONTINUE
200: END IF
201: END IF
202: END IF
203: IF (ALPHA.EQ.ZERO) RETURN
204: IF (LSAME(TRANS,'N')) THEN
205: *
206: * Form y := alpha*A*x + y.
207: *
208: JX = KX
209: IF (INCY.EQ.1) THEN
210: DO 60 J = 1,N
211: IF (X(JX).NE.ZERO) THEN
212: TEMP = ALPHA*X(JX)
213: DO 50 I = 1,M
214: Y(I) = Y(I) + TEMP*A(I,J)
215: 50 CONTINUE
216: END IF
217: JX = JX + INCX
218: 60 CONTINUE
219: ELSE
220: DO 80 J = 1,N
221: IF (X(JX).NE.ZERO) THEN
222: TEMP = ALPHA*X(JX)
223: IY = KY
224: DO 70 I = 1,M
225: Y(IY) = Y(IY) + TEMP*A(I,J)
226: IY = IY + INCY
227: 70 CONTINUE
228: END IF
229: JX = JX + INCX
230: 80 CONTINUE
231: END IF
232: ELSE
233: *
234: * Form y := alpha*A'*x + y.
235: *
236: JY = KY
237: IF (INCX.EQ.1) THEN
238: DO 100 J = 1,N
239: TEMP = ZERO
240: DO 90 I = 1,M
241: TEMP = TEMP + A(I,J)*X(I)
242: 90 CONTINUE
243: Y(JY) = Y(JY) + ALPHA*TEMP
244: JY = JY + INCY
245: 100 CONTINUE
246: ELSE
247: DO 120 J = 1,N
248: TEMP = ZERO
249: IX = KX
250: DO 110 I = 1,M
251: TEMP = TEMP + A(I,J)*X(IX)
252: IX = IX + INCX
253: 110 CONTINUE
254: Y(JY) = Y(JY) + ALPHA*TEMP
255: JY = JY + INCY
256: 120 CONTINUE
257: END IF
258: END IF
259: *
260: RETURN
261: *
262: * End of DGEMV .
263: *
264: END
CVSweb interface <joel.bertrand@systella.fr>