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> 72 template<
int>
int import_only(
const Epetra_CrsMatrix& M,
73 const Epetra_Map& targetMap,
74 CrsMatrixStruct& Mview,
75 const Epetra_Import * prototypeImporter);
77 template<
long long>
int import_only(
const Epetra_CrsMatrix& M,
78 const Epetra_Map& targetMap,
79 CrsMatrixStruct& Mview,
80 const Epetra_Import * prototypeImporter);
89 template<
typename int_type>
90 double sparsedot(
double* u, int_type* u_ind,
int u_len,
91 double* v, int_type* v_ind,
int v_len)
98 while(v_idx < v_len && u_idx < u_len) {
99 int_type ui = u_ind[u_idx];
100 int_type vi = v_ind[v_idx];
109 result += u[u_idx++]*v[v_idx++];
118 template<
typename int_type>
122 bool keep_all_hard_zeros)
128 for(i=0; i<Aview.
numRows; ++i) {
131 for(i=0; i<Bview.
numRows; ++i) {
141 int numBcols = Bview.
colMap->NumMyElements();
144 int iworklen = maxlen*2 + numBcols;
145 int_type* iwork =
new int_type[iworklen];
147 int_type * bcols = iwork+maxlen*2;
149 Bview.
colMap->MyGlobalElementsPtr(bgids);
150 double* bvals =
new double[maxlen*2];
151 double* avals = bvals+maxlen;
153 int_type max_all_b = (int_type) Bview.
colMap->MaxAllGID64();
154 int_type min_all_b = (int_type) Bview.
colMap->MinAllGID64();
158 for(i=0; i<numBcols; ++i) {
159 int blid = Bview.
colMap->LID(bgids[i]);
160 bcols[blid] = bgids[i];
167 int_type* b_firstcol =
new int_type[2*numBrows];
168 int_type* b_lastcol = b_firstcol+numBrows;
170 for(i=0; i<numBrows; ++i) {
171 b_firstcol[i] = max_all_b;
172 b_lastcol[i] = min_all_b;
175 if (Blen_i < 1)
continue;
176 int* Bindices_i = Bview.
indices[i];
179 for(k=0; k<Blen_i; ++k) {
180 temp = (int_type) Bview.
importColMap->GID64(Bindices_i[k]);
181 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
182 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
186 for(k=0; k<Blen_i; ++k) {
187 temp = bcols[Bindices_i[k]];
188 if (temp < b_firstcol[i]) b_firstcol[i] = temp;
189 if (temp > b_lastcol[i]) b_lastcol[i] = temp;
196 int_type* Aind = iwork;
197 int_type* Bind = iwork+maxlen;
199 bool C_filled = C.
Filled();
210 for(i=0; i<Aview.
numRows; ++i) {
215 int* Aindices_i = Aview.
indices[i];
216 double* Aval_i = Aview.
values[i];
222 for(k=0; k<A_len_i; ++k) {
223 Aind[k] = (int_type) Aview.
colMap->GID64(Aindices_i[k]);
224 avals[k] = Aval_i[k];
227 util.Sort<int_type>(
true, A_len_i, Aind, 1, &avals, 0, NULL, 0, 0);
229 int_type mina = Aind[0];
230 int_type maxa = Aind[A_len_i-1];
232 if (mina > max_all_b || maxa < min_all_b) {
236 int_type global_row = (int_type) Aview.
rowMap->GID64(i);
239 for(j=0; j<Bview.
numRows; ++j) {
240 if (b_firstcol[j] > maxa || b_lastcol[j] < mina) {
244 int* Bindices_j = Bview.
indices[j];
254 for(k=0; k<B_len_j; ++k) {
255 tmp = (int_type) Bview.
importColMap->GID64(Bindices_j[k]);
256 if (tmp < mina || tmp > maxa) {
260 bvals[Blen] = Bview.
values[j][k];
265 for(k=0; k<B_len_j; ++k) {
266 tmp = bcols[Bindices_j[k]];
267 if (tmp < mina || tmp > maxa) {
271 bvals[Blen] = Bview.
values[j][k];
280 util.Sort<int_type>(
true, Blen, Bind, 1, &bvals, 0, NULL, 0, 0);
282 double C_ij =
sparsedot(avals, Aind, A_len_i,
285 if (!keep_all_hard_zeros && C_ij == 0.0)
288 int_type global_col = (int_type) Bview.
rowMap->GID64(j);
303 std::cerr <<
"EpetraExt::MatrixMatrix::Multiply Warning: failed " 304 <<
"to insert value in result matrix at position "<<global_row
305 <<
"," <<global_col<<
", possibly because result matrix has a " 306 <<
"column-map that doesn't include column "<<global_col
307 <<
" on this proc." <<std::endl;
316 delete [] b_firstcol;
322 CrsMatrixStruct& Bview,
324 bool keep_all_hard_zeros)
326 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 327 if(Aview.rowMap->GlobalIndicesInt() &&
328 Aview.colMap->GlobalIndicesInt() &&
329 Aview.rowMap->GlobalIndicesInt() &&
330 Aview.colMap->GlobalIndicesInt()) {
331 return mult_A_Btrans<int>(Aview, Bview, C, keep_all_hard_zeros);
335 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 336 if(Aview.rowMap->GlobalIndicesLongLong() &&
337 Aview.colMap->GlobalIndicesLongLong() &&
338 Aview.rowMap->GlobalIndicesLongLong() &&
339 Aview.colMap->GlobalIndicesLongLong()) {
340 return mult_A_Btrans<long long>(Aview, Bview, C, keep_all_hard_zeros);
344 throw std::runtime_error(
"EpetraExt::mult_A_Btrans: GlobalIndices type unknown");
349 template<
typename int_type>
355 int C_firstCol = Bview.
colMap->MinLID();
356 int C_lastCol = Bview.
colMap->MaxLID();
358 int C_firstCol_import = 0;
359 int C_lastCol_import = -1;
366 int C_numCols = C_lastCol - C_firstCol + 1;
367 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
369 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
371 double* C_row_i =
new double[C_numCols];
372 int_type* C_colInds =
new int_type[C_numCols];
376 for(j=0; j<C_numCols; ++j) {
392 int localProc = Bview.
colMap->Comm().MyPID();
395 Aview.
rowMap->MyGlobalElementsPtr(Arows);
397 bool C_filled = C.
Filled();
400 for(i=0; i<Aview.
numRows; ++i) {
402 int* Aindices_i = Aview.
indices[i];
403 double* Aval_i = Aview.
values[i];
407 int Bi = Bview.
rowMap->LID(Arows[i]);
409 std::cout <<
"mult_Atrans_B ERROR, proc "<<localProc<<
" needs row " 410 <<Arows[i]<<
" of matrix B, but doesn't have it."<<std::endl;
414 int* Bcol_inds = Bview.
indices[Bi];
415 double* Bvals_i = Bview.
values[Bi];
427 for(j=0; j<Blen; ++j) {
428 C_colInds[j] = (int_type) Bview.
importColMap->GID64(Bcol_inds[j]);
432 for(j=0; j<Blen; ++j) {
433 C_colInds[j] = (int_type) Bview.
colMap->GID64(Bcol_inds[j]);
440 int Aj = Aindices_i[j];
441 double Aval = Aval_i[j];
448 global_row = (int_type) Aview.
colMap->GID64(Aj);
451 if (!C.
RowMap().MyGID(global_row)) {
455 for(k=0; k<Blen; ++k) {
456 C_row_i[k] = Aval*Bvals_i[k];
487 CrsMatrixStruct& Bview,
490 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 491 if(Aview.rowMap->GlobalIndicesInt() &&
492 Aview.colMap->GlobalIndicesInt() &&
493 Aview.rowMap->GlobalIndicesInt() &&
494 Aview.colMap->GlobalIndicesInt()) {
495 return mult_Atrans_B<int>(Aview, Bview, C);
499 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 500 if(Aview.rowMap->GlobalIndicesLongLong() &&
501 Aview.colMap->GlobalIndicesLongLong() &&
502 Aview.rowMap->GlobalIndicesLongLong() &&
503 Aview.colMap->GlobalIndicesLongLong()) {
504 return mult_Atrans_B<long long>(Aview, Bview, C);
508 throw std::runtime_error(
"EpetraExt::mult_Atrans_B: GlobalIndices type unknown");
512 template<
typename int_type>
516 bool keep_all_hard_zeros)
518 int C_firstCol = Aview.
rowMap->MinLID();
519 int C_lastCol = Aview.
rowMap->MaxLID();
521 int C_firstCol_import = 0;
522 int C_lastCol_import = -1;
529 int C_numCols = C_lastCol - C_firstCol + 1;
530 int C_numCols_import = C_lastCol_import - C_firstCol_import + 1;
532 if (C_numCols_import > C_numCols) C_numCols = C_numCols_import;
534 double* dwork =
new double[C_numCols];
535 int_type* iwork =
new int_type[C_numCols];
537 double* C_col_j = dwork;
538 int_type* C_inds = iwork;
542 for(j=0; j<C_numCols; ++j) {
547 int_type* A_col_inds = 0;
548 Aview.
colMap->MyGlobalElementsPtr(A_col_inds);
549 int_type* A_col_inds_import = 0;
551 Aview.
importColMap->MyGlobalElementsPtr(A_col_inds_import);
553 const Epetra_Map* Crowmap = &(C.
RowMap());
564 Bview.
rowMap->MyGlobalElementsPtr(Brows);
567 for(j=0; j<Bview.
numRows; ++j) {
568 int* Bindices_j = Bview.
indices[j];
569 double* Bvals_j = Bview.
values[j];
571 int_type global_col = Brows[j];
581 int bk = Bindices_j[k];
582 double Bval = Bvals_j[k];
589 global_k = (int_type) Bview.
colMap->GID64(bk);
593 int ak = Aview.
rowMap->LID(global_k);
598 int* Aindices_k = Aview.
indices[ak];
599 double* Avals_k = Aview.
values[ak];
605 C_col_j[C_len] = Avals_k[i]*Bval;
606 C_inds[C_len++] = A_col_inds_import[Aindices_k[i]];
611 C_col_j[C_len] = Avals_k[i]*Bval;
612 C_inds[C_len++] = A_col_inds[Aindices_k[i]];
618 for(i=0; i < C_len ; ++i) {
619 if (!keep_all_hard_zeros && C_col_j[i] == 0.0)
continue;
621 int_type global_row = C_inds[i];
622 if (!Crowmap->MyGID(global_row)) {
651 CrsMatrixStruct& Bview,
653 bool keep_all_hard_zeros)
655 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 656 if(Aview.rowMap->GlobalIndicesInt() &&
657 Aview.colMap->GlobalIndicesInt() &&
658 Aview.rowMap->GlobalIndicesInt() &&
659 Aview.colMap->GlobalIndicesInt()) {
660 return mult_Atrans_Btrans<int>(Aview, Bview, C, keep_all_hard_zeros);
664 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 665 if(Aview.rowMap->GlobalIndicesLongLong() &&
666 Aview.colMap->GlobalIndicesLongLong() &&
667 Aview.rowMap->GlobalIndicesLongLong() &&
668 Aview.colMap->GlobalIndicesLongLong()) {
669 return mult_Atrans_Btrans<long long>(Aview, Bview, C, keep_all_hard_zeros);
673 throw std::runtime_error(
"EpetraExt::mult_Atrans_Btrans: GlobalIndices type unknown");
677 template<
typename int_type>
679 const Epetra_Map& targetMap,
681 const Epetra_Import * prototypeImporter=0)
694 #ifdef ENABLE_MMM_TIMINGS 695 using Teuchos::TimeMonitor;
696 Teuchos::RCP<Teuchos::TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Alloc")));
701 const Epetra_Map& Mrowmap = M.RowMap();
702 int numProcs = Mrowmap.Comm().NumProc();
703 Mview.
numRows = targetMap.NumMyElements();
705 targetMap.MyGlobalElementsPtr(Mrows);
715 Mview.
rowMap = &targetMap;
716 Mview.
colMap = &(M.ColMap());
722 #ifdef ENABLE_MMM_TIMINGS 723 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Extract")));
726 int *rowptr=0, *colind=0;
729 EPETRA_CHK_ERR( M.ExtractCrsDataPointers(rowptr,colind,vals) );
731 if(Mrowmap.SameAs(targetMap)) {
733 for(
int i=0; i<Mview.
numRows; ++i) {
735 Mview.
indices[i] = colind + rowptr[i];
736 Mview.
values[i] = vals + rowptr[i];
741 else if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)){
743 int * PermuteToLIDs = prototypeImporter->PermuteToLIDs();
744 int * PermuteFromLIDs = prototypeImporter->PermuteFromLIDs();
745 int * RemoteLIDs = prototypeImporter->RemoteLIDs();
746 for(
int i=0; i<prototypeImporter->NumSameIDs();i++){
748 Mview.
indices[i] = colind + rowptr[i];
749 Mview.
values[i] = vals + rowptr[i];
752 for(
int i=0; i<prototypeImporter->NumPermuteIDs();i++){
753 int to = PermuteToLIDs[i];
754 int from = PermuteFromLIDs[i];
756 Mview.
indices[to] = colind + rowptr[from];
757 Mview.
values[to] = vals + rowptr[from];
760 for(
int i=0; i<prototypeImporter->NumRemoteIDs();i++){
761 Mview.
remote[RemoteLIDs[i]] =
true;
767 for(
int i=0; i<Mview.
numRows; ++i) {
768 int mlid = Mrowmap.LID(Mrows[i]);
775 Mview.
indices[i] = colind + rowptr[mlid];
776 Mview.
values[i] = vals + rowptr[mlid];
782 #ifdef ENABLE_MMM_TIMINGS 783 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Collective-0")));
788 std::cerr <<
"EpetraExt::MatrixMatrix::Multiply ERROR, numProcs < 2 but " 789 <<
"attempting to import remote matrix rows."<<std::endl;
802 int globalMaxNumRemote = 0;
803 Mrowmap.Comm().MaxAll(&Mview.
numRemote, &globalMaxNumRemote, 1);
805 if (globalMaxNumRemote > 0) {
806 #ifdef ENABLE_MMM_TIMINGS 807 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-1")));
813 for(
int i=0; i<Mview.
numRows; ++i) {
815 MremoteRows[offset++] = Mrows[i];
821 #ifdef ENABLE_MMM_TIMINGS 822 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-2")));
826 Epetra_Import * importer=0;
828 if(prototypeImporter && prototypeImporter->SourceMap().SameAs(M.RowMap()) && prototypeImporter->TargetMap().SameAs(targetMap)) {
831 else if(!prototypeImporter) {
832 Epetra_Map MremoteRowMap2((int_type) -1, Mview.
numRemote, MremoteRows, (int_type) Mrowmap.IndexBase64(), Mrowmap.Comm());
833 importer=
new Epetra_Import(MremoteRowMap2, Mrowmap);
836 throw std::runtime_error(
"prototypeImporter->SourceMap() does not match M.RowMap()!");
839 #ifdef ENABLE_MMM_TIMINGS 840 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-3")));
847 #ifdef ENABLE_MMM_TIMINGS 848 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM I&X Import-4")));
866 for(
int i=0; i<Mview.
numRows; ++i) {
868 int importLID = MremoteRowMap.LID(Mrows[i]);
870 Mview.
indices[i] = colind + rowptr[importLID];
871 Mview.
values[i] = vals + rowptr[importLID];
876 int_type * MyColGIDs = 0;
881 new Epetra_Map (static_cast<int_type> (-1),
886 delete [] MremoteRows;
887 #ifdef ENABLE_MMM_TIMINGS 904 template<
typename int_type>
906 const Epetra_Map* map2,
907 const Epetra_Map*& mapunion)
912 mapunion =
new Epetra_Map(*map2);
917 mapunion =
new Epetra_Map(*map1);
921 int map1_len = map1->NumMyElements();
922 int_type* map1_elements = 0;
923 map1->MyGlobalElementsPtr(map1_elements);
924 int map2_len = map2->NumMyElements();
925 int_type* map2_elements = 0;
926 map2->MyGlobalElementsPtr(map2_elements);
928 int_type* union_elements =
new int_type[map1_len+map2_len];
930 int map1_offset = 0, map2_offset = 0, union_offset = 0;
932 while(map1_offset < map1_len && map2_offset < map2_len) {
933 int_type map1_elem = map1_elements[map1_offset];
934 int_type map2_elem = map2_elements[map2_offset];
936 if (map1_elem < map2_elem) {
937 union_elements[union_offset++] = map1_elem;
940 else if (map1_elem > map2_elem) {
941 union_elements[union_offset++] = map2_elem;
945 union_elements[union_offset++] = map1_elem;
952 for(i=map1_offset; i<map1_len; ++i) {
953 union_elements[union_offset++] = map1_elements[i];
956 for(i=map2_offset; i<map2_len; ++i) {
957 union_elements[union_offset++] = map2_elements[i];
960 mapunion =
new Epetra_Map((int_type) -1, union_offset, union_elements,
961 (int_type) map1->IndexBase64(), map1->Comm());
963 delete [] union_elements;
969 template<
typename int_type>
971 const Epetra_Map& column_map)
978 std::pair<int_type,int_type> my_col_range = get_col_range<int_type>(column_map);
980 const Epetra_Comm& Comm = M.RowMap().Comm();
981 int num_procs = Comm.NumProc();
982 int my_proc = Comm.MyPID();
983 std::vector<int> send_procs;
984 send_procs.reserve(num_procs-1);
985 std::vector<int_type> col_ranges;
986 col_ranges.reserve(2*(num_procs-1));
987 for(
int p=0; p<num_procs; ++p) {
988 if (p == my_proc)
continue;
989 send_procs.push_back(p);
990 col_ranges.push_back(my_col_range.first);
991 col_ranges.push_back(my_col_range.second);
994 Epetra_Distributor* distor = Comm.CreateDistributor();
996 int num_recv_procs = 0;
997 int num_send_procs = send_procs.size();
998 bool deterministic =
true;
999 int* send_procs_ptr = send_procs.size()>0 ? &send_procs[0] : NULL;
1000 distor->CreateFromSends(num_send_procs, send_procs_ptr, deterministic, num_recv_procs);
1002 int len_import_chars = 0;
1003 char* import_chars = NULL;
1005 char* export_chars = col_ranges.size()>0 ?
reinterpret_cast<char*
>(&col_ranges[0]) : NULL;
1006 int err = distor->Do(export_chars, 2*
sizeof(int_type), len_import_chars, import_chars);
1008 std::cout <<
"ERROR returned from Distributor::Do."<<std::endl;
1011 int_type* IntImports =
reinterpret_cast<int_type*
>(import_chars);
1012 int num_import_pairs = len_import_chars/(2*
sizeof(int_type));
1014 std::vector<int> send_procs2;
1015 send_procs2.reserve(num_procs);
1016 std::vector<int_type> proc_col_ranges;
1017 std::pair<int_type,int_type> M_col_range = get_col_range<int_type>(M);
1018 for(
int i=0; i<num_import_pairs; ++i) {
1019 int_type first_col = IntImports[offset++];
1020 int_type last_col = IntImports[offset++];
1021 if (first_col <= M_col_range.second && last_col >= M_col_range.first) {
1022 send_procs2.push_back(send_procs[i]);
1023 proc_col_ranges.push_back(first_col);
1024 proc_col_ranges.push_back(last_col);
1028 std::vector<int_type> send_rows;
1029 std::vector<int> rows_per_send_proc;
1032 Epetra_Distributor* distor2 = Comm.CreateDistributor();
1034 int* send_procs2_ptr = send_procs2.size()>0 ? &send_procs2[0] : NULL;
1036 err = distor2->CreateFromSends(send_procs2.size(), send_procs2_ptr, deterministic, num_recv_procs);
1038 export_chars = send_rows.size()>0 ?
reinterpret_cast<char*
>(&send_rows[0]) : NULL;
1039 int* rows_per_send_proc_ptr = rows_per_send_proc.size()>0 ? &rows_per_send_proc[0] : NULL;
1040 len_import_chars = 0;
1041 err = distor2->Do(export_chars, (
int)
sizeof(int_type), rows_per_send_proc_ptr, len_import_chars, import_chars);
1043 int_type* recvd_row_ints =
reinterpret_cast<int_type*
>(import_chars);
1044 int num_recvd_rows = len_import_chars/
sizeof(int_type);
1046 Epetra_Map recvd_rows((int_type) -1, num_recvd_rows, recvd_row_ints, (int_type) column_map.IndexBase64(), Comm);
1050 delete [] import_chars;
1052 Epetra_Map* result_map = NULL;
1054 err = form_map_union<int_type>(&(M.RowMap()), &recvd_rows, (
const Epetra_Map*&)result_map);
1063 const Epetra_Map& column_map)
1065 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1066 if(M.RowMatrixRowMap().GlobalIndicesInt() && column_map.GlobalIndicesInt()) {
1067 return Tfind_rows_containing_cols<int>(M, column_map);
1071 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1072 if(M.RowMatrixRowMap().GlobalIndicesLongLong() && column_map.GlobalIndicesLongLong()) {
1073 return Tfind_rows_containing_cols<long long>(M, column_map);
1077 throw std::runtime_error(
"EpetraExt::find_rows_containing_cols: GlobalIndices type unknown");
1081 template<
typename int_type>
1082 int MatrixMatrix::TMultiply(
const Epetra_CrsMatrix& A,
1084 const Epetra_CrsMatrix& B,
1086 Epetra_CrsMatrix& C,
1087 bool call_FillComplete_on_result,
1088 bool keep_all_hard_zeros)
1092 #ifdef ENABLE_MMM_TIMINGS 1093 using Teuchos::TimeMonitor;
1094 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All Setup")));
1108 if (!A.Filled() || !B.Filled()) {
1113 bool NewFlag=!C.IndicesAreLocal() && !C.IndicesAreGlobal();
1119 if (transposeB && !transposeA) scenario = 2;
1120 if (transposeA && !transposeB) scenario = 3;
1121 if (transposeA && transposeB) scenario = 4;
1122 if(NewFlag && call_FillComplete_on_result && transposeA && !transposeB) scenario = 5;
1125 long long Aouter = transposeA ? A.NumGlobalCols64() : A.NumGlobalRows64();
1126 long long Bouter = transposeB ? B.NumGlobalRows64() : B.NumGlobalCols64();
1127 long long Ainner = transposeA ? A.NumGlobalRows64() : A.NumGlobalCols64();
1128 long long Binner = transposeB ? B.NumGlobalCols64() : B.NumGlobalRows64();
1129 if (Ainner != Binner) {
1130 std::cerr <<
"MatrixMatrix::Multiply: ERROR, inner dimensions of op(A) and op(B) " 1131 <<
"must match for matrix-matrix product. op(A) is " 1132 <<Aouter<<
"x"<<Ainner <<
", op(B) is "<<Binner<<
"x"<<Bouter<<std::endl;
1140 if (Aouter > C.NumGlobalRows64()) {
1141 std::cerr <<
"MatrixMatrix::Multiply: ERROR, dimensions of result C must " 1142 <<
"match dimensions of op(A) * op(B). C has "<<C.NumGlobalRows64()
1143 <<
" rows, should have at least "<<Aouter << std::endl;
1154 int numProcs = A.Comm().NumProc();
1159 const Epetra_Map* domainMap_A = &(A.DomainMap());
1160 const Epetra_Map* domainMap_B = &(B.DomainMap());
1162 const Epetra_Map* rowmap_A = &(A.RowMap());
1163 const Epetra_Map* rowmap_B = &(B.RowMap());
1167 const Epetra_Map* workmap1 = NULL;
1168 const Epetra_Map* workmap2 = NULL;
1169 const Epetra_Map* mapunion1 = NULL;
1173 CrsMatrixStruct Aview;
1174 CrsMatrixStruct Atransview;
1175 CrsMatrixStruct Bview;
1176 Teuchos::RCP<Epetra_CrsMatrix> Atrans;
1178 const Epetra_Map* targetMap_A = rowmap_A;
1179 const Epetra_Map* targetMap_B = rowmap_B;
1181 #ifdef ENABLE_MMM_TIMINGS 1182 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All I&X")));
1188 if (scenario == 3 || scenario == 4) {
1189 workmap1 = Tfind_rows_containing_cols<int_type>(A, *domainMap_A);
1190 targetMap_A = workmap1;
1193 if (scenario == 5) {
1194 targetMap_A = &(A.ColMap());
1198 if(scenario==1 && call_FillComplete_on_result) {
1199 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1201 else if (scenario == 5){
1204 Epetra_CrsMatrix * Anonconst =
const_cast<Epetra_CrsMatrix *
>(&A);
1205 Atrans = Teuchos::rcp(at.CreateTransposeLocal(*Anonconst));
1206 import_only<int_type>(*Atrans,*targetMap_A,Atransview);
1209 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1218 const Epetra_Map* colmap_op_A = NULL;
1219 if(scenario==1 || numProcs > 1){
1220 if (transposeA && scenario == 3) {
1221 colmap_op_A = targetMap_A;
1224 colmap_op_A = &(A.ColMap());
1226 targetMap_B = colmap_op_A;
1228 if(scenario==5) targetMap_B = &(B.RowMap());
1237 EPETRA_CHK_ERR( form_map_union<int_type>(colmap_op_A, domainMap_B, mapunion1) );
1238 workmap2 = Tfind_rows_containing_cols<int_type>(B, *mapunion1);
1239 targetMap_B = workmap2;
1244 if((scenario==1 && call_FillComplete_on_result) || scenario==5) {
1245 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1248 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1251 #ifdef ENABLE_MMM_TIMINGS 1252 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: MMM All Multiply")));
1256 if(C.Filled()) C.PutScalar(0.0);
1259 CrsWrapper_Epetra_CrsMatrix ecrsmat(C);
1262 case 1: EPETRA_CHK_ERR( mult_A_B(A,Aview,B,Bview,C,call_FillComplete_on_result, keep_all_hard_zeros) );
1264 case 2: EPETRA_CHK_ERR(
mult_A_Btrans(Aview, Bview, ecrsmat, keep_all_hard_zeros) );
1266 case 3: EPETRA_CHK_ERR(
mult_Atrans_B(Aview, Bview, ecrsmat) );
1268 case 4: EPETRA_CHK_ERR(
mult_Atrans_Btrans(Aview, Bview, ecrsmat, keep_all_hard_zeros) );
1270 case 5: EPETRA_CHK_ERR( mult_AT_B_newmatrix(Atransview, Bview, C, keep_all_hard_zeros) );
1275 if (scenario != 1 && call_FillComplete_on_result && scenario != 5) {
1283 const Epetra_Map* domainmap =
1284 transposeB ? &(B.RangeMap()) : &(B.DomainMap());
1286 const Epetra_Map* rangemap =
1287 transposeA ? &(A.DomainMap()) : &(A.RangeMap());
1290 EPETRA_CHK_ERR( C.FillComplete(*domainmap, *rangemap) );
1297 delete mapunion1; mapunion1 = NULL;
1298 delete workmap1; workmap1 = NULL;
1299 delete workmap2; workmap2 = NULL;
1306 const Epetra_CrsMatrix& B,
1308 Epetra_CrsMatrix& C,
1309 bool call_FillComplete_on_result,
1310 bool keep_all_hard_zeros)
1312 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1313 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1314 return TMultiply<int>(A, transposeA, B, transposeB, C, call_FillComplete_on_result, keep_all_hard_zeros);
1318 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1319 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1320 return TMultiply<long long>(A, transposeA, B, transposeB, C, call_FillComplete_on_result, keep_all_hard_zeros);
1324 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1328 template<
typename int_type>
1329 int MatrixMatrix::TAdd(
const Epetra_CrsMatrix& A,
1332 Epetra_CrsMatrix& B,
1344 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;
1349 Epetra_CrsMatrix * Aprime = 0;
1354 Aprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Atrans)(
const_cast<Epetra_CrsMatrix&
>(A)))));
1357 Aprime =
const_cast<Epetra_CrsMatrix*
>(&A);
1359 int MaxNumEntries = EPETRA_MAX( A.MaxNumEntries(), B.MaxNumEntries() );
1360 int A_NumEntries, B_NumEntries;
1361 int_type * A_Indices =
new int_type[MaxNumEntries];
1362 double * A_Values =
new double[MaxNumEntries];
1363 int* B_Indices_local;
1364 int_type* B_Indices_global;
1367 int NumMyRows = B.NumMyRows();
1375 for(
int i = 0; i < NumMyRows; ++i )
1377 Row = (int_type) B.GRID64(i);
1378 EPETRA_CHK_ERR( Aprime->ExtractGlobalRowCopy( Row, MaxNumEntries, A_NumEntries, A_Values, A_Indices ) );
1380 if (scalarB != 1.0) {
1382 EPETRA_CHK_ERR( B.ExtractGlobalRowView( Row, B_NumEntries,
1383 B_Values, B_Indices_global));
1386 EPETRA_CHK_ERR( B.ExtractMyRowView( i, B_NumEntries,
1387 B_Values, B_Indices_local));
1390 for(
int jj=0; jj<B_NumEntries; ++jj) {
1391 B_Values[jj] = scalarB*B_Values[jj];
1395 if( scalarA != 1.0 ) {
1396 for(
int j = 0; j < A_NumEntries; ++j ) A_Values[j] *= scalarA;
1400 err = B.SumIntoGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1402 if (err < 0) ierr = err;
1405 err = B.InsertGlobalValues( Row, A_NumEntries, A_Values, A_Indices );
1406 assert( err == 0 || err == 1 || err == 3 );
1407 if (err < 0) ierr = err;
1412 EPETRA_CHK_ERR( B.Scale(scalarB) );
1415 delete [] A_Indices;
1418 if( Atrans )
delete Atrans;
1426 Epetra_CrsMatrix& B,
1429 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1430 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1431 return TAdd<int>(A, transposeA, scalarA, B, scalarB);
1435 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1436 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1437 return TAdd<long long>(A, transposeA, scalarA, B, scalarB);
1441 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1444 template<
typename int_type>
1445 int MatrixMatrix::TAdd(
const Epetra_CrsMatrix& A,
1448 const Epetra_CrsMatrix & B,
1451 Epetra_CrsMatrix * & C)
1458 if (!A.Filled() || !B.Filled() ) {
1459 std::cerr <<
"EpetraExt::MatrixMatrix::Add ERROR, input matrix A.Filled() or B.Filled() is false," 1460 <<
"they are required to be true. (Result matrix C should be an empty pointer)" << std::endl;
1464 Epetra_CrsMatrix * Aprime = 0, * Bprime=0;
1470 Aprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Atrans)(
const_cast<Epetra_CrsMatrix&
>(A)))));
1473 Aprime =
const_cast<Epetra_CrsMatrix*
>(&A);
1478 Bprime = &(
dynamic_cast<Epetra_CrsMatrix&
>(((*Btrans)(
const_cast<Epetra_CrsMatrix&
>(B)))));
1481 Bprime =
const_cast<Epetra_CrsMatrix*
>(&B);
1487 C =
new Epetra_CrsMatrix(Copy,Aprime->RowMap(),0);
1491 Epetra_CrsMatrix * Mat[] = { Aprime,Bprime};
1492 double scalar[] = { scalarA, scalarB};
1495 for(
int k=0;k<2;k++) {
1496 int MaxNumEntries = Mat[k]->MaxNumEntries();
1498 int_type * Indices =
new int_type[MaxNumEntries];
1499 double * Values =
new double[MaxNumEntries];
1501 int NumMyRows = Mat[k]->NumMyRows();
1506 for(
int i = 0; i < NumMyRows; ++i ) {
1507 Row = (int_type) Mat[k]->GRID64(i);
1508 EPETRA_CHK_ERR( Mat[k]->ExtractGlobalRowCopy( Row, MaxNumEntries, NumEntries, Values, Indices));
1510 if( scalar[k] != 1.0 )
1511 for(
int j = 0; j < NumEntries; ++j ) Values[j] *= scalar[k];
1514 err = C->SumIntoGlobalValues( Row, NumEntries, Values, Indices );
1515 if (err < 0) ierr = err;
1517 err = C->InsertGlobalValues( Row, NumEntries, Values, Indices );
1518 if (err < 0) ierr = err;
1526 if( Atrans )
delete Atrans;
1527 if( Btrans )
delete Btrans;
1535 const Epetra_CrsMatrix & B,
1538 Epetra_CrsMatrix * & C)
1540 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1541 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1542 return TAdd<int>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1546 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1547 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1548 return TAdd<long long>(A, transposeA, scalarA, B, transposeB, scalarB, C);
1552 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Add: GlobalIndices type unknown");
1558 template<
typename int_type>
1559 int MatrixMatrix::TJacobi(
double omega,
1560 const Epetra_Vector & Dinv,
1561 const Epetra_CrsMatrix& A,
1562 const Epetra_CrsMatrix& B,
1563 Epetra_CrsMatrix& C,
1564 bool call_FillComplete_on_result)
1566 #ifdef ENABLE_MMM_TIMINGS 1567 using Teuchos::TimeMonitor;
1568 Teuchos::RCP<TimeMonitor> MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All Setup")));
1572 if (!A.Filled() || !B.Filled()) {
1577 long long Aouter = A.NumGlobalRows64();
1578 long long Bouter = B.NumGlobalCols64();
1579 long long Ainner = A.NumGlobalCols64();
1580 long long Binner = B.NumGlobalRows64();
1581 long long Dlen = Dinv.GlobalLength64();
1582 if (Ainner != Binner) {
1583 std::cerr <<
"MatrixMatrix::Jacobi: ERROR, inner dimensions of A and B " 1584 <<
"must match for matrix-matrix product. A is " 1585 <<Aouter<<
"x"<<Ainner <<
", B is "<<Binner<<
"x"<<Bouter<<std::endl;
1593 if (Aouter > C.NumGlobalRows64()) {
1594 std::cerr <<
"MatrixMatrix::Jacobi: ERROR, dimensions of result C must " 1595 <<
"match dimensions of A * B. C has "<<C.NumGlobalRows64()
1596 <<
" rows, should have at least "<<Aouter << std::endl;
1601 if(Dlen != Aouter) {
1602 std::cerr <<
"MatrixMatrix::Jacboi: ERROR, dimensions of result D must " 1603 <<
"match dimensions of A's rows. D has "<< Dlen
1604 <<
" rows, should have " << Aouter << std::endl;
1608 if(!A.RowMap().SameAs(B.RowMap()) || !A.RowMap().SameAs(Dinv.Map())) {
1609 std::cerr <<
"MatrixMatrix::Jacboi: ERROR, RowMap of A must match RowMap of B " 1610 <<
"and Map of D."<<std::endl;
1621 int numProcs = A.Comm().NumProc();
1624 const Epetra_Map* rowmap_A = &(A.RowMap());
1625 const Epetra_Map* rowmap_B = &(B.RowMap());
1631 const Epetra_Map* workmap1 = NULL;
1632 const Epetra_Map* workmap2 = NULL;
1633 const Epetra_Map* mapunion1 = NULL;
1637 CrsMatrixStruct Aview;
1638 CrsMatrixStruct Bview;
1640 const Epetra_Map* targetMap_A = rowmap_A;
1641 const Epetra_Map* targetMap_B = rowmap_B;
1643 #ifdef ENABLE_MMM_TIMINGS 1644 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All I&X")));
1648 if(call_FillComplete_on_result) {
1649 EPETRA_CHK_ERR(import_only<int_type>(A,*targetMap_A,Aview));
1652 EPETRA_CHK_ERR( import_and_extract_views<int_type>(A, *targetMap_A, Aview));
1659 const Epetra_Map* colmap_op_A = NULL;
1661 colmap_op_A = &(A.ColMap());
1662 targetMap_B = colmap_op_A;
1666 if(call_FillComplete_on_result) {
1667 EPETRA_CHK_ERR(import_only<int_type>(B,*targetMap_B,Bview,A.Importer()));
1670 EPETRA_CHK_ERR( import_and_extract_views<int_type>(B, *targetMap_B, Bview) );
1673 #ifdef ENABLE_MMM_TIMINGS 1674 MM = Teuchos::rcp(
new TimeMonitor(*TimeMonitor::getNewTimer(
"EpetraExt: Jacobi All Multiply")));
1678 if(C.Filled()) C.PutScalar(0.0);
1681 CrsWrapper_Epetra_CrsMatrix ecrsmat(C);
1682 EPETRA_CHK_ERR( jacobi_A_B(omega,Dinv,A,Aview,B,Bview,C,call_FillComplete_on_result) );
1686 delete mapunion1; mapunion1 = NULL;
1687 delete workmap1; workmap1 = NULL;
1688 delete workmap2; workmap2 = NULL;
1696 const Epetra_Vector & Dinv,
1697 const Epetra_CrsMatrix& A,
1698 const Epetra_CrsMatrix& B,
1699 Epetra_CrsMatrix& C,
1700 bool call_FillComplete_on_result)
1702 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1703 if(A.RowMap().GlobalIndicesInt() && B.RowMap().GlobalIndicesInt()) {
1704 return TJacobi<int>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1708 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1709 if(A.RowMap().GlobalIndicesLongLong() && B.RowMap().GlobalIndicesLongLong()) {
1710 return TJacobi<long long>(omega, Dinv, A, B, C, call_FillComplete_on_result);
1714 throw std::runtime_error(
"EpetraExt::MatrixMatrix::Jacobi: GlobalIndices type unknown");
LightweightCrsMatrix * importMatrix
Epetra_BlockMap * RowMapEP_
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double *Values, int *Indices)=0
LightweightMap * RowMapLW_
long long IndexBase64() const
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)
void MyGlobalElementsPtr(int *&MyGlobalElementList) const
const Epetra_Map * rowMap
static int Multiply(const Epetra_CrsMatrix &A, bool transposeA, const Epetra_CrsMatrix &B, bool transposeB, Epetra_CrsMatrix &C, bool call_FillComplete_on_result=true, bool keep_all_hard_zeros=false)
Given Epetra_CrsMatrix objects A, B and C, form the product C = A*B.
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_Atrans_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
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)
virtual const Epetra_Map & RowMap() const =0
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.
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...
Epetra_Map * Tfind_rows_containing_cols(const Epetra_CrsMatrix &M, const Epetra_Map &column_map)
std::vector< int > rowptr_
int mult_A_Btrans(CrsMatrixStruct &Aview, CrsMatrixStruct &Bview, CrsWrapper &C, bool keep_all_hard_zeros)
const Epetra_Map * colMap
const Epetra_BlockMap * importColMap
int NumMyElements() const
int import_only(const Epetra_CrsMatrix &M, const Epetra_Map &targetMap, CrsMatrixStruct &Mview, const Epetra_Import *prototypeImporter)