605 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin -> end";
608 auto get_a00_uids = [&]() {
609 auto get_field_bit = [&](
auto &name) {
610 return OP::getPtrFE()->mField.get_field_bit_number(name);
613 std::vector<std::pair<UId, UId>> a00_uids;
615 for (
auto ss = 0; ss !=
fieldsName.size(); ++ss) {
616 auto field_bit = get_field_bit(
fieldsName[ss]);
619 for (
auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
624 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
628 field_bit, get_id_for_min_type<MBVERTEX>());
630 field_bit, get_id_for_max_type<MBENTITYSET>());
631 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
638 auto get_field_name = [&](
auto uid) {
642 auto list_storage = [&](
auto &storage) {
645 for (
auto &p : storage) {
647 <<
"List schur storage: " <<
i <<
" " << p->iDX <<
": "
648 << get_field_name(p->uidRow) <<
" " << get_field_name(p->uidCol)
649 <<
" : " << p->getMat().size1() <<
" " << p->getMat().size2();
656 auto assemble_dense_blocks = [&]() {
657 using matrix_range = ublas::matrix_range<MatrixDouble>;
658 using range = ublas::range;
662 auto assemble_schur = [
this](
auto &
m,
auto &uid_row,
auto &uid_col,
663 auto *row_ind_ptr,
auto *col_ind_ptr) {
668 if (
auto ierr = this->assembleSchurMat(
670 schurMat, uid_row, *row_ind_ptr, uid_col, *col_ind_ptr,
m,
675 auto field_ents = OP::getPtrFE()->mField.get_field_ents();
676 auto row_ent_it = field_ents->find(uid_row);
677 auto col_ent_it = field_ents->find(uid_col);
679 if (row_ent_it != field_ents->end())
681 <<
"Assemble row entity: " << (*row_ent_it)->getName() <<
" "
682 << (*col_ent_it)->getEntTypeName() <<
" side "
683 << (*row_ent_it)->getSideNumber();
684 if (col_ent_it != field_ents->end())
686 <<
"Assemble col entity: " << (*col_ent_it)->getName() <<
" "
687 << (*col_ent_it)->getEntTypeName() <<
" side "
688 << (*col_ent_it)->getSideNumber();
697 auto a00_uids = get_a00_uids();
699 auto get_block_indexing = [&](
auto &a00_uids) {
701 std::vector<const SchurElemMats *> block_list;
702 block_list.reserve(storage.size());
703 for (
auto &rp_uid : a00_uids) {
704 auto [rlo_uid, rhi_uid] = rp_uid;
705 for (
auto &cp_uid : a00_uids) {
706 auto [clo_uid, chi_uid] = cp_uid;
709 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
710 boost::make_tuple(rlo_uid, clo_uid));
712 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
713 boost::make_tuple(rhi_uid, chi_uid));
715 for (; it != hi_it; ++it) {
716 if ((*it)->uidRow >= rlo_uid && (*it)->uidRow < rhi_uid &&
717 (*it)->uidCol >= clo_uid && (*it)->uidCol < chi_uid) {
718 block_list.push_back(*it);
725 std::map<UId, std::pair<size_t, const VectorInt *>>
727 for (
auto d : block_list) {
728 if (block_indexing.find(
d->uidRow) == block_indexing.end()) {
729 block_indexing[
d->uidRow] =
730 std::make_pair(
d->getRowInd().size(), &(
d->getRowInd()));
735 int mat_block_size = 0;
736 for (
auto &p_uid : a00_uids) {
737 auto [lo_uid, hi_uid] = p_uid;
738 auto lo = block_indexing.lower_bound(lo_uid);
739 auto up = block_indexing.upper_bound(hi_uid);
740 for (; lo != up; ++lo) {
741 lo->second.first = mat_block_size;
742 mat_block_size += lo->second.second->size();
746 return std::make_tuple(block_list, block_indexing, mat_block_size);
749 auto get_schur_block_list = [&](
auto &block_indexing) {
750 std::vector<const SchurElemMats *> block_list;
751 block_list.reserve(storage.size());
752 for (
auto &s : storage) {
753 if (block_indexing.find(s->uidRow) == block_indexing.end() &&
754 block_indexing.find(s->uidCol) == block_indexing.end()) {
755 block_list.push_back(s);
761 auto bi = get_block_indexing(a00_uids);
762 auto &[block_list, block_indexing, block_mat_size] = bi;
764 if (block_mat_size == 0) {
767 for (
auto &s : storage) {
769 &*s->getRowInd().begin());
771 &*s->getColInd().begin());
775 for (
auto &s : storage) {
776 auto &
m = s->getMat();
777 CHKERR assemble_schur(
m, s->uidRow, s->uidCol, &(s->getRowInd()),
783 blockMat.resize(block_mat_size, block_mat_size,
false);
786 auto get_range = [](
auto &bi) {
787 return range(bi.first, bi.first + bi.second->size());
790 for (
auto &s : block_list) {
791 auto &
m = s->getMat();
793 if (block_indexing.find(s->uidRow) == block_indexing.end())
795 if (block_indexing.find(s->uidCol) == block_indexing.end())
799 auto &rbi = block_indexing.at(s->uidRow);
800 auto &cbi = block_indexing.at(s->uidCol);
802 auto sub_mat = matrix_range(
blockMat, get_range(rbi), get_range(cbi));
806 auto get_zeroed_indices = [&](
auto extractor_uid,
auto extractor_ind) {
807 auto &[block_list, block_indexing, block_mat_size] = bi;
808 std::vector<int> zeroed_indices;
809 zeroed_indices.reserve(block_mat_size);
810 for (
auto &s : block_list) {
811 auto &bi = block_indexing.at(extractor_uid(s));
812 auto &
ind = extractor_ind(s);
813 for (
auto it =
ind.begin(); it !=
ind.end(); ++it) {
815 auto idx = bi.first + std::distance(
ind.begin(), it);
816 zeroed_indices.push_back(idx);
820 std::sort(zeroed_indices.begin(), zeroed_indices.end());
821 auto it = std::unique(zeroed_indices.begin(), zeroed_indices.end());
822 zeroed_indices.resize(std::distance(zeroed_indices.begin(), it));
823 return zeroed_indices;
825 auto zero_rows = get_zeroed_indices(
826 [](
auto &s) {
return s->uidRow; },
827 [](
auto &s) ->
VectorInt & {
return s->getRowInd(); });
828 auto zero_cols = get_zeroed_indices(
829 [](
auto &s) {
return s->uidCol; },
830 [](
auto &s) ->
VectorInt & {
return s->getColInd(); });
832 for (
auto i : zero_rows) {
837 for (
auto j : zero_cols) {
842 for (
auto i : zero_rows) {
849 for (
auto &s : block_list) {
850 auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
851 boost::make_tuple(s->uidRow, s->uidCol));
852 storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
857 for (
auto &s : storage) {
858 if (block_indexing.find(s->uidRow) == block_indexing.end()) {
860 &*s->getRowInd().begin());
862 if (block_indexing.find(s->uidCol) == block_indexing.end()) {
864 &*s->getColInd().begin());
869 auto schur_block_list = get_schur_block_list(block_indexing);
870 for (
auto &s : schur_block_list) {
871 auto &
m = s->getMat();
872 CHKERR assemble_schur(
m, s->uidRow, s->uidCol, &(s->getRowInd()),
875 for (
auto &s : schur_block_list) {
876 auto it = storage.template get<SchurElemMats::uid_mi_tag>().find(
877 boost::make_tuple(s->uidRow, s->uidCol));
878 storage.template get<SchurElemMats::uid_mi_tag>().erase(it);
880 schur_block_list.clear();
884 auto rp_uid_it = a00_uids.begin(); rp_uid_it != a00_uids.end();
888 auto [rlo_uid, rhi_uid] = *rp_uid_it;
889 for (
auto rm = block_indexing.lower_bound(rlo_uid);
890 rm != block_indexing.upper_bound(rhi_uid); ++rm) {
891 auto &rbi = rm->second;
894 storage.template get<SchurElemMats::col_mi_tag>().lower_bound(
897 storage.template get<SchurElemMats::col_mi_tag>().upper_bound(
902 auto cp_uid_it = (
symSchur) ? rp_uid_it : a00_uids.begin();
903 cp_uid_it != a00_uids.end(); ++cp_uid_it
906 auto [clo_uid, chi_uid] = *cp_uid_it;
907 for (
auto cm = block_indexing.lower_bound(clo_uid);
908 cm != block_indexing.upper_bound(chi_uid); ++cm) {
909 auto &cbi = cm->second;
912 storage.template get<SchurElemMats::uid_mi_tag>().lower_bound(
913 boost::make_tuple(cm->first, 0));
915 storage.template get<SchurElemMats::uid_mi_tag>().upper_bound(
920 matrix_range(
invMat, get_range(rbi), get_range(cbi));
921 bM.resize(sub_inv_mat.size1(), sub_inv_mat.size2());
922 noalias(
bM) = sub_inv_mat;
924 for (
auto a_lo = a_lo_tmp; a_lo != a_hi; ++a_lo) {
926 if (block_indexing.find((*a_lo)->uidRow) != block_indexing.end())
928 "Wrong a_lo->uidRow");
931 auto &
a = (*a_lo)->getMat();
932 abM.resize(
a.size1(),
bM.size2(),
false);
934 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
936 a.size1(),
bM.size2(),
a.size2(), 1.,
938 &*
a.data().begin(),
a.size2(),
940 &*
bM.data().begin(),
bM.size2(), 0.,
942 &*
abM.data().begin(),
abM.size2());
944 for (
auto c_lo = c_lo_tmp; c_lo != c_hi; ++c_lo) {
946 if (block_indexing.find((*c_lo)->uidCol) !=
947 block_indexing.end())
949 "Wrong a_lo->uidRow");
952 auto &
c = (*c_lo)->getMat();
954 abcM.resize(
abM.size1(),
c.size2(),
false);
957 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
958 abM.size1(),
c.size2(),
abM.size2(), -1.,
960 &*
abM.data().begin(),
abM.size2(),
962 &*
c.data().begin(),
c.size2(), 0.,
964 &*
abcM.data().begin(),
abcM.size2());
966 CHKERR assemble_schur(
abcM, (*a_lo)->uidRow, (*c_lo)->uidCol,
967 &((*a_lo)->getRowInd()),
968 &((*c_lo)->getColInd()));
970 if (
symSchur && rp_uid_it != cp_uid_it) {
972 CHKERR assemble_schur(
abcM, (*c_lo)->uidCol, (*a_lo)->uidRow,
973 &((*c_lo)->getColInd()),
974 &((*a_lo)->getRowInd()));
987 CHKERR assemble_dense_blocks();
991 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble done";