17 using namespace Eigen;
33 template <
typename MatrixType>
40 int saRepeats = stdRepeats * 4;
46 SquareMatrixType covMat =
a *
a.adjoint();
51 int r = internal::random<int>(0,covMat.rows()-1);
52 int c = internal::random<int>(0,covMat.cols()-1);
58 for (
int k=0; k<saRepeats; ++k)
72 for (
int k=0; k<stdRepeats; ++k)
81 if (MatrixType::RowsAtCompileTime==
Dynamic)
84 std::cout <<
"fixed ";
85 std::cout << covMat.rows() <<
" \t"
90 if (MatrixType::RowsAtCompileTime==
Dynamic)
95 gmm::dense_matrix<Scalar> gmmCovMat(covMat.rows(),covMat.cols());
96 gmm::dense_matrix<Scalar> eigvect(covMat.rows(),covMat.cols());
97 std::vector<Scalar> eigval(covMat.rows());
102 for (
int k=0; k<saRepeats; ++k)
104 gmm::symmetric_qr_algorithm(gmmCovMat, eigval, eigvect);
128 if (MatrixType::RowsAtCompileTime==
Dynamic)
133 gsl_matrix* gslCovMat = gsl_matrix_alloc(covMat.rows(),covMat.cols());
134 gsl_matrix* gslCopy = gsl_matrix_alloc(covMat.rows(),covMat.cols());
135 gsl_matrix* eigvect = gsl_matrix_alloc(covMat.rows(),covMat.cols());
136 gsl_vector* eigval = gsl_vector_alloc(covMat.rows());
137 gsl_eigen_symmv_workspace* eisymm = gsl_eigen_symmv_alloc(covMat.rows());
139 gsl_matrix_complex* eigvectz = gsl_matrix_complex_alloc(covMat.rows(),covMat.cols());
140 gsl_vector_complex* eigvalz = gsl_vector_complex_alloc(covMat.rows());
141 gsl_eigen_nonsymmv_workspace* einonsymm = gsl_eigen_nonsymmv_alloc(covMat.rows());
143 eiToGsl(covMat, &gslCovMat);
147 for (
int k=0; k<saRepeats; ++k)
149 gsl_matrix_memcpy(gslCopy,gslCovMat);
150 gsl_eigen_symmv(gslCopy, eigval, eigvect, eisymm);
151 acc += gsl_matrix_get(eigvect,r,
c);
158 for (
int k=0; k<stdRepeats; ++k)
160 gsl_matrix_memcpy(gslCopy,gslCovMat);
161 gsl_eigen_nonsymmv(gslCopy, eigvalz, eigvectz, einonsymm);
162 acc += GSL_REAL(gsl_matrix_complex_get(eigvectz,r,
c));
171 gsl_matrix_free(gslCovMat);
172 gsl_vector_free(gslCopy);
173 gsl_matrix_free(eigvect);
174 gsl_vector_free(eigval);
175 gsl_matrix_complex_free(eigvectz);
176 gsl_vector_complex_free(eigvalz);
177 gsl_eigen_symmv_free(eisymm);
178 gsl_eigen_nonsymmv_free(einonsymm);
189 int main(
int argc,
char* argv[])
191 const int dynsizes[] = {4,6,8,12,16,24,32,64,128,256,512,0};
192 std::cout <<
"size selfadjoint generic";
194 std::cout <<
" GMM++ ";
197 std::cout <<
" GSL (double + ATLAS) ";
200 for (uint
i=0; dynsizes[
i]>0; ++
i)