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