Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/blas/chpr.f
5191 views
1
SUBROUTINE CHPR ( UPLO, N, ALPHA, X, INCX, AP )
2
* .. Scalar Arguments ..
3
REAL ALPHA
4
INTEGER INCX, N
5
CHARACTER*1 UPLO
6
* .. Array Arguments ..
7
COMPLEX AP( * ), X( * )
8
* ..
9
*
10
* Purpose
11
* =======
12
*
13
* CHPR performs the hermitian rank 1 operation
14
*
15
* A := alpha*x*conjg( x' ) + A,
16
*
17
* where alpha is a real scalar, x is an n element vector and A is an
18
* n by n hermitian matrix, supplied in packed form.
19
*
20
* Parameters
21
* ==========
22
*
23
* UPLO - CHARACTER*1.
24
* On entry, UPLO specifies whether the upper or lower
25
* triangular part of the matrix A is supplied in the packed
26
* array AP as follows:
27
*
28
* UPLO = 'U' or 'u' The upper triangular part of A is
29
* supplied in AP.
30
*
31
* UPLO = 'L' or 'l' The lower triangular part of A is
32
* supplied in AP.
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 - REAL .
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
* AP - COMPLEX array of DIMENSION at least
57
* ( ( n*( n + 1 ) )/2 ).
58
* Before entry with UPLO = 'U' or 'u', the array AP must
59
* contain the upper triangular part of the hermitian matrix
60
* packed sequentially, column by column, so that AP( 1 )
61
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
62
* and a( 2, 2 ) respectively, and so on. On exit, the array
63
* AP is overwritten by the upper triangular part of the
64
* updated matrix.
65
* Before entry with UPLO = 'L' or 'l', the array AP must
66
* contain the lower triangular part of the hermitian matrix
67
* packed sequentially, column by column, so that AP( 1 )
68
* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
69
* and a( 3, 1 ) respectively, and so on. On exit, the array
70
* AP is overwritten by the lower triangular part of the
71
* updated matrix.
72
* Note that the imaginary parts of the diagonal elements need
73
* not be set, they are assumed to be zero, and on exit they
74
* are set to zero.
75
*
76
*
77
* Level 2 Blas routine.
78
*
79
* -- Written on 22-October-1986.
80
* Jack Dongarra, Argonne National Lab.
81
* Jeremy Du Croz, Nag Central Office.
82
* Sven Hammarling, Nag Central Office.
83
* Richard Hanson, Sandia National Labs.
84
*
85
*
86
* .. Parameters ..
87
COMPLEX ZERO
88
PARAMETER ( ZERO = ( 0.0E+0, 0.0E+0 ) )
89
* .. Local Scalars ..
90
COMPLEX TEMP
91
INTEGER I, INFO, IX, J, JX, K, KK, KX
92
* .. External Functions ..
93
LOGICAL LSAME
94
EXTERNAL LSAME
95
* .. External Subroutines ..
96
EXTERNAL XERBLA
97
* .. Intrinsic Functions ..
98
INTRINSIC CONJG, REAL
99
* ..
100
* .. Executable Statements ..
101
*
102
* Test the input parameters.
103
*
104
INFO = 0
105
IF ( .NOT.LSAME( UPLO, 'U' ).AND.
106
$ .NOT.LSAME( UPLO, 'L' ) )THEN
107
INFO = 1
108
ELSE IF( N.LT.0 )THEN
109
INFO = 2
110
ELSE IF( INCX.EQ.0 )THEN
111
INFO = 5
112
END IF
113
IF( INFO.NE.0 )THEN
114
CALL XERBLA( 'CHPR ', INFO )
115
RETURN
116
END IF
117
*
118
* Quick return if possible.
119
*
120
IF( ( N.EQ.0 ).OR.( ALPHA.EQ.REAL( ZERO ) ) )
121
$ RETURN
122
*
123
* Set the start point in X if the increment is not unity.
124
*
125
IF( INCX.LE.0 )THEN
126
KX = 1 - ( N - 1 )*INCX
127
ELSE IF( INCX.NE.1 )THEN
128
KX = 1
129
END IF
130
*
131
* Start the operations. In this version the elements of the array AP
132
* are accessed sequentially with one pass through AP.
133
*
134
KK = 1
135
IF( LSAME( UPLO, 'U' ) )THEN
136
*
137
* Form A when upper triangle is stored in AP.
138
*
139
IF( INCX.EQ.1 )THEN
140
DO 20, J = 1, N
141
IF( X( J ).NE.ZERO )THEN
142
TEMP = ALPHA*CONJG( X( J ) )
143
K = KK
144
DO 10, I = 1, J - 1
145
AP( K ) = AP( K ) + X( I )*TEMP
146
K = K + 1
147
10 CONTINUE
148
AP( KK + J - 1 ) = REAL( AP( KK + J - 1 ) )
149
$ + REAL( X( J )*TEMP )
150
ELSE
151
AP( KK + J - 1 ) = REAL( AP( KK + J - 1 ) )
152
END IF
153
KK = KK + J
154
20 CONTINUE
155
ELSE
156
JX = KX
157
DO 40, J = 1, N
158
IF( X( JX ).NE.ZERO )THEN
159
TEMP = ALPHA*CONJG( X( JX ) )
160
IX = KX
161
DO 30, K = KK, KK + J - 2
162
AP( K ) = AP( K ) + X( IX )*TEMP
163
IX = IX + INCX
164
30 CONTINUE
165
AP( KK + J - 1 ) = REAL( AP( KK + J - 1 ) )
166
$ + REAL( X( JX )*TEMP )
167
ELSE
168
AP( KK + J - 1 ) = REAL( AP( KK + J - 1 ) )
169
END IF
170
JX = JX + INCX
171
KK = KK + J
172
40 CONTINUE
173
END IF
174
ELSE
175
*
176
* Form A when lower triangle is stored in AP.
177
*
178
IF( INCX.EQ.1 )THEN
179
DO 60, J = 1, N
180
IF( X( J ).NE.ZERO )THEN
181
TEMP = ALPHA*CONJG( X( J ) )
182
AP( KK ) = REAL( AP( KK ) ) + REAL( TEMP*X( J ) )
183
K = KK + 1
184
DO 50, I = J + 1, N
185
AP( K ) = AP( K ) + X( I )*TEMP
186
K = K + 1
187
50 CONTINUE
188
ELSE
189
AP( KK ) = REAL( AP( KK ) )
190
END IF
191
KK = KK + N - J + 1
192
60 CONTINUE
193
ELSE
194
JX = KX
195
DO 80, J = 1, N
196
IF( X( JX ).NE.ZERO )THEN
197
TEMP = ALPHA*CONJG( X( JX ) )
198
AP( KK ) = REAL( AP( KK ) ) + REAL( TEMP*X( JX ) )
199
IX = JX
200
DO 70, K = KK + 1, KK + N - J
201
IX = IX + INCX
202
AP( K ) = AP( K ) + X( IX )*TEMP
203
70 CONTINUE
204
ELSE
205
AP( KK ) = REAL( AP( KK ) )
206
END IF
207
JX = JX + INCX
208
KK = KK + N - J + 1
209
80 CONTINUE
210
END IF
211
END IF
212
*
213
RETURN
214
*
215
* End of CHPR .
216
*
217
END
218
219