KTH framework for Nek5000 toolboxes; testing version  0.0.1
dsymv.f
Go to the documentation of this file.
1  SUBROUTINE dsymv ( UPLO, N, ALPHA, A, LDA, X, INCX,
2  $ BETA, Y, INCY )
3 * .. Scalar Arguments ..
4  DOUBLE PRECISION ALPHA, BETA
5  INTEGER INCX, INCY, LDA, N
6  CHARACTER*1 UPLO
7 * .. Array Arguments ..
8  DOUBLE PRECISION A( LDA, * ), X( * ), Y( * )
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DSYMV 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 matrix.
20 *
21 * Parameters
22 * ==========
23 *
24 * UPLO - CHARACTER*1.
25 * On entry, UPLO specifies whether the upper or lower
26 * triangular part of the array A is to be referenced as
27 * follows:
28 *
29 * UPLO = 'U' or 'u' Only the upper triangular part of A
30 * is to be referenced.
31 *
32 * UPLO = 'L' or 'l' Only the lower triangular part of A
33 * is to be referenced.
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 * ALPHA - DOUBLE PRECISION.
43 * On entry, ALPHA specifies the scalar alpha.
44 * Unchanged on exit.
45 *
46 * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
47 * Before entry with UPLO = 'U' or 'u', the leading n by n
48 * upper triangular part of the array A must contain the upper
49 * triangular part of the symmetric matrix and the strictly
50 * lower triangular part of A is not referenced.
51 * Before entry with UPLO = 'L' or 'l', the leading n by n
52 * lower triangular part of the array A must contain the lower
53 * triangular part of the symmetric matrix and the strictly
54 * upper triangular part of A is not referenced.
55 * Unchanged on exit.
56 *
57 * LDA - INTEGER.
58 * On entry, LDA specifies the first dimension of A as declared
59 * in the calling (sub) program. LDA must be at least
60 * max( 1, n ).
61 * Unchanged on exit.
62 *
63 * X - DOUBLE PRECISION array of dimension at least
64 * ( 1 + ( n - 1 )*abs( INCX ) ).
65 * Before entry, the incremented array X must contain the n
66 * element 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 + ( n - 1 )*abs( INCY ) ).
81 * Before entry, the incremented array Y must contain the n
82 * element vector y. On exit, Y is overwritten by the updated
83 * vector y.
84 *
85 * INCY - INTEGER.
86 * On entry, INCY specifies the increment for the elements of
87 * Y. INCY must not be zero.
88 * Unchanged on exit.
89 *
90 *
91 * Level 2 Blas routine.
92 *
93 * -- Written on 22-October-1986.
94 * Jack Dongarra, Argonne National Lab.
95 * Jeremy Du Croz, Nag Central Office.
96 * Sven Hammarling, Nag Central Office.
97 * Richard Hanson, Sandia National Labs.
98 *
99 *
100 * .. Parameters ..
101  DOUBLE PRECISION ONE , ZERO
102  parameter( one = 1.0d+0, zero = 0.0d+0 )
103 * .. Local Scalars ..
104  DOUBLE PRECISION TEMP1, TEMP2
105  INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY
106 * .. External Functions ..
107  LOGICAL LSAME
108  EXTERNAL lsame
109 * .. External Subroutines ..
110  EXTERNAL xerbla
111 * .. Intrinsic Functions ..
112  INTRINSIC max
113 * ..
114 * .. Executable Statements ..
115 *
116 * Test the input parameters.
117 *
118  info = 0
119  IF ( .NOT.lsame( uplo, 'U' ).AND.
120  $ .NOT.lsame( uplo, 'L' ) )THEN
121  info = 1
122  ELSE IF( n.LT.0 )THEN
123  info = 2
124  ELSE IF( lda.LT.max( 1, n ) )THEN
125  info = 5
126  ELSE IF( incx.EQ.0 )THEN
127  info = 7
128  ELSE IF( incy.EQ.0 )THEN
129  info = 10
130  END IF
131  IF( info.NE.0 )THEN
132  CALL xerbla( 'DSYMV ', info )
133  RETURN
134  END IF
135 *
136 * Quick return if possible.
137 *
138  IF( ( n.EQ.0 ).OR.( ( alpha.EQ.zero ).AND.( beta.EQ.one ) ) )
139  $ RETURN
140 *
141 * Set up the start points in X and Y.
142 *
143  IF( incx.GT.0 )THEN
144  kx = 1
145  ELSE
146  kx = 1 - ( n - 1 )*incx
147  END IF
148  IF( incy.GT.0 )THEN
149  ky = 1
150  ELSE
151  ky = 1 - ( n - 1 )*incy
152  END IF
153 *
154 * Start the operations. In this version the elements of A are
155 * accessed sequentially with one pass through the triangular part
156 * of A.
157 *
158 * First form y := beta*y.
159 *
160  IF( beta.NE.one )THEN
161  IF( incy.EQ.1 )THEN
162  IF( beta.EQ.zero )THEN
163  DO 10, i = 1, n
164  y( i ) = zero
165  10 CONTINUE
166  ELSE
167  DO 20, i = 1, n
168  y( i ) = beta*y( i )
169  20 CONTINUE
170  END IF
171  ELSE
172  iy = ky
173  IF( beta.EQ.zero )THEN
174  DO 30, i = 1, n
175  y( iy ) = zero
176  iy = iy + incy
177  30 CONTINUE
178  ELSE
179  DO 40, i = 1, n
180  y( iy ) = beta*y( iy )
181  iy = iy + incy
182  40 CONTINUE
183  END IF
184  END IF
185  END IF
186  IF( alpha.EQ.zero )
187  $ RETURN
188  IF( lsame( uplo, 'U' ) )THEN
189 *
190 * Form y when A is stored in upper triangle.
191 *
192  IF( ( incx.EQ.1 ).AND.( incy.EQ.1 ) )THEN
193  DO 60, j = 1, n
194  temp1 = alpha*x( j )
195  temp2 = zero
196  DO 50, i = 1, j - 1
197  y( i ) = y( i ) + temp1*a( i, j )
198  temp2 = temp2 + a( i, j )*x( i )
199  50 CONTINUE
200  y( j ) = y( j ) + temp1*a( j, j ) + alpha*temp2
201  60 CONTINUE
202  ELSE
203  jx = kx
204  jy = ky
205  DO 80, j = 1, n
206  temp1 = alpha*x( jx )
207  temp2 = zero
208  ix = kx
209  iy = ky
210  DO 70, i = 1, j - 1
211  y( iy ) = y( iy ) + temp1*a( i, j )
212  temp2 = temp2 + a( i, j )*x( ix )
213  ix = ix + incx
214  iy = iy + incy
215  70 CONTINUE
216  y( jy ) = y( jy ) + temp1*a( j, j ) + alpha*temp2
217  jx = jx + incx
218  jy = jy + incy
219  80 CONTINUE
220  END IF
221  ELSE
222 *
223 * Form y when A is stored in lower triangle.
224 *
225  IF( ( incx.EQ.1 ).AND.( incy.EQ.1 ) )THEN
226  DO 100, j = 1, n
227  temp1 = alpha*x( j )
228  temp2 = zero
229  y( j ) = y( j ) + temp1*a( j, j )
230  DO 90, i = j + 1, n
231  y( i ) = y( i ) + temp1*a( i, j )
232  temp2 = temp2 + a( i, j )*x( i )
233  90 CONTINUE
234  y( j ) = y( j ) + alpha*temp2
235  100 CONTINUE
236  ELSE
237  jx = kx
238  jy = ky
239  DO 120, j = 1, n
240  temp1 = alpha*x( jx )
241  temp2 = zero
242  y( jy ) = y( jy ) + temp1*a( j, j )
243  ix = jx
244  iy = jy
245  DO 110, i = j + 1, n
246  ix = ix + incx
247  iy = iy + incy
248  y( iy ) = y( iy ) + temp1*a( i, j )
249  temp2 = temp2 + a( i, j )*x( ix )
250  110 CONTINUE
251  y( jy ) = y( jy ) + alpha*temp2
252  jx = jx + incx
253  jy = jy + incy
254  120 CONTINUE
255  END IF
256  END IF
257 *
258  RETURN
259 *
260 * End of DSYMV .
261 *
262  END
subroutine dsymv(UPLO, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
Definition: dsymv.f:3
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2