zhbmv.c
Go to the documentation of this file.
1 /* zhbmv.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 zhbmv_(char *uplo, integer *n, integer *k, doublecomplex
16  *alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *
17  incx, doublecomplex *beta, doublecomplex *y, integer *incy, ftnlen
18  uplo_len)
19 {
20  /* System generated locals */
21  integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
22  doublereal d__1;
23  doublecomplex z__1, z__2, z__3, z__4;
24 
25  /* Builtin functions */
26  void d_cnjg(doublecomplex *, doublecomplex *);
27 
28  /* Local variables */
29  integer i__, j, l, ix, iy, jx, jy, kx, ky, info;
30  doublecomplex temp1, temp2;
31  extern logical lsame_(char *, char *, ftnlen, ftnlen);
32  integer kplus1;
33  extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
34 
35 /* .. Scalar Arguments .. */
36 /* .. */
37 /* .. Array Arguments .. */
38 /* .. */
39 
40 /* Purpose */
41 /* ======= */
42 
43 /* ZHBMV performs the matrix-vector operation */
44 
45 /* y := alpha*A*x + beta*y, */
46 
47 /* where alpha and beta are scalars, x and y are n element vectors and */
48 /* A is an n by n hermitian band matrix, with k super-diagonals. */
49 
50 /* Arguments */
51 /* ========== */
52 
53 /* UPLO - CHARACTER*1. */
54 /* On entry, UPLO specifies whether the upper or lower */
55 /* triangular part of the band matrix A is being supplied as */
56 /* follows: */
57 
58 /* UPLO = 'U' or 'u' The upper triangular part of A is */
59 /* being supplied. */
60 
61 /* UPLO = 'L' or 'l' The lower triangular part of A is */
62 /* being supplied. */
63 
64 /* Unchanged on exit. */
65 
66 /* N - INTEGER. */
67 /* On entry, N specifies the order of the matrix A. */
68 /* N must be at least zero. */
69 /* Unchanged on exit. */
70 
71 /* K - INTEGER. */
72 /* On entry, K specifies the number of super-diagonals of the */
73 /* matrix A. K must satisfy 0 .le. K. */
74 /* Unchanged on exit. */
75 
76 /* ALPHA - COMPLEX*16 . */
77 /* On entry, ALPHA specifies the scalar alpha. */
78 /* Unchanged on exit. */
79 
80 /* A - COMPLEX*16 array of DIMENSION ( LDA, n ). */
81 /* Before entry with UPLO = 'U' or 'u', the leading ( k + 1 ) */
82 /* by n part of the array A must contain the upper triangular */
83 /* band part of the hermitian matrix, supplied column by */
84 /* column, with the leading diagonal of the matrix in row */
85 /* ( k + 1 ) of the array, the first super-diagonal starting at */
86 /* position 2 in row k, and so on. The top left k by k triangle */
87 /* of the array A is not referenced. */
88 /* The following program segment will transfer the upper */
89 /* triangular part of a hermitian band matrix from conventional */
90 /* full matrix storage to band storage: */
91 
92 /* DO 20, J = 1, N */
93 /* M = K + 1 - J */
94 /* DO 10, I = MAX( 1, J - K ), J */
95 /* A( M + I, J ) = matrix( I, J ) */
96 /* 10 CONTINUE */
97 /* 20 CONTINUE */
98 
99 /* Before entry with UPLO = 'L' or 'l', the leading ( k + 1 ) */
100 /* by n part of the array A must contain the lower triangular */
101 /* band part of the hermitian matrix, supplied column by */
102 /* column, with the leading diagonal of the matrix in row 1 of */
103 /* the array, the first sub-diagonal starting at position 1 in */
104 /* row 2, and so on. The bottom right k by k triangle of the */
105 /* array A is not referenced. */
106 /* The following program segment will transfer the lower */
107 /* triangular part of a hermitian band matrix from conventional */
108 /* full matrix storage to band storage: */
109 
110 /* DO 20, J = 1, N */
111 /* M = 1 - J */
112 /* DO 10, I = J, MIN( N, J + K ) */
113 /* A( M + I, J ) = matrix( I, J ) */
114 /* 10 CONTINUE */
115 /* 20 CONTINUE */
116 
117 /* Note that the imaginary parts of the diagonal elements need */
118 /* not be set and are assumed to be zero. */
119 /* Unchanged on exit. */
120 
121 /* LDA - INTEGER. */
122 /* On entry, LDA specifies the first dimension of A as declared */
123 /* in the calling (sub) program. LDA must be at least */
124 /* ( k + 1 ). */
125 /* Unchanged on exit. */
126 
127 /* X - COMPLEX*16 array of DIMENSION at least */
128 /* ( 1 + ( n - 1 )*abs( INCX ) ). */
129 /* Before entry, the incremented array X must contain the */
130 /* vector x. */
131 /* Unchanged on exit. */
132 
133 /* INCX - INTEGER. */
134 /* On entry, INCX specifies the increment for the elements of */
135 /* X. INCX must not be zero. */
136 /* Unchanged on exit. */
137 
138 /* BETA - COMPLEX*16 . */
139 /* On entry, BETA specifies the scalar beta. */
140 /* Unchanged on exit. */
141 
142 /* Y - COMPLEX*16 array of DIMENSION at least */
143 /* ( 1 + ( n - 1 )*abs( INCY ) ). */
144 /* Before entry, the incremented array Y must contain the */
145 /* vector y. On exit, Y is overwritten by the updated vector y. */
146 
147 /* INCY - INTEGER. */
148 /* On entry, INCY specifies the increment for the elements of */
149 /* Y. INCY must not be zero. */
150 /* Unchanged on exit. */
151 
152 /* Further Details */
153 /* =============== */
154 
155 /* Level 2 Blas routine. */
156 
157 /* -- Written on 22-October-1986. */
158 /* Jack Dongarra, Argonne National Lab. */
159 /* Jeremy Du Croz, Nag Central Office. */
160 /* Sven Hammarling, Nag Central Office. */
161 /* Richard Hanson, Sandia National Labs. */
162 
163 /* ===================================================================== */
164 
165 /* .. Parameters .. */
166 /* .. */
167 /* .. Local Scalars .. */
168 /* .. */
169 /* .. External Functions .. */
170 /* .. */
171 /* .. External Subroutines .. */
172 /* .. */
173 /* .. Intrinsic Functions .. */
174 /* .. */
175 
176 /* Test the input parameters. */
177 
178  /* Parameter adjustments */
179  a_dim1 = *lda;
180  a_offset = 1 + a_dim1;
181  a -= a_offset;
182  --x;
183  --y;
184 
185  /* Function Body */
186  info = 0;
187  if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
188  ftnlen)1, (ftnlen)1)) {
189  info = 1;
190  } else if (*n < 0) {
191  info = 2;
192  } else if (*k < 0) {
193  info = 3;
194  } else if (*lda < *k + 1) {
195  info = 6;
196  } else if (*incx == 0) {
197  info = 8;
198  } else if (*incy == 0) {
199  info = 11;
200  }
201  if (info != 0) {
202  xerbla_("ZHBMV ", &info, (ftnlen)6);
203  return 0;
204  }
205 
206 /* Quick return if possible. */
207 
208  if (*n == 0 || (alpha->r == 0. && alpha->i == 0. && (beta->r == 1. &&
209  beta->i == 0.))) {
210  return 0;
211  }
212 
213 /* Set up the start points in X and Y. */
214 
215  if (*incx > 0) {
216  kx = 1;
217  } else {
218  kx = 1 - (*n - 1) * *incx;
219  }
220  if (*incy > 0) {
221  ky = 1;
222  } else {
223  ky = 1 - (*n - 1) * *incy;
224  }
225 
226 /* Start the operations. In this version the elements of the array A */
227 /* are accessed sequentially with one pass through A. */
228 
229 /* First form y := beta*y. */
230 
231  if (beta->r != 1. || beta->i != 0.) {
232  if (*incy == 1) {
233  if (beta->r == 0. && beta->i == 0.) {
234  i__1 = *n;
235  for (i__ = 1; i__ <= i__1; ++i__) {
236  i__2 = i__;
237  y[i__2].r = 0., y[i__2].i = 0.;
238 /* L10: */
239  }
240  } else {
241  i__1 = *n;
242  for (i__ = 1; i__ <= i__1; ++i__) {
243  i__2 = i__;
244  i__3 = i__;
245  z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
246  z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
247  .r;
248  y[i__2].r = z__1.r, y[i__2].i = z__1.i;
249 /* L20: */
250  }
251  }
252  } else {
253  iy = ky;
254  if (beta->r == 0. && beta->i == 0.) {
255  i__1 = *n;
256  for (i__ = 1; i__ <= i__1; ++i__) {
257  i__2 = iy;
258  y[i__2].r = 0., y[i__2].i = 0.;
259  iy += *incy;
260 /* L30: */
261  }
262  } else {
263  i__1 = *n;
264  for (i__ = 1; i__ <= i__1; ++i__) {
265  i__2 = iy;
266  i__3 = iy;
267  z__1.r = beta->r * y[i__3].r - beta->i * y[i__3].i,
268  z__1.i = beta->r * y[i__3].i + beta->i * y[i__3]
269  .r;
270  y[i__2].r = z__1.r, y[i__2].i = z__1.i;
271  iy += *incy;
272 /* L40: */
273  }
274  }
275  }
276  }
277  if (alpha->r == 0. && alpha->i == 0.) {
278  return 0;
279  }
280  if (lsame_(uplo, "U", (ftnlen)1, (ftnlen)1)) {
281 
282 /* Form y when upper triangle of A is stored. */
283 
284  kplus1 = *k + 1;
285  if (*incx == 1 && *incy == 1) {
286  i__1 = *n;
287  for (j = 1; j <= i__1; ++j) {
288  i__2 = j;
289  z__1.r = alpha->r * x[i__2].r - alpha->i * x[i__2].i, z__1.i =
290  alpha->r * x[i__2].i + alpha->i * x[i__2].r;
291  temp1.r = z__1.r, temp1.i = z__1.i;
292  temp2.r = 0., temp2.i = 0.;
293  l = kplus1 - j;
294 /* Computing MAX */
295  i__2 = 1, i__3 = j - *k;
296  i__4 = j - 1;
297  for (i__ = max(i__2,i__3); i__ <= i__4; ++i__) {
298  i__2 = i__;
299  i__3 = i__;
300  i__5 = l + i__ + j * a_dim1;
301  z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
302  z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
303  .r;
304  z__1.r = y[i__3].r + z__2.r, z__1.i = y[i__3].i + z__2.i;
305  y[i__2].r = z__1.r, y[i__2].i = z__1.i;
306  d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
307  i__2 = i__;
308  z__2.r = z__3.r * x[i__2].r - z__3.i * x[i__2].i, z__2.i =
309  z__3.r * x[i__2].i + z__3.i * x[i__2].r;
310  z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
311  temp2.r = z__1.r, temp2.i = z__1.i;
312 /* L50: */
313  }
314  i__4 = j;
315  i__2 = j;
316  i__3 = kplus1 + j * a_dim1;
317  d__1 = a[i__3].r;
318  z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
319  z__2.r = y[i__2].r + z__3.r, z__2.i = y[i__2].i + z__3.i;
320  z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i =
321  alpha->r * temp2.i + alpha->i * temp2.r;
322  z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
323  y[i__4].r = z__1.r, y[i__4].i = z__1.i;
324 /* L60: */
325  }
326  } else {
327  jx = kx;
328  jy = ky;
329  i__1 = *n;
330  for (j = 1; j <= i__1; ++j) {
331  i__4 = jx;
332  z__1.r = alpha->r * x[i__4].r - alpha->i * x[i__4].i, z__1.i =
333  alpha->r * x[i__4].i + alpha->i * x[i__4].r;
334  temp1.r = z__1.r, temp1.i = z__1.i;
335  temp2.r = 0., temp2.i = 0.;
336  ix = kx;
337  iy = ky;
338  l = kplus1 - j;
339 /* Computing MAX */
340  i__4 = 1, i__2 = j - *k;
341  i__3 = j - 1;
342  for (i__ = max(i__4,i__2); i__ <= i__3; ++i__) {
343  i__4 = iy;
344  i__2 = iy;
345  i__5 = l + i__ + j * a_dim1;
346  z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
347  z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
348  .r;
349  z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
350  y[i__4].r = z__1.r, y[i__4].i = z__1.i;
351  d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
352  i__4 = ix;
353  z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
354  z__3.r * x[i__4].i + z__3.i * x[i__4].r;
355  z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
356  temp2.r = z__1.r, temp2.i = z__1.i;
357  ix += *incx;
358  iy += *incy;
359 /* L70: */
360  }
361  i__3 = jy;
362  i__4 = jy;
363  i__2 = kplus1 + j * a_dim1;
364  d__1 = a[i__2].r;
365  z__3.r = d__1 * temp1.r, z__3.i = d__1 * temp1.i;
366  z__2.r = y[i__4].r + z__3.r, z__2.i = y[i__4].i + z__3.i;
367  z__4.r = alpha->r * temp2.r - alpha->i * temp2.i, z__4.i =
368  alpha->r * temp2.i + alpha->i * temp2.r;
369  z__1.r = z__2.r + z__4.r, z__1.i = z__2.i + z__4.i;
370  y[i__3].r = z__1.r, y[i__3].i = z__1.i;
371  jx += *incx;
372  jy += *incy;
373  if (j > *k) {
374  kx += *incx;
375  ky += *incy;
376  }
377 /* L80: */
378  }
379  }
380  } else {
381 
382 /* Form y when lower triangle of A is stored. */
383 
384  if (*incx == 1 && *incy == 1) {
385  i__1 = *n;
386  for (j = 1; j <= i__1; ++j) {
387  i__3 = j;
388  z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
389  alpha->r * x[i__3].i + alpha->i * x[i__3].r;
390  temp1.r = z__1.r, temp1.i = z__1.i;
391  temp2.r = 0., temp2.i = 0.;
392  i__3 = j;
393  i__4 = j;
394  i__2 = j * a_dim1 + 1;
395  d__1 = a[i__2].r;
396  z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
397  z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
398  y[i__3].r = z__1.r, y[i__3].i = z__1.i;
399  l = 1 - j;
400 /* Computing MIN */
401  i__4 = *n, i__2 = j + *k;
402  i__3 = min(i__4,i__2);
403  for (i__ = j + 1; i__ <= i__3; ++i__) {
404  i__4 = i__;
405  i__2 = i__;
406  i__5 = l + i__ + j * a_dim1;
407  z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
408  z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
409  .r;
410  z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
411  y[i__4].r = z__1.r, y[i__4].i = z__1.i;
412  d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
413  i__4 = i__;
414  z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
415  z__3.r * x[i__4].i + z__3.i * x[i__4].r;
416  z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
417  temp2.r = z__1.r, temp2.i = z__1.i;
418 /* L90: */
419  }
420  i__3 = j;
421  i__4 = j;
422  z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i =
423  alpha->r * temp2.i + alpha->i * temp2.r;
424  z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
425  y[i__3].r = z__1.r, y[i__3].i = z__1.i;
426 /* L100: */
427  }
428  } else {
429  jx = kx;
430  jy = ky;
431  i__1 = *n;
432  for (j = 1; j <= i__1; ++j) {
433  i__3 = jx;
434  z__1.r = alpha->r * x[i__3].r - alpha->i * x[i__3].i, z__1.i =
435  alpha->r * x[i__3].i + alpha->i * x[i__3].r;
436  temp1.r = z__1.r, temp1.i = z__1.i;
437  temp2.r = 0., temp2.i = 0.;
438  i__3 = jy;
439  i__4 = jy;
440  i__2 = j * a_dim1 + 1;
441  d__1 = a[i__2].r;
442  z__2.r = d__1 * temp1.r, z__2.i = d__1 * temp1.i;
443  z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
444  y[i__3].r = z__1.r, y[i__3].i = z__1.i;
445  l = 1 - j;
446  ix = jx;
447  iy = jy;
448 /* Computing MIN */
449  i__4 = *n, i__2 = j + *k;
450  i__3 = min(i__4,i__2);
451  for (i__ = j + 1; i__ <= i__3; ++i__) {
452  ix += *incx;
453  iy += *incy;
454  i__4 = iy;
455  i__2 = iy;
456  i__5 = l + i__ + j * a_dim1;
457  z__2.r = temp1.r * a[i__5].r - temp1.i * a[i__5].i,
458  z__2.i = temp1.r * a[i__5].i + temp1.i * a[i__5]
459  .r;
460  z__1.r = y[i__2].r + z__2.r, z__1.i = y[i__2].i + z__2.i;
461  y[i__4].r = z__1.r, y[i__4].i = z__1.i;
462  d_cnjg(&z__3, &a[l + i__ + j * a_dim1]);
463  i__4 = ix;
464  z__2.r = z__3.r * x[i__4].r - z__3.i * x[i__4].i, z__2.i =
465  z__3.r * x[i__4].i + z__3.i * x[i__4].r;
466  z__1.r = temp2.r + z__2.r, z__1.i = temp2.i + z__2.i;
467  temp2.r = z__1.r, temp2.i = z__1.i;
468 /* L110: */
469  }
470  i__3 = jy;
471  i__4 = jy;
472  z__2.r = alpha->r * temp2.r - alpha->i * temp2.i, z__2.i =
473  alpha->r * temp2.i + alpha->i * temp2.r;
474  z__1.r = y[i__4].r + z__2.r, z__1.i = y[i__4].i + z__2.i;
475  y[i__3].r = z__1.r, y[i__3].i = z__1.i;
476  jx += *incx;
477  jy += *incy;
478 /* L120: */
479  }
480  }
481  }
482 
483  return 0;
484 
485 /* End of ZHBMV . */
486 
487 } /* zhbmv_ */
488 
zhbmv_
int zhbmv_(char *uplo, integer *n, integer *k, doublecomplex *alpha, doublecomplex *a, integer *lda, doublecomplex *x, integer *incx, doublecomplex *beta, doublecomplex *y, integer *incy, ftnlen uplo_len)
Definition: zhbmv.c:15
max
#define max(a, b)
Definition: datatypes.h:20
y
Scalar * y
Definition: level1_cplx_impl.h:102
min
#define min(a, b)
Definition: datatypes.h:19
incx
RealScalar RealScalar int * incx
Definition: level1_cplx_impl.h:29
d_cnjg
void d_cnjg(doublecomplex *r, doublecomplex *z)
Definition: d_cnjg.c:3
n
int n
Definition: BiCGSTAB_simple.cpp:1
ftnlen
int ftnlen
Definition: datatypes.h:14
doublecomplex::r
doublereal r
Definition: datatypes.h:13
a
Array33i a
Definition: Cwise_product.cpp:1
info
else if n * info
Definition: 3rdparty/Eigen/lapack/cholesky.cpp:18
l
static const Line3 l(Rot3(), 1, 1)
xerbla_
EIGEN_WEAK_LINKING int xerbla_(const char *msg, int *info, int)
Definition: xerbla.cpp:15
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:125
lda
* lda
Definition: eigenvalues.cpp:59
doublecomplex::i
doublereal i
Definition: datatypes.h:13
logical
int logical
Definition: datatypes.h:15
lsame_
logical lsame_(char *ca, char *cb, ftnlen ca_len, ftnlen cb_len)
Definition: lsame.c:15
doublereal
double doublereal
Definition: datatypes.h:11
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:97

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:51:26