1 SUBROUTINE SSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2 * .. Scalar Arguments ..
3 REAL ALPHA,BETA
4 INTEGER INCX,INCY,K,LDA,N
5 CHARACTER UPLO
6 * ..
7 * .. Array Arguments ..
8 REAL A(LDA,*),X(*),Y(*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * SSBMV performs the matrix-vector operation
15 *
16 * y := alpha*A*x + beta*y,
17 *
18 * where alpha and beta are scalars, x and y are n element vectors and
19 * A is an n by n symmetric band matrix, with k super-diagonals.
20 *
21 * Arguments
22 * ==========
23 *
24 * UPLO - CHARACTER*1.
25 * On entry, UPLO specifies whether the upper or lower
26 * triangular part of the band matrix A is being supplied as
27 * follows:
28 *
29 * UPLO = 'U' or 'u' The upper triangular part of A is
30 * being supplied.
31 *
32 * UPLO = 'L' or 'l' The lower triangular part of A is
33 * being supplied.
34 *
35 * Unchanged on exit.
36 *
37 * N - INTEGER.
38 * On entry, N specifies the order of the matrix A.
39 * N must be at least zero.
40 * Unchanged on exit.
41 *
42 * K - INTEGER.
43 * On entry, K specifies the number of super-diagonals of the
44 * matrix A. K must satisfy 0 .le. K.
45 * Unchanged on exit.
46 *
47 * ALPHA - REAL .
48 * On entry, ALPHA specifies the scalar alpha.
49 * Unchanged on exit.
50 *
51 * A - REAL array of DIMENSION ( LDA, n ).
52 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
53 * by n part of the array A must contain the upper triangular
54 * band part of the symmetric matrix, supplied column by
55 * column, with the leading diagonal of the matrix in row
56 * ( k + 1 ) of the array, the first super-diagonal starting at
57 * position 2 in row k, and so on. The top left k by k triangle
58 * of the array A is not referenced.
59 * The following program segment will transfer the upper
60 * triangular part of a symmetric band matrix from conventional
61 * full matrix storage to band storage:
62 *
63 * DO 20, J = 1, N
64 * M = K + 1 - J
65 * DO 10, I = MAX( 1, J - K ), J
66 * A( M + I, J ) = matrix( I, J )
67 * 10 CONTINUE
68 * 20 CONTINUE
69 *
70 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
71 * by n part of the array A must contain the lower triangular
72 * band part of the symmetric matrix, supplied column by
73 * column, with the leading diagonal of the matrix in row 1 of
74 * the array, the first sub-diagonal starting at position 1 in
75 * row 2, and so on. The bottom right k by k triangle of the
76 * array A is not referenced.
77 * The following program segment will transfer the lower
78 * triangular part of a symmetric band matrix from conventional
79 * full matrix storage to band storage:
80 *
81 * DO 20, J = 1, N
82 * M = 1 - J
83 * DO 10, I = J, MIN( N, J + K )
84 * A( M + I, J ) = matrix( I, J )
85 * 10 CONTINUE
86 * 20 CONTINUE
87 *
88 * Unchanged on exit.
89 *
90 * LDA - INTEGER.
91 * On entry, LDA specifies the first dimension of A as declared
92 * in the calling (sub) program. LDA must be at least
93 * ( k + 1 ).
94 * Unchanged on exit.
95 *
96 * X - REAL array of DIMENSION at least
97 * ( 1 + ( n - 1 )*abs( INCX ) ).
98 * Before entry, the incremented array X must contain the
99 * vector x.
100 * Unchanged on exit.
101 *
102 * INCX - INTEGER.
103 * On entry, INCX specifies the increment for the elements of
104 * X. INCX must not be zero.
105 * Unchanged on exit.
106 *
107 * BETA - REAL .
108 * On entry, BETA specifies the scalar beta.
109 * Unchanged on exit.
110 *
111 * Y - REAL array of DIMENSION at least
112 * ( 1 + ( n - 1 )*abs( INCY ) ).
113 * Before entry, the incremented array Y must contain the
114 * vector y. On exit, Y is overwritten by the updated vector y.
115 *
116 * INCY - INTEGER.
117 * On entry, INCY specifies the increment for the elements of
118 * Y. INCY must not be zero.
119 * Unchanged on exit.
120 *
121 * Further Details
122 * ===============
123 *
124 * Level 2 Blas routine.
125 * The vector and matrix arguments are not referenced when N = 0, or M = 0
126 *
127 * -- Written on 22-October-1986.
128 * Jack Dongarra, Argonne National Lab.
129 * Jeremy Du Croz, Nag Central Office.
130 * Sven Hammarling, Nag Central Office.
131 * Richard Hanson, Sandia National Labs.
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136 REAL ONE,ZERO
137 PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
138 * ..
139 * .. Local Scalars ..
140 REAL TEMP1,TEMP2
141 INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
142 * ..
143 * .. External Functions ..
144 LOGICAL LSAME
145 EXTERNAL LSAME
146 * ..
147 * .. External Subroutines ..
148 EXTERNAL XERBLA
149 * ..
150 * .. Intrinsic Functions ..
151 INTRINSIC MAX,MIN
152 * ..
153 *
154 * Test the input parameters.
155 *
156 INFO = 0
157 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
158 INFO = 1
159 ELSE IF (N.LT.0) THEN
160 INFO = 2
161 ELSE IF (K.LT.0) THEN
162 INFO = 3
163 ELSE IF (LDA.LT. (K+1)) THEN
164 INFO = 6
165 ELSE IF (INCX.EQ.0) THEN
166 INFO = 8
167 ELSE IF (INCY.EQ.0) THEN
168 INFO = 11
169 END IF
170 IF (INFO.NE.0) THEN
171 CALL XERBLA('SSBMV ',INFO)
172 RETURN
173 END IF
174 *
175 * Quick return if possible.
176 *
177 IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
178 *
179 * Set up the start points in X and Y.
180 *
181 IF (INCX.GT.0) THEN
182 KX = 1
183 ELSE
184 KX = 1 - (N-1)*INCX
185 END IF
186 IF (INCY.GT.0) THEN
187 KY = 1
188 ELSE
189 KY = 1 - (N-1)*INCY
190 END IF
191 *
192 * Start the operations. In this version the elements of the array A
193 * are accessed sequentially with one pass through A.
194 *
195 * First form y := beta*y.
196 *
197 IF (BETA.NE.ONE) THEN
198 IF (INCY.EQ.1) THEN
199 IF (BETA.EQ.ZERO) THEN
200 DO 10 I = 1,N
201 Y(I) = ZERO
202 10 CONTINUE
203 ELSE
204 DO 20 I = 1,N
205 Y(I) = BETA*Y(I)
206 20 CONTINUE
207 END IF
208 ELSE
209 IY = KY
210 IF (BETA.EQ.ZERO) THEN
211 DO 30 I = 1,N
212 Y(IY) = ZERO
213 IY = IY + INCY
214 30 CONTINUE
215 ELSE
216 DO 40 I = 1,N
217 Y(IY) = BETA*Y(IY)
218 IY = IY + INCY
219 40 CONTINUE
220 END IF
221 END IF
222 END IF
223 IF (ALPHA.EQ.ZERO) RETURN
224 IF (LSAME(UPLO,'U')) THEN
225 *
226 * Form y when upper triangle of A is stored.
227 *
228 KPLUS1 = K + 1
229 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
230 DO 60 J = 1,N
231 TEMP1 = ALPHA*X(J)
232 TEMP2 = ZERO
233 L = KPLUS1 - J
234 DO 50 I = MAX(1,J-K),J - 1
235 Y(I) = Y(I) + TEMP1*A(L+I,J)
236 TEMP2 = TEMP2 + A(L+I,J)*X(I)
237 50 CONTINUE
238 Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
239 60 CONTINUE
240 ELSE
241 JX = KX
242 JY = KY
243 DO 80 J = 1,N
244 TEMP1 = ALPHA*X(JX)
245 TEMP2 = ZERO
246 IX = KX
247 IY = KY
248 L = KPLUS1 - J
249 DO 70 I = MAX(1,J-K),J - 1
250 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
251 TEMP2 = TEMP2 + A(L+I,J)*X(IX)
252 IX = IX + INCX
253 IY = IY + INCY
254 70 CONTINUE
255 Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
256 JX = JX + INCX
257 JY = JY + INCY
258 IF (J.GT.K) THEN
259 KX = KX + INCX
260 KY = KY + INCY
261 END IF
262 80 CONTINUE
263 END IF
264 ELSE
265 *
266 * Form y when lower triangle of A is stored.
267 *
268 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
269 DO 100 J = 1,N
270 TEMP1 = ALPHA*X(J)
271 TEMP2 = ZERO
272 Y(J) = Y(J) + TEMP1*A(1,J)
273 L = 1 - J
274 DO 90 I = J + 1,MIN(N,J+K)
275 Y(I) = Y(I) + TEMP1*A(L+I,J)
276 TEMP2 = TEMP2 + A(L+I,J)*X(I)
277 90 CONTINUE
278 Y(J) = Y(J) + ALPHA*TEMP2
279 100 CONTINUE
280 ELSE
281 JX = KX
282 JY = KY
283 DO 120 J = 1,N
284 TEMP1 = ALPHA*X(JX)
285 TEMP2 = ZERO
286 Y(JY) = Y(JY) + TEMP1*A(1,J)
287 L = 1 - J
288 IX = JX
289 IY = JY
290 DO 110 I = J + 1,MIN(N,J+K)
291 IX = IX + INCX
292 IY = IY + INCY
293 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
294 TEMP2 = TEMP2 + A(L+I,J)*X(IX)
295 110 CONTINUE
296 Y(JY) = Y(JY) + ALPHA*TEMP2
297 JX = JX + INCX
298 JY = JY + INCY
299 120 CONTINUE
300 END IF
301 END IF
302 *
303 RETURN
304 *
305 * End of SSBMV .
306 *
307 END
2 * .. Scalar Arguments ..
3 REAL ALPHA,BETA
4 INTEGER INCX,INCY,K,LDA,N
5 CHARACTER UPLO
6 * ..
7 * .. Array Arguments ..
8 REAL A(LDA,*),X(*),Y(*)
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * SSBMV performs the matrix-vector operation
15 *
16 * y := alpha*A*x + beta*y,
17 *
18 * where alpha and beta are scalars, x and y are n element vectors and
19 * A is an n by n symmetric band matrix, with k super-diagonals.
20 *
21 * Arguments
22 * ==========
23 *
24 * UPLO - CHARACTER*1.
25 * On entry, UPLO specifies whether the upper or lower
26 * triangular part of the band matrix A is being supplied as
27 * follows:
28 *
29 * UPLO = 'U' or 'u' The upper triangular part of A is
30 * being supplied.
31 *
32 * UPLO = 'L' or 'l' The lower triangular part of A is
33 * being supplied.
34 *
35 * Unchanged on exit.
36 *
37 * N - INTEGER.
38 * On entry, N specifies the order of the matrix A.
39 * N must be at least zero.
40 * Unchanged on exit.
41 *
42 * K - INTEGER.
43 * On entry, K specifies the number of super-diagonals of the
44 * matrix A. K must satisfy 0 .le. K.
45 * Unchanged on exit.
46 *
47 * ALPHA - REAL .
48 * On entry, ALPHA specifies the scalar alpha.
49 * Unchanged on exit.
50 *
51 * A - REAL array of DIMENSION ( LDA, n ).
52 * Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
53 * by n part of the array A must contain the upper triangular
54 * band part of the symmetric matrix, supplied column by
55 * column, with the leading diagonal of the matrix in row
56 * ( k + 1 ) of the array, the first super-diagonal starting at
57 * position 2 in row k, and so on. The top left k by k triangle
58 * of the array A is not referenced.
59 * The following program segment will transfer the upper
60 * triangular part of a symmetric band matrix from conventional
61 * full matrix storage to band storage:
62 *
63 * DO 20, J = 1, N
64 * M = K + 1 - J
65 * DO 10, I = MAX( 1, J - K ), J
66 * A( M + I, J ) = matrix( I, J )
67 * 10 CONTINUE
68 * 20 CONTINUE
69 *
70 * Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
71 * by n part of the array A must contain the lower triangular
72 * band part of the symmetric matrix, supplied column by
73 * column, with the leading diagonal of the matrix in row 1 of
74 * the array, the first sub-diagonal starting at position 1 in
75 * row 2, and so on. The bottom right k by k triangle of the
76 * array A is not referenced.
77 * The following program segment will transfer the lower
78 * triangular part of a symmetric band matrix from conventional
79 * full matrix storage to band storage:
80 *
81 * DO 20, J = 1, N
82 * M = 1 - J
83 * DO 10, I = J, MIN( N, J + K )
84 * A( M + I, J ) = matrix( I, J )
85 * 10 CONTINUE
86 * 20 CONTINUE
87 *
88 * Unchanged on exit.
89 *
90 * LDA - INTEGER.
91 * On entry, LDA specifies the first dimension of A as declared
92 * in the calling (sub) program. LDA must be at least
93 * ( k + 1 ).
94 * Unchanged on exit.
95 *
96 * X - REAL array of DIMENSION at least
97 * ( 1 + ( n - 1 )*abs( INCX ) ).
98 * Before entry, the incremented array X must contain the
99 * vector x.
100 * Unchanged on exit.
101 *
102 * INCX - INTEGER.
103 * On entry, INCX specifies the increment for the elements of
104 * X. INCX must not be zero.
105 * Unchanged on exit.
106 *
107 * BETA - REAL .
108 * On entry, BETA specifies the scalar beta.
109 * Unchanged on exit.
110 *
111 * Y - REAL array of DIMENSION at least
112 * ( 1 + ( n - 1 )*abs( INCY ) ).
113 * Before entry, the incremented array Y must contain the
114 * vector y. On exit, Y is overwritten by the updated vector y.
115 *
116 * INCY - INTEGER.
117 * On entry, INCY specifies the increment for the elements of
118 * Y. INCY must not be zero.
119 * Unchanged on exit.
120 *
121 * Further Details
122 * ===============
123 *
124 * Level 2 Blas routine.
125 * The vector and matrix arguments are not referenced when N = 0, or M = 0
126 *
127 * -- Written on 22-October-1986.
128 * Jack Dongarra, Argonne National Lab.
129 * Jeremy Du Croz, Nag Central Office.
130 * Sven Hammarling, Nag Central Office.
131 * Richard Hanson, Sandia National Labs.
132 *
133 * =====================================================================
134 *
135 * .. Parameters ..
136 REAL ONE,ZERO
137 PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
138 * ..
139 * .. Local Scalars ..
140 REAL TEMP1,TEMP2
141 INTEGER I,INFO,IX,IY,J,JX,JY,KPLUS1,KX,KY,L
142 * ..
143 * .. External Functions ..
144 LOGICAL LSAME
145 EXTERNAL LSAME
146 * ..
147 * .. External Subroutines ..
148 EXTERNAL XERBLA
149 * ..
150 * .. Intrinsic Functions ..
151 INTRINSIC MAX,MIN
152 * ..
153 *
154 * Test the input parameters.
155 *
156 INFO = 0
157 IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
158 INFO = 1
159 ELSE IF (N.LT.0) THEN
160 INFO = 2
161 ELSE IF (K.LT.0) THEN
162 INFO = 3
163 ELSE IF (LDA.LT. (K+1)) THEN
164 INFO = 6
165 ELSE IF (INCX.EQ.0) THEN
166 INFO = 8
167 ELSE IF (INCY.EQ.0) THEN
168 INFO = 11
169 END IF
170 IF (INFO.NE.0) THEN
171 CALL XERBLA('SSBMV ',INFO)
172 RETURN
173 END IF
174 *
175 * Quick return if possible.
176 *
177 IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
178 *
179 * Set up the start points in X and Y.
180 *
181 IF (INCX.GT.0) THEN
182 KX = 1
183 ELSE
184 KX = 1 - (N-1)*INCX
185 END IF
186 IF (INCY.GT.0) THEN
187 KY = 1
188 ELSE
189 KY = 1 - (N-1)*INCY
190 END IF
191 *
192 * Start the operations. In this version the elements of the array A
193 * are accessed sequentially with one pass through A.
194 *
195 * First form y := beta*y.
196 *
197 IF (BETA.NE.ONE) THEN
198 IF (INCY.EQ.1) THEN
199 IF (BETA.EQ.ZERO) THEN
200 DO 10 I = 1,N
201 Y(I) = ZERO
202 10 CONTINUE
203 ELSE
204 DO 20 I = 1,N
205 Y(I) = BETA*Y(I)
206 20 CONTINUE
207 END IF
208 ELSE
209 IY = KY
210 IF (BETA.EQ.ZERO) THEN
211 DO 30 I = 1,N
212 Y(IY) = ZERO
213 IY = IY + INCY
214 30 CONTINUE
215 ELSE
216 DO 40 I = 1,N
217 Y(IY) = BETA*Y(IY)
218 IY = IY + INCY
219 40 CONTINUE
220 END IF
221 END IF
222 END IF
223 IF (ALPHA.EQ.ZERO) RETURN
224 IF (LSAME(UPLO,'U')) THEN
225 *
226 * Form y when upper triangle of A is stored.
227 *
228 KPLUS1 = K + 1
229 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
230 DO 60 J = 1,N
231 TEMP1 = ALPHA*X(J)
232 TEMP2 = ZERO
233 L = KPLUS1 - J
234 DO 50 I = MAX(1,J-K),J - 1
235 Y(I) = Y(I) + TEMP1*A(L+I,J)
236 TEMP2 = TEMP2 + A(L+I,J)*X(I)
237 50 CONTINUE
238 Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
239 60 CONTINUE
240 ELSE
241 JX = KX
242 JY = KY
243 DO 80 J = 1,N
244 TEMP1 = ALPHA*X(JX)
245 TEMP2 = ZERO
246 IX = KX
247 IY = KY
248 L = KPLUS1 - J
249 DO 70 I = MAX(1,J-K),J - 1
250 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
251 TEMP2 = TEMP2 + A(L+I,J)*X(IX)
252 IX = IX + INCX
253 IY = IY + INCY
254 70 CONTINUE
255 Y(JY) = Y(JY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2
256 JX = JX + INCX
257 JY = JY + INCY
258 IF (J.GT.K) THEN
259 KX = KX + INCX
260 KY = KY + INCY
261 END IF
262 80 CONTINUE
263 END IF
264 ELSE
265 *
266 * Form y when lower triangle of A is stored.
267 *
268 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
269 DO 100 J = 1,N
270 TEMP1 = ALPHA*X(J)
271 TEMP2 = ZERO
272 Y(J) = Y(J) + TEMP1*A(1,J)
273 L = 1 - J
274 DO 90 I = J + 1,MIN(N,J+K)
275 Y(I) = Y(I) + TEMP1*A(L+I,J)
276 TEMP2 = TEMP2 + A(L+I,J)*X(I)
277 90 CONTINUE
278 Y(J) = Y(J) + ALPHA*TEMP2
279 100 CONTINUE
280 ELSE
281 JX = KX
282 JY = KY
283 DO 120 J = 1,N
284 TEMP1 = ALPHA*X(JX)
285 TEMP2 = ZERO
286 Y(JY) = Y(JY) + TEMP1*A(1,J)
287 L = 1 - J
288 IX = JX
289 IY = JY
290 DO 110 I = J + 1,MIN(N,J+K)
291 IX = IX + INCX
292 IY = IY + INCY
293 Y(IY) = Y(IY) + TEMP1*A(L+I,J)
294 TEMP2 = TEMP2 + A(L+I,J)*X(IX)
295 110 CONTINUE
296 Y(JY) = Y(JY) + ALPHA*TEMP2
297 JX = JX + INCX
298 JY = JY + INCY
299 120 CONTINUE
300 END IF
301 END IF
302 *
303 RETURN
304 *
305 * End of SSBMV .
306 *
307 END