← Back to branch summary

~ubuntu-branches/debian/sid/lammps/sid

« back to all changes in this revision

Viewing changes to lib/linalg/dsymv.f

  • Committer: Package Import Robot
  • Author(s): Anton Gladky
  • Date: 2015-04-29 23:44:49 UTC
  • mfrom: (5.1.3 experimental)
  • Revision ID: package-import@ubuntu.com-20150429234449-mbhy9utku6hp6oq8
Tags: 0~20150313.gitfa668e1-1
Upload into unstable.

Show diffs side-by-side

added added

removed removed

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