$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 namespace Eigen { 00002 00003 namespace internal { 00004 00005 template <typename Scalar> 00006 void r1updt( 00007 Matrix< Scalar, Dynamic, Dynamic > &s, 00008 const Matrix< Scalar, Dynamic, 1> &u, 00009 std::vector<JacobiRotation<Scalar> > &v_givens, 00010 std::vector<JacobiRotation<Scalar> > &w_givens, 00011 Matrix< Scalar, Dynamic, 1> &v, 00012 Matrix< Scalar, Dynamic, 1> &w, 00013 bool *sing) 00014 { 00015 typedef DenseIndex Index; 00016 const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1,0); 00017 00018 /* Local variables */ 00019 const Index m = s.rows(); 00020 const Index n = s.cols(); 00021 Index i, j=1; 00022 Scalar temp; 00023 JacobiRotation<Scalar> givens; 00024 00025 // r1updt had a broader usecase, but we dont use it here. And, more 00026 // importantly, we can not test it. 00027 eigen_assert(m==n); 00028 eigen_assert(u.size()==m); 00029 eigen_assert(v.size()==n); 00030 eigen_assert(w.size()==n); 00031 00032 /* move the nontrivial part of the last column of s into w. */ 00033 w[n-1] = s(n-1,n-1); 00034 00035 /* rotate the vector v into a multiple of the n-th unit vector */ 00036 /* in such a way that a spike is introduced into w. */ 00037 for (j=n-2; j>=0; --j) { 00038 w[j] = 0.; 00039 if (v[j] != 0.) { 00040 /* determine a givens rotation which eliminates the */ 00041 /* j-th element of v. */ 00042 givens.makeGivens(-v[n-1], v[j]); 00043 00044 /* apply the transformation to v and store the information */ 00045 /* necessary to recover the givens rotation. */ 00046 v[n-1] = givens.s() * v[j] + givens.c() * v[n-1]; 00047 v_givens[j] = givens; 00048 00049 /* apply the transformation to s and extend the spike in w. */ 00050 for (i = j; i < m; ++i) { 00051 temp = givens.c() * s(j,i) - givens.s() * w[i]; 00052 w[i] = givens.s() * s(j,i) + givens.c() * w[i]; 00053 s(j,i) = temp; 00054 } 00055 } else 00056 v_givens[j] = IdentityRotation; 00057 } 00058 00059 /* add the spike from the rank 1 update to w. */ 00060 w += v[n-1] * u; 00061 00062 /* eliminate the spike. */ 00063 *sing = false; 00064 for (j = 0; j < n-1; ++j) { 00065 if (w[j] != 0.) { 00066 /* determine a givens rotation which eliminates the */ 00067 /* j-th element of the spike. */ 00068 givens.makeGivens(-s(j,j), w[j]); 00069 00070 /* apply the transformation to s and reduce the spike in w. */ 00071 for (i = j; i < m; ++i) { 00072 temp = givens.c() * s(j,i) + givens.s() * w[i]; 00073 w[i] = -givens.s() * s(j,i) + givens.c() * w[i]; 00074 s(j,i) = temp; 00075 } 00076 00077 /* store the information necessary to recover the */ 00078 /* givens rotation. */ 00079 w_givens[j] = givens; 00080 } else 00081 v_givens[j] = IdentityRotation; 00082 00083 /* test for zero diagonal elements in the output s. */ 00084 if (s(j,j) == 0.) { 00085 *sing = true; 00086 } 00087 } 00088 /* move w back into the last column of the output s. */ 00089 s(n-1,n-1) = w[n-1]; 00090 00091 if (s(j,j) == 0.) { 00092 *sing = true; 00093 } 00094 return; 00095 } 00096 00097 } // end namespace internal 00098 00099 } // end namespace Eigen