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**T*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**T*x + beta*y.
31 *
32 * TRANS = 'C' or 'c' y := alpha*A**T*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 * The vector and matrix arguments are not referenced when N = 0, or M = 0
97 *
98 * -- Written on 22-October-1986.
99 * Jack Dongarra, Argonne National Lab.
100 * Jeremy Du Croz, Nag Central Office.
101 * Sven Hammarling, Nag Central Office.
102 * Richard Hanson, Sandia National Labs.
103 *
104 * =====================================================================
105 *
106 * .. Parameters ..
107 DOUBLE PRECISION ONE,ZERO
108 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
109 * ..
110 * .. Local Scalars ..
111 DOUBLE PRECISION TEMP
112 INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
113 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 EXTERNAL LSAME
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL XERBLA
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC MAX
123 * ..
124 *
125 * Test the input parameters.
126 *
127 INFO = 0
128 IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
129 + .NOT.LSAME(TRANS,'C')) THEN
130 INFO = 1
131 ELSE IF (M.LT.0) THEN
132 INFO = 2
133 ELSE IF (N.LT.0) THEN
134 INFO = 3
135 ELSE IF (LDA.LT.MAX(1,M)) THEN
136 INFO = 6
137 ELSE IF (INCX.EQ.0) THEN
138 INFO = 8
139 ELSE IF (INCY.EQ.0) THEN
140 INFO = 11
141 END IF
142 IF (INFO.NE.0) THEN
143 CALL XERBLA('DGEMV ',INFO)
144 RETURN
145 END IF
146 *
147 * Quick return if possible.
148 *
149 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
150 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
151 *
152 * Set LENX and LENY, the lengths of the vectors x and y, and set
153 * up the start points in X and Y.
154 *
155 IF (LSAME(TRANS,'N')) THEN
156 LENX = N
157 LENY = M
158 ELSE
159 LENX = M
160 LENY = N
161 END IF
162 IF (INCX.GT.0) THEN
163 KX = 1
164 ELSE
165 KX = 1 - (LENX-1)*INCX
166 END IF
167 IF (INCY.GT.0) THEN
168 KY = 1
169 ELSE
170 KY = 1 - (LENY-1)*INCY
171 END IF
172 *
173 * Start the operations. In this version the elements of A are
174 * accessed sequentially with one pass through A.
175 *
176 * First form y := beta*y.
177 *
178 IF (BETA.NE.ONE) THEN
179 IF (INCY.EQ.1) THEN
180 IF (BETA.EQ.ZERO) THEN
181 DO 10 I = 1,LENY
182 Y(I) = ZERO
183 10 CONTINUE
184 ELSE
185 DO 20 I = 1,LENY
186 Y(I) = BETA*Y(I)
187 20 CONTINUE
188 END IF
189 ELSE
190 IY = KY
191 IF (BETA.EQ.ZERO) THEN
192 DO 30 I = 1,LENY
193 Y(IY) = ZERO
194 IY = IY + INCY
195 30 CONTINUE
196 ELSE
197 DO 40 I = 1,LENY
198 Y(IY) = BETA*Y(IY)
199 IY = IY + INCY
200 40 CONTINUE
201 END IF
202 END IF
203 END IF
204 IF (ALPHA.EQ.ZERO) RETURN
205 IF (LSAME(TRANS,'N')) THEN
206 *
207 * Form y := alpha*A*x + y.
208 *
209 JX = KX
210 IF (INCY.EQ.1) THEN
211 DO 60 J = 1,N
212 IF (X(JX).NE.ZERO) THEN
213 TEMP = ALPHA*X(JX)
214 DO 50 I = 1,M
215 Y(I) = Y(I) + TEMP*A(I,J)
216 50 CONTINUE
217 END IF
218 JX = JX + INCX
219 60 CONTINUE
220 ELSE
221 DO 80 J = 1,N
222 IF (X(JX).NE.ZERO) THEN
223 TEMP = ALPHA*X(JX)
224 IY = KY
225 DO 70 I = 1,M
226 Y(IY) = Y(IY) + TEMP*A(I,J)
227 IY = IY + INCY
228 70 CONTINUE
229 END IF
230 JX = JX + INCX
231 80 CONTINUE
232 END IF
233 ELSE
234 *
235 * Form y := alpha*A**T*x + y.
236 *
237 JY = KY
238 IF (INCX.EQ.1) THEN
239 DO 100 J = 1,N
240 TEMP = ZERO
241 DO 90 I = 1,M
242 TEMP = TEMP + A(I,J)*X(I)
243 90 CONTINUE
244 Y(JY) = Y(JY) + ALPHA*TEMP
245 JY = JY + INCY
246 100 CONTINUE
247 ELSE
248 DO 120 J = 1,N
249 TEMP = ZERO
250 IX = KX
251 DO 110 I = 1,M
252 TEMP = TEMP + A(I,J)*X(IX)
253 IX = IX + INCX
254 110 CONTINUE
255 Y(JY) = Y(JY) + ALPHA*TEMP
256 JY = JY + INCY
257 120 CONTINUE
258 END IF
259 END IF
260 *
261 RETURN
262 *
263 * End of DGEMV .
264 *
265 END
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**T*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**T*x + beta*y.
31 *
32 * TRANS = 'C' or 'c' y := alpha*A**T*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 * The vector and matrix arguments are not referenced when N = 0, or M = 0
97 *
98 * -- Written on 22-October-1986.
99 * Jack Dongarra, Argonne National Lab.
100 * Jeremy Du Croz, Nag Central Office.
101 * Sven Hammarling, Nag Central Office.
102 * Richard Hanson, Sandia National Labs.
103 *
104 * =====================================================================
105 *
106 * .. Parameters ..
107 DOUBLE PRECISION ONE,ZERO
108 PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
109 * ..
110 * .. Local Scalars ..
111 DOUBLE PRECISION TEMP
112 INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
113 * ..
114 * .. External Functions ..
115 LOGICAL LSAME
116 EXTERNAL LSAME
117 * ..
118 * .. External Subroutines ..
119 EXTERNAL XERBLA
120 * ..
121 * .. Intrinsic Functions ..
122 INTRINSIC MAX
123 * ..
124 *
125 * Test the input parameters.
126 *
127 INFO = 0
128 IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
129 + .NOT.LSAME(TRANS,'C')) THEN
130 INFO = 1
131 ELSE IF (M.LT.0) THEN
132 INFO = 2
133 ELSE IF (N.LT.0) THEN
134 INFO = 3
135 ELSE IF (LDA.LT.MAX(1,M)) THEN
136 INFO = 6
137 ELSE IF (INCX.EQ.0) THEN
138 INFO = 8
139 ELSE IF (INCY.EQ.0) THEN
140 INFO = 11
141 END IF
142 IF (INFO.NE.0) THEN
143 CALL XERBLA('DGEMV ',INFO)
144 RETURN
145 END IF
146 *
147 * Quick return if possible.
148 *
149 IF ((M.EQ.0) .OR. (N.EQ.0) .OR.
150 + ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
151 *
152 * Set LENX and LENY, the lengths of the vectors x and y, and set
153 * up the start points in X and Y.
154 *
155 IF (LSAME(TRANS,'N')) THEN
156 LENX = N
157 LENY = M
158 ELSE
159 LENX = M
160 LENY = N
161 END IF
162 IF (INCX.GT.0) THEN
163 KX = 1
164 ELSE
165 KX = 1 - (LENX-1)*INCX
166 END IF
167 IF (INCY.GT.0) THEN
168 KY = 1
169 ELSE
170 KY = 1 - (LENY-1)*INCY
171 END IF
172 *
173 * Start the operations. In this version the elements of A are
174 * accessed sequentially with one pass through A.
175 *
176 * First form y := beta*y.
177 *
178 IF (BETA.NE.ONE) THEN
179 IF (INCY.EQ.1) THEN
180 IF (BETA.EQ.ZERO) THEN
181 DO 10 I = 1,LENY
182 Y(I) = ZERO
183 10 CONTINUE
184 ELSE
185 DO 20 I = 1,LENY
186 Y(I) = BETA*Y(I)
187 20 CONTINUE
188 END IF
189 ELSE
190 IY = KY
191 IF (BETA.EQ.ZERO) THEN
192 DO 30 I = 1,LENY
193 Y(IY) = ZERO
194 IY = IY + INCY
195 30 CONTINUE
196 ELSE
197 DO 40 I = 1,LENY
198 Y(IY) = BETA*Y(IY)
199 IY = IY + INCY
200 40 CONTINUE
201 END IF
202 END IF
203 END IF
204 IF (ALPHA.EQ.ZERO) RETURN
205 IF (LSAME(TRANS,'N')) THEN
206 *
207 * Form y := alpha*A*x + y.
208 *
209 JX = KX
210 IF (INCY.EQ.1) THEN
211 DO 60 J = 1,N
212 IF (X(JX).NE.ZERO) THEN
213 TEMP = ALPHA*X(JX)
214 DO 50 I = 1,M
215 Y(I) = Y(I) + TEMP*A(I,J)
216 50 CONTINUE
217 END IF
218 JX = JX + INCX
219 60 CONTINUE
220 ELSE
221 DO 80 J = 1,N
222 IF (X(JX).NE.ZERO) THEN
223 TEMP = ALPHA*X(JX)
224 IY = KY
225 DO 70 I = 1,M
226 Y(IY) = Y(IY) + TEMP*A(I,J)
227 IY = IY + INCY
228 70 CONTINUE
229 END IF
230 JX = JX + INCX
231 80 CONTINUE
232 END IF
233 ELSE
234 *
235 * Form y := alpha*A**T*x + y.
236 *
237 JY = KY
238 IF (INCX.EQ.1) THEN
239 DO 100 J = 1,N
240 TEMP = ZERO
241 DO 90 I = 1,M
242 TEMP = TEMP + A(I,J)*X(I)
243 90 CONTINUE
244 Y(JY) = Y(JY) + ALPHA*TEMP
245 JY = JY + INCY
246 100 CONTINUE
247 ELSE
248 DO 120 J = 1,N
249 TEMP = ZERO
250 IX = KX
251 DO 110 I = 1,M
252 TEMP = TEMP + A(I,J)*X(IX)
253 IX = IX + INCX
254 110 CONTINUE
255 Y(JY) = Y(JY) + ALPHA*TEMP
256 JY = JY + INCY
257 120 CONTINUE
258 END IF
259 END IF
260 *
261 RETURN
262 *
263 * End of DGEMV .
264 *
265 END