51 S, row_ind.size(), &*row_ind.begin(), col_ind.size(), &*col_ind.begin(),
52 &*
m.data().begin(), iora
72 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
79 boost::shared_ptr<std::vector<unsigned char>>
markerPtr;
87 template <
typename OP_SCHUR_ASSEMBLE_BASE>
105 std::vector<std::string> fields_name,
106 std::vector<boost::shared_ptr<Range>> field_ents,
110 bool sym_schur,
bool symm_op);
113 template <
typename I>
161 struct SchurElemMats :
public boost::enable_shared_from_this<SchurElemMats> {
187 template <
typename OP_SCHUR_ASSEMBLE_BASE>
202 member<SchurElemMats, const UId, &SchurElemMats::uidRow>,
203 member<SchurElemMats, const UId, &SchurElemMats::uidCol>
212 static boost::ptr_vector<MatrixDouble>
locMats;
279 const_mem_fun<Indexes, UId, &Indexes::getRowUId>,
280 const_mem_fun<Indexes, UId, &Indexes::getColUId>>
286 const_mem_fun<Indexes, UId, &Indexes::getFEUId>
292 const_mem_fun<Indexes, UId, &Indexes::getRowUId>
298 const_mem_fun<Indexes, UId, &Indexes::getColUId>
344 : iDX(idx), uidRow(uid_row), uidCol(uid_col) {}
358 auto get_row_indices = [&]() ->
const VectorInt & {
360 if (
auto stored_data_ptr =
362 return stored_data_ptr->entityIndices;
368 const auto &row_ind = get_row_indices();
371 const auto nb_rows = row_ind.size();
372 const auto nb_cols = col_ind.size();
375 if (mat.size1() != nb_rows) {
377 "Wrong mat size %d != %d", mat.size1(), nb_rows);
379 if (mat.size2() != nb_cols) {
381 "Wrong mat size %d != %d", mat.size2(), nb_cols);
386 auto get_uid = [](
auto &data) {
387 if (data.getFieldEntities().size() == 1) {
389 return data.getFieldEntities()[0]->getLocalUniqueId();
398 auto &uid0 = data.getFieldEntities()[0]->getLocalUniqueId();
404 for (
auto i = 1;
i < data.getFieldEntities().size(); ++
i) {
411 data.getFieldEntities()[
i]->getLocalUniqueId())
425 auto uid_row = get_uid(row_data);
426 auto uid_col = get_uid(col_data);
430 .find(boost::make_tuple(uid_row, uid_col));
462 auto asmb = [&](
auto &sm) {
463 sm.resize(nb_rows, nb_cols,
false);
467 asmb((*p.first)->getMat());
469 auto add_indices = [](
auto &storage,
auto &
ind) {
470 storage.resize(
ind.size(),
false);
471 noalias(storage) =
ind;
474 add_indices((*p.first)->getRowInd(), row_ind);
475 add_indices((*p.first)->getColInd(), col_ind);
480 auto asmb = [&](
auto &sm) {
484 if (sm.size1() != nb_rows) {
486 "Wrong mat or storage size %d != %d", sm.size1(), nb_rows);
488 if (sm.size2() != nb_cols) {
490 "Wrong mat or storage size %d != %d", sm.size2(), nb_cols);
503 "Assembly type not implemented");
508 CHKERR asmb((*it)->getMat());
528 MOFEM_LOG(
"SELF", Sev::noisy) <<
"Schur assemble begin";
537 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
double diag_val)
547 auto prb_ptr = fe_ptr->problemPtr;
550 auto find_loc_dof_index = [&](
auto glob) {
551 auto it = dof.find(glob);
552 if (it != dof.end()) {
553 return (*it)->getPetscLocalDofIdx();
555 MOFEM_LOG(
"SELF", Sev::error) <<
"Global dof " << glob;
558 "Global dof index not found in local numered dofs storage");
564 auto zero_ind = [&](
auto &
ind) {
565 for (
auto i = 0;
i <
ind.size(); ++
i) {
574 zero_ind(s->getRowInd());
575 zero_ind(s->getColInd());
582 template <
typename OP_SCHUR_ASSEMBLE_BASE>
584 std::vector<std::string> fields_name,
585 std::vector<boost::shared_ptr<Range>> field_ents,
588 :
OP(
NOSPACE,
OP::OPSPACE, symm_op), fieldsName(fields_name),
589 fieldEnts(field_ents), schurAO(schur_ao), schurMat(schur_mat),
590 symSchur(sym_schur) {}
592 template <
typename OP_SCHUR_ASSEMBLE_BASE>
593 template <
typename I>
600 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
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;
614 a00_uids.reserve(fieldsName.size());
615 for (
auto ss = 0; ss != fieldsName.size(); ++ss) {
616 auto field_bit = get_field_bit(fieldsName[ss]);
617 auto row_ents = fieldEnts[ss];
619 for (
auto p = row_ents->pair_begin(); p != row_ents->pair_end(); ++p) {
621 FieldEntity::getLoLocalEntityBitNumber(field_bit, p->first);
623 FieldEntity::getHiLocalEntityBitNumber(field_bit, p->second);
624 a00_uids.push_back(std::make_pair(lo_uid, hi_uid));
627 auto lo_uid = FieldEntity::getLoLocalEntityBitNumber(
628 field_bit, get_id_for_min_type<MBVERTEX>());
629 auto hi_uid = FieldEntity::getHiLocalEntityBitNumber(
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) {
640 FieldEntity::getFieldBitNumberFromUniqueId(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;
660 auto &storage = SchurElemMats::schurL2Storage;
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) {
768 CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
769 &*s->getRowInd().begin());
770 CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
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) {
833 for (
auto j = 0;
j != blockMat.size2(); ++
j) {
837 for (
auto j : zero_cols) {
838 for (
auto i = 0;
i != blockMat.size1(); ++
i) {
842 for (
auto i : zero_rows) {
846 CHKERR I::invertMat(blockMat, invMat);
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()) {
859 CHKERR AOApplicationToPetsc(schurAO, s->getRowInd().size(),
860 &*s->getRowInd().begin());
862 if (block_indexing.find(s->uidCol) == block_indexing.end()) {
863 CHKERR AOApplicationToPetsc(schurAO, s->getColInd().size(),
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(
916 boost::make_tuple(cm->first, FieldEntity::getHiBitNumberUId(
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";
995 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_opSchurAssembleEnd, 0, 0, 0, 0);
1006 const auto nb =
m.size1();
1008 if (nb !=
m.size2()) {
1010 "It should be square matrix %d != %d", nb,
m.size2());
1014 inv.resize(nb, nb,
false);
1016 m.resize(2 * nb, 2 * nb,
false);
1023 info =
lapack_dgetrf(nb, nb, &*inv.data().begin(), nb, &*ipiv.begin());
1026 "lapack error info = %d", info);
1027 info =
lapack_dgetri(nb, &*inv.data().begin(), nb, &*ipiv.begin(),
1028 &*
m.data().begin(),
m.data().size());
1031 "lapack error info = %d", info);
1042 const auto nb =
m.size1();
1044 if (nb !=
m.size2()) {
1046 "It should be square matrix %d != %d", nb,
m.size2());
1050 inv.resize(nb, nb,
false);
1056 info =
lapack_dgetrf(nb, nb, &*inv.data().begin(), nb, &*ipiv.begin());
1059 "lapack error info = %d", info);
1060 info =
lapack_dgetri(nb, &*inv.data().begin(), nb, &*ipiv.begin(),
1061 &*
m.data().begin(),
m.data().size());
1064 "lapack error info = %d", info);
1073 return doWorkImpl<SchurDSYSV>(side,
type, data);
1079 return doWorkImpl<SchurDGESV>(side,
type, data);
1089 auto cmp_uid_lo = [](
const boost::weak_ptr<FieldEntity> &
a,
const UId &b) {
1090 if (
auto a_ptr =
a.lock()) {
1091 if (a_ptr->getLocalUniqueId() < b)
1100 auto cmp_uid_hi = [](
const UId &b,
const boost::weak_ptr<FieldEntity> &
a) {
1101 if (
auto a_ptr =
a.lock()) {
1102 if (b < a_ptr->getLocalUniqueId())
1112 auto get_uid_pair = [](
const auto &field_id) {
1113 auto lo_uid = FieldEntity::getLocalUniqueIdCalculate(
1114 field_id, get_id_for_min_type<MBVERTEX>());
1115 auto hi_uid = FieldEntity::getLocalUniqueIdCalculate(
1116 field_id, get_id_for_max_type<MBENTITYSET>());
1117 return std::make_pair(lo_uid, hi_uid);
1122 auto lo = std::lower_bound(field_ents.begin(), field_ents.end(),
1124 auto hi = std::upper_bound(field_ents.begin(), field_ents.end(),
1126 return std::make_pair(lo, hi);
1131 auto row_extractor = [](
auto &e) {
return e->entityCacheRowDofs; };
1132 auto col_extractor = [](
auto &e) {
return e->entityCacheColDofs; };
1134 auto extract_data = [](
auto &&its,
auto extractor) {
1135 std::vector<std::tuple<UId, int, int, int>> data;
1136 data.reserve(std::distance(its.first, its.second));
1138 for (; its.first != its.second; ++its.first) {
1139 if (
auto e = its.first->lock()) {
1140 if (
auto cache = extractor(e).lock()) {
1142 for (
auto dof = cache->loHi[0]; dof != cache->loHi[1]; ++dof) {
1143 if ((*dof)->getPetscGlobalDofIdx() != -1)
1146 auto uid = e->getLocalUniqueId();
1149 auto glob = (*cache->loHi[0])->getPetscGlobalDofIdx();
1150 auto loc = (*cache->loHi[0])->getPetscLocalDofIdx();
1151 while (glob == -1 && cache->loHi[0] != cache->loHi[1]) {
1153 glob = (*cache->loHi[0])->getPetscGlobalDofIdx();
1154 loc = (*cache->loHi[0])->getPetscLocalDofIdx();
1156 data.emplace_back(uid, glob, nb_dofs, loc);
1160 for (
auto lo = cache->loHi[0]; lo != cache->loHi[1]; ++lo) {
1161 auto glob = (*lo)->getPetscGlobalDofIdx();
1164 "Wrong global index");
1170 data.emplace_back(uid, -1, 0, -1);
1178 auto data_ptr = boost::make_shared<BlockStructure>();
1182 auto fe_method = boost::shared_ptr<MoFEM::FEMethod>(
new MoFEM::FEMethod());
1184 for (
auto &
d : schur_fe_op_vec) {
1187 auto get_bit_numbers = [&
d](
auto op) {
1188 std::vector<FieldBitNumber> bit_numbers(
d.second.size());
1194 auto row_bit_numbers = get_bit_numbers(
1195 [&](
auto &p) {
return m_field_ptr->get_field_bit_number(p.first); });
1197 auto col_bit_numbers = get_bit_numbers(
1198 [&](
auto &p) {
return m_field_ptr->get_field_bit_number(p.second); });
1200 fe_method->preProcessHook = []() {
return 0; };
1201 fe_method->postProcessHook = []() {
return 0; };
1202 fe_method->operatorHook = [&]() {
1205 auto fe_uid = fe_method->numeredEntFiniteElementPtr->getLocalUniqueId();
1207 for (
auto f = 0;
f != row_bit_numbers.size(); ++
f) {
1210 extract_data(get_it_pair(fe_method->getRowFieldEnts(),
1211 get_uid_pair(row_bit_numbers[
f])),
1214 extract_data(get_it_pair(fe_method->getColFieldEnts(),
1215 get_uid_pair(col_bit_numbers[
f])),
1218 for (
auto &
r : row_data) {
1219 auto [r_uid, r_glob, r_nb_dofs, r_loc] =
r;
1220 for (
auto &
c : col_data) {
1221 auto [c_uid, c_glob, c_nb_dofs, c_loc] =
c;
1223 r_uid, c_uid, fe_uid, r_glob, c_glob, r_nb_dofs, c_nb_dofs,
1233 m_field_ptr->get_comm_size());
1238 for (
auto &
v : data_ptr->blockIndex.get<0>()) {
1239 v.getMatShift() = mem_size;
1240 mem_size +=
v.getNbCols() *
v.getNbRows();
1243 std::vector<double> tmp;
1244 if (mem_size > tmp.max_size())
1247 data_ptr->dataBlocksPtr =
1248 boost::make_shared<std::vector<double>>(mem_size, 0);
1252 CHKERR VecSetDM(ghost_x, PETSC_NULL);
1253 CHKERR VecSetDM(ghost_y, PETSC_NULL);
1255 data_ptr->ghostX = ghost_x;
1256 data_ptr->ghostY = ghost_y;
1262 Mat mat,
Vec x,
Vec y, InsertMode iora,
1264 int(DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator)>
1266 boost::shared_ptr<std::vector<double>> data_blocks_ptr,
1267 bool multiply_by_preconditioner);
1274 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1276 mat, x, y, INSERT_VALUES,
1278 [](DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator it) {
1279 return it->getMatShift();
1282 ctx->dataBlocksPtr,
true);
1286 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1288 mat, x, y, ADD_VALUES,
1290 [](DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator it) {
1291 return it->getMatShift();
1294 ctx->dataBlocksPtr,
true);
1304 const PetscInt rows[], PetscScalar diag,
1309 CHKERR MatShellGetContext(
A, (
void **)&ctx);
1311 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_zeroRowsAndCols, 0, 0, 0, 0);
1313 using BlockIndexView = multi_index_container<
1322 &DiagBlockIndex::Indexes::rowShift>
1329 &DiagBlockIndex::Indexes::colShift>
1335 BlockIndexView view;
1336 auto hint = view.get<0>().end();
1337 for (
auto &
v : ctx->blockIndex) {
1338 hint = view.insert(hint, &
v);
1341 const int *ptr = &rows[0];
1346 CHKERR PetscObjectGetComm((PetscObject)
A, &comm);
1348 MPI_Comm_size(comm, &size);
1354 CHKERR ISGetSize(is_local, &is_nb_rows);
1357 CHKERR ISView(is_local, PETSC_VIEWER_STDOUT_WORLD);
1360 CHKERR ISGetIndices(is_local, &ptr);
1364 CHKERR MatGetLocalSize(
A, &loc_m, &loc_n);
1366 for (
auto n = 0;
n != is_nb_rows; ++
n) {
1368 auto rlo = view.get<0>().lower_bound(row);
1369 auto rhi = view.get<0>().end();
1370 for (; rlo != rhi; ++rlo) {
1371 auto r_shift = row - (*rlo)->getRow();
1372 if (r_shift >= 0 && r_shift < (*rlo)->getNbRows()) {
1373 auto *ptr = &(*ctx->dataBlocksPtr)[(*rlo)->getMatShift()];
1374 for (
auto i = 0;
i != (*rlo)->getNbCols(); ++
i) {
1375 ptr[
i + r_shift * (*rlo)->getNbCols()] = 0;
1377 }
else if ((*rlo)->getRow() + (*rlo)->getNbRows() > row) {
1383 for (
auto n = 0;
n != is_nb_rows; ++
n) {
1385 auto clo = view.get<1>().lower_bound(col);
1386 auto chi = view.get<1>().end();
1387 for (; clo != chi; ++clo) {
1388 auto c_shift = col - (*clo)->getCol();
1389 if (c_shift >= 0 && c_shift < (*clo)->getNbCols()) {
1391 auto *ptr = &(*ctx->dataBlocksPtr)[(*clo)->getMatShift()];
1392 for (
auto i = 0;
i != (*clo)->getNbRows(); ++
i) {
1393 ptr[c_shift +
i * (*clo)->getNbCols()] = 0;
1399 (*clo)->getRow() == (*clo)->getCol() &&
1400 (*clo)->getLocRow() < loc_m && (*clo)->getLocCol() < loc_n
1403 auto r_shift = col - (*clo)->getCol();
1404 if (r_shift >= 0 && r_shift < (*clo)->getNbRows()) {
1405 ptr[c_shift + r_shift * (*clo)->getNbCols()] = diag;
1408 }
else if ((*clo)->getCol() + (*clo)->getNbCols() > col) {
1415 CHKERR ISRestoreIndices(is_local, &ptr);
1418 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_zeroRowsAndCols, 0, 0, 0, 0);
1426 CHKERR MatShellGetContext(
m, (
void **)&ctx);
1427 if (ctx->dataBlocksPtr)
1428 std::fill(ctx->dataBlocksPtr->begin(), ctx->dataBlocksPtr->end(), 0.);
1429 if (ctx->preconditionerBlocksPtr)
1430 std::fill(ctx->preconditionerBlocksPtr->begin(),
1431 ctx->preconditionerBlocksPtr->end(), 0.);
1437 CHKERR MatShellSetManageScalingShifts(mat_raw);
1438 CHKERR MatShellSetOperation(mat_raw, MATOP_MULT, (
void (*)(
void))
mult);
1439 CHKERR MatShellSetOperation(mat_raw, MATOP_MULT_ADD,
1441 CHKERR MatShellSetOperation(mat_raw, MATOP_SOLVE, (
void (*)(
void))
solve);
1442 CHKERR MatShellSetOperation(mat_raw, MATOP_SOLVE_ADD,
1444 CHKERR MatShellSetOperation(mat_raw, MATOP_ZERO_ENTRIES,
1446 CHKERR MatShellSetOperation(mat_raw, MATOP_ZERO_ROWS_COLUMNS,
1453 boost::shared_ptr<BlockStructure> data) {
1456 auto nb_local = problem_ptr->nbLocDofsRow;
1457 auto nb_global = problem_ptr->nbDofsRow;
1460 if (nb_local != problem_ptr->nbLocDofsCol) {
1462 <<
"Wrong size " << nb_local <<
" != " << problem_ptr->nbLocDofsCol;
1464 "nb. cols is inconsistent with nb. rows");
1466 if (nb_global != problem_ptr->nbDofsCol) {
1468 <<
"Wrong size " << nb_global <<
" != " << problem_ptr->nbDofsCol;
1470 "nb. cols is inconsistent with nb. rows");
1475 CHKERR PetscObjectGetComm((PetscObject)dm, &comm);
1478 CHKERR MatCreateShell(comm, nb_local, nb_local, nb_global, nb_global,
1479 (
void *)data.get(), &mat_raw);
1489 Mat mat,
Vec x,
Vec y, InsertMode iora,
1492 int(DiagBlockIndex::BlockIndex::nth_index<0>::type::iterator)>
1495 boost::shared_ptr<std::vector<double>> data_blocks_ptr,
1496 bool multiply_by_preconditioner) {
1499 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1501 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_BlockStructureMult, 0, 0, 0, 0);
1504 CHKERR VecGetLocalSize(x, &x_loc_size);
1506 CHKERR VecGetLocalSize(y, &y_loc_size);
1508 Vec ghost_x = ctx->ghostX;
1509 Vec ghost_y = ctx->ghostY;
1511 CHKERR VecCopy(x, ghost_x);
1515 CHKERR VecGhostGetLocalForm(ghost_y, &loc_ghost_y);
1517 CHKERR VecGetLocalSize(loc_ghost_y, &nb_y);
1518 CHKERR VecGetArray(loc_ghost_y, &y_array);
1519 for (
auto i = 0;
i != nb_y; ++
i)
1521 CHKERR VecRestoreArray(loc_ghost_y, &y_array);
1522 CHKERR VecGhostRestoreLocalForm(ghost_y, &loc_ghost_y);
1524 auto mult = [&](
int low_x,
int hi_x,
int low_y,
int hi_y) {
1529 CHKERR VecGhostGetLocalForm(ghost_x, &loc_ghost_x);
1530 CHKERR VecGetArray(loc_ghost_x, &x_array);
1534 CHKERR VecGhostGetLocalForm(ghost_y, &loc_ghost_y);
1536 CHKERR VecGetLocalSize(loc_ghost_y, &nb_y);
1537 CHKERR VecGetArray(loc_ghost_y, &y_array);
1539 double *block_ptr = &*data_blocks_ptr->begin();
1540 auto it = ctx->blockIndex.get<0>().begin();
1541 auto hi = ctx->blockIndex.get<0>().end();
1544 if (it->getLocRow() < low_y || it->getLocRow() >= hi_y ||
1545 it->getLocCol() < low_x || it->getLocCol() >= hi_x) {
1550 auto nb_rows = it->getNbRows();
1551 auto nb_cols = it->getNbCols();
1552 auto x_ptr = &x_array[it->getLocCol()];
1553 auto y_ptr = &y_array[it->getLocRow()];
1554 auto ptr = &block_ptr[shift_extractor(it)];
1557 cblas_dgemv(CblasRowMajor, CblasNoTrans, nb_rows, nb_cols, 1.0, ptr,
1558 nb_cols, x_ptr, 1, 1.0, y_ptr, 1);
1560 for (
auto r = 0;
r != nb_rows; ++
r) {
1561 for (
auto c = 0;
c != nb_cols; ++
c) {
1562 y_ptr[
r] += ptr[
r * nb_cols +
c] * x_ptr[
c];
1569 if (multiply_by_preconditioner && ctx->multiplyByPreconditioner) {
1571 if (!ctx->preconditionerBlocksPtr)
1573 "No parentBlockStructurePtr");
1575 auto preconditioner_ptr = &*ctx->preconditionerBlocksPtr->begin();
1577 auto it = ctx->blockIndex.get<0>().begin();
1578 auto hi = ctx->blockIndex.get<0>().end();
1581 if (it->getLocRow() < low_y || it->getLocRow() >= hi_y ||
1582 it->getLocCol() < low_x || it->getLocCol() >= hi_x) {
1587 if (it->getMatShift() != -1) {
1588 auto nb_rows = it->getNbRows();
1589 auto nb_cols = it->getNbCols();
1590 auto x_ptr = &x_array[it->getLocCol()];
1591 auto y_ptr = &y_array[it->getLocRow()];
1592 auto ptr = &preconditioner_ptr[it->getMatShift()];
1594 cblas_dgemv(CblasRowMajor, CblasNoTrans, nb_rows, nb_cols, 1.0, ptr,
1595 nb_cols, x_ptr, 1, 1.0, y_ptr, 1);
1597 for (
auto r = 0;
r != nb_rows; ++
r) {
1598 for (
auto c = 0;
c != nb_cols; ++
c) {
1599 y_ptr[
r] += ptr[
r * nb_cols +
c] * x_ptr[
c];
1609 CHKERR VecRestoreArray(loc_ghost_x, &x_array);
1610 CHKERR VecRestoreArray(loc_ghost_y, &y_array);
1611 CHKERR VecGhostRestoreLocalForm(ghost_x, &loc_ghost_x);
1612 CHKERR VecGhostRestoreLocalForm(ghost_y, &loc_ghost_y);
1616 constexpr
auto max_int = std::numeric_limits<int>::max();
1617 CHKERR VecGhostUpdateBegin(ghost_x, INSERT_VALUES, SCATTER_FORWARD);
1619 CHKERR VecGhostUpdateEnd(ghost_x, INSERT_VALUES, SCATTER_FORWARD);
1620 CHKERR mult(x_loc_size, max_int, 0, max_int);
1622 CHKERR VecGhostUpdateBegin(ghost_y, ADD_VALUES, SCATTER_REVERSE);
1623 CHKERR VecGhostUpdateEnd(ghost_y, ADD_VALUES, SCATTER_REVERSE);
1627 CHKERR VecCopy(ghost_y, y);
1630 CHKERR VecAXPY(y, 1., ghost_y);
1638 auto print_norm = [&](
auto msg,
auto y) {
1641 CHKERR VecNorm(y, NORM_2, &norm);
1643 CHKERR VecGetLocalSize(y, &nb_loc_y);
1645 CHKERR VecGetSize(y, &nb_y);
1647 << msg <<
" " << nb_y <<
" " << nb_loc_y <<
" norm " << norm;
1653 print_norm(
"mult_schur_block_shell insert x", x);
1654 print_norm(
"mult_schur_block_shell insert y", y);
1657 print_norm(
"mult_schur_block_shell add x", x);
1658 print_norm(
"mult_schur_block_shell add y", y);
1667 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_BlockStructureMult, 0, 0, 0, 0);
1674 using matrix_range = ublas::matrix_range<MatrixDouble>;
1675 using range = ublas::range;
1678 CHKERR MatShellGetContext(mat, (
void **)&ctx);
1680 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_BlockStructureSolve, 0, 0, 0, 0);
1682 if (ctx->multiplyByPreconditioner) {
1683 if (!ctx->preconditionerBlocksPtr)
1685 "No preconditionerBlocksPtr");
1688 if (iora == INSERT_VALUES) {
1689 CHKERR VecZeroEntries(x);
1693 CHKERR VecGetArray(x, &x_array);
1695 CHKERR VecGetArray(y, &y_array);
1697 std::vector<const DiagBlockIndex::Indexes *> blocks;
1702 auto &block_index = ctx->blockIndex.get<1>();
1703 for (
auto it = block_index.begin(); it != block_index.end();) {
1707 auto last_uid = it->getFEUId();
1708 while (it != block_index.end() && it->getFEUId() == last_uid) {
1709 blocks.push_back(&*it);
1714 std::map<UId, std::tuple<int, int, int>> block_indexing;
1715 for (
auto b : blocks) {
1716 if (block_indexing.find(b->getRowUId()) == block_indexing.end()) {
1717 block_indexing[b->getRowUId()] =
1718 std::make_tuple(b->getNbRows(), b->getNbRows(), b->getLocRow());
1721 std::sort(blocks.begin(), blocks.end(), [](
auto a,
auto b) {
1722 if (a->getRowUId() == b->getRowUId())
1723 return a->getColUId() < b->getColUId();
1725 return a->getRowUId() < b->getRowUId();
1729 int mat_block_size = 0;
1730 for (
auto &b : block_indexing) {
1731 auto &[index, size, loc] = b.second;
1732 index = mat_block_size;
1733 mat_block_size += size;
1736 std::vector<std::tuple<int, int, int, int, int>> block_data;
1737 block_data.reserve(blocks.size());
1738 for (
auto s : blocks) {
1739 auto ruid = s->getRowUId();
1740 auto cuid = s->getColUId();
1741 auto &rbi = block_indexing.at(ruid);
1742 auto &cbi = block_indexing.at(cuid);
1743 if (
auto shift = s->getMatShift(); shift != -1) {
1744 block_data.push_back(std::make_tuple(
1748 s->getNbRows(), s->getNbCols(),
1750 get<0>(rbi), get<0>(cbi))
1755 block_mat.resize(mat_block_size, mat_block_size,
false);
1757 for (
auto &bd : block_data) {
1758 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
1759 auto ptr = &(*(ctx->dataBlocksPtr))[shift];
1760 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
1762 cblas_dcopy(nb_cols, ptr, 1, sub_ptr, 1);
1765 if (ctx->multiplyByPreconditioner) {
1766 for (
auto &bd : block_data) {
1767 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
1768 auto ptr = &(*(ctx->preconditionerBlocksPtr))[shift];
1769 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
1771 cblas_daxpy(nb_cols, 1., ptr, 1, sub_ptr, 1);
1777 block_y.resize(mat_block_size,
false);
1780 for (
auto &b : block_indexing) {
1781 auto [index, size, loc] = b.second;
1782 auto ptr = &y_array[loc];
1783 cblas_dcopy(size, ptr, 1, &block_y[index], 1);
1787 constexpr
int nrhs = 1;
1788 ipiv.resize(mat_block_size,
false);
1790 mat_block_size, &*ipiv.data().begin(),
1791 &*block_y.data().begin(), mat_block_size);
1794 <<
"error lapack solve dgesv info = " << info;
1797 "error lapack solve dgesv info = %d", info);
1800 for (
auto &b : block_indexing) {
1801 auto [index, size, loc] = b.second;
1802 auto ptr = &x_array[loc];
1803 cblas_daxpy(size, 1.0, &block_y[index], 1, ptr, 1);
1807 CHKERR VecRestoreArray(x, &x_array);
1808 CHKERR VecRestoreArray(y, &y_array);
1811 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_BlockStructureSolve, 0, 0, 0, 0);
1818 std::vector<std::string> fields_name,
1819 std::vector<boost::shared_ptr<Range>> field_ents,
1821 using matrix_range = ublas::matrix_range<MatrixDouble>;
1822 using range = ublas::range;
1825 CHKERR MatShellGetContext(B, (
void **)&ctx);
1827 constexpr
bool debug =
false;
1829 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_AssembleSchurMat, 0, 0, 0, 0);
1835 if (ctx->multiplyByPreconditioner) {
1836 if (!ctx->preconditionerBlocksPtr)
1838 "No preconditionerBlocksPtr");
1841 std::vector<std::pair<UId, UId>> a00_fields_uids;
1842 a00_fields_uids.reserve(fields_name.size());
1843 auto it_field_ents = field_ents.begin();
1844 if (fields_name.size() != field_ents.size()) {
1846 "fields_name.size() != field_ents.size()");
1848 for (
auto &
f : fields_name) {
1850 if (*it_field_ents) {
1851 for (
auto p = (*it_field_ents)->pair_begin();
1852 p != (*it_field_ents)->pair_end(); ++p) {
1853 a00_fields_uids.emplace_back(
1854 DofEntity::getLoFieldEntityUId(bn, p->first),
1855 DofEntity::getHiFieldEntityUId(bn, p->second));
1858 a00_fields_uids.emplace_back(FieldEntity::getLoBitNumberUId(bn),
1859 FieldEntity::getHiBitNumberUId(bn));
1864 std::vector<const DiagBlockIndex::Indexes *> fe_blocks;
1865 std::vector<const DiagBlockIndex::Indexes *> a00_blocks;
1866 std::vector<const DiagBlockIndex::Indexes *> a01_blocks;
1867 std::vector<const DiagBlockIndex::Indexes *> a10_blocks;
1868 std::vector<const DiagBlockIndex::Indexes *> a11_blocks;
1874 std::vector<int> glob_idx;
1875 std::vector<int> row_glob_idx, col_glob_idx;
1877 auto &block_index = ctx->blockIndex.get<1>();
1878 for (
auto it = block_index.begin(); it != block_index.end();) {
1887 auto last_uid = it->getFEUId();
1888 while (it != block_index.end() && it->getFEUId() == last_uid) {
1889 fe_blocks.push_back(&*it);
1893 for (
auto b_ptr : fe_blocks) {
1894 auto check_id = [&](
auto uid) {
1895 for (
auto &uid_pair : a00_fields_uids) {
1896 if (uid >= uid_pair.first && uid < uid_pair.second) {
1902 auto r00 = check_id(b_ptr->getRowUId());
1903 auto c00 = check_id(b_ptr->getColUId());
1905 a00_blocks.push_back(b_ptr);
1906 }
else if (r00 && !c00) {
1907 a01_blocks.push_back(b_ptr);
1908 }
else if (!r00 && c00) {
1909 a10_blocks.push_back(b_ptr);
1911 a11_blocks.push_back(b_ptr);
1915 for (
auto b : a11_blocks) {
1917 row_glob_idx.resize(b->getNbRows());
1918 col_glob_idx.resize(b->getNbCols());
1920 std::iota(row_glob_idx.begin(), row_glob_idx.end(), b->getRow());
1921 std::iota(col_glob_idx.begin(), col_glob_idx.end(), b->getCol());
1922 CHKERR AOApplicationToPetsc(ao, col_glob_idx.size(), col_glob_idx.data());
1923 CHKERR AOApplicationToPetsc(ao, row_glob_idx.size(), row_glob_idx.data());
1925 auto shift = b->getMatShift();
1927 auto ptr = &(*(ctx->dataBlocksPtr))[shift];
1929 col_glob_idx.size(), col_glob_idx.data(), ptr,
1931 if (ctx->multiplyByPreconditioner) {
1932 auto ptr = &(*(ctx->preconditionerBlocksPtr))[shift];
1934 col_glob_idx.size(), col_glob_idx.data(), ptr,
1942 <<
"a00_blocks.size() " << a00_blocks.size();
1944 <<
"a01_blocks.size() " << a01_blocks.size();
1946 <<
"a10_blocks.size() " << a10_blocks.size();
1948 <<
"a11_blocks.size() " << a11_blocks.size();
1951 if(a00_blocks.size()) {
1953 for (
auto r : a00_blocks) {
1954 auto r_uid =
r->getRowUId();
1955 auto range = ctx->blockIndex.get<2>().equal_range(r_uid);
1956 for (
auto it = range.first; it != range.second; ++it) {
1957 if (it->getFEUId() != last_uid && it->getColUId() != r_uid) {
1958 a01_blocks.push_back(&*it);
1963 for (
auto c : a00_blocks) {
1964 auto c_uid =
c->getColUId();
1965 auto range = ctx->blockIndex.get<3>().equal_range(c_uid);
1966 for (
auto it = range.first; it != range.second; ++it) {
1967 if (it->getFEUId() != last_uid && it->getRowUId() != c_uid) {
1968 a10_blocks.push_back(&*it);
1973 auto sort = [](
auto &blocks) {
1974 std::sort(blocks.begin(), blocks.end(), [](
auto a,
auto b) {
1975 if (a->getRowUId() == b->getRowUId())
1976 return a->getColUId() < b->getColUId();
1978 return a->getRowUId() < b->getRowUId();
1988 <<
"a01_blocks.size() " << a01_blocks.size();
1990 <<
"a10_blocks.size() " << a10_blocks.size();
1994 auto set_local_indices = [](
auto &row_blocks,
auto &col_blocks) {
1996 std::map<UId, std::tuple<int, int, int>>
1998 for (
auto b : row_blocks) {
1999 if (block_indexing.find(b->getRowUId()) == block_indexing.end()) {
2000 block_indexing[b->getRowUId()] =
2001 std::make_tuple(b->getNbRows(), b->getNbRows(), b->getRow());
2004 for (
auto b : col_blocks) {
2005 if (block_indexing.find(b->getColUId()) == block_indexing.end()) {
2006 block_indexing[b->getColUId()] =
2007 std::make_tuple(b->getNbCols(), b->getNbCols(), b->getCol());
2011 int mat_block_size = 0;
2012 for (
auto &b : block_indexing) {
2014 auto &[index, size, glob] = b.second;
2015 index = mat_block_size;
2016 mat_block_size += size;
2018 return std::make_pair(block_indexing, mat_block_size);
2021 auto set_block = [&](
auto &blocks,
auto &row_indexing,
auto &col_indexing,
2023 int col_block_size) {
2024 std::vector<std::tuple<int, int, int, int, int>> block_data;
2025 block_data.reserve(blocks.size());
2026 for (
auto s : blocks) {
2027 auto ruid = s->getRowUId();
2028 auto cuid = s->getColUId();
2029 auto &rbi = row_indexing.at(ruid);
2030 auto &cbi = col_indexing.at(cuid);
2031 if (
auto shift = s->getMatShift(); shift != -1) {
2032 block_data.push_back(std::make_tuple(
2036 s->getNbRows(), s->getNbCols(),
2038 get<0>(rbi), get<0>(cbi))
2043 block_mat.resize(row_block_size, col_block_size,
false);
2045 for (
auto &bd : block_data) {
2046 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
2047 auto ptr = &(*(ctx->dataBlocksPtr))[shift];
2048 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
2050 cblas_dcopy(nb_cols, ptr, 1, sub_ptr, 1);
2053 if (ctx->multiplyByPreconditioner) {
2054 for (
auto &bd : block_data) {
2055 auto &[shift, nb_rows, nb_cols, ridx, cidx] = bd;
2056 auto ptr = &(*(ctx->preconditionerBlocksPtr))[shift];
2057 for (
auto i = 0;
i != nb_rows; ++
i, ptr += nb_cols) {
2059 cblas_daxpy(nb_cols, 1., ptr, 1, sub_ptr, 1);
2065 auto [a00_indexing, a00_size] = set_local_indices(a00_blocks, a00_blocks);
2066 auto [a11_indexing, a11_size] = set_local_indices(a10_blocks, a01_blocks);
2070 <<
"a00_indexing.size() " << a00_indexing.size() <<
" a00_size "
2073 <<
"a11_indexing.size() " << a11_indexing.size() <<
" a11_size "
2077 set_block(a00_blocks, a00_indexing, a00_indexing, block_mat_a00, a00_size,
2079 set_block(a01_blocks, a00_indexing, a11_indexing, block_mat_a01, a00_size,
2081 set_block(a10_blocks, a11_indexing, a00_indexing, block_mat_a10, a11_size,
2083 block_mat_a11.resize(a11_size, a11_size,
false);
2084 block_mat_a11.clear();
2086 block_mat_a00 = trans(block_mat_a00);
2087 block_mat_a01 = trans(block_mat_a01);
2089 ipiv.resize(block_mat_a00.size1(),
false);
2091 lapack_dgesv(block_mat_a00.size1(), block_mat_a01.size1(),
2092 block_mat_a00.data().data(), block_mat_a00.size2(),
2093 &*ipiv.data().begin(), block_mat_a01.data().data(),
2094 block_mat_a01.size2());
2097 "lapack error info = %d", info);
2100 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
2101 block_mat_a11.size1(), block_mat_a11.size2(),
2102 block_mat_a00.size1(), -1.,
2104 block_mat_a10.data().data(), block_mat_a10.size2(),
2106 block_mat_a01.data().data(), block_mat_a01.size2(),
2108 0., block_mat_a11.data().data(), block_mat_a11.size2()
2113 std::vector<int> glob_idx(block_mat_a11.size1());
2114 for (
auto &
r : a11_indexing) {
2115 auto [r_index, r_size, r_glob] =
r.second;
2116 std::iota(&glob_idx[idx], &glob_idx[idx + r_size], r_glob);
2119 CHKERR AOApplicationToPetsc(ao, glob_idx.size(), glob_idx.data());
2121 glob_idx.data(), block_mat_a11.data().data(),
2127 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_AssembleSchurMat, 0, 0, 0, 0);
2134 CHKERR MatShellGetContext(B, (
void **)&ctx);
2135 return ctx->dataBlocksPtr;
2138 boost::shared_ptr<std::vector<double>>
2141 CHKERR MatShellGetContext(B, (
void **)&ctx);
2142 return ctx->dataBlocksPtr;
2150 boost::shared_ptr<std::vector<double>> data_blocks_ptr) {
2161 PetscLogEventBegin(SchurEvents::MOFEM_EVENT_BlockStructureSetValues, 0, 0, 0,
2174 auto get_row_data = [&]() -> std::pair<bool, VectorInt> {
2176 if (
auto stored_data_ptr =
2178 MOFEM_LOG(
"SELF", Sev::warning) <<
"Can lead to unhandled behaviour";
2179 return std::make_pair(
true, stored_data_ptr->entityIndices);
2182 return std::make_pair(
false, row_data.
getIndices());
2185 auto row_indices = get_row_data();
2187 std::vector<int> ent_row_indices;
2188 std::vector<int> ent_col_indices;
2191 if (
auto r_cache = rent->entityCacheRowDofs.lock()) {
2193 auto &row_uid = rent->getLocalUniqueId();
2194 auto &row_ind = row_indices.second;
2197 if (
auto c_cache = cent->entityCacheColDofs.lock()) {
2199 auto &col_uid = cent->getLocalUniqueId();
2204 boost::make_tuple(row_uid, col_uid)
2213 <<
"missing block: row "
2217 "Block not allocated");
2222 auto ri = row_ind.begin();
2223 auto rie = row_ind.end();
2225 auto shift = shift_extractor(&*it);
2226 auto s_mat = &(*data_blocks_ptr)[shift];
2228 auto get_ent_indices = [](
auto &cache,
auto &
ind) {
2230 ind.reserve(std::distance(cache->loHi[0], cache->loHi[1]));
2231 for (
auto e = cache->loHi[0]; e != cache->loHi[1]; ++e) {
2232 auto glob = (*e)->getPetscGlobalDofIdx();
2234 ind.push_back(glob);
2238 get_ent_indices(r_cache, ent_row_indices);
2239 if (ent_row_indices.empty())
2241 get_ent_indices(c_cache, ent_col_indices);
2242 if (ent_col_indices.empty())
2245 if (mat.size1() == ent_row_indices.size() &&
2246 mat.size2() == ent_col_indices.size()) {
2248 if (iora == ADD_VALUES) {
2249 cblas_daxpy(mat.data().size(), 1.0, mat.data().data(), 1, s_mat,
2252 cblas_dcopy(mat.data().size(), mat.data().data(), 1, s_mat, 1);
2258 for (
auto re : ent_row_indices) {
2259 ri = std::find(ri, rie, re);
2260 if (!(ri == rie && *ri != -1)) {
2262 auto ci = col_ind.begin();
2263 auto cie = col_ind.end();
2266 for (
auto ce : ent_col_indices) {
2267 ci = std::find(ci, cie, ce);
2268 if (!(ci == cie && *ci != -1)) {
2269 auto &
m = s_mat[row * ent_col_indices.size() + col];
2270 if (iora == ADD_VALUES) {
2271 m += mat(std::distance(row_ind.begin(), ri),
2272 std::distance(col_ind.begin(), ci));
2274 m = mat(std::distance(row_ind.begin(), ri),
2275 std::distance(col_ind.begin(), ci));
2291 PetscLogEventEnd(SchurEvents::MOFEM_EVENT_BlockStructureSetValues, 0, 0, 0,
2303 PetscBool is_mat_shell = PETSC_FALSE;
2304 PetscObjectTypeCompare((PetscObject)
M, MATSHELL, &is_mat_shell);
2307 CHKERR MatShellGetContext(
M, (
void **)&ctx);
2309 ctx, row_data, col_data, mat, iora,
2311 ctx->dataBlocksPtr);
2313 CHKERR MatSetValues<AssemblyTypeSelector<PETSC>>(
M, row_data, col_data, mat,
2324 PetscBool is_mat_shell = PETSC_FALSE;
2325 PetscObjectTypeCompare((PetscObject)
M, MATSHELL, &is_mat_shell);
2328 CHKERR MatShellGetContext(
M, (
void **)&ctx);
2329 if (!ctx->preconditionerBlocksPtr)
2331 "No preconditionerBlocksPtr");
2333 ctx, row_data, col_data, mat, iora,
2335 ctx->preconditionerBlocksPtr);
2337 CHKERR MatSetValues<AssemblyTypeSelector<PETSC>>(
M, row_data, col_data, mat,
2346 boost::shared_ptr<BlockStructure> block_mat_data_ptr,
2348 std::vector<std::string> fields_names,
2349 std::vector<boost::shared_ptr<Range>>
2351 bool add_preconditioner_block) {
2353 if (!block_mat_data_ptr) {
2357 if (fields_names.size() != field_ents.size())
2359 "fields_names.size() != field_ents.size()");
2361 auto [schur_dm, block_dm] = dms;
2366 auto schur_dofs_row = schur_prb->getNumeredRowDofsPtr();
2367 auto schur_dofs_col = schur_prb->getNumeredColDofsPtr();
2368 auto block_dofs_row = block_prb->getNumeredRowDofsPtr();
2369 auto block_dofs_col = block_prb->getNumeredColDofsPtr();
2371 auto ao_schur_row = schur_prb->getSubData()->getSmartRowMap();
2372 auto ao_schur_col = schur_prb->getSubData()->getSmartColMap();
2373 auto ao_block_row = block_prb->getSubData()->getSmartRowMap();
2374 auto ao_block_col = block_prb->getSubData()->getSmartColMap();
2380 CHKERR VecSetDM(schur_vec_x, PETSC_NULL);
2381 CHKERR VecSetDM(block_vec_x, PETSC_NULL);
2382 CHKERR VecSetDM(schur_vec_y, PETSC_NULL);
2383 CHKERR VecSetDM(block_vec_y, PETSC_NULL);
2385 auto find_field_ent = [&](
auto uid,
auto prb,
auto rc) {
2386 boost::shared_ptr<NumeredDofEntity_multiIndex> dofs;
2390 dofs = prb->getNumeredRowDofsPtr();
2393 dofs = prb->getNumeredColDofsPtr();
2402 return boost::shared_ptr<NumeredDofEntity>();
2408 return boost::shared_ptr<NumeredDofEntity>();
2411 std::array<boost::shared_ptr<BlockStructure>, 4> data_ptrs;
2413 for (
auto r = 0;
r != 3; ++
r) {
2414 data_ptrs[
r] = boost::make_shared<BlockStructure>();
2415 data_ptrs[
r]->dataBlocksPtr = block_mat_data_ptr->dataBlocksPtr;
2417 data_ptrs[3] = boost::make_shared<BlockStructure>();
2418 data_ptrs[3]->dataBlocksPtr = block_mat_data_ptr->dataBlocksPtr;
2420 data_ptrs[0]->ghostX = schur_vec_x;
2421 data_ptrs[0]->ghostY = schur_vec_y;
2422 data_ptrs[1]->ghostX = block_vec_x;
2423 data_ptrs[1]->ghostY = schur_vec_y;
2424 data_ptrs[2]->ghostX = schur_vec_x;
2425 data_ptrs[2]->ghostY = block_vec_y;
2426 data_ptrs[3]->ghostX = block_vec_x;
2427 data_ptrs[3]->ghostY = block_vec_y;
2430 for (
auto &
d : block_mat_data_ptr->blockIndex.get<0>()) {
2432 auto insert = [&](
auto &
m,
auto &dof_r,
auto &dof_c,
auto &
d) {
2436 d.getRowUId(),
d.getColUId(),
d.getFEUId(),
2438 dof_r->getPetscGlobalDofIdx(), dof_c->getPetscGlobalDofIdx(),
2440 d.getNbRows(),
d.getNbCols(),
2442 dof_r->getPetscLocalDofIdx(), dof_c->getPetscLocalDofIdx(),
2449 auto dof_schur_row = find_field_ent(
d.getRowUId(), schur_prb,
ROW);
2450 auto dof_schur_col = find_field_ent(
d.getColUId(), schur_prb,
COL);
2451 auto dof_block_row = find_field_ent(
d.getRowUId(), block_prb,
ROW);
2452 auto dof_block_col = find_field_ent(
d.getColUId(), block_prb,
COL);
2454 if (dof_schur_row && dof_schur_col) {
2455 insert(data_ptrs[0]->blockIndex, dof_schur_row, dof_schur_col,
d);
2458 if (dof_schur_row && dof_block_col) {
2459 insert(data_ptrs[1]->blockIndex, dof_schur_row, dof_block_col,
d);
2462 if (dof_block_row && dof_schur_col) {
2463 insert(data_ptrs[2]->blockIndex, dof_block_row, dof_schur_col,
d);
2466 if (dof_block_row && dof_block_col) {
2467 insert(data_ptrs[3]->blockIndex, dof_block_row, dof_block_col,
d);
2474 auto set_up_a00_data = [&](
auto inv_block_data) {
2477 if (add_preconditioner_block) {
2478 auto preconditioned_block = boost::make_shared<std::vector<double>>(
2479 inv_block_data->dataBlocksPtr->size(), 0);
2480 inv_block_data->preconditionerBlocksPtr = preconditioned_block;
2481 inv_block_data->multiplyByPreconditioner =
true;
2482 block_mat_data_ptr->preconditionerBlocksPtr =
2483 inv_block_data->preconditionerBlocksPtr;
2484 block_mat_data_ptr->multiplyByPreconditioner =
false;
2490 CHKERR set_up_a00_data(data_ptrs[3]);
2493 CHKERR PetscObjectGetComm((PetscObject)schur_dm, &comm);
2495 auto create_shell_mat = [&](
auto nb_r_loc,
auto nb_c_loc,
auto nb_r_glob,
2496 auto nb_c_glob,
auto data_ptr) {
2498 CHKERR MatCreateShell(comm, nb_r_loc, nb_c_loc, nb_r_glob, nb_c_glob,
2499 (
void *)data_ptr.get(), &mat_raw);
2504 auto schur_nb_global = schur_prb->getNbDofsRow();
2505 auto block_nb_global = block_prb->getNbDofsRow();
2506 auto schur_nb_local = schur_prb->getNbLocalDofsRow();
2507 auto block_nb_local = block_prb->getNbLocalDofsRow();
2509 std::array<SmartPetscObj<Mat>, 4> mats_array;
2511 create_shell_mat(schur_nb_local, schur_nb_local, schur_nb_global,
2512 schur_nb_global, data_ptrs[0]);
2514 create_shell_mat(schur_nb_local, block_nb_local, schur_nb_global,
2515 block_nb_global, data_ptrs[1]);
2517 create_shell_mat(block_nb_local, schur_nb_local, block_nb_global,
2518 schur_nb_global, data_ptrs[2]);
2520 create_shell_mat(block_nb_local, block_nb_local, block_nb_global,
2521 block_nb_global, data_ptrs[3]);
2524 <<
"(0, 0) " << schur_nb_local <<
" " << schur_nb_global <<
" "
2525 << data_ptrs[0]->blockIndex.size();
2527 <<
"(0, 1) " << schur_nb_local <<
" " << block_nb_local <<
" "
2528 << schur_nb_global <<
" " << block_nb_global <<
" "
2529 << data_ptrs[1]->blockIndex.size();
2531 <<
"(1, 0) " << block_nb_local <<
" " << schur_nb_local <<
" "
2532 << block_nb_global <<
" " << schur_nb_global <<
" "
2533 << data_ptrs[2]->blockIndex.size();
2535 <<
"(1, 1) " << block_nb_local <<
" " << block_nb_global <<
" "
2536 << data_ptrs[3]->blockIndex.size();
2540 auto schur_is = schur_prb->getSubData()->getSmartRowIs();
2541 auto block_is = block_prb->getSubData()->getSmartRowIs();
2543 return boost::make_shared<NestSchurData>(
2545 mats_array, data_ptrs, block_mat_data_ptr,
2546 std::make_pair(schur_is, block_is)
2551 std::pair<SmartPetscObj<Mat>, boost::shared_ptr<NestSchurData>>
2554 if (!schur_net_data_ptr)
2557 auto [mat_arrays, data_ptrs, block_mat_data_ptr, iss] = *schur_net_data_ptr;
2558 auto [schur_is, block_is] = iss;
2560 std::array<IS, 2> is_a = {schur_is, block_is};
2561 std::array<Mat, 4> mats_a = {
2563 mat_arrays[0], mat_arrays[1], mat_arrays[2], mat_arrays[3]
2568 CHKERR PetscObjectGetComm((PetscObject)mat_arrays[0], &comm);
2573 comm, 2, is_a.data(), 2, is_a.data(), mats_a.data(), &mat_raw
2584 OpSchurAssembleBase *
2586 std::vector<boost::shared_ptr<Range>> field_ents,
2588 bool sym_schur,
bool symm_op) {
2591 schur, sym_schur, symm_op);
2594 schur, sym_schur, symm_op);
2598 boost::shared_ptr<std::vector<unsigned char>> marker_ptr,
double diag_val) {
2602 OpSchurAssembleBase *
2604 std::vector<boost::shared_ptr<Range>> field_ents,
2607 std::vector<bool> sym_schur,
bool symm_op,
2608 boost::shared_ptr<BlockStructure> diag_blocks) {
2610 fields_name, field_ents, sequence_of_aos.back(), sequence_of_mats.back(),
2611 sym_schur.back(), symm_op);
2614 OpSchurAssembleBase *
2616 std::vector<boost::shared_ptr<Range>> field_ents,
2619 std::vector<bool> sym_schur,
2620 std::vector<double> diag_eps,
bool symm_op,
2621 boost::shared_ptr<BlockStructure> diag_blocks) {
2623 fields_name, field_ents, sequence_of_aos.back(), sequence_of_mats.back(),
2624 sym_schur.back(), symm_op);
2630 auto apply = [](PC pc,
Vec f,
Vec x) {
2638 CHKERR PCSetType(pc, PCSHELL);
2639 CHKERR PCShellSetApply(pc, apply);
2640 CHKERR PCShellSetName(pc,
"MoFEMSchurBlockPC");
2663 return MatSetValues<AssemblyTypeSelector<PETSC>>(
M, row_data, col_data, mat,
2668 SchurBackendMatSetValuesPtr::MatSetValuesPtr
2678 CHKERR assembleStorage(row_data, col_data, mat, iora);
2679 CHKERR SchurBackendMatSetValuesPtr::matSetValuesPtr(
M, row_data, col_data,
2712 SchurBackendMatSetValuesPtr::matSetValuesBlockPtr =
2721 CHKERR assembleStorage(row_data, col_data, mat, iora);
2722 CHKERR SchurBackendMatSetValuesPtr::matSetValuesBlockPtr(
M, row_data,
2723 col_data, mat, iora);
2745 SchurBackendMatSetValuesPtr::matSetValuesPreconditionedBlockPtr =
2753 CHKERR assembleStorage(row_data, col_data, mat, iora);
2754 CHKERR SchurBackendMatSetValuesPtr::matSetValuesPreconditionedBlockPtr(
2755 M, row_data, col_data, mat, iora);
2771 if (block_mat_data->multiplyByPreconditioner) {
2772 block_mat_data->multiplyByPreconditioner =
false;
2774 if (!block_mat_data->preconditionerBlocksPtr) {
2776 "preconditionerBlocksPtr not set");
2778 block_mat_data->multiplyByPreconditioner =
true;
2785 std::string filename) {
2791 ReadUtilIface *iface;
2792 CHKERR moab.query_interface(iface);
2794 auto size = block_mat_data->blockIndex.size();
2797 vector<double *> arrays_coord;
2798 CHKERR iface->get_node_coords(3, 4 * size, 0, starting_vertex, arrays_coord);
2801 CHKERR iface->get_element_connect(size, 4, MBQUAD, 0, starting_handle, array);
2803 double def_val_nrm2 = 0;
2805 CHKERR moab.tag_get_handle(
"nrm2", 1, MB_TYPE_DOUBLE, th_nrm2,
2806 MB_TAG_CREAT | MB_TAG_DENSE, &def_val_nrm2);
2808 int def_val_mat_shift = 0;
2810 CHKERR moab.tag_get_handle(
"mat_shift", 1, MB_TYPE_INTEGER, th_mat_shift,
2811 MB_TAG_CREAT | MB_TAG_DENSE, &def_val_mat_shift);
2814 for (
auto &
d : block_mat_data->blockIndex) {
2815 auto row =
d.getRow();
2816 auto col =
d.getCol();
2817 auto nb_rows =
d.getNbRows();
2818 auto nb_cols =
d.getNbCols();
2819 auto mat_shift =
d.getMatShift();
2822 arrays_coord[0][4 *
i + 0] = row;
2823 arrays_coord[1][4 *
i + 0] = col;
2824 arrays_coord[2][4 *
i + 0] = 0;
2827 arrays_coord[0][4 *
i + 1] = row + nb_rows;
2828 arrays_coord[1][4 *
i + 1] = col;
2829 arrays_coord[2][4 *
i + 1] = 0;
2832 arrays_coord[0][4 *
i + 2] = row + nb_rows;
2833 arrays_coord[1][4 *
i + 2] = col + nb_cols;
2834 arrays_coord[2][4 *
i + 2] = 0;
2837 arrays_coord[0][4 *
i + 3] = row;
2838 arrays_coord[1][4 *
i + 3] = col + nb_cols;
2839 arrays_coord[2][4 *
i + 3] = 0;
2842 array[4 *
i + 0] = starting_vertex + 4 *
i + 0;
2843 array[4 *
i + 1] = starting_vertex + 4 *
i + 1;
2844 array[4 *
i + 2] = starting_vertex + 4 *
i + 2;
2845 array[4 *
i + 3] = starting_vertex + 4 *
i + 3;
2847 auto prt = &(*block_mat_data->dataBlocksPtr)[
d.getMatShift()];
2848 auto nrm2 = cblas_dnrm2(nb_rows * nb_cols, prt, 1);
2850 CHKERR moab.tag_set_data(th_nrm2, &ele, 1, &nrm2);
2851 CHKERR moab.tag_set_data(th_mat_shift, &ele, 1, &mat_shift);
2856 CHKERR iface->update_adjacencies(starting_handle, size, 4, array);
2857 CHKERR moab.write_file(filename.c_str(),
"VTK",
"");