• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /* chpmv.f -- translated by f2c (version 20100827).
2    You must link the resulting object file with libf2c:
3 	on Microsoft Windows system, link with libf2c.lib;
4 	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 	or, if you install libf2c.a in a standard place, with -lf2c -lm
6 	-- in that order, at the end of the command line, as in
7 		cc *.o -lf2c -lm
8 	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9 
10 		http://www.netlib.org/f2c/libf2c.zip
11 */
12 
13 #include "datatypes.h"
14 
chpmv_(char * uplo,integer * n,complex * alpha,complex * ap,complex * x,integer * incx,complex * beta,complex * y,integer * incy,ftnlen uplo_len)15 /* Subroutine */ int chpmv_(char *uplo, integer *n, complex *alpha, complex *
16 	ap, complex *x, integer *incx, complex *beta, complex *y, integer *
17 	incy, ftnlen uplo_len)
18 {
19     /* System generated locals */
20     integer i__1, i__2, i__3, i__4, i__5;
21     real r__1;
22     complex q__1, q__2, q__3, q__4;
23 
24     /* Builtin functions */
25     void r_cnjg(complex *, complex *);
26 
27     /* Local variables */
28     integer i__, j, k, kk, ix, iy, jx, jy, kx, ky, info;
29     complex temp1, temp2;
30     extern logical lsame_(char *, char *, ftnlen, ftnlen);
31     extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
32 
33 /*     .. Scalar Arguments .. */
34 /*     .. */
35 /*     .. Array Arguments .. */
36 /*     .. */
37 
38 /*  Purpose */
39 /*  ======= */
40 
41 /*  CHPMV  performs the matrix-vector operation */
42 
43 /*     y := alpha*A*x + beta*y, */
44 
45 /*  where alpha and beta are scalars, x and y are n element vectors and */
46 /*  A is an n by n hermitian matrix, supplied in packed form. */
47 
48 /*  Arguments */
49 /*  ========== */
50 
51 /*  UPLO   - CHARACTER*1. */
52 /*           On entry, UPLO specifies whether the upper or lower */
53 /*           triangular part of the matrix A is supplied in the packed */
54 /*           array AP as follows: */
55 
56 /*              UPLO = 'U' or 'u'   The upper triangular part of A is */
57 /*                                  supplied in AP. */
58 
59 /*              UPLO = 'L' or 'l'   The lower triangular part of A is */
60 /*                                  supplied in AP. */
61 
62 /*           Unchanged on exit. */
63 
64 /*  N      - INTEGER. */
65 /*           On entry, N specifies the order of the matrix A. */
66 /*           N must be at least zero. */
67 /*           Unchanged on exit. */
68 
69 /*  ALPHA  - COMPLEX         . */
70 /*           On entry, ALPHA specifies the scalar alpha. */
71 /*           Unchanged on exit. */
72 
73 /*  AP     - COMPLEX          array of DIMENSION at least */
74 /*           ( ( n*( n + 1 ) )/2 ). */
75 /*           Before entry with UPLO = 'U' or 'u', the array AP must */
76 /*           contain the upper triangular part of the hermitian matrix */
77 /*           packed sequentially, column by column, so that AP( 1 ) */
78 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) */
79 /*           and a( 2, 2 ) respectively, and so on. */
80 /*           Before entry with UPLO = 'L' or 'l', the array AP must */
81 /*           contain the lower triangular part of the hermitian matrix */
82 /*           packed sequentially, column by column, so that AP( 1 ) */
83 /*           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) */
84 /*           and a( 3, 1 ) respectively, and so on. */
85 /*           Note that the imaginary parts of the diagonal elements need */
86 /*           not be set and are assumed to be zero. */
87 /*           Unchanged on exit. */
88 
89 /*  X      - COMPLEX          array of dimension at least */
90 /*           ( 1 + ( n - 1 )*abs( INCX ) ). */
91 /*           Before entry, the incremented array X must contain the n */
92 /*           element vector x. */
93 /*           Unchanged on exit. */
94 
95 /*  INCX   - INTEGER. */
96 /*           On entry, INCX specifies the increment for the elements of */
97 /*           X. INCX must not be zero. */
98 /*           Unchanged on exit. */
99 
100 /*  BETA   - COMPLEX         . */
101 /*           On entry, BETA specifies the scalar beta. When BETA is */
102 /*           supplied as zero then Y need not be set on input. */
103 /*           Unchanged on exit. */
104 
105 /*  Y      - COMPLEX          array of dimension at least */
106 /*           ( 1 + ( n - 1 )*abs( INCY ) ). */
107 /*           Before entry, the incremented array Y must contain the n */
108 /*           element vector y. On exit, Y is overwritten by the updated */
109 /*           vector y. */
110 
111 /*  INCY   - INTEGER. */
112 /*           On entry, INCY specifies the increment for the elements of */
113 /*           Y. INCY must not be zero. */
114 /*           Unchanged on exit. */
115 
116 /*  Further Details */
117 /*  =============== */
118 
119 /*  Level 2 Blas routine. */
120 
121 /*  -- Written on 22-October-1986. */
122 /*     Jack Dongarra, Argonne National Lab. */
123 /*     Jeremy Du Croz, Nag Central Office. */
124 /*     Sven Hammarling, Nag Central Office. */
125 /*     Richard Hanson, Sandia National Labs. */
126 
127 /*  ===================================================================== */
128 
129 /*     .. Parameters .. */
130 /*     .. */
131 /*     .. Local Scalars .. */
132 /*     .. */
133 /*     .. External Functions .. */
134 /*     .. */
135 /*     .. External Subroutines .. */
136 /*     .. */
137 /*     .. Intrinsic Functions .. */
138 /*     .. */
139 
140 /*     Test the input parameters. */
141 
142     /* Parameter adjustments */
143     --y;
144     --x;
145     --ap;
146 
147     /* Function Body */
148     info = 0;
149     if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
150 	    ftnlen)1, (ftnlen)1)) {
151 	info = 1;
152     } else if (*n < 0) {
153 	info = 2;
154     } else if (*incx == 0) {
155 	info = 6;
156     } else if (*incy == 0) {
157 	info = 9;
158     }
159     if (info != 0) {
160 	xerbla_("CHPMV ", &info, (ftnlen)6);
161 	return 0;
162     }
163 
164 /*     Quick return if possible. */
165 
166     if (*n == 0 || (alpha->r == 0.f && alpha->i == 0.f && (beta->r == 1.f &&
167                                                            beta->i == 0.f))) {
168 	return 0;
169     }
170 
171 /*     Set up the start points in  X  and  Y. */
172 
173     if (*incx > 0) {
174 	kx = 1;
175     } else {
176 	kx = 1 - (*n - 1) * *incx;
177     }
178     if (*incy > 0) {
179 	ky = 1;
180     } else {
181 	ky = 1 - (*n - 1) * *incy;
182     }
183 
184 /*     Start the operations. In this version the elements of the array AP */
185 /*     are accessed sequentially with one pass through AP. */
186 
187 /*     First form  y := beta*y. */
188 
189     if (beta->r != 1.f || beta->i != 0.f) {
190 	if (*incy == 1) {
191 	    if (beta->r == 0.f && beta->i == 0.f) {
192 		i__1 = *n;
193 		for (i__ = 1; i__ <= i__1; ++i__) {
194 		    i__2 = i__;
195 		    y[i__2].r = 0.f, y[i__2].i = 0.f;
196 /* L10: */
197 		}
198 	    } else {
199 		i__1 = *n;
200 		for (i__ = 1; i__ <= i__1; ++i__) {
201 		    i__2 = i__;
202 		    i__3 = i__;
203 		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
204 			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
205 			    .r;
206 		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
207 /* L20: */
208 		}
209 	    }
210 	} else {
211 	    iy = ky;
212 	    if (beta->r == 0.f && beta->i == 0.f) {
213 		i__1 = *n;
214 		for (i__ = 1; i__ <= i__1; ++i__) {
215 		    i__2 = iy;
216 		    y[i__2].r = 0.f, y[i__2].i = 0.f;
217 		    iy += *incy;
218 /* L30: */
219 		}
220 	    } else {
221 		i__1 = *n;
222 		for (i__ = 1; i__ <= i__1; ++i__) {
223 		    i__2 = iy;
224 		    i__3 = iy;
225 		    q__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
226 			    q__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
227 			    .r;
228 		    y[i__2].r = q__1.r, y[i__2].i = q__1.i;
229 		    iy += *incy;
230 /* L40: */
231 		}
232 	    }
233 	}
234     }
235     if (alpha->r == 0.f && alpha->i == 0.f) {
236 	return 0;
237     }
238     kk = 1;
239     if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
240 
241 /*        Form  y  when AP contains the upper triangle. */
242 
243 	if (*incx == 1 && *incy == 1) {
244 	    i__1 = *n;
245 	    for (j = 1; j <= i__1; ++j) {
246 		i__2 = j;
247 		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
248 			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
249 		temp1.r = q__1.r, temp1.i = q__1.i;
250 		temp2.r = 0.f, temp2.i = 0.f;
251 		k = kk;
252 		i__2 = j - 1;
253 		for (i__ = 1; i__ <= i__2; ++i__) {
254 		    i__3 = i__;
255 		    i__4 = i__;
256 		    i__5 = k;
257 		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i,
258 			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
259 			    .r;
260 		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
261 		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
262 		    r_cnjg(&q__3, &ap[k]);
263 		    i__3 = i__;
264 		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
265 			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
266 		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
267 		    temp2.r = q__1.r, temp2.i = q__1.i;
268 		    ++k;
269 /* L50: */
270 		}
271 		i__2 = j;
272 		i__3 = j;
273 		i__4 = kk + j - 1;
274 		r__1 = ap[i__4].r;
275 		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
276 		q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i;
277 		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i =
278 			alpha->r * temp2.i + alpha->i * temp2.r;
279 		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
280 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
281 		kk += j;
282 /* L60: */
283 	    }
284 	} else {
285 	    jx = kx;
286 	    jy = ky;
287 	    i__1 = *n;
288 	    for (j = 1; j <= i__1; ++j) {
289 		i__2 = jx;
290 		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
291 			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
292 		temp1.r = q__1.r, temp1.i = q__1.i;
293 		temp2.r = 0.f, temp2.i = 0.f;
294 		ix = kx;
295 		iy = ky;
296 		i__2 = kk + j - 2;
297 		for (k = kk; k <= i__2; ++k) {
298 		    i__3 = iy;
299 		    i__4 = iy;
300 		    i__5 = k;
301 		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i,
302 			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
303 			    .r;
304 		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
305 		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
306 		    r_cnjg(&q__3, &ap[k]);
307 		    i__3 = ix;
308 		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
309 			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
310 		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
311 		    temp2.r = q__1.r, temp2.i = q__1.i;
312 		    ix += *incx;
313 		    iy += *incy;
314 /* L70: */
315 		}
316 		i__2 = jy;
317 		i__3 = jy;
318 		i__4 = kk + j - 1;
319 		r__1 = ap[i__4].r;
320 		q__3.r = r__1 * temp1.r, q__3.i = r__1 * temp1.i;
321 		q__2.r = y[i__3].r + q__3.r, q__2.i = y[i__3].i + q__3.i;
322 		q__4.r = alpha->r * temp2.r - alpha->i * temp2.i, q__4.i =
323 			alpha->r * temp2.i + alpha->i * temp2.r;
324 		q__1.r = q__2.r + q__4.r, q__1.i = q__2.i + q__4.i;
325 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
326 		jx += *incx;
327 		jy += *incy;
328 		kk += j;
329 /* L80: */
330 	    }
331 	}
332     } else {
333 
334 /*        Form  y  when AP contains the lower triangle. */
335 
336 	if (*incx == 1 && *incy == 1) {
337 	    i__1 = *n;
338 	    for (j = 1; j <= i__1; ++j) {
339 		i__2 = j;
340 		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
341 			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
342 		temp1.r = q__1.r, temp1.i = q__1.i;
343 		temp2.r = 0.f, temp2.i = 0.f;
344 		i__2 = j;
345 		i__3 = j;
346 		i__4 = kk;
347 		r__1 = ap[i__4].r;
348 		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
349 		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
350 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
351 		k = kk + 1;
352 		i__2 = *n;
353 		for (i__ = j + 1; i__ <= i__2; ++i__) {
354 		    i__3 = i__;
355 		    i__4 = i__;
356 		    i__5 = k;
357 		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i,
358 			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
359 			    .r;
360 		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
361 		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
362 		    r_cnjg(&q__3, &ap[k]);
363 		    i__3 = i__;
364 		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
365 			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
366 		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
367 		    temp2.r = q__1.r, temp2.i = q__1.i;
368 		    ++k;
369 /* L90: */
370 		}
371 		i__2 = j;
372 		i__3 = j;
373 		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i =
374 			alpha->r * temp2.i + alpha->i * temp2.r;
375 		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
376 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
377 		kk += *n - j + 1;
378 /* L100: */
379 	    }
380 	} else {
381 	    jx = kx;
382 	    jy = ky;
383 	    i__1 = *n;
384 	    for (j = 1; j <= i__1; ++j) {
385 		i__2 = jx;
386 		q__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, q__1.i =
387 			 alpha->r * x[i__2].i + alpha->i * x[i__2].r;
388 		temp1.r = q__1.r, temp1.i = q__1.i;
389 		temp2.r = 0.f, temp2.i = 0.f;
390 		i__2 = jy;
391 		i__3 = jy;
392 		i__4 = kk;
393 		r__1 = ap[i__4].r;
394 		q__2.r = r__1 * temp1.r, q__2.i = r__1 * temp1.i;
395 		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
396 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
397 		ix = jx;
398 		iy = jy;
399 		i__2 = kk + *n - j;
400 		for (k = kk + 1; k <= i__2; ++k) {
401 		    ix += *incx;
402 		    iy += *incy;
403 		    i__3 = iy;
404 		    i__4 = iy;
405 		    i__5 = k;
406 		    q__2.r = temp1.r * ap[i__5].r - temp1.i * ap[i__5].i,
407 			    q__2.i = temp1.r * ap[i__5].i + temp1.i * ap[i__5]
408 			    .r;
409 		    q__1.r = y[i__4].r + q__2.r, q__1.i = y[i__4].i + q__2.i;
410 		    y[i__3].r = q__1.r, y[i__3].i = q__1.i;
411 		    r_cnjg(&q__3, &ap[k]);
412 		    i__3 = ix;
413 		    q__2.r = q__3.r * x[i__3].r - q__3.i * x[i__3].i, q__2.i =
414 			     q__3.r * x[i__3].i + q__3.i * x[i__3].r;
415 		    q__1.r = temp2.r + q__2.r, q__1.i = temp2.i + q__2.i;
416 		    temp2.r = q__1.r, temp2.i = q__1.i;
417 /* L110: */
418 		}
419 		i__2 = jy;
420 		i__3 = jy;
421 		q__2.r = alpha->r * temp2.r - alpha->i * temp2.i, q__2.i =
422 			alpha->r * temp2.i + alpha->i * temp2.r;
423 		q__1.r = y[i__3].r + q__2.r, q__1.i = y[i__3].i + q__2.i;
424 		y[i__2].r = q__1.r, y[i__2].i = q__1.i;
425 		jx += *incx;
426 		jy += *incy;
427 		kk += *n - j + 1;
428 /* L120: */
429 	    }
430 	}
431     }
432 
433     return 0;
434 
435 /*     End of CHPMV . */
436 
437 } /* chpmv_ */
438 
439