← Back to branch summary

~ubuntu-branches/ubuntu/vivid/atlas/vivid

« back to all changes in this revision

Viewing changes to src/blas/f77reference/cher2.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.3 upstream)
  • mto: (2.2.21 experimental)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-b72z8f621tuhbzn0
Tags: upstream-3.10.1
Import upstream version 3.10.1

Show diffs side-by-side

added added

removed removed

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