1 /*
3 BLIS
4 An object-based framework for developing high-performance BLAS-like
5 libraries.
7 Copyright (C) 2014, The University of Texas at Austin
9 Redistribution and use in source and binary forms, with or without
10 modification, are permitted provided that the following conditions are
11 met:
12 - Redistributions of source code must retain the above copyright
13 notice, this list of conditions and the following disclaimer.
14 - Redistributions in binary form must reproduce the above copyright
15 notice, this list of conditions and the following disclaimer in the
16 documentation and/or other materials provided with the distribution.
17 - Neither the name of The University of Texas at Austin nor the names
18 of its contributors may be used to endorse or promote products
19 derived from this software without specific prior written permission.
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25 HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 */
35 #include "blis.h"
37 #ifdef BLIS_ENABLE_BLAS2BLIS
39 /* chpr2.f -- translated by f2c (version 19991025).
40 You must link the resulting object file with the libraries:
41 -lf2c -lm (in that order)
42 */
44 /* Subroutine */ int PASTEF77(c,hpr2)(character *uplo, integer *n, singlecomplex *alpha, singlecomplex *x, integer *incx, singlecomplex *y, integer *incy, singlecomplex *ap)
45 {
46 /* System generated locals */
47 integer i__1, i__2, i__3, i__4, i__5, i__6;
48 real r__1;
49 singlecomplex q__1, q__2, q__3, q__4;
51 /* Builtin functions */
52 void bla_r_cnjg(singlecomplex *, singlecomplex *);
54 /* Local variables */
55 integer info;
56 singlecomplex temp1, temp2;
57 integer i__, j, k;
58 extern logical PASTEF770(lsame)(character *, character *, ftnlen, ftnlen);
59 integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
60 extern /* Subroutine */ int PASTEF770(xerbla)(character *, integer *, ftnlen);
62 /* .. Scalar Arguments .. */
63 /* .. Array Arguments .. */
64 /* .. */
66 /* Purpose */
67 /* ======= */
69 /* CHPR2 performs the hermitian rank 2 operation */
71 /* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
73 /* where alpha is a scalar, x and y are n element vectors and A is an */
74 /* n by n hermitian matrix, supplied in packed form. */
76 /* Parameters */
77 /* ========== */
79 /* UPLO - CHARACTER*1. */
80 /* On entry, UPLO specifies whether the upper or lower */
81 /* triangular part of the matrix A is supplied in the packed */
82 /* array AP as follows: */
84 /* UPLO = 'U' or 'u' The upper triangular part of A is */
85 /* supplied in AP. */
87 /* UPLO = 'L' or 'l' The lower triangular part of A is */
88 /* supplied in AP. */
90 /* Unchanged on exit. */
92 /* N - INTEGER. */
93 /* On entry, N specifies the order of the matrix A. */
94 /* N must be at least zero. */
95 /* Unchanged on exit. */
97 /* ALPHA - COMPLEX . */
98 /* On entry, ALPHA specifies the scalar alpha. */
99 /* Unchanged on exit. */
101 /* X - COMPLEX array of dimension at least */
102 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
103 /* Before entry, the incremented array X must contain the n */
104 /* element vector x. */
105 /* Unchanged on exit. */
107 /* INCX - INTEGER. */
108 /* On entry, INCX specifies the increment for the elements of */
109 /* X. INCX must not be zero. */
110 /* Unchanged on exit. */
112 /* Y - COMPLEX array of dimension at least */
113 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
114 /* Before entry, the incremented array Y must contain the n */
115 /* element vector y. */
116 /* Unchanged on exit. */
118 /* INCY - INTEGER. */
119 /* On entry, INCY specifies the increment for the elements of */
120 /* Y. INCY must not be zero. */
121 /* Unchanged on exit. */
123 /* AP - COMPLEX array of DIMENSION at least */
124 /* ( ( n*( n + 1 ) )/2 ). */
125 /* Before entry with UPLO = 'U' or 'u', the array AP must */
126 /* contain the upper triangular part of the hermitian matrix */
127 /* packed sequentially, column by column, so that AP( 1 ) */
128 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
129 /* and a( 2, 2 ) respectively, and so on. On exit, the array */
130 /* AP is overwritten by the upper triangular part of the */
131 /* updated matrix. */
132 /* Before entry with UPLO = 'L' or 'l', the array AP must */
133 /* contain the lower triangular part of the hermitian matrix */
134 /* packed sequentially, column by column, so that AP( 1 ) */
135 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
136 /* and a( 3, 1 ) respectively, and so on. On exit, the array */
137 /* AP is overwritten by the lower triangular part of the */
138 /* updated matrix. */
139 /* Note that the imaginary parts of the diagonal elements need */
140 /* not be set, they are assumed to be zero, and on exit they */
141 /* are set to zero. */
144 /* Level 2 Blas routine. */
146 /* -- Written on 22-October-1986. */
147 /* Jack Dongarra, Argonne National Lab. */
148 /* Jeremy Du Croz, Nag Central Office. */
149 /* Sven Hammarling, Nag Central Office. */
150 /* Richard Hanson, Sandia National Labs. */
153 /* .. Parameters .. */
154 /* .. Local Scalars .. */
155 /* .. External Functions .. */
156 /* .. External Subroutines .. */
157 /* .. Intrinsic Functions .. */
158 /* .. */
159 /* .. Executable Statements .. */
161 /* Test the input parameters. */
163 /* Parameter adjustments */
164 --ap;
165 --y;
166 --x;
168 /* Function Body */
169 info = 0;
170 if (! PASTEF770(lsame)(uplo, "U", (ftnlen)1, (ftnlen)1) && ! PASTEF770(lsame)(uplo, "L", (
171 ftnlen)1, (ftnlen)1)) {
172 info = 1;
173 } else if (*n < 0) {
174 info = 2;
175 } else if (*incx == 0) {
176 info = 5;
177 } else if (*incy == 0) {
178 info = 7;
179 }
180 if (info != 0) {
181 PASTEF770(xerbla)("CHPR2 ", &info, (ftnlen)6);
182 return 0;
183 }
185 /* Quick return if possible. */
187 if (*n == 0 || (bli_creal(*alpha) == 0.f && bli_cimag(*alpha) == 0.f)) {
188 return 0;
189 }
191 /* Set up the start points in X and Y if the increments are not both */
192 /* unity. */
194 if (*incx != 1 || *incy != 1) {
195 if (*incx > 0) {
196 kx = 1;
197 } else {
198 kx = 1 - (*n - 1) * *incx;
199 }
200 if (*incy > 0) {
201 ky = 1;
202 } else {
203 ky = 1 - (*n - 1) * *incy;
204 }
205 jx = kx;
206 jy = ky;
207 }
209 /* Start the operations. In this version the elements of the array AP */
210 /* are accessed sequentially with one pass through AP. */
212 kk = 1;
213 if (PASTEF770(lsame)(uplo, "U", (ftnlen)1, (ftnlen)1)) {
215 /* Form A when upper triangle is stored in AP. */
217 if (*incx == 1 && *incy == 1) {
218 i__1 = *n;
219 for (j = 1; j <= i__1; ++j) {
220 i__2 = j;
221 i__3 = j;
222 if (bli_creal(x[i__2]) != 0.f || bli_cimag(x[i__2]) != 0.f || (bli_creal(y[i__3]) != 0.f
223 || bli_cimag(y[i__3]) != 0.f)) {
224 bla_r_cnjg(&q__2, &y[j]);
225 bli_csets( (bli_creal(*alpha) * bli_creal(q__2) - bli_cimag(*alpha) * bli_cimag(q__2)), (bli_creal(*alpha) * bli_cimag(q__2) + bli_cimag(*alpha) * bli_creal(q__2)), q__1 );
226 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp1 );
227 i__2 = j;
228 bli_csets( (bli_creal(*alpha) * bli_creal(x[i__2]) - bli_cimag(*alpha) * bli_cimag(x[i__2])), (bli_creal(*alpha) * bli_cimag(x[i__2]) + bli_cimag(*alpha) * bli_creal(x[i__2])), q__2 );
229 bla_r_cnjg(&q__1, &q__2);
230 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp2 );
231 k = kk;
232 i__2 = j - 1;
233 for (i__ = 1; i__ <= i__2; ++i__) {
234 i__3 = k;
235 i__4 = k;
236 i__5 = i__;
237 bli_csets( (bli_creal(x[i__5]) * bli_creal(temp1) - bli_cimag(x[i__5]) * bli_cimag(temp1)), (bli_creal(x[i__5]) * bli_cimag(temp1) + bli_cimag(x[i__5]) * bli_creal(temp1)), q__3 );
238 bli_csets( (bli_creal(ap[i__4]) + bli_creal(q__3)), (bli_cimag(ap[i__4]) + bli_cimag(q__3)), q__2 );
239 i__6 = i__;
240 bli_csets( (bli_creal(y[i__6]) * bli_creal(temp2) - bli_cimag(y[i__6]) * bli_cimag(temp2)), (bli_creal(y[i__6]) * bli_cimag(temp2) + bli_cimag(y[i__6]) * bli_creal(temp2)), q__4 );
241 bli_csets( (bli_creal(q__2) + bli_creal(q__4)), (bli_cimag(q__2) + bli_cimag(q__4)), q__1 );
242 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), ap[i__3] );
243 ++k;
244 /* L10: */
245 }
246 i__2 = kk + j - 1;
247 i__3 = kk + j - 1;
248 i__4 = j;
249 bli_csets( (bli_creal(x[i__4]) * bli_creal(temp1) - bli_cimag(x[i__4]) * bli_cimag(temp1)), (bli_creal(x[i__4]) * bli_cimag(temp1) + bli_cimag(x[i__4]) * bli_creal(temp1)), q__2 );
250 i__5 = j;
251 bli_csets( (bli_creal(y[i__5]) * bli_creal(temp2) - bli_cimag(y[i__5]) * bli_cimag(temp2)), (bli_creal(y[i__5]) * bli_cimag(temp2) + bli_cimag(y[i__5]) * bli_creal(temp2)), q__3 );
252 bli_csets( (bli_creal(q__2) + bli_creal(q__3)), (bli_cimag(q__2) + bli_cimag(q__3)), q__1 );
253 r__1 = bli_creal(ap[i__3]) + bli_creal(q__1);
254 bli_csets( (r__1), (0.f), ap[i__2] );
255 } else {
256 i__2 = kk + j - 1;
257 i__3 = kk + j - 1;
258 r__1 = bli_creal(ap[i__3]);
259 bli_csets( (r__1), (0.f), ap[i__2] );
260 }
261 kk += j;
262 /* L20: */
263 }
264 } else {
265 i__1 = *n;
266 for (j = 1; j <= i__1; ++j) {
267 i__2 = jx;
268 i__3 = jy;
269 if (bli_creal(x[i__2]) != 0.f || bli_cimag(x[i__2]) != 0.f || (bli_creal(y[i__3]) != 0.f
270 || bli_cimag(y[i__3]) != 0.f)) {
271 bla_r_cnjg(&q__2, &y[jy]);
272 bli_csets( (bli_creal(*alpha) * bli_creal(q__2) - bli_cimag(*alpha) * bli_cimag(q__2)), (bli_creal(*alpha) * bli_cimag(q__2) + bli_cimag(*alpha) * bli_creal(q__2)), q__1 );
273 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp1 );
274 i__2 = jx;
275 bli_csets( (bli_creal(*alpha) * bli_creal(x[i__2]) - bli_cimag(*alpha) * bli_cimag(x[i__2])), (bli_creal(*alpha) * bli_cimag(x[i__2]) + bli_cimag(*alpha) * bli_creal(x[i__2])), q__2 );
276 bla_r_cnjg(&q__1, &q__2);
277 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp2 );
278 ix = kx;
279 iy = ky;
280 i__2 = kk + j - 2;
281 for (k = kk; k <= i__2; ++k) {
282 i__3 = k;
283 i__4 = k;
284 i__5 = ix;
285 bli_csets( (bli_creal(x[i__5]) * bli_creal(temp1) - bli_cimag(x[i__5]) * bli_cimag(temp1)), (bli_creal(x[i__5]) * bli_cimag(temp1) + bli_cimag(x[i__5]) * bli_creal(temp1)), q__3 );
286 bli_csets( (bli_creal(ap[i__4]) + bli_creal(q__3)), (bli_cimag(ap[i__4]) + bli_cimag(q__3)), q__2 );
287 i__6 = iy;
288 bli_csets( (bli_creal(y[i__6]) * bli_creal(temp2) - bli_cimag(y[i__6]) * bli_cimag(temp2)), (bli_creal(y[i__6]) * bli_cimag(temp2) + bli_cimag(y[i__6]) * bli_creal(temp2)), q__4 );
289 bli_csets( (bli_creal(q__2) + bli_creal(q__4)), (bli_cimag(q__2) + bli_cimag(q__4)), q__1 );
290 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), ap[i__3] );
291 ix += *incx;
292 iy += *incy;
293 /* L30: */
294 }
295 i__2 = kk + j - 1;
296 i__3 = kk + j - 1;
297 i__4 = jx;
298 bli_csets( (bli_creal(x[i__4]) * bli_creal(temp1) - bli_cimag(x[i__4]) * bli_cimag(temp1)), (bli_creal(x[i__4]) * bli_cimag(temp1) + bli_cimag(x[i__4]) * bli_creal(temp1)), q__2 );
299 i__5 = jy;
300 bli_csets( (bli_creal(y[i__5]) * bli_creal(temp2) - bli_cimag(y[i__5]) * bli_cimag(temp2)), (bli_creal(y[i__5]) * bli_cimag(temp2) + bli_cimag(y[i__5]) * bli_creal(temp2)), q__3 );
301 bli_csets( (bli_creal(q__2) + bli_creal(q__3)), (bli_cimag(q__2) + bli_cimag(q__3)), q__1 );
302 r__1 = bli_creal(ap[i__3]) + bli_creal(q__1);
303 bli_csets( (r__1), (0.f), ap[i__2] );
304 } else {
305 i__2 = kk + j - 1;
306 i__3 = kk + j - 1;
307 r__1 = bli_creal(ap[i__3]);
308 bli_csets( (r__1), (0.f), ap[i__2] );
309 }
310 jx += *incx;
311 jy += *incy;
312 kk += j;
313 /* L40: */
314 }
315 }
316 } else {
318 /* Form A when lower triangle is stored in AP. */
320 if (*incx == 1 && *incy == 1) {
321 i__1 = *n;
322 for (j = 1; j <= i__1; ++j) {
323 i__2 = j;
324 i__3 = j;
325 if (bli_creal(x[i__2]) != 0.f || bli_cimag(x[i__2]) != 0.f || (bli_creal(y[i__3]) != 0.f
326 || bli_cimag(y[i__3]) != 0.f)) {
327 bla_r_cnjg(&q__2, &y[j]);
328 bli_csets( (bli_creal(*alpha) * bli_creal(q__2) - bli_cimag(*alpha) * bli_cimag(q__2)), (bli_creal(*alpha) * bli_cimag(q__2) + bli_cimag(*alpha) * bli_creal(q__2)), q__1 );
329 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp1 );
330 i__2 = j;
331 bli_csets( (bli_creal(*alpha) * bli_creal(x[i__2]) - bli_cimag(*alpha) * bli_cimag(x[i__2])), (bli_creal(*alpha) * bli_cimag(x[i__2]) + bli_cimag(*alpha) * bli_creal(x[i__2])), q__2 );
332 bla_r_cnjg(&q__1, &q__2);
333 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp2 );
334 i__2 = kk;
335 i__3 = kk;
336 i__4 = j;
337 bli_csets( (bli_creal(x[i__4]) * bli_creal(temp1) - bli_cimag(x[i__4]) * bli_cimag(temp1)), (bli_creal(x[i__4]) * bli_cimag(temp1) + bli_cimag(x[i__4]) * bli_creal(temp1)), q__2 );
338 i__5 = j;
339 bli_csets( (bli_creal(y[i__5]) * bli_creal(temp2) - bli_cimag(y[i__5]) * bli_cimag(temp2)), (bli_creal(y[i__5]) * bli_cimag(temp2) + bli_cimag(y[i__5]) * bli_creal(temp2)), q__3 );
340 bli_csets( (bli_creal(q__2) + bli_creal(q__3)), (bli_cimag(q__2) + bli_cimag(q__3)), q__1 );
341 r__1 = bli_creal(ap[i__3]) + bli_creal(q__1);
342 bli_csets( (r__1), (0.f), ap[i__2] );
343 k = kk + 1;
344 i__2 = *n;
345 for (i__ = j + 1; i__ <= i__2; ++i__) {
346 i__3 = k;
347 i__4 = k;
348 i__5 = i__;
349 bli_csets( (bli_creal(x[i__5]) * bli_creal(temp1) - bli_cimag(x[i__5]) * bli_cimag(temp1)), (bli_creal(x[i__5]) * bli_cimag(temp1) + bli_cimag(x[i__5]) * bli_creal(temp1)), q__3 );
350 bli_csets( (bli_creal(ap[i__4]) + bli_creal(q__3)), (bli_cimag(ap[i__4]) + bli_cimag(q__3)), q__2 );
351 i__6 = i__;
352 bli_csets( (bli_creal(y[i__6]) * bli_creal(temp2) - bli_cimag(y[i__6]) * bli_cimag(temp2)), (bli_creal(y[i__6]) * bli_cimag(temp2) + bli_cimag(y[i__6]) * bli_creal(temp2)), q__4 );
353 bli_csets( (bli_creal(q__2) + bli_creal(q__4)), (bli_cimag(q__2) + bli_cimag(q__4)), q__1 );
354 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), ap[i__3] );
355 ++k;
356 /* L50: */
357 }
358 } else {
359 i__2 = kk;
360 i__3 = kk;
361 r__1 = bli_creal(ap[i__3]);
362 bli_csets( (r__1), (0.f), ap[i__2] );
363 }
364 kk = kk + *n - j + 1;
365 /* L60: */
366 }
367 } else {
368 i__1 = *n;
369 for (j = 1; j <= i__1; ++j) {
370 i__2 = jx;
371 i__3 = jy;
372 if (bli_creal(x[i__2]) != 0.f || bli_cimag(x[i__2]) != 0.f || (bli_creal(y[i__3]) != 0.f
373 || bli_cimag(y[i__3]) != 0.f)) {
374 bla_r_cnjg(&q__2, &y[jy]);
375 bli_csets( (bli_creal(*alpha) * bli_creal(q__2) - bli_cimag(*alpha) * bli_cimag(q__2)), (bli_creal(*alpha) * bli_cimag(q__2) + bli_cimag(*alpha) * bli_creal(q__2)), q__1 );
376 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp1 );
377 i__2 = jx;
378 bli_csets( (bli_creal(*alpha) * bli_creal(x[i__2]) - bli_cimag(*alpha) * bli_cimag(x[i__2])), (bli_creal(*alpha) * bli_cimag(x[i__2]) + bli_cimag(*alpha) * bli_creal(x[i__2])), q__2 );
379 bla_r_cnjg(&q__1, &q__2);
380 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), temp2 );
381 i__2 = kk;
382 i__3 = kk;
383 i__4 = jx;
384 bli_csets( (bli_creal(x[i__4]) * bli_creal(temp1) - bli_cimag(x[i__4]) * bli_cimag(temp1)), (bli_creal(x[i__4]) * bli_cimag(temp1) + bli_cimag(x[i__4]) * bli_creal(temp1)), q__2 );
385 i__5 = jy;
386 bli_csets( (bli_creal(y[i__5]) * bli_creal(temp2) - bli_cimag(y[i__5]) * bli_cimag(temp2)), (bli_creal(y[i__5]) * bli_cimag(temp2) + bli_cimag(y[i__5]) * bli_creal(temp2)), q__3 );
387 bli_csets( (bli_creal(q__2) + bli_creal(q__3)), (bli_cimag(q__2) + bli_cimag(q__3)), q__1 );
388 r__1 = bli_creal(ap[i__3]) + bli_creal(q__1);
389 bli_csets( (r__1), (0.f), ap[i__2] );
390 ix = jx;
391 iy = jy;
392 i__2 = kk + *n - j;
393 for (k = kk + 1; k <= i__2; ++k) {
394 ix += *incx;
395 iy += *incy;
396 i__3 = k;
397 i__4 = k;
398 i__5 = ix;
399 bli_csets( (bli_creal(x[i__5]) * bli_creal(temp1) - bli_cimag(x[i__5]) * bli_cimag(temp1)), (bli_creal(x[i__5]) * bli_cimag(temp1) + bli_cimag(x[i__5]) * bli_creal(temp1)), q__3 );
400 bli_csets( (bli_creal(ap[i__4]) + bli_creal(q__3)), (bli_cimag(ap[i__4]) + bli_cimag(q__3)), q__2 );
401 i__6 = iy;
402 bli_csets( (bli_creal(y[i__6]) * bli_creal(temp2) - bli_cimag(y[i__6]) * bli_cimag(temp2)), (bli_creal(y[i__6]) * bli_cimag(temp2) + bli_cimag(y[i__6]) * bli_creal(temp2)), q__4 );
403 bli_csets( (bli_creal(q__2) + bli_creal(q__4)), (bli_cimag(q__2) + bli_cimag(q__4)), q__1 );
404 bli_csets( (bli_creal(q__1)), (bli_cimag(q__1)), ap[i__3] );
405 /* L70: */
406 }
407 } else {
408 i__2 = kk;
409 i__3 = kk;
410 r__1 = bli_creal(ap[i__3]);
411 bli_csets( (r__1), (0.f), ap[i__2] );
412 }
413 jx += *incx;
414 jy += *incy;
415 kk = kk + *n - j + 1;
416 /* L80: */
417 }
418 }
419 }
421 return 0;
423 /* End of CHPR2 . */
425 } /* chpr2_ */
427 /* zhpr2.f -- translated by f2c (version 19991025).
428 You must link the resulting object file with the libraries:
429 -lf2c -lm (in that order)
430 */
432 /* Subroutine */ int PASTEF77(z,hpr2)(character *uplo, integer *n, doublecomplex *alpha, doublecomplex *x, integer *incx, doublecomplex *y, integer *incy, doublecomplex *ap)
433 {
434 /* System generated locals */
435 integer i__1, i__2, i__3, i__4, i__5, i__6;
436 doublereal d__1;
437 doublecomplex z__1, z__2, z__3, z__4;
439 /* Builtin functions */
440 void bla_d_cnjg(doublecomplex *, doublecomplex *);
442 /* Local variables */
443 integer info;
444 doublecomplex temp1, temp2;
445 integer i__, j, k;
446 extern logical PASTEF770(lsame)(character *, character *, ftnlen, ftnlen);
447 integer kk, ix, iy, jx = 0, jy = 0, kx = 0, ky = 0;
448 extern /* Subroutine */ int PASTEF770(xerbla)(character *, integer *, ftnlen);
450 /* .. Scalar Arguments .. */
451 /* .. Array Arguments .. */
452 /* .. */
454 /* Purpose */
455 /* ======= */
457 /* ZHPR2 performs the hermitian rank 2 operation */
459 /* A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A, */
461 /* where alpha is a scalar, x and y are n element vectors and A is an */
462 /* n by n hermitian matrix, supplied in packed form. */
464 /* Parameters */
465 /* ========== */
467 /* UPLO - CHARACTER*1. */
468 /* On entry, UPLO specifies whether the upper or lower */
469 /* triangular part of the matrix A is supplied in the packed */
470 /* array AP as follows: */
472 /* UPLO = 'U' or 'u' The upper triangular part of A is */
473 /* supplied in AP. */
475 /* UPLO = 'L' or 'l' The lower triangular part of A is */
476 /* supplied in AP. */
478 /* Unchanged on exit. */
480 /* N - INTEGER. */
481 /* On entry, N specifies the order of the matrix A. */
482 /* N must be at least zero. */
483 /* Unchanged on exit. */
485 /* ALPHA - COMPLEX*16 . */
486 /* On entry, ALPHA specifies the scalar alpha. */
487 /* Unchanged on exit. */
489 /* X - COMPLEX*16 array of dimension at least */
490 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
491 /* Before entry, the incremented array X must contain the n */
492 /* element vector x. */
493 /* Unchanged on exit. */
495 /* INCX - INTEGER. */
496 /* On entry, INCX specifies the increment for the elements of */
497 /* X. INCX must not be zero. */
498 /* Unchanged on exit. */
500 /* Y - COMPLEX*16 array of dimension at least */
501 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
502 /* Before entry, the incremented array Y must contain the n */
503 /* element vector y. */
504 /* Unchanged on exit. */
506 /* INCY - INTEGER. */
507 /* On entry, INCY specifies the increment for the elements of */
508 /* Y. INCY must not be zero. */
509 /* Unchanged on exit. */
511 /* AP - COMPLEX*16 array of DIMENSION at least */
512 /* ( ( n*( n + 1 ) )/2 ). */
513 /* Before entry with UPLO = 'U' or 'u', the array AP must */
514 /* contain the upper triangular part of the hermitian matrix */
515 /* packed sequentially, column by column, so that AP( 1 ) */
516 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
517 /* and a( 2, 2 ) respectively, and so on. On exit, the array */
518 /* AP is overwritten by the upper triangular part of the */
519 /* updated matrix. */
520 /* Before entry with UPLO = 'L' or 'l', the array AP must */
521 /* contain the lower triangular part of the hermitian matrix */
522 /* packed sequentially, column by column, so that AP( 1 ) */
523 /* contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
524 /* and a( 3, 1 ) respectively, and so on. On exit, the array */
525 /* AP is overwritten by the lower triangular part of the */
526 /* updated matrix. */
527 /* Note that the imaginary parts of the diagonal elements need */
528 /* not be set, they are assumed to be zero, and on exit they */
529 /* are set to zero. */
532 /* Level 2 Blas routine. */
534 /* -- Written on 22-October-1986. */
535 /* Jack Dongarra, Argonne National Lab. */
536 /* Jeremy Du Croz, Nag Central Office. */
537 /* Sven Hammarling, Nag Central Office. */
538 /* Richard Hanson, Sandia National Labs. */
541 /* .. Parameters .. */
542 /* .. Local Scalars .. */
543 /* .. External Functions .. */
544 /* .. External Subroutines .. */
545 /* .. Intrinsic Functions .. */
546 /* .. */
547 /* .. Executable Statements .. */
549 /* Test the input parameters. */
551 /* Parameter adjustments */
552 --ap;
553 --y;
554 --x;
556 /* Function Body */
557 info = 0;
558 if (! PASTEF770(lsame)(uplo, "U", (ftnlen)1, (ftnlen)1) && ! PASTEF770(lsame)(uplo, "L", (
559 ftnlen)1, (ftnlen)1)) {
560 info = 1;
561 } else if (*n < 0) {
562 info = 2;
563 } else if (*incx == 0) {
564 info = 5;
565 } else if (*incy == 0) {
566 info = 7;
567 }
568 if (info != 0) {
569 PASTEF770(xerbla)("ZHPR2 ", &info, (ftnlen)6);
570 return 0;
571 }
573 /* Quick return if possible. */
575 if (*n == 0 || (bli_zreal(*alpha) == 0. && bli_zimag(*alpha) == 0.)) {
576 return 0;
577 }
579 /* Set up the start points in X and Y if the increments are not both */
580 /* unity. */
582 if (*incx != 1 || *incy != 1) {
583 if (*incx > 0) {
584 kx = 1;
585 } else {
586 kx = 1 - (*n - 1) * *incx;
587 }
588 if (*incy > 0) {
589 ky = 1;
590 } else {
591 ky = 1 - (*n - 1) * *incy;
592 }
593 jx = kx;
594 jy = ky;
595 }
597 /* Start the operations. In this version the elements of the array AP */
598 /* are accessed sequentially with one pass through AP. */
600 kk = 1;
601 if (PASTEF770(lsame)(uplo, "U", (ftnlen)1, (ftnlen)1)) {
603 /* Form A when upper triangle is stored in AP. */
605 if (*incx == 1 && *incy == 1) {
606 i__1 = *n;
607 for (j = 1; j <= i__1; ++j) {
608 i__2 = j;
609 i__3 = j;
610 if (bli_zreal(x[i__2]) != 0. || bli_zimag(x[i__2]) != 0. || (bli_zreal(y[i__3]) != 0. ||
611 bli_zimag(y[i__3]) != 0.)) {
612 bla_d_cnjg(&z__2, &y[j]);
613 bli_zsets( (bli_zreal(*alpha) * bli_zreal(z__2) - bli_zimag(*alpha) * bli_zimag(z__2)), (bli_zreal(*alpha) * bli_zimag(z__2) + bli_zimag(*alpha) * bli_zreal(z__2)), z__1 );
614 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp1 );
615 i__2 = j;
616 bli_zsets( (bli_zreal(*alpha) * bli_zreal(x[i__2]) - bli_zimag(*alpha) * bli_zimag(x[i__2])), (bli_zreal(*alpha) * bli_zimag(x[i__2]) + bli_zimag(*alpha) * bli_zreal(x[i__2])), z__2 );
617 bla_d_cnjg(&z__1, &z__2);
618 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp2 );
619 k = kk;
620 i__2 = j - 1;
621 for (i__ = 1; i__ <= i__2; ++i__) {
622 i__3 = k;
623 i__4 = k;
624 i__5 = i__;
625 bli_zsets( (bli_zreal(x[i__5]) * bli_zreal(temp1) - bli_zimag(x[i__5]) * bli_zimag(temp1)), (bli_zreal(x[i__5]) * bli_zimag(temp1) + bli_zimag(x[i__5]) * bli_zreal(temp1)), z__3 );
626 bli_zsets( (bli_zreal(ap[i__4]) + bli_zreal(z__3)), (bli_zimag(ap[i__4]) + bli_zimag(z__3)), z__2 );
627 i__6 = i__;
628 bli_zsets( (bli_zreal(y[i__6]) * bli_zreal(temp2) - bli_zimag(y[i__6]) * bli_zimag(temp2)), (bli_zreal(y[i__6]) * bli_zimag(temp2) + bli_zimag(y[i__6]) * bli_zreal(temp2)), z__4 );
629 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__4)), (bli_zimag(z__2) + bli_zimag(z__4)), z__1 );
630 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), ap[i__3] );
631 ++k;
632 /* L10: */
633 }
634 i__2 = kk + j - 1;
635 i__3 = kk + j - 1;
636 i__4 = j;
637 bli_zsets( (bli_zreal(x[i__4]) * bli_zreal(temp1) - bli_zimag(x[i__4]) * bli_zimag(temp1)), (bli_zreal(x[i__4]) * bli_zimag(temp1) + bli_zimag(x[i__4]) * bli_zreal(temp1)), z__2 );
638 i__5 = j;
639 bli_zsets( (bli_zreal(y[i__5]) * bli_zreal(temp2) - bli_zimag(y[i__5]) * bli_zimag(temp2)), (bli_zreal(y[i__5]) * bli_zimag(temp2) + bli_zimag(y[i__5]) * bli_zreal(temp2)), z__3 );
640 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__3)), (bli_zimag(z__2) + bli_zimag(z__3)), z__1 );
641 d__1 = bli_zreal(ap[i__3]) + bli_zreal(z__1);
642 bli_zsets( (d__1), (0.), ap[i__2] );
643 } else {
644 i__2 = kk + j - 1;
645 i__3 = kk + j - 1;
646 d__1 = bli_zreal(ap[i__3]);
647 bli_zsets( (d__1), (0.), ap[i__2] );
648 }
649 kk += j;
650 /* L20: */
651 }
652 } else {
653 i__1 = *n;
654 for (j = 1; j <= i__1; ++j) {
655 i__2 = jx;
656 i__3 = jy;
657 if (bli_zreal(x[i__2]) != 0. || bli_zimag(x[i__2]) != 0. || (bli_zreal(y[i__3]) != 0. ||
658 bli_zimag(y[i__3]) != 0.)) {
659 bla_d_cnjg(&z__2, &y[jy]);
660 bli_zsets( (bli_zreal(*alpha) * bli_zreal(z__2) - bli_zimag(*alpha) * bli_zimag(z__2)), (bli_zreal(*alpha) * bli_zimag(z__2) + bli_zimag(*alpha) * bli_zreal(z__2)), z__1 );
661 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp1 );
662 i__2 = jx;
663 bli_zsets( (bli_zreal(*alpha) * bli_zreal(x[i__2]) - bli_zimag(*alpha) * bli_zimag(x[i__2])), (bli_zreal(*alpha) * bli_zimag(x[i__2]) + bli_zimag(*alpha) * bli_zreal(x[i__2])), z__2 );
664 bla_d_cnjg(&z__1, &z__2);
665 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp2 );
666 ix = kx;
667 iy = ky;
668 i__2 = kk + j - 2;
669 for (k = kk; k <= i__2; ++k) {
670 i__3 = k;
671 i__4 = k;
672 i__5 = ix;
673 bli_zsets( (bli_zreal(x[i__5]) * bli_zreal(temp1) - bli_zimag(x[i__5]) * bli_zimag(temp1)), (bli_zreal(x[i__5]) * bli_zimag(temp1) + bli_zimag(x[i__5]) * bli_zreal(temp1)), z__3 );
674 bli_zsets( (bli_zreal(ap[i__4]) + bli_zreal(z__3)), (bli_zimag(ap[i__4]) + bli_zimag(z__3)), z__2 );
675 i__6 = iy;
676 bli_zsets( (bli_zreal(y[i__6]) * bli_zreal(temp2) - bli_zimag(y[i__6]) * bli_zimag(temp2)), (bli_zreal(y[i__6]) * bli_zimag(temp2) + bli_zimag(y[i__6]) * bli_zreal(temp2)), z__4 );
677 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__4)), (bli_zimag(z__2) + bli_zimag(z__4)), z__1 );
678 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), ap[i__3] );
679 ix += *incx;
680 iy += *incy;
681 /* L30: */
682 }
683 i__2 = kk + j - 1;
684 i__3 = kk + j - 1;
685 i__4 = jx;
686 bli_zsets( (bli_zreal(x[i__4]) * bli_zreal(temp1) - bli_zimag(x[i__4]) * bli_zimag(temp1)), (bli_zreal(x[i__4]) * bli_zimag(temp1) + bli_zimag(x[i__4]) * bli_zreal(temp1)), z__2 );
687 i__5 = jy;
688 bli_zsets( (bli_zreal(y[i__5]) * bli_zreal(temp2) - bli_zimag(y[i__5]) * bli_zimag(temp2)), (bli_zreal(y[i__5]) * bli_zimag(temp2) + bli_zimag(y[i__5]) * bli_zreal(temp2)), z__3 );
689 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__3)), (bli_zimag(z__2) + bli_zimag(z__3)), z__1 );
690 d__1 = bli_zreal(ap[i__3]) + bli_zreal(z__1);
691 bli_zsets( (d__1), (0.), ap[i__2] );
692 } else {
693 i__2 = kk + j - 1;
694 i__3 = kk + j - 1;
695 d__1 = bli_zreal(ap[i__3]);
696 bli_zsets( (d__1), (0.), ap[i__2] );
697 }
698 jx += *incx;
699 jy += *incy;
700 kk += j;
701 /* L40: */
702 }
703 }
704 } else {
706 /* Form A when lower triangle is stored in AP. */
708 if (*incx == 1 && *incy == 1) {
709 i__1 = *n;
710 for (j = 1; j <= i__1; ++j) {
711 i__2 = j;
712 i__3 = j;
713 if (bli_zreal(x[i__2]) != 0. || bli_zimag(x[i__2]) != 0. || (bli_zreal(y[i__3]) != 0. ||
714 bli_zimag(y[i__3]) != 0.)) {
715 bla_d_cnjg(&z__2, &y[j]);
716 bli_zsets( (bli_zreal(*alpha) * bli_zreal(z__2) - bli_zimag(*alpha) * bli_zimag(z__2)), (bli_zreal(*alpha) * bli_zimag(z__2) + bli_zimag(*alpha) * bli_zreal(z__2)), z__1 );
717 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp1 );
718 i__2 = j;
719 bli_zsets( (bli_zreal(*alpha) * bli_zreal(x[i__2]) - bli_zimag(*alpha) * bli_zimag(x[i__2])), (bli_zreal(*alpha) * bli_zimag(x[i__2]) + bli_zimag(*alpha) * bli_zreal(x[i__2])), z__2 );
720 bla_d_cnjg(&z__1, &z__2);
721 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp2 );
722 i__2 = kk;
723 i__3 = kk;
724 i__4 = j;
725 bli_zsets( (bli_zreal(x[i__4]) * bli_zreal(temp1) - bli_zimag(x[i__4]) * bli_zimag(temp1)), (bli_zreal(x[i__4]) * bli_zimag(temp1) + bli_zimag(x[i__4]) * bli_zreal(temp1)), z__2 );
726 i__5 = j;
727 bli_zsets( (bli_zreal(y[i__5]) * bli_zreal(temp2) - bli_zimag(y[i__5]) * bli_zimag(temp2)), (bli_zreal(y[i__5]) * bli_zimag(temp2) + bli_zimag(y[i__5]) * bli_zreal(temp2)), z__3 );
728 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__3)), (bli_zimag(z__2) + bli_zimag(z__3)), z__1 );
729 d__1 = bli_zreal(ap[i__3]) + bli_zreal(z__1);
730 bli_zsets( (d__1), (0.), ap[i__2] );
731 k = kk + 1;
732 i__2 = *n;
733 for (i__ = j + 1; i__ <= i__2; ++i__) {
734 i__3 = k;
735 i__4 = k;
736 i__5 = i__;
737 bli_zsets( (bli_zreal(x[i__5]) * bli_zreal(temp1) - bli_zimag(x[i__5]) * bli_zimag(temp1)), (bli_zreal(x[i__5]) * bli_zimag(temp1) + bli_zimag(x[i__5]) * bli_zreal(temp1)), z__3 );
738 bli_zsets( (bli_zreal(ap[i__4]) + bli_zreal(z__3)), (bli_zimag(ap[i__4]) + bli_zimag(z__3)), z__2 );
739 i__6 = i__;
740 bli_zsets( (bli_zreal(y[i__6]) * bli_zreal(temp2) - bli_zimag(y[i__6]) * bli_zimag(temp2)), (bli_zreal(y[i__6]) * bli_zimag(temp2) + bli_zimag(y[i__6]) * bli_zreal(temp2)), z__4 );
741 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__4)), (bli_zimag(z__2) + bli_zimag(z__4)), z__1 );
742 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), ap[i__3] );
743 ++k;
744 /* L50: */
745 }
746 } else {
747 i__2 = kk;
748 i__3 = kk;
749 d__1 = bli_zreal(ap[i__3]);
750 bli_zsets( (d__1), (0.), ap[i__2] );
751 }
752 kk = kk + *n - j + 1;
753 /* L60: */
754 }
755 } else {
756 i__1 = *n;
757 for (j = 1; j <= i__1; ++j) {
758 i__2 = jx;
759 i__3 = jy;
760 if (bli_zreal(x[i__2]) != 0. || bli_zimag(x[i__2]) != 0. || (bli_zreal(y[i__3]) != 0. ||
761 bli_zimag(y[i__3]) != 0.)) {
762 bla_d_cnjg(&z__2, &y[jy]);
763 bli_zsets( (bli_zreal(*alpha) * bli_zreal(z__2) - bli_zimag(*alpha) * bli_zimag(z__2)), (bli_zreal(*alpha) * bli_zimag(z__2) + bli_zimag(*alpha) * bli_zreal(z__2)), z__1 );
764 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp1 );
765 i__2 = jx;
766 bli_zsets( (bli_zreal(*alpha) * bli_zreal(x[i__2]) - bli_zimag(*alpha) * bli_zimag(x[i__2])), (bli_zreal(*alpha) * bli_zimag(x[i__2]) + bli_zimag(*alpha) * bli_zreal(x[i__2])), z__2 );
767 bla_d_cnjg(&z__1, &z__2);
768 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), temp2 );
769 i__2 = kk;
770 i__3 = kk;
771 i__4 = jx;
772 bli_zsets( (bli_zreal(x[i__4]) * bli_zreal(temp1) - bli_zimag(x[i__4]) * bli_zimag(temp1)), (bli_zreal(x[i__4]) * bli_zimag(temp1) + bli_zimag(x[i__4]) * bli_zreal(temp1)), z__2 );
773 i__5 = jy;
774 bli_zsets( (bli_zreal(y[i__5]) * bli_zreal(temp2) - bli_zimag(y[i__5]) * bli_zimag(temp2)), (bli_zreal(y[i__5]) * bli_zimag(temp2) + bli_zimag(y[i__5]) * bli_zreal(temp2)), z__3 );
775 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__3)), (bli_zimag(z__2) + bli_zimag(z__3)), z__1 );
776 d__1 = bli_zreal(ap[i__3]) + bli_zreal(z__1);
777 bli_zsets( (d__1), (0.), ap[i__2] );
778 ix = jx;
779 iy = jy;
780 i__2 = kk + *n - j;
781 for (k = kk + 1; k <= i__2; ++k) {
782 ix += *incx;
783 iy += *incy;
784 i__3 = k;
785 i__4 = k;
786 i__5 = ix;
787 bli_zsets( (bli_zreal(x[i__5]) * bli_zreal(temp1) - bli_zimag(x[i__5]) * bli_zimag(temp1)), (bli_zreal(x[i__5]) * bli_zimag(temp1) + bli_zimag(x[i__5]) * bli_zreal(temp1)), z__3 );
788 bli_zsets( (bli_zreal(ap[i__4]) + bli_zreal(z__3)), (bli_zimag(ap[i__4]) + bli_zimag(z__3)), z__2 );
789 i__6 = iy;
790 bli_zsets( (bli_zreal(y[i__6]) * bli_zreal(temp2) - bli_zimag(y[i__6]) * bli_zimag(temp2)), (bli_zreal(y[i__6]) * bli_zimag(temp2) + bli_zimag(y[i__6]) * bli_zreal(temp2)), z__4 );
791 bli_zsets( (bli_zreal(z__2) + bli_zreal(z__4)), (bli_zimag(z__2) + bli_zimag(z__4)), z__1 );
792 bli_zsets( (bli_zreal(z__1)), (bli_zimag(z__1)), ap[i__3] );
793 /* L70: */
794 }
795 } else {
796 i__2 = kk;
797 i__3 = kk;
798 d__1 = bli_zreal(ap[i__3]);
799 bli_zsets( (d__1), (0.), ap[i__2] );
800 }
801 jx += *incx;
802 jy += *incy;
803 kk = kk + *n - j + 1;
804 /* L80: */
805 }
806 }
807 }
809 return 0;
811 /* End of ZHPR2 . */
813 } /* zhpr2_ */
815 #endif