| 1443 | } |
| 1444 | |
| 1445 | MatFreeDavidsonResult generalized_davidson_eigensolver( |
| 1446 | const HamiltonianOp& h_op, |
| 1447 | const HamiltonianOp& s_op, |
| 1448 | const std::vector<double>& kinetic_diag, |
| 1449 | int n_pw, |
| 1450 | double ecut_ha, |
| 1451 | const MatFreeDavidsonOptions& opts, |
| 1452 | const std::vector<Complex>& X0 |
| 1453 | ) { |
| 1454 | if (n_pw <= 0) { |
| 1455 | throw std::invalid_argument("generalized_davidson_eigensolver: n_pw must be positive"); |
| 1456 | } |
| 1457 | if (static_cast<int>(kinetic_diag.size()) != n_pw) { |
| 1458 | throw std::invalid_argument( |
| 1459 | "generalized_davidson_eigensolver: kinetic_diag size must match n_pw"); |
| 1460 | } |
| 1461 | |
| 1462 | const int n_bands = (opts.n_bands > 0) ? std::min(opts.n_bands, n_pw) : n_pw; |
| 1463 | if (n_pw <= opts.dense_fallback_threshold) { |
| 1464 | return dense_generalized_fallback(h_op, s_op, n_pw, n_bands); |
| 1465 | } |
| 1466 | |
| 1467 | const int subspace_factor = std::max(2, opts.max_subspace_factor); |
| 1468 | const int max_subspace = std::min(n_pw, std::max(n_bands, n_bands * subspace_factor)); |
| 1469 | const int max_iter = opts.max_iter; |
| 1470 | const double tol = opts.tol; |
| 1471 | |
| 1472 | std::vector<int> sorted_idx(static_cast<std::size_t>(n_pw)); |
| 1473 | std::iota(sorted_idx.begin(), sorted_idx.end(), 0); |
| 1474 | std::sort(sorted_idx.begin(), sorted_idx.end(), [&](int a, int b) { |
| 1475 | return kinetic_diag[static_cast<std::size_t>(a)] < kinetic_diag[static_cast<std::size_t>(b)]; |
| 1476 | }); |
| 1477 | |
| 1478 | int n_sub = 0; |
| 1479 | std::vector<Complex> V(static_cast<std::size_t>(n_pw) * max_subspace, Complex{0.0, 0.0}); |
| 1480 | std::vector<Complex> SV(static_cast<std::size_t>(n_pw) * max_subspace, Complex{0.0, 0.0}); |
| 1481 | std::vector<Complex> HV(static_cast<std::size_t>(n_pw) * max_subspace, Complex{0.0, 0.0}); |
| 1482 | |
| 1483 | if (!X0.empty() && static_cast<int>(X0.size()) == n_pw * n_bands) { |
| 1484 | for (int j = 0; j < n_bands; ++j) { |
| 1485 | for (int i = 0; i < n_pw; ++i) { |
| 1486 | V[static_cast<std::size_t>(i + n_pw * j)] = |
| 1487 | X0[static_cast<std::size_t>(i + n_pw * j)]; |
| 1488 | } |
| 1489 | } |
| 1490 | n_sub = s_orthonormalize_columns(V.data(), SV.data(), n_pw, n_bands, s_op); |
| 1491 | } else { |
| 1492 | for (int j = 0; j < n_bands; ++j) { |
| 1493 | V[static_cast<std::size_t>(sorted_idx[static_cast<std::size_t>(j)] + n_pw * j)] = |
| 1494 | Complex{1.0, 0.0}; |
| 1495 | } |
| 1496 | n_sub = s_orthonormalize_columns(V.data(), SV.data(), n_pw, n_bands, s_op); |
| 1497 | } |
| 1498 | |
| 1499 | if (n_sub < n_bands) { |
| 1500 | throw std::runtime_error( |
| 1501 | "generalized_davidson_eigensolver: could not form enough initial vectors"); |
| 1502 | } |