From cd0fbe7c2836028ef9700c9bd1a1a6f26f003e25 Mon Sep 17 00:00:00 2001 From: Zhentao Wang Date: Sun, 17 Sep 2017 01:12:08 -0400 Subject: [PATCH] updated examples. --- README.md | 25 +- examples/trans_absent/latt_chain/chain_tJ.cc | 7 +- .../trans_absent/latt_kagome/kagome_tJ.cc | 2 - .../trans_symmetric/latt_chain/chain_tJ.cc | 107 ++++++++ .../trans_symmetric/latt_kagome/kagome_tJ.cc | 237 ++++++++++++++++++ examples/trans_symmetric/object_list.txt | 3 + src/model.cc | 1 + 7 files changed, 366 insertions(+), 16 deletions(-) create mode 100644 examples/trans_symmetric/latt_chain/chain_tJ.cc create mode 100644 examples/trans_symmetric/latt_kagome/kagome_tJ.cc diff --git a/README.md b/README.md index 860c9a2..6d0c9f1 100644 --- a/README.md +++ b/README.md @@ -3,16 +3,21 @@ Basis of condensed matter quantum lattice problems, for usage in exact diagonali ## Examples To learn how to use this library to design ED code for your own models, please refer to the folder "examples": -- Heisenberg spin-1/2 chain -- Heisenberg spin-1 chain -- Heisenberg spin-1/2 triangular lattice -- Heisenberg spin-1/2 kagome lattice -- Spinless fermion on honeycomb lattice -- Fermi-Hubbard Model on square lattice -- Bose-Hubbard Model on square lattice -- Kondo Lattice Model on a chain -- tJ Model on chain lattice (E0 checked, but E1 inconsistent with ALPS) -- tJ Model on kagome lattice (E0 checked, but E1 inconsistent with ALPS) +- Chain + - Heisenberg spin-1/2 + - Heisenberg spin-1 + - Kondo Lattice model + - t-J model +- Honeycomb lattice + - Spinless fermion +- Kagome lattice + - Heisenberg spin-1/2 + - t-J model +- Square lattice + - Bose-Hubbard model + - Fermi-Hubbard model +- Triangular lattice + - Heisenberg spin-1/2 ## Dependencies: - [Boost](http://www.boost.org/) diff --git a/examples/trans_absent/latt_chain/chain_tJ.cc b/examples/trans_absent/latt_chain/chain_tJ.cc index 702767f..7ddf765 100644 --- a/examples/trans_absent/latt_chain/chain_tJ.cc +++ b/examples/trans_absent/latt_chain/chain_tJ.cc @@ -3,10 +3,8 @@ #include "qbasis.h" // tJ model on chain lattice -// E0 checked. Excitations seems not correct. int main() { qbasis::enable_ckpt = true; - std::cout << "ground state energy checked correct. 1st excited state inconsistent with ALPS!" << std::endl; std::cout << std::setprecision(10); // parameters double t = 1.0; @@ -93,10 +91,11 @@ int main() { // obtaining the eigenvals of the matrix - tJ.locate_E0_lanczos(0); + tJ.locate_E0_full(4,8); std::cout << std::endl; // for the parameters considered, we should obtain: - //assert(std::abs(tJ.eigenvals_full[0] + 15.41931496) < 1e-8); + assert(std::abs(tJ.eigenvals_full[0] + 9.762087307) < 1e-8); + assert(std::abs(tJ.eigenvals_full[1] + 9.762087307) < 1e-8); } diff --git a/examples/trans_absent/latt_kagome/kagome_tJ.cc b/examples/trans_absent/latt_kagome/kagome_tJ.cc index 28422d4..73188ba 100644 --- a/examples/trans_absent/latt_kagome/kagome_tJ.cc +++ b/examples/trans_absent/latt_kagome/kagome_tJ.cc @@ -3,10 +3,8 @@ #include "qbasis.h" // tJ model on kagome lattice -// Ground state energy consistent with ALPS, 1st excitation NOT agree. NEED further check. int main() { qbasis::enable_ckpt = true; - std::cout << "ground state energy checked correct. 1st excited state inconsistent with ALPS!" << std::endl; std::cout << std::setprecision(10); // parameters double t = 1.0; diff --git a/examples/trans_symmetric/latt_chain/chain_tJ.cc b/examples/trans_symmetric/latt_chain/chain_tJ.cc new file mode 100644 index 0000000..716cc7e --- /dev/null +++ b/examples/trans_symmetric/latt_chain/chain_tJ.cc @@ -0,0 +1,107 @@ +#include +#include +#include "qbasis.h" + +// tJ model on chain lattice +int main() { + qbasis::enable_ckpt = false; + std::cout << std::setprecision(10); + // parameters + double t = 1.0; + double J = 1.0; + int Lx = 12; + double N_total_val = 8; + double Sz_total_val = 0; + + std::cout << "Lx = " << Lx << std::endl; + std::cout << "t = " << t << std::endl; + std::cout << "J = " << J << std::endl; + std::cout << "N = " << N_total_val << std::endl; + std::cout << "Sz = " << Sz_total_val << std::endl << std::endl; + + + // lattice object + std::vector bc{"pbc"}; + qbasis::lattice lattice("chain",std::vector{static_cast(Lx)},bc); + + + // local matrix representation + auto c_up = std::vector>>(3,std::vector>(3, 0.0)); + auto c_dn = std::vector>>(3,std::vector>(3, 0.0)); + c_up[0][1] = std::complex(1.0,0.0); + c_dn[0][2] = std::complex(1.0,0.0); + + + // constructing the Hamiltonian in operator representation + qbasis::model> tJ; + tJ.add_orbital(lattice.total_sites(), "tJ"); + qbasis::mopr> Sz_total; // operators representating total Sz + qbasis::mopr> N_total; // operators representating total N + + for (int m = 0; m < Lx; m++) { + uint32_t site_i, site_j; + std::vector work(lattice.dimension()); + lattice.coor2site(std::vector{m}, 0, site_i, work); // obtain site label of (x,y) + lattice.coor2site(std::vector{m+1}, 0, site_j, work); + // construct operators on each site + auto c_up_i = qbasis::opr>(site_i,0,true,c_up); + auto c_dn_i = qbasis::opr>(site_i,0,true,c_dn); + auto c_up_dg_i = c_up_i; c_up_dg_i.dagger(); + auto c_dn_dg_i = c_dn_i; c_dn_dg_i.dagger(); + auto Splus_i = c_up_dg_i * c_dn_i; + auto Sminus_i = c_dn_dg_i * c_up_i; + auto Sz_i = std::complex(0.5,0.0) * (c_up_dg_i * c_up_i - c_dn_dg_i * c_dn_i); + auto N_i = (c_up_dg_i * c_up_i + c_dn_dg_i * c_dn_i); + + auto c_up_j = qbasis::opr>(site_j,0,true,c_up); + auto c_dn_j = qbasis::opr>(site_j,0,true,c_dn); + auto c_up_dg_j = c_up_j; c_up_dg_j.dagger(); + auto c_dn_dg_j = c_dn_j; c_dn_dg_j.dagger(); + auto Splus_j = c_up_dg_j * c_dn_j; + auto Sminus_j = c_dn_dg_j * c_up_j; + auto Sz_j = std::complex(0.5,0.0) * (c_up_dg_j * c_up_j - c_dn_dg_j * c_dn_j); + auto N_j = (c_up_dg_j * c_up_j + c_dn_dg_j * c_dn_j); + + if (bc[0] == "pbc" || (bc[0] == "obc" && m < Lx - 1)) { + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i * c_up_j )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_j * c_up_i )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i * c_dn_j )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_j * c_dn_i )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i * Sminus_j + Sminus_i * Splus_j)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i * Sz_j )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i * N_j )); + } + + + // total Sz operator + Sz_total += Sz_i; + // total N operator + N_total += N_i; + + } + + // to use translational symmetry, we first fill the Weisse tables + tJ.fill_Weisse_table(lattice); + + std::vector E0_list; + for (int momentum = 0; momentum < Lx; momentum++) { + // constructing the Hilbert space basis + tJ.enumerate_basis_repr({momentum}, {Sz_total, N_total}, {Sz_total_val, N_total_val}); + + // optional in future, will use more memory and give higher speed + // generating matrix of the Hamiltonian in the full Hilbert space + tJ.generate_Ham_sparse_repr(); + std::cout << std::endl; + + // obtaining the eigenvals of the matrix + tJ.locate_E0_lanczos(1); + std::cout << std::endl; + + E0_list.push_back(tJ.eigenvals_repr[0]); + + } + + for (int momentum = 0; momentum < Lx; momentum++) { + std::cout << "E0(k=" << momentum << ")=\t" << E0_list[momentum] << std::endl; + } +} diff --git a/examples/trans_symmetric/latt_kagome/kagome_tJ.cc b/examples/trans_symmetric/latt_kagome/kagome_tJ.cc new file mode 100644 index 0000000..da79954 --- /dev/null +++ b/examples/trans_symmetric/latt_kagome/kagome_tJ.cc @@ -0,0 +1,237 @@ +#include +#include +#include "qbasis.h" + +// tJ model on kagome lattice +int main() { + qbasis::enable_ckpt = false; + std::cout << std::setprecision(10); + // parameters + double t = 1.0; + double J = 1.0; + int Lx = 2; + int Ly = 2; + double N_total_val = Lx * Ly * 2; + double Sz_total_val = 0; + + std::cout << "Lx = " << Lx << std::endl; + std::cout << "Ly = " << Ly << std::endl; + std::cout << "t = " << t << std::endl; + std::cout << "J = " << J << std::endl; + std::cout << "N = " << N_total_val << std::endl; + std::cout << "Sz = " << Sz_total_val << std::endl << std::endl; + + + // lattice object + std::vector bc{"pbc", "pbc"}; + assert(bc[0] == "pbc" && bc[1] == "pbc"); // "obc" requires more careful job in the following + qbasis::lattice lattice("kagome",std::vector{static_cast(Lx), static_cast(Ly)},bc); + + + // local matrix representation + auto c_up = std::vector>>(3,std::vector>(3, 0.0)); + auto c_dn = std::vector>>(3,std::vector>(3, 0.0)); + c_up[0][1] = std::complex(1.0,0.0); + c_dn[0][2] = std::complex(1.0,0.0); + + + // constructing the Hamiltonian in operator representation + qbasis::model> tJ; + tJ.add_orbital(lattice.total_sites(), "tJ"); + qbasis::mopr> Sz_total; // operators representating total Sz + qbasis::mopr> N_total; // operators representating total N + + for (int m = 0; m < Lx; m++) { + for (int n = 0; n < Ly; n++) { + uint32_t site_i0, site_i1, site_i2; + std::vector work(lattice.dimension()); + lattice.coor2site(std::vector{m,n}, 0, site_i0, work); // obtain site label of (x,y) + lattice.coor2site(std::vector{m,n}, 1, site_i1, work); + lattice.coor2site(std::vector{m,n}, 2, site_i2, work); + // construct operators on each site + auto c_up_i0 = qbasis::opr>(site_i0,0,true,c_up); + auto c_dn_i0 = qbasis::opr>(site_i0,0,true,c_dn); + auto c_up_dg_i0 = c_up_i0; c_up_dg_i0.dagger(); + auto c_dn_dg_i0 = c_dn_i0; c_dn_dg_i0.dagger(); + auto Splus_i0 = c_up_dg_i0 * c_dn_i0; + auto Sminus_i0 = c_dn_dg_i0 * c_up_i0; + auto Sz_i0 = std::complex(0.5,0.0) * (c_up_dg_i0 * c_up_i0 - c_dn_dg_i0 * c_dn_i0); + auto N_i0 = (c_up_dg_i0 * c_up_i0 + c_dn_dg_i0 * c_dn_i0); + + auto c_up_i1 = qbasis::opr>(site_i1,0,true,c_up); + auto c_dn_i1 = qbasis::opr>(site_i1,0,true,c_dn); + auto c_up_dg_i1 = c_up_i1; c_up_dg_i1.dagger(); + auto c_dn_dg_i1 = c_dn_i1; c_dn_dg_i1.dagger(); + auto Splus_i1 = c_up_dg_i1 * c_dn_i1; + auto Sminus_i1 = c_dn_dg_i1 * c_up_i1; + auto Sz_i1 = std::complex(0.5,0.0) * (c_up_dg_i1 * c_up_i1 - c_dn_dg_i1 * c_dn_i1); + auto N_i1 = (c_up_dg_i1 * c_up_i1 + c_dn_dg_i1 * c_dn_i1); + + auto c_up_i2 = qbasis::opr>(site_i2,0,true,c_up); + auto c_dn_i2 = qbasis::opr>(site_i2,0,true,c_dn); + auto c_up_dg_i2 = c_up_i2; c_up_dg_i2.dagger(); + auto c_dn_dg_i2 = c_dn_i2; c_dn_dg_i2.dagger(); + auto Splus_i2 = c_up_dg_i2 * c_dn_i2; + auto Sminus_i2 = c_dn_dg_i2 * c_up_i2; + auto Sz_i2 = std::complex(0.5,0.0) * (c_up_dg_i2 * c_up_i2 - c_dn_dg_i2 * c_dn_i2); + auto N_i2 = (c_up_dg_i2 * c_up_i2 + c_dn_dg_i2 * c_dn_i2); + + // 1 -- 0 -> 1 + { + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i0 * c_up_i1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i1 * c_up_i0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i0 * c_dn_i1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i1 * c_dn_i0 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i0 * Sminus_i1 + Sminus_i0 * Splus_i1)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i0 * Sz_i1 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i0 * N_i1 )); + } + + + // 1 <- 0 -- 1 + { + uint32_t site_j1; + lattice.coor2site(std::vector{m-1,n}, 1, site_j1, work); + auto c_up_j1 = qbasis::opr>(site_j1,0,true,c_up); + auto c_dn_j1 = qbasis::opr>(site_j1,0,true,c_dn); + auto c_up_dg_j1 = c_up_j1; c_up_dg_j1.dagger(); + auto c_dn_dg_j1 = c_dn_j1; c_dn_dg_j1.dagger(); + auto Splus_j1 = c_up_dg_j1 * c_dn_j1; + auto Sminus_j1 = c_dn_dg_j1 * c_up_j1; + auto Sz_j1 = std::complex(0.5,0.0) * (c_up_dg_j1 * c_up_j1 - c_dn_dg_j1 * c_dn_j1); + auto N_j1 = (c_up_dg_j1 * c_up_j1 + c_dn_dg_j1 * c_dn_j1); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i0 * c_up_j1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_j1 * c_up_i0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i0 * c_dn_j1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_j1 * c_dn_i0 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i0 * Sminus_j1 + Sminus_i0 * Splus_j1)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i0 * Sz_j1 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i0 * N_j1 )); + } + + + // 2 + // ^ + // \ + // 1 + // \ + // \ + // 2 + { + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i1 * c_up_i2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i2 * c_up_i1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i1 * c_dn_i2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i2 * c_dn_i1 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i1 * Sminus_i2 + Sminus_i1 * Splus_i2)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i1 * Sz_i2 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i1 * N_i2 )); + } + + + // 2 + // \ + // \ + // 1 + // \ + // v + // 2 + { + uint32_t site_j2; + lattice.coor2site(std::vector{m+1,n-1}, 2, site_j2, work); + auto c_up_j2 = qbasis::opr>(site_j2,0,true,c_up); + auto c_dn_j2 = qbasis::opr>(site_j2,0,true,c_dn); + auto c_up_dg_j2 = c_up_j2; c_up_dg_j2.dagger(); + auto c_dn_dg_j2 = c_dn_j2; c_dn_dg_j2.dagger(); + auto Splus_j2 = c_up_dg_j2 * c_dn_j2; + auto Sminus_j2 = c_dn_dg_j2 * c_up_j2; + auto Sz_j2 = std::complex(0.5,0.0) * (c_up_dg_j2 * c_up_j2 - c_dn_dg_j2 * c_dn_j2); + auto N_j2 = (c_up_dg_j2 * c_up_j2 + c_dn_dg_j2 * c_dn_j2); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i1 * c_up_j2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_j2 * c_up_i1 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i1 * c_dn_j2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_j2 * c_dn_i1 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i1 * Sminus_j2 + Sminus_i1 * Splus_j2)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i1 * Sz_j2 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i1 * N_j2 )); + } + + + // 0 + // / + // / + // 2 + // / + // v + // 0 + { + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i2 * c_up_i0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i0 * c_up_i2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i2 * c_dn_i0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i0 * c_dn_i2 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i2 * Sminus_i0 + Sminus_i2 * Splus_i0)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i2 * Sz_i0 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i2 * N_i0 )); + } + + + // 0 + // ^ + // / + // 2 + // / + // / + // 0 + { + uint32_t site_j0; + lattice.coor2site(std::vector{m,n+1}, 0, site_j0, work); + auto c_up_j0 = qbasis::opr>(site_j0,0,true,c_up); + auto c_dn_j0 = qbasis::opr>(site_j0,0,true,c_dn); + auto c_up_dg_j0 = c_up_j0; c_up_dg_j0.dagger(); + auto c_dn_dg_j0 = c_dn_j0; c_dn_dg_j0.dagger(); + auto Splus_j0 = c_up_dg_j0 * c_dn_j0; + auto Sminus_j0 = c_dn_dg_j0 * c_up_j0; + auto Sz_j0 = std::complex(0.5,0.0) * (c_up_dg_j0 * c_up_j0 - c_dn_dg_j0 * c_dn_j0); + auto N_j0 = (c_up_dg_j0 * c_up_j0 + c_dn_dg_j0 * c_dn_j0); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_i2 * c_up_j0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_up_dg_j0 * c_up_i2 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_i2 * c_dn_j0 )); + tJ.add_offdiagonal_Ham(std::complex(-t,0.0) * ( c_dn_dg_j0 * c_dn_i2 )); + tJ.add_offdiagonal_Ham(std::complex(0.5*J,0.0) * (Splus_i2 * Sminus_j0 + Sminus_i2 * Splus_j0)); + tJ.add_diagonal_Ham(std::complex(J,0.0) * ( Sz_i2 * Sz_j0 )); + tJ.add_diagonal_Ham(std::complex(-0.25*J,0.0) * ( N_i2 * N_j0 )); + } + + // total Sz operator + Sz_total += (Sz_i0 + Sz_i1 + Sz_i2); + // total N operator + N_total += (N_i0 + N_i1 + N_i2); + } + } + + // to use translational symmetry, we first fill the Weisse tables + tJ.fill_Weisse_table(lattice); + + std::vector E0_list; + for (int m = 0; m < Lx; m++) { + for (int n = 0; n < Ly; n++) { + // constructing the Hilbert space basis + tJ.enumerate_basis_repr({m,n}, {Sz_total, N_total}, {Sz_total_val, N_total_val}); + + // generating matrix of the Hamiltonian in the subspace + tJ.generate_Ham_sparse_repr(); + std::cout << std::endl; + + // obtaining the lowest eigenvals of the matrix + tJ.locate_E0_lanczos(1); + std::cout << std::endl; + + E0_list.push_back(tJ.eigenvals_repr[0]); + } + } + + // for the parameters considered, we should obtain: + assert(std::abs(E0_list[0] + 15.41931496) < 1e-8); + assert(std::abs(E0_list[1] + 14.40277723) < 1e-8); + assert(std::abs(E0_list[2] + 14.40277723) < 1e-8); + assert(std::abs(E0_list[3] + 14.40277723) < 1e-8); +} diff --git a/examples/trans_symmetric/object_list.txt b/examples/trans_symmetric/object_list.txt index a8f233f..ba95d00 100644 --- a/examples/trans_symmetric/object_list.txt +++ b/examples/trans_symmetric/object_list.txt @@ -1,11 +1,14 @@ vpath %.cc ../latt_chain \ ../latt_honeycomb \ + ../latt_kagome \ ../latt_square \ ../latt_triangular EXEC = chain_Heisenberg_spin_half.x \ chain_Heisenberg_spin_one.x \ + chain_tJ.x \ chain_Kondo.x \ honeycomb_Spinless_Fermion.x \ + kagome_tJ.x \ square_Fermi_Hubbard.x \ triangular_Heisenberg_spin_half.x diff --git a/src/model.cc b/src/model.cc index 10065f8..571ed79 100644 --- a/src/model.cc +++ b/src/model.cc @@ -808,6 +808,7 @@ namespace qbasis { if (nev == 2 && ! E1_done) { start = std::chrono::system_clock::now(); std::cout << "Calculating 1st excited state energy (simple Lanczos)..." << std::endl; + std::cout << "Lanczos may/maynot miss degenerate states, BE CAREFUL!" << std::endl; vec_randomize(dim, v.data(), seed); auto alpha = dotc(dim, v.data() + 2 * dim, 1, v.data(), 1); // (phi0, v0) axpy(dim, -alpha, v.data() + 2 * dim, 1, v.data(), 1); // v0 -= alpha * phi0