|
1 | 1 | #include "atom_pair.h" |
2 | 2 | #include <complex> |
3 | 3 | #include <cassert> |
| 4 | +#include "source_base/module_external/blas_connector.h" |
4 | 5 |
|
5 | 6 | namespace hamilt |
6 | 7 | { |
@@ -585,6 +586,261 @@ void AtomPair<T>::merge_to_gamma() |
585 | 586 | this->values.push_back(tmp); |
586 | 587 | } |
587 | 588 |
|
| 589 | +template <typename T> |
| 590 | +void AtomPair<T>::add_from_matrix(const std::complex<T>* hk, |
| 591 | + const int ld_hk, |
| 592 | + const std::complex<T>& kphase, |
| 593 | + const int hk_type) |
| 594 | +{ |
| 595 | + // Only works for R_index[0], which is (0,0,0) |
| 596 | + const BaseMatrix<T>& matrix = values[0]; |
| 597 | + T* hr_tmp = matrix.get_pointer(); |
| 598 | + const std::complex<T>* hk_tmp = hk; |
| 599 | + const T* hk_real_pointer = nullptr; |
| 600 | + const T* hk_imag_pointer = nullptr; |
| 601 | + const int ld_hk_2 = ld_hk * 2; |
| 602 | + // row major |
| 603 | + if (hk_type == 0) |
| 604 | + { |
| 605 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 606 | + for (int mu = 0; mu < this->row_size; mu++) |
| 607 | + { |
| 608 | + hk_real_pointer = (T*)hk_tmp; |
| 609 | + hk_imag_pointer = hk_real_pointer+1; |
| 610 | + BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, 2, hr_tmp, 1); |
| 611 | + BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, 2, hr_tmp, 1); |
| 612 | + hk_tmp += ld_hk; |
| 613 | + hr_tmp += this->col_size; |
| 614 | + } |
| 615 | + } |
| 616 | + // column major |
| 617 | + else if (hk_type == 1) |
| 618 | + { |
| 619 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 620 | + for (int mu = 0; mu < this->row_size; mu++) |
| 621 | + { |
| 622 | + hk_real_pointer = (T*)hk_tmp; |
| 623 | + hk_imag_pointer = hk_real_pointer+1; |
| 624 | + BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, ld_hk_2, hr_tmp, 1); |
| 625 | + BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, ld_hk_2, hr_tmp, 1); |
| 626 | + hk_tmp ++; |
| 627 | + hr_tmp += this->col_size; |
| 628 | + } |
| 629 | + } |
| 630 | +} |
| 631 | + |
| 632 | +// add_from_matrix |
| 633 | +template <typename T> |
| 634 | +void AtomPair<T>::add_from_matrix(const T* hk, const int ld_hk, const T& kphase, const int hk_type) |
| 635 | +{ |
| 636 | + // Only works for R_index[0], which is (0,0,0) |
| 637 | + const BaseMatrix<T>& matrix = values[0]; |
| 638 | + T* hr_tmp = matrix.get_pointer(); |
| 639 | + const T* hk_tmp = hk; |
| 640 | + // row major |
| 641 | + if (hk_type == 0) |
| 642 | + { |
| 643 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 644 | + for (int mu = 0; mu < this->row_size; mu++) |
| 645 | + { |
| 646 | + BlasConnector::axpy(this->col_size, kphase, hk_tmp, 1, hr_tmp, 1); |
| 647 | + hk_tmp += ld_hk; |
| 648 | + hr_tmp += this->col_size; |
| 649 | + } |
| 650 | + } |
| 651 | + // column major |
| 652 | + else if (hk_type == 1) |
| 653 | + { |
| 654 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 655 | + for (int mu = 0; mu < this->row_size; mu++) |
| 656 | + { |
| 657 | + BlasConnector::axpy(this->col_size, kphase, hk_tmp, ld_hk, hr_tmp, 1); |
| 658 | + ++hk_tmp; |
| 659 | + hr_tmp += this->col_size; |
| 660 | + } |
| 661 | + } |
| 662 | +} |
| 663 | + |
| 664 | +// add_from_matrix with explicit R_index - thread-safe version |
| 665 | +template <typename T> |
| 666 | +void AtomPair<T>::add_from_matrix(const int R_index, |
| 667 | + const std::complex<T>* hk, |
| 668 | + const int ld_hk, |
| 669 | + const std::complex<T>& kphase, |
| 670 | + const int hk_type) |
| 671 | +{ |
| 672 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 673 | + T* hr_tmp = matrix.get_pointer(); |
| 674 | + const std::complex<T>* hk_tmp = hk; |
| 675 | + const T* hk_real_pointer = nullptr; |
| 676 | + const T* hk_imag_pointer = nullptr; |
| 677 | + const int ld_hk_2 = ld_hk * 2; |
| 678 | + // row major |
| 679 | + if (hk_type == 0) |
| 680 | + { |
| 681 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 682 | + for (int mu = 0; mu < this->row_size; mu++) |
| 683 | + { |
| 684 | + hk_real_pointer = (T*)hk_tmp; |
| 685 | + hk_imag_pointer = hk_real_pointer+1; |
| 686 | + BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, 2, hr_tmp, 1); |
| 687 | + BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, 2, hr_tmp, 1); |
| 688 | + hk_tmp += ld_hk; |
| 689 | + hr_tmp += this->col_size; |
| 690 | + } |
| 691 | + } |
| 692 | + // column major |
| 693 | + else if (hk_type == 1) |
| 694 | + { |
| 695 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 696 | + for (int mu = 0; mu < this->row_size; mu++) |
| 697 | + { |
| 698 | + hk_real_pointer = (T*)hk_tmp; |
| 699 | + hk_imag_pointer = hk_real_pointer+1; |
| 700 | + BlasConnector::axpy(this->col_size, kphase.real(), hk_real_pointer, ld_hk_2, hr_tmp, 1); |
| 701 | + BlasConnector::axpy(this->col_size, -kphase.imag(), hk_imag_pointer, ld_hk_2, hr_tmp, 1); |
| 702 | + hk_tmp ++; |
| 703 | + hr_tmp += this->col_size; |
| 704 | + } |
| 705 | + } |
| 706 | +} |
| 707 | + |
| 708 | +// add_from_matrix with explicit R_index - thread-safe version |
| 709 | +template <typename T> |
| 710 | +void AtomPair<T>::add_from_matrix(const int R_index, |
| 711 | + const T* hk, |
| 712 | + const int ld_hk, |
| 713 | + const T& kphase, |
| 714 | + const int hk_type) |
| 715 | +{ |
| 716 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 717 | + T* hr_tmp = matrix.get_pointer(); |
| 718 | + const T* hk_tmp = hk; |
| 719 | + // row major |
| 720 | + if (hk_type == 0) |
| 721 | + { |
| 722 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 723 | + for (int mu = 0; mu < this->row_size; mu++) |
| 724 | + { |
| 725 | + BlasConnector::axpy(this->col_size, kphase, hk_tmp, 1, hr_tmp, 1); |
| 726 | + hk_tmp += ld_hk; |
| 727 | + hr_tmp += this->col_size; |
| 728 | + } |
| 729 | + } |
| 730 | + // column major |
| 731 | + else if (hk_type == 1) |
| 732 | + { |
| 733 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 734 | + for (int mu = 0; mu < this->row_size; mu++) |
| 735 | + { |
| 736 | + BlasConnector::axpy(this->col_size, kphase, hk_tmp, ld_hk, hr_tmp, 1); |
| 737 | + ++hk_tmp; |
| 738 | + hr_tmp += this->col_size; |
| 739 | + } |
| 740 | + } |
| 741 | +} |
| 742 | + |
| 743 | +// add_to_matrix with explicit R_index - thread-safe version |
| 744 | +template <typename T> |
| 745 | +void AtomPair<T>::add_to_matrix(const int R_index, |
| 746 | + std::complex<T>* hk, |
| 747 | + const int ld_hk, |
| 748 | + const std::complex<T>& kphase, |
| 749 | + const int hk_type) const |
| 750 | +{ |
| 751 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 752 | + T* hr_tmp = matrix.get_pointer(); |
| 753 | + std::complex<T>* hk_tmp = hk; |
| 754 | + T* hk_real_pointer = nullptr; |
| 755 | + T* hk_imag_pointer = nullptr; |
| 756 | + const int ld_hk_2 = ld_hk * 2; |
| 757 | + // row major |
| 758 | + if (hk_type == 0) |
| 759 | + { |
| 760 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 761 | + for (int mu = 0; mu < this->row_size; mu++) |
| 762 | + { |
| 763 | + hk_real_pointer = (T*)hk_tmp; |
| 764 | + hk_imag_pointer = hk_real_pointer+1; |
| 765 | + BlasConnector::axpy(this->col_size, kphase.real(), hr_tmp, 1, hk_real_pointer, 2); |
| 766 | + BlasConnector::axpy(this->col_size, kphase.imag(), hr_tmp, 1, hk_imag_pointer, 2); |
| 767 | + hk_tmp += ld_hk; |
| 768 | + hr_tmp += this->col_size; |
| 769 | + } |
| 770 | + } |
| 771 | + // column major |
| 772 | + else if (hk_type == 1) |
| 773 | + { |
| 774 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 775 | + for (int mu = 0; mu < this->row_size; mu++) |
| 776 | + { |
| 777 | + hk_real_pointer = (T*)hk_tmp; |
| 778 | + hk_imag_pointer = hk_real_pointer+1; |
| 779 | + BlasConnector::axpy(this->col_size, kphase.real(), hr_tmp, 1, hk_real_pointer, ld_hk_2); |
| 780 | + BlasConnector::axpy(this->col_size, kphase.imag(), hr_tmp, 1, hk_imag_pointer, ld_hk_2); |
| 781 | + hk_tmp ++; |
| 782 | + hr_tmp += this->col_size; |
| 783 | + } |
| 784 | + } |
| 785 | +} |
| 786 | + |
| 787 | +// add_to_matrix with explicit R_index - thread-safe version |
| 788 | +template <typename T> |
| 789 | +void AtomPair<T>::add_to_matrix(const int R_index, |
| 790 | + T* hk, |
| 791 | + const int ld_hk, |
| 792 | + const T& kphase, |
| 793 | + const int hk_type) const |
| 794 | +{ |
| 795 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 796 | + T* hr_tmp = matrix.get_pointer(); |
| 797 | + T* hk_tmp = hk; |
| 798 | + // row major |
| 799 | + if (hk_type == 0) |
| 800 | + { |
| 801 | + hk_tmp += this->row_ap * ld_hk + this->col_ap; |
| 802 | + for (int mu = 0; mu < this->row_size; mu++) |
| 803 | + { |
| 804 | + BlasConnector::axpy(this->col_size, kphase, hr_tmp, 1, hk_tmp, 1); |
| 805 | + hk_tmp += ld_hk; |
| 806 | + hr_tmp += this->col_size; |
| 807 | + } |
| 808 | + } |
| 809 | + // column major |
| 810 | + else if (hk_type == 1) |
| 811 | + { |
| 812 | + hk_tmp += this->col_ap * ld_hk + this->row_ap; |
| 813 | + for (int mu = 0; mu < this->row_size; mu++) |
| 814 | + { |
| 815 | + BlasConnector::axpy(this->col_size, kphase, hr_tmp, 1, hk_tmp, ld_hk); |
| 816 | + ++hk_tmp; |
| 817 | + hr_tmp += this->col_size; |
| 818 | + } |
| 819 | + } |
| 820 | +} |
| 821 | + |
| 822 | +// add_to_array with explicit R_index - thread-safe version |
| 823 | +template <typename T> |
| 824 | +void AtomPair<T>::add_to_array(const int R_index, T* array, const T& kphase) const |
| 825 | +{ |
| 826 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 827 | + for (int i = 0; i < this->row_size * this->col_size; i++) |
| 828 | + { |
| 829 | + array[i] += matrix.get_pointer()[i] * kphase; |
| 830 | + } |
| 831 | +} |
| 832 | + |
| 833 | +// add_to_array with explicit R_index - thread-safe version |
| 834 | +template <typename T> |
| 835 | +void AtomPair<T>::add_to_array(const int R_index, std::complex<T>* array, const std::complex<T>& kphase) const |
| 836 | +{ |
| 837 | + const BaseMatrix<T>& matrix = values[R_index]; |
| 838 | + for (int i = 0; i < this->row_size * this->col_size; i++) |
| 839 | + { |
| 840 | + array[i] += matrix.get_pointer()[i] * kphase; |
| 841 | + } |
| 842 | +} |
| 843 | + |
588 | 844 | template <typename T> |
589 | 845 | std::tuple<std::vector<int>, T*> AtomPair<T>::get_matrix_values(int ir) const |
590 | 846 | { |
@@ -657,10 +913,6 @@ size_t AtomPair<T>::get_memory_size() const |
657 | 913 | } |
658 | 914 |
|
659 | 915 | // T of AtomPair can be double, float, or std::complex<double> |
660 | | -// explicit class-level instantiation for the remaining methods. |
661 | | -// Note: add_from_matrix, add_to_matrix, add_to_array methods are |
662 | | -// separately instantiated in atom_pair_kernels.cpp and atom_pair_output.cpp |
663 | | -// to keep their definitions there (out-of-line). |
664 | 916 | template class AtomPair<float>; |
665 | 917 | template class AtomPair<double>; |
666 | 918 | template class AtomPair<std::complex<double>>; |
|
0 commit comments