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,numMyRows,IndexOffset)
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,
2208 int_type * Indices = 0;
2209 double * values = 0;
2212 Indices =
new int_type[maxNumEntries];
2213 values =
new double[maxNumEntries];
2221 for (i=0; i<NumSameIDs; i++) {
2222 Row = (int_type)
GRID64(i);
2232 for (i=0; i<NumSameIDs; i++) {
2233 Row = (int_type)
GRID64(i);
2245 for (i=0; i<NumSameIDs; i++) {
2246 Row = (int_type)
GRID64(i);
2253 for (i=0; i<NumSameIDs; i++) {
2254 Row = (int_type)
GRID64(i);
2265 for (i=0; i<NumSameIDs; i++) {
2266 Row = (int_type)
GRID64(i);
2276 for (i=0; i<NumSameIDs; i++) {
2277 Row = (int_type)
GRID64(i);
2289 for (i=0; i<NumSameIDs; i++) {
2290 Row = (int_type)
GRID64(i);
2297 for (i=0; i<NumSameIDs; i++) {
2298 Row = (int_type)
GRID64(i);
2309 int_type FromRow, ToRow;
2310 if (NumPermuteIDs>0) {
2314 for (i=0; i<NumPermuteIDs; i++) {
2315 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2316 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2326 for (i=0; i<NumPermuteIDs; i++) {
2327 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2328 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2340 for (i=0; i<NumPermuteIDs; i++) {
2341 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2342 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2349 for (i=0; i<NumPermuteIDs; i++) {
2350 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2351 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2362 for (i=0; i<NumPermuteIDs; i++) {
2363 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2364 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2374 for (i=0; i<NumPermuteIDs; i++) {
2375 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2376 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2388 for (i=0; i<NumPermuteIDs; i++) {
2389 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2390 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2397 for (i=0; i<NumPermuteIDs; i++) {
2398 FromRow = (int_type) A.
GRID64(PermuteFromLIDs[i]);
2399 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2420 int * PermuteToLIDs,
2421 int *PermuteFromLIDs,
2426 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteCrsMatrix: Incoming global index type does not match the one for *this",-1);
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);
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;
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);
2492 for (i=0; i<NumSameIDs; i++) {
2493 Row = (int_type)
GRID64(i);
2494 int AlocalRow = rowMap.
LID(Row);
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++) {
2508 Row = (int_type)
GRID64(i);
2514 for (i=0; i<NumSameIDs; i++) {
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];
2537 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2543 for (i=0; i<NumPermuteIDs; i++) {
2544 FromRow = PermuteFromLIDs[i];
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];
2560 ToRow = (int_type)
GRID64(PermuteToLIDs[i]);
2566 for (i=0; i<NumPermuteIDs; i++) {
2567 FromRow = PermuteFromLIDs[i];
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,
2605 throw ReportError(
"Epetra_CrsMatrix::CopyAndPermuteRowMatrix: Incoming global index type does not match the one for *this",-1);
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);
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 )
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 )
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;
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;
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);
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
3115 out = Teuchos::VerboseObjectBase::getDefaultOStream();
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,numMyRows)
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,numMyRows)
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,numMyRows)
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);
3519 #ifdef EPETRA_HAVE_OMP
3520 #pragma omp parallel for default(none) shared(IndexOffset,Indices,values,NumVectors,numMyRows,xp,yp)
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();
3548 #ifdef EPETRA_HAVE_OMP
3549 #pragma omp parallel for default(none) shared(NumEntriesPerRow,Indices,srcValues,NumVectors,numMyRows,xp,yp)
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]];
3569 #ifdef EPETRA_HAVE_OMP
3570 #pragma omp parallel for default(none) shared(NumVectors,numMyRows,xp,yp)
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;
3788 for (i=NumMyRows_-1; i >= 0; i--) {
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;
3929 for(i = NumMyRows_ - 1; i >= 0; i--) {
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;
4029 for(i = NumMyRows_ - 1; i >= 0; i--) {
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;
4564 for(i = NumMyRows_ - 1; i >= 0; i--) {
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;
4677 for(
int i=0;i<m;i++){
4718 for(
int i=0; i<m; i++){
4723 if(NumIndices > 0) {
4726 int jl_0 = col_indices[0];
4727 int jl_n = col_indices[NumIndices-1];
4733 if(numMyDiagonals == -1) {
4737 int insertPoint = -1;
4739 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
4748 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
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 GlobalConstantsComputed_
int NumMyRowEntries(int MyRow, int &NumEntries) const
Return the current number of values stored for the specified local row.
const Epetra_Export * Exporter_
int TransformToLocal()
Use FillComplete() instead.
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector &x, Epetra_Vector &y) const
const Epetra_BlockMap & RangeMap() const
Returns the RangeMap associated with this graph.
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 ReplaceDiagonalValues(const Epetra_Vector &Diagonal)
Replaces diagonal values of the matrix with those in the user-provided vector.
bool HaveColMap() const
Returns true if we have a well-defined ColMap, and returns false otherwise.
bool DistributedGlobal() const
Returns true if this multi-vector is distributed global, i.e., not local replicated.
Epetra_Map: A class for partitioning vectors and matrices.
bool DistributedGlobal() const
Returns true if map is defined across more than one processor.
const Epetra_Map & RangeMap() const
Returns the Epetra_Map object associated with the range of this matrix operator.
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Epetra_BlockMap RangeMap_
Epetra_IntSerialDenseVector: A class for constructing and using dense vectors.
int FillComplete()
Tranform to local index space. Perform other operations to allow optimal matrix operations.
void GeneralMM(double **X, int LDX, double **Y, int LDY, int NumVectors) const
int GlobalMaxNumNonzeros_
int ElementSize() const
Returns the size of elements in the map; only valid if map has constant element size.
bool StorageOptimized() const
If OptimizeStorage() has been called, this query returns true, otherwise it returns false...
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 MaxValue(double *Result) const
Compute maximum value of each vector in multi-vector.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.
long long NumGlobalRows64() const
int Resize(int Length_in)
Resize a Epetra_IntSerialDenseVector object.
int SumIntoOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices)
bool IndicesAreContiguous_
int NumMyEntries(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
bool GlobalIndicesLongLong() const
Returns true if map create with long long NumGlobalElements.
Epetra_Distributor: The Epetra Gather/Scatter Setup Base Class.
bool SameAs(const Epetra_BlockMap &Map) const
Returns true if this and Map are identical maps.
int * Values_alloc_lengths_
int ExtractDiagonalCopy(Epetra_Vector &Diagonal) const
Returns a copy of the main diagonal in a user-provided vector.
Epetra_OffsetIndex: This class builds index for efficient mapping of data from one Epetra_CrsGraph ba...
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 ReplaceOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices)
int * IndexOffset() const
int SortEntries()
Sort column entries, row-by-row, in ascending order.
const Epetra_BlockMap & TargetMap() const
Returns the TargetMap used to construct this exporter.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer...
void DecrementReferenceCount()
Decrement reference count.
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
virtual void Print(std::ostream &os) const
Print object to an output stream.
long long NumGlobalEntries_
const Epetra_BlockMap & SourceMap() const
Returns the SourceMap used to construct this exporter.
double * All_Values() const
void GeneralMTV(double *x, double *y) const
long long NumGlobalElements64() const
long long NumGlobalNonzeros_
int TCopyAndPermuteRowMatrix(const Epetra_RowMatrix &A, int NumSameIDs, int NumPermuteIDs, int *PermuteToLIDs, int *PermuteFromLIDs, const Epetra_OffsetIndex *Indexor, Epetra_CombineMode CombineMode)
long long NumGlobalCols64() const
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 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 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 NumMyEntries() const
Returns the number of entries on this processor.
double NormInf() const
Returns the infinity norm of the global matrix.
int MyLength() const
Returns the local vector length on the calling processor of vectors in the multi-vector.