Mercurial > hg > octave-nkf

comparison libcruft/blas/zher.f @ 5164:57077d0ddc8e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2005-02-25 19:55:24 by jwe]
author jwe
date Fri, 25 Feb 2005 19:55:28 +0000
parents
children
comparison
equal deleted inserted replaced
5163:9f3299378193 5164:57077d0ddc8e
1 SUBROUTINE ZHER ( UPLO, N, ALPHA, X, INCX, A, LDA )
2 * .. Scalar Arguments ..
3 DOUBLE PRECISION ALPHA
4 INTEGER INCX, LDA, N
5 CHARACTER*1 UPLO
6 * .. Array Arguments ..
7 COMPLEX*16 A( LDA, * ), X( * )
8 * ..
9 *
10 * Purpose
11 * =======
12 *
13 * ZHER 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.
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 - DOUBLE PRECISION.
42 * On entry, ALPHA specifies the scalar alpha.
43 * Unchanged on exit.
44 *
45 * X - COMPLEX*16 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 * A - COMPLEX*16 array of DIMENSION ( LDA, n ).
57 * Before entry with UPLO = 'U' or 'u', the leading n by n
58 * upper triangular part of the array A must contain the upper
59 * triangular part of the hermitian matrix and the strictly
60 * lower triangular part of A is not referenced. On exit, the
61 * upper triangular part of the array A is overwritten by the
62 * upper triangular part of the updated matrix.
63 * Before entry with UPLO = 'L' or 'l', the leading n by n
64 * lower triangular part of the array A must contain the lower
65 * triangular part of the hermitian matrix and the strictly
66 * upper triangular part of A is not referenced. On exit, the
67 * lower triangular part of the array A is overwritten by the
68 * lower triangular part of the updated matrix.
69 * Note that the imaginary parts of the diagonal elements need
70 * not be set, they are assumed to be zero, and on exit they
71 * are set to zero.
72 *
73 * LDA - INTEGER.
74 * On entry, LDA specifies the first dimension of A as declared
75 * in the calling (sub) program. LDA must be at least
76 * max( 1, n ).
77 * Unchanged on exit.
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 * .. Parameters ..
90 COMPLEX*16 ZERO
91 PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ) )
92 * .. Local Scalars ..
93 COMPLEX*16 TEMP
94 INTEGER I, INFO, IX, J, JX, KX
95 * .. External Functions ..
96 LOGICAL LSAME
97 EXTERNAL LSAME
98 * .. External Subroutines ..
99 EXTERNAL XERBLA
100 * .. Intrinsic Functions ..
101 INTRINSIC DCONJG, MAX, DBLE
102 * ..
103 * .. Executable Statements ..
104 *
105 * Test the input parameters.
106 *
107 INFO = 0
108 IF ( .NOT.LSAME( UPLO, 'U' ).AND.
109 $ .NOT.LSAME( UPLO, 'L' ) )THEN
110 INFO = 1
111 ELSE IF( N.LT.0 )THEN
112 INFO = 2
113 ELSE IF( INCX.EQ.0 )THEN
114 INFO = 5
115 ELSE IF( LDA.LT.MAX( 1, N ) )THEN
116 INFO = 7
117 END IF
118 IF( INFO.NE.0 )THEN
119 CALL XERBLA( 'ZHER ', INFO )
120 RETURN
121 END IF
122 *
123 * Quick return if possible.
124 *
125 IF( ( N.EQ.0 ).OR.( ALPHA.EQ.DBLE( ZERO ) ) )
126 $ 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 A are
137 * accessed sequentially with one pass through the triangular part
138 * of A.
139 *
140 IF( LSAME( UPLO, 'U' ) )THEN
141 *
142 * Form A when A is stored in upper triangle.
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 DO 10, I = 1, J - 1
149 A( I, J ) = A( I, J ) + X( I )*TEMP
150 10 CONTINUE
151 A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( J )*TEMP )
152 ELSE
153 A( J, J ) = DBLE( A( J, J ) )
154 END IF
155 20 CONTINUE
156 ELSE
157 JX = KX
158 DO 40, J = 1, N
159 IF( X( JX ).NE.ZERO )THEN
160 TEMP = ALPHA*DCONJG( X( JX ) )
161 IX = KX
162 DO 30, I = 1, J - 1
163 A( I, J ) = A( I, J ) + X( IX )*TEMP
164 IX = IX + INCX
165 30 CONTINUE
166 A( J, J ) = DBLE( A( J, J ) ) + DBLE( X( JX )*TEMP )
167 ELSE
168 A( J, J ) = DBLE( A( J, J ) )
169 END IF
170 JX = JX + INCX
171 40 CONTINUE
172 END IF
173 ELSE
174 *
175 * Form A when A is stored in lower triangle.
176 *
177 IF( INCX.EQ.1 )THEN
178 DO 60, J = 1, N
179 IF( X( J ).NE.ZERO )THEN
180 TEMP = ALPHA*DCONJG( X( J ) )
181 A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( J ) )
182 DO 50, I = J + 1, N
183 A( I, J ) = A( I, J ) + X( I )*TEMP
184 50 CONTINUE
185 ELSE
186 A( J, J ) = DBLE( A( J, J ) )
187 END IF
188 60 CONTINUE
189 ELSE
190 JX = KX
191 DO 80, J = 1, N
192 IF( X( JX ).NE.ZERO )THEN
193 TEMP = ALPHA*DCONJG( X( JX ) )
194 A( J, J ) = DBLE( A( J, J ) ) + DBLE( TEMP*X( JX ) )
195 IX = JX
196 DO 70, I = J + 1, N
197 IX = IX + INCX
198 A( I, J ) = A( I, J ) + X( IX )*TEMP
199 70 CONTINUE
200 ELSE
201 A( J, J ) = DBLE( A( J, J ) )
202 END IF
203 JX = JX + INCX
204 80 CONTINUE
205 END IF
206 END IF
207 *
208 RETURN
209 *
210 * End of ZHER .
211 *
212 END