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

gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:14