Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2789 added Finbarr O'Sullivan smoothing spline code

File size: 7.8 KB
Line 
1      SUBROUTINE CHEMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
2*     .. Scalar Arguments ..
3      COMPLEX ALPHA,BETA
4      INTEGER INCX,INCY,LDA,N
5      CHARACTER UPLO
6*     ..
7*     .. Array Arguments ..
8      COMPLEX A(LDA,*),X(*),Y(*)
9*     ..
10*
11*  Purpose
12*  =======
13*
14*  CHEMV  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 hermitian 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  - COMPLEX         .
43*           On entry, ALPHA specifies the scalar alpha.
44*           Unchanged on exit.
45*
46*  A      - COMPLEX          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 hermitian 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 hermitian matrix and the strictly
54*           upper triangular part of A is not referenced.
55*           Note that the imaginary parts of the diagonal elements need
56*           not be set and are assumed to be zero.
57*           Unchanged on exit.
58*
59*  LDA    - INTEGER.
60*           On entry, LDA specifies the first dimension of A as declared
61*           in the calling (sub) program. LDA must be at least
62*           max( 1, n ).
63*           Unchanged on exit.
64*
65*  X      - COMPLEX          array of dimension at least
66*           ( 1 + ( n - 1 )*abs( INCX ) ).
67*           Before entry, the incremented array X must contain the n
68*           element vector x.
69*           Unchanged on exit.
70*
71*  INCX   - INTEGER.
72*           On entry, INCX specifies the increment for the elements of
73*           X. INCX must not be zero.
74*           Unchanged on exit.
75*
76*  BETA   - COMPLEX         .
77*           On entry, BETA specifies the scalar beta. When BETA is
78*           supplied as zero then Y need not be set on input.
79*           Unchanged on exit.
80*
81*  Y      - COMPLEX          array of dimension at least
82*           ( 1 + ( n - 1 )*abs( INCY ) ).
83*           Before entry, the incremented array Y must contain the n
84*           element vector y. On exit, Y is overwritten by the updated
85*           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*
93*  Level 2 Blas routine.
94*
95*  -- Written on 22-October-1986.
96*     Jack Dongarra, Argonne National Lab.
97*     Jeremy Du Croz, Nag Central Office.
98*     Sven Hammarling, Nag Central Office.
99*     Richard Hanson, Sandia National Labs.
100*
101*
102*     .. Parameters ..
103      COMPLEX ONE
104      PARAMETER (ONE= (1.0E+0,0.0E+0))
105      COMPLEX ZERO
106      PARAMETER (ZERO= (0.0E+0,0.0E+0))
107*     ..
108*     .. Local Scalars ..
109      COMPLEX TEMP1,TEMP2
110      INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
111*     ..
112*     .. External Functions ..
113      LOGICAL LSAME
114      EXTERNAL LSAME
115*     ..
116*     .. External Subroutines ..
117      EXTERNAL XERBLA
118*     ..
119*     .. Intrinsic Functions ..
120      INTRINSIC CONJG,MAX,REAL
121*     ..
122*
123*     Test the input parameters.
124*
125      INFO = 0
126      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
127          INFO = 1
128      ELSE IF (N.LT.0) THEN
129          INFO = 2
130      ELSE IF (LDA.LT.MAX(1,N)) THEN
131          INFO = 5
132      ELSE IF (INCX.EQ.0) THEN
133          INFO = 7
134      ELSE IF (INCY.EQ.0) THEN
135          INFO = 10
136      END IF
137      IF (INFO.NE.0) THEN
138          CALL XERBLA('CHEMV ',INFO)
139          RETURN
140      END IF
141*
142*     Quick return if possible.
143*
144      IF ((N.EQ.0) .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN
145*
146*     Set up the start points in  X  and  Y.
147*
148      IF (INCX.GT.0) THEN
149          KX = 1
150      ELSE
151          KX = 1 - (N-1)*INCX
152      END IF
153      IF (INCY.GT.0) THEN
154          KY = 1
155      ELSE
156          KY = 1 - (N-1)*INCY
157      END IF
158*
159*     Start the operations. In this version the elements of A are
160*     accessed sequentially with one pass through the triangular part
161*     of A.
162*
163*     First form  y := beta*y.
164*
165      IF (BETA.NE.ONE) THEN
166          IF (INCY.EQ.1) THEN
167              IF (BETA.EQ.ZERO) THEN
168                  DO 10 I = 1,N
169                      Y(I) = ZERO
170   10             CONTINUE
171              ELSE
172                  DO 20 I = 1,N
173                      Y(I) = BETA*Y(I)
174   20             CONTINUE
175              END IF
176          ELSE
177              IY = KY
178              IF (BETA.EQ.ZERO) THEN
179                  DO 30 I = 1,N
180                      Y(IY) = ZERO
181                      IY = IY + INCY
182   30             CONTINUE
183              ELSE
184                  DO 40 I = 1,N
185                      Y(IY) = BETA*Y(IY)
186                      IY = IY + INCY
187   40             CONTINUE
188              END IF
189          END IF
190      END IF
191      IF (ALPHA.EQ.ZERO) RETURN
192      IF (LSAME(UPLO,'U')) THEN
193*
194*        Form  y  when A is stored in upper triangle.
195*
196          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
197              DO 60 J = 1,N
198                  TEMP1 = ALPHA*X(J)
199                  TEMP2 = ZERO
200                  DO 50 I = 1,J - 1
201                      Y(I) = Y(I) + TEMP1*A(I,J)
202                      TEMP2 = TEMP2 + CONJG(A(I,J))*X(I)
203   50             CONTINUE
204                  Y(J) = Y(J) + TEMP1*REAL(A(J,J)) + ALPHA*TEMP2
205   60         CONTINUE
206          ELSE
207              JX = KX
208              JY = KY
209              DO 80 J = 1,N
210                  TEMP1 = ALPHA*X(JX)
211                  TEMP2 = ZERO
212                  IX = KX
213                  IY = KY
214                  DO 70 I = 1,J - 1
215                      Y(IY) = Y(IY) + TEMP1*A(I,J)
216                      TEMP2 = TEMP2 + CONJG(A(I,J))*X(IX)
217                      IX = IX + INCX
218                      IY = IY + INCY
219   70             CONTINUE
220                  Y(JY) = Y(JY) + TEMP1*REAL(A(J,J)) + ALPHA*TEMP2
221                  JX = JX + INCX
222                  JY = JY + INCY
223   80         CONTINUE
224          END IF
225      ELSE
226*
227*        Form  y  when A is stored in lower triangle.
228*
229          IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
230              DO 100 J = 1,N
231                  TEMP1 = ALPHA*X(J)
232                  TEMP2 = ZERO
233                  Y(J) = Y(J) + TEMP1*REAL(A(J,J))
234                  DO 90 I = J + 1,N
235                      Y(I) = Y(I) + TEMP1*A(I,J)
236                      TEMP2 = TEMP2 + CONJG(A(I,J))*X(I)
237   90             CONTINUE
238                  Y(J) = Y(J) + ALPHA*TEMP2
239  100         CONTINUE
240          ELSE
241              JX = KX
242              JY = KY
243              DO 120 J = 1,N
244                  TEMP1 = ALPHA*X(JX)
245                  TEMP2 = ZERO
246                  Y(JY) = Y(JY) + TEMP1*REAL(A(J,J))
247                  IX = JX
248                  IY = JY
249                  DO 110 I = J + 1,N
250                      IX = IX + INCX
251                      IY = IY + INCY
252                      Y(IY) = Y(IY) + TEMP1*A(I,J)
253                      TEMP2 = TEMP2 + CONJG(A(I,J))*X(IX)
254  110             CONTINUE
255                  Y(JY) = Y(JY) + ALPHA*TEMP2
256                  JX = JX + INCX
257                  JY = JY + INCY
258  120         CONTINUE
259          END IF
260      END IF
261*
262      RETURN
263*
264*     End of CHEMV .
265*
266      END
Note: See TracBrowser for help on using the repository browser.