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

gtsam
Author(s):
autogenerated on Sat May 8 2021 02:41:46