Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/ssymv.f @ 17185

Last change on this file since 17185 was 15457, checked in by gkronber, 7 years ago

#2789 added Finbarr O'Sullivan smoothing spline code

File size: 7.6 KB
Line 
1      SUBROUTINE SSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2*     .. Scalar Arguments ..
3      REAL ALPHA,BETA
4      INTEGER INCX,INCY,LDA,N
5      CHARACTER UPLO
6*     ..
7*     .. Array Arguments ..
8      REAL A(LDA,*),X(*),Y(*)
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  SSYMV  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*  Arguments
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  - REAL            .
43*           On entry, ALPHA specifies the scalar alpha.
44*           Unchanged on exit.
45*
46*  A      - REAL             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      - REAL             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   - REAL            .
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      - REAL             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      REAL ONE,ZERO
102      PARAMETER (ONE=1.0E+0,ZERO=0.0E+0)
103*     ..
104*     .. Local Scalars ..
105      REAL TEMP1,TEMP2
106      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
107*     ..
108*     .. External Functions ..
109      LOGICAL LSAME
110      EXTERNAL LSAME
111*     ..
112*     .. External Subroutines ..
113      EXTERNAL XERBLA
114*     ..
115*     .. Intrinsic Functions ..
116      INTRINSIC MAX
117*     ..
118*
119*     Test the input parameters.
120*
121      INFO = 0
122      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
123          INFO = 1
124      ELSE IF (N.LT.0) THEN
125          INFO = 2
126      ELSE IF (LDA.LT.MAX(1,N)) THEN
127          INFO = 5
128      ELSE IF (INCX.EQ.0) THEN
129          INFO = 7
130      ELSE IF (INCY.EQ.0) THEN
131          INFO = 10
132      END IF
133      IF (INFO.NE.0) THEN
134          CALL XERBLA('SSYMV ',INFO)
135          RETURN
136      END IF
137*
138*     Quick return if possible.
139*
140      IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
141*
142*     Set up the start points in  X  and  Y.
143*
144      IF (INCX.GT.0) THEN
145          KX = 1
146      ELSE
147          KX = 1 - (N-1)*INCX
148      END IF
149      IF (INCY.GT.0) THEN
150          KY = 1
151      ELSE
152          KY = 1 - (N-1)*INCY
153      END IF
154*
155*     Start the operations. In this version the elements of A are
156*     accessed sequentially with one pass through the triangular part
157*     of A.
158*
159*     First form  y := beta*y.
160*
161      IF (BETA.NE.ONE) THEN
162          IF (INCY.EQ.1) THEN
163              IF (BETA.EQ.ZERO) THEN
164                  DO 10 I = 1,N
165                      Y(I) = ZERO
166   10             CONTINUE
167              ELSE
168                  DO 20 I = 1,N
169                      Y(I) = BETA*Y(I)
170   20             CONTINUE
171              END IF
172          ELSE
173              IY = KY
174              IF (BETA.EQ.ZERO) THEN
175                  DO 30 I = 1,N
176                      Y(IY) = ZERO
177                      IY = IY + INCY
178   30             CONTINUE
179              ELSE
180                  DO 40 I = 1,N
181                      Y(IY) = BETA*Y(IY)
182                      IY = IY + INCY
183   40             CONTINUE
184              END IF
185          END IF
186      END IF
187      IF (ALPHA.EQ.ZERO) 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 SSYMV .
261*
262      END
Note: See TracBrowser for help on using the repository browser.