Go to the documentation of this file.00001 namespace internal {
00002
00003 template <typename Scalar>
00004 void r1updt(
00005 Matrix< Scalar, Dynamic, Dynamic > &s,
00006 const Matrix< Scalar, Dynamic, 1> &u,
00007 std::vector<JacobiRotation<Scalar> > &v_givens,
00008 std::vector<JacobiRotation<Scalar> > &w_givens,
00009 Matrix< Scalar, Dynamic, 1> &v,
00010 Matrix< Scalar, Dynamic, 1> &w,
00011 bool *sing)
00012 {
00013 typedef DenseIndex Index;
00014
00015
00016 const Index m = s.rows();
00017 const Index n = s.cols();
00018 Index i, j=1;
00019 Scalar temp;
00020 JacobiRotation<Scalar> givens;
00021
00022
00023
00024 assert(m==n);
00025 assert(u.size()==m);
00026 assert(v.size()==n);
00027 assert(w.size()==n);
00028
00029
00030 w[n-1] = s(n-1,n-1);
00031
00032
00033
00034 for (j=n-2; j>=0; --j) {
00035 w[j] = 0.;
00036 if (v[j] != 0.) {
00037
00038
00039 givens.makeGivens(-v[n-1], v[j]);
00040
00041
00042
00043 v[n-1] = givens.s() * v[j] + givens.c() * v[n-1];
00044 v_givens[j] = givens;
00045
00046
00047 for (i = j; i < m; ++i) {
00048 temp = givens.c() * s(j,i) - givens.s() * w[i];
00049 w[i] = givens.s() * s(j,i) + givens.c() * w[i];
00050 s(j,i) = temp;
00051 }
00052 }
00053 }
00054
00055
00056 w += v[n-1] * u;
00057
00058
00059 *sing = false;
00060 for (j = 0; j < n-1; ++j) {
00061 if (w[j] != 0.) {
00062
00063
00064 givens.makeGivens(-s(j,j), w[j]);
00065
00066
00067 for (i = j; i < m; ++i) {
00068 temp = givens.c() * s(j,i) + givens.s() * w[i];
00069 w[i] = -givens.s() * s(j,i) + givens.c() * w[i];
00070 s(j,i) = temp;
00071 }
00072
00073
00074
00075 w_givens[j] = givens;
00076 }
00077
00078
00079 if (s(j,j) == 0.) {
00080 *sing = true;
00081 }
00082 }
00083
00084 s(n-1,n-1) = w[n-1];
00085
00086 if (s(j,j) == 0.) {
00087 *sing = true;
00088 }
00089 return;
00090 }
00091
00092 }