Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/cher2.f
5191 views
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*1 UPLO
6
* .. Array Arguments ..
7
COMPLEX A( LDA, * ), X( * ), Y( * )
8
* ..
9
*
10
* Purpose
11
* =======
12
*
13
* CHER2 performs the hermitian rank 2 operation
14
*
15
* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
16
*
17
* where alpha is a scalar, x and y are n element vectors and A is an n
18
* by n hermitian matrix.
19
*
20
* Parameters
21
* ==========
22
*
23
* UPLO - CHARACTER*1.
24
* On entry, UPLO specifies whether the upper or lower
25
* triangular part of the array A is to be referenced as
26
* follows:
27
*
28
* UPLO = 'U' or 'u' Only the upper triangular part of A
29
* is to be referenced.
30
*
31
* UPLO = 'L' or 'l' Only the lower triangular part of A
32
* is to be referenced.
33
*
34
* Unchanged on exit.
35
*
36
* N - INTEGER.
37
* On entry, N specifies the order of the matrix A.
38
* N must be at least zero.
39
* Unchanged on exit.
40
*
41
* ALPHA - COMPLEX .
42
* On entry, ALPHA specifies the scalar alpha.
43
* Unchanged on exit.
44
*
45
* X - COMPLEX array of dimension at least
46
* ( 1 + ( n - 1 )*abs( INCX ) ).
47
* Before entry, the incremented array X must contain the n
48
* element vector x.
49
* Unchanged on exit.
50
*
51
* INCX - INTEGER.
52
* On entry, INCX specifies the increment for the elements of
53
* X. INCX must not be zero.
54
* Unchanged on exit.
55
*
56
* Y - COMPLEX array of dimension at least
57
* ( 1 + ( n - 1 )*abs( INCY ) ).
58
* Before entry, the incremented array Y must contain the n
59
* element vector y.
60
* Unchanged on exit.
61
*
62
* INCY - INTEGER.
63
* On entry, INCY specifies the increment for the elements of
64
* Y. INCY must not be zero.
65
* Unchanged on exit.
66
*
67
* A - COMPLEX array of DIMENSION ( LDA, n ).
68
* Before entry with UPLO = 'U' or 'u', the leading n by n
69
* upper triangular part of the array A must contain the upper
70
* triangular part of the hermitian matrix and the strictly
71
* lower triangular part of A is not referenced. On exit, the
72
* upper triangular part of the array A is overwritten by the
73
* upper triangular part of the updated matrix.
74
* Before entry with UPLO = 'L' or 'l', the leading n by n
75
* lower triangular part of the array A must contain the lower
76
* triangular part of the hermitian matrix and the strictly
77
* upper triangular part of A is not referenced. On exit, the
78
* lower triangular part of the array A is overwritten by the
79
* lower triangular part of the updated matrix.
80
* Note that the imaginary parts of the diagonal elements need
81
* not be set, they are assumed to be zero, and on exit they
82
* are set to zero.
83
*
84
* LDA - INTEGER.
85
* On entry, LDA specifies the first dimension of A as declared
86
* in the calling (sub) program. LDA must be at least
87
* max( 1, n ).
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
COMPLEX ZERO
102
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) )
103
* .. Local Scalars ..
104
COMPLEX 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 CONJG, MAX, REAL
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( INCX.EQ.0 )THEN
125
INFO = 5
126
ELSE IF( INCY.EQ.0 )THEN
127
INFO = 7
128
ELSE IF( LDA.LT.MAX( 1, N ) )THEN
129
INFO = 9
130
END IF
131
IF( INFO.NE.0 )THEN
132
CALL XERBLA( 'CHER2 ', INFO )
133
RETURN
134
END IF
135
*
136
* Quick return if possible.
137
*
138
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
139
$ RETURN
140
*
141
* Set up the start points in X and Y if the increments are not both
142
* unity.
143
*
144
IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
145
IF( INCX.GT.0 )THEN
146
KX = 1
147
ELSE
148
KX = 1 - ( N - 1 )*INCX
149
END IF
150
IF( INCY.GT.0 )THEN
151
KY = 1
152
ELSE
153
KY = 1 - ( N - 1 )*INCY
154
END IF
155
JX = KX
156
JY = KY
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
IF( LSAME( UPLO, 'U' ) )THEN
164
*
165
* Form A when A is stored in the upper triangle.
166
*
167
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
168
DO 20, J = 1, N
169
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
170
TEMP1 = ALPHA*CONJG( Y( J ) )
171
TEMP2 = CONJG( ALPHA*X( J ) )
172
DO 10, I = 1, J - 1
173
A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
174
10 CONTINUE
175
A( J, J ) = REAL( A( J, J ) ) +
176
$ REAL( X( J )*TEMP1 + Y( J )*TEMP2 )
177
ELSE
178
A( J, J ) = REAL( A( J, J ) )
179
END IF
180
20 CONTINUE
181
ELSE
182
DO 40, J = 1, N
183
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
184
TEMP1 = ALPHA*CONJG( Y( JY ) )
185
TEMP2 = CONJG( ALPHA*X( JX ) )
186
IX = KX
187
IY = KY
188
DO 30, I = 1, J - 1
189
A( I, J ) = A( I, J ) + X( IX )*TEMP1
190
$ + Y( IY )*TEMP2
191
IX = IX + INCX
192
IY = IY + INCY
193
30 CONTINUE
194
A( J, J ) = REAL( A( J, J ) ) +
195
$ REAL( X( JX )*TEMP1 + Y( JY )*TEMP2 )
196
ELSE
197
A( J, J ) = REAL( A( J, J ) )
198
END IF
199
JX = JX + INCX
200
JY = JY + INCY
201
40 CONTINUE
202
END IF
203
ELSE
204
*
205
* Form A when A is stored in the lower triangle.
206
*
207
IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
208
DO 60, J = 1, N
209
IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
210
TEMP1 = ALPHA*CONJG( Y( J ) )
211
TEMP2 = CONJG( ALPHA*X( J ) )
212
A( J, J ) = REAL( A( J, J ) ) +
213
$ REAL( X( J )*TEMP1 + Y( J )*TEMP2 )
214
DO 50, I = J + 1, N
215
A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
216
50 CONTINUE
217
ELSE
218
A( J, J ) = REAL( A( J, J ) )
219
END IF
220
60 CONTINUE
221
ELSE
222
DO 80, J = 1, N
223
IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
224
TEMP1 = ALPHA*CONJG( Y( JY ) )
225
TEMP2 = CONJG( ALPHA*X( JX ) )
226
A( J, J ) = REAL( A( J, J ) ) +
227
$ REAL( X( JX )*TEMP1 + Y( JY )*TEMP2 )
228
IX = JX
229
IY = JY
230
DO 70, I = J + 1, N
231
IX = IX + INCX
232
IY = IY + INCY
233
A( I, J ) = A( I, J ) + X( IX )*TEMP1
234
$ + 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
250
251