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);
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);
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;