chpmv.c
Go to the documentation of this file.
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 
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 
complex::i
real i
Definition: datatypes.h:12
y
Scalar * y
Definition: level1_cplx_impl.h:124
r_cnjg
void r_cnjg(complex *r, complex *z)
Definition: r_cnjg.c:3
chpmv_
int chpmv_(char *uplo, integer *n, complex *alpha, complex *ap, complex *x, integer *incx, complex *beta, complex *y, integer *incy, ftnlen uplo_len)
Definition: chpmv.c:15
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
n
int n
Definition: BiCGSTAB_simple.cpp:1
ftnlen
int ftnlen
Definition: datatypes.h:14
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
real
Definition: main.h:100
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
logical
int logical
Definition: datatypes.h:15
lsame_
logical lsame_(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: lsame.c:15
integer
int integer
Definition: datatypes.h:8
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
incy
int RealScalar int RealScalar int * incy
Definition: level1_cplx_impl.h:119
complex::r
real r
Definition: datatypes.h:12

gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:01