67 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 69 bool Epetra_CrsMatrixTraceDumpMultiply =
false;
70 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 72 #ifdef HAVE_EPETRA_TEUCHOS 74 # define EPETRA_CRSMATRIX_TEUCHOS_TIMERS 77 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 81 #if defined(Epetra_ENABLE_MKL_SPARSE) && defined(Epetra_ENABLE_CASK) 82 #error Error: Epetra_ENABLE_MKL_SPARSE and Epetra_ENABLE_CASK both defined. 85 #ifdef Epetra_ENABLE_MKL_SPARSE 86 #include "mkl_spblas.h" 94 Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
98 constructedWithFilledGraph_(false),
99 matrixFillCompleteCalled_(false),
100 StorageOptimized_(false),
102 Values_alloc_lengths_(0),
107 NumMyRows_(rowMap.NumMyPoints()),
111 squareFillCompleteCalled_(false)
122 Graph_(CV, rowMap, NumEntriesPerRow, StaticProfile),
125 UseTranspose_(false),
126 constructedWithFilledGraph_(false),
127 matrixFillCompleteCalled_(false),
128 StorageOptimized_(false),
130 Values_alloc_lengths_(0),
135 NumMyRows_(rowMap.NumMyPoints()),
139 squareFillCompleteCalled_(false)
146 const Epetra_Map& colMap,
const int* NumEntriesPerRow,
bool StaticProfile)
150 Graph_(CV, rowMap, colMap, NumEntriesPerRow, StaticProfile),
153 UseTranspose_(false),
154 constructedWithFilledGraph_(false),
155 matrixFillCompleteCalled_(false),
156 StorageOptimized_(false),
158 Values_alloc_lengths_(0),
163 NumMyRows_(rowMap.NumMyPoints()),
167 squareFillCompleteCalled_(false)
175 const Epetra_Map& colMap,
int NumEntriesPerRow,
bool StaticProfile)
179 Graph_(CV, rowMap, colMap, NumEntriesPerRow, StaticProfile),
182 UseTranspose_(false),
183 constructedWithFilledGraph_(false),
184 matrixFillCompleteCalled_(false),
185 StorageOptimized_(false),
187 Values_alloc_lengths_(0),
192 NumMyRows_(rowMap.NumMyPoints()),
196 squareFillCompleteCalled_(false)
199 throw ReportError(
"Epetra_CrsMatrix::Epetra_CrsMatrix: cannot be called with different indices types for rowMap and colMap", -1);
212 UseTranspose_(false),
213 constructedWithFilledGraph_(false),
214 matrixFillCompleteCalled_(false),
215 StorageOptimized_(false),
217 Values_alloc_lengths_(0),
222 NumMyRows_(graph.NumMyRows()),
226 squareFillCompleteCalled_(false)
238 Graph_(Matrix.Graph()),
241 UseTranspose_(Matrix.UseTranspose_),
242 constructedWithFilledGraph_(false),
243 matrixFillCompleteCalled_(false),
244 StorageOptimized_(false),
246 Values_alloc_lengths_(0),
251 NumMyRows_(Matrix.NumMyRows()),
255 squareFillCompleteCalled_(false)
268 if (!src.
Filled())
throw ReportError(
"Copying an Epetra_CrsMatrix requires source matrix to have Filled()==true", -1);
294 if (numMyNonzeros>0)
All_Values_ =
new double[numMyNonzeros];
296 for (
int i=0; i<numMyNonzeros; ++i)
All_Values_[i] = srcValues[i];
298 #ifdef Epetra_ENABLE_CASK 300 cask = cask_handler_copy(src.cask);
309 double *
const srcValues = src.
Values(i);
310 double * targValues =
Values(i);
311 for (
int j=0; j< NumEntries; j++) targValues[j] = srcValues[j];
350 if (numMyNonzeros>0)
All_Values_ =
new double[numMyNonzeros];
359 if (NumAllocatedEntries > 0) {
362 all_values += NumAllocatedEntries;
365 Values_[i] =
new double[NumAllocatedEntries];
372 for(j=0; j< NumAllocatedEntries; j++)
382 #ifdef Epetra_ENABLE_CASK 403 if (
Graph().NumAllocatedMyIndices(i) >0)
423 #ifdef Epetra_ENABLE_CASK 427 cask_handler_destroy(cask);
491 for (
int i=0; i<length; ++i)
All_Values_[i] = ScalarConstant;
496 double * targValues =
Values(i);
497 for(
int j=0; j< NumEntries; j++)
498 targValues[j] = ScalarConstant;
508 for (
int i=0; i<length; ++i)
All_Values_[i] *= ScalarConstant;
513 double * targValues =
Values(i);
514 for(
int j=0; j< NumEntries; j++)
515 targValues[j] *= ScalarConstant;
522 template<
typename int_type>
524 const double* values,
525 const int_type* Indices)
539 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 541 const double* values,
544 if(
RowMap().GlobalIndicesInt())
545 return TInsertGlobalValues<int>(Row, NumEntries, values, Indices);
547 throw ReportError(
"Epetra_CrsMatrix::InsertGlobalValues int version called for a matrix that is not int.", -1);
550 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 552 const double* values,
553 const long long* Indices)
555 if(
RowMap().GlobalIndicesLongLong())
556 return TInsertGlobalValues<long long>(Row, NumEntries, values, Indices);
558 throw ReportError(
"Epetra_CrsMatrix::InsertGlobalValues long long version called for a matrix that is not long long.", -1);
563 template<
typename int_type>
580 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 585 if(
RowMap().GlobalIndicesInt())
586 return TInsertGlobalValues<int>(Row, NumEntries, values, Indices);
588 throw ReportError(
"Epetra_CrsMatrix::InsertGlobalValues int version called for a matrix that is not int.", -1);
591 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 596 if(
RowMap().GlobalIndicesLongLong()) {
597 return TInsertGlobalValues<long long>(Row, NumEntries, values, Indices);
600 throw ReportError(
"Epetra_CrsMatrix::InsertGlobalValues long long version called for a matrix that is not long long.", -1);
605 const double* values,
614 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 617 throw ReportError(
"Epetra_CrsMatrix::InsertMyValues: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
635 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 638 throw ReportError(
"Epetra_CrsMatrix::InsertMyValues: Failure because neither 32 bit nor 64 bit indices insertable.", -1);
646 template<
typename int_type>
648 const double* values,
649 const int_type* Indices)
659 return(
InsertValues(Row, NumEntries, const_cast<double*>(values), const_cast<int_type*>(Indices)));
664 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 666 const double* values,
670 return InsertValues<int>(Row, NumEntries, values, Indices);
672 throw ReportError(
"Epetra_CrsMatrix::InsertValues int version called for a matrix that is not int.", -1);
675 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 677 const double* values,
678 const long long* Indices)
680 if(
RowMap().GlobalIndicesLongLong())
681 return InsertValues<long long>(Row, NumEntries, values, Indices);
683 throw ReportError(
"Epetra_CrsMatrix::InsertValues long long version called for a matrix that is not long long.", -1);
687 template<
typename int_type>
709 if(NumEntries != testNumEntries)
711 for(
int i = 0; i < NumEntries; ++i)
712 match = match && (Indices[i]==testIndices[i]);
726 int tmpNumEntries = NumEntries;
729 const double* tmpValues = values;
730 values =
new double[NumEntries];
733 for(
int i = 0; i < NumEntries; ++i)
735 values[loc++] = tmpValues[i];
738 for(
int i = 0; i < NumEntries; ++i)
740 values[loc++] = tmpValues[i];
742 if(NumEntries != loc)
748 int stop = start + NumEntries;
750 if(stop > NumAllocatedEntries) {
751 if (
Graph().StaticProfile() && stop >
Graph().NumAllocatedMyIndices(Row)) {
754 if(NumAllocatedEntries == 0) {
755 Values_[Row] =
new double[NumEntries];
760 double* tmp_Values =
new double[stop];
761 for(j = 0; j < start; j++)
762 tmp_Values[j] =
Values_[Row][j];
769 for(j = start; j < stop; j++)
770 Values_[Row][j] = values[j-start];
773 NumEntries = tmpNumEntries;
790 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 796 return InsertValues<int>(Row, NumEntries, values, Indices);
798 throw ReportError(
"Epetra_CrsMatrix::InsertValues int version called for a matrix that is not int.", -1);
801 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 806 if(
RowMap().GlobalIndicesLongLong())
807 return InsertValues<long long>(Row, NumEntries, values, Indices);
809 throw ReportError(
"Epetra_CrsMatrix::InsertValues long long version called for a matrix that is not long long.", -1);
822 template<
typename int_type>
834 double * targValues =
Values(locRow);
835 for (j=0; j<NumEntries; j++) {
836 int_type Index = Indices[j];
838 targValues[Loc] = srcValues[j];
851 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 853 if(
RowMap().GlobalIndicesInt())
854 return TReplaceGlobalValues<int>(Row, NumEntries, srcValues, Indices);
856 throw ReportError(
"Epetra_CrsMatrix::ReplaceGlobalValues int version called for a matrix that is not int.", -1);
859 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 861 if(
RowMap().GlobalIndicesLongLong())
862 return TReplaceGlobalValues<long long>(Row, NumEntries, srcValues, Indices);
864 throw ReportError(
"Epetra_CrsMatrix::ReplaceGlobalValues long long version called for a matrix that is not long long.", -1);
882 double* RowValues =
Values(Row);
883 for (j=0; j<NumEntries; j++) {
884 int Index = Indices[j];
886 RowValues[Loc] = srcValues[j];
901 const double * srcValues,
const int *Offsets)
907 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES 908 int locRow =
LRID((
int) Row);
910 int locRow =
LRID(Row);
917 double* RowValues =
Values(locRow);
918 for(j=0; j<NumEntries; j++) {
919 if( Offsets[j] != -1 )
920 RowValues[Offsets[j]] = srcValues[j];
932 template<
typename int_type>
935 const double * srcValues,
936 const int_type *Indices)
953 double * RowValues =
Values(locRow);
956 for (j=0; j<NumEntries; j++) {
957 int_type Index = Indices[j];
959 #ifdef EPETRA_HAVE_OMP
960 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
964 RowValues[Loc] += srcValues[j];
976 for (j=0; j<NumEntries; j++) {
977 int Index = colmap.
LID(Indices[j]);
981 if (Loc < NumColIndices && Index == ColIndices[Loc])
982 #ifdef EPETRA_HAVE_OMP 983 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 987 RowValues[Loc] += srcValues[j];
991 #ifdef EPETRA_HAVE_OMP 992 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE 996 RowValues[Loc] += srcValues[j];
1004 for (j=0; j<NumEntries; j++) {
1005 int Index = colmap.
LID(Indices[j]);
1007 #ifdef EPETRA_HAVE_OMP
1008 #ifdef EPETRA_HAVE_OMP_NONASSOCIATIVE
1012 RowValues[Loc] += srcValues[j];
1027 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1030 const double * srcValues,
1033 if(
RowMap().GlobalIndicesInt())
1034 return TSumIntoGlobalValues<int>(Row, NumEntries, srcValues, Indices);
1036 throw ReportError(
"Epetra_CrsMatrix::SumIntoGlobalValues int version called for a matrix that is not int.", -1);
1039 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1042 const double * srcValues,
1043 const long long *Indices)
1045 if(
RowMap().GlobalIndicesLongLong())
1046 return TSumIntoGlobalValues<long long>(Row, NumEntries, srcValues, Indices);
1048 throw ReportError(
"Epetra_CrsMatrix::SumIntoGlobalValues long long version called for a matrix that is not long long.", -1);
1067 double* RowValues =
Values(Row);
1071 for (j=0; j<NumEntries; j++) {
1072 int Index = Indices[j];
1076 if (Loc < NumColIndices && Index == ColIndices[Loc])
1077 RowValues[Loc] += srcValues[j];
1081 RowValues[Loc] += srcValues[j];
1089 for (j=0; j<NumEntries; j++) {
1090 int Index = Indices[j];
1092 RowValues[Loc] += srcValues[j];
1114 #ifdef EPETRA_NO_64BIT_GLOBAL_INDICES 1115 int locRow =
LRID((
int) Row);
1117 int locRow =
LRID(Row);
1124 double* RowValues =
Values(locRow);
1125 for (j=0; j<NumEntries; j++) {
1126 if( Offsets[j] != -1 )
1127 RowValues[Offsets[j]] += srcValues[j];
1148 const Epetra_Map& range_map,
bool OptimizeDataStorage)
1150 int returnValue = 0;
1183 return(returnValue);
1212 double* locValues =
Values(i);
1221 for(
int j = 0; j < max; j++) {
1222 for(
int k = j; k >= 0; k-=m) {
1223 if(locIndices[k+m] >= locIndices[k])
1225 double dtemp = locValues[k+m];
1226 locValues[k+m] = locValues[k];
1227 locValues[k] = dtemp;
1228 int itemp = locIndices[k+m];
1229 locIndices[k+m] = locIndices[k];
1230 locIndices[k] = itemp;
1256 if(NumEntries > 1) {
1257 double*
const locValues =
Values(i);
1260 double curValue = locValues[0];
1261 for(
int k = 1; k < NumEntries; k++) {
1262 if(locIndices[k] == locIndices[k-1])
1263 curValue += locValues[k];
1265 locValues[curEntry++] = curValue;
1266 curValue = locValues[k];
1269 locValues[curEntry] = curValue;
1290 bool Contiguous =
true;
1305 if ((
CV_==
View) && !Contiguous)
1323 #ifdef EPETRA_HAVE_OMP 1324 #pragma omp parallel for default(none) shared(Values_s,All_Values_s) 1326 for (
int i=0; i<numMyRows; i++) {
1328 int curOffset = IndexOffset[i];
1329 double * values = Values_s[i];
1330 double * newValues = All_Values_s+curOffset;
1331 for (
int j=0; j<NumEntries; j++) newValues[j] = values[j];
1340 for (
int j=0; j<NumEntries; j++) tmp[j] = values[j];
1362 #ifdef Epetra_ENABLE_CASK 1367 cask_handler_initialize(&cask);
1368 cask_csr_analysis(
NumMyRows_, NumMyCols_, IndexOffset, Indices,
1381 template<
typename int_type>
1383 int_type * Indices)
const 1394 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1396 int * Indices)
const 1398 if(
RowMap().GlobalIndicesInt())
1399 return ExtractGlobalRowCopy<int>(Row, Length, NumEntries, values, Indices);
1401 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowCopy int version called for a matrix that is not int.", -1);
1404 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1406 long long * Indices)
const 1408 if(
RowMap().GlobalIndicesLongLong())
1409 return ExtractGlobalRowCopy<long long>(Row, Length, NumEntries, values, Indices);
1411 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowCopy long long version called for a matrix that is not long long.", -1);
1417 int * Indices)
const 1438 template<
typename int_type>
1447 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1450 if(
RowMap().GlobalIndicesInt())
1451 return ExtractGlobalRowCopy<int>(Row, Length, NumEntries, values);
1453 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowCopy int version called for a matrix that is not int.", -1);
1456 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1459 if(
RowMap().GlobalIndicesLongLong())
1460 return ExtractGlobalRowCopy<long long>(Row, Length, NumEntries, values);
1462 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowCopy long long version called for a matrix that is not long long.", -1);
1475 if (Length < NumEntries)
1478 double * srcValues =
Values(Row);
1480 for(j=0; j<NumEntries; j++)
1481 targValues[j] = srcValues[j];
1495 long long ii =
GRID64(i);
1498 double * srcValues =
Values(i);
1501 for(
int j = 0; j < NumEntries; j++) {
1502 if(ii ==
GCID64(Indices[j])) {
1503 Diagonal[i] = srcValues[j];
1520 long long ii =
GRID64(i);
1523 double * targValues =
Values(i);
1524 bool DiagMissing =
true;
1525 for(
int j = 0; j < NumEntries; j++) {
1526 if(ii ==
GCID64(Indices[j])) {
1527 targValues[j] = Diagonal[i];
1528 DiagMissing =
false;
1546 template<
typename int_type>
1557 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1560 if(
RowMap().GlobalIndicesInt())
1561 return ExtractGlobalRowView<int>(Row, NumEntries, values, Indices);
1563 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowView int version called for a matrix that is not int.", -1);
1566 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1569 if(
RowMap().GlobalIndicesLongLong())
1570 return ExtractGlobalRowView<long long>(Row, NumEntries, values, Indices);
1572 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowView long long version called for a matrix that is not long long.", -1);
1587 template<
typename int_type>
1596 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 1599 if(
RowMap().GlobalIndicesInt())
1600 return ExtractGlobalRowView<int>(Row, NumEntries, values);
1602 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowView int version called for a matrix that is not int.", -1);
1605 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 1608 if(
RowMap().GlobalIndicesLongLong())
1609 return ExtractGlobalRowView<long long>(Row, NumEntries, values);
1611 throw ReportError(
"Epetra_CrsMatrix::ExtractGlobalRowView long long version called for a matrix that is not long long.", -1);
1623 targValues =
Values(Row);
1633 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 1660 double *xp = (
double*)x.
Values();
1661 double *yp = (
double*)y.
Values();
1664 GeneralSV(Upper, Trans, UnitDiagonal, xp, yp);
1673 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 1698 double** Xp = (
double**) X.
Pointers();
1699 double** Yp = (
double**) Y.
Pointers();
1707 GeneralSV(Upper, Trans, UnitDiagonal, *Xp, *Yp);
1709 GeneralSM(Upper, Trans, UnitDiagonal, Xp, LDX, Yp, LDY, NumVectors);
1725 double * xp = (
double*)x.
Values();
1729 double * x_tmp_p = (
double*)x_tmp.
Values();
1732 double * RowValues =
Values(i);
1733 for (j=0; j < NumEntries; j++) x_tmp_p[i] += std::abs(RowValues[j]);
1737 for (i=0; i<myLength; i++) {
1739 if (xp[i]==0.0) ierr = 1;
1740 else if (ierr!=1) ierr = 2;
1750 double * RowValues =
Values(i);
1752 for (j=0; j < NumEntries; j++) scale += std::abs(RowValues[j]);
1754 if (scale==0.0) ierr = 1;
1755 else if (ierr!=1) ierr = 2;
1779 bool needExport =
false;
1780 double * xp = (
double*)x.
Values();
1786 xp = (
double*)x_tmp->
Values();
1794 double * RowValues =
Values(i);
1796 for (j=0; j < NumEntries; j++) scale =
EPETRA_MAX(std::abs(RowValues[j]),scale);
1798 if (scale==0.0) ierr = 1;
1799 else if (ierr!=1) ierr = 2;
1826 double* xp = (
double*)x.
Values();
1830 double * x_tmp_p = (
double*)x_tmp.
Values();
1834 double* RowValues =
Values(i);
1835 for(j = 0; j < NumEntries; j++)
1836 x_tmp_p[ColIndices[j]] += std::abs(RowValues[j]);
1844 double* RowValues =
Values(i);
1845 for(j = 0; j < NumEntries; j++)
1846 xp[ColIndices[j]] += std::abs(RowValues[j]);
1854 for(i = 0; i < MapNumMyElements; i++) {
1855 double scale = xp[i];
1864 xp[i] = 1.0 / scale;
1883 double* xp = (
double*)x.
Values();
1887 double * x_tmp_p = (
double*)x_tmp.
Values();
1891 double* RowValues =
Values(i);
1892 for(j = 0; j < NumEntries; j++)
1893 x_tmp_p[ColIndices[j]] =
EPETRA_MAX(std::abs(RowValues[j]),x_tmp_p[ColIndices[j]]);
1901 double* RowValues =
Values(i);
1902 for(j = 0; j < NumEntries; j++)
1903 xp[ColIndices[j]] =
EPETRA_MAX(std::abs(RowValues[j]),xp[ColIndices[j]]);
1911 for(i = 0; i < MapNumMyElements; i++) {
1912 double scale = xp[i];
1921 xp[i] = 1.0 / scale;
1948 xp = (
double*)x.
Values();
1950 xp = (
double*)x.
Values();
1958 double* RowValues =
Values(i);
1959 double scale = xp[i];
1960 for(j = 0; j < NumEntries; j++)
1961 RowValues[j] *= scale;
1990 xp = (
double*)x.
Values();
1992 xp = (
double*)x.
Values();
2000 double* RowValues =
Values(i);
2001 for(j = 0; j < NumEntries; j++)
2002 RowValues[j] *= xp[ColIndices[j]];
2032 double* xp = (
double*)x.
Values();
2038 xp = (
double*)x_tmp->
Values();
2045 double* RowValues =
Values(i);
2046 for(j = 0; j < NumEntries; j++)
2047 xp[i] += std::abs(RowValues[j]);
2080 double* xp = (
double*)x.
Values();
2088 xp = (
double*)x_tmp->
Values();
2092 for(i = 0; i < NumCols; i++)
2098 double* RowValues =
Values(i);
2099 for(j = 0; j < NumEntries; j++)
2100 xp[ColIndices[j]] += std::abs(RowValues[j]);
2131 double local_sum = 0.0;
2135 double* RowValues =
Values(i);
2136 for(
int j = 0; j < NumEntries; j++) {
2137 local_sum += RowValues[j]*RowValues[j];
2141 double global_sum = 0.0;
2166 int * PermuteToLIDs,
2167 int *PermuteFromLIDs,
2172 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermute: Incoming global index type does not match the one for *this",-1);
2177 PermuteFromLIDs,Indexor,CombineMode));
2183 PermuteFromLIDs,Indexor,CombineMode));
2194 template<
typename int_type>
2198 int * PermuteToLIDs,
2199 int *PermuteFromLIDs,
2207 int maxNumEntries =
A.MaxNumEntries();
2208 int_type * Indices = 0;
2209 double * values = 0;
2211 if (maxNumEntries>0 &&
A.IndicesAreLocal() ) {
2212 Indices =
new int_type[maxNumEntries];
2213 values =
new double[maxNumEntries];
2218 if (
A.IndicesAreLocal()) {
2221 for (i=0; i<NumSameIDs; i++) {
2222 Row = (int_type)
GRID64(i);
2223 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices));
2232 for (i=0; i<NumSameIDs; i++) {
2233 Row = (int_type)
GRID64(i);
2234 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices));
2245 for (i=0; i<NumSameIDs; i++) {
2246 Row = (int_type)
GRID64(i);
2247 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices));
2253 for (i=0; i<NumSameIDs; i++) {
2254 Row = (int_type)
GRID64(i);
2255 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(Row, maxNumEntries, NumEntries, values, Indices));
2265 for (i=0; i<NumSameIDs; i++) {
2266 Row = (int_type)
GRID64(i);
2267 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(Row, NumEntries, values, Indices));
2276 for (i=0; i<NumSameIDs; i++) {
2277 Row = (int_type)
GRID64(i);
2278 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(Row, NumEntries, values, Indices));
2289 for (i=0; i<NumSameIDs; i++) {
2290 Row = (int_type)
GRID64(i);
2291 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(Row, NumEntries, values, Indices));
2297 for (i=0; i<NumSameIDs; i++) {
2298 Row = (int_type)
GRID64(i);
2299 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(Row, NumEntries, values, Indices));
2309 int_type FromRow, ToRow;
2310 if (NumPermuteIDs>0) {
2311 if (
A.IndicesAreLocal()) {
2314 for (i=0; i<NumPermuteIDs; i++) {
2315 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2316 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2317 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
2326 for (i=0; i<NumPermuteIDs; i++) {
2327 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2328 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2329 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
2340 for (i=0; i<NumPermuteIDs; i++) {
2341 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2342 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2343 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
2349 for (i=0; i<NumPermuteIDs; i++) {
2350 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2351 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2352 EPETRA_CHK_ERR(
A.ExtractGlobalRowCopy(FromRow, maxNumEntries, NumEntries, values, Indices));
2362 for (i=0; i<NumPermuteIDs; i++) {
2363 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2364 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2365 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices));
2374 for (i=0; i<NumPermuteIDs; i++) {
2375 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2376 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2377 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices));
2388 for (i=0; i<NumPermuteIDs; i++) {
2389 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2390 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2391 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices));
2397 for (i=0; i<NumPermuteIDs; i++) {
2398 FromRow = (int_type)
A.GRID64(PermuteFromLIDs[i]);
2399 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2400 EPETRA_CHK_ERR(
A.ExtractGlobalRowView(FromRow, NumEntries, values, Indices));
2409 if (maxNumEntries>0 &&
A.IndicesAreLocal() ) {
2420 int * PermuteToLIDs,
2421 int *PermuteFromLIDs,
2425 if(!
A.RowMap().GlobalIndicesTypeMatch(
RowMap()))
2426 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: Incoming global index type does not match the one for *this",-1);
2428 if(
A.RowMap().GlobalIndicesInt())
2429 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2430 return TCopyAndPermuteCrsMatrix<int>(
A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2432 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: ERROR, GlobalIndicesInt but no API for it.",-1);
2435 if(
A.RowMap().GlobalIndicesLongLong())
2436 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2437 return TCopyAndPermuteCrsMatrix<long long>(
A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2439 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2442 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: Internal error, unable to determine global index type of maps", -1);
2446 template<
typename int_type>
2450 int * PermuteToLIDs,
2451 int *PermuteFromLIDs,
2456 int_type Row, ToRow;
2459 int maxNumEntries =
A.MaxNumEntries();
2460 int_type * gen_Indices = 0;
2461 int * int_Indices = 0;
2462 double * values = 0;
2464 if (maxNumEntries>0) {
2466 int_Indices =
new int[maxNumEntries];
2469 gen_Indices =
new int_type[maxNumEntries];
2470 int_Indices =
reinterpret_cast<int*
>(gen_Indices);
2473 values =
new double[maxNumEntries];
2483 for (i=0; i<NumSameIDs; i++) {
2484 Row = (int_type)
GRID64(i);
2485 int AlocalRow = rowMap.
LID(Row);
2486 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, int_Indices));
2492 for (i=0; i<NumSameIDs; i++) {
2493 Row = (int_type)
GRID64(i);
2494 int AlocalRow = rowMap.
LID(Row);
2495 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(AlocalRow, maxNumEntries, NumEntries, values, int_Indices));
2496 for(j=0; j<NumEntries; ++j) {
2497 int_Indices[j] =
LCID((int_type) colMap.
GID64(int_Indices[j]));
2506 for (i=0; i<NumSameIDs; i++) {
2507 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, int_Indices));
2508 Row = (int_type)
GRID64(i);
2514 for (i=0; i<NumSameIDs; i++) {
2515 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(i, maxNumEntries, NumEntries, values, int_Indices));
2516 Row = (int_type)
GRID64(i);
2519 for(j = NumEntries; j > 0;) {
2521 gen_Indices[j] = (int_type) colMap.
GID64(int_Indices[j]);
2531 if (NumPermuteIDs>0) {
2534 for (i=0; i<NumPermuteIDs; i++) {
2535 FromRow = PermuteFromLIDs[i];
2536 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
2537 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2543 for (i=0; i<NumPermuteIDs; i++) {
2544 FromRow = PermuteFromLIDs[i];
2545 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
2546 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2547 for(j=0; j<NumEntries; ++j) {
2548 int_Indices[j] =
LCID((int_type) colMap.
GID64(int_Indices[j]));
2557 for (i=0; i<NumPermuteIDs; i++) {
2558 FromRow = PermuteFromLIDs[i];
2559 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
2560 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2566 for (i=0; i<NumPermuteIDs; i++) {
2567 FromRow = PermuteFromLIDs[i];
2568 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(FromRow, maxNumEntries, NumEntries, values, int_Indices));
2571 for(j = NumEntries; j > 0;) {
2573 gen_Indices[j] = (int_type) colMap.
GID64(int_Indices[j]);
2576 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2584 if (maxNumEntries>0) {
2587 delete [] int_Indices;
2590 delete [] gen_Indices;
2599 int * PermuteToLIDs,
2600 int *PermuteFromLIDs,
2604 if(!
A.Map().GlobalIndicesTypeMatch(
RowMap()))
2605 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
2607 if(
A.RowMatrixRowMap().GlobalIndicesInt())
2608 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2609 return TCopyAndPermuteRowMatrix<int>(
A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2611 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteRowMatrix: ERROR, GlobalIndicesInt but no API for it.",-1);
2614 if(
A.RowMatrixRowMap().GlobalIndicesLongLong())
2615 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2616 return TCopyAndPermuteRowMatrix<long long>(
A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs, Indexor, CombineMode);
2618 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteRowMatrix: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2621 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteRowMatrix: Internal error, unable to determine global index type of maps", -1);
2636 throw ReportError(
"Epetra_CrsMatrix::PackAndPrepare: Incoming global index type does not match the one for *this",-1);
2644 int TotalSendLength = 0;
2646 if( NumExportIDs>0 ) IntSizes =
new int[NumExportIDs];
2648 int SizeofIntType = -1;
2650 SizeofIntType = (
int)
sizeof(int);
2652 SizeofIntType = (
int)
sizeof(
long long);
2654 throw ReportError(
"Epetra_CrsMatrix::PackAndPrepare: Unable to determine source global index type",-1);
2656 for(
int i = 0; i < NumExportIDs; ++i )
2659 A.NumMyRowEntries( ExportLIDs[i], NumEntries );
2661 Sizes[i] = NumEntries;
2662 IntSizes[i] = 1 + (((NumEntries+2)*SizeofIntType)/(int)
sizeof(
double));
2663 TotalSendLength += (Sizes[i]+IntSizes[i]);
2666 double * DoubleExports = 0;
2667 SizeOfPacket = (int)
sizeof(
double);
2670 if( TotalSendLength*SizeOfPacket > LenExports )
2672 if( LenExports > 0 )
delete [] Exports;
2673 LenExports = TotalSendLength*SizeOfPacket;
2674 DoubleExports =
new double[TotalSendLength];
2675 for(
int i = 0; i < TotalSendLength; ++i ) DoubleExports[i] = 0.0;
2676 Exports = (
char *) DoubleExports;
2681 double * valptr, * dintptr;
2691 if( NumExportIDs > 0 )
2698 int maxNumEntries =
A.MaxNumEntries();
2699 dintptr = (
double *) Exports;
2700 valptr = dintptr + IntSizes[0];
2701 intptr = (
int *) dintptr;
2702 for (
int i=0; i<NumExportIDs; i++)
2704 FromRow = (int) rowMap.
GID64(ExportLIDs[i]);
2705 intptr[0] = FromRow;
2707 Indices = intptr + 2;
2708 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, Indices));
2709 for (
int j=0; j<NumEntries; j++) Indices[j] = (
int) colMap.
GID64(Indices[j]);
2710 intptr[1] = NumEntries;
2711 if( i < (NumExportIDs-1) )
2713 dintptr += (IntSizes[i]+Sizes[i]);
2714 valptr = dintptr + IntSizes[i+1];
2715 intptr = (
int *) dintptr;
2720 long long * LL_Indices;
2724 int maxNumEntries =
A.MaxNumEntries();
2725 dintptr = (
double *) Exports;
2726 valptr = dintptr + IntSizes[0];
2727 LLptr = (
long long *) dintptr;
2728 for (
int i=0; i<NumExportIDs; i++)
2730 FromRow = rowMap.
GID64(ExportLIDs[i]);
2733 LL_Indices = LLptr + 2;
2734 int * int_indices =
reinterpret_cast<int*
>(LL_Indices);
2735 EPETRA_CHK_ERR(
A.ExtractMyRowCopy(ExportLIDs[i], maxNumEntries, NumEntries, values, int_indices));
2738 for(
int j = NumEntries; j > 0;) {
2740 LL_Indices[j] = colMap.
GID64(int_indices[j]);
2743 LLptr[1] = NumEntries;
2744 if( i < (NumExportIDs-1) )
2746 dintptr += (IntSizes[i]+Sizes[i]);
2747 valptr = dintptr + IntSizes[i+1];
2748 LLptr = (
long long *) dintptr;
2753 for(
int i = 0; i < NumExportIDs; ++i )
2754 Sizes[i] += IntSizes[i];
2757 if( IntSizes )
delete [] IntSizes;
2763 template<
typename int_type>
2778 if (NumImportIDs<=0)
return(0);
2780 if ( CombineMode !=
Add 2782 && CombineMode !=
Zero )
2792 double * valptr, *dintptr;
2800 dintptr = (
double *) Imports;
2801 intptr = (int_type *) dintptr;
2802 NumEntries = (int) intptr[1];
2803 IntSize = 1 + (((NumEntries+2)*(
int)
sizeof(int_type))/(int)
sizeof(
double));
2804 valptr = dintptr + IntSize;
2806 for (i=0; i<NumImportIDs; i++)
2808 ToRow = (int_type)
GRID64(ImportLIDs[i]);
2809 assert((intptr[0])==ToRow);
2811 Indices = intptr + 2;
2813 if (CombineMode==
Add) {
2828 else if (CombineMode==
Insert) {
2844 if( i < (NumImportIDs-1) )
2846 dintptr += IntSize + NumEntries;
2847 intptr = (int_type *) dintptr;
2848 NumEntries = (int) intptr[1];
2849 IntSize = 1 + (((NumEntries+2)*(
int)
sizeof(int_type))/(int)
sizeof(
double));
2850 valptr = dintptr + IntSize;
2868 throw ReportError(
"Epetra_CrsMatrix::UnpackAndCombine: Incoming global index type does not match the one for *this",-1);
2871 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
2872 return TUnpackAndCombine<int>(Source, NumImportIDs, ImportLIDs, LenImports,
2873 Imports, SizeOfPacket, Distor, CombineMode, Indexor);
2875 throw ReportError(
"Epetra_CrsMatrix::UnpackAndCombine: ERROR, GlobalIndicesInt but no API for it.",-1);
2879 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
2880 return TUnpackAndCombine<long long>(Source, NumImportIDs, ImportLIDs, LenImports,
2881 Imports, SizeOfPacket, Distor, CombineMode, Indexor);
2883 throw ReportError(
"Epetra_CrsMatrix::UnpackAndCombine: ERROR, GlobalIndicesLongLong but no API for it.",-1);
2886 throw ReportError(
"Epetra_CrsMatrix::UnpackAndCombine: Internal error, unable to determine global index type of maps", -1);
2895 for (
int iproc=0; iproc < NumProc; iproc++) {
2901 os <<
"\nNumber of Global Rows = "; os <<
NumGlobalRows64(); os << std::endl;
2902 os <<
"Number of Global Cols = "; os <<
NumGlobalCols64(); os << std::endl;
2906 if (
LowerTriangular()) { os <<
" ** Matrix is Lower Triangular **"; os << std::endl; }
2907 if (
UpperTriangular()) { os <<
" ** Matrix is Upper Triangular **"; os << std::endl; }
2908 if (
NoDiagonal()) { os <<
" ** Matrix has no diagonal **"; os << std::endl; os << std::endl; }
2911 os <<
"\nNumber of My Rows = "; os <<
NumMyRows(); os << std::endl;
2912 os <<
"Number of My Cols = "; os <<
NumMyCols(); os << std::endl;
2913 os <<
"Number of My Diagonals = "; os <<
NumMyDiagonals(); os << std::endl;
2914 os <<
"Number of My Nonzeros = "; os <<
NumMyNonzeros(); os << std::endl;
2915 os <<
"My Maximum Num Entries = "; os <<
MaxNumEntries(); os << std::endl; os << std::endl;
2931 {
for (
int iproc=0; iproc < NumProc; iproc++) {
2936 int * Indices_int = 0;
2937 long long * Indices_LL = 0;
2938 if(
RowMap().GlobalIndicesInt()) {
2939 Indices_int =
new int[MaxNumIndices];
2941 else if(
RowMap().GlobalIndicesLongLong()) {
2942 Indices_LL =
new long long[MaxNumIndices];
2945 throw ReportError(
"Epetra_CrsGraph::Print: Unable to determine source global index type",-1);
2947 double * values =
new double[MaxNumIndices];
2948 #if !defined(EPETRA_NO_32BIT_GLOBAL_INDICES) || !defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 2956 os <<
" Processor ";
2958 os <<
" Row Index ";
2960 os <<
" Col Index ";
2965 for (i=0; i<NumMyRows1; i++) {
2966 if(
RowMap().GlobalIndicesInt()) {
2967 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 2968 int Row = (int)
GRID64(i);
2971 for (j = 0; j < NumIndices ; j++) {
2973 os << MyPID ; os <<
" ";
2975 os << Row ; os <<
" ";
2977 os << Indices_int[j]; os <<
" ";
2979 os << values[j]; os <<
" ";
2983 throw ReportError(
"Epetra_CrsMatrix::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
2986 else if(
RowMap().GlobalIndicesLongLong()) {
2987 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 2988 long long Row =
GRID64(i);
2991 for (j = 0; j < NumIndices ; j++) {
2993 os << MyPID ; os <<
" ";
2995 os << Row ; os <<
" ";
2997 os << Indices_LL[j]; os <<
" ";
2999 os << values[j]; os <<
" ";
3003 throw ReportError(
"Epetra_CrsMatrix::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
3008 if(
RowMap().GlobalIndicesInt()) {
3009 delete [] Indices_int;
3011 else if(
RowMap().GlobalIndicesLongLong()) {
3012 delete [] Indices_LL;
3030 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 3041 double* xp = (
double*) x.
Values();
3042 double* yp = (
double*) y.
Values();
3047 xp = (
double *) xcopy->
Values();
3109 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 3113 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3117 if(Epetra_CrsMatrixTraceDumpMultiply) {
3118 *out << std::boolalpha;
3119 *out <<
"\nEntering Epetra_CrsMatrix::Multipy("<<TransA<<
",X,Y) ...\n";
3121 *out <<
"\nDomainMap =\n";
3125 *out <<
"\nRangeMap =\n";
3128 *out <<
"\nInitial input X with " << ( TransA ?
"RangeMap" :
"DomainMap" ) <<
" =\n\n";
3131 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3145 double** Xp = (
double**) X.
Pointers();
3146 double** Yp = (
double**) Y.
Pointers();
3154 Xp = (
double **) Xcopy->
Pointers();
3167 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3168 if(Epetra_CrsMatrixTraceDumpMultiply) {
3169 *out <<
"\nColMap =\n";
3171 *out <<
"\nX after import from DomainMap to ColMap =\n\n";
3174 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3187 GeneralMM(Xp, LDX, Yp, LDY, NumVectors);
3188 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3189 if(Epetra_CrsMatrixTraceDumpMultiply) {
3190 *out <<
"\nRowMap =\n";
3192 *out <<
"\nY after local mat-vec where Y has RowMap =\n\n";
3198 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3202 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3203 if(Epetra_CrsMatrixTraceDumpMultiply) {
3204 *out <<
"\nRangeMap =\n";
3206 *out <<
"\nY after export from RowMap to RangeMap = \n\n";
3209 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3222 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3223 if(Epetra_CrsMatrixTraceDumpMultiply) {
3224 *out <<
"\nRowMap =\n";
3226 *out <<
"\nX after import from RangeMap to RowMap =\n\n";
3229 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3243 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3244 if(Epetra_CrsMatrixTraceDumpMultiply) {
3245 *out <<
"\nColMap =\n";
3247 *out <<
"\nY after local transpose mat-vec where Y has ColMap =\n\n";
3253 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3257 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3258 if(Epetra_CrsMatrixTraceDumpMultiply) {
3259 *out <<
"\nDomainMap =\n";
3261 *out <<
"\nY after export from ColMap to DomainMap =\n\n";
3264 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3277 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3278 if(Epetra_CrsMatrixTraceDumpMultiply) {
3279 *out <<
"\nFinal output Y is the last Y printed above!\n";
3280 *out <<
"\nLeaving Epetra_CrsMatrix::Multipy("<<TransA<<
",X,Y) ...\n";
3282 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY 3322 #if defined(Epetra_ENABLE_MKL_SPARSE) 3326 double alpha = 1, beta = 0;
3328 char matdescra[6] =
"G//C/";
3329 mkl_dcsrmv(&transa, &m, &NumCols, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, x, &beta, y);
3330 #elif defined(EPETRA_HAVE_OMP) 3332 #pragma omp parallel for default(none) shared(IndexOffset,values,Indices,y,x) 3333 for (
int row=0; row<numMyRows; ++row)
3335 const int curOffset = IndexOffset[row];
3336 const double *val_ptr = values+curOffset;
3337 const int *colnum_ptr = Indices+curOffset;
3339 const double *
const val_end_of_row = &values[IndexOffset[row+1]];
3340 while (val_ptr != val_end_of_row)
3341 s += *val_ptr++ * x[*colnum_ptr++];
3344 #elif defined(Epetra_ENABLE_CASK) 3345 cask_csr_dax_new(
NumMyRows_, IndexOffset, Indices,
3346 values, x, y, cask);
3347 #elif !defined(FORTRAN_DISABLED) 3353 const double *val_ptr = values;
3354 const int *colnum_ptr = Indices;
3355 double * dst_ptr = y;
3359 const double *
const val_end_of_row = &values[IndexOffset[row+1]];
3360 while (val_ptr != val_end_of_row)
3361 s += *val_ptr++ * x[*colnum_ptr++];
3364 #endif // Epetra_ENABLE_CASK or EPETRA_HAVE_OMP or Epetra_ENABLE_MKL_SPARSE 3373 double** srcValues =
Values();
3377 #ifdef EPETRA_HAVE_OMP 3378 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,y,x) 3380 for(
int i = 0; i < numMyRows; i++) {
3381 int NumEntries = NumEntriesPerRow[i];
3382 int* RowIndices = Indices[i];
3383 double* RowValues = srcValues[i];
3385 for(
int j = 0; j < NumEntries; j++)
3386 sum += *RowValues++ * x[*RowIndices++];
3397 #ifdef EPETRA_HAVE_OMP 3398 #pragma omp parallel for default(none) shared(x,y) 3400 for(
int i = 0; i < numMyRows; i++) {
3403 double* RowValues =
Values(i);
3405 for(
int j = 0; j < NumEntries; j++)
3406 sum += *RowValues++ * x[*RowIndices++];
3417 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE) 3423 #if defined(Epetra_ENABLE_MKL_SPARSE) 3426 double alpha = 1, beta = 0;
3428 char matdescra[6] =
"G//C/";
3429 mkl_dcsrmv(&transa, &m, &NumCols, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, x, &beta, y);
3430 #elif defined(Epetra_ENABLE_CASK) 3431 cask_csr_datx(
NumMyRows_, NumCols, IndexOffset, Indices, values,x ,y );
3439 #endif // FORTRAN_DISABLED etc 3441 for(
int i = 0; i < NumCols; i++)
3449 int prevOffset = *IndexOffset++;
3450 int NumEntries = *IndexOffset - prevOffset;
3452 for(
int j = 0; j < NumEntries; j++)
3453 y[*Indices++] += *values++ * xi;
3460 double** srcValues =
Values();
3463 int NumEntries = *NumEntriesPerRow++;
3464 int* RowIndices = *Indices++;
3465 double* RowValues = *srcValues++;
3467 for(
int j = 0; j < NumEntries; j++)
3468 y[*RowIndices++] += *RowValues++ * xi;
3476 double* RowValues =
Values(i);
3478 for(
int j = 0; j < NumEntries; j++)
3479 y[*RowIndices++] += *RowValues++ * xi;
3488 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)) 3494 if (LDX!=0 && LDY!=0) {
3495 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM) 3496 int * IndicesPlus1 =
Graph().All_IndicesPlus1();
3500 double alpha = 1, beta = 0;
3504 char matdescra[6] =
"G//F/";
3505 mkl_dcsrmm(&transa, &NumRows, &NumVectors, &NumCols, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *X, &LDX, &beta, *Y, &LDY);
3506 #elif defined(Epetra_ENABLE_CASK) 3508 IndexOffset, Indices, values, *X, LDX, 0.0, *Y, LDY,cask);
3511 EPETRA_DCRSMM_F77(&izero, &
NumMyRows_, &
NumMyRows_, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
3516 double **
const xp = X;
3517 double **
const yp = Y;
3519 #ifdef EPETRA_HAVE_OMP 3520 #pragma omp parallel for default(none) shared(IndexOffset,Indices,values,NumVectors) 3522 for (
int i=0; i < numMyRows; i++) {
3523 int prevOffset = IndexOffset[i];
3524 int NumEntries = IndexOffset[i+1] - prevOffset;
3525 int * RowIndices = Indices+prevOffset;
3526 double * RowValues = values+prevOffset;
3527 for (
int k=0; k<NumVectors; k++) {
3529 const double *
const x = xp[k];
3530 double *
const y = yp[k];
3531 for (
int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
3539 #endif // FORTRAN_DISABLED etc 3543 double** srcValues =
Values();
3544 double **
const xp = X;
3545 double **
const yp = Y;
3548 #ifdef EPETRA_HAVE_OMP 3549 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,NumVectors) 3551 for (
int i=0; i < numMyRows; i++) {
3552 int NumEntries = NumEntriesPerRow[i];
3553 int * RowIndices = Indices[i];
3554 double * RowValues = srcValues[i];
3555 for (
int k=0; k<NumVectors; k++) {
3557 const double *
const x = xp[k];
3558 double *
const y = yp[k];
3559 for (
int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
3566 double **
const xp = X;
3567 double **
const yp = Y;
3569 #ifdef EPETRA_HAVE_OMP 3570 #pragma omp parallel for default(none) shared(NumVectors) 3572 for (
int i=0; i < numMyRows; i++) {
3575 double* RowValues =
Values(i);
3576 for (
int k=0; k<NumVectors; k++) {
3579 for (
int j=0; j < NumEntries; j++) sum += RowValues[j] * x[RowIndices[j]];
3590 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)) 3592 if (LDX!=0 && LDY!=0) {
3596 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM) 3597 int * IndicesPlus1 =
Graph().All_IndicesPlus1();
3601 double alpha = 1, beta = 0;
3605 char matdescra[6] =
"G//F/";
3606 mkl_dcsrmm(&transa, &NumRows, &NumVectors, &NumCols, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *X, &LDX, &beta, *Y, &LDY);
3607 #elif defined(Epetra_ENABLE_CASK) 3608 cask_csr_dgesmm_new(1, 1.0,
NumMyRows_, NumCols, NumVectors,
3609 IndexOffset, Indices, values, *X, LDX, 0.0,
3613 EPETRA_DCRSMM_F77(&ione, &
NumMyRows_, &NumCols, values, Indices, IndexOffset, *X, &LDX, *Y, &LDY, &NumVectors);
3618 #endif // FORTRAN_DISABLED etc 3619 for (
int k=0; k<NumVectors; k++)
3620 for (
int i=0; i < NumCols; i++)
3628 int prevOffset = *IndexOffset++;
3629 int NumEntries = *IndexOffset - prevOffset;
3630 int * RowIndices = Indices+prevOffset;
3631 double * RowValues = values+prevOffset;
3633 for (
int k=0; k<NumVectors; k++) {
3636 for (
int j=0; j < NumEntries; j++)
3637 y[RowIndices[j]] += RowValues[j] * x[i];
3645 double** srcValues =
Values();
3648 int NumEntries = *NumEntriesPerRow++;
3649 int * RowIndices = *Indices++;
3650 double * RowValues = *srcValues++;
3651 for (
int k=0; k<NumVectors; k++) {
3654 for (
int j=0; j < NumEntries; j++)
3655 y[RowIndices[j]] += RowValues[j] * x[i];
3664 double* RowValues =
Values(i);
3665 for (
int k=0; k<NumVectors; k++) {
3668 for (
int j=0; j < NumEntries; j++)
3669 y[RowIndices[j]] += RowValues[j] * x[i];
3681 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE) 3687 #if !defined(Epetra_ENABLE_MKL_SPARSE) 3688 int iupper = Upper ? 1:0;
3689 int itrans = Trans ? 1:0;
3690 int udiag = UnitDiagonal ? 1:0;
3692 int xysame = (xp==yp) ? 1:0;
3695 #if defined(Epetra_ENABLE_MKL_SPARSE) 3696 char transa = Trans ?
't' :
'n';
3700 char matdescra[6] = {
'T', Upper ?
'U' :
'L', UnitDiagonal ?
'U' :
'N',
'C',
'/',
'\0'};
3701 mkl_dcsrsv(&transa, &m, &alpha, matdescra, values, Indices, IndexOffset, IndexOffset + 1, xp, yp);
3702 #elif defined(Epetra_ENABLE_CASK) 3703 cask_csr_dtrsv_new( iupper, itrans, udiag, nodiag, 0, xysame,
NumMyRows_,
3704 NumMyRows_, IndexOffset, Indices, values, xp, yp, cask);
3706 EPETRA_DCRSSV_F77( &iupper, &itrans, &udiag, &nodiag, &
NumMyRows_, &
NumMyRows_, values, Indices, IndexOffset, xp, yp, &xysame);
3724 double * RowValues =
Values(i);
3726 for (j=j0; j < NumEntries; j++)
3727 sum += RowValues[j] * yp[RowIndices[j]];
3730 yp[i] = xp[i] - sum;
3732 yp[i] = (xp[i] - sum)/RowValues[0];
3743 double * RowValues =
Values(i);
3745 for (j=0; j < NumEntries; j++)
3746 sum += RowValues[j] * yp[RowIndices[j]];
3749 yp[i] = xp[i] - sum;
3751 yp[i] = (xp[i] - sum)/RowValues[NumEntries];
3774 double * RowValues =
Values(i);
3776 yp[i] = yp[i]/RowValues[0];
3777 double ytmp = yp[i];
3778 for (j=j0; j < NumEntries; j++)
3779 yp[RowIndices[j]] -= RowValues[j] * ytmp;
3791 double * RowValues =
Values(i);
3793 yp[i] = yp[i]/RowValues[NumEntries];
3794 double ytmp = yp[i];
3795 for (j=0; j < NumEntries; j++)
3796 yp[RowIndices[j]] -= RowValues[j] * ytmp;
3801 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || defined(Epetra_ENABLE_MKL_SPARSE) 3816 #if !defined(FORTRAN_DISABLED) || defined(Epetra_ENABLE_CASK) || (defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM)) 3817 if (LDX!=0 && LDY!=0 && ((UnitDiagonal &&
NoDiagonal()) || (!UnitDiagonal && !
NoDiagonal()))) {
3819 #if !defined(Epetra_ENABLE_MKL_SPARSE) || defined(Epetra_DISABLE_MKL_SPARSE_MM) 3820 int iupper = Upper ? 1:0;
3821 int itrans = Trans ? 1:0;
3822 int udiag = UnitDiagonal ? 1:0;
3824 int xysame = (Xp==Yp) ? 1:0;
3827 #if defined(Epetra_ENABLE_MKL_SPARSE) && !defined(Epetra_DISABLE_MKL_SPARSE_MM) 3828 int * IndicesPlus1 =
Graph().All_IndicesPlus1();
3829 char transa = Trans ?
't' :
'n';
3835 char matdescra[6] = {
'T', Upper ?
'U' :
'L', UnitDiagonal ?
'U' :
'N',
'F',
'/',
'\0'};
3836 mkl_dcsrsm(&transa, &NumRows, &NumVectors, &alpha, matdescra, values, IndicesPlus1, IndexOffset, IndexOffset + 1, *Xp, &LDX, *Yp, &LDY);
3837 #elif defined(Epetra_ENABLE_CASK) 3838 cask_csr_dtrsm( iupper, itrans, udiag, nodiag, 0, xysame,
NumMyRows_,
3839 NumMyRows_, NumVectors, IndexOffset, Indices, values,
3840 *Xp, LDX, *Yp, LDY);
3843 *Xp, &LDX, *Yp, &LDY, &xysame, &NumVectors);
3847 #endif // FORTRAN_DISABLED etc 3854 int Offset = IndexOffset[i];
3855 int NumEntries = IndexOffset[i+1]-Offset;
3856 int * RowIndices = Indices+Offset;
3857 double * RowValues = values+Offset;
3859 diag = 1.0/RowValues[0];
3860 for(k = 0; k < NumVectors; k++) {
3862 for(j = j0; j < NumEntries; j++)
3863 sum += RowValues[j] * Yp[k][RowIndices[j]];
3866 Yp[k][i] = Xp[k][i] - sum;
3868 Yp[k][i] = (Xp[k][i] - sum) * diag;
3877 int Offset = IndexOffset[i];
3878 int NumEntries = IndexOffset[i+1]-Offset - j0;
3879 int * RowIndices = Indices+Offset;
3880 double * RowValues = values+Offset;
3882 diag = 1.0/RowValues[NumEntries];
3883 for(k = 0; k < NumVectors; k++) {
3885 for(j = 0; j < NumEntries; j++)
3886 sum += RowValues[j] * Yp[k][RowIndices[j]];
3889 Yp[k][i] = Xp[k][i] - sum;
3891 Yp[k][i] = (Xp[k][i] - sum)*diag;
3899 for(k = 0; k < NumVectors; k++)
3902 Yp[k][i] = Xp[k][i];
3910 int Offset = IndexOffset[i];
3911 int NumEntries = IndexOffset[i+1]-Offset;
3912 int * RowIndices = Indices+Offset;
3913 double * RowValues = values+Offset;
3915 diag = 1.0/RowValues[0];
3916 for(k = 0; k < NumVectors; k++) {
3918 Yp[k][i] = Yp[k][i]*diag;
3919 double ytmp = Yp[k][i];
3920 for(j = j0; j < NumEntries; j++)
3921 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
3930 int Offset = IndexOffset[i];
3931 int NumEntries = IndexOffset[i+1]-Offset - j0;
3932 int * RowIndices = Indices+Offset;
3933 double * RowValues = values+Offset;
3935 diag = 1.0/RowValues[NumEntries];
3936 for(k = 0; k < NumVectors; k++) {
3938 Yp[k][i] = Yp[k][i]*diag;
3939 double ytmp = Yp[k][i];
3940 for(j = 0; j < NumEntries; j++)
3941 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
3959 double* RowValues =
Values(i);
3961 diag = 1.0/RowValues[0];
3962 for(k = 0; k < NumVectors; k++) {
3964 for(j = j0; j < NumEntries; j++)
3965 sum += RowValues[j] * Yp[k][RowIndices[j]];
3968 Yp[k][i] = Xp[k][i] - sum;
3970 Yp[k][i] = (Xp[k][i] - sum) * diag;
3981 double* RowValues =
Values(i);
3983 diag = 1.0/RowValues[NumEntries];
3984 for(k = 0; k < NumVectors; k++) {
3986 for(j = 0; j < NumEntries; j++)
3987 sum += RowValues[j] * Yp[k][RowIndices[j]];
3990 Yp[k][i] = Xp[k][i] - sum;
3992 Yp[k][i] = (Xp[k][i] - sum)*diag;
4000 for(k = 0; k < NumVectors; k++)
4003 Yp[k][i] = Xp[k][i];
4013 double* RowValues =
Values(i);
4015 diag = 1.0/RowValues[0];
4016 for(k = 0; k < NumVectors; k++) {
4018 Yp[k][i] = Yp[k][i]*diag;
4019 double ytmp = Yp[k][i];
4020 for(j = j0; j < NumEntries; j++)
4021 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
4032 double* RowValues =
Values(i);
4034 diag = 1.0/RowValues[NumEntries];
4035 for(k = 0; k < NumVectors; k++) {
4037 Yp[k][i] = Yp[k][i]*diag;
4038 double ytmp = Yp[k][i];
4039 for(j = 0; j < NumEntries; j++)
4040 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
4052 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 4064 double* xp = (
double*) x.
Values();
4065 double* yp = (
double*) y.
Values();
4101 double* RowValues =
Values(i);
4103 for(j = 0; j < NumEntries; j++)
4104 sum += RowValues[j] * xp[RowIndices[j]];
4147 for(i = 0; i < NumMyCols_; i++)
4153 double* RowValues =
Values(i);
4154 for(j = 0; j < NumEntries; j++)
4155 yp[RowIndices[j]] += RowValues[j] * xp[i];
4171 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 4179 double* xp = (
double*) X[0];
4180 double* yp = (
double*) Y[0];
4192 double** Xp = (
double**) X.
Pointers();
4193 double** Yp = (
double**) Y.
Pointers();
4234 double * RowValues =
Values(i);
4235 for (k=0; k<NumVectors; k++) {
4237 for (j=0; j < NumEntries; j++) sum += RowValues[j] * Xp[k][RowIndices[j]];
4279 for (k=0; k<NumVectors; k++)
4280 for (i=0; i < NumMyCols_; i++)
4286 double * RowValues =
Values(i);
4287 for (k=0; k<NumVectors; k++) {
4288 for (j=0; j < NumEntries; j++)
4289 Yp[k][RowIndices[j]] += RowValues[j] * Xp[k][i];
4309 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 4334 double ** Vals =
Values();
4338 if ((Upper && !Trans) || (!Upper && Trans)) {
4344 double *xp = (
double*)x.
Values();
4345 double *yp = (
double*)y.
Values();
4355 int NumEntries = *NumEntriesPerRow--;
4356 int * RowIndices = *Indices--;
4357 double * RowValues = *Vals--;
4359 for (j=j0; j < NumEntries; j++)
4360 sum += RowValues[j] * yp[RowIndices[j]];
4363 yp[i] = xp[i] - sum;
4365 yp[i] = (xp[i] - sum)/RowValues[0];
4374 int NumEntries = *NumEntriesPerRow++ - j0;
4375 int * RowIndices = *Indices++;
4376 double * RowValues = *Vals++;
4378 for (j=0; j < NumEntries; j++)
4379 sum += RowValues[j] * yp[RowIndices[j]];
4382 yp[i] = xp[i] - sum;
4384 yp[i] = (xp[i] - sum)/RowValues[NumEntries];
4395 for (i=0; i < NumMyCols_; i++)
4405 int NumEntries = *NumEntriesPerRow++;
4406 int * RowIndices = *Indices++;
4407 double * RowValues = *Vals++;
4409 yp[i] = yp[i]/RowValues[0];
4410 double ytmp = yp[i];
4411 for (j=j0; j < NumEntries; j++)
4412 yp[RowIndices[j]] -= RowValues[j] * ytmp;
4422 int NumEntries = *NumEntriesPerRow-- - j0;
4423 int * RowIndices = *Indices--;
4424 double * RowValues = *Vals--;
4426 yp[i] = yp[i]/RowValues[NumEntries];
4427 double ytmp = yp[i];
4428 for (j=0; j < NumEntries; j++)
4429 yp[RowIndices[j]] -= RowValues[j] * ytmp;
4441 #ifdef EPETRA_CRSMATRIX_TEUCHOS_TIMERS 4449 double* xp = (
double*) X[0];
4450 double* yp = (
double*) Y[0];
4471 double** Vals =
Values();
4475 if((Upper && !Trans) || (!Upper && Trans)) {
4481 double** Xp = (
double**) X.
Pointers();
4482 double** Yp = (
double**) Y.
Pointers();
4492 int NumEntries = *NumEntriesPerRow--;
4493 int* RowIndices = *Indices--;
4494 double* RowValues = *Vals--;
4496 diag = 1.0/RowValues[0];
4497 for(k = 0; k < NumVectors; k++) {
4499 for(j = j0; j < NumEntries; j++)
4500 sum += RowValues[j] * Yp[k][RowIndices[j]];
4503 Yp[k][i] = Xp[k][i] - sum;
4505 Yp[k][i] = (Xp[k][i] - sum) * diag;
4514 int NumEntries = *NumEntriesPerRow++ - j0;
4515 int* RowIndices = *Indices++;
4516 double* RowValues = *Vals++;
4518 diag = 1.0/RowValues[NumEntries];
4519 for(k = 0; k < NumVectors; k++) {
4521 for(j = 0; j < NumEntries; j++)
4522 sum += RowValues[j] * Yp[k][RowIndices[j]];
4525 Yp[k][i] = Xp[k][i] - sum;
4527 Yp[k][i] = (Xp[k][i] - sum)*diag;
4535 for(k = 0; k < NumVectors; k++)
4538 Yp[k][i] = Xp[k][i];
4546 int NumEntries = *NumEntriesPerRow++;
4547 int* RowIndices = *Indices++;
4548 double* RowValues = *Vals++;
4550 diag = 1.0/RowValues[0];
4551 for(k = 0; k < NumVectors; k++) {
4553 Yp[k][i] = Yp[k][i]*diag;
4554 double ytmp = Yp[k][i];
4555 for(j = j0; j < NumEntries; j++)
4556 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
4565 int NumEntries = *NumEntriesPerRow-- - j0;
4566 int* RowIndices = *Indices--;
4567 double* RowValues = *Vals--;
4569 diag = 1.0/RowValues[NumEntries];
4570 for(k = 0; k < NumVectors; k++) {
4572 Yp[k][i] = Yp[k][i]*diag;
4573 double ytmp = Yp[k][i];
4574 for(j = 0; j < NumEntries; j++)
4575 Yp[k][RowIndices[j]] -= RowValues[j] * ytmp;
4610 int m=
D.RowMap_.NumMyElements();
4613 if(
D.RowMap_.GlobalIndicesLongLong()) UseLL=
true;
4616 if(!
D.RowMap_.ConstantElementSize() ||
D.RowMap_.ElementSize()!=1 ||
4617 !
D.ColMap_.ConstantElementSize() ||
D.ColMap_.ElementSize()!=1)
4621 D.DomainMap_ = theDomainMap;
4622 D.RangeMap_ = theRangeMap;
4625 if (!
D.ColMap_.SameAs(
D.DomainMap_)) {
4626 if (
D.Importer_ != 0) {
4631 D.Importer_=theImporter;
4638 if (
D.Importer_ != 0) {
4645 if (!
D.RowMap_.SameAs(
D.RangeMap_)) {
4646 if (
D.Exporter_ != 0) {
4651 D.Exporter_=theExporter;
4658 if (
D.Exporter_ != 0) {
4677 for(
int i=0;i<m;i++){
4680 D.data->SortedEntries_.clear();
4684 delete []
D.data->Indices_;
D.data->Indices_=0;
4685 D.NumAllocatedIndicesPerRow_.Resize(0);
4690 D.Allocated_ =
true;
4692 D.StorageOptimized_ =
true;
4693 D.NoRedundancies_ =
true;
4694 D.IndicesAreGlobal_ =
false;
4695 D.IndicesAreLocal_ =
true;
4696 D.IndicesAreContiguous_ =
true;
4697 D.GlobalConstantsComputed_ =
true;
4698 D.StaticProfile_ =
true;
4699 D.SortGhostsAssociatedWithEachProcessor_ =
true;
4700 D.HaveColMap_ =
true;
4705 int nnz=
D.IndexOffset_[m]-
D.IndexOffset_[0];
4706 D.NumMyRows_ =
D.NumMyBlockRows_ = m;
4707 D.NumMyCols_ =
D.NumMyBlockCols_ =
D.ColMap_.NumMyElements();
4708 D.NumMyNonzeros_ =
D.NumMyEntries_ = nnz;
4713 D.GlobalMaxRowDim_ = 1;
4714 D.GlobalMaxColDim_ = 1;
4718 for(
int i=0; i<m; i++){
4719 int NumIndices=
D.IndexOffset_[i+1]-
D.IndexOffset_[i];
4720 D.MaxNumIndices_=
EPETRA_MAX(
D.MaxNumIndices_,NumIndices);
4723 if(NumIndices > 0) {
4724 const int* col_indices = &
D.data->All_Indices_[
D.IndexOffset_[i]];
4726 int jl_0 = col_indices[0];
4727 int jl_n = col_indices[NumIndices-1];
4729 if(jl_n > i)
D.LowerTriangular_ =
false;
4730 if(jl_0 < i)
D.UpperTriangular_ =
false;
4733 if(numMyDiagonals == -1) {
4737 int insertPoint = -1;
4739 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 4740 long long jg =
D.RowMap_.GID64(i);
4742 D.NumMyBlockDiagonals_++;
4743 D.NumMyDiagonals_ ++;
4748 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 4749 int jg =
D.RowMap_.GID(i);
4751 D.NumMyBlockDiagonals_++;
4752 D.NumMyDiagonals_ ++;
4760 if(numMyDiagonals > -1)
D.NumMyDiagonals_ =
D.NumMyBlockDiagonals_ = numMyDiagonals;
4762 D.MaxNumNonzeros_=
D.MaxNumIndices_;
4766 tempvec[0] =
D.NumMyEntries_;
4767 tempvec[1] =
D.NumMyBlockDiagonals_;
4768 tempvec[2] =
D.NumMyDiagonals_;
4769 tempvec[3] =
D.NumMyNonzeros_;
4773 D.NumGlobalEntries_ = tempvec[4];
4774 D.NumGlobalBlockDiagonals_ = tempvec[5];
4775 D.NumGlobalDiagonals_ = tempvec[6];
4776 D.NumGlobalNonzeros_ = tempvec[7];
4778 tempvec[0] =
D.MaxNumIndices_;
4779 tempvec[1] =
D.MaxNumNonzeros_;
4783 D.GlobalMaxNumIndices_ = tempvec[2];
4784 D.GlobalMaxNumNonzeros_ = tempvec[3];
4790 D.NumGlobalRows_ =
D.RangeMap_.NumGlobalPoints64();
4791 D.NumGlobalCols_ =
D.DomainMap_.NumGlobalPoints64();
4794 D.NumGlobalRows_ = (int)
D.RangeMap_.NumGlobalPoints64();
4795 D.NumGlobalCols_ = (int)
D.DomainMap_.NumGlobalPoints64();
4801 template<
class TransferType>
4804 const TransferType& RowTransfer,
4805 const TransferType* DomainTransfer,
4808 const bool RestrictCommunicator)
4813 bool communication_needed = RowTransfer.SourceMap().DistributedGlobal();
4823 throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor unable to determine source global index type",-1);
4833 if (!SourceMatrix.
RowMap().
SameAs(RowTransfer.SourceMap()))
4834 throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor requires RowTransfer.SourceMap() to match SourceMatrix.RowMap()",-2);
4837 std::vector<int> SourcePids;
4838 std::vector<int> TargetPids;
4846 const Epetra_Map* BaseDomainMap = MyDomainMap;
4849 Epetra_Map *ReducedRowMap=0, *ReducedDomainMap=0, *ReducedRangeMap=0;
4855 if(RestrictCommunicator) {
4857 ReducedComm = ReducedRowMap ? &ReducedRowMap->
Comm() : 0;
4864 MyRowMap = ReducedRowMap;
4865 MyDomainMap = ReducedDomainMap;
4866 MyRangeMap = ReducedRangeMap;
4869 if(ReducedComm) MyPID = ReducedComm->
MyPID();
4882 bool bSameDomainMap = BaseDomainMap->
SameAs(SourceMatrix.
DomainMap());
4885 if(!RestrictCommunicator && MyImporter && bSameDomainMap) {
4892 else if (RestrictCommunicator && MyImporter && bSameDomainMap) {
4902 SourceCol_pids.Import(SourceDomain_pids,*MyImporter,
Insert);
4904 else if(!MyImporter && bSameDomainMap) {
4911 else if (MyImporter && DomainTransfer) {
4929 if(
typeid(TransferType)==
typeid(
Epetra_Import) && DomainTransfer->TargetMap().SameBlockMapDataAs(*theDomainMap))
4930 SourceDomain_pids.
Export(TargetDomain_pids,*DomainTransfer,
Insert);
4931 else if(
typeid(TransferType)==
typeid(
Epetra_Export) && DomainTransfer->SourceMap().SameBlockMapDataAs(*theDomainMap))
4932 SourceDomain_pids.Import(TargetDomain_pids,*DomainTransfer,
Insert);
4934 throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor TransferType must be Epetra_Import or Epetra_Export",-31);
4936 SourceCol_pids.Import(SourceDomain_pids,*MyImporter,
Insert);
4950 SourceRow_pids.Export(TargetRow_pids,RowTransfer,
Insert);
4952 SourceRow_pids.Import(TargetRow_pids,RowTransfer,
Insert);
4954 throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor TransferType must be Epetra_Import or Epetra_Export",-31);
4955 SourceCol_pids.Import(SourceRow_pids,*MyImporter,
Insert);
4959 std::cout <<
"Error in FusedTransfer:" << std::endl;
4962 std::cout <<
"BaseRowMap " << BaseRowMap->
NumMyElements() << std::endl;
4963 std::cout <<
"BaseDomainMap" << BaseDomainMap->
NumMyElements() << std::endl;
4964 if(DomainTransfer == NULL) std::cout <<
"DomainTransfer = NULL" << std::endl;
4965 else std::cout <<
"DomainTransfer is not NULL" << std::endl;
4966 throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor only supports *theDomainMap==SourceMatrix.DomainMap() || DomainTransfer != NULL || *theDomainMap==RowTransfer.TargetMap() && SourceMatrix.DomainMap() == SourceMatrix.RowMap().",-4);
4970 int NumSameIDs = RowTransfer.NumSameIDs();
4971 int NumPermuteIDs = RowTransfer.NumPermuteIDs();
4972 int NumRemoteIDs = RowTransfer.NumRemoteIDs();
4973 int NumExportIDs = RowTransfer.NumExportIDs();
4974 int* ExportLIDs = RowTransfer.ExportLIDs();
4975 int* RemoteLIDs = RowTransfer.RemoteLIDs();
4976 int* PermuteToLIDs = RowTransfer.PermuteToLIDs();
4977 int* PermuteFromLIDs = RowTransfer.PermuteFromLIDs();
4987 std::vector<long long> CSR_colind_LL;
4990 if(rv)
throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor failed in CheckSizes()",-3);
4993 bool VarSizes =
false;
4994 if( NumExportIDs > 0) {
4996 Sizes_ =
new int[NumExportIDs];
5001 NumExportIDs,ExportLIDs,
5003 Sizes_,VarSizes,SourcePids);
5004 if(rv)
throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor failed in PackAndPrepareWithOwningPIDs()",-5);
5006 if (communication_needed) {
5013 if(rv)
throw ReportError(
"Epetra_CrsMatrix: Fused import/export constructor failed in Distor.Do",-6);
5023 CSR_colind.
Resize(mynnz);
5024 if(UseLL) CSR_colind_LL.resize(mynnz);
5026 CSR_vals =
new double[mynnz];
5030 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,
LenImports_,
Imports_,
NumMyRows(),mynnz,MyPID,CSR_rowptr.
Values(),
Epetra_Util_data_ptr(CSR_colind_LL),CSR_vals,SourcePids,TargetPids);
5032 Epetra_Import_Util::UnpackAndCombineIntoCrsArrays(SourceMatrix,NumSameIDs,NumRemoteIDs,RemoteLIDs,NumPermuteIDs,PermuteToLIDs,PermuteFromLIDs,
LenImports_,
Imports_,
NumMyRows(),mynnz,MyPID,CSR_rowptr.
Values(),CSR_colind.
Values(),CSR_vals,SourcePids,TargetPids);
5037 if(RestrictCommunicator) {
5050 #if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES) 5069 std::vector<int> RemotePIDs;
5072 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES 5075 *MyDomainMap,pids_ptr,
5082 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES 5105 int NumRemotePIDs = RemotePIDs.size();
5119 if(ReducedDomainMap!=ReducedRowMap)
delete ReducedDomainMap;
5120 if(ReducedRangeMap !=ReducedRowMap)
delete ReducedRangeMap;
5121 delete ReducedRowMap;
5131 Graph_(
Copy, RowImporter.TargetMap(), 0, false),
5135 NumMyRows_(RowImporter.TargetMap().NumMyPoints()),
5140 FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,NULL,theDomainMap,theRangeMap,RestrictCommunicator);
5149 Graph_(
Copy, RowExporter.TargetMap(), 0, false),
5151 Values_alloc_lengths_(0),
5153 NumMyRows_(RowExporter.TargetMap().NumMyPoints()),
5158 FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,NULL,theDomainMap,theRangeMap,RestrictCommunicator);
5166 Graph_(
Copy, RowImporter.TargetMap(), 0, false),
5168 Values_alloc_lengths_(0),
5170 NumMyRows_(RowImporter.TargetMap().NumMyPoints()),
5175 FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,DomainImporter,theDomainMap,theRangeMap,RestrictCommunicator);
5184 Graph_(
Copy, RowExporter.TargetMap(), 0, false),
5186 Values_alloc_lengths_(0),
5188 NumMyRows_(RowExporter.TargetMap().NumMyPoints()),
5193 FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,DomainExporter,theDomainMap,theRangeMap,RestrictCommunicator);
5201 bool RestrictCommunicator) {
5202 FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,NULL,theDomainMap,theRangeMap,RestrictCommunicator);
5210 bool RestrictCommunicator) {
5211 FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,NULL,theDomainMap,theRangeMap,RestrictCommunicator);
5219 bool RestrictCommunicator) {
5220 FusedTransfer<Epetra_Import>(SourceMatrix,RowImporter,DomainImporter,theDomainMap,theRangeMap,RestrictCommunicator);
5229 bool RestrictCommunicator) {
5230 FusedTransfer<Epetra_Export>(SourceMatrix,RowExporter,DomainExporter,theDomainMap,theRangeMap,RestrictCommunicator);
bool Sorted() const
If SortEntries() has been called, this query returns true, otherwise it returns false.
const Epetra_Map & RowMap() const
Returns the Epetra_Map object associated with the rows of this matrix.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
int TransformToLocal()
Use FillComplete() instead.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
bool FindMyIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
bool HaveColMap() const
Returns true if we have a well-defined ColMap, and returns false otherwise.
double * Values() const
Get pointer to MultiVector values.
int GlobalMaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on all processors.
const Epetra_BlockMap & Map() const
Returns the address of the Epetra_BlockMap for this multi-vector.
int MaxNumEntries() const
Returns the maximum number of nonzero entries across all rows on this processor.
int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified local row in user-provided arrays.
Epetra_Map: A class for partitioning vectors and matrices.
bool MyLID(int lid) const
Returns true if the LID passed in belongs to the calling processor in this map, otherwise returns fal...
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Epetra_BlockMap RangeMap_
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
void UpdateFlops(int Flops_in) const
Increment Flop count for this object.
void GeneralMTM(double **X, int LDX, double **Y, int LDY, int NumVectors) const
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
void FusedTransfer(const Epetra_CrsMatrix &SourceMatrix, const TransferType &RowTransfer, const TransferType *DomainTransfer, const Epetra_Map *DomainMap, const Epetra_Map *RangeMap, bool RestrictCommunicator)
void InitializeDefaults()
int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int Resize(int Length_in)
Resize a Epetra_IntSerialDenseVector object.
const Epetra_BlockMapData * DataPtr() const
Returns a pointer to the BlockMapData instance this BlockMap uses.
int ReferenceCount() const
Get reference count.
int SumIntoOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices)
int Multiply1(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
bool FindGlobalIndexLoc(int LocalRow, int Index, int Start, int &Loc) const
long long NumGlobalNonzeros64() const
int Multiply(bool TransA, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a Epetra_CrsMatrix multiplied by a Epetra_Vector x in y.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.
int NumMyDiagonals() const
Returns the number of local nonzero diagonal entries, based on global row/column index comparisons...
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
bool DistributedGlobal() const
Returns true if map is defined across more than one processor.
virtual void Print(std::ostream &os) const
Print method.
int * Values_alloc_lengths_
int ExtractGlobalRowView(int GlobalRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified global row values via pointers to internal data.
long long NumGlobalRows64() const
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
int ReplaceOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices)
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
int SortEntries()
Sort column entries, row-by-row, in ascending order.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
void GeneralMV(double *x, double *y) const
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
void DecrementReferenceCount()
Decrement reference count.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int TCopyAndPermuteRowMatrix(const Epetra_RowMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
double NormOne() const
Returns the one norm of the global matrix.
int PackAndPrepare(const Epetra_SrcDistObject &Source, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, Epetra_Distributor &Distor)
Perform any packing or preparation required for call to DoTransfer().
int NumMyNonzeros() const
Returns the number of indices in the local graph.
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this importer.
int PutScalar(double ScalarConstant)
Initialize all values in a multi-vector with constant value.
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given global row of the matrix.
int NumAllocatedMyIndices(int Row) const
Returns the allocated number of nonzero entries in specified local row on this processor.
bool SortGhostsAssociatedWithEachProcessor_
int InsertOffsetValues(long long GlobalRow, int NumEntries, double *Values, int *Indices)
double * All_Values() const
virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Add this list of entries to existing values for a given global row of the matrix. ...
T * Epetra_Util_data_ptr(std::vector< T > &vec)
Function that returns either a pointer to the first entry in the vector or, if the vector is empty...
Epetra_BlockMap DomainMap_
#define EPETRA_CHK_ERR(a)
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true)...
int MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
int LCID(int GCID_in) const
Returns the local column index for given global column index, returns -1 if no local column for this ...
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this importer.
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
int UnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
Perform any unpacking and combining after call to DoTransfer().
bool constructedWithFilledGraph_
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements...
long long NumGlobalDiagonals64() const
int PutValue(int Value)
Set all elements of the vector to Value.
int MakeIndicesLocal(const Epetra_BlockMap &DomainMap, const Epetra_BlockMap &RangeMap)
int ReplaceRowMap(const Epetra_BlockMap &newmap)
Replaces the current RowMap with the user-specified map object, but only if currentmap->PointSameAs(n...
Epetra_Vector: A class for constructing and using dense vectors on a parallel computer.
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
bool IndicesAreLocal() const
If matrix indices has been transformed to local, this query returns true, otherwise it returns false...
const Epetra_Comm & Comm() const
Returns a pointer to the Epetra_Comm communicator associated with this matrix.
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
virtual int ReportError(const std::string Message, int ErrorCode) const
Error reporting method.
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
Epetra_Comm Global Min function.
long long NumGlobalElements64() const
int ExtractGlobalRowView(int GlobalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified global row of the graph.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements...
virtual void Barrier() const =0
Epetra_Comm Barrier function.
const Epetra_BlockMap & RowMap() const
Returns the RowMap associated with this graph.
int ReplaceDomainMapAndImporter(const Epetra_Map &NewDomainMap, const Epetra_Import *NewImporter)
Replaces the current DomainMap & Importer with the user-specified map object.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Returns a view of the specified local row values via pointers to internal data.
int PackAndPrepareWithOwningPIDs(const Epetra_CrsMatrix &SourceMatrix, int NumExportIDs, int *ExportLIDs, int &LenExports, char *&Exports, int &SizeOfPacket, int *Sizes, bool &VarSizes, std::vector< int > &SourcePids)
PackAndPrepareWithOwningPIDs.
static int SortAndMergeCrsEntries(int NumRows, int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortAndMergeCrsEntries function.
void PREFIX EPETRA_DCRSSM_F77(const int *, const int *, const int *, const int *, const int *, const int *, const double *, const int *, const int *, double *, const int *, double *, const int *, const int *, const int *)
int ReplaceDomainMapAndImporter(const Epetra_BlockMap &NewDomainMap, const Epetra_Import *NewImporter)
Replaces the current DomainMap & Importer with the user-specified map object.
void FusedImport(const Epetra_CrsMatrix &SourceMatrix, const Epetra_Import &RowImporter, const Epetra_Map *DomainMap, const Epetra_Map *RangeMap, bool RestrictCommunicator)
static int GetPids(const Epetra_Import &Importer, std::vector< int > &pids, bool use_minus_one_for_local)
Epetra_Util GetPids function.
Epetra_CrsGraphData * CrsGraphData_
int InvRowMaxs(Epetra_Vector &x) const
Computes the inverse of the max of absolute values of the rows of the Epetra_CrsMatrix, results returned in x.
virtual int MyPID() const =0
Return my process ID.
bool squareFillCompleteCalled_
const Epetra_Map & RangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space. ...
int RemoveRedundantIndices()
Removes any redundant column indices in the rows of the graph.
const Epetra_Export * Exporter() const
Returns the Epetra_Export object that contains the export operations for distributed operations...
int ExpertMakeUniqueCrsGraphData()
Makes sure this matrix has a unique CrsGraphData object.
virtual int MaxAll(double *PartialMaxs, double *GlobalMaxs, int Count) const =0
Epetra_Comm Global Max function.
int * IndexOffset() const
Epetra_CrsGraphData: The Epetra CrsGraph Data Class.
bool StaticGraph()
Returns true if the graph associated with this matrix was pre-constructed and therefore not changeabl...
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
void GeneralMM(double **X, int LDX, double **Y, int LDY, int NumVectors) const
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true...
const Epetra_BlockMap & ColMap() const
Returns the Column Map associated with this graph.
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
Epetra_BLAS: The Epetra BLAS Wrapper Class.
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
Epetra_Comm Global Sum function.
int UnpackWithOwningPIDsCount(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports)
UnpackWithOwningPIDsCount.
int InvColMaxs(Epetra_Vector &x) const
Computes the max of absolute values of the columns of the Epetra_CrsMatrix, results returned in x...
int ** PermuteOffsets() const
Accessor.
static int SortCrsEntries(int NumRows, const int *CRS_rowptr, int *CRS_colind, double *CRS_vals)
Epetra_Util SortCrsEntries function.
bool MyGID(int GID_in) const
Returns true if the GID passed in belongs to the calling processor in this map, otherwise returns fal...
Epetra_MultiVector * ExportVector_
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
int InsertIndices(int Row, int NumIndices, int *Indices)
void PREFIX EPETRA_DCRSSV_F77(const int *, const int *, const int *, const int *, const int *, const int *, const double *, const int *, const int *, double *, double *, const int *)
int UnpackAndCombineIntoCrsArrays(const Epetra_CrsMatrix &SourceMatrix, int NumSameIDs, int NumRemoteIDs, const int *RemoteLIDs, int NumPermuteIDs, const int *PermuteToLIDs, const int *PermuteFromLIDs, int LenImports, char *Imports, int TargetNumRows, int TargetNumNonzeros, int MyTargetPID, int *CSR_rowptr, int *CSR_colind, double *CSR_values, const std::vector< int > &SourcePids, std::vector< int > &TargetPids)
UnpackAndCombineIntoCrsArrays.
Epetra_Comm: The Epetra Communication Abstract Base Class.
long long GID64(int LID) const
int CopyAndPermute(const Epetra_SrcDistObject &Source, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode=Zero)
Perform ID copies and permutations that are on processor.
void SetIndicesAreGlobal(bool Flag)
virtual ~Epetra_CrsMatrix()
Epetra_CrsMatrix Destructor.
void SetIndicesAreLocal(bool Flag)
void UpdateImportVector(int NumVectors) const
const Epetra_BlockMap & RangeMap() const
Returns the RangeMap associated with this graph.
int ** SameOffsets() const
Accessor.
int InvRowSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the rows of the Epetra_CrsMatrix, results returned in x.
const Epetra_Comm * Comm_
int ReplaceRowMap(const Epetra_BlockMap &newmap)
Replaces the current RowMap with the user-specified map object.
int NumVectors() const
Returns the number of vectors in the multi-vector.
int * All_Indices() const
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
double NormFrobenius() const
Returns the frobenius norm of the global matrix.
int InsertValues(int LocalRow, int NumEntries, double *Values, int *Indices)
virtual void Print(std::ostream &os) const
Print object to an output stream.
int OptimizeStorage()
Make consecutive row index sections contiguous, minimize internal storage used for constructing graph...
int TInsertGlobalValues(int_type Row, int NumEntries, const double *values, const int_type *Indices)
int TSumIntoGlobalValues(int_type Row, int NumEntries, const double *srcValues, const int_type *Indices)
bool GlobalIndicesInt() const
Returns true if map create with int NumGlobalElements.
virtual const Epetra_BlockMap & Map() const =0
Returns a reference to the Epetra_BlockMap for this object.
int NumMyElements() const
Number of elements on the calling processor.
const double Epetra_MaxDouble
int * NumIndicesPerRow() const
static RCP< FancyOStream > getDefaultOStream()
void FusedExport(const Epetra_CrsMatrix &SourceMatrix, const Epetra_Export &RowExporter, const Epetra_Map *DomainMap, const Epetra_Map *RangeMap, bool RestrictCommunicator)
Epetra_CompObject: Functionality and data that is common to all computational classes.
Epetra_BlockMap: A class for partitioning block element vectors and matrices.
int RemoveEmptyProcessesInPlace(const Epetra_BlockMap *NewMap)
Remove processes owning zero rows from the Maps and their communicator.
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int Epetra_Util_binary_search_aux(T item, const int *list, const T *aux_list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
long long GRID64(int LRID_in) const
int ** RemoteOffsets() const
Accessor.
int LowCommunicationMakeColMapAndReindex(int N, const int *rowptr, int *colind, const Epetra_Map &domainMap, const int *owningPIDs, bool SortGhostsAssociatedWithEachProcessor, std::vector< int > &RemotePIDs, Epetra_BlockMap &NewColMap)
LowCommunicationMakeColMapAndReindex.
int Epetra_Util_binary_search(T item, const T *list, int len, int &insertPoint)
Utility function to perform a binary-search on a list of data.
int LeftScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the left with a Epetra_Vector x.
bool LowerTriangular() const
If matrix is lower triangular in local index space, this query returns true, otherwise it returns fal...
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
const Epetra_Map & ColMap() const
Returns the Epetra_Map object that describes the set of column-indices that appear in each processor'...
Epetra_IntSerialDenseVector & ExpertExtractIndices()
Returns a reference to the Epetra_IntSerialDenseVector used to hold the local All_Indices (CRS colind...
Epetra_IntSerialDenseVector All_Indices_
int TReplaceGlobalValues(int_type Row, int NumEntries, const double *srcValues, const int_type *Indices)
virtual int Do(char *export_objs, int obj_size, int &len_import_objs, char *&import_objs)=0
Execute plan on buffer of export objects in a single step.
Epetra_SerialComm: The Epetra Serial Communication Class.
int LID(int GID) const
Returns local ID of global ID, return -1 if not found on this processor.
virtual void Print(std::ostream &os) const
Print method.
void PREFIX EPETRA_DCRSMM_F77(const int *, const int *, const int *, const double *, const int *, const int *, double *, int *, double *, int *, int *)
bool matrixFillCompleteCalled_
int SetAllocated(bool Flag)
int SumIntoMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Add this list of entries to existing values for a given local row of the matrix.
Epetra_IntSerialDenseVector & ExpertExtractIndexOffset()
Returns a reference to the Epetra_IntSerialDenseVector used to hold the local IndexOffsets (CRS rowpt...
const double Epetra_MinDouble
long long NumGlobalCols64() const
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor's portion of the matrix.
void SetSorted(bool Flag)
int ReplaceColMap(const Epetra_BlockMap &newmap)
Replaces the current ColMap with the user-specified map object, but only if no entries have been inse...
int ExtractGlobalRowCopy(int GlobalRow, int Length, int &NumEntries, double *Values, int *Indices) const
Returns a copy of the specified global row in user-provided arrays.
int TUnpackAndCombine(const Epetra_SrcDistObject &Source, int NumImportIDs, int *ImportLIDs, int LenImports, char *Imports, int &SizeOfPacket, Epetra_Distributor &Distor, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor)
int CopyAndPermuteRowMatrix(const Epetra_RowMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
bool MyLRID(int LRID_in) const
Returns true if the LRID passed in belongs to the calling processor in this map, otherwise returns fa...
virtual int NumProc() const =0
Returns total number of processes.
bool IndicesAreGlobal() const
If matrix indices has not been transformed to local, this query returns true, otherwise it returns fa...
void UpdateExportVector(int NumVectors) const
Epetra_Map * RemoveEmptyProcesses() const
Return a new BlockMap with processes with zero elements removed.
const Epetra_BlockMap & DomainMap() const
Returns the DomainMap associated with this graph.
bool UpperTriangular() const
If matrix is upper triangular in local index space, this query returns true, otherwise it returns fal...
Epetra_IntSerialDenseVector IndexOffset_
bool GlobalIndicesTypeMatch(const Epetra_BlockMap &other) const
bool NoDiagonal() const
If matrix has no diagonal entries in global index space, this query returns true, otherwise it return...
Epetra_Map * ReplaceCommWithSubset(const Epetra_Comm *Comm) const
Replace this Map's communicator with a subset communicator.
int ExtractMyRowCopy(int LocalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified local row of the graph. Put into storage allocated by calli...
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int InvColSums(Epetra_Vector &x) const
Computes the inverse of the sum of absolute values of the columns of the Epetra_CrsMatrix, results returned in x.
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
Epetra_CrsMatrix & operator=(const Epetra_CrsMatrix &src)
Assignment operator.
double NormInf() const
Returns the infinity norm of the global matrix.
void GeneralSV(bool Upper, bool Trans, bool UnitDiagonal, double *x, double *y) const
int NumMyEntries() const
Returns the number of entries on this processor.
int ExtractMyRowView(int LocalRow, int &NumIndices, int *&Indices) const
Get a view of the elements in a specified local row of the graph.
int ExtractGlobalRowCopy(int GlobalRow, int LenOfIndices, int &NumIndices, int *Indices) const
Extract a list of elements in a specified global row of the graph. Put into storage allocated by call...
void GeneralSM(bool Upper, bool Trans, bool UnitDiagonal, double **X, int LDX, double **Y, int LDY, int NumVectors) const
int CopyAndPermuteCrsMatrix(const Epetra_CrsMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this exporter.
Epetra_SrcDistObject: A class for supporting flexible source distributed objects for import/export op...
const Epetra_Import * Importer() const
Returns the Epetra_Import object that contains the import operations for distributed operations...
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices...
void PREFIX EPETRA_DCRSMV_F77(const int *, const int *, const int *, const double *, const int *, const int *, double *, double *)
int Scale(double ScalarConstant)
Multiply all values in the matrix by a constant value (in place: A <- ScalarConstant * A)...
long long GCID64(int LCID_in) const
int * Values()
Returns pointer to the values in vector.
int MergeRedundantEntries()
Add entries that have the same column index. Remove redundant entries from list.
int ReplaceMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Replace current values with this list of entries for a given local row of the matrix.
int InsertMyValues(int MyRow, int NumEntries, const double *Values, const int *Indices)
Insert a list of elements in a given local row of the matrix.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
int TCopyAndPermuteCrsMatrix(const Epetra_CrsMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
int RightScale(const Epetra_Vector &x)
Scales the Epetra_CrsMatrix on the right with a Epetra_Vector x.
Epetra_MultiVector * ImportVector_
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Replace specified existing values with this list of entries for a given global row of the matrix...
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this exporter.
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
Returns the result of a local solve using the Epetra_CrsMatrix on a Epetra_Vector x in y...
int ExpertStaticFillComplete(const Epetra_Map &DomainMap, const Epetra_Map &RangeMap, const Epetra_Import *Importer=0, const Epetra_Export *Exporter=0, int NumMyDiagonals=-1)
Performs a FillComplete on an object that aready has filled CRS data.
Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map &RowMap, const int *NumEntriesPerRow, bool StaticProfile=false)
Epetra_CrsMatrix constructor with variable number of indices per row.
double ** Pointers() const
Get pointer to individual vector pointers.
bool Filled() const
If FillComplete() has been called, this query returns true, otherwise it returns false.
int CheckSizes(const Epetra_SrcDistObject &A)
Allows the source and target (this) objects to be compared for compatibility, return nonzero if not...
double *& ExpertExtractValues()
Returns a reference to the double* used to hold the values array.
int LRID(int GRID_in) const
Returns the local row index for given global row index, returns -1 if no local row for this global ro...
bool NoRedundancies() const
If MergeRedundantEntries() has been called, this query returns true, otherwise it returns false...
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs. ...
void GeneralMTV(double *x, double *y) const
int ReplaceColMap(const Epetra_BlockMap &newmap)
Replaces the current ColMap with the user-specified map object.
const Epetra_Map & DomainMap() const
Returns the Epetra_Map object associated with the domain of this matrix operator. ...
Epetra_DistObject: A class for constructing and using dense multi-vectors, vectors and matrices in pa...
bool Sorted() const
If SortIndices() has been called, this query returns true, otherwise it returns false.