From f082838c5f0482189e5c97ed14755b3f3929e5a5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 21 Jan 2020 15:29:39 +0200 Subject: [PATCH] Beyn algorithm "cherry-pick" from 'newbeyn_unitcell' - Add rank_min_sel argument to beyn_solve() and beyn_solve_gsl() - Fix order of K and K_coarse evaluation (K_coarse should probably be removed). Former-commit-id: c0a241f8712439ba84e7c907658ebb6071528482 --- qpms/beyn.c | 22 +++++++++++----------- qpms/beyn.h | 6 ++++-- qpms/qpms_cdefs.pxd | 4 ++-- 3 files changed, 17 insertions(+), 15 deletions(-) diff --git a/qpms/beyn.c b/qpms/beyn.c index 28911be..b04649f 100644 --- a/qpms/beyn.c +++ b/qpms/beyn.c @@ -345,7 +345,7 @@ void BeynSolver_srandom(BeynSolver *solver, unsigned int RandSeed) static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_function, void *Params, gsl_matrix_complex *A0, gsl_matrix_complex *A1, double complex z0, - gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, const double res_tol) + gsl_vector_complex *eigenvalues, gsl_matrix_complex *eigenvectors, const double rank_tol, size_t rank_sel_min, const double res_tol) { const size_t m = solver->M; const size_t l = solver->L; @@ -379,7 +379,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun int K=0; for (int k=0; ksize /* this is l, actually */; k++) { if (verbose) printf("Beyn: SV(%d)=%e\n",k,gsl_vector_get(Sigma, k)); - if (gsl_vector_get(Sigma, k) > rank_tol) + if (k < rank_sel_min || gsl_vector_get(Sigma, k) > rank_tol) K++; } if (verbose)printf(" Beyn: %d/%zd relevant singular values\n",K,l); @@ -494,7 +494,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun gsl_vector_complex_set(eigenvalues, KRetained, zgsl); if(eigenvectors) { - gsl_matrix_complex_set_col(eigenvectors, KRetained, &(Vk.vector)); + gsl_matrix_complex_set_row(eigenvectors, KRetained, &(Vk.vector)); gsl_vector_set(solver->residuals, KRetained, residual); } ++KRetained; @@ -513,7 +513,7 @@ static int beyn_process_matrices(BeynSolver *solver, beyn_function_M_gsl_t M_fun beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, beyn_function_M_gsl_t M_function, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat_function, void *params, const beyn_contour_t *contour, - double rank_tol, double res_tol) + double rank_tol, size_t rank_sel_min, double res_tol) { BeynSolver *solver = BeynSolver_create(m, l); @@ -583,14 +583,14 @@ beyn_result_gsl_t *beyn_solve_gsl(const size_t m, const size_t l, gsl_vector_complex *eigenvalue_errors = solver->eigenvalue_errors; gsl_matrix_complex *eigenvectors = solver->eigenvectors; - // Beyn Steps 3–6 - int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, res_tol); // Repeat Steps 3–6 with rougher contour approximation to get an error estimate. - int K_coarse = beyn_process_matrices(solver, M_function, params, A0_coarse, A1_coarse, z0, eigenvalue_errors, eigenvectors, rank_tol, res_tol); - - gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); + int K_coarse = beyn_process_matrices(solver, M_function, params, A0_coarse, A1_coarse, z0, eigenvalue_errors, /*eigenvectors_coarse*/ NULL, rank_tol, rank_sel_min, res_tol); // Reid did also fabs on the complex and real parts ^^^. + // Beyn Steps 3–6 + int K = beyn_process_matrices(solver, M_function, params, A0, A1, z0, eigenvalues, eigenvectors, rank_tol, rank_sel_min, res_tol); + gsl_blas_zaxpy(gsl_complex_rect(-1,0), eigenvalues, eigenvalue_errors); + beyn_result_gsl_t *result; QPMS_CRASHING_MALLOC(result, sizeof(beyn_result_gsl_t)); result->eigval = solver->eigenvalues; @@ -636,12 +636,12 @@ static int beyn_function_M_inv_Vhat_carr2gsl(gsl_matrix_complex *target, } beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat, - void *params, const beyn_contour_t *contour, double rank_tol, double res_tol) { + void *params, const beyn_contour_t *contour, double rank_tol, size_t rank_sel_min, double res_tol) { struct beyn_function_M_carr2gsl_param p = {M, M_inv_Vhat, params}; return beyn_result_from_beyn_result_gsl( beyn_solve_gsl(m, l, beyn_function_M_carr2gsl, (p.M_inv_Vhat_function) ? beyn_function_M_inv_Vhat_carr2gsl : NULL, - (void *) &p, contour, rank_tol, res_tol) + (void *) &p, contour, rank_tol, rank_sel_min, res_tol) ); } diff --git a/qpms/beyn.h b/qpms/beyn.h index 25bc1c7..86f7702 100644 --- a/qpms/beyn.h +++ b/qpms/beyn.h @@ -72,7 +72,7 @@ typedef struct beyn_result_gsl_t { gsl_vector_complex *eigval; gsl_vector_complex *eigval_err; gsl_vector *residuals; - gsl_matrix_complex *eigvec; + gsl_matrix_complex *eigvec; // Rows are the eigenvectors gsl_vector *ranktest_SV; } beyn_result_gsl_t; @@ -85,7 +85,7 @@ typedef struct beyn_result_t { complex double *eigval; complex double *eigval_err; double *residuals; - complex double *eigvec; + complex double *eigvec; // Rows are the eigenvectors double *ranktest_SV; /// Private, we wrap it around beyn_result_gsl_t for now. @@ -109,6 +109,7 @@ beyn_result_gsl_t *beyn_solve_gsl( void *params, ///< Parameter pointer passed to M() and M_inv_Vhat(). const beyn_contour_t *contour, ///< Integration contour. double rank_tol, ///< (default: `1e-4`) TODO DOC. + size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol. double res_tol ///< (default: `0.0`) TODO DOC. ); @@ -120,6 +121,7 @@ beyn_result_t *beyn_solve( void *params, ///< Parameter pointer passed to M() and M_inv_Vhat(). const beyn_contour_t *contour, ///< Integration contour. double rank_tol, ///< (default: `1e-4`) TODO DOC. + size_t rank_min_sel, ///< Minimum number of eigenvalue candidates, even if they don't pass \a rank_tol. double res_tol ///< (default: `0.0`) TODO DOC. ); diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 6df0435..a0d892a 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -674,11 +674,11 @@ cdef extern from "beyn.h": beyn_result_gsl_t *beyn_solve_gsl(size_t m, size_t l, beyn_function_M_gsl_t M, beyn_function_M_inv_Vhat_gsl_t M_inv_Vhat, void *params, const beyn_contour_t *contour, - double rank_tol, double res_tol) + double rank_tol, size_t rank_min_sel, double res_tol) beyn_result_t *beyn_solve(size_t m, size_t l, beyn_function_M_t M, beyn_function_M_inv_Vhat_t M_inv_Vhat, void *params, const beyn_contour_t *contour, - double rank_tol, double res_tol) + double rank_tol, size_t rank_min_sel, double res_tol) beyn_contour_t *beyn_contour_ellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints) beyn_contour_t *beyn_contour_halfellipse(cdouble centre, double halfax_re, double halfax_im, size_t npoints,