r1updt.h
Go to the documentation of this file.
1 namespace Eigen {
2 
3 namespace internal {
4 
5 template <typename Scalar>
6 void r1updt(
9  std::vector<JacobiRotation<Scalar> > &v_givens,
10  std::vector<JacobiRotation<Scalar> > &w_givens,
13  bool *sing)
14 {
15  typedef DenseIndex Index;
16  const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0);
17 
18  /* Local variables */
19  const Index m = s.rows();
20  const Index n = s.cols();
21  Index i, j=1;
22  Scalar temp;
24 
25  // r1updt had a broader usecase, but we dont use it here. And, more
26  // importantly, we can not test it.
27  eigen_assert(m==n);
28  eigen_assert(u.size()==m);
29  eigen_assert(v.size()==n);
30  eigen_assert(w.size()==n);
31 
32  /* move the nontrivial part of the last column of s into w. */
33  w[n-1] = s(n-1,n-1);
34 
35  /* rotate the vector v into a multiple of the n-th unit vector */
36  /* in such a way that a spike is introduced into w. */
37  for (j=n-2; j>=0; --j) {
38  w[j] = 0.;
39  if (v[j] != 0.) {
40  /* determine a givens rotation which eliminates the */
41  /* j-th element of v. */
42  givens.makeGivens(-v[n-1], v[j]);
43 
44  /* apply the transformation to v and store the information */
45  /* necessary to recover the givens rotation. */
46  v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
47  v_givens[j] = givens;
48 
49  /* apply the transformation to s and extend the spike in w. */
50  for (i = j; i < m; ++i) {
51  temp = givens.c() * s(j,i) - givens.s() * w[i];
52  w[i] = givens.s() * s(j,i) + givens.c() * w[i];
53  s(j,i) = temp;
54  }
55  } else
56  v_givens[j] = IdentityRotation;
57  }
58 
59  /* add the spike from the rank 1 update to w. */
60  w += v[n-1] * u;
61 
62  /* eliminate the spike. */
63  *sing = false;
64  for (j = 0; j < n-1; ++j) {
65  if (w[j] != 0.) {
66  /* determine a givens rotation which eliminates the */
67  /* j-th element of the spike. */
68  givens.makeGivens(-s(j,j), w[j]);
69 
70  /* apply the transformation to s and reduce the spike in w. */
71  for (i = j; i < m; ++i) {
72  temp = givens.c() * s(j,i) + givens.s() * w[i];
73  w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
74  s(j,i) = temp;
75  }
76 
77  /* store the information necessary to recover the */
78  /* givens rotation. */
79  w_givens[j] = givens;
80  } else
81  v_givens[j] = IdentityRotation;
82 
83  /* test for zero diagonal elements in the output s. */
84  if (s(j,j) == 0.) {
85  *sing = true;
86  }
87  }
88  /* move w back into the last column of the output s. */
89  s(n-1,n-1) = w[n-1];
90 
91  if (s(j,j) == 0.) {
92  *sing = true;
93  }
94  return;
95 }
96 
97 } // end namespace internal
98 
99 } // end namespace Eigen
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index cols() const
Matrix3f m
Scalar & c()
Definition: Jacobi.h:45
SCALAR Scalar
Definition: bench_gemm.cpp:33
void r1updt(Matrix< Scalar, Dynamic, Dynamic > &s, const Matrix< Scalar, Dynamic, 1 > &u, std::vector< JacobiRotation< Scalar > > &v_givens, std::vector< JacobiRotation< Scalar > > &w_givens, Matrix< Scalar, Dynamic, 1 > &v, Matrix< Scalar, Dynamic, 1 > &w, bool *sing)
Definition: r1updt.h:6
ArrayXcf v
Definition: Cwise_arg.cpp:1
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Index rows() const
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Rotation given by a cosine-sine pair.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
RealScalar s
RowVector3d w
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: Meta.h:25
Scalar & s()
Definition: Jacobi.h:47
std::ptrdiff_t j
void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:148


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:43:47