49 #include <Epetra_Export.h> 50 #include <Epetra_Import.h> 51 #include <Epetra_Util.h> 52 #include <Epetra_Map.h> 53 #include <Epetra_Comm.h> 54 #include <Epetra_CrsMatrix.h> 55 #include <Epetra_Vector.h> 56 #include <Epetra_Directory.h> 57 #include <Epetra_HashTable.h> 58 #include <Epetra_Distributor.h> 60 #include <Teuchos_TimeMonitor.hpp> 75 template<
typename int_type>
76 double sparsedot(
double* u, int_type* u_ind,
int u_len,
77 double* v, int_type* v_ind,
int v_len)
84 while(v_idx < v_len && u_idx < u_len) {
85 int_type ui = u_ind[u_idx];
86 int_type vi = v_ind[v_idx];
95 result += u[u_idx++]*v[v_idx++];
104 template<
typename int_type>
113 for(i=0; i<Aview.
numRows; ++i) {
116 for(i=0; i<Bview.
numRows; ++i) {
126 int numBcols = Bview.
colMap->NumMyElements();
129 int iworklen = maxlen*2 + numBcols;
130 int_type* iwork =
new int_type[iworklen];
132 int_type * bcols = iwork+maxlen*2;
134 Bview.
colMap->MyGlobalElementsPtr(bgids);
135 double* bvals =
new double[maxlen*2];
136 double* avals = bvals+maxlen;
138 int_type max_all_b = (int_type) Bview.
colMap->MaxAllGID64();
139 int_type min_all_b = (int_type) Bview.
colMap->MinAllGID64();
143 for(i=0; i<numBcols; ++i) {
144 int blid = Bview.
colMap->LID(bgids[i]);
145 bcols[blid] = bgids[i];
152 int_type* b_firstcol =
new int_type[2*numBrows];
153 int_type* b_lastcol = b_firstcol+numBrows;
155 for(i=0; i<numBrows; ++i) {
156 b_firstcol[i] = max_all_b;
157 b_lastcol[i] = min_all_b;
160 if (Blen_i < 1)
continue;
161 int* Bindices_i = Bview.
indices[i];
164 for(k=0; k<Blen_i; ++k) {
165 temp = (int_type) Bview.
importColMap->GID64(Bindices_i[k]);
166 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
167 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
171 for(k=0; k<Blen_i; ++k) {
172 temp = bcols[Bindices_i[k]];
173 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
174 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
181 int_type* Aind = iwork;
182 int_type* Bind = iwork+maxlen;
184 bool C_filled = C.
Filled();
195 for(i=0; i<Aview.
numRows; ++i) {
200 int* Aindices_i = Aview.
indices[i];
201 double* Aval_i = Aview.
values[i];
207 for(k=0; k<A_len_i; ++k) {
208 Aind[k] = (int_type) Aview.
colMap->GID64(Aindices_i[k]);
209 avals[k] = Aval_i[k];
212 util.Sort<int_type>(
true, A_len_i, Aind, 1, &avals, 0, NULL, 0, 0);
214 int_type mina = Aind[0];
215 int_type maxa = Aind[A_len_i-1];
217 if (mina > max_all_b || maxa < min_all_b) {
221 int_type global_row = (int_type) Aview.
rowMap->GID64(i);
224 for(j=0; j<Bview.
numRows; ++j) {
225 if (b_firstcol[j] > maxa || b_lastcol[j] < mina) {
229 int* Bindices_j = Bview.
indices[j];
239 for(k=0; k<B_len_j; ++k) {
240 tmp = (int_type) Bview.
importColMap->GID64(Bindices_j[k]);
241 if (tmp < mina || tmp > maxa) {
245 bvals[Blen] = Bview.
values[j][k];
250 for(k=0; k<B_len_j; ++k) {
251 tmp = bcols[Bindices_j[k]];
252 if (tmp < mina || tmp > maxa) {
256 bvals[Blen] = Bview.
values[j][k];
265 util.Sort<int_type>(
true, Blen, Bind, 1, &bvals, 0, NULL, 0, 0);
267 double C_ij =
sparsedot(avals, Aind, A_len_i,
273 int_type global_col = (int_type) Bview.
rowMap->GID64(j);
288 std::cerr <<
"EpetraExt::MatrixMatrix::Multiply Warning: failed " 289 <<
"to insert value in result matrix at position "<<global_row
290 <<
"," <<global_col<<
", possibly because result matrix has a " 291 <<
"column-map that doesn't include column "<<global_col
292 <<
" on this proc." <<std::endl;
301 delete [] b_firstcol;
310 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 311 if(Aview.
rowMap->GlobalIndicesInt() &&
312 Aview.
colMap->GlobalIndicesInt() &&
313 Aview.
rowMap->GlobalIndicesInt() &&
314 Aview.
colMap->GlobalIndicesInt()) {
315 return mult_A_Btrans<int>(Aview, Bview, C);
319 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 320 if(Aview.
rowMap->GlobalIndicesLongLong() &&
321 Aview.
colMap->GlobalIndicesLongLong() &&
322 Aview.
rowMap->GlobalIndicesLongLong() &&
323 Aview.
colMap->GlobalIndicesLongLong()) {
324 return mult_A_Btrans<long long>(Aview, Bview, C);
328 throw std::runtime_error(
"EpetraExt::mult_A_Btrans: GlobalIndices type unknown");
333 template<
typename int_type>
338 int C_firstCol = Bview.
colMap->MinLID();
339 int C_lastCol = Bview.
colMap->MaxLID();
341 int C_firstCol_import = 0;
342 int C_lastCol_import = -1;
349 int C_numCols = C_lastCol - C_firstCol + 1;
350 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
352 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
354 double* C_row_i =
new double[C_numCols];
355 int_type* C_colInds =
new int_type[C_numCols];
359 for(j=0; j<C_numCols; ++j) {
375 int localProc = Bview.
colMap->Comm().MyPID();
378 Aview.
rowMap->MyGlobalElementsPtr(Arows);
380 bool C_filled = C.
Filled();
383 for(i=0; i<Aview.
numRows; ++i) {
385 int* Aindices_i = Aview.
indices[i];
386 double* Aval_i = Aview.
values[i];
390 int Bi = Bview.
rowMap->LID(Arows[i]);
392 std::cout <<
"mult_Atrans_B ERROR, proc "<<localProc<<
" needs row " 393 <<Arows[i]<<
" of matrix B, but doesn't have it."<<std::endl;
397 int* Bcol_inds = Bview.
indices[Bi];
398 double* Bvals_i = Bview.
values[Bi];
410 for(j=0; j<Blen; ++j) {
411 C_colInds[j] = (int_type) Bview.
importColMap->GID64(Bcol_inds[j]);
415 for(j=0; j<Blen; ++j) {
416 C_colInds[j] = (int_type) Bview.
colMap->GID64(Bcol_inds[j]);
423 int Aj = Aindices_i[j];
424 double Aval = Aval_i[j];
431 global_row = (int_type) Aview.
colMap->GID64(Aj);
434 if (!C.
RowMap().MyGID(global_row)) {
438 for(k=0; k<Blen; ++k) {
439 C_row_i[k] = Aval*Bvals_i[k];
473 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 474 if(Aview.
rowMap->GlobalIndicesInt() &&
475 Aview.
colMap->GlobalIndicesInt() &&
476 Aview.
rowMap->GlobalIndicesInt() &&
477 Aview.
colMap->GlobalIndicesInt()) {
478 return mult_Atrans_B<int>(Aview, Bview, C);
482 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 483 if(Aview.
rowMap->GlobalIndicesLongLong() &&
484 Aview.
colMap->GlobalIndicesLongLong() &&
485 Aview.
rowMap->GlobalIndicesLongLong() &&
486 Aview.
colMap->GlobalIndicesLongLong()) {
487 return mult_Atrans_B<long long>(Aview, Bview, C);
491 throw std::runtime_error(
"EpetraExt::mult_Atrans_B: GlobalIndices type unknown");
495 template<
typename int_type>
500 int C_firstCol = Aview.
rowMap->MinLID();
501 int C_lastCol = Aview.
rowMap->MaxLID();
503 int C_firstCol_import = 0;
504 int C_lastCol_import = -1;
511 int C_numCols = C_lastCol - C_firstCol + 1;
512 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
514 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
516 double* dwork =
new double[C_numCols];
517 int_type* iwork =
new int_type[C_numCols];
519 double* C_col_j = dwork;
520 int_type* C_inds = iwork;
531 for(j=0; j<C_numCols; ++j) {
536 int_type* A_col_inds = 0;
537 Aview.
colMap->MyGlobalElementsPtr(A_col_inds);
538 int_type* A_col_inds_import = 0;
540 Aview.
importColMap->MyGlobalElementsPtr(A_col_inds_import);
542 const Epetra_Map* Crowmap = &(C.
RowMap());
553 Bview.
rowMap->MyGlobalElementsPtr(Brows);
556 for(j=0; j<Bview.
numRows; ++j) {
557 int* Bindices_j = Bview.
indices[j];
558 double* Bvals_j = Bview.
values[j];
560 int_type global_col = Brows[j];
570 int bk = Bindices_j[k];
571 double Bval = Bvals_j[k];
578 global_k = (int_type) Bview.
colMap->GID64(bk);
582 int ak = Aview.
rowMap->LID(global_k);
587 int* Aindices_k = Aview.
indices[ak];
588 double* Avals_k = Aview.
values[ak];
594 C_col_j[C_len] = Avals_k[i]*Bval;
595 C_inds[C_len++] = A_col_inds_import[Aindices_k[i]];
600 C_col_j[C_len] = Avals_k[i]*Bval;
601 C_inds[C_len++] = A_col_inds[Aindices_k[i]];
607 for(i=0; i < C_len ; ++i) {
608 if (C_col_j[i] == 0.0)
continue;
610 int_type global_row = C_inds[i];
611 if (!Crowmap->MyGID(global_row)) {
643 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 644 if(Aview.
rowMap->GlobalIndicesInt() &&
645 Aview.
colMap->GlobalIndicesInt() &&
646 Aview.
rowMap->GlobalIndicesInt() &&
647 Aview.
colMap->GlobalIndicesInt()) {
648 return mult_Atrans_Btrans<int>(Aview, Bview, C);
652 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 653 if(Aview.
rowMap->GlobalIndicesLongLong() &&
654 Aview.
colMap->GlobalIndicesLongLong() &&
655 Aview.
rowMap->GlobalIndicesLongLong() &&
656 Aview.
colMap->GlobalIndicesLongLong()) {
657 return mult_Atrans_Btrans<long long>(Aview, Bview, C);
661 throw std::runtime_error(
"EpetraExt::mult_Atrans_Btrans: GlobalIndices type unknown");
665 template<
typename int_type>
667 const Epetra_Map& targetMap,
669 const Epetra_Import * prototypeImporter=0)
682 #ifdef ENABLE_MMM_TIMINGS 683 using Teuchos::TimeMonitor;
684 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Alloc")));
689 const Epetra_Map& Mrowmap = M.RowMap();
690 int numProcs = Mrowmap.Comm().NumProc();
691 Mview.
numRows = targetMap.NumMyElements();
693 targetMap.MyGlobalElementsPtr(Mrows);
703 Mview.
rowMap = &targetMap;
704 Mview.
colMap = &(M.ColMap());
710 #ifdef ENABLE_MMM_TIMINGS 711 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Extract")));
714 int *rowptr=0, *colind=0;
717 EPETRA_CHK_ERR( M.ExtractCrsDataPointers(rowptr,colind,vals) );
719 if(Mrowmap.SameAs(targetMap)) {
721 for(
int i=0; i<Mview.
numRows; ++i) {
723 Mview.
indices[i] = colind + rowptr[i];
724 Mview.
values[i] = vals + rowptr[i];
729 else if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)){
731 int * PermuteToLIDs = prototypeImporter->PermuteToLIDs();
732 int * PermuteFromLIDs = prototypeImporter->PermuteFromLIDs();
733 int * RemoteLIDs = prototypeImporter->RemoteLIDs();
734 for(
int i=0; i<prototypeImporter->NumSameIDs();i++){
736 Mview.
indices[i] = colind + rowptr[i];
737 Mview.
values[i] = vals + rowptr[i];
740 for(
int i=0; i<prototypeImporter->NumPermuteIDs();i++){
741 int to = PermuteToLIDs[i];
742 int from = PermuteFromLIDs[i];
744 Mview.
indices[to] = colind + rowptr[from];
745 Mview.
values[to] = vals + rowptr[from];
748 for(
int i=0; i<prototypeImporter->NumRemoteIDs();i++){
749 Mview.
remote[RemoteLIDs[i]] =
true;
755 for(
int i=0; i<Mview.
numRows; ++i) {
756 int mlid = Mrowmap.LID(Mrows[i]);
763 Mview.
indices[i] = colind + rowptr[mlid];
764 Mview.
values[i] = vals + rowptr[mlid];
770 #ifdef ENABLE_MMM_TIMINGS 771 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Collective-0")));
776 std::cerr <<
"EpetraExt::MatrixMatrix::Multiply ERROR, numProcs < 2 but " 777 <<
"attempting to import remote matrix rows."<<std::endl;
790 int globalMaxNumRemote = 0;
791 Mrowmap.Comm().MaxAll(&Mview.
numRemote, &globalMaxNumRemote, 1);
793 if (globalMaxNumRemote > 0) {
794 #ifdef ENABLE_MMM_TIMINGS 795 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-1")));
801 for(
int i=0; i<Mview.
numRows; ++i) {
803 MremoteRows[offset++] = Mrows[i];
809 #ifdef ENABLE_MMM_TIMINGS 810 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-2")));
814 Epetra_Import * importer=0;
816 if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)) {
819 else if(!prototypeImporter) {
820 Epetra_Map MremoteRowMap2((int_type) -1, Mview.
numRemote, MremoteRows, (int_type) Mrowmap.IndexBase64(), Mrowmap.Comm());
821 importer=
new Epetra_Import(MremoteRowMap2, Mrowmap);
824 throw std::runtime_error(
"prototypeImporter->SourceMap() does not match M.RowMap()!");
827 #ifdef ENABLE_MMM_TIMINGS 828 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-3")));
835 #ifdef ENABLE_MMM_TIMINGS 836 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-4")));
854 for(
int i=0; i<Mview.
numRows; ++i) {
856 int importLID = MremoteRowMap.LID(Mrows[i]);
858 Mview.
indices[i] = colind + rowptr[importLID];
859 Mview.
values[i] = vals + rowptr[importLID];
864 int_type * MyColGIDs = 0;
869 new Epetra_Map (static_cast<int_type> (-1),
874 delete [] MremoteRows;
875 #ifdef ENABLE_MMM_TIMINGS 890 template<
typename int_type>
892 const Epetra_Map& targetMap,
894 const Epetra_Import * prototypeImporter=0)
905 #ifdef ENABLE_MMM_TIMINGS 906 using Teuchos::TimeMonitor;
907 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM Ionly Setup")));
913 const Epetra_Map& Mrowmap = M.RowMap();
914 int numProcs = Mrowmap.Comm().NumProc();
915 Mview.
numRows = targetMap.NumMyElements();
918 Mview.
rowMap = &targetMap;
919 Mview.
colMap = &(M.ColMap());
925 int N = Mview.
rowMap->NumMyElements();
927 if(Mrowmap.SameAs(targetMap)) {
933 else if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)){
934 numRemote = prototypeImporter->NumRemoteIDs();
939 const int * PermuteToLIDs = prototypeImporter->PermuteToLIDs();
940 const int * PermuteFromLIDs = prototypeImporter->PermuteFromLIDs();
941 const int * RemoteLIDs = prototypeImporter->RemoteLIDs();
943 for(i=0; i<prototypeImporter->NumSameIDs();i++)
947 for(i=0; i<prototypeImporter->NumPermuteIDs();i++)
950 for(i=0; i<prototypeImporter->NumRemoteIDs();i++)
955 throw std::runtime_error(
"import_only: This routine only works if you either have the right map or no prototypeImporter");
959 std::cerr <<
"EpetraExt::MatrixMatrix::Multiply ERROR, numProcs < 2 but " 960 <<
"attempting to import remote matrix rows."<<std::endl;
968 #ifdef ENABLE_MMM_TIMINGS 969 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM Ionly Import-1")));
971 const int * RemoteLIDs = prototypeImporter->RemoteLIDs();
974 int_type* MremoteRows = numRemote>0 ?
new int_type[prototypeImporter->NumRemoteIDs()] : 0;
975 for(i=0; i<prototypeImporter->NumRemoteIDs(); i++)
976 MremoteRows[i] = (int_type) targetMap.GID64(RemoteLIDs[i]);
978 LightweightMap MremoteRowMap((int_type) -1, numRemote, MremoteRows, (int_type)Mrowmap.IndexBase64());
980 #ifdef ENABLE_MMM_TIMINGS 981 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM Ionly Import-2")));
988 #ifdef ENABLE_MMM_TIMINGS 989 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM Ionly Import-3")));
994 #ifdef ENABLE_MMM_TIMINGS 995 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM Ionly Import-4")));
1000 delete [] MremoteRows;
1009 template<
typename int_type>
1011 const Epetra_Map* map2,
1012 const Epetra_Map*& mapunion)
1017 mapunion =
new Epetra_Map(*map2);
1022 mapunion =
new Epetra_Map(*map1);
1026 int map1_len = map1->NumMyElements();
1027 int_type* map1_elements = 0;
1028 map1->MyGlobalElementsPtr(map1_elements);
1029 int map2_len = map2->NumMyElements();
1030 int_type* map2_elements = 0;
1031 map2->MyGlobalElementsPtr(map2_elements);
1033 int_type* union_elements =
new int_type[map1_len+map2_len];
1035 int map1_offset = 0, map2_offset = 0, union_offset = 0;
1037 while(map1_offset < map1_len && map2_offset < map2_len) {
1038 int_type map1_elem = map1_elements[map1_offset];
1039 int_type map2_elem = map2_elements[map2_offset];
1041 if (map1_elem < map2_elem) {
1042 union_elements[union_offset++] = map1_elem;
1045 else if (map1_elem > map2_elem) {
1046 union_elements[union_offset++] = map2_elem;
1050 union_elements[union_offset++] = map1_elem;
1057 for(i=map1_offset; i<map1_len; ++i) {
1058 union_elements[union_offset++] = map1_elements[i];
1061 for(i=map2_offset; i<map2_len; ++i) {
1062 union_elements[union_offset++] = map2_elements[i];
1065 mapunion =
new Epetra_Map((int_type) -1, union_offset, union_elements,
1066 (int_type) map1->IndexBase64(), map1->Comm());
1068 delete [] union_elements;
1074 template<
typename int_type>
1076 const Epetra_Map& column_map)
1083 std::pair<int_type,int_type> my_col_range = get_col_range<int_type>(column_map);
1085 const Epetra_Comm& Comm = M.RowMap().Comm();
1086 int num_procs = Comm.NumProc();
1087 int my_proc = Comm.MyPID();
1088 std::vector<int> send_procs;
1089 send_procs.reserve(num_procs-1);
1090 std::vector<int_type> col_ranges;
1091 col_ranges.reserve(2*(num_procs-1));
1092 for(
int p=0; p<num_procs; ++p) {
1093 if (p == my_proc)
continue;
1094 send_procs.push_back(p);
1095 col_ranges.push_back(my_col_range.first);
1096 col_ranges.push_back(my_col_range.second);
1099 Epetra_Distributor* distor = Comm.CreateDistributor();
1101 int num_recv_procs = 0;
1102 int num_send_procs = send_procs.size();
1103 bool deterministic =
true;
1104 int* send_procs_ptr = send_procs.size()>0 ? &send_procs[0] : NULL;
1105 distor->CreateFromSends(num_send_procs, send_procs_ptr, deterministic, num_recv_procs);
1107 int len_import_chars = 0;
1108 char* import_chars = NULL;
1110 char* export_chars = col_ranges.size()>0 ?
reinterpret_cast<char*
>(&col_ranges[0]) : NULL;
1111 int err = distor->Do(export_chars, 2*
sizeof(int_type), len_import_chars, import_chars);
1113 std::cout <<
"ERROR returned from Distributor::Do."<<std::endl;
1116 int_type* IntImports =
reinterpret_cast<int_type*
>(import_chars);
1117 int num_import_pairs = len_import_chars/(2*
sizeof(int_type));
1119 std::vector<int> send_procs2;
1120 send_procs2.reserve(num_procs);
1121 std::vector<int_type> proc_col_ranges;
1122 std::pair<int_type,int_type> M_col_range = get_col_range<int_type>(M);
1123 for(
int i=0; i<num_import_pairs; ++i) {
1124 int_type first_col = IntImports[offset++];
1125 int_type last_col = IntImports[offset++];
1126 if (first_col <= M_col_range.second && last_col >= M_col_range.first) {
1127 send_procs2.push_back(send_procs[i]);
1128 proc_col_ranges.push_back(first_col);
1129 proc_col_ranges.push_back(last_col);
1133 std::vector<int_type> send_rows;
1134 std::vector<int> rows_per_send_proc;
1137 Epetra_Distributor* distor2 = Comm.CreateDistributor();
1139 int* send_procs2_ptr = send_procs2.size()>0 ? &send_procs2[0] : NULL;
1141 err = distor2->CreateFromSends(send_procs2.size(), send_procs2_ptr, deterministic, num_recv_procs);
1143 export_chars = send_rows.size()>0 ?
reinterpret_cast<char*
>(&send_rows[0]) : NULL;
1144 int* rows_per_send_proc_ptr = rows_per_send_proc.size()>0 ? &rows_per_send_proc[0] : NULL;
1145 len_import_chars = 0;
1146 err = distor2->Do(export_chars, (
int)
sizeof(int_type), rows_per_send_proc_ptr, len_import_chars, import_chars);
1148 int_type* recvd_row_ints =
reinterpret_cast<int_type*
>(import_chars);
1149 int num_recvd_rows = len_import_chars/
sizeof(int_type);
1151 Epetra_Map recvd_rows((int_type) -1, num_recvd_rows, recvd_row_ints, (int_type) column_map.IndexBase64(), Comm);
1155 delete [] import_chars;
1157 Epetra_Map* result_map = NULL;
1159 err = form_map_union<int_type>(&(M.RowMap()), &recvd_rows, (
const Epetra_Map*&)result_map);
1168 const Epetra_Map& column_map)
1170 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1171 if(M.RowMatrixRowMap().GlobalIndicesInt() && column_map.GlobalIndicesInt()) {
1172 return Tfind_rows_containing_cols<int>(M, column_map);
1176 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1177 if(M.RowMatrixRowMap().GlobalIndicesLongLong() && column_map.GlobalIndicesLongLong()) {
1178 return Tfind_rows_containing_cols<long long>(M, column_map);
1182 throw std::runtime_error(
"EpetraExt::find_rows_containing_cols: GlobalIndices type unknown");
1186 template<
typename int_type>
1187 int MatrixMatrix::TMultiply(
const Epetra_CrsMatrix& A,
1189 const Epetra_CrsMatrix& B,
1191 Epetra_CrsMatrix& C,
1192 bool call_FillComplete_on_result)
1196 #ifdef ENABLE_MMM_TIMINGS 1197 using Teuchos::TimeMonitor;
1198 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All Setup")));
1212 if (!A.Filled() || !B.Filled()) {
1217 bool NewFlag=!C.IndicesAreLocal() && !C.IndicesAreGlobal();
1223 if (transposeB && !transposeA) scenario = 2;
1224 if (transposeA && !transposeB) scenario = 3;
1225 if (transposeA && transposeB) scenario = 4;
1226 if(NewFlag && transposeA && !transposeB) scenario = 5;
1229 long long Aouter = transposeA ? A.NumGlobalCols64() : A.NumGlobalRows64();
1230 long long Bouter = transposeB ? B.NumGlobalRows64() : B.NumGlobalCols64();
1231 long long Ainner = transposeA ? A.NumGlobalRows64() : A.NumGlobalCols64();
1232 long long Binner = transposeB ? B.NumGlobalCols64() : B.NumGlobalRows64();
1233 if (Ainner != Binner) {
1234 std::cerr <<
"MatrixMatrix::Multiply: ERROR, inner dimensions of op(A) and op(B) " 1235 <<
"must match for matrix-matrix product. op(A) is " 1236 <<Aouter<<
"x"<<Ainner <<
", op(B) is "<<Binner<<
"x"<<Bouter<<std::endl;
1244 if (Aouter > C.NumGlobalRows64()) {
1245 std::cerr <<
"MatrixMatrix::Multiply: ERROR, dimensions of result C must " 1246 <<
"match dimensions of op(A) * op(B). C has "<<C.NumGlobalRows64()
1247 <<
" rows, should have at least "<<Aouter << std::endl;
1258 int numProcs = A.Comm().NumProc();
1263 const Epetra_Map* domainMap_A = &(A.DomainMap());
1264 const Epetra_Map* domainMap_B = &(B.DomainMap());
1266 const Epetra_Map* rowmap_A = &(A.RowMap());
1267 const Epetra_Map* rowmap_B = &(B.RowMap());
1271 const Epetra_Map* workmap1 = NULL;
1272 const Epetra_Map* workmap2 = NULL;
1273 const Epetra_Map* mapunion1 = NULL;
1280 Teuchos::RCP<Epetra_CrsMatrix> Atrans;
1282 const Epetra_Map* targetMap_A = rowmap_A;
1283 const Epetra_Map* targetMap_B = rowmap_B;
1285 #ifdef ENABLE_MMM_TIMINGS 1286 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All I&X")));
1292 if (scenario == 3 || scenario == 4) {
1293 workmap1 = Tfind_rows_containing_cols<int_type>(A, *domainMap_A);
1294 targetMap_A = workmap1;
1297 if (scenario == 5) {
1298 targetMap_A = &(A.ColMap());
1302 if(scenario==1 && call_FillComplete_on_result) {
1303 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1305 else if (scenario == 5){
1308 Epetra_CrsMatrix * Anonconst =
const_cast<Epetra_CrsMatrix *
>(&A);
1310 import_only<int_type>(*Atrans,*targetMap_A,Atransview);
1313 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1322 const Epetra_Map* colmap_op_A = NULL;
1323 if(scenario==1 || numProcs > 1){
1324 if (transposeA && scenario == 3) {
1325 colmap_op_A = targetMap_A;
1328 colmap_op_A = &(A.ColMap());
1330 targetMap_B = colmap_op_A;
1332 if(scenario==5) targetMap_B = &(B.RowMap());
1341 EPETRA_CHK_ERR( form_map_union<int_type>(colmap_op_A, domainMap_B, mapunion1) );
1342 workmap2 = Tfind_rows_containing_cols<int_type>(B, *mapunion1);
1343 targetMap_B = workmap2;
1348 if((scenario==1 && call_FillComplete_on_result) || scenario==5) {
1349 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1352 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1355 #ifdef ENABLE_MMM_TIMINGS 1356 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All Multiply")));
1360 if(C.Filled()) C.PutScalar(0.0);
1366 case 1: EPETRA_CHK_ERR( mult_A_B(A,Aview,B,Bview,C,call_FillComplete_on_result) );
1368 case 2: EPETRA_CHK_ERR(
mult_A_Btrans(Aview, Bview, ecrsmat) );
1370 case 3: EPETRA_CHK_ERR(
mult_Atrans_B(Aview, Bview, ecrsmat) );
1374 case 5: EPETRA_CHK_ERR( mult_AT_B_newmatrix(Atransview, Bview, C) );
1379 if (scenario != 1 && call_FillComplete_on_result && scenario != 5) {
1387 const Epetra_Map* domainmap =
1388 transposeB ? &(B.RangeMap()) : &(B.DomainMap());
1390 const Epetra_Map* rangemap =
1391 transposeA ? &(A.DomainMap()) : &(A.RangeMap());
1394 EPETRA_CHK_ERR( C.FillComplete(*domainmap, *rangemap) );
1401 delete mapunion1; mapunion1 = NULL;
1402 delete workmap1; workmap1 = NULL;
1403 delete workmap2; workmap2 = NULL;
1410 const Epetra_CrsMatrix& B,
1412 Epetra_CrsMatrix& C,
1413 bool call_FillComplete_on_result)
1415 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1416 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1417 return TMultiply<int>(A, transposeA, B, transposeB, C, call_FillComplete_on_result);
1421 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1422 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1423 return TMultiply<long long>(A, transposeA, B, transposeB, C, call_FillComplete_on_result);
1427 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1431 template<
typename int_type>
1432 int MatrixMatrix::TAdd(
const Epetra_CrsMatrix& A,
1435 Epetra_CrsMatrix& B,
1447 std::cerr <<
"EpetraExt::MatrixMatrix::Add ERROR, input matrix A.Filled() is false, it is required to be true. (Result matrix B is not required to be Filled())."<<std::endl;
1452 Epetra_CrsMatrix * Aprime = 0;
1457 Aprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Atrans)(
const_cast<Epetra_CrsMatrix&
>(A)))));
1460 Aprime =
const_cast<Epetra_CrsMatrix*
>(&A);
1462 int MaxNumEntries = EPETRA_MAX( A.MaxNumEntries(), B.MaxNumEntries() );
1463 int A_NumEntries, B_NumEntries;
1464 int_type * A_Indices =
new int_type[MaxNumEntries];
1465 double * A_Values =
new double[MaxNumEntries];
1466 int* B_Indices_local;
1467 int_type* B_Indices_global;
1470 int NumMyRows = B.NumMyRows();
1478 for(
int i = 0; i < NumMyRows; ++i )
1480 Row = (int_type) B.GRID64(i);
1481 EPETRA_CHK_ERR( Aprime->ExtractGlobalRowCopy( Row, MaxNumEntries, A_NumEntries, A_Values, A_Indices ) );
1483 if (scalarB != 1.0) {
1485 EPETRA_CHK_ERR( B.ExtractGlobalRowView( Row, B_NumEntries,
1486 B_Values, B_Indices_global));
1489 EPETRA_CHK_ERR( B.ExtractMyRowView( i, B_NumEntries,
1490 B_Values, B_Indices_local));
1493 for(
int jj=0; jj<B_NumEntries; ++jj) {
1494 B_Values[jj] = scalarB*B_Values[jj];
1498 if( scalarA != 1.0 ) {
1499 for(
int j = 0; j < A_NumEntries; ++j ) A_Values[j] *= scalarA;
1503 err = B.SumIntoGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1505 if (err < 0) ierr = err;
1508 err = B.InsertGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1509 assert( err == 0 || err == 1 || err == 3 );
1510 if (err < 0) ierr = err;
1515 EPETRA_CHK_ERR( B.Scale(scalarB) );
1518 delete [] A_Indices;
1521 if( Atrans )
delete Atrans;
1529 Epetra_CrsMatrix& B,
1532 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1533 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1534 return TAdd<int>(A, transposeA, scalarA, B, scalarB);
1538 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1539 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1540 return TAdd<long long>(A, transposeA, scalarA, B, scalarB);
1544 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1547 template<
typename int_type>
1548 int MatrixMatrix::TAdd(
const Epetra_CrsMatrix& A,
1551 const Epetra_CrsMatrix & B,
1554 Epetra_CrsMatrix * & C)
1561 if (!A.Filled() || !B.Filled() ) {
1562 std::cerr <<
"EpetraExt::MatrixMatrix::Add ERROR, input matrix A.Filled() or B.Filled() is false," 1563 <<
"they are required to be true. (Result matrix C should be an empty pointer)" << std::endl;
1567 Epetra_CrsMatrix * Aprime = 0, * Bprime=0;
1573 Aprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Atrans)(
const_cast<Epetra_CrsMatrix&
>(A)))));
1576 Aprime =
const_cast<Epetra_CrsMatrix*
>(&A);
1581 Bprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Btrans)(
const_cast<Epetra_CrsMatrix&
>(B)))));
1584 Bprime =
const_cast<Epetra_CrsMatrix*
>(&B);
1590 C =
new Epetra_CrsMatrix(Copy,Aprime->RowMap(),0);
1594 Epetra_CrsMatrix * Mat[] = { Aprime,Bprime};
1595 double scalar[] = { scalarA, scalarB};
1598 for(
int k=0;k<2;k++) {
1599 int MaxNumEntries = Mat[k]->MaxNumEntries();
1601 int_type * Indices =
new int_type[MaxNumEntries];
1602 double * Values =
new double[MaxNumEntries];
1604 int NumMyRows = Mat[k]->NumMyRows();
1609 for(
int i = 0; i < NumMyRows; ++i ) {
1610 Row = (int_type) Mat[k]->GRID64(i);
1611 EPETRA_CHK_ERR( Mat[k]->ExtractGlobalRowCopy( Row, MaxNumEntries, NumEntries, Values, Indices));
1613 if( scalar[k] != 1.0 )
1614 for(
int j = 0; j < NumEntries; ++j ) Values[j] *= scalar[k];
1617 err = C->SumIntoGlobalValues( Row, NumEntries, Values, Indices );
1618 if (err < 0) ierr = err;
1620 err = C->InsertGlobalValues( Row, NumEntries, Values, Indices );
1621 if (err < 0) ierr = err;
1629 if( Atrans )
delete Atrans;
1630 if( Btrans )
delete Btrans;
1638 const Epetra_CrsMatrix & B,
1641 Epetra_CrsMatrix * & C)
1643 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1644 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1645 return TAdd<int>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1649 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1650 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1651 return TAdd<long long>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1655 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1661 template<
typename int_type>
1662 int MatrixMatrix::TJacobi(
double omega,
1663 const Epetra_Vector & Dinv,
1664 const Epetra_CrsMatrix& A,
1665 const Epetra_CrsMatrix& B,
1666 Epetra_CrsMatrix& C,
1667 bool call_FillComplete_on_result)
1669 #ifdef ENABLE_MMM_TIMINGS 1670 using Teuchos::TimeMonitor;
1671 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All Setup")));
1675 if (!A.Filled() || !B.Filled()) {
1680 long long Aouter = A.NumGlobalRows64();
1681 long long Bouter = B.NumGlobalCols64();
1682 long long Ainner = A.NumGlobalCols64();
1683 long long Binner = B.NumGlobalRows64();
1684 long long Dlen = Dinv.GlobalLength64();
1685 if (Ainner != Binner) {
1686 std::cerr <<
"MatrixMatrix::Jacobi: ERROR, inner dimensions of A and B " 1687 <<
"must match for matrix-matrix product. A is " 1688 <<Aouter<<
"x"<<Ainner <<
", B is "<<Binner<<
"x"<<Bouter<<std::endl;
1696 if (Aouter > C.NumGlobalRows64()) {
1697 std::cerr <<
"MatrixMatrix::Jacobi: ERROR, dimensions of result C must " 1698 <<
"match dimensions of A * B. C has "<<C.NumGlobalRows64()
1699 <<
" rows, should have at least "<<Aouter << std::endl;
1704 if(Dlen != Aouter) {
1705 std::cerr <<
"MatrixMatrix::Jacboi: ERROR, dimensions of result D must " 1706 <<
"match dimensions of A's rows. D has "<< Dlen
1707 <<
" rows, should have " << Aouter << std::endl;
1711 if(!A.RowMap().SameAs(B.RowMap()) || !A.RowMap().SameAs(Dinv.Map())) {
1712 std::cerr <<
"MatrixMatrix::Jacboi: ERROR, RowMap of A must match RowMap of B " 1713 <<
"and Map of D."<<std::endl;
1724 int numProcs = A.Comm().NumProc();
1727 const Epetra_Map* rowmap_A = &(A.RowMap());
1728 const Epetra_Map* rowmap_B = &(B.RowMap());
1734 const Epetra_Map* workmap1 = NULL;
1735 const Epetra_Map* workmap2 = NULL;
1736 const Epetra_Map* mapunion1 = NULL;
1743 const Epetra_Map* targetMap_A = rowmap_A;
1744 const Epetra_Map* targetMap_B = rowmap_B;
1746 #ifdef ENABLE_MMM_TIMINGS 1747 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All I&X")));
1751 if(call_FillComplete_on_result) {
1752 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1755 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1762 const Epetra_Map* colmap_op_A = NULL;
1764 colmap_op_A = &(A.ColMap());
1765 targetMap_B = colmap_op_A;
1769 if(call_FillComplete_on_result) {
1770 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1773 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1776 #ifdef ENABLE_MMM_TIMINGS 1777 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All Multiply")));
1781 if(C.Filled()) C.PutScalar(0.0);
1785 EPETRA_CHK_ERR( jacobi_A_B(omega,Dinv,A,Aview,B,Bview,C,call_FillComplete_on_result) );
1789 delete mapunion1; mapunion1 = NULL;
1790 delete workmap1; workmap1 = NULL;
1791 delete workmap2; workmap2 = NULL;
1799 const Epetra_Vector & Dinv,
1800 const Epetra_CrsMatrix& A,
1801 const Epetra_CrsMatrix& B,
1802 Epetra_CrsMatrix& C,
1803 bool call_FillComplete_on_result)
1805 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1806 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1807 return TJacobi<int>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1811 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1812 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1813 return TJacobi<long long>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1817 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Jacobi: GlobalIndices type unknown");
LightweightCrsMatrix * importMatrix
int mult_Atrans_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C)
std::vector< int > targetMapToOrigRow
Epetra_BlockMap * RowMapEP_
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
static int Multiply(const Epetra_CrsMatrix &A, bool transposeA, const Epetra_CrsMatrix &B, bool transposeB, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true)
Given Epetra_CrsMatrix objects A, B and C, form the product C = A*B.
LightweightMap * RowMapLW_
void pack_outgoing_rows(const Epetra_CrsMatrix &mtx, const std::vector< int > &proc_col_ranges, std::vector< int > &send_rows, std::vector< int > &rows_per_send_proc)
const Epetra_Map * rowMap
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
std::vector< int > colind_
const Epetra_CrsMatrix * origMatrix
Transform to form the explicit transpose of a Epetra_RowMatrix.
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
Epetra_Map * find_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
int mult_A_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C)
int form_map_union(const Epetra_Map *map1, const Epetra_Map *map2, const Epetra_Map *&mapunion)
const Epetra_Map * origRowMap
int import_and_extract_views(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter=0)
std::vector< int > targetMapToImportRow
virtual const Epetra_Map & RowMap() const =0
int NumMyElements() const
int mult_Atrans_B(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C)
const Epetra_Map * domainMap
std::vector< double > vals_
static int Add(const Epetra_CrsMatrix &A, bool transposeA, double scalarA, Epetra_CrsMatrix &B, double scalarB)
Given Epetra_CrsMatrix objects A and B, form the sum B = a*A + b*B.
long long IndexBase64() const
double sparsedot(double *u, int_type *u_ind, int u_len, double *v, int_type *v_ind, int v_len)
Method for internal use...
static int Jacobi(double omega, const Epetra_Vector &Dinv, const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true)
Given Epetra_CrsMatrix objects A, B and C, and Epetra_Vector Dinv, form the product C = (I-omega * Di...
int import_only(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter=0)
Epetra_Map * Tfind_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
Epetra_CrsMatrix * CreateTransposeLocal(OriginalTypeRef orig)
Local-only transpose operator. Don't use this unless you're sure you know what you're doing...
std::vector< int > rowptr_
const Epetra_Map * colMap
const Epetra_BlockMap * importColMap