25 int svd_eigen_Macie(
const Eigen::MatrixXd& A,Eigen::MatrixXd& U,Eigen::VectorXd& S, Eigen::MatrixXd& V,
26 Eigen::MatrixXd& B, Eigen::VectorXd& tempi,
27 double threshold,
bool toggle)
30 unsigned int sweeps=0;
31 unsigned int rotations=0;
39 for(
unsigned int i=0;i<B.cols();i++){
40 for(
unsigned int j=i+1;j<B.cols();j++){
42 double p = B.col(i).dot(B.col(j));
43 double qi =B.col(i).dot(B.col(i));
44 double qj = B.col(j).dot(B.col(j));
46 double alpha =
pow(p,2.0)/(qi*qj);
55 cos=
sqrt((c+q)/(2*c));
59 sin=-
sqrt((c-q)/(2*c));
61 sin=
sqrt((c-q)/(2*c));
66 tempi = cos*B.col(i) + sin*B.col(j);
67 B.col(j) = - sin*B.col(i) + cos*B.col(j);
71 tempi.head(V.rows()) = cos*V.col(i) + sin*V.col(j);
72 V.col(j) = - sin*V.col(i) + cos*V.col(j);
73 V.col(i) = tempi.head(V.rows());
80 for(
unsigned int i=0;i<U.rows();i++) {
82 double si=
sqrt(B.col(i).dot(B.col(i)));
86 U.col(i) = B.col(i)/si;
98 B = U.transpose().lazyProduct(A);
103 for(
unsigned int i=0;i<B.cols();i++){
104 for(
unsigned int j=i+1;j<B.cols();j++){
106 double p = B.row(i).dot(B.row(j));
107 double qi = B.row(i).dot(B.row(i));
108 double qj = B.row(j).dot(B.row(j));
111 double alpha =
pow(p,2.0)/(qi*qj);
121 cos=
sqrt((c+q)/(2*c));
125 sin=-
sqrt((c-q)/(2*c));
127 sin=
sqrt((c-q)/(2*c));
132 tempi.head(B.cols()) = cos*B.row(i) + sin*B.row(j);
133 B.row(j) = - sin*B.row(i) + cos*B.row(j);
134 B.row(i) = tempi.head(B.cols());
138 tempi.head(U.rows()) = cos*U.col(i) + sin*U.col(j);
139 U.col(j) = - sin*U.col(i) + cos*U.col(j);
140 U.col(i) = tempi.head(U.rows());
148 for(
unsigned int i=0;i<V.rows();i++) {
149 double si=
sqrt(B.row(i).dot(B.row(i)));
153 V.col(i) = B.row(i)/si;
int svd_eigen_Macie(const Eigen::MatrixXd &A, Eigen::MatrixXd &U, Eigen::VectorXd &S, Eigen::MatrixXd &V, Eigen::MatrixXd &B, Eigen::VectorXd &tempi, double threshold, bool toggle)
INLINE Rall1d< T, V, S > sqrt(const Rall1d< T, V, S > &arg)
INLINE Rall1d< T, V, S > pow(const Rall1d< T, V, S > &arg, double m)
INLINE Rall1d< T, V, S > cos(const Rall1d< T, V, S > &arg)
INLINE Rall1d< T, V, S > sin(const Rall1d< T, V, S > &arg)